Molecular dynamics on quantum annealers.

In this work we demonstrate a practical prospect of using quantum annealers for simulation of molecular dynamics. A methodology developed for this goal, dubbed Quantum Differential Equations (QDE), is applied to propagate classical trajectories for the vibration of the hydrogen molecule in several regimes: nearly harmonic, highly anharmonic, and dissociative motion. The results obtained using the D-Wave 2000Q quantum annealer are all consistent and quickly converge to the analytical reference solution. Several alternative strategies for such calculations are explored and it was found that the most accurate results and the best efficiency are obtained by combining the quantum annealer with classical post-processing (greedy algorithm). Importantly, the QDE framework developed here is entirely general and can be applied to solve any system of first-order ordinary nonlinear differential equations using a quantum annealer.

Introduction

Full-fledged quantum computers are expected to offer significant computational advantages compared to the traditional (classical) computers but, realistically, this is a distant future[1,2]. Present-day quantum computers have relatively modest computational capabilities and are not yet free of technical issues such as noise and decoherence that cause errors[3]. Still, these early machines can be used as a testbed for the development of quantum algorithms and for the small-scale proof-of-principle simulations, paving the way to practical applications on the next generation hardware. One of the long-thought applications of quantum computers is to simulate quantum systems, such as molecules and materials at the atomistic scale[4-11]. Several algorithms were recently proposed for the calculations of molecular electronic structure (quantum chemistry) using either gate-based quantum computers with only a few tens of qubits, such as Sycamore[12,13], IBM Q-machines[14-16], or using quantum annealers with a larger number of qubits, such as the D-Wave systems[17-21].

Another essential component of chemical modeling is the propagation of the equations of motion for constituent atoms in time and space (e.g., along the reaction path), called molecular dynamics simulations. For over 50 years, molecular dynamics simulations have played a pivotal role not only in guiding and analyzing experiments, but also as an interdisciplinary computational tool that is able to reach far beyond experimental conditions. Thus, molecular dynamics simulations are an indispensable tool for a wide variety of research topics and several attempts to harness the power of quantum computing for its benefit have already been made, including a method for solution of the quantum vibrational eigenvalue problem on a quantum annealer[22,23], and the trajectory simulation of molecular vibrations on the IBM quantum devices[24,25]. These pioneering calculations were restricted to small molecules, such as diatomic and triatomic systems, and were often done in a hybrid quantum/classical fashion, when only a part of the overall algorithm is run on a quantum device. For example, in the trajectory simulations of Ref.[24], a quantum computer was used only to obtain the gradient of the molecular potential energy (the force determined by electrons), whereas the motion of the atoms (the actual trajectory) along this potential energy surface was propagated on a classical computer.

The goal of the present paper is to demonstrate that quantum architectures, in our case the D-Wave quantum annealer, can be used to simulate the molecular dynamics component pertaining to any computational chemistry problem, with possible future applications in chemical dynamics, material science or drug design. For this purpose, we developed a new method named Quantum Differential Equations (QDE) and carried out the first ever calculations of this sort on an actual quantum annealer (D-Wave 2000Q) to propagate trajectories for the motion of atoms in the hydrogen molecule, H2 (i.e., vibrations). Although this is the simplest diatomic molecule, our method and code are general and suitable for simulation of polyatomic molecules with multiple vibrational modes, which will be accessible on the next generation of D-Wave annealers, such as Advantage[26]. Moreover, our algorithm paves the way to the full-quantum molecular dynamics simulations within the time-dependent framework, including such methods as quantum trajectory, path-integral or thawed Gaussian[27].

Theory

In general, the goal of molecular dynamics is to determine how a system of interacting atoms evolves in time. Here we will focus on the vibration of one molecular bond, when the trajectory is obtained by solving the following system of Hamilton’s equations: where is time, is the bond length, is its associated momentum, is a reduced mass, and is the force field acting on atoms due to the local gradient of potential energy. In polyatomic molecules, equations similar to Eq. (1) can be written for every degree of freedom, coupled through the overall force field . Therefore, a general analogue of Eq. (1) can be written as the following system of differential equations: where is a vector of unknown real functions , and is a set of known arbitrary real functions. Without loss of generality, we can consider a system of 1st-order equations, since higher-order equations can always be rewritten as a system of 1st-order equations. To represent the functions and numerically, we introduce an equidistant grid of points with step size over the range of problem-specific values of the variable . The functions and are then represented as arrays of their values at the grid points: and . Our goal is to formulate this problem as a binary minimization problem, which will enable its solution on a quantum annealer[28]. While other methods to solve systems of differential equations on quantum annealers have also been considered in the literature[29,30], we believe our approach is more general.

