Molecular Structure Optimization Based on Electrons-Nuclei Quantum Dynamics Computation.

A new concept of the molecular structure optimization method based on quantum dynamics computations is presented. Nuclei are treated as quantum mechanical particles, as are electrons, and the many-body wave function of the system is optimized by the imaginary time evolution method. The numerical demonstrations with a two-dimensional H2 + system and a H-C-N system exemplify two possible advantages of our proposed method: (1) the optimized nuclear positions can be specified with a small number of observations (quantum measurements) and (2) the global minimum structure of nuclei can be obtained without starting from any sophisticated initial structure and getting stuck in the local minima. This method is considered to be suitable for quantum computers, the development of which will realize its application as a powerful method.

Introduction

Accurate quantum chemistry computations remain a challenge on classical computers, especially for molecules with industrially relevant sizes, despite the significant efforts and developments of recent years. The computational cost of exact methods for quantum chemistry on classical computers grows exponentially with the molecular size,[1,2] whereas the cost can be suppressed in polynomial scaling on quantum computers.[3] For this reason, quantum chemistry computations have been considered to be a promising application of quantum computers. By the manipulation of quantum states of matter and with the advantage of their unique features, such as superposition and entanglement, quantum computers promise to efficiently deliver accurate results for many important problems in quantum chemistry, such as the electronic structure of molecules.[4,5]

In quantum chemistry computations, the Born–Oppenheimer (BO) approximation[6] has typically been used to speed up the computation of molecular wave functions and other properties for large molecules. The BO approximation is the assumption that the wave functions of atomic nuclei and electrons in molecular systems can be treated separately, on the basis of the fact that the nuclei are much heavier than the electrons. In many cases of quantum chemistry computations by classical computers, nuclei have been treated as point charges (classical particles). Many quantum algorithms for quantum chemistry computations on quantum computers are also typically based on the BO approximation: for example, the quantum phase estimation (QPE) method[7,8] for a fault-tolerant quantum computer (FTQC) and the variational quantum eigensolver (VQE) method[9,10] for a noisy intermediate-scale quantum (NISQ)[11] device are often employed to solve systems under the BO approximation. On the other hand, it is also possible to think of methods beyond the BO approximation with quantum computers, and such methods have been proposed in the literature.[4,12−17] Kassal et al.,[13] for instance, reported that a completely nonadiabatic grid-based method on FTQCs, where nuclei are treated as quantum particles and correlated with electrons, is not only more accurate but also faster and more efficient than the methods based on the BO approximation. Moreover, the molecular orbital theory beyond the Born–Oppenheimer approximation[18] has been extended to the NISQ algorithms such as VQE.[17]

Optimization of the geometry of molecules is an important process to obtain the equilibrium molecular structures in quantum chemistry computations.[19] Since the physical and chemical properties of molecules are dependent on their specific geometrical structures, elucidation of the optimized structures enables the prediction of properties and identification of chemical products. In the conventional molecular structure optimization methods, which are based on the BO approximation and treat nuclei as classical particles, the electronic states should be computed to obtain the forces acting on nuclei and update the molecular structure for each optimization step. These conventional methods can be accelerated through the use of quantum computers to perform the electronic state calculations.

However, a large number of circuit executions and observations (measurements for qubits) are required to obtain expectation values of observables (the forces acting on nuclei in this case) because we collapse the state of qubits and destroy a great deal of the information when we measure the qubits. Moreover, it is necessary to repeat the above repetition (force computations) at each iteration in the geometry optimization process. Furthermore, as gradient-driven descent routines are often used to optimize the molecular structure, the system tends to be relaxed to the closest local minimum from an initial structure for the geometry optimization. It is thus a difficult task to find the global minimum structure in systems that have multiple local minima (i.e., many isomers), such as alloy cluster systems.[19]

