Communication: Constant uncertainty molecular dynamics: A simple and efficient algorithm to incorporate quantum nature into a real-time molecular dynamics simulation Free

Quantum nature of nuclei governs the tunneling of the protons and many chemical reactions. Furthermore, the vibrational spectra are highly affected by the nuclear quantum effects. Nevertheless, most of the simulations have been conducted under the approximation that a nucleus can be dealt with as a classical particle because the simulation of real-time nuclear quantum dynamics of large and complex molecular system in condensed phase is still computationally and theoretically challenging. A major route to simulate the real-time nuclear quantum dynamics in condensed phase is to use the imaginary-time path integral molecular dynamics (PIMD).^1–5 Although it was formulated to predict static thermodynamic properties, the PIMD formalism can be used to predict dynamical properties. The most successful methods are the imaginary-time path integral centroid MD (CMD) and the ring-polymer MD (RPMD).^3,5,6 These methods utilize the fact that the correlation function of the average of imaginary-time path beads is an approximation of the Kubo-transformed quantum real-time correlation function.^3,5,7 Another method is the linearized semi-classical initial value representation (LSC-IVR), or the classical Wigner approximation, in which the system evolves along real-time.^8–11 While these trajectory-based methods are practical for simulating the real-time quantum dynamics in condensed phase, they suffer from a common drawback in their nature, the absence of quantum nature in real-time evolution. This is because all the quantum effects considered in the above methods arise from the quantum Boltzmann operator. In this communication, we call these effects quantum imaginary-time statistical effects as distinguished from the quantum nature in real-time evolution. The lack of the quantum nature in real-time evolution imposes some limitations. For example, the LSC-IVR method employs classical dynamics in real-time evolution, so that it suffers from the so-called zero-point energy (ZPE) leakage problem.^12,13 Apart from the trajectory-based methods, wave packet/wave function-based methods were used to calculate the real-time quantum dynamics.^14,15 The quantized Hamiltonian dynamics (QHD) method,^16,17 which utilizes the quantum expectation values as the basic dynamical variables, was also developed to simulate the real-time quantum dynamics of molecular systems. The quantum nature in real-time evolution is inherently included in these methods. However, the wave packet/wave function-based methods and the QHD method are computationally too expensive to evaluate potential functions for complex and large molecules in condensed phase. A simple and efficient MD algorithm is therefore desirable for dealing with the quantum effect in real-time evolution in complex molecular systems. The ZPE conserving-trajectory approaches^18–27 and the entangled trajectories MD approach^28,29 have been shown to capture the quantum effect in real-time evolution to some extent. While the applications of these approaches are limited to rather small systems, the adequacy of the strategies of these approaches encourages us to develop the semi-classical algorithm which utilizes basic quantum properties to regulate the dynamics of the classical systems. Then we developed a novel formalism to deal with the quantum effect in real-time evolution. Our formalism focuses on the basic property of the quantum mechanics, that is, the uncertainty principle between position and momentum. The idea is that performing the MD simulation of the ensemble of classical systems under the constant uncertainty constraint on its dynamics so that the quantum nature in real-time evolution can be taken into account. In this communication, we report that the real-time quantum dynamics can approximately be handled by the simple and efficient new algorithm, which we call the constant uncertainty molecular dynamics (CUMD).

Consider an ensemble consisting of N copies of a one-dimensional system which has a classical particle with mass m on a potential energy V(q). The time evolutions of the particles in the ensemble are described by a set of Hamilton’s equations,

where i is the index of the ith system in the ensemble. We then define the time-dependent ensemble average of the function A(p, q) at time t as

where we define $p 0 = ( p 1 0 , p 2 0 , … , p N 0 )$ and $q 0 = ( q 1 0 , q 2 0 , … , q N 0 )$ as the vectors of momenta and positions of the classical particles in the ensemble at t = 0. Since the time evolutions of the systems in the ensemble are classical, the momenta p_i(t) and positions q_i(t) obey Eq. (1). Assuming that this ensemble average formally corresponds to the Wigner function of the operator $A ˆ$⁠,^30,31 we can write the uncertainty of the ensemble of the classical particles as

The expectation value $p q$ formally corresponds to the symmetrized quantum operator $p ˆ q ˆ + q ˆ p ˆ /2$ according to the definition of the Wigner function. Schrödinger’s uncertainty relation is employed instead of Heisenberg’s one to derive Eq. (3) because former is a more general principle and includes the correlation effect between position and momentum.³² Note that the quantum uncertainty condition g(t) ≥ ħ²/4 is not necessary to be satisfied for the ensemble of the classical particles obeying Eq. (1), which gives rise to the difficulty in treating quantum dynamics with the ensemble of the classical trajectories. We therefore consider the modified classical dynamics where the uncertainty constraint is imposed on the dynamics of the classical systems. Since the exact time dependence of the quantum uncertainty is not tractable, we instead assume the ensemble of the systems which have the time-independent uncertainty. This constant uncertainty assumption is justified by the fact that the lowest-order quantum correction to the classical trajectory appears as g(t) = ħ²/4 in perturbative analysis.³³ The constant uncertainty constraint is written as