The only problem that a quantum annealer can handle is the minimization of a user-defined functional (subject to certain limitations, as discussed below). Therefore, to make system (2) solvable on a quantum annealer, we need to define a functional such that the correct solution vector minimizes its value. One way to do this is to introduce a cost function defined as the total squared difference between left- and right-hand sides of all equations in the system, at all points of the grid, namely:

The value of the derivative can be approximated by the first order finite difference scheme:

Substituting this into Eq. (3) we obtain:

Note, that Eq. (5) treats all values of on an equal footing, which may not be desirable in practice, for example, if the relevant ranges of their values are substantially different. In this case, data rescaling or introduction of a penalty factor might be helpful to achieve better results.

The first practical consideration for Eq. (5) is that the variation of the solution at all points of the grid at once (one run with maximum value of ) may not be computationally feasible using present-day quantum annealers. Instead, one can consider only a subset of points in a given run of the annealer, and then use the final point obtained in this run as an initial condition for the next subset of points (smaller but multiple runs). In the limiting case when only one point of the grid is considered at a time (i.e. , maximum number of consecutive runs), the values of are all known from the previous step, which decouples all Eq. (5) and makes it possible to solve them one by one. This procedure permits us to solve, piece by piece, even large systems, even on annealers with a small number of qubits. As more qubits become available in the future, they can be effectively utilized to consider more grid points at the same time, and thus reduce the total number of runs.

Another practical problem is that is an arbitrary function, which makes Eq. (5) non-quadratic in general, whereas the current generation of quantum annealers does not have native support for non-quadratic functions. To keep Eq. (5) quadratic in , the has to be no more than linear in . One way to solve this issue is to split into linear segments and find the solution for each segment separately, using the last solution point in each segment as an initial condition for the next one. This separation into multiple runs is similar to what has been discussed in the previous paragraph and both of these procedures can be combined. Namely, when solving for only a few grid points per run (so that the solution does not advance too far), the local behavior of around the initial point of this run can be approximated with a linear segment reasonably well. In the case when only one point is considered at a time, the values of are constant, so the issue of separating into linear segments does not arise at all. Also, in special cases (e.g. harmonic oscillator) the forces have naturally linear dependence, and no separation is necessary as well.

Now let us obtain an expression for the functional at each segment. The right-hand side function is linear with respect to within a given segment, so it can be written as: where the index labels individual terms in the linear expansion of . Plugging this into Eq. (5) and expanding squared terms, one obtains the following expression (quadratic in ):

From this expanded form, it is easy to see the individual coefficients for all unknown variables and to compose the so-called Quadratic Programming (QP)[31] matrices and for this function, such that Eq. (7) can be written in matrix form as:

The last step that we need to do before we can submit this task to a quantum annealer is to convert the QP-matrices into the corresponding QUBO matrix (for Quadratic Unconstrained Binary Optimization)[28], which has the same meaning as QP-matrices, but, instead of the continuous variables of Eqs. (7) and (8), it uses binary variables that represent the states of the qubits after the read operation. This conversion can be done as follows. First, a given continuous variable can be approximated by a set of binary variables via a regular signed fixed-point number representation[22]: where and are discretization parameters, equal to the number of qubits used to represent integer and decimal parts of , respectively. Equation (9) can approximate numbers in the range from to with an error of up to . The floating-point representation, while being superior to the fixed-point representation on classical computers, is not used here because this representation is not quadratic in , which makes it unsuitable for use on quantum annealers, as discussed above.