In this study, we propose a method for optimizing molecular structures based on quantum dynamics computations with working on an FTQC in mind. In our method, the many-body wave functions of nuclei and electrons are directly treated as wavepackets and optimized by the imaginary time evolution method. The many-body wave function on the lattice can be mapped to a state on a quantum computer with the number of qubits growing only logarithmically with respect to the number of lattice points.[5,12] Although the imaginary time evolution is a nonunitary operation that is not straightforward to implement on a quantum computer, which is based on the unitary evolution of quantum systems, probabilistic methods[20−22] and measurement-based methods[23] for FTQC have been proposed. The method of linear combination of unitaries could also be useful.[24,25] For NISQ devices, a variational method has also been proposed,[26] as we will review later. Since nuclei are treated as quantum particles, the nuclear quantum effects are naturally included in our method. More importantly, our method is expected to have the following favorable characteristics. First, it is theoretically expected that the global minimum of the molecular structure can be obtained just by starting with an initial many-body wave function that covers all molecular structures. This is because the imaginary time evolution is not a gradient-based optimization and in principle will not suffer from being trapped in the local minima of the molecular structures. Second, the optimized many-body wave functions have large stochastic amplitudes at the most stable structure of nuclei, so that the determination of the optimized nuclei positions, i.e., the molecular structure optimization, can be performed with a small number of observations (quantum measurements). We demonstrate our concept of the molecular structure optimization method based on quantum dynamics computations for a two-dimensional H2+ system and a H–C–N system with numerical simulations by a classical computer.

Methods

In this section, we explain our proposal for the molecular structure optimization.

Quantum Imaginary Time Evolution

In nonrelativistic quantum dynamics computations, the time evolution of quantum systems can be described by the time-dependent Schrödinger equationwhere H is the Hamiltonian that describes the system. The solution of the time-dependent Schrödinger equation is formally written by

That is, if an appropriate initial wave function ψ(0) can be prepared, then the time evolution of the system can be obtained by application of the time evolution operator

If ψ(0) is the eigenstate of H, i.e. ψ(0) = ϕ and Hϕ = Eϕ, then eq becomes

If not (in general cases), then the exponential operator must be applied to the wave function to compute the time evolution of the system. A number of methods have been developed to apply the exponential operator to the wave function.[27,28] Here we explain the method based on the second-order Suzuki–Trotter decomposition.[4,5,29] It should be noted that the first-order decomposition is also used frequently in quantum computing, though it is subject to errors due to the noncommutativity of the terms in the Hamiltonian. The second-order decomposition compensates for noncommutativity and eliminates some error terms. One might think that it would be better to use a third- or higher-order decomposition in order to increase accuracy, but this is not practical because the number of terms for the decomposition increases exponentially with the order.[30] We have to minimize the number of quantum gates to reduce errors in exectuting quantum gates to obtain meaningful results in quantum computing, especially when NISQ devices are used. Therefore, lower-order decompositions are preferred. In this study, the second-order decomposition is used because of its balance between accuracy and number of terms.

Denoting the kinetic energy term of H by T and the potential energy term by V, the second-order Suzuki–Trotter decomposition is

The error arises because T and V are noncommutative. To reduce the error, we set t = Nstepdt and express the time evolution operator as the product of Nstep operators in time increments of dt

Because T is diagonal in wavenumber space, when we apply in the second-order Suzuki–Trotter decomposition, the wave function is expressed in wavenumber space by performing a fast Fourier transformation (FFT) for the wave function in real space, when the calculation is run on classical computers. Similarly, V is diagonal in real space; thus, the operator is applied to the wave function in real space after application of the inverse FFT to the wave function in wavenumber space, . Note that if a matrix is diagonalthen its exponential can be obtained by exponentiating each entry on the main diagonal

Using these procedures, the application of the time evolution operator U(dt) can be computed without expanding the exponential function operator. When these calculations are repeated Nstep times, the wave function at the desired time t can be obtained.

Now let us convert time into an imaginary number, it → τ. The solution of the time-dependent Schrödinger equation then becomesWhen we formally expand ψ(0) by the eigenfunctions of H, {ϕ}, we havewhere is a coefficient for the expansion ψ(0) = ∑cϕ and E is the eigenvalue of ϕ. If we increase τ, the factor attenuates more quickly for large E, and the ground state ϕ0, which has the smallest eigenvalue E0, remains until the end as long as the coefficient (overlap between the initial state ψ(0) and the ground state) c0 is nonzero. This is how we obtain the ground state by an imaginary time evolution. Note that the norm of the wave function decreases during the imaginary time evolution; thus, it is necessary to renormalize the wave function after the evolution.

Geometry Optimization with Fully Quantum Wave Function for Nuclei and Electrons