where g(0) is the initial uncertainty of the ensemble. Then the time derivative of the constraint can be written as

Hereafter we drop the argument t for simplicity. Since Eq. (5) is linear with respect to the momenta, the Lagrange multiplier method can be used to regulate the dynamics of the classical systems. Following the standard procedure to treat a nonholonomic constraint,^34–36 one can derive the equations of motion of the constant uncertainty dynamics as

where the superscript CUMD denotes that the time evolution of the variable obeys Eqs. (6) and (7), i.e., it is calculated by the CUMD simulation, λ is the Lagrange multiplier which is determined to satisfy the constant uncertainty condition g^CUMD(t) − g(0) = 0, and $A CUMD$ is the time-dependent expectation value in the CUMD simulation, which is defined as

where $q i CUMD ( p 0 , q 0 , t )$ and $p i CUMD ( p 0 , q 0 , t )$ are the position and momentum at time t, which are calculated by CUMD from the initial positions q⁰ and initial momenta p⁰. While the explicit expression of the multiplier can be available, we employed the velocity-Verlet integrator with the RATTLE-type algorithm to perform the CUMD simulation in this study.³⁷ 

Before carrying out a CUMD simulation, one needs to prepare the appropriate initial positions q⁰ and momenta p⁰ of the classical particles representing a quantum state as the constant uncertainty (CU) ensemble. We can prepare an initial condition which satisfies g(0) = ħ²/4 by choosing q⁰ and p⁰ both of which represent a specific wave function. However, no thermostatting algorithm to equilibrate the prepared initial condition is available for CUMD at present time. Thus we employ a simple ad hoc scheme to generate them at finite temperature by the RPMD simulation with N imaginary-time path beads in this study. First, the initial position of ith copy of the classical system in the CU ensemble is sampled from the RPMD trajectory as

where $q i RPMD$ is the position of the ith bead of the ring-polymer. Then the initial momentum of ith copy of the classical system is calculated as

where $K CL = p RPMD 2 /2m$⁠, p^RPMD represents the momentum of the ring-polymer beads, and we conveniently choose $s i = q i − 1 RPMD − q i + 1 RPMD$ because $s$ and $s q$ become zero. The cyclic boundary conditions $q N + 1 RPMD = q 1 RPMD$ and $p N + 1 RPMD = p 1 RPMD$ are applied. α is determined for the kinetic energy of the initial momenta to be equal to the sum of the instantaneous classical kinetic energy K^CL and the quantum kinetic energy K^QM, namely, $p 0 2 /2m= K CL + K QM$⁠. That is, the values of α are the pair of solutions of the quadratic equation,

where s = (s₁, s₂, …, s_N). We use the centroid virial estimator to calculate the quantum kinetic energy, which is given by^38,39

While there are several options for the preparation of the value of K^QM, a time-averaged value of K^QM,virial over the total RPMD trajectory is used for K^QM in this study. The initial uncertainty $g 0$ is automatically determined by the initial momenta and positions in the ensemble.

In this approach, the Kubo-transformed quantum correlation function of operators $A ˆ$ and $B ˆ$ can be written as

where β_n = β/N, β is the inverse temperature, and

and

are the ring-polymer Hamiltonian and the partition function, respectively.⁵ Note that, although this expression is related to the ring-polymer representation of the LSC-IVR correlation function,⁴⁰ there are two differences between LSC-IVR and our CUMD-based approaches. First, the quantum interference effect between initial positions and initial momenta which is accounted for in the LSC-IVR method is neglected by using the time-averaged quantum kinetic energy. Second, the time evolutions of the particles are governed by the constant uncertainty dynamics instead of the classical dynamics which is used in the LSC-IVR method.

In order to test the performance of the CUMD method, we calculated the Kubo-transformed quantum position autocorrelation functions of a quartic oscillator which has the potential $V q = q 4 /4$ and the mass m = 1 in atomic unit (a.u.). The simulations were carried out with the classical MD, RPMD, and CUMD methods with and without the uncertainty constraint at high (β = 8 a.u.) and low (β = 1 a.u.) temperatures. We used 64 imaginary-time path beads for the high temperature case and 128 beads for the low temperature case. The velocity-Verlet integrator with the time step of 0.0025 a.u. was used to solve the equations.