Next, to express the functional in terms of the binary variables , one could plug Eq. (9) into Eq. (7), and compose the QUBO-matrix in the same way the QP-matrices and were composed. But this is tedious, and the resulting expression is rather long. Therefore, we will not derive it explicitly here. Instead, the QUBO-matrix can also be obtained directly from the QP-matrices by simply replacing each element of the matrix with a block (), where the value of element is given by: Here and are indices of rows and columns of the original QP-matrices and , while and are local indices of rows and columns within each block. The range of in each block is the same as in Eq. (9): from to . Note that since there is no need for a separate vector , as the relevant contributions are simply added to the diagonal of . Thus, the minimized functional can be written in QUBO formalism as:

The QUBO-matrix defined by Eq. (10) can be used as a direct input to a quantum annealer to obtain the values of that minimize Eq. (11), which, in turn, can be used to reconstruct the continuous solution to the problem of Eq. (2).

A schematic representation of the workflow of QDE is shown in Fig. 1. First, the potential energy surface is used to calculate the forces and compose a QUBO matrix for the next trajectory points. The QUBO matrix is then given to D-Wave, which finds the solution and updates the trajectory. Finally, the forces at the new trajectory points are used to compose a new QUBO matrix and the cycle repeats until the trajectory propagates for the specified time.

Figure 1

Schematic representation of the workflow of QDE.

Results

Practical implementations of this method were carried out using the D-Wave 2000Q quantum annealer at Los Alamos National Laboratory. The values of molecular parameters were chosen to represent H2 with an analytical Morse potential as described in the Supplemetary Information (in principle, an arbitrary potential energy function could be used). Discretization parameters and were set to 6 and 15, respectively. This choice can represent the numbers from – 32 to + 32 with rounding errors up to (in atomic units for both coordinate and momentum), using qubits per number.

The D-Wave 2000Q has 2048 physical qubits, but it is important to emphasize that these qubits are different from those used in gate-based quantum computers in the sense that there is only one operation that can be performed on them (minimization), and each qubit can only interact with 6 neighbors. In cases when full connectivity is required (as in this work), a number of physical qubits have to be “merged” into a single fully connected logical qubit representing a single variable. The maximum number of such fully connected logical qubits on the D-Wave 2000Q is only 64[26], which is nevertheless sufficient to either propagate through one time step using both equations of the system of Eq. (1) at once, or to propagate the equations for position and momentum individually one after another. We tested both strategies. All trajectories in this work were propagated through the time interval close to 10 femtoseconds (400 a.u. of time, or about 1.25 vibrational periods of the ground state H2). The only exception are the trajectories in Fig. 3B, where the propagation time was doubled.

Figure 3

Different types of trajectories for the vibration of the H2 molecule. All trajectories are plotted in the phase space (momentum vs coordinate): (A) in the low-energy regime of nearly harmonic motion, (B) in the high-energy regime of strongly anharmonic motion, and (C) at an energy above the dissociation threshold that leads to bond breaking. The exact analytical solution (dashed black), a classical QP solution (dashed magenta) and solutions obtained using the D-Wave quantum annealer with (blue) and without (red) greedy postprocessing are presented for each case.

Figure 2 shows several solutions for a low energy trajectory which starts at Bohr with zero initial momentum (), all of which are computed with 1000 equal time steps, i.e. attoseconds. This bond length is close to the equilibrium position, so the vibrational motion is expected to be fairly harmonic. The black dashed line in Fig. 2 shows the exact analytical solution available for this problem (see the Supplemetary Information for details). The red line shows a direct application of our method, using the D-Wave quantum annealer to propagate one equation at a time through one time step (). As one can see, the accuracy of this trajectory is rather poor.

Figure 2

Trajectory for the vibration of the H2 molecule in the low-energy regime of nearly harmonic motion. The trajectory is plotted in the phase space (momentum vs. coordinate). The points are connected in the order of time and start from Bohr with . The exact analytical solution (dashed black) and three solutions obtained on a quantum annealer (color lines) are presented.

There are multiple sources of error that can contribute to this result. First, the QUBO matrix translates into the physical parameters of the qubits with limited accuracy[22], which can be regarded as errors in the QUBO matrix itself. The limited precision and range of the weights when mapped from the QUBO impacts performance. Second, D-Wave’s limited sensitivity to small differences in energies corresponding to similar solutions together with the probabilistic nature of quantum measurements makes it impossible to guarantee selection of the true lowest energy solution. Finally, additional errors may be contributed by imperfect hardware implementation (e.g. noise and decoherence).