In the previous section, we showed that the most stable quantum state of the system can be obtained by the imaginary time evolution method. In this section, we show that this property can be used to obtain the most stable structure of a molecule where the atomic nuclei are regarded as quantum mechanical particles as well as the electrons, in contrast to the conventional method based on the BO approximation. We consider the many-body wave function ψ(R, r, τ = 0), including the degrees of freedom of nuclei and electrons, where R represents the coordinates of Nnuc nuclei, , and r represents the coordinates of n electrons, r = r1, r2, ..., r. From any initial wave function as the starting point (ultimately, it may be a constant over the entire simulation region in real space but should have nonzero overlap with the desired ground state), the most stable state (ground state) of the nuclei and the electrons can be obtained by the imaginary time evolution of the initial wave function under the Hamiltonianwhere Q is a electric charge of the ith nucleus and qe (>0) is a charge of an electron. The optimized wave function is expected to have a large norm at the state of the most stable molecular structure; thus, the stable structure can be obtained with a high probability by performing the measurement to the state after imaginary time evolution. This calculation includes quantum effects such as the zero-point oscillation effect and the nonadiabatic effects that are ignored under the Born–Oppenheimer approximation.[31] If we take an initial state (for the nuclei part of the wave function) as the superposition state for all possible molecular structures, the most stable structure among all isomers of the molecule would be obtained. Because the search for the most stable isomer generally becomes difficult as the number of atoms increases,[19] it is an appealing feature of this method that treats the atomic nuclei as quantum mechanical particles.

Such an approach could not be taken to date due to the difficulty of describing many-body wave functions with a classical computer. In contrast, a many-body wave function can be described with a reasonable number of qubits on a quantum computer. This is discussed in more detail later.

Computational Setups for Numerical Demonstrations of the Proof of Concept

We present numerical demonstrations of our method by taking a two-dimensional (2D) H2+ molecule and a H–C–N system as examples to show that the molecular structure optimization is possible with a small number of observations and illustrate the feature that the most stable molecular structure can be found without being trapped in local minima.

Computational Setup for the 2D H2+ Molecule

A H2+ molecule is a three-particle system with two protons and an electron that is naively described as a six-dimensional system, which will be difficult to calculate with a classical computer. For example, if we take a grid of 25 = 32 per dimension in that space, the wave function has 230 = 1073741824 elements and the Hamiltonian has 260 = 1152921504606846976 elements in the naive and worst estimate. On the other hand, in a quantum computer, a wave function can be expressed with 30 qubits because N qubits can represent 2 variables: i.e., 2 wave function elements on grids. In addition, the number of qubits scales linearly with the number of particles.

In numerical demonstrations, we perform a dimensional reduction to the system so that the calculation can be performed on a classical computer. First, the center of mass of the system is fixed at the origin (0, 0). As a result, the degrees of freedom of translational motion of the molecule can be ignored. Now that the purpose is to optimize the molecular structure, we only need to focus on the internal degrees of freedom. Classically, the rotational motion around the center of mass has also been frozen, but quantum mechanically, the rotational motion and the vibration mode are coupled, so that the rotational motion cannot be separated. Even if the rotational motion is left, observing the optimized wave function gives only one molecular structure in a specific orientation of the rotational degrees of freedom. When we obtain multiple results for the measurements of the positions of the nuclei, the obtained structures are described in terms of internal coordinates, which allows a rotational-free description of molecules. The mean of each nucleus position is then calculated, which gives the molecular structure. Next, the relative coordinates R⃗ between the two nuclei are introduced to further reduce the dimension. The relationship between the reduced coordinates and the original coordinates is given aswhere , , and mp are the coordinates and mass of the protons, respectively, and and me are the coordinates and mass of the electron, respectively. are the coordinates of the electron as seen from the center of mass of the nuclei, and is the center of mass of the system, which is fixed at the origin in the simulation. We also define the total mass of the system M = 2mp + me and that of the nuclei Mn = 2mp. The reduced coordinates have four degrees of freedom, R⃗ = (R, R) and . The reduced masses for R⃗ and are μR–1=mp–1 + mp–1 and , respectively. By using the reduced coordinates, the Hamiltonian of the system becomes

Here the atomic unit is introduced for simplicity (ℏ = 1, Qp = 1, qe = 1). We take an evenly spaced isotropic grid of 25 elements per dimension in each two-dimensional spaces for . See the Supporting Information for more details on the computational setup.

Computational Setup for the H–C–N System

