{'fraction_non_alphanumeric': 0.032130396292246004, 'fraction_numerical': 0.009998437744102483, 'mean_word_length': 5.112030553787396, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Incorporating Nuclear Quantum Effects in Molecular Dynamics

The accurate incorporation of nuclear quantum effects in large-scale molecular dynamics simulations remains a significant challenge. Here, we develop a new formulation of the equations of motion for molecular dynamics that utilizes an effective energy surface obtained from a constrained energy minimization with given quantum nuclear expectation positions, thus enabling the incorporation of nuclear quantum effects in molecular dynamics simulations. With a series of model systems, we demonstrate that this new approach is significantly more accurate in describing vibrations than classical molecular dynamics and can also accurately describe tunneling, which classical molecular dynamics is unable to do. It is comparable to or more accurate than centroid molecular dynamics and ring-polymer molecular dynamics. Through a combination with constrained nuclear electronic orbital density functional theory, this new approach is promising for practical molecular simulations with nuclear quantum effects incorporated.

poration of nuclear quantum effects in molecular dynamics simulations. With a series of model systems, we demonstrate that this new approach is significantly more accurate in describing vibrations than classical molecular dynamics and can also accurately describe tunneling, which classical molecular dynamics is unable to do. It is comparable to or more accurate than centroid molecular dynamics and ring-polymer molecular dynamics. Through a combination with constrained nuclear electronic orbital density functional theory, this new approach is promising for practical molecular simulations with nuclear quantum effects incorporated.
Nuclear quantum effects (NQEs) [1] have a great impact on the structural, thermodynamical, and kinetic properties of a wide range of chemical and biological systems [2]. They usually include zero-point and tunneling effects and are significant when light nuclei, such as hydrogen, are present. The accurate incorporation of NQEs in molecular simulations is important for understanding many fundamental properties but remains a significant challenge for large-scale molecular simulations. For example, the anomalous properties of water are closely related to the NQEs of the complex hydrogen bond network [3,4] and thus cannot be fully explained with classical molecular dynamics (MD) without an accurate inclusion of NQEs [5].
There have been many theoretical developments on the incorporation of NQEs in molecular simulations. Some empirical force fields [6] have been used to include NQEs implicitly, and have been able to treat protons and deuteriums differently in water [7]. Quantum wave packet dynamics is based on the exact time evolution of a quantum system according to the time-dependent Schrödinger equation and can give theoretical predictions that accurately match to experiments [8][9][10][11][12]. Quantum trajectory methods [13] are based on the de Broglie-Bohm formulation of quantum mechanics [14][15][16][17][18], which attributes quantum effects to the quantum potential. With a reasonable approximation to the quantum potential, quantum trajectory methods have been applied to many model systems and give accurate results [19,20]. Multicomponent quantum theories [21][22][23][24][25][26] also include NQEs in real-time simulations [27][28][29][30][31]. They usually simultaneously treat both electrons and key nuclei quantum mechanically and therefore do not rely on conventional Born-Oppenheimer potential energy surfaces (PESs). Dynamics simulations can be performed through quantum time evolution of multicomponent wave functions or density matrices [29,30], and practical problems, such as proton transfer processes, have been studied with multicomponent quantum theory [32].
Although the aforementioned methods are highly accurate in describing NQEs, they are often hindered by their high computational costs in large molecular or bulk systems. This challenge can be partially addressed using methods based on the path integral formulation of quantum mechanics [33,34]. By simultaneously simulating a set of coupled replicas for a system, path integral molecular dynamics (PIMD) [35,36] is able to capture NQEs and accurately describe static properties of the system [37]. Its extensions such as centroid molecular dynamics (CMD) [38][39][40][41][42][43] and ring-polymer molecular dynamics (RPMD) [44] can describe dynamical properties using approximate correlation functions. However, while dynamical properties from CMD and RPMD are considerably more accurate than those from classical MD, challenges still exist with the curvature problem in CMD and spurious frequencies in RPMD. Both of these problems can lead to unreliable vibrational spectra [45], although several recent developments can mitigate them to some extent, including thermostatted RPMD [46,47], Matsubara dynamics [48], and quasicentroid molecular dynamics [49,50].
In this Letter, we present an alternative formulation of the equations of motion for classical MD that incorporates NQEs. With this formulation, NQEs can be accurately described with MD on an effective PES, which in practice is approximated by a constrained minimized energy surface (CMES). We first analytically show that CMES-MD remains exact for the harmonic oscillator model. Then, with numerical analyses of a Morse oscillator model and a quartic double-well potential model, we show that CMES-MD is generally much more accurate in describing vibrations and tunneling effects than conventional MD and is comparable to or better than CMD and RPMD.
We start with the polar representation of a time-dependent wave function ψ(x, t) = A(x, t) exp(iS(x, t)/h), where the amplitude A and the phase S are both real functions. For simplicity, we assume that there is only one quantum particle throughout our derivation.
However, the formulation can be easily generalized to multiple quantum particle cases if the particles can be assumed to be distinguishable, such as nuclei in regular molecular and bulk systems. With this polar representation, the kinetic energy can be decomposed into two terms A(x, t) The first term is the kinetic energy evaluated with the amplitude function A only. Since A is associated with the real space probability density distribution with ρ( this term can be perceived as the kinetic energy due to quantum delocalization, or the zero-point kinetic energy. In the second term, the key quantity ∇S is associated with the observable momentum and is in fact the momentum field in Bohmian mechanics [14][15][16][17][18] with the definition p(x, t) = ∇S(x, t). Therefore, after multiplying [∇S(x, t)] 2 /2m with the probability density A 2 (x, t) and integrating over the spatial variable, the second term can be viewed as the kinetic energy associated with the observable momentum p.
With the variance of the observable momentum defined as the variance of the momentum field in Bohmian mechanics: the kinetic energy can be further expressed as The terms in Eq. 3 correspond to, respectively, the zero-point kinetic energy, the classical kinetic energy associated with the expectation value of the observable momentum, and an energy contribution from the variance of the observable momentum.
Another way of expressing the kinetic energy is simply T (t) = Ĥ (t) − V (t), which can be plugged into the left side of Eq. 3. Then we can take the time derivative on both sides of the equation and simplify the equation into Note that we have used the relationship d Ĥ (t)/ dt = ∂V /∂t in the simplification procedure. Eq. 4 relates the time dependence of momentum to the time dependence of energetic terms. It is exact without any approximation. However, in order to make it into an equation of motion that can be practically used in MD simulations, we next proceed with an approximation that builds a connection between quantum states and the classical phase space.
Conventionally, when assuming the potential is slowly varying in space, the Ehrenfest theorem provides such a connection between the classical Newtonian dynamics in the phase space (X, P ) and the evolution of quantum expectation position and momentum ( x , p ) Here we build on the same mapping philosophy but instead of assuming the behavior of the potential, we approximate the quantum state as the energy-minimized state for a given phase space point. That is, when the system is at a particular phase space point given by an expectation position and momentum combination, i.e., ( x , p ) = (X, P ), the quantum state |A always adapts to the energy-minimized state for that phase space point. We note that this approximation is essentially an adiabatic approximation during the dynamics process, and it is not trivially justifiable, however, to keep the flow of the derivation, we leave discussions of its applicability as well as limitations for the later part of this Letter.
Under this approximation, the quantum state |ψ becomes an explicit function of (X, P ) but an implicit function of time t, i.e., |ψ (X(t), P (t)). At a particular phase space point (X(t), P (t)) the system evolves to at time t, this state can be obtained from a constrained minimization for the total energy with the following Lagrangian: where f is the Lagrange multiplier associated with the expectation position, v is the Lagrange multiplier associated with the expectation momentum, andẼ is the Lagrange multiplier associated with the wave function normalization. This Lagrangian can be further expressed in terms of A and S by Making the Lagrangian function stationary with respect to the variation of ∇S and A leads Further combining these equations with the expectation position constraint, the expectation momentum constraint, and the normalization constraint gives v = P (t)/m, ∇S(x) = mv = P (t), and the eigenvalue equation can be simplified to The eigenvalueẼ + P 2 /2m, eigenstate |A , and the Lagrange multiplier f can be solved under the expectation position and normalization constraints for |A . Note that the solution of the amplitude function A only depends on the expectation position constraint, and the expectation momentum constraint only affects the phase function S.
The fact that the constrained minimization requires ∇S to agree with the expectation momentum (∇S(x) = P (t)) naturally leads to σ 2 p = 0 according to the definition in Eq. 2, and with the quantum state |A obtained as an explicit function of X and thus an implicit Note that here we have used dX/ dt = d x / dt = p /m. According to classical mechanics, it is reasonable to assume that the change of p should have an opposite direction to the energy gradient term (∇ X A|Ĥ(t)|A ), the common prefactor | p |/m can be dropped, and we arrive at the final expression where V CMES is a constrained minimized energy surface (CMES) associated with the amplitude part |A . It can also be viewed as an effective potential energy surface that includes not only the potential energy but also the quantum delocalization kinetic energy. Eq. 12, together with d x / dt = p /m, forms the equations of motion for CMES-MD. These equations of motion represent a new way of performing MD simulations that incorporates NQEs.
They are highly similar in structure to Newton's equations used in conventional MD simulations, with the difference that the time evolution is now on the quantum expectation positions and momenta rather than the classical ones.
motion within the framework of Feynman's path-integral formulation of quantum mechanics. Within this framework, the effective potential guides the motion of the ring-polymer centroid, and is claimed to be equal to the zero-temperature limit of the centroid potential for CMD. While this set of equations was used to gain insight into the behavior of CMD, its broader potential for practical MD applications was not realized. In this work, we derived this set of equations starting from the time-dependent Schrödinger equation with the effective potential utilized to guide the motion of the quantum wave function expectation values in the classical phase space. Furthermore, this work is the theoretical foundation for our recently published work [54], in which we demonstrated that, this framework is especially powerful in real systems with the effective potential obtained from constrained nuclear electronic orbital density functional theory (cNEO-DFT) [55][56][57].
In this Letter, we investigate the model systems of a harmonic oscillator, a Morse oscillator and a quartic double-well potential, which have easily accessible exact quantum solutions and can be conveniently compared with CMD and RPMD. For the Morse oscillator and the quartic double-well models, classical MD and CMES-MD are performed with an inhouse python script, and 1-dimensional CMD and RPMD simulations are performed with a modified i-PI package [58]. We note that solving for the wave function A in each step for these models is computationally unfavorable, and therefore we pre-calculate the CMES on This effective potential universally shifts the original harmonic potential upwards byhω/2, which is the well-known zero-point energy for a harmonic oscillator (see Fig. S1). This result may seem counterintuitive since conventionally ZPE is considered to be a property of the whole energy surface rather than a point-wise property, however, we note that here the ZPE should be more accurately considered as a quantum delocalization energy, which always exists as the quantum wave packet travels through space. Since classical MD produces the exact trajectory on the harmonic potential, the trajectory produced by CMES-MD on V CMES will also be exact without any need for numerical tests.  zero-temperature limit of CMD centroid potential [51,52], we observe that CMES-MD at finite temperatures is in fact more accurate than RPMD and CMD with better agreement with the exact quantum correlation functions and more accurate vibrational frequencies.
As the temperature increases, CMES-MD and other simulation methods start to have broader vibrational peaks. For classical MD and CMES-MD, this peak broadening is accompanied by a red shift in the peak position, while RPMD and CMD have relatively more stable peak positions. For this Morse potential, we observe that CMES-MD performs well for a temperature range between 50 K and 1500 K (see Fig. S2). Note that although 1500 K is considered high temperature for most chemical and biophysical processes, it is still low temperature for this Morse oscillator that mimics the stretch of the O -H radical molecule, since 3700 cm −1 roughly corresponds to 5400 K. These results suggest CMES-MD is a reliable method in a large temperature range that most chemical and biophysical reactions are performed at.
Next we investigate a more challenging double-well potential model, in which quantum tunneling is expected to occur. We use a quartic double-well potential with a 0.125 eV barrier height and a 0.5Å separation between the potential minima, which can roughly represent the potential energy landscape for a practical proton transfer reaction. As shown in Fig. S3, the CMES of this double-well potential is a single well with the minimum located at X = 0. This huge difference between the original potential and the effective potential has also been observed in previous literatures [51,52]. The reason is that the ground state for this double-well potential has a symmetrical wave function with two peaks, whose expectation position is at x = 0, and therefore the CMES reach its minimum at X = 0. As the constrained expectation position deviates from the center, the constrained minimized wave function becomes less symmetrical with more and more excited state characters mixed in, thus increasing the energy and forming a single-well effective potential. behave very similarly with broad peaks that maximize around 1300-1500 cm −1 (see Fig. S4).
All of these results from CMES-MD show that CMES-MD is accurate and comparable to CMD in describing the dynamics in the double-well potential system.
We note that in principle, the effective potential energy surfaces in MD simulations should be temperature-dependent to fully account for nuclear quantum effects. However, although CMES-MD relies on CMES, which is a temperature-independent effective potential, we observe good performance of CMES-MD over a relatively large temperature range. This suggests that the adiabatic approximation that effectively assumes that the quantum state adapts its real space distribution to the lowest-energy state for a particular phase space point is reasonable. Nevertheless, it is possible that when the temperature is very high and the particle is moving fast, the probability distribution may not adapt fast enough to the constrained minimized wave function, thus breaking the adiabatic approximation.
Therefore, we expect CMES-MD to be accurate at low temperatures relative to the mode frequency. Fortunately, for most vibrational modes, the room temperature is still considered low temperature, and therefore CMES-MD can be accurate in a good range of temperatures typically investigated by chemical physicists and biophysicists.
Similarly to conventional MD, the classical treatment brings not only efficiency but also some limitations. For example, quantum coherence is missing, which is reflected by a decreasing amplitude of the correlation function (Fig. 2). Heat capacities will not approach zero when T → 0 K due to the loss of the energy quantization picture. Furthermore, classical dynamics with distinguishable particles is incapable of capturing the exchange effect, which is important in systems with heavily packed particles, such as in a Bose-Einstein condensate [61,62]. Although detailed studies of these possible limitations due to the approximations made during the derivation are beyond the scope of the current work, they are important topics for our future research for better understanding the applicabilities and limitations of CMES-MD. We finally note that due to the similarity between CMES-MD and classical MD, in practical systems, we can expect analytical force field models or even machine-learning force fields (ML-FFs) [63] to be built based on the CMES, which will allow for an even more efficient incorporation of NQEs in MD simulations.
In summary, we provide a formal derivation for incorporating NQEs in the framework of MD. This is achieved through the calculation of a CMES, which serves as the effective potential for MD simulations. In CMES-MD, quantum delocalization and tunneling effects are inherently included and therefore CMES-MD is capable of accurately describing dynamical vibrational frequencies that are comparable to or better than CMD and RPMD. In future simulations of practical systems, the cNEO approximation can be applied. The resulting cNEO-MD will be vastly more computationally efficient than conventional ways of including    1) and high energy initial configurations from the Maxwell-Boltzmann distribution at 1500 K must be discarded. Therefore, the CMES-MD spectrum is not able to produce lower frequencies as one can see from the classical MD spectrum.