One way to improve the quality of D-Wave’s results is to simply restart the calculation of each trajectory point multiple times and select the best solution among these runs (i.e., the one with the smallest ) as the answer. The effect of this strategy is shown with the green line, where each task was allowed to be solved up to 10 times. Comparing these two trajectories, one can see that restarting improves the quality of the D-Wave’s results quite dramatically. The values of , and other settings in this case are the same as for the red line, which confirms that the relatively poor results of direct application of the D-Wave are primarily due to random internal errors of the annealer.

Another way to improve the accuracy is to allow some classical postprocessing. The blue line in Fig. 2 shows an example of a trajectory where each solution found by D-Wave was used as a starting point for a classical greedy algorithm (available as a part of D-Wave’s Ocean tools)[32], which then tried to improve the solution found by the quantum annealer. Both equations were propagated at once in this case (, which is strictly harder than ) and only one attempt per task was allowed, but this was sufficient to obtain an even better trajectory. We conclude that this hybrid quantum–classical approach is the most efficient and accurate across the three scenarios considered here.

The next set of results is presented in Fig. 3, where we increased the time step ten-fold (i.e., attoseconds). In this picture each point of the time-grid is shown by a symbol, to show the discrete nature of the numerical solution. In addition to the analytical and the D-Wave’s results, we also present a trajectory obtained with the QP version of our algorithm, by minimizing the target functional in the form of Eq. (8) on a classical computer. This method is free of the technical issues that are still present in the D-Wave annealer. In this case it serves as the best solution that one can achieve within the framework of our QDE algorithm at a chosen level of discretization, if all minimizations are carried out correctly. Figure 3A represents the same trajectory as Fig. 2 and shows that increasing the time step by a factor of 10 leads to a noticeable deviation of all numeric solutions (including QP) from the exact analytic solution, but the results obtained by the D-Wave quantum annealer remain very close to the QP result. We also tried to propagate trajectories with the high level of vibrational excitation, when the motion is expected to be highly anharmonic. Figure 3B illustrates one example of such a trajectory, where the energy of H2 was close to 50% of its dissociation energy ( Bohr, which corresponds to a compressed bond, ). Figure 3C gives another example, where the energy of H2 is chosen above the dissociation threshold and this leads to a relatively fast bond-breaking, or dissociation of the molecule into atoms ( Bohr in the repulsive range, ). The level of discretization and the methods/colors are the same as in Fig. 3A. We see that the results of a combined quantum–classical approach (D-Wave + greedy) overlap with the QP results in all three frames of Fig. 3.

When comparing Figs. 2 and 3, one can see how multiple components of the overall error influence the quality of solutions. In Fig. 3, the overall error is dominated by the propagation error of QDE discretization (i.e., large time step), so the internal D-Wave’s error (i.e., error of minimization of the target functional) is not noticeable and both QP and D-Wave’s solutions look similar. In contrast, in Fig. 2 the time step is small, so the propagation error becomes negligible, and the internal D-Wave’s error becomes apparent, which is responsible for the “ragged” look of the D-Wave’s solution in this case.

To provide more insight into the performance of the QDE method itself, and of its execution on the D-Wave annealer, we have analyzed these trajectories on a more quantitative level. In Fig. 4 we demonstrate convergence properties for several approaches we tried, by plotting their values of Residual Mean Squared Error (RMSE, see the Supplemetary Information) with respect to the analytical trajectory, as a function of the total number of grid points or time steps (larger number of grid points corresponds to smaller time step ). All data in Fig. 4 were obtained using the simplest version of the QDE algorithm (, when one equation is propagated at a time through one time step). In these and further tests, the value of was set to 1.3 Bohr and , the same as in Fig. 2.

Figure 4

RMSE of solutions obtained with different methods as a function of the total number of grid points (or time steps). All methods propagate one equation at a time.

