BQ-Chem : A Quantum Software Program for Chemistry Simulation Based on the Full Quantum Eigensolver Algorithm

,


Introduction
Quantum chemistry is also called molecular quantum mechanics, where quantum mechanics is applied to the feld of chemistry to obtain the chemical properties of molecules at the atomic level. A central challenge in the feld of quantum chemistry is to determine the low-lying molecular energies and electronic structures of a chemical system with high precision, and its electronic Hamiltonian usually can be obtained with a set of nuclear geometries under the Born-Oppenheimer approximation. Moreover, by varying nuclear positions of a chemical system, a potential energy surface (e.g., ground-state energy) can be determined as a function of the bond length or bond angle of interest, which is key to understanding chemical reactivity, product distributions, and reaction rates.
Conventionally, the number of classical bits needed to simulate a multiatom molecule increases exponentially as the number of atoms. And in the worst case, quantum simulation of a chemical system is exponentially hard on classical computers. Despite the great success of approximation methods, tackling the problem accurately with conventional computers is still a difcult task. On the other hand, the chemistry of molecules is accurately described with quantum mechanics; thus, the quantum simulation of a chemical system with a quantum simulator or quantum computer has natural advantages in the sense of computational resources needed [1,2].
However, the problem of estimating the ground state and its energy of Hamiltonian is in general QMA-hard [3] and even quantum computers are not expected to efciently solve it. In the past few decades, various quantum algorithms are harnessed to address this issue. Te main idea of these methods is to prepare a state which extremely approximates to the desired ground state instead of preparing the exact ground state which is computationally expensive. Based on the well-known adiabatic theorem [4,5] and quantum phase estimation algorithm [6], Aspuru-Guzik et al. [7,8] constructed an adequate quantum circuit to evolve the system adiabatically in the ground state, i.e., adiabatic state preparation. Te quantum imaginary time evolution algorithm [9,10], which uses the linear combination of unitaries (LCU) [11] formalism, is another approach as a heuristic method to approximate the ground state.
Despite developments in quantum algorithms and optimization of resource requirements, many of the algorithms have hardware requirements far beyond the capability of near-term quantum computers. In 2014, a quantum-classical hybrid algorithm called variational quantum eigensolver (VQE) [12,13] was proposed by Peruzzo et al. VQE can prepare a class of ansatz states with a variational quantum circuit and optimize its parameters, leading the state approximating to the ground state step by step. VQE requires fewer quantum gates and a shallower circuit, making it more feasible on noisy intermediate-scale quantum (NISQ) [14] devices. In 2020, Wei et al. developed a full quantum eigensolver (FQE) [15] scheme which attracted much attention in recent years [16][17][18][19][20][21]. Compared with VQE where its optimization step is completed on a classical computer, FQE, using the quantum gradient descent algorithm in the LCU paradigm, is a fully quantum algorithm where the calculation of the expectation value of Hamiltonian and the optimization process is fully conducted on quantum computers, enabling it an advantage of providing a smooth transition to future large-scale quantum computers.