Since there are no isomers in the H2+ system, we cannot illustrate one of the promising features of our method that the most stable molecular structure can be found without being trapped in local minima. Therefore, we present a numerical demonstration of our method in a H–C–N system, which has two stable linear isomers: i.e., HCN and HNC. These isomers have been intensely studied[32−35] because this system is one of the simplest systems in which isomers exist and is important as a tracer of moderate gas densities in molecular clouds and galaxies.[36] From these studies, HNC is known to be unstable by about 0.6 eV over HCN.[35] In principle, we can treat the many-body electron–nucleus Hamiltonian and wave function for this system as we conducted for the H2+ system, though antisymmetricity of electrons must be taken into account in multielectron systems. How to impose antisymmetry will be discussed later. Due to the limitations of classical computers, however, it is difficult to simulate the quantum dynamics for systems with large degrees of freedom. To save computational resources, one can use ab initio density functional theory for the electronic degrees of freedom or empirical potentials to completely omit the calculation of electronic degrees of freedom. Because the purpose of the demonstration here is to show the behavior of the method in a system with isomers, we employ the Lennard–Jones potential for nuclei to avoid treating electronic degrees of freedom and consider only nuclear degrees of freedom. Moreover, we treat the problem as a one-dimensional problem along the molecular axis since the two isomers, HCN and HNC, have linear geometries. The Hamiltonian of this system is then expressed aswhere MI and RI are the mass and the coordinate of atom I, respectively, and VIJLJ is the Lennard–Jones potential for the atomic pair I and J. The parameters of the Lennard–Jones potential were determined on the basis of the two-body binding energy and bond length between each atom. With the Lennard–Jones potential, HCN was found to be the most stable structure and HNC was 0.184 eV higher in energy. We take an evenly spaced isotropic grid of 27 elements for each one-dimensional space of RH, RC, and RN. Details of the computational setup, including the parameters of the Lennard–Jones potential, are given in the Supporting Information.

Results and Discussion

2D H2+ Molecule

We present results of numerical demonstrations for the method using the 2D H2+ molecule. The obtained probability density distribution for all coordinates is multidimensional, making it difficult to display on paper. We instead plot the conditional probability density distribution of the relative distance R⃗ between the protons when the electron is near the origin (0.0375, 0.0375), i.e., |ψ(R, R, rc = 0.0375, rc = 0.0375)|2, in Figure (note that the grid used in this study does not contain the origin (0, 0)).

Figure 1

Conditional probability density distribution of the relative coordinates of protons R⃗ for the H2+ molecule when an electron is near the origin point (0.0375, 0.0375).

We find that the probability density distribution of the relative coordinates R⃗ between the protons has a donut-shaped structure. This means that the H2+ molecule can be oriented in any direction because it has rotational degrees of freedom. This probability density distribution must be donut-shaped and uniform for all directions; however, the probability is large and small in specific directions in Figure . The reason for this is considered to be the following. The roughness of the grid space used in this study cannot represent the exact internuclear distance for the equilibrium structure in specific directions. Therefore, the probability density distributions for such directions would be small and the probability density distributions for the specific directions where the grid is well-matched with the equilibrium internuclear distance would be large. A better shaped “donut” will be obtained by using a polar coordinate grid. It should be noted that it is necessary to repeat the calculations and observations many times to obtain such a probability density distribution with a quantum computer since the observation makes the wave function (superposed state) collapses to a single selected component. However, it is not necessary to get the probability density distribution to obtain the optimized molecular structure because we are interested only in the states with high probabilities corresponding to the stable structures.

We also calculate the expectation value (mean) of the proton–proton distance |R⃗|, defined asfor our simulation data in Figure . The result is 0.41 bohr, which is slightly larger than the value of 0.37 bohr, where the protons were treated as point charges (see the Supporting Information). This larger expectation value is considered to be due to the spread of the nuclear wave function and its binding potential, which is very tight for smaller |R⃗| but loose for larger |R⃗|. It should be noted that these bond distances differ significantly from those of an ordinary three-dimensional H2+ molecule. This is because not only the atomic configuration but also the electronic degrees of freedom were treated as a two-dimensional system.