First, we solved Eq. (1) classically, as a continuous QP-problem (dashed magenta line in Fig. 4). Once again, the QP-solution provides a helpful reference because it is exact and shows the minimum error that can be achieved for a given number of grid points, if the function of Eq. (8) is minimized correctly. One can see that our method converges quickly with respect to the total number of grid points, which is the only convergence parameter that we have here. Next, we found solutions in binary variables (similar to qubits) using the QBSolv software tool also available as a part of D-Wave’s Ocean tools[32], which uses a classical heuristic probabilistic algorithm to find the minimum of Eq. (11) in binary variables. Finding the global minimum in discrete variables is harder than in continuous variables, therefore, in contrast to the QP solver, QBSolv is not guaranteed to find the correct solution. However, as one can see from Fig. 4, the solutions obtained with QBSolv (dashed line) nearly coincide with the ideal QP-solutions, so this algorithm is working well for small QUBO-matrices (). This test demonstrates that the correct solution to Eq. (1) can be obtained in binary variables if minimization of Eq. (11) is done sufficiently accurately.

As a next step, we used the D-Wave to minimize Eq. (11), instead of the classical solver QBSolv. As one can see from the red line in Fig. 4 (which also corresponds to the red line in Fig. 2 for grid points), the D-Wave’s result is close to the results of QP and QBSolv for small grids, but deviates further for larger grids, unable to achieve low values of RMSE. Loosely speaking, we define an “accurate” trajectory as one with RMSE < 10–2 Bohr, which means that in this test the D-Wave failed to provide any acceptable solution. From the QBSolv’s result we know that better solutions exist, but D-Wave’s annealer was not able to find them. We can conclude that the D-Wave has internal errors or noise that results in an RMSE on the order of 10–2 Bohr (for our problem) and is responsible for the deviation of the red line observed in Fig. 4. We tried to improve this result by varying annealing parameters, such as annealing time and chain strength[28], but this did not lead to any significant changes in the errors of trajectories. The only parameter that reliably improved the quality of solutions in this test was the number of reads, which was set to the maximum allowed value of 10,000 for all D-Wave’s results in this work.

Finally, the green line in Fig. 4 corresponds to the green line in Fig. 2, i.e., to the case when each point was restarted multiple times (up to 10 attempts). This may be similar to increasing the number of reads beyond the maximum value, depending on the D-Wave’s internal implementation of it. We see that this simple fix reduces the RMSE of the solutions by almost an order of magnitude, and produces a sufficiently accurate trajectory.

In the next computational experiment, reported in Fig. 5, we tried to propagate the equations for position and momentum together in one run (). Once again, the dashed magenta line shows perfect results of the QP solver, which can be used as a reference line. The results of QBSolv with 1 attempt per task (cyan line) in this case are significantly worse, but this is expected since the algorithm is classical, and the size of the search space increases exponentially with the size of QUBO matrix. However, similar to what we saw for D-Wave in Fig. 4, the results of QBSolv can be significantly improved simply by restarting each task several times (yellow line) and improved even further if one rescales the ranges of two variables in the two equations (brown line) to equalize contributions of position and momentum to the minimization functional (see the Supplemetary Information). In this final form, the accuracy of QBSolv approaches that of QP.

Figure 5

RMSE of solutions obtained with different methods as a function of the total number of points in the grid. All methods propagate both equations at once.

Returning back to the performance of the D-Wave annealer, we see that even with 10 attempts per task (green line in Fig. 5) the results are much worse compared to the analogous results for the case when only one equation was propagated at a time (green line in Fig. 4). This is in sharp contrast with classical QBSolv and could indicate that the overall quality of the quantum annealing hardware is reduced when the number of qubits is increased for a given problem. We suspect that such dependence could exist due to the fact that only a few other qubits are directly connected to any given qubit on the D-Wave and involvement of a larger number of qubits translates into longer qubit connection chains, which are imperfect. Surprisingly, rescaling the variables in the two equations did not help in this case either.