Figure 1 shows $C q q K ( t )$ for a quartic oscillator at a low-temperature (β = 8 a.u.) and at a high-temperature (β = 1 a.u.). In the low-temperature case (Fig. 1(a)), while the classical result far deviates from the exact quantum one even at t = 0 a.u. and fails to reproduce the quantum oscillation period, the CUMD result is in good agreement with the exact quantum one especially in short time. If the constraint is removed from the CUMD simulation, no quantum oscillation is observed at t > 10 a.u. This indicates that, although the initial CU ensemble possess the quantum imaginary-time statistical effect, it fades out along the real-time evolution and totally disappears at t > 10 a.u. due to violation of the uncertain principle of the ensemble of the classical trajectories. On the other hand, the oscillation of the CUMD result at t > 10 a.u. is successfully represented by the thermally averaged CU ensemble. In the high-temperature case (Fig. 1(b)), since the systems have more classical nature when their temperatures increase, all the results, including the classical result, are in better agreement with the exact quantum one than those in the low-temperature at t < 5 a.u. However, only the CUMD and the exact quantum results retain their oscillating behaviors at t > 10 a.u.

Next we compare the CUMD and the RPMD results. At the low temperature (Fig. 1(a)), CUMD reproduces the exact quantum result better than RPMD does because of the inclusion of the quantum uncertainty in real-time evolution of the CU ensemble. At the high temperature (Fig. 1(b)), RPMD does not reproduce the quantum oscillation after the first oscillation period and the result is similar to the classical one. This failure is explained by the fact that, at the high temperature, the quantum effects in real-time evolution, which is neglected in RPMD, are more important than the quantum imaginary-time statistical effect.⁵ Since CUMD can approximately describe the quantum real-time dynamics through the conserved quantum uncertainty of the ensemble, the CUMD result shows the oscillation even at t > 10 a.u. where the RPMD result shows no oscillations.

Finally, we check how accurately CUMD reproduces the vibrational transition frequency because the accurate description of ZPE is critical to reproduce the correct quantum fundamental transition frequency of an anharmonic oscillator.^41,42 The Fourier transformed Kubo-transformed position autocorrelation functions $C q q K ( ω ) = ∫ − ∞ ∞ e − i ω t C q q K ( t ) dt$ of a Morse oscillator were calculated. The potential parameters and mass were chosen so that the quantum fundamental transition frequency becomes 3500 cm⁻¹ and the frequency difference between the 0-1 vibrational transition and the 1-2 vibrational transition, so called the anharmonicity, was to be 150 cm⁻¹. These parameters led to the classical fundamental frequency of the Morse oscillator to be 3650 cm⁻¹. The velocity-Verlet integrator with the time step of 0.05 fs was used. The spectra were calculated by the fast Fourier transformation with the Hann window. In Fig. 2, the calculated $C q q K ( t )$ are shown. The classical spectrum shows the peak at ∼3650 cm⁻¹ because the classical simulation totally discounts ZPE. The RPMD result shows a peak at ∼3580 cm⁻¹, which is a slightly higher frequency than the exact quantum transition frequency. This disagreement is explained by the inaccurate evaluation of the instantaneous quantum kinetic energy in the RPMD simulation.⁴¹ The both CUMD results show a peak at ∼3500 cm⁻¹, which is in excellent agreement with the exact quantum transition frequency. This is not surprising because the initial quantum kinetic energy is accurately estimated by the centroid virial estimator from the total RPMD trajectory and preserves its quantum nature during the real-time evolution by CUMD. This relevant role of the constant uncertainty constraint can be seen in the result without the uncertainty constraint. The broad peak indicates that the initial quantum imaginary-time statistical effect gradually disappears and hence the oscillation frequency becomes gradually blue-shifted during the time evolution. We note that the CUMD results show an additional broad peak around ∼3850 cm⁻¹. Similar to the case of the Gaussian wave packet approach,⁴³ this artificial peak may be attributed to the vibration of 〈q²〉 around 〈q〉 due to the uncertainty constraint. However, further detailed investigation should be needed for this point.

In this communication, we reported the new MD algorithm for approximately dealing with the real-time quantum dynamics by means of the constant uncertainty constraint on the dynamics of classical systems. The algorithm was implemented to the conventional MD with the initial conditions which were prepared by the RPMD simulation. The CUMD method can describe the quantum correlation function better than the RPMD method due to the inclusion of the quantum effect in real-time evolution as well as the quantum imaginary-time statistical effect in the initial condition. We also showed that the transition frequency of the Morse oscillator was correctly reproduced by the CUMD simulation because ZPE was preserved under the constant uncertainty dynamics. The new method is easy-to-implement and computationally not so expensive. Because the method can deal with the ZPE effect accurately in real-time evolution, versatile applications of the CUMD simulation are expected for the future studies, such as the vibrational energy relaxation process and the vibrational spectrum calculation in condensed phase.

T.H. would like to thank S. Hayashi, Y. Nagata, and K. Abe for valuable comments on the manuscript. This research was supported by MEXT as “Priority Issue on Post-K computer” (Building Innovative Drug Discovery Infrastructure through Functional Control of Biomolecular Systems).