Next, we estimate how many measurements for the wave function would be required to optimize the molecular structure within an acceptable error. gives the simultaneous probability of finding relative proton-to-proton coordinates in R⃗ and electrons in . The grid space for (220 points) is represented by 20 qubits using the binary encoding (0···000, 0···010, ..., 1···111). Each binary combination corresponds to the specific coordinates of R⃗ and . If we perform the projective measurement on these qubits, one of the state (the coordinates R⃗, ) is chosen from (0···000, 0···010, ..., 1···111) according to the probability . We conduct the measurement (observation) simulations by randomly sampling the qubit states with the pseudorandom numbers generated by the Mersenne Twister method.[37]

Figure shows how the mean of |R⃗| for Nobs observations converges as we increase Nobs. The exact mean and the standard deviation (corresponding to the zero-point vibration of the ground state) calculated from p(R⃗) are also shown by dotted lines. Note that the optimized |ψ(R)|2 (even for τ → ∞) has a distribution with nonzero width by zero-point vibration. We see that the mean of the sampled |R⃗| gets close to the exact mean with sufficient accuracy after as small as 200 observations. The mean |R⃗| obtained from 200 observations was 0.4082 bohr, while the exact value is 0.41 bohr as explained earlier. The standard deviation for 200 observations was 0.0611 bohr, and the mode (the value that was obtained most frequently) was 0.4142 bohr. The reason why the converged bond length R⃗ was obtained with such a small number of observations is that the wave function has a peak (cusp) at the state corresponding to the most stable structure and that the measurements will pick it up with high probability. It should be noted that, as the size of the molecule increases, the number of vibrational modes also increases. As a result, the fluctuation width of the zero-point vibration also increases, which may affect the required number of measurements to determine the molecular structure with a certain precision. As discussed in the Supporting Information, the standard deviation of the zero-point fluctuation scales as , where N is the number of nuclei. The magnitude of the error after M measurements under the fluctuation is given by . Therefore, the required number of measurements for a certain precision will scale as O(N).

Figure 2

Mean of the sampled bond distance |R⃗| drawn from the optimized wave function. We plot the number of observations Nobs versus the mean up to Nobs observations, , where is an ith sample of the observed R⃗. The red (blue) dotted line represents the exact mean (standard deviation) of |R⃗|.

Finally, we illustrate how the presented result above changes when we consider heavier nuclei. It is expected that, as a nucleus becomes heavier, the fewer wavelike properties it has and the closer it is to a classical particle. It is predicted that the peaks at the nuclear positions of the wave function will be even sharper and that converged structures of molecules will be obtained with a smaller number of observations. To see this, we consider the deuterium molecular ion D2+, in which the protons in the hydrogen molecular ion are replaced with deuterons (nucleus of deuterium; the mass is 2 times that of a proton), and the tritium molecular ion T2+, in which tritons (nucleus of tritium; the mass is 3 times that of a proton) are substituted. We perform the same calculations as for H2+ for these isotopes. Figure shows the conditional probability density distributions of the relative coordinates R⃗ for deuterons and tritons, respectively, when the electron is near the origin (0.0375, 0.0375).

Figure 3

Conditional probability density distributions of the relative coordinates of nuclei for (a) D2+ and (b) T2+ molecules when the electron is near the origin point (0.0375, 0.0375).

The absolute values of the peaks are higher than those for H2+, and the spreads of the wave functions are narrowed. This tendency was more significant for the heavier T2+. In addition, the peak bias increases for specific directions with the heavier isotopes D2+ and T2+ because, as the wave functions become narrower, it becomes difficult for the grid points to exist at the positions where the peaks are in the continuous limit (infinite number of meshes of the gird) because the grid size is not fine enough. Subsequently, as in the case of H2+, we simulate how many observations were required to optimize the molecular structure of D2+ and T2+ within an acceptable error. Figure plots the same quantity as Figure , the mean of the sampled |R⃗| from the optimized wave function, for D2+ and T2+ molecules. The exact mean and the standard deviation are also shown as dotted lines.

Figure 4

Mean of the sampled bond distance |R⃗| drawn from the optimized wave function for (a) D2+ and (b) T2+ molecules. We plot the number of observations Nobs versus the mean up to Nobs observations, , where R⃗ is an ith sample of the observed R⃗. The red (blue) dotted line represents the exact mean (standard deviation) of |R⃗|.

In addition, the means, standard deviations, and modes obtained from 200 observations are summarized in Table along with the result for H2+ for comparison.

Table 1

Mean, Standard Deviation, and Mode of the 200 Observations for |R⃗| in bohr