However, we found that a hybrid strategy (when each solution found by D-Wave was used as a starting point for the classical greedy algorithm) works really well for this harder problem (blue line in Fig. 5, which also corresponds to blue line in Fig. 2). This result matches the quality of QP solutions even for a large number of grid points (i.e., small time steps and high accuracy). One might think that the greedy algorithm is doing all the work here and that the involvement of D-Wave is unnecessary. In order to rule out this possibility, we tried to lower the quality of D-Wave’s solutions by significantly reducing the number of reads, which basically gave a near-random starting point for the greedy algorithm. We found that in this case greedy algorithm does not perform well, which indicates that the D-Wave’s initial guess is important here. Notably, while the purely classical algorithm (QBSolv) was also able to achieve the same results, the hybrid quantum–classical strategy (D-Wave + greedy) did it much more easily with only one attempt per task and even without rescaling, which is quite optimistic. One explanation of the success of this hybrid approach is that D-Wave is good at identifying the global minimum, but it cannot descend all the way there with a high precision due to the presence of noise and related errors. On the other hand, greedy algorithm is good at descending to the very bottom of the minimization functional but needs to be placed in the vicinity of the global minimum first. Thus, the two methods complement each other and together converge to a very good result.

Practicality of presented numerical algorithms is obviously defined by the scaling of hardware requirements with the system size. Overall, the resources required by QDE grow as , where is the total number of points in the grid, which depends on the desired accuracy, and is the total number of equations in the system, which increases linearly with the dimensionality of the problem. For example, the new D-Wave Advantage quantum annealer with 177 fully connected logical qubits[26] could be used to propagate up to 8 equations in one run, which is sufficient to describe the vibrational motion of any triatomic and some tetratomic molecules.

For more complicated systems (large number of equations due to large number of atoms), the trajectories could be propagated by splitting the overall workload into multiple runs. As long as the number of available qubits is sufficient to propagate at least one time step of one equation, we can propagate each equation and each time step sequentially, one after another, since all equations are decoupled in the limit when only one time step is propagated at a time (see Theory section for more details). This way, the number of qubits necessary for propagation can be constant regardless of system size or propagation time, but the total number of runs would grow instead. Having more qubits would certainly help by allowing to solve more equations or more time steps in a single run of a quantum annealer, thus decreasing the overall wall clock time, but is ultimately not a requirement.

In cases when the right-hand side functions of Eq. (2) are non-linear, our method will not be able to obtain all points of a solution at once, even on annealers with sufficiently large number of qubits, since only quadratic functions can be minimized on such devices. As explained above, this can be easily circumvented by running the same algorithm multiple times. Ultimately, the minimum required number of consecutive runs will be dictated by the number of segments approximating , which depends on desired accuracy of a solution and the degree of non-linearity of the . However, we can still make use of all available resources of a given (large, if available) annealer by increasing the number of grid points in each segment. Note that this restriction only concerns the dependence of on ; the dependence on can still be arbitrary and cause no problems. Although the methods to simulate higher-order polynomial behavior on quantum annealers exist and can be used to approximate with segments of higher order to reduce the total number of necessary runs, this does not necessarily work better for arbitrary functions and was not discussed here.

An important advantage of the presented methodology over classical propagation techniques is the ability to obtain multiple trajectory points simultaneously rather than sequentially. This opens a potentially exciting route of calculating the entire trajectory in just a few runs, given a sufficiently large number of qubits.

Conclusions

To summarize, in this work we reported first-ever formulation of molecular dynamics as a Quadratic Unconstrained Binary Optimization (QUBO) problem and its solution on a D-Wave quantum annealer. Importantly, QUBO problems can be solved not only on quantum annealers but also on existing universal gate-based quantum computers, and on exotic classical devices (e.g. Fujitsu digital annealer). The new methodology developed here (the QDE method) was applied to propagate classical trajectories simulating vibrational motions of the hydrogen molecule (H2) in three different energy regimes (nearly harmonic, highly anharmonic, and dissociative limits). The results obtained using the D-Wave 2000Q quantum annealer are mutually consistent and quickly converge to the analytical solution. Several alternative strategies for such calculations were explored and it was found that the most accurate results and the best efficiency are obtained by combining the quantum annealer with classical post-processing (greedy algorithm). Importantly, the QDE framework developed here is entirely general and can be applied to solve any system of first-order ordinary nonlinear differential equations using a quantum annealer. The new generation of quantum annealers, such as the D-Wave Advantage with more qubits and better connectivity, could be used to either compute multiple time-grid points at once () or to explore more complicated molecules with many degrees of freedom (). Higher order finite difference schemes can also be useful in some applications. Development of trajectory based quantum algorithms, suitable for execution on a quantum annealer, is another promising potential future direction.

Supplementary Information 1.

Supplementary Information 2.

Supplementary Information 3.