Abstract
We describe a quantum chemistry simulation software program BQ-Chem, which can calculate the low-energy spectrum and potential energy surface of molecules on a quantum computer. BQ-Chem is based on the full quantum eigensolver (FQE), which is implemented with a quantum gradient descent algorithm. Benefiting from FQE, BQ-Chem can perform all the calculations on a quantum computer. Compared with the classical optimization methods which encounter the optimization difficulty of high-dimensional and multivariable functions in dealing with multielectron orbitals of macromolecules, FQE provides an exponential speedup. FQE works fully on a quantum computer; thus, BQ-Chem can be smoothly transited to future large-scale quantum computers.
1. Introduction
Quantum chemistry is also called molecular quantum mechanics, where quantum mechanics is applied to the field of chemistry to obtain the chemical properties of molecules at the atomic level. A central challenge in the field 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 difficult 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 efficiently solve it. In the past few decades, various quantum algorithms are harnessed to address this issue. The 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. The 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–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.
2. 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. This quantum chemistry simulation module supports multitask parallel computing, and users can edit their specific contents after adding projects on the main page.
2.1. Specifying a Molecule and Its Parameters
In order to simplify the process and improve the user’s experience of the program, we preinstalled the specific molecular configures of about 50 common molecules for users as default in BQ-Chem, including the atomic type, the , , and coordinate of each atomic nucleus, and the bonding properties between them. The 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 flexibly by adding more atoms and changing atomic coordinates and corresponding bonding properties, as shown in Figure 1, where we take molecule as an example of specifying a molecule and its parameters in BQ-Chem.
2.2. 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. The 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 fixed 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 first set the preprocessing information, including basis set choice, fermion-spin transformation model choice, and the initial state choice.
2.2.1. 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 efficiency of integral evaluation. More details are described in Appendix A.(2)Choosing the Fermi-spin transformation method from Jordan–Wigner or Bravyi–Kitaev transformation. The 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 efficient 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. The Hartree–Fock approximation state usually is a better choice for faster convergence. More details are described in Appendix B.
We take molecule as an example; its configuration 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.
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.
2.2.2. 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 field, 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 modified. 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 molecule as an example again. Its configuration and setting of the PES mode are shown in Figure 4, i.e., “bond length” as the scanning option, all chemical bonds of are specified, and 11 scanning points of chemical bond length set ranging from to 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 field of biomedicine, materials, and chemical industry. In the example of calculating PES of the 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 molecule under the fixed chemical bond angles.
3. 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. The 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]. 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. The ground state of the Hamiltonian can be obtained by applying the FQE. In a specific 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 final quantum state to obtain the information that the next loop is required. In the process of continuous iteration of the quantum circuit, the final 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 difficult in quantum chemistry and requires massive classical computation. Similar to solving the ground state energy problem, FQE can also be used to find 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. Through 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 denote the qubit on which the operators act and Greek indices refer to the type of Pauli operators; for example, means Pauli matrix acting on the -th qubit. Apparently, 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.
The following methods for finding 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 . We start from the initial state 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 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 under the operator .where .
We define a new Hamiltonian as follows:where is the number of Pauli terms of operator . The process of the quantum gradient descent algorithm can be expressed as follows:
s 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. The entire quantum circuit diagram described by equation (7) is divided into four parts as shown in Figure 7: Wave division: The entire quantum system is composed of a work system and an ancillary system. Firstly, we encode the initial vector into the initial state of the work system. In the field of quantum chemistry, the Hartree–Fock (HF) state is often adopted as an initial state. The ancillary system is initialized from , where , to a specific superposition as follows: Here, is a normalization constant and is the computational basis. We denote the composite state of the whole system as . 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 Hadamard gates on the ancillary register, the state of the whole system in the subspace where the ancillary system in state is Measurement: Here, we make measurements on the ancillary register. If is detected, our algorithm succeeds and we obtain the state , which is proportional to the next iterative state. The probability of success in detecting the ancillary state in is
After successfully obtaining in the ancillary system, we can continue the gradient descent process by repeating the above four steps. We can preset a threshold of as the criterion for stopping the iteration. After enough iterations, the final state of work system will converge to our interested Hamiltonian ’s ground state and is the corresponding ground state energy.
The successful probability after -time measurements is , which is an exponential function of . The number of measurements is , 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.
4. Summary
In this paper, we present a quantum chemical simulation software program based on FQE and BQ-Chem. With BQ-Chem, we can calculate the single-point energy spectrum (SES) and potential energy surface (PES) after specifying the configuration and setting of a molecule. In the SES mode, we can obtain the ground-state energy directly with FQE. In addition, we can calculate the excited states and energies by replacing the interested Hamiltonian with the original 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 different distances or angles.
Appendix
A. Basis Set
In quantum chemical calculations based on the self-consistent field 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 , its single-electron wave function takes the formwhere the radial functionand [23]. Since the mass of the nucleus can be regarded as infinite relative to the mass of the electron, the reduced mass can be approximately taken as the mass of the electron, i.e., ; then, is the corresponding Bohr radius. A generalized Laguerre function can be expressed as follows:
Besides, the angular function in equation (A.1) iswhere is the imaginary symbol, and the adjoint Legendre polynomial has the form
A.1. Slater-Type Orbitals (sto)
Generally speaking, the single-electron radial wave function of a hydrogen-like atom of equation (A.2) is extremely complex and difficult to handle in numerical calculations. For this reason, John C. Slater simplified the radial function part of equation (A.1) approximately [24, 25] as follows:which becomes the basis set of Slater-type orbitals function (sto). We have the following:(i) is the normalization coefficient(ii) is the principal quantum number(iii) is the radial distance between the electron and the atomic kernel(iv) is a constant related to the effective 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)
The Gaussian-type orbitals (gto) set is another typical atomic electron orbital basis set. Boys [26, 27] first reduced the radial wave function in (A.2) towhere , and take the same definition as equation (A.6) and is a constant number related to the effective nuclear charge.
According to the Gaussian Product Theorem, the product of any two Gaussian functions (gtos) centered at two different positions can be replaced by a finite 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 -atom molecule can be expressed as the product of “gto” atomic electronic states at 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 centers and again further expressed as the linear superposition of multiple Gaussian integral function productions with the same center. In other words, the computational efficiency 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 are quite different, which means that although the computational efficiency 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]. These “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 coefficients are obtained and fixed by prior optimization. It remains unchanged in the subsequent variational optimization process for finding the ground state of the molecule; i.e., the superposition coefficient is not used as an optimization parameter.
Generally, the approximate “sto” atomic orbital basis set obtained by the linear superposition of Gaussian basis vectors is denoted by “sto-ng.” For example, as depicted in Figure 9, “sto-3g” is a linear superposition of three Gaussian-type atomic orbitals to approximate a Slater-type atomic orbital.
B. Hartree–Fock Self-Consistent Field Method
In Appendix A, we get the basis set , where represents the -th single-electron atomic orbital of the -th atom. A single-electron molecular orbital can be obtained by linearly superposing different one-electron atomic orbitals. After the direct product of different single-electron molecular orbitals and the commutative antisymmetry operation, a multielectron molecular orbital that satisfies fermion statistics can be obtained.where is the commutation operation, and equation (B.1) is the Slater determinant.
Under the Born–Oppenheimer approximation, we regard the nucleus as stationary and only provide a background electrostatic field for the electrons; then, the multielectron Hamiltonian of the molecule can be expressed as follows:
is the potential energy of the kinetic energy of the electron under the background electric field, that is, the single energy, and describes the Coulomb potential energy between the electrons. The energy of the corresponding many-body quantum state is expressed as follows:
For the zero-order approximation , since the direct product state is an eigenstate of , can be directly described as the summation of single-particle energies.
is the interaction Hamiltonian, expressed as
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 , and then, the variation of equation (B.3) is actually the variation of . 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 and taking the extreme , there should bewhich is the Hartree–Fock self-consistent field equation [30, 31]. For the Hermitian matrix , there always exists a unitary transformation that can diagonalize it, and then, if we take as the new , a regular form Hartree–Fock equation can be expressed aswhere is the -th row of the Fock operator that has a dependency on the single-electron state , is the -th eigenvalue of , andis the monomer energy of the -th single-electron state,is the Coulomb repulsion energy of the -th single-electron state corresponding to other electrons, andis the Fermi exchange energy between the -th single-electron state and other electronic states.
The variational optimization process of the Hartree–Fock method can be depicted as the flowchart in Figure 10:(1)Start from the preset single-electron molecular state (2)Calculate the corresponding Fock operator according to equations (B.9)–(B.11)(3)Diagonalize the Fock operator to get the corresponding eigenvectors and eigenvalues(4)Continue to substitute the new single-electron molecular state as in (2) to obtain a new 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 different 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 satisfies commutative symmetry, i.e., SU (2) spin, which satisfies
By comparing equation (C.1) with equation (C.2) and combining the occupation properties of the fermion state and the state directivity of the spin, we can form an equation between the single-fermion operators and and the − SU (2) spin operators , , as follows:
This shows that if only the particles on the same lattice are considered, the Fock state basis vector 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 different 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 based on the algebras of fermion particles. In order to unify the sign of Fermi-spin operators on the -th lattice, we can just calculate the parity of all fermions occupied before the -th lattice and then assign it to the local operators of the -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–36] is a global unitary transformation.
Because of the existence of string operators and 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 .
C.2. Bravyi–Kitaev Transformation
When encoding the fermionic operators with spin operators, two pieces of information should be included, the fermionic occupation on -th lattice and the fermionic parity of -th lattice (i.e., the parity of the fermionic occupation numbers of all lattices preceding -th lattice). The Jordan–Wigner transformation encodes the fermionic occupation information with the local spin operator and encodes the fermionic occupation number parity information with the nonlocal string operator . Another protocol called “Parity encoding” is in the opposite encoding way, where the local spin operator at -th lattice encodes into the parity of the fermionic occupation numbers at all lattices up to and including lattice number and encodes the spin string operator into the occupation state at the -th fermionic lattice. Similar to the Jordan–Wigner transformation, the parity encoding of the Bravyi–Kitaev transformation [37–39] also has a linear averaged complexity of .
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. The 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 -th (i.e., even ) lattice on the left branch stores the local occupation information of the fermion chain at the same -th lattice and the right branch at the -th (i.e., odd ) lattice stores the fermion occupation number parity information for some finite 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 -th lattice.
Under the Bravyi–Kitaev transformation, the Fock state basis vector can be connected with the spin state by a transformation matrix .
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, means that 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. This means the averaged encoding complexity of the fermionic creation and annihilation operators under the Bravyi–Kitaev transformation is .
Data Availability
The data sets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 12005015 and 12004206).