system	mean	standard deviation	mode
H2+	0.4082	0.0611	0.4142
D2+	0.3968	0.0494	0.3712
T2+	0.3778	0.0355	0.3712

These results in Figure and Table indicate that the variation of the observed R⃗ becomes smaller and the convergence becomes faster as the mass increases. This is because of the smaller fluctuation widths of the zero-point vibrations for heavier nuclei. Since all the nuclei of atoms after the He atom in the periodic table are heavier than a triton, it is expected that molecular structures obtained by measurements will converge even faster in terms of the number of measurements. We briefly comment on the bond length |R⃗| before ending the explanation of numerical demonstrations. The mean of the bond length approaches the value for point charges (0.37 bohr) as the mass increases. The tendency for the bond length to become shorter as the mass increases has also been confirmed for the 3D hydrogen molecule and its isotope molecules.[38] The modes of the sampled |R⃗| are the same for D2+ and T2+ because the values between 0.3712 and 0.4142 cannot be taken due to the rough grid used in our calculations. This result indicates that the step size of the grid must be sufficiently smaller than the desired resolution of the molecular structures. For calculations on a classical computer, making the grid finer leads to a exponential increase in computational cost, but calculations on a quantum computer will scale polynomially.[13]

H–C–N System

The time evolution of the probability density distribution during the imaginary time evolution for the H–C–N system is presented in Figure , where the conditional probability density distributions when the C atom is located at the center of the coordinate (RC = 7 bohr) are shown. Figure indicates that four peaks appear in the conditional probability density distribution immediately after the start of the imaginary time evolution. Two of the four peaks represent the HCN molecule and its oppositely oriented molecule (NCH). The remaining two peaks represent the HNC molecule and its oppositely oriented molecule (CNH). As the imaginary time evolution proceeds further, the peaks corresponding to the HNC molecule decay and only the peaks corresponding to HCN molecule, the most stable state, remain. However, since the two isomers are energetically close, there are still small peaks corresponding to the HNC molecule even after 2000 time steps.

Figure 5

Imaginary time evolution of the conditional probability density distributions of the H and N atomic coordinates for the HCN molecule when the C atom is located at the center of the coordinate (RC = 7 bohr).

Next, a numerical simulation of the measurement (observation) of the nuclear positions was conducted for the H–C–N system as we did with H2+ molecule. The histogram of the coordinates of each atom obtained from 100 observations is shown in Figure a.

Figure 6

Histogram of observed H–C, C–N, and H–N distances for the HCN molecule when (a) the actual masses are used and (b) the artificial masses (M = 184000 au) are used as a classical limit.

Note that the optimized |ψ(R)|2 has a distribution with nonzero width by zero-point vibration. Figure shows that the most stable state HCN (RHC ≈ 2.2 bohr, RCN ≈ 2.2 bohr, and RHN ≈ 4.4 bohr) has been the most frequently observed. However, the HNC isomer (RHN ≈ 2 bohr, RCN ≈ 2.2 bohr, and RHC ≈ 4.2 bohr) has also been observed, albeit very infrequently. The reason why RCN has only one peak is because the distance between C and N does not change between HCN and HNC. These results indicate that, by deliberately choosing the number of steps in the imaginary time evolution, the structures of metastable isomers can also be determined.

Finally, we present the same result obtained by using the artificial masses for the nuclei (the masses of all atoms are set to 184000 au) in Figure b to see the classical limit. Only one peak corresponding to HCN can be seen in Figure b. The possible reasons for this are as follows. When actual (natural) masses were used, the probability density distribution peaked at two states due to the effect of the zero-point vibrational energy, which is greater than the energy difference between the two isomers. However, when the masses are increased, the contribution of the zero-point vibration energy disappears, and only the peak of the most stable state remains. We also conducted a simulation without the kinetic energy terms of the atoms (corresponding to M → ∞), and got the same result as in Figure b. What these results tell us is that we can ignore the kinetic energy term of the nuclei and reduce the quantum gates for that if we only want to know the most stable structure in the classical limit, ignoring the quantum effects of the nuclei. This is valid not only for our demonstration where the wave function of electrons is not directly treated but also for general cases where it is explicitly considered. We note that the electron kinetic energy term cannot be neglected even in those cases.

Implementation on Quantum Devices