BQ-Chem
Here, we give an introduction to a newly designed quantum chemistry simulation software program BQ-Chem (https:// biwonq.baqis.ac.cn/#/pages/chemistry), which is based on FQE. Tis quantum chemistry simulation module supports multitask parallel computing, and users can edit their specifc contents after adding projects on the main page.

Specifying a Molecule and Its Parameters.
In order to simplify the process and improve the user's experience of the program, we preinstalled the specifc molecular confgures of about 50 common molecules for users as default in BQ-Chem, including the atomic type, the x, y, and z coordinate of each atomic nucleus, and the bonding properties between them. Te default data of chemical compositions and their parameters about our preinstalled molecules are extracted from the National Institute of Standards and Technology (NIST) Chemistry WebBook (https://webbook.nist.gov/chemistry/) of the United States. In addition, users are also allowed to build their own molecular models fexibly by adding more atoms and changing atomic coordinates and corresponding bonding properties, as shown in Figure 1, where we take NH 3 molecule as an example of specifying a molecule and its parameters in BQ-Chem.

Two Calculation Modes.
After specifying a molecule and its parameters from the last step, users can calculate the chemical properties of the molecule of their interest. Te BQ-Chem simulation module supports two calculation modes for users to choose from: the static single-point energy spectrum (SES) calculation mode where all molecular geometry parameters given are fxed and the dynamic potential energy surface (PES) calculation mode where some molecular parameters of the user's interest range around the given geometry. As the following shows, the SES mode supports calculating the ground and excited energies sequentially, while the PES mode supports calculating the ground energy around the given geometry by varying the bond length or bond angle of a molecule. In both modes, one should frst set the preprocessing information, including basis set choice, fermion-spin transformation model choice, and the initial state choice.

Single-Point Energy Spectrum (SES) Mode.
In the mode of SES, one can choose to calculate the ground-state energy with FQE and up to 15 excited-level energies one by one by subtracting the ground-state energy part from the system Hamiltonian and calculating the new ground-state energy of the subtracted Hamiltonian.
BQ-Chem provides several classical methods for preprocessing molecular electronic structures for users. As depicted in Figure 2, the following steps are taken: (1) Choosing the single-electron atomic orbital basis set from the option of "sto-3g" and "sto-6g." Here, we take "sto-ng" atomic orbitals instead of the original "sto" or "gto" basis, which means the basis functions are expressed as sums of Gaussian functions rather than the original Slater-type orbitals to enhance the efciency of integral evaluation. More details are described in Appendix A. (2) Choosing the Fermi-spin transformation method from Jordan-Wigner or Bravyi-Kitaev transformation. Te Jordan-Wigner transformation directly maps the fermionic occupation state of a particular atomic orbital to a state of qubits, while the Bravyi-Kitaev transformation is a compromise between the Jordan-Wigner transformation and the parity encoding transformation, which is more efcient in operational complexity. More details are described in Appendix C. (3) Choosing the initial state of the full quantum eigensolver solution process from the option of the Hartree-Fock approximation state or the uniform superposition state. Te Hartree-Fock approximation state usually is a better choice for faster convergence. More details are described in Appendix B.
We take NH 3 molecule as an example; its confguration and setting of the SES mode are shown in Figure 2, i.e., Hartree-Fock state as the initial state, "sto-6g" as the basis, Jordan-Wigner transformation method as the chosen Fermi-spin transformation, and the goal is to obtain the ground energy and up to three excitation energies.

Quantum Engineering
After executing the process, a resultant spectrum of 4 energy levels is presented in Figure 3, where we can see as the iteration of the quantum gradient descent algorithm increases, the simulated ground-state and up to three excitation energies have a convergence to the theoretical energies of the interested Hamiltonian, respectively. To be noted, the ground-state energy can be obtained directly from the original Hamiltonian with FQE, while the 1st excitation energy is acquired from a new Hamiltonian by reducing the ground-state energy part of the original Hamiltonian with the same FQE process. By repeating the above process, one can reach any excited-level energy one after the other.

Potential Energy Surface (PES) Mode.
In the mode of PES, one can choose the scanning mode from "bond length" and "bond angle"; moreover, one can set the scanning number and scanning feld, as depicted in Figure 4. If the "bond length" option is selected, one is further required to provide which set of chemical bonds is going to be modifed. Otherwise, if the "bond angle" option is selected, one is further required to specify two chemical bonds which share a common atom (as the vertex of two connected edges).
We take NH 3 molecule as an example again. Its confguration and setting of the PES mode are shown in Figure 4, i.e., "bond length" as the scanning option, all chemical bonds of NH 3 are specifed, and 11 scanning points of chemical bond length set ranging from 50% to 150% of the original chemical bond length and its setting are maintained the same with the SES mode.
PES can help search for the most stable structure of the interested molecule by ranging the length and angles of the chemical bonds, which is very useful in the feld of biomedicine, materials, and chemical industry. In the example of calculating PES of the NH 3 molecule, the resultant PES as presented in Figure 5 shows that the molecule with the original set of chemical bond lengths has the lowest energy which is the most stable structure of NH 3 molecule under the fxed chemical bond angles.

Full Quantum Eigensolver
In this section, we will give a detailed description of the underlying algorithm in the BQ-Chem module, FQE. A diagram of FQE for solving the ground-state energy of a molecule is depicted in Figure 6. As discussed above, the word "full" for FQE [15] means that all computations of FQE are realized on quantum computers, while other eigensolvers require a sequence of data transformation between quantum computers and classical computers. Te reason why FQE can be fully completed on a quantum computer is that it takes advantage of the technology of linear combination of unitary operators (LCU) which is proposed by Long in 2002 [11]. : Input interface of the NH 3 molecular data as an example in the BQ-Chem simulation module. One can select a molecule structure from a preset database and modify the structure by changing the coordinates of its element. One can also construct another molecule structure by adding or deleting an atom from the selected molecular structure.
LCU can realize the four arithmetic operations, addition, subtraction, multiplication, and division of unitary operators.
After the secondary quantization, the fermion Hamiltonian of the molecule is transformed into the qubit Hamiltonian in the Hilbert space through mathematical mapping. Te ground state of the Hamiltonian can be obtained by applying the FQE. In a specifc quantum computing system, we shall choose a suitable initial quantum state, then construct a quantum circuit realizing the quantum gradient descent process, and measure the fnal quantum state to obtain the information that the next loop is required. In the process of continuous iteration of the quantum circuit, the fnal state converges to the ground state of our interested Hamiltonian. On the other hand, searching for the most stable structure of a molecule is very difcult in quantum chemistry and requires massive classical computation. Similar to solving the ground state energy problem, FQE can also be used to fnd the state with the lowest energy by changing the distance between atoms or angles between chemical bonds in a molecule and obtain the most stable structure of the molecule.
A molecular system containing nuclei and electrons can be described by its molecular Hamiltonian. Trough Jordan-Wigner (JW) or Bravyi-Kitaev (BK) transformation, the molecular Hamiltonian can be mapped to the Hamiltonian in a qubit form.
where the Roman indices i, j denote the qubit on which the operators act and Greek indices α, β refer to the type of Pauli operators; for example, σ i x means Pauli matrix σ x acting on the i -th qubit. Apparently, H in equation (1) is a linear combination of unitary Pauli matrices. We can calculate the ground-state energy by minimizing the expected value of the Hamiltonian.
Te following methods for fnding the molecular ground state and its energy are all based on it.
In classical computing, we usually use the gradient descent algorithm to obtain the minimum value of the objective function f(X → ). We start from the initial state X →(0) ∈ R N and evolve along the object function's gradient direction to the next state: In the case of solving the ground-state energy of the Hamiltonian problem with FQE, the gradient of f(X → ) can be expressed as follows: And, the iterative process of classical gradient descent can be mapped into a quantum version, which can be considered as the evolution of the quantum state X → under the operator H.
We defne a new Hamiltonian as follows: where M is the number of Pauli terms of operator H g . Te process of the quantum gradient descent algorithm can be expressed as follows: s H g is an LCU operator, which is a nonunitary evolution and can be realized on a quantum computer by adding auxiliary qubits to convert it into a unitary evolution in a larger Hilbert space. Te entire quantum circuit diagram described by equation (7) is divided into four parts as shown in Figure 7:  Quantum Engineering Wave division: Te entire quantum system is composed of a work system and an ancillary system. Firstly, we encode the initial vector X � (x 1 , . . . , x N ) T into the initial state | x (t) 〉 of the work system. In the feld of quantum chemistry, the Hartree-Fock (HF) state is often adopted as an initial state. Te ancillary system is initialized from |0 m , where m � log 2 M, to a specifc superposition as follows:  Quantum Engineering Here, C � ������� M− 1 i�0 β 2 i is a normalization constant and | i〉 is the computational basis. We denote the composite state of the whole system as |Φ � |ψ s 〉| x →(t) 〉.

Entanglement: A series of ancilla-controlled operations
on the work qubits are implemented, the work qubits and the ancillary qubits are now entangled, and the state is transformed into Wave combination: After performing m Hadamard gates on the ancillary register, the state of the whole system in the subspace where the ancillary system in state | 0〉 is Measurement: Here, we make measurements on the ancillary register. If | 0〉 is detected, our algorithm succeeds and we obtain the state is proportional to the next iterative state. Te probability of success in detecting the ancillary state in | 0〉 is After successfully obtaining | 0〉 in the ancillary system, we can continue the gradient descent process by repeating the above four steps. We can preset a threshold of → t 〉 as the criterion for stopping the iteration. After enough iterations, the fnal state of work system | X〉 will converge to our interested Hamiltonian H's ground state and 〈X|H|X〉 is the corresponding ground state energy.
Te successful probability after n-time measurements is 1 − ((C 2 M − ‖H g | x →(t) 〉‖ 2 )/C 2 M) n , which is an exponential function of n. Te number of measurements is C 2 M/‖H g | x →(t) 〉‖ 2 , which shows its complexity will grow exponentially with respect to the number of iteration steps [22], which is a relatively large limitation and should be solved in future work.

Summary
In   Hamiltonian after reducing the part corresponding to its ground state following the same procedure. In PES mode, by varying the distances between atoms or angles between the chemical bonds in a molecule, we can search for the lowest energy and the corresponding most stable structure of a molecule with a set of distances or angles, from its potential energy surface corresponding to diferent distances or angles.

A. Basis Set
In quantum chemical calculations based on the self-consistent feld method (SCF), the single-electron orbital basis vectors of molecules are generally approximated by the linear combination of the single-electron orbital basis vectors (basis set) of each constituent atom. We know that for a hydrogen-like atom with nuclear electron number Z, its single-electron wave function takes the form where the radial function and a μ � 4πε 0 ℏ 2 /μe 2 [23]. Since the mass of the nucleus can be regarded as infnite relative to the mass of the electron, the reduced mass can be approximately taken as the mass of the electron, i.e., μ � m n m e /(m n + m e ) ≈ m e ; then, a μ is the corresponding Bohr radius. A generalized Laguerre function can be expressed as follows:

Besides, the angular function in equation (A.1) is
is the imaginary symbol, and the adjoint Legendre polynomial has the form which becomes the basis set of Slater-type orbitals function (sto). We have the following: (i) N � (2ζ) n ������� 2ζ/(2n)! is the normalization coefcient (ii) n is the principal quantum number (iii) r is the radial distance between the electron and the atomic kernel (iv) ζ is a constant related to the efective nuclear charge determined by Slater's rule After substituting equations (A.4) and (A.6) into equation (A.1), the "sto" basis set is obtained. As we can see, "sto" has a natural physical meaning; thus, the approximation is better.

A.2. Gaussian-Type Orbitals (gto)
. Te Gaussian-type orbitals (gto) set is another typical atomic electron orbital basis set. Boys [26,27] frst reduced the radial wave function in (A.2) to where N, n, and r take the same defnition as equation (A.6) and α is a constant number related to the efective nuclear charge.
According to the Gaussian Product Teorem, the product of any two Gaussian functions (gtos) centered at two diferent positions can be replaced by a fnite sum of gtos centered on a point along the axis connecting them. When we select the "gto" basis set, the molecular electronic state of the n-atom molecule can be expressed as the product of n "gto" atomic electronic states at n nuclei centers, respectively. According to the Gaussian Reproduction Lemma, the product can be further expressed as the linear superposition of Gaussian integral function production with n/2 centers and again further expressed as the linear superposition of multiple Gaussian integral function productions with the same center. In other words, the computational efciency of the "gto" basis set is much higher than that of the sto basis set, and this advantage can sometimes reach 4-5 orders of magnitude. However, by comparing equation (A.6) with (A.7), as shown in Figure 8, it can be found that the two function properties near r � 0 are quite diferent, which means that although the computational efciency of "gto" is higher than that of "sto," its computational accuracy is lower than that of "sto." A.3. "sto-ng" Atomic Orbitals. In order to solve the problem that the Gaussian function does not match the linearity of real atomic orbitals, a "gto" orbit can be approximated by the linear superposition of multiple Gaussian orbits with different parameters α [28,29]. Tese "sto" orbitals participating in the superposition are called the primitive Gaussian basis set, and their linear superposition is called the contracted Gaussian basis set, where the superposition coefcients are obtained and fxed by prior optimization. It remains unchanged in the subsequent variational 8 Quantum Engineering optimization process for fnding the ground state of the molecule; i.e., the superposition coefcient is not used as an optimization parameter. Generally, the approximate "sto" atomic orbital basis set obtained by the linear superposition of n Gaussian basis vectors is denoted by "sto-ng." For example, as depicted in Figure 9, "sto-3g" is a linear superposition of three Gaussiantype atomic orbitals to approximate a Slater-type atomic orbital.

B. Hartree-Fock Self-Consistent Field Method
In Appendix A, we get the basis set ϕ i k (r) , where ϕ i k (r) represents the i-th single-electron atomic orbital of the k-th atom.
A single-electron molecular orbital ϕ j (r) � k c j k,i φ i k (r) can be obtained by linearly superposing diferent one-electron atomic orbitals. After the direct product of diferent single-electron molecular orbitals and the commutative antisymmetry operation, a multielectron molecular orbital that satisfes fermion statistics can be obtained.
Under the Born-Oppenheimer approximation, we regard the nucleus as stationary and only provide a background electrostatic feld for the electrons; then, the multielectron Hamiltonian of the molecule can be expressed as follows: H 0 is the potential energy of the kinetic energy of the electron under the background electric feld, that is, the single energy, and H I describes the Coulomb potential energy between the electrons. Te energy of the corresponding manybody quantum state Ψ is expressed as follows: For the zero-order approximation E 0 , since the direct product state |φ 1 |φ 2 · · · |φ n 〉 is an eigenstate of H 0 , E 0 can be directly described as the summation of single-particle energies.
H I is the interaction Hamiltonian, expressed as

(B.5)
In principle, the corresponding many-body ground state and ground state energy can be obtained by calculating the minimum of equation (B.3) with calculus of variations. We know that the many-body state | Ψ〉 is obtained by the Slater symmetrization of the single-electron molecular state | φ i 〉, and then, the variation of equation (B.3) is actually the variation of | φ i 〉. Considering the normalization condition of the one-electron molecular state, equations (B.5) and (B.5) can be substituted into equation (B.3), and the Lagrangian multiplier term can be matched at the same time: After a variation of the whole system on a certain state 〈φ i | and taking the extreme δE/δ〈φ i | � 0, there should be which is the Hartree-Fock self-consistent feld equation [30,31]. For the Hermitian matrix ϵ i,j , there always exists a unitary transformation U that can diagonalize it, and then, if we take n j�1 U i,j |ϕ j 〉 as the new φ i , a regular form Hartree-Fock equation can be expressed as Slater-type 1s orbital Gaussian-type 1s orbital Quantum Engineering 9 where F i is the i-th row of the Fock operator F that has a dependency on the single-electron state |φ i 〉 , ε i is the i-th eigenvalue of F, and is the monomer energy of the i-th single-electron state, is the Coulomb repulsion energy of the i-th single-electron state corresponding to other electrons, and is the Fermi exchange energy between the i-th single-electron state and other electronic states. Te variational optimization process of the Hartree-Fock method can be depicted as the fowchart in Figure 10: (1) Start from the preset single-electron molecular state φ i (2) Calculate the corresponding Fock operator F( φ i ) according to equations (B.9)-(B.11) (3) Diagonalize the Fock operator F to get the corresponding eigenvectors φ i′ and eigenvalues (4) Continue to substitute the new single-electron molecular state as in (2) to obtain a new F ′ until convergence

C. Fermi-Spin Transformation Method
Quantum simulation of chemistry is always dealing with fermions that satisfy the Pauli exclusion principle (i.e., electrons) [32,33]. We can simply handle the electrons of spin ↑↓ as handling two diferent kinds of spin-free fermions without spin-orbit coupling. In this section, spin-free fermions are used as an example to introduce related calculation methods. Since fermions satisfy commutative antisymmetry, However, as the basic storage unit of a quantum computer, the qubit is a boson that satisfes commutative symmetry, i.e., SU (2) spin, which satisfes Tis shows that if only the particles on the same lattice are considered, the Fock state basis vector |0, |1〉 of the Fermi system and the spin basis vector |↓, |↑〉 can be one-to-one correspondence, and furthermore, its Hamiltonian can be related to each other with the same transformation.
However, the situation could be quite diferent when dealing with multiple-lattice Fermi systems, the structure obtained from equation (C.3) loses the commutative antisymmetry property of fermions, which will cause confusion of sign in the calculation, leading to completely wrong results. In that sense, a unitary operation that faithfully connects the fermionic Fock space with the spin space is required. Obviously, such a unitary transformation cannot be local.

C.1. Jordan-Wigner Transformation.
Considering the fermions as a one-dimensional chain, we can denote the fermion lattice as i � 0, 1, . . . , n − 1 based on the algebras of fermion particles. In order to unify the sign of Fermi-spin  operators on the j-th lattice, we can just calculate the parity of all fermions occupied before the j-th lattice and then assign it to the local operators of the j-th lattice, which is exactly the Jordan-Wigner transformation and its inverse is shown in equation (C.4). We can simply verify that the Jordan-Wigner transformation [34][35][36] is a global unitary transformation.
Because of the existence of string operators ( i<j (1 − 2f † i f i )) and ( i<j σ z j ) in equation (C.4), the averaged number of spin operators required to encode a fermion operator is linear to the system size under the Jordan-Wigner transformation; i.e., the Jordan-Wigner transformation has an encoding complexity of O(n).

C.2. Bravyi-Kitaev Transformation.
When encoding the fermionic operators with spin operators, two pieces of information should be included, the fermionic occupation on i -th lattice and the fermionic parity of i-th lattice (i.e., the parity of the fermionic occupation numbers of all lattices preceding i-th lattice). Te Jordan-Wigner transformation encodes the fermionic occupation information with the local spin operator σ z j and encodes the fermionic occupation number parity information with the nonlocal string operator i<j σ z j . Another protocol called "Parity encoding" is in the opposite encoding way, where the local spin operator at j-th lattice encodes into the parity of the fermionic occupation numbers at all lattices up to and including lattice number j and encodes the spin string operator into the occupation state at the j-th fermionic lattice. Similar to the Jordan-Wigner transformation, the parity encoding of the Bravyi-Kitaev transformation [37][38][39] also has a linear averaged complexity of O(n).
In order to reduce the encoding complexity of fermion operators in the spin formula, Bravyi and Kitaev [40] proposed another Fermi-spin transformation method in 2002, namely, the Bravyi-Kitaev transformation, which is a compromise between the Jordan-Wigner transformation and the parity encoding transformation. Te Bravyi-Kitaev method draws on the binary tree structure to encode the information of the Fermi subsystem on the spin system, as shown in Figure 11, where the solid gray lines represent the left branches, the dashed black lines represent the right branches, and the solid red circles are the roots. At the bottom layer of Figure 11, the spin state at the j-th (i.e., even j) lattice on the left branch stores the local occupation information of the fermion chain at the same j-th lattice and the right branch at the j-th (i.e., odd j) lattice stores the fermion occupation number parity information for some fnite sets of the lattices, which conclude all lattices under the root of the highest level that can be reached along the right branch route starting from the j-th lattice.
We take 8 lattices as an example in Figure 11; the Bravyi-Kitaev transformation can be expressed as follows: where ⊕ is a symbol of modulo-2 binary addition. Obviously, β 2 n β † 2 n ≠ I means that β 2 n is not a unitary transformation; nevertheless, it is always linearly invertible.
By applying the conjugate of the transformation matrix β to the fermionic Hamiltonian, the spin Hamiltonian with the same energy spectrum can be obtained and thus solved on a quantum computer. Under the Bravyi-Kitaev transformation, the fermionic single lattice occupation number information and the parity information are stored neither locally in the single lattice of spins nor in the string operator with linear complexity, but in the string operator with logarithmic length. Tis means the averaged encoding complexity of the fermionic creation and annihilation operators under the Bravyi-Kitaev transformation is O(log n).

Data Availability
Te data sets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest
Te authors declare that they have no conficts of interest.