In previous sections, we have demonstrated the concept of molecular structure optimization using the imaginary time evolution method of an electron–nucleus wave function. However, the method is difficult to apply to general systems, except for very small systems as used in the demonstrations. This is because of the difficulty of describing high-dimensional wave functions with a classical computer. In contrast, high-dimensional wave functions can be described with a reasonable number of qubits on a quantum computer.

In this section, we provide a brief idea on how to implement the imaginary time evolution of wave functions on NISQ devices. We consider a single-particle Hamiltonian in a one-dimensional systemwhere we expand the Hamiltonian on the position basis. H is mapped to a qubit model in the following way: we first discretize the position operator and the squared momentum operator as and , respectively, with a being the size of each mesh. The latter replacement is derived from the formula . These replacements lead to the discretized modelwhere t = ℏ2/(2ma). The discretized position basis |k⟩ is expressed in the computational basis, which yields the qubit model of Hlat. Any one-particle state ψ(x) is approximated aswhich can be prepared or realized on quantum computers. The energy expectation value ⟨ψ|Hlat|ψ⟩ can be efficiently computed with a circuit of polynomial size; it can be shown that the increment operator ∑|k⟩⟨k + 1| and its conjugate are constructed by Toffoli and CNOT gates (see e.g. refs (39 and 40)), and also the diagonal part can be expanded as a sum of unitary operations,[41,42] so that the expectation value for each term can be evaluated using the Hadamard test, whose output can be measured with a low-depth circuit.[43]

As was explained in the Introduction, one can use various methods to perform the imaginary time evolution. To make our proposal more concrete, we shall review one such method, the variatonal quantum simulation (VQS).[44] VQS is based on the McLachlan variational principle,[45] which seeks the best parameter θ(t) of a given ansatz at each t, to describe the time evolution. The time derivative is approximated by , and is obtained by minimizing the cost function . The imaginary time evolution is also simulated in a very similar manner with the cost function .[26] should be introduced so as to preserve the norm of the state.

Moreover, we explain one technical aspect of our method: the symmetry of many-body wave function. For systems with identical particles, the wave function should be (anti)symmetrical with respect to a particle exchange. In such cases, one can (anti)symmetrize the wave function by adding a penalty term in the Hamiltonian.[46] For example, let us consider a system with two fermions which we call Htwo. We also denote the discretized position the basis of two-particle states by |k1,k2⟩ ≡ |k1⟩⊗|k2⟩. A projector which projects out antisymmetric states is written as follows:

Note that a corresponding projector for boson is (1 – Oswap)/2. By adding gP with a sufficiently large coefficient g, the ground state of Htwo + gP coincides with the lowest energy state in antisymmetric wave functions of Htwo. Since the (unitary) operator Oswap is expressed as a sequence of SWAP gates between qubits for the positions of two particles, the expectation value of Oswap can be evaluated by a Hadamard-test type circuit of those gates. The extension to a multiparticle system is straightforward by adding more penalty terms corresponding to an exchange of all pairs of particles, which amounts to M(M – 1)/2 terms, with M being the number of identical particles. We also note that another quantum algorithm for preparing antisymmetrical wave functions has been proposed starting with a Hartree product of molecular orbitals.[47]

Conclusions

The concept of a molecular structure optimization method using quantum dynamics computation has been presented. The nuclei have been treated as quantum mechanical particles as were the electrons, and the many-body wave function of the system was optimized using the imaginary time evolution method. The optimized wave function has a large probability amplitude at the most stable structure of the nuclei, which allows us to determine the optimized nuclear positions with a small number of observations (quantum measurements). Our method has a favorable feature that the most stable isomer structure can be obtained for complex systems without being trapped in local minima by virtue of the imaginary time evolution. This method may be particularly promising in the search for the most stable structure for systems with many isomers, such as metal alloy clusters. Another aspect of our method is that it includes nuclear quantum effects (zero-point oscillation, nonadiabatic correction, etc.) that are typically ignored by conventional methods. The Coulomb interactions between multiple particles are explicitly incorporated, and there is no approximation in nucleus–nucleus, electron–nucleus, and electron–electron interactions. Although executing our method for industrially interesting large molecules with classical computers is difficult because of huge computational costs in treating the fully quantum wave function for nuclei and electrons, quantum computers (possibly FTQC) can be appealing candidates to run our method. Our proposal can give a new insight into quantum chemistry computations on quantum computers.