From classical to quantum and back: Hamiltonian adaptive resolution path integral, ring polymer, and centroid molecular dynamics

Kreis, Karsten; Kremer, Kurt; Potestio, Raffaello; Tuckerman, Mark E.

doi:10.1063/1.5000701

Path integral-based methodologies play a crucial role for the investigation of nuclear quantum effects by means of computer simulations. However, these techniques are significantly more demanding than corresponding classical simulations. To reduce this numerical effort, we recently proposed a method, based on a rigorous Hamiltonian formulation, which restricts the quantum modeling to a small but relevant spatial region within a larger reservoir where particles are treated classically. In this work, we extend this idea and show how it can be implemented along with state-of-the-art path integral simulation techniques, including path-integral molecular dynamics, which allows for the calculation of quantum statistical properties, and ring-polymer and centroid molecular dynamics, which allow the calculation of approximate quantum dynamical properties. To this end, we derive a new integration algorithm that also makes use of multiple time-stepping. The scheme is validated via adaptive classical–path-integral simulations of liquid water. Potential applications of the proposed multiresolution method are diverse and include efficient quantum simulations of interfaces as well as complex biomolecular systems such as membranes and proteins.

I. INTRODUCTION

Quantum delocalization of light atomic nuclei plays an important role in many aqueous and soft matter systems, ranging from low temperature helium or hydrogen^1–5 to complex biological systems at room temperature. Examples include proton transfer in biomolecules and membranes,^6–12 thermodynamics of ice,¹³ water adlayers on catalysts,^14,15 aqueous proton and hydroxide transport,^16–21 and even the structure and dynamics of bulk water.^22–26

In computer simulations, nuclear quantum effects are typically modeled using Feynman’s path integral (PI) formulation of quantum statistical mechanics.^2,27,28 The atomic nuclei are mapped onto classical-like (CL) ring polymers, whose beads correspond to the imaginary time slices of the PI. Based on this approach, various techniques have been developed to compute approximate quantum mechanical properties. Path integral molecular dynamics (PIMD)^28–34 and path integral Monte Carlo (PIMC)^28,33,35 directly sample the quantum canonical distribution obtained after path integral quantization and can be employed to calculate time-independent quantum statistical properties. Centroid molecular dynamics (CMD),^36–47 which follows the artificial dynamics of the ring-polymer centroids, and ring polymer molecular dynamics (RPMD),^{45,46,48–51} which is based on the evolution of the individual PI beads, additionally enable the calculation of approximate quantum dynamical properties.

However, PI-based methods are significantly more expensive than corresponding classical simulations. To overcome this, different techniques have been proposed. For example, ring polymer contraction (RPC)^52–55 makes use of the fact that long-ranged and non-bonded interactions typically do not need to be evaluated on as many PI beads as bonded interactions. A related technique is the mixed time slicing scheme,⁵⁶ in which different particles are described with a different number of imaginary time slices. Other approaches include higher-order Trotter factorization^57–61 and advanced thermostating procedures based on generalized Langevin equations (GLEs).^62–65 Additionally, multiple time-stepping (MTS) techniques are frequently employed in PI simulations to decouple the computation of the expensive but slowly varying non-bonded forces and the high frequency internal motion of the ring polymers.^33,34,66

Most of these methods correspond to a modification of the path integral computation itself. For example, in RPC, one evaluates the long-range forces using fewer beads than are used for short-range and bonded interactions. A different approach is provided by adaptive resolution methods, which do not modify the path integral calculation itself but rather aim for a gain in numerical efficiency by restricting the PI description to a small subregion within the simulation box and coupling it with a classical model. The available computational resources can then be concentrated on the quantum mechanical (QM) subregion leading to an overall speedup compared with full QM simulations. This strategy is useful when only a small part of the overall large system actually needs to be described taking into account quantum delocalization effects, which might be the case, for example, in simulations of surfaces, membranes, or the active site of a protein. One such method, based on the adaptive resolution simulation scheme (AdResS),^67–69 is the direct spatial interpolation of a classical force field with the PI-based forces obtained after quantization.^70–74 This method is, however, not compatible with an overall Hamiltonian description and, thus, inconsistent with the PI formalism.⁷⁵ Nevertheless, it can be used to simulate open quantum systems,^76,77 for example. For a detailed account of these methods, see Ref. 78.

Recently, we proposed a related multiresolution quantum-classical method that, instead of interpolating forces, directly changes the quantum character, degree of quantum delocalization, or “quantumness” of the particles themselves.⁷⁹ In the QM region, the ring polymers are defined as usual, while in the classical region, they collapse to point-like particles, thereby effectively behaving classically. When diffusing between the different regions, the particles change their resolution on the fly. Furthermore, the number of particles in the QM region is not fixed but is allowed to fluctuate. Hence, the scheme can, for example, be used to simulate a quantum grand canonical ensemble efficiently, provided the high- and low-resolution regions are sufficiently large.^76,78,80 The approach is derived in a rigorous fashion from the bottom up and is also compatible with a Hamiltonian description. When restricting the QM part to a small but relevant region in space, the scheme leads to a significant computational speedup. In the example of our previous paper, a liquid parahydrogen system, the calculation of the particle pair interactions was accelerated by a factor of ≈10. Furthermore, the approach can be combined with the previously mentioned methods for efficient PI simulations, such as RPC or GLE thermostating. Therefore, approaching the problem from a different perspective and reducing the number of PI-based interactions in the system, our method is complementary to techniques that make PI computations themselves more efficient.

In our previous paper,⁷⁹ we proposed the general Hamiltonian adaptive quantum-classical scheme, performed a simple validation of the method using a Monte Carlo algorithm to sample the hybrid (HY) Hamiltonian, and demonstrated that the approach can speed up PI-based simulations. In this follow-up article, we show how the method can be extended to perform Hamiltonian multiresolution quantum-classical CMD and RPMD simulations. To this end, we derive an MTS integration protocol suited for the proposed methodology and validate the method by adaptive quantum-classical simulations of liquid water.

Our scheme enables efficient simulations of complex systems by locally taking into account QM delocalization effects. As we demonstrate, this includes the local calculation of approximate quantum statistical and dynamical quantities, such as particle correlation functions, order parameters, and vibrational spectra. We also show in test simulations of liquid water that the technique enables quantum-classical adaptive resolution simulations that are computationally more efficient than full-quantum simulations of the same size. The method can be useful, for example, for interface systems and in simulations of biological objects such as membranes or proteins. Additionally, as previously noted, it allows an efficient implementation of the QM grand canonical ensemble and can, in principle, also be combined with quantum mechanics/molecular mechanics (QM/MM) approaches, in particular those which are based on a similar Hamiltonian interpolation scheme.⁸¹

The paper is organized as follows: In Sec. II, we review the adaptive quantum-classical scheme proposed in our previous work, and in Sec. III we present its implementation within the PIMD framework. In Sec. IV, we discuss how to use the methodology to calculate approximate quantum dynamical quantities in the context of adaptive RPMD and CMD simulations. We describe the details of the simulations we performed for validation in Sec. V and the results are discussed in Sec. VI. In Sec. VII, we summarize the article and conclude.

II. QUANTUM-CLASSICAL PATH INTEGRALS

In quantum statistical mechanics, the canonical partition function of a system of N interacting particles with indices α, momenta ${\hat{p}}_{α}$ ⁠, masses m_α, kinetic energy $\hat{K} = \sum_{α = 1}^{N} {\hat{p}}_{α}^{2} / 2 m_{α}$ ⁠, and potential energy $\hat{V}$ is $Q = Tr [exp (- β \hat{H})]$ with the inverse temperature β = 1/k_bT and the Hamiltonian $\hat{H} = \hat{K} + \hat{V}$ ⁠. Assuming Boltzmann statistics, when performing PI quantization using P imaginary time slices (Trotter number P), the kinetic energy term gives rise to a configurational energy that is equivalent to that of a classical-like ring polymer with P beads that are coupled to their nearest neighbors via harmonic springs (for a detailed derivation see, for example, Ref. 28). This mapping from a quantum particle onto a classical-like polymer ring is exact in the limit P → ∞. In practice, however, well converged results can be obtained for finite values of P, which typically range from 16 to 48 beads for standard PI simulations aimed at treating nuclear quantum effects under ambient conditions.^{10,25,82–86}

The strength of the spring constants between the beads of the ring polymers is $m ω_{P}^{2}$ with $ω_{P} = \sqrt{P} / β ℏ$ ⁠. It is proportional to the temperature as well as the particles’ masses. In other words, the higher the temperature and/or mass, the more collapsed are the ring polymers. The extension of the ring polymers can be interpreted as a measure of the quantum character of the QM particles, with classical behavior corresponding to fully collapsed and therefore localized ring polymers.

The previously proposed method for quantum-classical adaptive resolution simulations⁷⁹ is based on the following observation: In a PI-based formulation of quantum statistical mechanics, the only role of a particle’s mass is to determine the spring constant. It determines the extent to which a particle exhibits its quantum character. The scheme is as follows: For each particle α, we define a resolution parameter $λ_{α} = λ ({\hat{r}}_{α})$ that is a function of the particle’s position ${\hat{r}}_{α}$ ⁠. It smoothly changes from 1 in a spatially predefined QM region to 0 in a classical (CL) region via an intermediate hybrid (HY) transition region (see Fig. 1). Based on this resolution function, we then define a variable mass for particle α as m_α → μ_α(λ_α) = λ_αm_α + (1 − λ_α)M_α. Therefore, in the QM region, μ_α(1) = m_α where m_α is the real mass of the particles while in the CL region, μ_α(0) = M_α $≫$ m_α. The mass M_α must be chosen large enough that particles with μ_α(0) = M_α behave essentially classically (in our previous work, we used M_α = 100 m_α). In this way, particles in the QM region exhibit proper QM behavior, while in the CL region, the ring polymers are forced to collapse to nearly point-like particles that behave essentially classically.

FIG. 1.

View large Download slide

Simulation snapshots from the adaptive quantum-classical simulations. (a) A complete box of quantum-classical water: The blue line shows the resolution function λ switching smoothly from 1 in the QM region to 0 in the CL region. In the QM region, the ring polymers, which correspond to the water atoms, are extended, modeling their quantum mechanical delocalization. In the CL region, they are collapsed to points. (b) Interacting atoms in the CL and QM regions: In the QM region, interactions between different atoms are given as the average over the P time-slices. In the CL region, in contrast, only a single calculation is required for the interaction between a pair of atoms due to the point-like structure of the particles. This alleviates the numerical effort and reproduces normal classical computational efficiency in the CL subsystem.

In addition to variable masses, we also use the Hamiltonian adaptive resolution simulation (H-AdResS) formalism,^87–90 which can be employed to couple different force fields via interpolation of potential energies. A classical H-AdResS Hamiltonian H of a system of N interacting molecules reads

H = K + \sum_{α = 1}^{N} [λ_{α} V_{α}^{1} + (1 - λ_{α}) V_{α}^{0} + V_{α}^{i n t} - Δ H (R_{α})],

(1)

where $K$ is the kinetic energy, α indexes the N particles, and λ_α = λ(R_α) is the previously defined resolution function (when assigning single resolution values λ_α to whole molecules, one typically uses the molecular center of mass R_α as a reference coordinate to determine λ_α). The single-particle potentials $V_{α}^{R e s}$ (with Res = 0, 1) are the sums of all intermolecular potentials acting on particle α, properly normalized so that double counting is avoided. In the present study, the term $V_{α}^{i n t}$ represents all intramolecular interactions, such as bond and angle potentials, which are not subject to interpolation. The term ΔH, referred to as the Free Energy Compensation (FEC),^87,88 is an external field acting in the HY transition region to eliminate the density imbalance that naturally occurs in such dual-resolution systems. Different models of the same physical system exhibit a free energy difference that needs to be neutralized in order to enforce identical thermodynamical and/or structural properties (e.g., density) everywhere in the simulation domain. The FEC levels off these free energy imbalances. H-AdResS has been used mainly to couple atomistic and coarse-grained two-body force fields.^87,88,91,92 In general, however, the potentials can refer to any non-bonded interaction $V_{α}^{R e s}$ ⁠.

We have combined the mass-based quantum-classical interpolation with the H-AdResS scheme in order to enable the use of different force fields within the QM and the CL regions. This can be advantageous, for example, when the CL region’s only role is that of a particle reservoir, in which case a very simple force field can be used in the CL region. On the other hand, when simulating a protein, for example, and only a small part of it needs to treat nuclear quantum effects explicitly, it would be possible to resort to the same force-field everywhere in the system and just “add” the nuclear quantum effects in the relevant region with the present quantum-classical multiresolution scheme.

Combining the quantum-classical mass interpolation and the H-AdResS scheme, we can write the Hamiltonian operator of N interacting Boltzmann particles in three dimensions as⁷⁹

\hat{H} = \sum_{α = 1}^{N} (\frac{1}{2} {\hat{p}}_{α} μ^{- 1} ({\hat{r}}_{α}) {\hat{p}}_{α} + V^{H - A d R e s S} ({\hat{r}}_{α})),

(2)

where $μ^{- 1} ({\hat{r}}_{α})$ is the inverse mass operator. The potential energy term $V^{H - A d R e s S} ({\hat{r}}_{α})$ corresponds to the interpolated H-AdResS potential energy, i.e., the term within the sum in Eq. (1). We have shown that PI quantization then leads to the following approximate expression for the partition function⁷⁹ (the exact nature of the approximation will be discussed shortly):

Q = lim_{P \to \infty} [\prod_{k = 1}^{P} \prod_{α = 1}^{N} \int d r_{α, k} {(\frac{μ_{α, k} P}{2 π β ℏ^{2}})}^{\frac{3}{2}}] e^{- β V_{P}^{μ}}

(3)

with

\begin{matrix} V_{P}^{μ} & = & \sum_{k = 1}^{P} \sum_{α = 1}^{N} \{\frac{μ_{α, k} ω_{P}^{2}}{2} | r_{α, k} - r_{α, k + 1} |^{2}) \\ + (\frac{1}{P} [λ_{α, k} V_{α, k}^{Q M} + (1 - λ_{α, k}) V_{α, k}^{C L} + V_{α, k}^{i n t} - Δ H (r_{α, k})]\}, \end{matrix}

(4)

where α indexes the different particles, k the individual Trotter beads for each of them, and $μ_{α, k} \equiv μ (r_{α, k})$ ⁠. A resolution value $λ_{α, k} = λ (r_{α, k})$ is associated with each bead. We have renamed $V_{α, k}^{0,1}$ to $V_{α, k}^{C L, Q M}$ to emphasize that $V_{α, k}^{C L}$ (⁠ $V_{α, k}^{Q M}$ ⁠) is the intermolecular potential acting in the CL (QM) region. Note that interactions between different ring polymers are computed only between beads having the same Trotter index k, which is why the ring polymers are referred to as “classical-like”. Note, in addition, that the normalization term in Eq. (3) depends on the position of the particles via μ_α,k. To obtain a constant normalization factor, one can transform this position dependent term to a potential energy in $V_{P}^{μ}$ ⁠, as was done in our previous article,⁷⁹ and then treat it as a constant field in the hybrid region. In this work, however, we will deal with it in a different way, which we will discuss in detail later. We would like to highlight that the key point in this section is the introduction of the adaptive mass. Instead of a potential energy term that is based on the spatial interpolation of two classical potentials as in H-AdResS, one could also use other potential energies such as a regular classical potential, an ab initio-based potential, or a QM/MM potential energy.⁸¹ All these choices are equivalent with respect to the central element of a varying effective quantum character of the particles.

The above expression [Eqs. (3) and (4)] is consistent with a rigorous PI quantization if

|\frac{d μ (x)}{d x}| ≪ \frac{2 μ (x)}{Λ_{μ} (x)}

(5)

with $Λ_{μ} (x) \equiv \sqrt{β ℏ^{2} / (P μ (x))}$ ⁠. A derivation of this condition was presented in Ref. 79. This criterion indicates that the interpolation between the QM and the CL parts of the system must be sufficiently smooth such that the mass difference between two neighboring beads on a ring polymer in the HY region is negligible. This requirement can always be satisfied by choosing a sufficiently large HY region. Note, however, that even if Eq. (5) is not fulfilled, Eqs. (3) and (4) correspond to a well-defined quantum-classical simulation protocol.

The ring polymers described by the energy function $V_{P}^{μ}$ [Eq. (4)] are expanded in the region where the mass is small and collapse to nearly classical point-like particles in the large-mass region. Therefore, in the CL region, the interactions between different ring polymers do not need to be computed as an average over the P bead pairs as done in the QM region. Instead, due to their nearly point-like structure, one can use only the centroid with negligible error (see Fig. 1). In this manner, classical computational efficiency is regained in the CL region.

III. QUANTUM-CLASSICAL PATH INTEGRAL MOLECULAR DYNAMICS

In our previous paper,⁷⁹ we validated the scheme introduced above via simulations of liquid parahydrogen using a basic Monte Carlo algorithm to directly sample the phase space defined by $V_{P}^{μ}$ ⁠. A more state-of-the-art approach to the numerical evaluation of PIs is provided by PIMD, in which a thermostated dynamics is generated in phase space to sample the quantum canonical ensemble.^28–34 PIMD is more easily parallelizable compared to Monte Carlo methods and therefore significantly more efficient for typical simulation setups and on multicore computer architectures. Moreover, PIMD is the method of choice for ab initio path integrals.^93,94 In the following, we show how the proposed quantum-classical multiresolution method can be implemented in PIMD as well as CMD and RPMD.

A. Evaluation of the adaptive mass and the resolution function on the centroids

In order to decouple the modes of the cyclic ring polymers from each other, PIMD is typically performed using staging variables^33,35 or, more popularly, normal modes.^33,39 In our case, we cannot transform smoothly into normal mode space because the beads of the individual ring polymers have different masses μ_α,k within the HY region. For typical systems like liquid water at room temperature, however, this mass difference is small as the extension of the ring polymers, measured by the root-mean-square radius of gyration r^g, is short, even in the QM region. This suggests that we associate a single resolution value λ_α and a single adaptive mass value μ_α with each atomic or molecular particle α instead of with every single bead k. λ_α and μ_α can then be determined using the ring polymers’ centroid positions.

Although this corresponds only to a minor modification in Eq. (4), we can ask to what extent the configurational energy $V_{P}^{μ}$ is then still compatible with formal PI quantization. To this end, we first consider the adaptive mass μ_α,k of an individual Trotter bead within the HY region which we can approximate as

μ_{α, k} \approx μ_{α}^{c} + δ x_{α, k} \frac{d μ (x_{α}^{c})}{d x},

(6)

where $x_{α}^{c}$ is the centroid coordinate along the direction of resolution change of the ring polymer α. For a setup where the resolution changes along the x-direction, this would be the centroid’s x-coordinate and for a system where the QM region is spherical and the resolution changes radially, this would be the radial distance from the center. $μ_{α}^{c}$ is the mass function evaluated at the centroid $x_{α}^{c}$ and δx_α,k is the distance along the direction of resolution change between the kth bead of ring α and its centroid $x_{α}^{c}$ ⁠. We have $δ x_{α, k} ≲ r^{g} (x_{α}^{c})$ ⁠, where $r^{g} (x) = \sqrt{β ℏ^{2} / (4 μ (x))} \cdot \sqrt{1 - 1 / P^{2}} \approx \sqrt{β ℏ^{2} / (4 μ (x))}$ is the radius of gyration of a free ring with mass μ(x) for large P. Hence, we can approximate $μ_{α, k} \approx μ_{α}^{c}$ if

|\frac{d μ (x_{α}^{c})}{d x}| ≪ \frac{μ_{α}^{c}}{r^{g} (x_{α}^{c})}

(7)

or simply

|\frac{d μ (x)}{d x}| ≪ \frac{μ (x)}{r^{g} (x)},

(8)

for general x. Here, we used the free ring-polymer radius of gyration. However, for the ring polymers in typical systems, the average radius of gyration differs only slightly from the free-particle radius.

Equation (8) is a slightly stronger criterion than that in Eq. (5). This makes sense since Eq. (5) essentially provides the condition under which the mass can be considered as constant between two neighboring beads, while the new criterion, Eq. (8), gives the condition for treating the mass as constant over an entire ring polymer. For liquid water at room temperature, Eq. (8) is satisfied by a hybrid region wider than ≈1 nm. An even smaller hybrid region would not be desirable anyway since the interaction cutoffs of typical interaction potentials are also of the order ≈1 nm.

Next, we consider the resolution λ itself, which we also seek to treat as constant over a whole ring polymer so that it can be approximated as a function of the centroid position only. On the one hand, λ varies between 0 and 1 and changes most steeply in the center of the HY region (see Fig. 1). On the other hand, an upper bound on the extension of the ring polymers is provided by the radius of gyration of the ring polymers in the QM region, which will be denoted as $r_{QM}^{g}$ ⁠. Therefore, if the change in λ, corresponding to its gradient in the center of the HY region, over a distance $r_{QM}^{g}$ is much smaller than ≈1, we can approximate $λ_{α, k} \approx λ_{α}^{c}$ everywhere (as for the adaptive mass before, $λ_{α}^{c}$ denotes the resolution function of ring α evaluated at its centroid $x_{α}^{c}$ ⁠). This corresponds to

|(\frac{d λ (x)}{d x} {|)}_{x_{=} d_{HY} / 2})| r_{QM}^{g} ≪ 1,

(9)

where $d λ (x) / d x |_{x = d_{HY} / 2}$ denotes the gradient of λ in the center of the HY region and d_HY denotes the width of the HY region. The criterion in Eq. (9) can be easily fulfilled for typical systems such as water at room temperature, with a HY region of width d_HY ≈ 1 nm.

To summarize, in the following, we will treat both the mass and the resolution as a constant over entire ring polymers and write, for simplicity, μ_α,k(x) → μ_α(x) and λ_α,k(x) → λ_α(x), where we assume that the mass and the resolution functions for an atom α have been evaluated using its centroid coordinate $x_{α}^{c}$ and that the resulting parameters have been assigned to all beads belonging to ring α. Then we can write the partition function as

Q = lim_{P \to \infty} [\prod_{k = 1}^{P} \prod_{α = 1}^{N} \int d r_{α, k} {(\frac{μ_{α} P}{2 π β ℏ^{2}})}^{\frac{3}{2}}] e^{- β {\tilde{V}}_{P}^{μ}}

(10)

with

\begin{matrix} {\tilde{V}}_{P}^{μ} & = & \sum_{k = 1}^{P} \sum_{α = 1}^{N} \{\frac{μ_{α} ω_{P}^{2}}{2} | r_{α, k} - r_{α, k + 1} |^{2}) \\ + (\frac{1}{P} [λ_{α} V_{α, k}^{Q M} + (1 - λ_{α}) V_{α, k}^{C L} + V_{α, k}^{i n t} - Δ H (r_{α}^{c})]\}, \end{matrix}

(11)

where $r_{α}^{c}$ denotes the centroid of ring polymer α. Note that the FEC is now also applied at the single-atom, i.e., centroid level. The criteria in Eqs. (8) and (9) quantify the degree to which this partition function is still compatible with a formal, bottom-up PI quantization. The inequalities can always be fulfilled by choosing a sufficiently wide hybrid region. However, even if they are not met, the final partition function, Eqs. (10) and (11), still represents a well-defined Hamiltonian multiresolution quantum-classical simulation scheme.

B. Introducing normal modes

Now that the different beads of each ring polymer α all have the same adaptive mass μ_α, we can proceed with transforming the Cartesian coordinates into normal modes u_α,k via

u_{α, k} = \sum_{j = 1}^{P} r_{α, j} C_{j k},

(12)

where, for even P, the orthogonal transformation matrix is⁶³

C_{j k} = \{\begin{matrix} \sqrt{1 / P} & i f k = 1 \\ \sqrt{2 / P} c o s (2 π j k / P) & i f 2 \leq k \leq P / 2 \\ \sqrt{1 / P} {(- 1)}^{j} & i f k = P / 2 + 1 \\ \sqrt{2 / P} s i n (2 π j k / P) & i f P / 2 + 2 \leq k \leq P \end{matrix})

(13)

such that for a given ring polymer at position x,

μ_{α} (x) \sum_{k = 1}^{P} | r_{α, k} - r_{α, k + 1} |^{2} = μ_{α} (x) \sum_{k = 1}^{P} ξ_{k} u_{α, k}^{2}

(14)

with

ξ_{k} = 4 {sin}^{2} (\frac{(k - 1) π}{P}) .

(15)

In the normal mode representation, the centroid of a ring polymer α is given by the rescaled first mode coordinate

r_{α}^{c} = \frac{1}{P} \sum_{k = 1}^{P} r_{α, k} = \frac{1}{\sqrt{P}} u_{α, 1} .

(16)

Therefore, the adaptive mass μ_α and the resolution function λ_α, being evaluated at the centroids, are functions of the first mode only, i.e., $μ_{α} (r_{α}^{c}) = μ_{α} (1 / \sqrt{P} u_{α, 1})$ and $λ_{α} (r_{α}^{c}) = λ_{α} (1 / \sqrt{P} u_{α, 1})$ ⁠. To lighten the notation, though, we will drop the $1 / \sqrt{P}$ factor and write simply μ_α = μ_α(u_α,1) and λ_α = λ_α(u_α,1). We then obtain the following partition function:

Q = lim_{P \to \infty} [\prod_{k = 1}^{P} \prod_{α = 1}^{N} \int d u_{α, k} {(\frac{μ_{α} (u_{α, 1}) P}{2 π β ℏ^{2}})}^{\frac{3}{2}}] e^{- β {\tilde{V}}_{P}^{μ}}

(17)

with

\begin{matrix} {\tilde{V}}_{P}^{μ} & = & \sum_{k = 1}^{P} \sum_{α = 1}^{N} \{\frac{1}{2} ν_{α, k} (u_{α, 1}) ω_{P}^{2} u_{α, k}^{2} + \frac{1}{P} [λ_{α} (u_{α, 1}) V_{α, k}^{Q M} (u))) \\ ((+ (1 - λ_{α} (u_{α, 1})) V_{α, k}^{C L} (u) + V_{α, k}^{i n t} (u) - Δ H (u_{α, 1})]\} \end{matrix}

(18)

and the rescaled adaptive mass ν_α,k(u_α,1) = μ_α(u_α,1) ξ_k. For the centroid, i.e., k = 1, we have ν_α,k(u_α,1) = 0. In Eqs. (17) and (18), we have explicitly indicated the dependencies of the different terms on the normal modes. The notation u without any indices is a shorthand for the compound set of all coordinates u_α,k.

C. Introducing momenta

In PIMD, one usually recasts the prefactor of the partition function as Gaussian integrals over a set of variables that can be interpreted as momenta conjugate to the coordinates u_α,k.²⁸ The energy term in the exponential can then be interpreted as a classical Hamiltonian and sampled via thermostated molecular dynamics. In our case, however, the prefactor is position dependent and, therefore, we have different options to proceed.

1. Constant kinetic masses (CKMs)

As was done in our previous article,⁷⁹ we can write

\begin{matrix} {(\frac{μ_{α} (u_{α, 1}) P}{2 π β ℏ^{2}})}^{\frac{3}{2}} & = & {(\frac{\tilde{m} P}{2 π β ℏ^{2}})}^{\frac{3}{2}} \\ \times exp \{- β (- \frac{3}{2 β} log (\frac{μ_{α} (u_{α, 1})}{\tilde{m}}))\}, \end{matrix}

(19)

where we used an arbitrary mass $\tilde{m}$ as the reference mass scale. Then we can pull the term $- (3 / 2 β) log (μ_{α} / \tilde{m})$ into ${\tilde{V}}_{P}^{μ}$ as an external field in the HY region and use the now constant prefactor to introduce a set of momenta via rephrasing the prefactor as Gaussian integrals. This yields

Q = lim_{P \to \infty} [\prod_{k = 1}^{P} \prod_{α = 1}^{N} \int d u_{α, k} \int d p_{α, k}] e^{- β H_{P}^{c k m}}

(20)

with the Hamiltonian

\begin{matrix} H_{P}^{c k m} & = & \sum_{k = 1}^{P} \sum_{α = 1}^{N} \{\frac{p_{α, k}^{2}}{2 m_{α, k}^{'}} + \frac{1}{2} ν_{α, k} ω_{P}^{2} u_{α, k}^{2} - \frac{3}{2 β} log \frac{μ_{α} (u_{α, 1})}{\tilde{m}}) \\ (+ \frac{1}{P} [λ_{α} (u_{α, 1}) V_{α, k}^{Q M} (u) + (1 - λ_{α} (u_{α, 1})) V_{α, k}^{C L} (u))) \\ ((+ V_{α, k}^{i n t} (u) - Δ H (u_{α, 1})]\} . \end{matrix}

(21)

As the additional logarithmic term only acts as an external field in the HY region, it can be exactly removed via the FEC function ΔH. We will denote the fictitious masses $m_{α, k}^{'}$ in the following as “kinetic” masses in contrast to ν_α,k, which we will refer to as spring masses. In principle, the set of $m_{α, k}^{'}$ can be chosen freely, as their rescaling does not affect thermodynamic averages.²⁸

$H_{P}^{c k m}$ in Eq. (21) defines a classical Hamiltonian system composed of ring polymers representing the delocalized atoms. In the QM region where λ_α = 1 and μ_α = m, the ring polymers are extended and the regular quantum mechanical behaviour is recovered. In the CL region where λ_α = 0 and μ_α = M, the ring polymers are collapsed to essentially point-like particles, thereby reproducing classical mechanics. The Hamiltonian gives rise to regular equations of motion that can be integrated by a symplectic integrator such as the velocity Verlet algorithm,²⁸ with the possibility of employing multiple time-stepping.

However, let us take a closer look at the different masses in the system. While the spring masses ν_α,k change between the CL and QM subregions of the system, the kinetic masses $m_{α, k}^{'}$ do not. We choose $m_{α, k}^{'} = m_{α} / P$ with m_α being the real mass of atom α since this corresponds to a realistic bead-wise approximate quantum dynamical behavior in the QM region similar to RPMD [in RPMD, one usually chooses $m_{α, k}^{'} = m_{α}$ without 1/P but rescales the potential energy terms by P and runs the simulation at a P-times higher temperature.^{45,46,48–51} Here, the factor of 1/P in the mass is equivalent to this procedure, as we perform the simulations at the actual temperature and use a Hamiltonian, Eq. (21), without rescaling potential energies]. Therefore, in the QM region, the modes oscillate with vibration frequencies $ω_{P} \sqrt{ξ_{k} P}$ ⁠. In the CL region, however, where ν_α,k is significantly larger than in the QM subsystem, the modes oscillate faster than in the QM region by a factor of $\sqrt{M_{α} / m_{α}}$ ⁠. For the case of M_α = 100 m_α, this results in 10 times higher frequencies. This would require a 10 times smaller time step in the integration algorithm compared to a normal quantum simulation or compared to what would be required in the QM subregion. Although this poses no fundamental hurdle, it may slow down the simulations notably.

2. Adaptive kinetic masses (AKMs)

The previous observation suggests an alternative approach: We can also directly recast the prefactor as a Gaussian integral, which includes the position dependent mass μ_α,

\begin{matrix} {(\frac{μ_{α} (u_{α, 1}) P}{2 π β ℏ^{2}})}^{\frac{3}{2}} = {(\frac{P^{2}}{4 π^{2} ℏ^{2}})}^{\frac{3}{2}} \int d p_{α, k} exp \{- β P \frac{p_{α, k}^{2}}{2 μ_{α} (u_{α, 1})}\} . \end{matrix}

(22)

In this way, we can introduce a kinetic energy term that has adaptive kinetic masses. This leads to the construction of a Hamiltonian in which both the spring and the kinetic masses vary in the same fashion such that the modes oscillate with the same frequencies everywhere in the quantum-classical adaptive resolution setup.

Specifically, we propose the following: Overall, we have N × P prefactors of the form

{(\frac{μ_{α} (u_{α, 1}) P}{2 π β ℏ^{2}})}^{\frac{3}{2}},

(23)

one for each atom and mode. For all higher modes with k > 1, we transform the prefactors according to Eq. (22) and introduce the momentum terms in the Hamiltonian with a variable mass in the denominator. The remaining N prefactors are then treated via Eq. (19), and the kinetic masses for the centroid modes, which are not associated with springs since ν_α,1 = 0, are chosen constant. We then obtain

Q = lim_{P \to \infty} [\prod_{k = 1}^{P} \prod_{α = 1}^{N} \int d u_{α, k} \int d p_{α, k}] e^{- β H_{P}^{a k m}}

(24)

with the Hamiltonian

\begin{matrix} H_{P}^{a k m} & = & \sum_{k = 1}^{P} \sum_{α = 1}^{N} \{\frac{p_{α, k}^{2}}{2 ν_{α, k}^{'} (u_{α, 1})} + \frac{1}{2} ν_{α, k} ω_{P}^{2} u_{α, k}^{2}) \\ (- \frac{3}{2 P β} log \frac{μ_{α} (u_{α, 1})}{\tilde{m}} + \frac{1}{P} [λ_{α} (u_{α, 1}) V_{α, k}^{Q M} (u))) \\ + (((1 - λ_{α} (u_{α, 1})) V_{α, k}^{C L} (u) + V_{α, k}^{i n t} (u) - Δ H (u_{α, 1})]\}, \end{matrix}

(25)

where $ν_{α, k}^{'}$ is the kinetic mass of bead k of atom α. Note that a 1/P factor appears in front of the logarithmic term in $H_{P}^{a k m}$ because the term still appears in the sum over all P, although we obtain the logarithmic term only for the centroid modes. Choosing appropriate prefactors, the parameters $ν_{α, k}^{'}$ are

ν_{α, 1}^{'} = m_{α} / P, k = 1, c e n t r o i d m o d e,

(26)

ν_{α, k}^{'} = {\tilde{ν}}_{α} = μ_{α} / P, k > 1, h i g h e r m o d e s .

(27)

Since the kinetic masses for the higher modes are all equal, we introduced a new abbreviation, ${\tilde{ν}}_{α}$ ⁠, for them without the index k. We choose a factor 1/P to ensure that the approximate quantum dynamical time evolution of the centroids proceeds on the real time scale, that is, the same as in corresponding classical Newtonian dynamics. As already pointed out, this choice is equivalent to the temperature rescaling often done in ring polymer molecular dynamics,^{45,46,48–51} which we do not perform here. A further rescaling of the kinetic masses would be allowed when sampling only canonical averages.^28,33

We can interpret the choice of the kinetic masses in the following way: While moving from the QM to the CL via the HY region, the spring constants of the higher modes become stronger and the ring polymers collapse. Simultaneously, however, these higher modes become heavier such that they do not vibrate faster in the CL region than in the QM region, despite the stiffer springs. Their oscillation frequencies are the same everywhere in the system such that their configurations can be sampled efficiently throughout the whole system with the same time step. The centroid modes do not undergo such oscillations, as they represent only the displacements of the ring polymers as a whole. Hence, their masses do not need to change across the transition from the QM to the CL part of the system and, therefore, they are chosen to be the real masses.

In the following, we will refer to the constant kinetic mass (CKM) approach, defined by the Hamiltonian $H_{P}^{c k m}$ in Eq. (21), as the CKM approach and to the adaptive kinetic mass (AKM) scheme, defined by the Hamiltonian $H_{P}^{a k m}$ in Eq. (25), as the AKM approach.

D. Equations of motion

We will not discuss the equations of motion obtained in the CKM approach, as they resemble a regular structure that can be integrated, for example, by a regular velocity Verlet algorithm.²⁸ Instead, we focus on the AKM scheme, from which the CKM approach can be obtained as a special case.

In the following, we will assume that the logarithmic term in the Hamiltonian in Eq. (25) has been exactly canceled by an appropriately chosen FEC function ΔH(λ) and therefore omit it to simplify the notation.

The equations of motion then read as follows:

\begin{matrix} {\dot{u}}_{α, 1} = \frac{p_{α, 1}}{ν_{α, 1}^{'} (u_{α, 1})} = P \frac{p_{α, 1}}{m_{α}}, k = 1, c e n t r o i d m o d e, \end{matrix}

(28)

\begin{matrix} {\dot{u}}_{α, k} = \frac{p_{α, k}}{{\tilde{ν}}_{α} (u_{α, 1})} = P \frac{p_{α, k}}{μ_{α} (u_{α, 1})}, k > 1, h i g h e r m o d e s, \end{matrix}

(29)

\begin{matrix} {\dot{p}}_{α, 1} & = & F_{α, 1}^{H - A d R e s S} (u) \\ - \frac{1}{P} [\sum_{i = 1}^{P} \{V_{α, i}^{Q M} (u) - V_{α, i}^{C L} (u)\}] \nabla_{u_{α, 1}} λ (u_{α, 1}) (F_{1}^{d r i f t}) \\ - [\sum_{i = 2}^{P} \frac{1}{2} \frac{(M_{α} - m_{α})}{P {\tilde{ν}}_{α} {(u_{α, 1})}^{2}} p_{α, i}^{2}] \nabla_{u_{α, 1}} λ (u_{α, 1}) (F_{2}^{d r i f t}) \\ - [\sum_{i = 2}^{P} \frac{1}{2} ξ_{i} (m_{α} - M_{α}) ω_{P}^{2} u_{α, i}^{2}] \nabla_{u_{α, 1}} λ (u_{α, 1}) (F_{3}^{d r i f t}) \\ + \nabla_{u_{α, 1}} Δ H (u_{α, 1}), \\ k = 1, c e n t r o i d m o d e, \end{matrix}

(30)

\begin{matrix} {\dot{p}}_{α, k} & = & - ν_{α, k} (u_{α, 1}) ω_{P}^{2} u_{α, k} + F_{α, k}^{H - A d R e s S} (u), \\ k > 1, h i g h e r m o d e s, \end{matrix}

(31)

where

\begin{matrix} F_{α, k}^{H - A d R e s S} (u) & = & - \frac{1}{P} \sum_{j = 1}^{P} \sum_{β = 1}^{N} [λ_{β} \nabla_{u_{α, k}} V_{β, j}^{Q M} (u)) \\ + ((1 - λ_{β}) \nabla_{u_{α, k}} V_{β, j}^{C L} (u) + \nabla_{u_{α, k}} V_{β, j}^{i n t} (u)] . \end{matrix}

(32)

The terms in lines 2-5 in Eq. (30) stem from the application of the derivative on the position dependent resolution function. The terms $F_{i}^{d r i f t}$ are undesired forces that act only in the hybrid region, can lead to thermodynamic imbalances in the system, and, for example, artificially push particles from one subregion of the system to the other. In accordance with earlier studies using the H-AdResS scheme, we will refer to these forces as drift forces.⁸⁷ They need to be compensated, which can be achieved via the FEC, line 5 in Eq. (30). In fact, the latter is typically constructed to cancel their average effect.^87,88,90

The drift force $F_{1}^{d r i f t}$ comes from the potential energy interpolation and would not be present if we changed only the masses of the atoms but not the force field. $F_{2}^{d r i f t}$ is a result of choosing the kinetic masses of the higher modes to be adaptive and would be absent in the CKM approach. Finally, $F_{3}^{d r i f t}$ corresponds to the adaptive spring masses. As the resolution λ_α and the adaptive mass μ_α depend only on the centroid positions of the ring polymers, drift forces only occur in the equations of motion for the centroid. Therefore, the internal motion of the ring polymers is not disturbed by any drift forces and can smoothly change from quantum to classical and vice versa when the atoms move through the system. Only the translation of the ring polymers is affected by drift forces, which can be corrected via the FEC. Note that the sums in $F_{2}^{d r i f t}$ and $F_{3}^{d r i f t}$ only run over non-centroid modes since the centroid mode neither appears in the spring constant term in the Hamiltonian nor it is associated with a variable mass in the kinetic energy onto which the position derivative could act.

An interesting observation is that the drift forces $F_{2}^{d r i f t}$ and $F_{3}^{d r i f t}$ have strictly opposite signs. Furthermore $F_{2}^{d r i f t}$ is proportional to $p_{α, i}^{2}$ ⁠, which is related to the kinetic energy of mode i, while $F_{3}^{d r i f t}$ is proportional to $u_{α, i}^{2}$ ⁠, which corresponds to the potential energy in mode i. Since the higher modes behave largely like harmonic oscillators, the energy fluctuates back and forth between the potential and kinetic energy. Therefore, $F_{2}^{d r i f t}$ and $F_{3}^{d r i f t}$ partially cancel out each other. This effect would not be present in the CKM approach, where $F_{2}^{d r i f t}$ is missing. Therefore, we expect the overall drift forces in the CKM approach to be stronger than in the AKM method and this is indeed observed in our simulations.

E. Integration

To devise a suitable integration scheme for the equations of motion, Eqs. (28)–(31), we use the Liouville operator formalism. To simplify the notation, in the following, we will drop the atom index α, as the Liouville operators for different ring polymers, commute, and do not require further discussion.

The Liouville operator iL can be written as

i L = i L^{(1)} + i L^{(2)},

(33)

where iL⁽¹⁾ propagates positions and iL⁽²⁾ propagates momenta. We further decompose

i L^{(1)} = i L_{1}^{(1)} + i L_{{k}}^{(1)},

(34)

i L^{(2)} = i L_{1}^{(2)} + i L_{{k}}^{(2)},

(35)

where $i L_{{k}}^{(i)} = \sum_{k > 1}^{P} i L_{k}^{(i)}$ ⁠, operators $i L_{1}^{(i)}$ propagate the first mode, and $i L_{k}^{(i)}$ act on modes $k$ with $k > 1$ ⁠. They are

i L_{1}^{(1)} = P \frac{p_{1}}{m} \frac{\partial}{\partial u_{1}},

(36)

i L_{k}^{(1)} = \frac{p_{k}}{\tilde{ν} (u_{1})} \frac{\partial}{\partial u_{k}},

(37)

i L_{1}^{(2)} = (F_{1} (u) + η (u_{1}) \sum_{k = 2}^{P} p_{k}^{2}) \frac{\partial}{\partial p_{1}},

(38)

i L_{k}^{(2)} = F_{k} (u) \frac{\partial}{\partial p_{k}},

(39)

with

\begin{matrix} F_{1} (u) & = & F_{1}^{H - A d R e s S} (u) \\ - \frac{1}{P} [\sum_{i = 1}^{P} \{V_{i}^{Q M} (u) - V_{i}^{C L} (u)\}] \nabla_{u_{1}} λ (u_{1}) \\ - [\sum_{i = 2}^{P} \frac{1}{2} ξ_{i} (m - M) ω_{P}^{2} u_{i}^{2}] \nabla_{u_{1}} λ (u_{1}) + \nabla_{u_{1}} Δ H (u_{1}), \end{matrix}

(40)

and

F_{k} (u) = - ν_{k} (u_{1}) ω_{P}^{2} u_{k} + F_{k}^{H - A d R e s S} (u),

(41)

η (u_{1}) = - \frac{1}{2} \frac{(M - m)}{P \tilde{ν} {(u_{1})}^{2}} \nabla_{u_{1}} λ (u_{1}) .

(42)

The vector notation denotes that each Liouville operator in Eqs. (36)–(39) represents a set of three operators for each direction, which commute and can therefore be applied in arbitrary order.

In the first step, we decompose the classical propagator exp{iLt} as

e^{i L t} = e^{i L^{(1)} t + i L^{(2)} t} = lim_{M \to \infty} {(e^{i L^{(2)} \frac{t}{2 M}} e^{i L^{(1)} \frac{t}{M}} e^{i L^{(2)} \frac{t}{2 M}})}^{M}

(43)

using the symmetric Trotter theorem.^95,96 Defining the time step Δt = t/M, this yields the velocity Verlet integrator

e^{i L Δ t} \approx e^{i L^{(2)} \frac{Δ t}{2}} e^{i L^{(1)} Δ t} e^{i L^{(2)} \frac{Δ t}{2}} .

(44)

However, recalling the definitions of iL⁽¹⁾ and iL⁽²⁾ we recognize that their constituents $i L_{1}^{(1)}$ and $i L_{{k}}^{(1)}$ as well as $i L_{1}^{(2)}$ and $i L_{{k}}^{(2)}$ do not commute. Therefore, using again the Trotter theorem, we decompose each of the propagators in Eq. (44) further to

\begin{matrix} e^{i L Δ t} & \approx & e^{i L^{(2)} \frac{Δ t}{2}} e^{i L^{(1)} Δ t} e^{i L^{(2)} \frac{Δ t}{2}} \\ \approx & (e^{i L_{{k}}^{(2)} \frac{Δ t}{4}} e^{i L_{1}^{(2)} \frac{Δ t}{2}} e^{i L_{{k}}^{(2)} \frac{Δ t}{4}}) (e^{i L_{1}^{(1)} \frac{Δ t}{2}} e^{i L_{{k}}^{(1)} Δ t} e^{i L_{1}^{(1)} \frac{Δ t}{2}}) \\ \times (e^{i L_{{k}}^{(2)} \frac{Δ t}{4}} e^{i L_{1}^{(2)} \frac{Δ t}{2}} e^{i L_{{k}}^{(2)} \frac{Δ t}{4}}), \end{matrix}

(45)

which is correct up to second order in Δt and is suited for the integration of the system’s equations of motion. Note that for the CKM approach, the adaptive mass $\tilde{ν} (u_{1})$ in $i L_{k}^{(1)}$ ⁠, Eq. (37), would be constant and also the $p_{k}^{2}$ term in $i L_{1}^{(2)}$ ⁠, Eq. (38), would be missing. Therefore, no second decomposition level would be required and we could stick with velocity Verlet.

The interpretation of the integration scheme in Eq. (45) is straightforward. The first and the last term in brackets (⋯) correspond to the propagation of momenta, while the term in the middle propagates coordinates. The next decomposition level tells us how to update the momenta and the coordinates. The momenta are integrated in the following way: We first propagate all higher mode momenta by a quarter step, then we update the $p_{k}^{2}$ term in $i L_{1}^{(2)}$ using these new momenta and propagate the first mode’s momentum by a half step. Next, we perform the other quarter step for the higher modes. To get the full time step, the procedure is repeated after the position update in the center of the scheme. The coordinates are updated as follows: We first propagate the first mode by a half step using $i L_{1}^{(1)}$ ⁠. Then we update the adaptive masses in the Liouville operators $i L_{k}^{(1)}$ and propagate the higher modes by a full step. Finally, we integrate the first mode by another half step. The scheme requires little additional computational overhead compared to a regular velocity Verlet scheme. The number of additional operations scales only linearly with the number of particles, and the force computation, usually the numerically most demanding part of a simulation, has to be performed as usual only once after all coordinates are fully propagated.

Finally, we want to address the symplecticity of the new integrator. The Hamiltonian $H_{P}^{a k m}$ ⁠, Eq. (25), is not trivially separable into two parts, one depending only on coordinates and one containing only momenta. However, it has no term in which both the momentum and the corresponding conjugate coordinate of the same mode appear together. This is a result of our choice of the adaptive masses: The higher mode masses are position dependent, but they do not depend on their own mode coordinates, i.e., they only depend on the centroid coordinate. However, the centroid itself, which determines the resolution and the masses of the ring polymers, is associated with a constant mass. Therefore, the Hamiltonian $H_{P}^{a k m}$ and the corresponding equations of motion still define a symplectic structure. As a consequence, the integration scheme in Eq. (45) is also symplectic, as it is constructed in a rigorous fashion using the Liouville operator associated with the equations of motion derived from $H_{P}^{a k m}$ ⁠. This can also be understood considering the Liouville operators themselves: The operators $i L_{k}^{(1)}$ ⁠, Eq. (37), which propagate the higher mode coordinates u_k do have an additional position dependence but it is only on the centroid mode. Hence, they can be applied as usual in a well-defined manner. Similarly, the operator $i L_{1}^{(2)}$ ⁠, Eq. (38), propagating the momentum of the centroid mode has an additional dependence on momenta, but it is only on the higher mode momenta. Consequently, the determinant J of the time evolution matrix is 1. The symplecticity has the practical advantage that we are able to derive an energy conserving integrator, which, in our case, is exact up to second order in time, similar to standard velocity Verlet.

It is worth pointing out that the previous observation is in contrast with earlier studies using adaptive masses.^97,98 There, both momenta and the corresponding conjugate coordinates appear together in the same terms in the Hamiltonian. Hence, in those cases, Liouville’s theorem no longer holds, the Liouville operator formalism breaks down, and symplecticity is lost.

F. Multiple time-stepping

In typical complex soft matter systems, non-bonded interactions as well as bonds, angles, and dihedrals generate motion on different time scales. In PIMD, we additionally have the springs between the beads of the ring polymers onto which the quantum particles are mapped. If the kinetic masses for the higher modes are small, they vibrate strongly, which requires a small integration time step. When only statistical averages are sought, the kinetic masses can be chosen freely, for example, such that all higher modes vibrate with the same frequency. When calculating approximate quantum dynamical quantities, however, the kinetic mode masses must either correspond to the real ones, as in RPMD, or must be significantly decreased, as in adiabatic CMD. This leads to an internal ring polymer dynamics which is significantly faster than the motion due to typical interatomic non-bonded or bonded potentials. Furthermore, in the CKM approach, the modes’ oscillation frequencies are increased in the CL region, as the kinetic mass will be small there compared to the increased spring mass. This strongly motivates the introduction of multiple time-stepping into our integrator.

We employ the RESPA scheme⁶⁶ and decompose the force computation into three parts: one for non-bonded forces, a second for the bonds, and a third for the internal ring polymer motion. The first drift term, $F_{1}^{d r i f t}$ ⁠, depends only on the energies associated with the non-bonded potentials and is therefore evaluated together with the rest of the non-bonded forces. The second and the third drift terms $F_{2}^{d r i f t}$ and $F_{3}^{d r i f t}$ ⁠, however, depend directly on the motion of the higher modes and therefore need to be evaluated together with them. Hence, we define

i L_{1}^{(1)} = P \frac{p_{1}}{m} \frac{\partial}{\partial u_{1}},

(46)

i L_{k}^{(1)} = \frac{p_{k}}{\tilde{ν} (u_{1})} \frac{\partial}{\partial u_{k}},

(47)

i L_{1}^{(2)} = (F_{1}^{m o d e} (u) + η (u_{1}) \sum_{k = 2}^{P} p_{k}^{2}) \frac{\partial}{\partial p_{1}},

(48)

i L_{k}^{(2)} = F_{k}^{m o d e} (u) \frac{\partial}{\partial p_{k}},

(49)

i L_{q}^{(3)} = F_{q}^{i n t} (u) \frac{\partial}{\partial p_{q}},

(50)

i L_{q}^{(4)} = F_{q}^{n b} (u) \frac{\partial}{\partial p_{q}},

(51)

with

F_{1}^{m o d e} (u) = - [\sum_{i = 2}^{P} \frac{1}{2} ξ_{i} (m - M) ω_{P}^{2} u_{i}^{2}] \nabla_{u_{1}} λ (u_{1}),

(52)

F_{k}^{m o d e} (u) = - ν_{k} (u_{1}) ω_{P}^{2} u_{k},

(53)

F_{q}^{i n t} (u) = - \frac{1}{P} \sum_{j = 1}^{P} \sum_{β = 1}^{N} \nabla_{u_{q}} V_{β, j}^{i n t} (u),

(54)

F_{q}^{n b} (u) = - \frac{1}{P} \sum_{j = 1}^{P} \sum_{β = 1}^{N} [λ_{β} \nabla_{u_{q}} V_{β, j}^{Q M} (u)) + ((1 - λ_{β}) \nabla_{u_{q}} V_{β, j}^{C L} (u)]

- \frac{1}{P} [\sum_{i = 1}^{P} \{V_{i}^{Q M} (u) - V_{i}^{C L} (u)\}] \nabla_{u_{q}} λ (u_{1})

+ \nabla_{u_{q}} Δ H (u_{1}) .

(55)

Note that the Liouville operators $i L_{q}^{(3)}$ and $i L_{q}^{(4)}$ do not need to be split into centroid and higher terms, as these commute in this case. Hence, for the sake of brevity, we have subsumed both parts and changed to the index q, which includes all 1 ≤ q ≤ P. We then define Liouville operators acting on all modes q as $i L_{{q}}^{(i)} = \sum_{q = 1}^{P} {i L}_{q}^{(i)}$ ⁠. Finally, we obtain the following RESPA multiple time-stepping scheme:

\begin{matrix} e^{i L Δ t} & \approx & e^{i L_{{q}}^{(4)} \frac{Δ t}{2}} \{e^{i L_{{q}}^{(3)} \frac{δ t}{2}} [(e^{i L_{{k}}^{(2)} \frac{d t}{4}} e^{i L_{1}^{(2)} \frac{d t}{2}} e^{i L_{{k}}^{(2)} \frac{d t}{4}}))) \\ \times (e^{i L_{1}^{(1)} \frac{d t}{2}} e^{i L_{{k}}^{(1)} d t} e^{i L_{1}^{(1)} \frac{d t}{2}}) \\ {({(\times (e^{i L_{{k}}^{(2)} \frac{d t}{4}} e^{i L_{1}^{(2)} \frac{d t}{2}} e^{i L_{{k}}^{(2)} \frac{d t}{4}})]}^{n} e^{i L_{{q}}^{(3)} \frac{δ t}{2}}\}}^{N} e^{i L_{{q}}^{(4)} \frac{Δ t}{2}}, \end{matrix}

(56)

where Δt = N · δt = N · n · dt. The internal ring vibrations as well as the drift terms depending on these higher ring modes are integrated with the shortest time step dt. The intramolecular bonds and angles are integrated with a distinct, intermediate time step δt, and the intermolecular non-bonded interactions as well as the corresponding drift terms are integrated with the largest time step Δt. The whole integration may be carried out in normal mode space, although in practice, the interatomic forces are computed in real space and then transformed into mode space.

In conventional normal mode PIMD, the integration of the higher modes is typically carried out analytically via the exact propagator, which makes the integrator very stable.⁶³ In our case, this is not possible in a straightforward manner due to the adaptive masses in the spring constant and kinetic energy terms, which lead to drift forces and require an additional decomposition of the integrator. As will be demonstrated in our simulations, we nevertheless obtain a stable integration scheme.

G. Langevin thermostating

To generate a canonical ensemble, we need to couple the system to a thermostat. We resort to a Langevin thermostat, as Langevin equation-based frameworks have been shown to be favorable in PIMD and RPMD simulations and can be used to optimize sampling efficiency.^62–65 As the focus of this work is not advanced thermostating, however, we use a simple white noise Langevin thermostat without memory instead of, for example, a GLE approach. The implementation follows the BAOAB method by Leimkuhler and Matthews,⁹⁹ which provides high configurational sampling accuracy (a detailed comparison of different integrator splittings for the implementation of Langevin thermostats in the context of PIMD can be found in Ref. 100). Within the proposed multiple time-stepping scheme, this yields

\begin{matrix} e^{i L Δ t} & \approx & e^{i L_{{q}}^{(4)} \frac{Δ t}{2}} \{e^{i L_{{q}}^{(3)} \frac{δ t}{2}} [(e^{i L_{{k}}^{(2)} \frac{d t}{4}} e^{i L_{1}^{(2)} \frac{d t}{2}} e^{i L_{{k}}^{(2)} \frac{d t}{4}}))) \\ \times (e^{i L_{1}^{(1)} \frac{d t}{4}} e^{i L_{{k}}^{(1)} \frac{d t}{2}} e^{i L_{1}^{(1)} \frac{d t}{4}}) e^{i L_{{q}}^{L a n g e v i n} d t} (e^{i L_{1}^{(1)} \frac{d t}{4}} e^{i L_{{k}}^{(1)} \frac{d t}{2}} e^{i L_{1}^{(1)} \frac{d t}{4}}) \\ {({(\times (e^{i L_{{k}}^{(2)} \frac{d t}{4}} e^{i L_{1}^{(2)} \frac{d t}{2}} e^{i L_{{k}}^{(2)} \frac{d t}{4}})]}^{n} e^{i L_{{q}}^{(3)} \frac{δ t}{2}}\}}^{N} e^{i L_{{q}}^{(4)} \frac{Δ t}{2}}, \end{matrix}

(57)

where $i L_{{q}}^{L a n g e v i n} = \sum_{q = 1}^{P} i L_{q}^{L a n g e v i n}$ and with the action of the Langevin Liouville operator $i L_{{q}}^{L a n g e v i n}$ on mode q,

e^{i L^{L a n g e v i n} d t} u_{q} = u_{q},

(58)

e^{i L^{L a n g e v i n} d t} p_{q} = p_{q} e^{- γ d t} + \sqrt{ν_{q}^{'} (u_{1}) k_{b} T (1 - e^{- 2 γ d t})} R (t),

(59)

where γ is the friction parameter, T is the temperature, k_B is Boltzmann’s constant, and R(t) are independent and identically distributed normal random numbers with mean 0, variance 1, and 〈R(t)R(t′)〉 = δ(t − t′). This thermostating method can also be adapted such that each mode is thermostated with a different optimized friction constant, as done in the path integral Langevin equation (PILE) scheme by Ceriotti et al.^51,63

Using the integration scheme of Eq. (57), we can perform efficient adaptive quantum-classical PIMD simulations with either the AKM or the CKM approach. It is derived in a rigorous fashion from a symplectic Hamiltonian and is also consistent with PI quantization, provided that the criteria in Eqs. (8) and (9) are satisfied. It is computationally advantageous over full-quantum simulations because in the CL region, all forces between interacting ring polymers can be approximated by a single calculation between the centroids. Furthermore, because the ring polymers are collapsed in the CL region and interact classically, the integration of the internal motion, i.e., of the higher modes, can be stopped and the ring polymers can be frozen in this part of the system. In the CKM case or for full-quantum systems, the algorithm reduces to a regular velocity Verlet scheme with multiple time-stepping and Langevin thermostating.

The derived quantum-classical multiresolution scheme can be combined with other optimization techniques for PI simulations. For example, the non-bonded forces in the QM region could also be evaluated using fewer than P beads, via the RPC scheme by Markland and Manolopoulos.^52,53 Alternatively, instead of white noise Langevin thermostating, we could make use of a colored noise thermostat. It was shown by Ceriotti et al. that a carefully parametrized PI GLE can lower the number of Trotter beads required for converged quantum behavior.^62–65 As our approach reduces the overall computational effort of a PI simulation by restricting the QM region of the system, it is complementary to these methods, which reduce the numerical complexity of the PI interactions themselves.

IV. APPROXIMATE QUANTUM DYNAMICS

In PIMD, we only measure quantum statistical properties. In the following, we will discuss how our integration scheme can be extended to allow for multiresolution quantum-classical CMD and RPMD with only minor changes.

A. Quantum-classical centroid molecular dynamics

Centroid molecular dynamics (CMD) is a method for the calculation of real-time quantum correlation functions in the short-time limit.^36–47 It is based on the notion that approximate quantum dynamical properties can be calculated from the time evolution of the centroid subject to the potential of mean force generated by the non-centroid modes of the ring polymer. Formally, this potential is obtained by integration over all possible ring configurations with constrained centroid position. This would not only be computationally expensive but practically intractable. The idea of CMD is to adiabatically decouple the internal fluctuations of the ring polymers from the centroid motion. By rescaling the higher mode kinetic masses (k > 1) with a sufficiently small adiabaticity parameter $0 < γ_{C M D}^{2} < 1$ such that $ν_{k}^{'} \to γ_{C M D}^{2} ν_{k}^{'}$ ⁠, the higher modes can be forced to evolve significantly faster than the centroid. Thereby, the centroid potential of mean force of the ring polymer is generated “on the fly” during the simulation.¹³² It has been shown, however, that in practice a partial adiabatic decoupling is sufficient for most applications.⁴⁵ In addition to the mass rescaling, the higher modes alone are coupled to thermostats such that the centroid dynamics remains Newtonian.

The previously described implementation of CMD, i.e., the removal of the thermostat from the centroid, and the kinetic mass rescaling, can also be done easily in our quantum-classical multiresolution scheme. In practice, one would typically be interested in the quantum dynamics only in the QM region. Hence, it suffices to only remove the centroid thermostat in this region. Then we could measure approximate quantum dynamical properties in the QM region while the classical part would still behave as in the canonical ensemble and could serve, for example, as a particle reservoir for the QM region. Another relevant scenario is the simulation of a complex biomolecule like a protein. In this case, an overall large simulation box would be required to preserve the structure and the solvating environment of the system, although we may want to probe the dynamics only in a specific subregion, such as near the protein’s active site.¹⁰¹

B. Quantum-classical ring polymer molecular dynamics

An alternative approach to the calculation of approximate quantum dynamics is provided by ring polymer molecular dynamics (RPMD).^{45,46,48–50} In the normal RPMD approach, the kinetic masses are chosen to be the real physical masses (as was done above for PIMD) and no thermostats are used so that the ring polymer evolution is completely Newtonian. In comparison to CMD, RPMD uses the entire chain to approximate quantum time correlation functions. However, the internal ring fluctuations can lead to artifacts when measuring, for example, vibrational spectra.¹⁰² To overcome this deficiency, Rossi et al. recently proposed a thermostated ring polymer molecular dynamics (TRPMD) approach⁵¹ that shares some characteristics of both normal RPMD and CMD. In TRPMD, the kinetic masses are also chosen to be the real physical masses and measurements are performed using all of the beads. However, as in CMD, Langevin thermostats are attached to all higher modes for which k > 1. Provided the thermostats are adjusted carefully, TRPMD avoids both the spurious resonances in the vibrational spectra and also the curvature problem of CMD¹⁰² while retaining the appealing properties of RPMD. An ideal choice for the Langevin friction parameters in TRPMD is given by the PILE scheme.^51,63 In the PILE method, each higher mode k > 1 is thermostated with a different optimized coupling constant γ_k based on the mode vibration frequency as $γ_{k} = ω_{P} \sqrt{ξ_{k} P}$ ⁠. Note that it was recently shown how both CMD and (T)RPMD can be derived as different approximations to the so-called Matsubara dynamics.¹⁰³

Just as with CMD, TRPMD simulations can be easily run with our quantum-classical PI scheme and the corresponding integrator, Eq. (57). We only need to adapt the thermostats on the different modes accordingly and remove the thermostat on the centroid.

Note that in the AKM approach, the kinetic masses will change in the CL region. However, in this part of the simulation box, the spring masses also have different values and we are typically not interested in the dynamics anyway. Thus, one may also reintroduce the centroid thermostat in this outer region.

V. SIMULATIONS

To validate the proposed adaptive resolution PIMD approach, we implemented it in the ESPResSo++ molecular simulation package¹⁰⁴ and performed adaptive resolution simulations of liquid water. Nuclear quantum effects in liquid water have been thoroughly investigated and shown to be important for an accurate description of its structure and dynamics.^22–25 Hence, water is an ideal test case for the method.

A. Water system

We consider a system of 918 water molecules in a slab-shaped box with dimensions L_x = 6.92 nm, L_y = L_z = 2.0 nm (33.168 molecules/nm³), and periodic boundary conditions in all directions. The resolution changes along the X-direction, the full width of the QM region is set to d_QM = 2.0 nm, and the width of the adjacent HY regions to d_HY = 1.5 nm. The resolution function is given by a squared cosine, commonly used in adaptive resolution simulations.^{72,79,90,92,105–107} We perform the simulations at a temperature of 300 K and use a Trotter number P = 32, as this has been shown to provide well-converged results for most dynamical and structural properties of water.^22,23,25 A simulation snapshot of the system is presented in Fig. 1(a).

We employed a water force field that was recently developed by Fritsch et al. specifically for PI simulations of bulk liquid water.²⁵ It is parametrized from ab initio density functional theory calculations using the force matching^108–110 and iterative Boltzmann inversion methods.¹¹¹ All interactions are mapped onto a set of short-ranged tabulated potentials and no explicit charges are present. Separate potentials are provided for the non-bonded O–O, O–H, H–H interactions, for the O–H bond and for the H–O–H angle. An additional bonded potential is applied between the two H-atoms of the same molecule. This force field describes the structural and dynamical properties of liquid water at 300 K and at a density of 1.1 g/cm² with reasonable accuracy. Furthermore, it is very efficient to simulate since it is purely short-ranged with an interaction cutoff of 7.8 Å. We have chosen this potential for its numerical efficiency and its suitability for PI simulations, and we note that the present adaptive resolution methodology can also be applied to analytical potentials that have short-range and/or long-range contributions.

In order to collapse the ring polymers in the CL region, we choose M_α = 100 m_α for all particles α. Because of their point-like structure, we only use the centroids to calculate non-bonded and bonded interactions between atoms in the CL region [see Fig. 1(b)]. Furthermore, we stop the integration of the higher modes in the CL region, i.e., we freeze the internal degrees of freedom of the ring polymers. Note that the setup satisfies the criteria in Eqs. (8) and (9) and can therefore be considered to be consistent with formal path integral quantization.

We run simulations using both the AKM and CKM approaches, although we focus on the AKM method, for which we have derived a non-standard integrator and which allows larger time steps. In all simulations, the kinetic mass of the centroids is given by Eq. (26), which corresponds to the real mass. For the AKM simulations, we choose the kinetic masses of the higher modes according to Eq. (27). As already argued, this corresponds to the real masses in the QM region, which facilitates a realistic ring polymer time evolution and therefore allows the calculation of approximate quantum dynamical properties from RPMD simulations. In the CL region, the higher mode masses are increased as explained previously. For the CKM simulations, we choose the masses in a similar way although they remain constant over all simulation domains. For CMD simulations, we introduce an additional rescaling of the higher modes’ kinetic masses with an adiabaticity parameter $γ_{C M D}^{2} = 0.05$ ⁠. Note that we do not rescale the kinetic masses of the higher modes with the eigenvalues of the normal mode transformation, as it is often done in CMD.⁴⁵ We also keep the additional 1/P factor in the kinetic masses, which we introduced earlier to ensure dynamics on the correct time scale. Therefore, the ring polymers’ higher modes vibrate with frequencies $ω_{P} \sqrt{ξ_{k} P / γ_{C M D}^{2}}$ ⁠.

To enforce the correct temperature, we couple all modes to white noise Langevin thermostats. The centroid mode is thermostated with a friction constant γ = 2.0 ps⁻¹, except in CMD and TRPMD simulations, where no thermostat is applied on it. For the higher modes k, we employed the path-integral Langevin equation scheme by Ceriotti et al.^51,63 and used friction coefficients $γ_{k} = ω_{P} \sqrt{ξ_{k} P}$ ⁠, which are proportional to the normal-mode frequencies (for CMD simulations, we used $γ_{k} = ω_{P} \sqrt{ξ_{k} P / γ_{C M D}^{2}}$ ⁠). The path-integral Langevin equation method leads to optimal sampling and can also be applied in the context of TRPMD simulations.

The present adaptive quantum-classical simulation method allows us to change not only the quantum delocalization of the particles but also their non-bonded interaction potentials. This was demonstrated in our previous paper in simulations of liquid parahydrogen.⁷⁹ Here, we perform the majority of the validation simulations using the same interaction potential in both the QM and the CL regions. We also test the scheme using a different potential in the CL region, specifically a purely repulsive Weeks-Chandler-Andersen (WCA) potential¹¹² of the form

V^{C L} (r) = \{\begin{matrix} 4 ϵ [{(\frac{σ}{r})}^{12} - {(\frac{σ}{r})}^{6} + \frac{1}{4}] & : r \leq R_{c} \\ 0 & : r > R_{c} \end{matrix})

(60)

with ϵ = k_BT, σ = 0.25 nm, and $R_{c} = 2^{\frac{1}{6}} σ = 0.28 nm$ ⁠. The potential acts only between the oxygen atoms. Note that the intramolecular bonded interactions are kept in the CL region to prevent the molecules from disintegrating.

B. Setups

We perform simulations employing the following setups:

Setup 1: We use adaptive kinetic masses and the same interaction potentials in both regions. Applying thermostats to all modes, we calculate various structural properties of the water in the QM region. Additionally, we remove the thermostat from the centroid and use TRPMD to calculate several dynamical quantities.
Setup 2: This setup is the same as setup 1 except the kinetic masses of the higher modes are rescaled with the adiabaticity parameter $γ_{C M D}^{2}$ ⁠. We then calculate the dynamical properties via CMD.
Setup 3: The AKM method is applied as in setup 1, but the WCA potential is employed to model intermolecular interactions in the CL region. In this scenario, we only validate the coupling by calculating density profiles as well as profiles of the atomistic radii of gyration.
Setup 4: We switch to constant kinetic masses and employ the same interaction potentials in both regions. As in setup 3, we only validate the coupling and calculate density profiles as well as profiles of the atomistic radii of gyration.

Additionally, we perform full-quantum and full-classical (P = 1) reference simulations without any interpolation between different particle masses or interaction potentials. All simulation parameters, including the box dimensions, are the same as for the adaptive simulations. The only exception is the friction constant of the Langevin thermostat, which is set to 10.0 ps⁻¹ in the full-classical simulations.

The time steps used in the simulation setups are presented in Table I. For the full-classical and the full-quantum reference simulations, we use the same time steps as in the corresponding adaptive resolution setups. The time steps in the table refer to those used in equilibration simulations, during the derivation of the free energy correction and the thermodynamic force (see Sec. V C), as well as during all other simulations sampling statistical averages. For the calculation of dynamic quantities in setup 1, we reduce all time steps to the same ones as used in the CMD simulations in setup 2. We do this for two reasons: On the one hand, our implementation of the integration scheme allows one to print out positions or velocities only after a full step Δt. Hence, this large time step needs to be short enough to allow a fine sampling when calculating, for example, velocity autocorrelation functions. On the other hand, we wish to avoid artifacts resulting from the use of different time steps when comparing CMD to TRPMD. Note, however, that only few and very short simulations need to be run with this modification. The majority of simulations use the time steps in Table I.

TABLE I.

Time steps for the quantum-classical adaptive resolution PIMD, RPMD, and CMD simulations and for the reference calculations.

Setup	Δt (fs)	δt (fs)	dt (fs)
#1: AKM method, same potentials, TRPMD	2.0	0.5	0.05
#2: AKM method, same potentials, CMD	0.4	0.1	0.01
#3: AKM method, WCA potential	2.0	0.5	0.05
in CL region, TRPMD
#4: CKM method, same potentials, TRPMD	1.0	0.1	0.006 25

Setup	Δt (fs)	δt (fs)	dt (fs)
#1: AKM method, same potentials, TRPMD	2.0	0.5	0.05
#2: AKM method, same potentials, CMD	0.4	0.1	0.01
#3: AKM method, WCA potential	2.0	0.5	0.05
in CL region, TRPMD
#4: CKM method, same potentials, TRPMD	1.0	0.1	0.006 25

In general, all time steps are chosen to be as large as possible but still sufficiently small to accurately sample phase space, retain an acceptable level of energy-conservation in microcanonical test simulations, and generate the correct temperature in simulations in the canonical ensemble. The time steps we find to work well seem reasonable: In classical simulations, updating the regular non-bonded forces every 1-2 fs is a frequent choice,^{92,106,113–116} while the vibration frequency for the bonds and angles in water demands a time step of around 0.5 fs.^25,117 The vibrational frequency of the springs between the PI beads is yet higher, requiring an even smaller time step. Furthermore, CMD simulations are known to require particularly small time steps, as the internal ring polymer motion is strongly accelerated. A similar effect is observed in simulations with the CKM approach. In this case, the internal motion of collapsed ring polymers is also significantly enhanced (we mentioned that in the CL region, the ring polymers are frozen. However, the ring polymers are already strongly collapsed at the outer parts of the HY region, where a full integration of the internal motion is still necessary to accommodate the gradual collapse and extension of the ring polymers). Therefore, it becomes clear that the AKM scheme is better suited for the proposed adaptive quantum-classical simulation protocol than the naive CKM method. We stress, however, that finding optimal time steps is not the primary goal of this work; clearly, there is room for further fine-tuning.

For the calculation of all structural quantities and statistical averages, we run simulations of duration 200 ps, if not otherwise indicated. Additionally, we perform short 2 ps runs during which we calculate velocity autocorrelation functions and vibrational spectra. We also measure hydrogen bond population fluctuations, which is done in simulations of duration 32 ps. In all cases, we start from equilibrated configurations, run 10 independent simulations, and average the results.

Because averaging is performed over independent runs as well as over the different imaginary time slices and because roughly 200 molecules [see Fig. 7(b)] are used for the calculation of the quantities of interest in the QM region of adaptive simulations or in full-quantum reference simulations, the chosen lengths of the simulations are sufficient to obtain well-converged results. Furthermore, the durations of the simulations for the calculation of dynamical quantities were chosen to cover the relevant time scales of the studied processes. If not otherwise indicated, the statistical error is within line thickness.

C. Free energy corrections

To correct for the thermodynamic imbalance between the low-mass QM and the high-mass CL region, we apply a free energy correction (FEC) ΔH.^{79,87–90,92} We derive the FEC via Kirkwood thermodynamic integration (KTI)¹¹⁸ between the fully CL (λ = 0) and the fully QM (λ = 1) systems, and we calculate the averages

1 . {⟨\frac{1}{N P} \sum_{α = 1}^{N} \sum_{i = 1}^{P} \{V_{α, i}^{Q M} (u) - V_{α, i}^{C L} (u)\}⟩}_{λ},

(61)

2 . {⟨\frac{1}{N} \sum_{α = 1}^{N} \sum_{i = 2}^{P} \frac{1}{2} \frac{(M_{α} - m_{α})}{P {\tilde{ν}}_{α} {(u_{α, 1})}^{2}} p_{α, i}^{2}⟩}_{λ},

(62)

3 . {⟨\frac{1}{N} \sum_{α = 1}^{N} \sum_{i = 2}^{P} \frac{1}{2} ξ_{i} (m_{α} - M_{α}) ω_{P}^{2} u_{α, i}^{2}⟩}_{λ},

(63)

as well as the pressure p(λ) for a set of 101 λ’s along the integration path from λ = 0 to λ = 1. The KTI is run in a smaller box of dimensions L_x = L_y = L_z = 3.0 nm. All other simulation parameters are as explained above except for the thermostat friction of the centroid mode, which was set to 10 ps⁻¹ to achieve rapid equilibration after changing λ. We start the KTI from an equilibrated system at λ = 0 and we perform for each λ a short 0.3 ps equilibration run (for setup 2, only 0.12 ps due to the short time step) and another 1.5 ps run (for setup 2, only 0.6 ps) during which we take measurements. From these results, we construct the FEC ΔH to cancel the averages of the drift forces, Eq. (30), and the pressure difference between the subsystems. Since the calculation of the FEC via KTI is an approximate method to correct for the thermodynamic imbalance, we refine the FEC using the thermodynamic force (TF) scheme.^69,79 The TF is an iterative approach that directly constructs a correction force in the HY region from the distorted density profile along the direction of resolution change in order to flatten the density throughout the system. Each TF iteration consists of a 50 ps equilibration run (10 ps for setup 2 and 15 ps for setup 4) and a 150 ps production run (30 ps for setup 2 and 20 ps for setup 4) during which we sample the density. We perform 20 iterations for each setup.

Although we calculate the quantities in Eqs. (61)–(63) separately for the oxygen and hydrogen atoms, we determine global molecular pressures instead of species-wise partial pressures. After constructing the correction force from the pressure p(λ), we distribute it between oxygen and hydrogens proportionally to their masses. Finally, the FEC is applied on the atomistic level based on the atom’s centroid positions.

The calculation of the FEC via KTI and the iterative correction are both well-established methods for achieving a smooth coupling in adaptive resolution simulations. See, for example, Refs. 69, 79, 87–90, and 92 for further technical details.

VI. RESULTS

A. Structure

We first investigate the structural properties of the adaptive quantum-classical water systems. Figure 2 shows the density profiles along the x-direction of the four setups without correcting for the thermodynamic imbalance (green curves), with FEC but without iterative refinement (blue curves), and with FEC including the iterative refinement via TF (red curves). Without any corrections, the density is strongly distorted. Applying the non-iterative FEC significantly improves the coupling between the regions, although the density in the QM region is still slightly too low for setups 1-3 and much too low for setup 4. This can be expected, as the non-iterative FEC is approximate, and statistical inaccuracies can occur when it is generated via KTI. Refining the FEC with the iterative TF technique, we obtain flat density profiles for all setups, except setup 4, for which significant deviations in the HY region remain. Note that for setup 4, which uses the CKM scheme, we were not able to run stable simulations without any compensation. In this case, the drift forces are so strong that all molecules are immediately pushed to one subregion. In comparison, the AKM approach works better and requires a smaller FEC. Furthermore, it is interesting to note that for setup 3, which features a WCA potential in the CL region, the density profile without any corrections is significantly flatter than the uncorrected profiles for setups 1 and 2. This suggests that the different pressure of the WCA potential together with the additional drift term leads to an overall more balanced equilibrium state between the QM and CL regions. Since these effects are concurrent, it is difficult to disentangle one from the other.

FIG. 2.

View large Download slide

Normalized density profiles ρ in the quantum-classical adaptive resolution simulations of liquid water for the four different setups without FEC (green), with KTI-based FEC but without iterative refinement (blue), and with FEC including iterative refinement via TF (red). For setup 4, we are not able to run stable simulations without any FEC (setup 1: AKM, same potentials, TRPMD. Setup 2: AKM, same potentials, CMD. Setup 3: AKM, WCA potential in CL region, TRPMD. Setup 4: CKM, same potentials, TRPMD).

The derivation of the FEC via KTI and several iterations of TF may seem cumbersome. However, both the KTI and the TF iterations can be run using simulation setups that are much smaller than the actual system. For large applications, this step will likely take significantly less time than the simulation of the complete system. Additionally, more advanced approaches have recently been developed that efficiently calculate the FEC on the fly during the simulation of the full system or a representative smaller one.¹¹⁹ Note that a FEC or a similar compensation force is required in all adaptive resolution methods that allow a free exchange of particles between subregions that feature different thermodynamics.^{69,73,74,87,88,119} All results reported below are calculated in setups in which the refined FEC is applied.

Figure 3 presents the radii of gyration of the ring polymers corresponding to the water’s oxygen and hydrogen atoms as a function of their position along the x-direction. In the QM region, the radii of gyration perfectly match with those from full-quantum reference simulations, while in the CL region, they drop by ≈90% (also see Fig. 1). Therefore, the molecules exhibit their full quantum character in the QM region, while in the CL region, the ring polymers shrink to nearly point-like particles and behave classically. To collapse the ring polymers even further, one would simply need to choose a heavier particle mass M_α in the CL region. The radius of gyration in the CL region is approximately proportional to $1 / \sqrt{M_{α}}$ ⁠. Note that the data in Fig. 3 correspond to system 1 and that the other setups show the exact same behavior.

FIG. 3.

View large Download slide

Radii of gyration of the ring polymers corresponding to oxygen and hydrogen atoms as a function of the ring polymers’ position along the x-direction. The magenta and cyan lines correspond to the radii of gyration of atoms in full-quantum reference simulations.

Using setup 1, we also calculate the radial distribution functions (RDFs) and the tetrahedral order parameter q_tet within the QM region (Fig. 4). We left a small buffer of 0.25 nm at the interface to the HY region and considered the inner 1.5 nm of the QM region in order to avoid artifacts by molecules at the outer edges of the QM region that interact strongly with molecules in the HY region. For a molecule i, q_tet is given by

q_{t e t} = 1 - \frac{3}{8} \sum_{j = 1}^{3} \sum_{k = j + 1}^{4} {(c o s (θ_{j, k}) + \frac{1}{3})}^{2} .

(64)

The indices j and k run over i’s four nearest neighbor molecules and the angle θ_j,k is formed by the oxygen atoms of molecules i, j, and k with i in the center. The order parameter q_tet is defined such that it is 1 when the molecule forms a perfect tetrahedron with its four nearest neighbors and on average 0 for an ideal gas. The RDFs and q_tet in the QM region of the adaptive quantum-classical water systems match nearly perfectly the results from full-quantum reference simulations. While we do not find any significant impact of quantum effects on the tetrahedral order parameter q_tet, the H–H and O–H RDFs show a significant dependence on the PI description. g_OO(r), on the other hand, remains largely unaffected. These observations are consistent with previous work.²⁵ We conclude that the PI-based water structure in the QM region is undisturbed by the coupling to the CL particle reservoir.

FIG. 4.

View large Download slide

(a) (Bead-bead) RDFs of the quantum-classical adaptive resolution simulations calculated in the QM region and of full-quantum and full-classical (P = 1) reference simulations. (b) Same for the tetrahedral order parameter q_tet.

B. Dynamics

We also probe the dynamics in the inner QM region of the adaptive quantum-classical water systems. First, we calculate the vibrational spectrum from the water molecules’ velocity autocorrelation function. We do this both via TRPMD (setup 1) and CMD (setup 2) and compare the results to full-quantum and full-classical (P = 1) reference simulations (Fig. 5). The vibrational dynamics in the QM region perfectly reproduces the full-quantum reference data, both for CMD and TRPMD. While CMD and TRPMD give similar results, the classical system shows blue shifts in the H–O–H bending and O–H stretching modes. The spectra also agree with the results from Fritsch et al.²⁵

FIG. 5.

View large Download slide

Vibrational density of states in the QM region of the quantum-classical adaptive resolution simulations of liquid water and in full-quantum and full-classical reference simulations, calculated using TRPMD and CMD. The dashed vertical lines indicate the diffusion mode, the H–O–H bending mode, and the O–H stretching mode.

As hydrogen bonds play a critical role in the behavior of water,^120–122 we additionally assess the hydrogen bonding kinetics of water in the QM region. The breaking and forming of hydrogen bonds can be characterized by the correlation function

C (t) = \frac{⟨ h (0) h (t) ⟩}{⟨ h ⟩},

(65)

which measures fluctuations in the hydrogen bond populations throughout the system.^123,124 The hydrogen bond population operator h(t) is 1, if a particular pair of molecules is hydrogen bonded and 0 otherwise [〈h〉 denotes the average of h(t)]. We consider two molecules to be hydrogen bonded if the distance between their oxygen atoms is <3.5 Å and the angle between the O–O axis and one of the O–H bonds is <30°. Based on Eq. (65), we can determine the hydrogen bond relaxation rate k(t) as

k (t) = - \frac{d C (t)}{d t} .

(66)

The quantity −k(t) can be interpreted as the average rate of change of hydrogen bonds that are broken at time t later. It has been widely used in studies of the hydrogen bond kinetics in liquid water.^123–128

We employ the centroids for measuring the hydrogen bonds and calculate C(t) via CMD using setup 2. We also perform full-quantum and full-classical reference simulations. The results are shown in Fig. 6. The hydrogen bonding kinetics in the QM region of the adaptive system reproduces the full-quantum reference within the statistical error (see shaded regions in Fig. 6). We conclude that the hydrogen bond kinetics in the QM region of the adaptive simulations is well preserved. Furthermore, we observe no quantum effects. The classical and the quantum systems behave the same in their hydrogen bonding dynamics.

FIG. 6.

View large Download slide

Hydrogen bond population fluctuations characterized via the correlation function C(t) in the QM region of the quantum-classical adaptive resolution simulations and in full-quantum and full-classical reference simulations. The shaded regions indicate the standard deviations of the data corresponding to the full-quantum and adaptive simulations. Inset: average hydrogen bond relaxation rates k(t) in a semi-log plot.

We conclude that both the water structure and the PI-based dynamics in the QM region are unaffected by the coupling to the CL domain. Importantly, we have shown that one can apply both CMD and TRPMD in the proposed quantum-classical adaptive resolution simulation scheme.

C. Particle fluctuations

It is important that the proposed method allows for a free flow of particles through the HY region without any barriers. The QM region must behave as if embedded in an overall QM environment. To test this, we label all molecules that reside at the beginning of a simulation in the inner QM region, leaving a buffer of 0.25 nm. We then track how many of the labeled particles remain in the inner QM region after time t and compare this to a full-quantum reference simulation in which we label and track all molecules in a similar subregion of the system. Note that we keep the thermostat on the ring polymers’ centroid modes for these simulations. The results are presented in Fig. 7(a) and show that the particles diffuse out of the QM region in the adaptive setup in a similar fashion as in the full-quantum system (for a more detailed study of the diffusion properties of the water model employed in this work, see Fritsch et al.²⁵).

FIG. 7.

View large Download slide

(a) Fraction of molecules that resided at the beginning of the simulation in the inner quantum region and that still stayed there after time t for an adaptive system (setup 1) and for a full-quantum reference system, considering the same subregion. The thick lines denote the average from 10 simulations (thin lines). (b) Particle number probability distribution of the inner part of the quantum region in the adaptive simulations and of the same volume in a reference full-quantum simulation. The black curve is a Gaussian fit to the latter (μ = 199.2, σ = 3.4). The error bars denote the standard deviation of the adaptive simulation data. The statistical error of the full-quantum data is similar.

Additionally, we measure the particle number fluctuations in the inner QM region [Fig. 7(b)]. The fluctuations match the full-quantum reference nearly perfectly. Note that the data for the adaptive system in Fig. 7 correspond to setup 1. All other setups show similar behavior.

The results indicate that the HY region allows a free exchange of molecules between the CL and QM regions and that the QM region exchanges particles with its environment as if embedded in a full-quantum environment. Considering the complexity of the setup, the application of a correction force in the HY region and the different structures and thermodynamics in the CL and QM subsystems, this is non-trivial. Because of the free flow of particles and the correct particle number fluctuations in the QM region, the scheme can, for example, be used for efficient simulations of open quantum systems.

D. Computational speedup

In practice, the proposed scheme will be useful in applications where the QM region only occupies a small area in the overall system compared to the CL region. To estimate its speedup compared to full-quantum simulations, we therefore set up a box with a CL domain 10 times as large as the QM and HY regions together. The system is the same as setup 1, but the simulation box was extended in the x-direction to L_x = 55.0 nm. The full width of the QM region is again set to d_QM = 2.0 nm and the width of the adjacent HY regions to d_HY = 1.5 nm. Hence, the central 5 nm of the box is occupied by the QM and HY regions and they are surrounded by 50 nm of the CL region. Starting from equilibrated configurations, we run short simulations of 1 ps and measure the average time required for one full step of the integrator [cf. Eq. (57)] as well as the time spent by our code only for interatomic force computations. Note that because of the MTS scheme, one full step constitutes 40 iterations of the integrator’s innermost loop and 4 iterations of the central loop (N = 4, n = 40). Additionally, we perform full-quantum reference simulations and estimate the speedup as the ratio t_QM/t_adaptive, where t_QM and t_adaptive denote the time spent by the full-quantum and the adaptive simulations. All simulations are run using the same implementation of the method and we use a single CPU to provide parallelization independent results. The data are shown in Table II.

TABLE II.

Time spent for a full integration step and for the interatomic force computation only in adaptive and full-quantum reference simulations.

	Adaptive	Full-quantum	Speedup
Total time per full step	10.4 s	26.5 s	2.5
Interatomic force computation time per full step	1.1 s	9.2 s	8.4

	Adaptive	Full-quantum	Speedup
Total time per full step	10.4 s	26.5 s	2.5
Interatomic force computation time per full step	1.1 s	9.2 s	8.4

Since we can evaluate the interatomic forces in the CL region in the adaptive setup using only the centroids due to the point-like structure of the ring polymers, the cost of computing the interatomic forces is reduced by a factor of about 8.4. This is close to the theoretical upper bound of 11, which one can obtain by assuming that the force computation in the CL region becomes completely negligible and comparing the volume of the combined QM/HY region with the total simulation domain in the full-quantum simulation. However, in our MTS integration scheme, the expensive interatomic forces are evaluated infrequently compared to the internal ring polymer vibrations. Furthermore, the employed interatomic force field has a very short cut off and is highly efficient. Therefore, the overall speedup is significantly lower and the integration scheme’s overhead plays a significant role in the overall time spent by the code on the calculations. Nevertheless, we find that the overall simulation speedup is still 2.5. Much larger speedup factors would be expected if the recently introduced stochastic isokinetic approach of Ref. 133 had been employed.

These results suggest that the proposed method is indeed practically useful for large-scale systems. In realistic applications, the QM region will likely constitute only a tiny fraction of the simulation box. Considering the example of a protein in which only the active site is modeled using nuclear quantum effects via a small spherical QM region, the CL domain will be larger than the QM domain by much more than 10, the factor used here. Furthermore, many practical force fields are longer ranged and computationally more demanding to evaluate than the one demonstrated here. In such a situation, the evaluation of the interatomic potentials will likely dominate the overall computation and the integrator’s overhead will become negligible. Therefore, we expect the speedup in large-scale applications to be even greater than that demonstrated here. For this work, a significant amount of time was spent in the generation of the FEC from scratch. In general, however, the FEC can be derived in a much smaller box than the one actually used for production simulations later. Additionally, there exist recently developed and more efficient techniques that enable a quick calculation of the FEC on-the-fly without any lengthy iterative procedure, nearly completely eliminating the computational cost of separately deriving it.¹¹⁹ Lastly, a standard library of FEC’s can be compiled to be employed as a starting point for iterative refinement to be then employed in specific cases differing, e.g., by temperature, density etc.

Finally, it needs to be mentioned that, in practice, the method’s computational efficiency strongly depends on its implementation and, in particular, on the parallelization scheme used when making use of many CPUs. For real-life applications, it would thus make sense to employ an appropriate load-balancing scheme that can concentrate the computational resources in the QM region.

VII. DISCUSSION AND CONCLUSIONS

We have proposed and validated a concurrent multiscale method for Hamiltonian adaptive resolution molecular dynamics simulations using the PI formalism. The scheme is based on a position-dependent particle mass, which controls the extension and collapse of the ring polymers. In the QM region, where the particles have their real masses, the ring polymers are extended, while in the CL region, where the mass is increased, the ring polymers collapse to point-like particles. Therefore, the interaction becomes classical and the dynamics obeys classical Newtonian mechanics in the CL region. The particles freely diffuse between the two regions and change their description on the fly. The method allows a more efficient evaluation of forces and energies in the CL domain, which leads to a speedup compared to full PI simulations. Importantly, we provide criteria that quantify to what extent such an adaptive PI setup is consistent with a bottom-up PI quantization. We want to point out that this differentiates our approach from related methodologies which are based on a direct interpolation of the forces corresponding to a classical and a PI system.^70–74 These techniques do not allow a Hamiltonian description of the system,⁷⁵ which, however, is the basis for a bottom-up PI treatment in the first place. Our scheme aims at overcoming this limitation. It allows both adaptive PIMD simulations sampling quantum statistical averages and quantum-classical RPMD and CMD, which enable us to calculate approximate quantum dynamical quantities and time correlation functions. Finally, the method allows one not only to selectively turn on and off nuclear quantum effects in different regions but also to change the intermolecular interaction potential. In this way, one can use a more efficient, possibly coarse-grained model in the CL region. This would be useful, for example, when the CL domain only serves as a particle reservoir.

To implement our methodology in a molecular dynamics framework, a kinetic energy term needs to be introduced into the configurational energy obtained from PI quantization. At this point, our approach can be implemented in two different ways: The kinetic masses in this kinetic energy term can be chosen either to be constant throughout the whole system (CKM) or to vary in a way similar to the particle masses that control the springs between the PI beads (AKM). The CKM approach results in a simpler scheme that can be integrated with a standard velocity Verlet integrator. However, it leads to strong thermodynamic imbalances between the CL and the QM regions and requires very small time steps due to the accelerated vibrations of the ring polymers in the CL region. On the other hand, the AKM method requires a more sophisticated integration scheme due to the position-dependent kinetic masses. We derived an integrator that is tailored to the problem, employs multiple time-stepping, and allows a symplectic integration of the equations of motion. This scheme also facilitates time steps that are much larger than those in the CKM protocol and that are similar to those used in normal PI simulations. The AKM method also enables a smoother connection of the CL and QM systems, requiring a milder correction force in the HY coupling region.

The new integrator may appear complicated, but it requires little additional overhead in practice. In molecular dynamics simulations, most time is typically spent for non-bonded force calculations and for inter-processor communication. However, these two tasks do not need to be performed more often than in a standard velocity Verlet integrator. In essence, our methodology elegantly decouples the change of the particles’ quantum character, which is connected only to the higher modes and requires an additional decomposition step in the inner loop of the integrator, from the interatomic and intermolecular interactions, which are related to the more expensive bonded and non-bonded force calculations. In fact, it is only the masses of the higher modes that are position-dependent, while the masses of the centroid modes are constant. In our implementation, the integration is performed in normal mode space. Only before the calculation of the bonded and non-bonded forces, the particles’ real positions are updated and the force calculation is performed in real space. Afterwards, the forces are transformed back to normal modes. Therefore, the innermost loop of our integrator, in which the additional decomposition step occurs, does not require any inter-processor communication.

As we demonstrated, the new integration method, including the efficient evaluation of interatomic forces in the CL region due to the point-like ring polymer structure in this domain, indeed leads to a large computational speedup over corresponding full-quantum simulations. While in our benchmark simulations, the overhead of the integrator was still significant, in systems in which only a tiny part of the simulation domain is modeled quantum mechanically, and where the interatomic potentials are computationally more demanding than the highly efficient force field employed in this work, the additional overhead of the method will become negligible and the gain in computational efficiency over a similar full QM system will be even higher. In general, the speedup depends on various factors such as the system at hand, the scheme’s implementation, the parallelization methodology, and the load balancing protocol (in highly parallelized simulations that employ many CPUs, a suitable load balancing method that allows to concentrate computational resources in the QM region is crucial). Note that the present adaptive resolution method can, of course, also be used in setups with different geometrical arrangements of the QM and CL regions compared to that used in this article. A typical example would be a small spherical QM domain positioned at an area of particular interest within a large CL system.

The proposed adaptive PI simulation scheme gains its efficiency by restricting the QM region to a small but relevant region in space and treating the rest of the system with a more efficient classical model. This is in contrast to other approaches that aim to alleviate the computational cost of PI simulations by modifying the PI calculations themselves, such as RPC,^52,53 higher order Trotter factorizations,⁵⁷ or advanced thermostating techniques.^62–65 Our method is complementary to these approaches and could be combined with them. For example, one could apply RPC to further reduce the numerical effort in the QM region or use a colored noise instead of a white-noise Langevin thermostat, which would allow one to employ less Trotter beads. One could also further improve the multiple time-stepping and tailor it to the investigated systems.

The applications of the proposed methodology are diverse. The scheme is useful whenever only a small subdomain of an overall large system needs to be treated using path integrals. This could be the case, for example, in biomolecular systems, in which the study of nuclear quantum effects has gained significant interest.^{6–13,16–21,101} Biological systems are often complicated and quantum delocalization plays an important role usually only in a small part of the system, such as the active site of proteins.¹⁰¹ Our multiscale method could be used to describe the active site quantum mechanically and an efficient classical model could be employed for the rest of the system. This would allow for an extension of the accessible length and time scales compared to full path integral simulations. Similar applications of the scheme are simulations of interfaces or membranes. The possibility to selectively switch on and off the nuclear quantum effects in different regions also allows one to investigate the locality of quantum properties. One can ask, for example, how much quantum mechanically modeled environment is required to support the quantum mechanical features in a certain subregion.^{115,129–131} This would not be possible in bulk path integral simulations in a straightforward way. Furthermore, the method enables an efficient simulation of a quantum grand canonical ensemble: a QM region can be coupled to a large particle reservoir, which itself is described classically. Yet another interesting possible application of our methodology is its combination with QM/MM techniques, which concurrently couple ab initio and classical empirical force fields. Recently, Boereboom et al.⁸¹ proposed an adaptive QM/MM method based on the Hamiltonian adaptive resolution scheme. The latter is also used in the PI-based adaptive resolution scheme presented in this article. In fact, although the interatomic potentials employed in this work are empirical force fields, the forces and energies could also come from ab initio calculations. Therefore, one could combine our approach with the one from Boereboom et al. and construct a Hamiltonian adaptive QM/MM scheme that also incorporates a multiscale treatment of PIs.

Finally, we would like to point out that the present concurrent multiscale PI simulation methodology has been implemented in the ESPResSo++ package¹⁰⁴ and is publicly available.

ACKNOWLEDGMENTS

K. Kreis is the recipient of a fellowship funded through the Excellence Initiative (DFG/GSC 266). M.E.T. acknowledges support from the National Science Foundation, Grant No. CHE-1566085. The authors thank Joe Rudzinski for careful reading of the manuscript and his helpful suggestions.

REFERENCES

1.

D.

Scharf

,

G. J.

Martyna

, and

M. L.

Klein

,

Fiz. Nizk. Temp.

19

,

516

(

1993

).

Google Scholar

2.

D. M.

Ceperley

,

Rev. Mod. Phys.

67

,

279

(

1995

).

https://doi.org/10.1103/revmodphys.67.279

Google Scholar

Crossref

3.

R. T.

Azuah

,

W. G.

Stirling

,

H. R.

Glyde

,

M.

Boninsegni

,

P. E.

Sokol

, and

S. M.

Bennington

,

Phys. Rev. B

56

,

14620

(

1997

).

https://doi.org/10.1103/physrevb.56.14620

Google Scholar

Crossref

4.

P.

Nozières

and

D.

Pines

,

The Theory of Quantum Liquids

, Advanced Book Classics (

Perseus

,

Cambridge, MA

,

1999

).

Google Scholar

5.

T.

Lindenau

,

M. L.

Ristig

,

K. A.

Gernoth

,

J.

Dawidowski

, and

F. J.

Bermejo

,

Int. J. Mod. Phys. B

20

,

5035

(

2006

).

https://doi.org/10.1142/s0217979206036077

Google Scholar

Crossref

6.

P.-O.

Löwdin

,

Rev. Mod. Phys.

35

,

724

(

1963

).

https://doi.org/10.1103/revmodphys.35.724

Google Scholar

Crossref

7.

R.

Rein

and

F. E.

Harris

,

Science

146

,

649

(

1964

).

https://doi.org/10.1126/science.146.3644.649

Google Scholar

Crossref

PubMed

8.

T. H.

Haines

,

Prog. Lipid Res.

40

,

299

(

2001

).

https://doi.org/10.1016/s0163-7827(01)00009-1

Google Scholar

Crossref

PubMed

9.

L. J. C.

Jeuken

,

R. J.

Bushby

, and

S. D.

Evans

,

Electrochem. Commun.

9

,

610

(

2007

).

https://doi.org/10.1016/j.elecom.2006.10.045

Google Scholar

Crossref

10.

A.

Pérez

,

M. E.

Tuckerman

,

H. P.

Hjalmarson

, and

O. A.

Von Lilienfeld

,

J. Am. Chem. Soc.

132

,

11510

(

2010

).

https://doi.org/10.1021/ja102004b

Google Scholar

Crossref

PubMed

11.

A. Y.

Smirnov

,

L. G.

Mourokh

, and

F.

Nori

,

J. Phys.: Condens. Matter

23

,

234101

(

2011

).

https://doi.org/10.1088/0953-8984/23/23/234101

Google Scholar

Crossref

PubMed

12.

D.

Jacquemin

,

J.

Zúñiga

,

A.

Requena

, and

J. P.

Céron-Carrasco

,

Acc. Chem. Res.

47

,

2467

(

2014

).

https://doi.org/10.1021/ar500148c

Google Scholar

Crossref

PubMed

13.

B.

Pamuk

,

J. M.

Soler

,

R.

Ramírez

,

C. P.

Herrero

,

P. W.

Stephens

,

P. B.

Allen

, and

M. V.

Fernández-Serra

,

Phys. Rev. Lett.

108

,

193003

(

2012

).

https://doi.org/10.1103/physrevlett.108.193003

Google Scholar

Crossref

PubMed

14.

W.

Li

,

H.

Yoshii

,

N.

Hori

,

T.

Kameda

, and

S.

Takada

,

Methods

52

,

106

(

2010

).

https://doi.org/10.1016/j.ymeth.2010.04.014

Google Scholar

Crossref

PubMed

15.

Y.

Nagata

,

R. E.

Pool

,

E. H. G.

Backus

, and

M.

Bonn

,

Phys. Rev. Lett.

109

,

226101

(

2012

).

https://doi.org/10.1103/physrevlett.109.226101

Google Scholar

Crossref

PubMed

16.

D.

Marx

,

M. E.

Tuckerman

,

J.

Hutter

, and

M.

Parrinello

,

Nature

397

,

601

(

1999

).

https://doi.org/10.1038/17579

Google Scholar

Crossref

17.

U. W.

Schmitt

and

G. A.

Voth

,

J. Chem. Phys.

111

,

9361

(

1999

).

https://doi.org/10.1063/1.480032

Google Scholar

Crossref

18.

M. E.

Tuckerman

,

D.

Marx

, and

M.

Parrinello

,

Nature

417

,

925

(

2002

).

https://doi.org/10.1038/nature00797

Google Scholar

Crossref

PubMed

19.

Y.

Wu

,

H.

Chen

,

F.

Wang

,

F.

Paesani

, and

G. A.

Voth

,

J. Phys. Chem. B

112

,

467

(

2008

).

https://doi.org/10.1021/jp076658h

Google Scholar

Crossref

PubMed

20.

T. C.

Berkelbach

,

H. S.

Lee

, and

M. E.

Tuckerman

,

Phys. Rev. Lett.

103

,

238302

(

2009

).

https://doi.org/10.1103/physrevlett.103.238302

Google Scholar

Crossref

PubMed

21.

D.

Marx

,

A.

Chandra

, and

M. E.

Tuckerman

,

Chem. Rev.

110

,

2174

(

2010

).

https://doi.org/10.1021/cr900233f

Google Scholar

Crossref

PubMed

22.

J.

Morrone

and

R.

Car

,

Phys. Rev. Lett.

101

,

017801

(

2008

).

https://doi.org/10.1103/physrevlett.101.017801

Google Scholar

Crossref

PubMed

23.

F.

Paesani

,

S.

Yoo

,

H. J.

Bakker

, and

S. S.

Xantheas

,

J. Phys. Chem. Lett.

1

,

2316

(

2010

).

https://doi.org/10.1021/jz100734w

Google Scholar

Crossref

24.

M.

Ceriotti

,

J.

Cuny

,

M.

Parrinello

, and

D. E.

Manolopoulos

,

Proc. Natl. Acad. Sci. U. S. A.

110

,

15591

(

2013

).

https://doi.org/10.1073/pnas.1308560110

Google Scholar

Crossref

PubMed

25.

S.

Fritsch

,

R.

Potestio

,

D.

Donadio

, and

K.

Kremer

,

J. Chem. Theory Comput.

10

,

816

(

2014

).

https://doi.org/10.1021/ct4010504

Google Scholar

Crossref

PubMed

26.

M.

Ceriotti

,

W.

Fang

,

P. G.

Kusalik

,

R. H.

McKenzie

,

A.

Michaelides

,

M. A.

Morales

, and

T. E.

Markland

,

Chem. Rev.

116

,

7529

(

2016

).

https://doi.org/10.1021/acs.chemrev.5b00674

Google Scholar

Crossref

PubMed

27.

R.

Feynman

and

A. R.

Hibbs

,

Quantum Mechanics and Path Integrals

(

McGraw-Hill

,

New York

,

1965

).

Google Scholar

28.

M. E.

Tuckerman

,

Statistical Mechanics: Theory and Molecular Simulation

(

Oxford University Press

,

New York

,

2010

).

Google Scholar

29.

D.

Chandler

and

P. G.

Wolynes

,

J. Chem. Phys.

74

,

4078

(

1981

).

https://doi.org/10.1063/1.441588

Google Scholar

Crossref

30.

M.

Parrinello

and

A.

Rahman

,

J. Chem. Phys.

80

,

860

(

1984

).

https://doi.org/10.1063/1.446740

Google Scholar

Crossref

31.

R. W.

Hall

and

B. J.

Berne

,

J. Chem. Phys.

81

,

3641

(

1984

).

https://doi.org/10.1063/1.448112

Google Scholar

Crossref

32.

B. J.

Berne

and

D.

Thirumalai

,

Annu. Rev. Phys. Chem.

37

,

401

(

1986

).

https://doi.org/10.1146/annurev.pc.37.100186.002153

Google Scholar

Crossref

33.

M. E.

Tuckerman

,

B. J.

Berne

,

G. J.

Martyna

, and

M. L.

Klein

,

J. Chem. Phys.

99

,

2796

(

1993

).

https://doi.org/10.1063/1.465188

Google Scholar

Crossref

34.

G. J.

Martyna

,

A.

Hughes

, and

M. E.

Tuckerman

,

J. Chem. Phys.

110

,

3275

(

1999

).

https://doi.org/10.1063/1.478193

Google Scholar

Crossref

35.

E. L.

Pollock

and

D. M.

Ceperley

,

Phys. Rev. B

30

,

2555

(

1984

).

https://doi.org/10.1103/physrevb.30.2555

Google Scholar

Crossref

36.

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

100

,

5093

(

1994

).

https://doi.org/10.1063/1.467175

Google Scholar

Crossref

37.

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

100

,

5106

(

1994

).

https://doi.org/10.1063/1.467176

Google Scholar

Crossref

38.

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

101

,

6157

(

1994

).

https://doi.org/10.1063/1.468503

Google Scholar

Crossref

39.

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

101

,

6168

(

1994

).

https://doi.org/10.1063/1.468399

Google Scholar

Crossref

40.

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

101

,

6184

(

1994

).

https://doi.org/10.1063/1.468400

Google Scholar

Crossref

41.

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

104

,

273

(

1996

).

https://doi.org/10.1063/1.470898

Google Scholar

Crossref

42.

G. J.

Martyna

,

J. Chem. Phys.

104

,

2018

(

1996

).

https://doi.org/10.1063/1.470958

Google Scholar

Crossref

43.

J.

Cao

and

G. J.

Martyna

,

J. Chem. Phys.

104

,

2028

(

1996

).

https://doi.org/10.1063/1.470959

Google Scholar

Crossref

44.

S.

Jang

and

G. A.

Voth

,

J. Chem. Phys.

111

,

2371

(

1999

).

https://doi.org/10.1063/1.479515

Google Scholar

Crossref

45.

T. D.

Hone

,

P. J.

Rossky

, and

G. A.

Voth

,

J. Chem. Phys.

124

,

154103

(

2006

).

https://doi.org/10.1063/1.2186636

Google Scholar

Crossref

PubMed

46.

A.

Pérez

,

M. E.

Tuckerman

, and

M. H.

Müser

,

J. Chem. Phys.

130

,

184105

(

2009

).

https://doi.org/10.1063/1.3126950

Google Scholar

Crossref

PubMed

47.

E. A.

Polyakov

,

A. P.

Lyubartsev

, and

P. N.

Vorontsov-Velyaminov

,

J. Chem. Phys.

133

,

194103

(

2010

).

https://doi.org/10.1063/1.3484490

Google Scholar

Crossref

PubMed

48.

I. R.

Craig

and

D. E.

Manolopoulos

,

J. Chem. Phys.

121

,

3368

(

2004

).

https://doi.org/10.1063/1.1777575

Google Scholar

Crossref

PubMed

49.

B. J.

Braams

and

D. E.

Manolopoulos

,

J. Chem. Phys.

125

,

124105

(

2006

).

https://doi.org/10.1063/1.2357599

Google Scholar

Crossref

PubMed

50.

S.

Habershon

,

D. E.

Manolopoulos

,

T. E.

Markland

, and

T. F.

Miller

III,

Annu. Rev. Phys. Chem.

64

,

387

(

2013

).

https://doi.org/10.1146/annurev-physchem-040412-110122

Google Scholar

Crossref

PubMed

51.

M.

Rossi

,

M.

Ceriotti

, and

D. E.

Manolopoulos

,

J. Chem. Phys.

140

,

234116

(

2014

).

https://doi.org/10.1063/1.4883861

Google Scholar

Crossref

PubMed

52.

T. E.

Markland

and

D. E.

Manolopoulos

,

J. Chem. Phys.

129

,

024105

(

2008

).

https://doi.org/10.1063/1.2953308

Google Scholar

Crossref

PubMed

53.

T. E.

Markland

and

D. E.

Manolopoulos

,

Chem. Phys. Lett.

464

,

256

(

2008

).

https://doi.org/10.1016/j.cplett.2008.09.019

Google Scholar

Crossref

54.

O.

Marsalek

and

T. E.

Markland

,

J. Chem. Phys.

144

,

054112

(

2016

).

https://doi.org/10.1063/1.4941093

Google Scholar

Crossref

PubMed

55.

C.

John

,

T.

Spura

,

S.

Habershon

, and

T. D.

Kühne

,

Phys. Rev. E

93

,

043305

(

2016

).

https://doi.org/10.1103/physreve.93.043305

Google Scholar

Crossref

PubMed

56.

R. P.

Steele

,

J.

Zwickl

,

P.

Shushkov

, and

J. C.

Tully

,

J. Chem. Phys.

134

,

074112

(

2011

).

https://doi.org/10.1063/1.3518714

Google Scholar

Crossref

PubMed

57.

M.

Takahashi

and

M.

Imada

,

J. Phys. Soc. Jpn.

53

,

3765

(

1984

).

https://doi.org/10.1143/jpsj.53.3765

Google Scholar

Crossref

58.

S.

Jang

,

S.

Jang

, and

G. A.

Voth

,

J. Chem. Phys.

115

,

7832

(

2001

).

https://doi.org/10.1063/1.1410117

Google Scholar

Crossref

59.

K.

Suzuki

,

M.

Tachikawa

, and

M.

Shiga

,

J. Chem. Phys.

132

,

144108

(

2010

).

https://doi.org/10.1063/1.3367724

Google Scholar

Crossref

PubMed

60.

A.

Pérez

and

M. E.

Tuckerman

,

J. Chem. Phys.

135

,

064104

(

2011

).

https://doi.org/10.1063/1.3609120

Google Scholar

Crossref

PubMed

61.

I.

Poltavsky

and

A.

Tkatchenko

,

Chem. Sci.

7

,

1368

(

2016

).

https://doi.org/10.1039/c5sc03443d

Google Scholar

Crossref

PubMed

62.

M.

Ceriotti

,

G.

Bussi

, and

M.

Parrinello

,

Phys. Rev. Lett.

103

,

030603

(

2009

).

https://doi.org/10.1103/PhysRevLett.103.030603

Google Scholar

Crossref

PubMed

63.

M.

Ceriotti

,

M.

Parrinello

,

T. E.

Markland

, and

D. E.

Manolopoulos

,

J. Chem. Phys.

133

,

124104

(

2010

).

https://doi.org/10.1063/1.3489925

Google Scholar

Crossref

PubMed

64.

M.

Ceriotti

,

D. E.

Manolopoulos

, and

M.

Parrinello

,

J. Chem. Phys.

134

,

084104

(

2011

).

https://doi.org/10.1063/1.3556661

Google Scholar

Crossref

PubMed

65.

M.

Ceriotti

and

D. E.

Manolopoulos

,

Phys. Rev. Lett.

109

,

100604

(

2012

).

https://doi.org/10.1103/physrevlett.109.100604

Google Scholar

Crossref

PubMed

66.

M. E.

Tuckerman

,

B. J.

Berne

, and

G. J.

Martyna

,

J. Chem. Phys.

97

,

1990

(

1992

).

https://doi.org/10.1063/1.463137

Google Scholar

Crossref

67.

M.

Praprotnik

,

L.

Delle Site

, and

K.

Kremer

,

J. Chem. Phys.

123

,

224106

(

2005

).

https://doi.org/10.1063/1.2132286

Google Scholar

Crossref

PubMed

68.

M.

Praprotnik

,

L.

Delle Site

, and

K.

Kremer

,

Annu. Rev. Phys. Chem.

59

,

545

(

2008

).

https://doi.org/10.1146/annurev.physchem.59.032607.093707

Google Scholar

Crossref

PubMed

69.

S.

Fritsch

,

S.

Poblete

,

C.

Junghans

,

G.

Ciccotti

,

L.

Delle Site

, and

K.

Kremer

,

Phys. Rev. Lett.

108

,

170602

(

2012

).

https://doi.org/10.1103/physrevlett.108.170602

Google Scholar

Crossref

PubMed

70.

A. B.

Poma

and

L.

Delle Site

,

Phys. Rev. Lett.

104

,

250201

(

2010

).

https://doi.org/10.1103/physrevlett.104.250201

Google Scholar

Crossref

PubMed

71.

A. B.

Poma

and

L.

Delle Site

,

Phys. Chem. Chem. Phys.

13

,

10510

(

2011

).

https://doi.org/10.1039/c0cp02865g

Google Scholar

Crossref

PubMed

72.

R.

Potestio

and

L.

Delle Site

,

J. Chem. Phys.

136

,

054101

(

2012

).

https://doi.org/10.1063/1.3678587

Google Scholar

Crossref

PubMed

73.

A.

Agarwal

and

L.

Delle Site

,

J. Chem. Phys.

143

,

094102

(

2015

).

https://doi.org/10.1063/1.4929738

Google Scholar

Crossref

PubMed

74.

A.

Agarwal

and

L.

Delle Site

,

Comput. Phys. Commun.

206

,

26

(

2016

).

https://doi.org/10.1016/j.cpc.2016.05.001

Google Scholar

Crossref

75.

L.

Delle Site

,

Phys. Rev. E

76

,

047701

(

2007

).

https://doi.org/10.1103/physreve.76.047701

Google Scholar

Crossref

76.

H.

Wang

,

C.

Hartmann

,

C.

Schütte

, and

L.

Delle Site

,

Phys. Rev. X

3

,

011018

(

2013

).

https://doi.org/10.1103/physrevx.3.011018

Google Scholar

Crossref

77.

A.

Agarwal

,

J.

Zhu

,

C.

Hartmann

,

H.

Wang

, and

L.

Delle Site

,

New J. Phys.

17

,

083042

(

2015

).

https://doi.org/10.1088/1367-2630/17/8/083042

Google Scholar

Crossref

78.

L.

Delle Site

and

M.

Praprotnik

,

Phys. Rep.

693

,

1

(

2017

).

https://doi.org/10.1016/j.physrep.2017.05.007

Google Scholar

Crossref

79.

K.

Kreis

,

M. E.

Tuckerman

,

D.

Donadio

,

K.

Kremer

, and

R.

Potestio

,

J. Chem. Theory Comput.

12

,

3030

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00242

Google Scholar

Crossref

PubMed

80.

R.

Cortes-Huerto

,

K.

Kremer

, and

R.

Potestio

,

J. Chem. Phys.

145

,

141103

(

2016

).

https://doi.org/10.1063/1.4964779

Google Scholar

Crossref

PubMed

81.

J. M.

Boereboom

,

R.

Potestio

,

D.

Donadio

, and

R. E.

Bulo

,

J. Chem. Theory Comput.

12

,

3441

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00205

Google Scholar

Crossref

PubMed

82.

M. E.

Tuckerman

and

D.

Marx

,

Phys. Rev. Lett.

86

,

4946

(

2001

).

https://doi.org/10.1103/physrevlett.86.4946

Google Scholar

Crossref

PubMed

83.

S.

Habershon

,

T. E.

Markland

, and

D. E.

Manolopoulos

,

J. Chem. Phys.

131

,

024501

(

2009

).

https://doi.org/10.1063/1.3167790

Google Scholar

Crossref

PubMed

84.

L.

Wang

,

M.

Ceriotti

, and

T. E.

Markland

,

J. Chem. Phys.

141

,

104502

(

2014

).

https://doi.org/10.1063/1.4894287

Google Scholar

Crossref

PubMed

85.

D. M.

Wilkins

,

D. E.

Manolopoulos

, and

L. X.

Dang

,

J. Chem. Phys.

142

,

064509

(

2015

).

https://doi.org/10.1063/1.4907554

Google Scholar

Crossref

PubMed

86.

T. F.

Miller

III and

D. E.

Manolopoulos

,

J. Chem. Phys.

122

,

184503

(

2005

).

https://doi.org/10.1063/1.1893956

Google Scholar

Crossref

PubMed

87.

R.

Potestio

,

S.

Fritsch

,

P.

Español

,

R.

Delgado-Buscalioni

,

K.

Kremer

,

R.

Everaers

, and

D.

Donadio

,

Phys. Rev. Lett.

110

,

108301

(

2013

).

https://doi.org/10.1103/physrevlett.110.108301

Google Scholar

Crossref

PubMed

88.

R.

Potestio

,

P.

Español

,

R.

Delgado-Buscalioni

,

R.

Everaers

,

K.

Kremer

, and

D.

Donadio

,

Phys. Rev. Lett.

111

,

060601

(

2013

).

https://doi.org/10.1103/physrevlett.111.060601

Google Scholar

Crossref

PubMed

89.

K.

Kreis

,

D.

Donadio

,

K.

Kremer

, and

R.

Potestio

,

Europhys. Lett.

108

,

30007

(

2014

).

https://doi.org/10.1209/0295-5075/108/30007

Google Scholar

Crossref

90.

P.

Español

,

R.

Delgado-Buscalioni

,

R.

Everaers

,

R.

Potestio

,

D.

Donadio

, and

K.

Kremer

,

J. Chem. Phys.

142

,

064115

(

2015

).

https://doi.org/10.1063/1.4907006

Google Scholar

Crossref

PubMed

91.

K.

Kreis

,

A. C.

Fogarty

,

K.

Kremer

, and

R.

Potestio

,

Eur. Phys. J.: Spec. Top.

224

,

2289

(

2015

).

https://doi.org/10.1140/epjst/e2015-02412-1

Google Scholar

Crossref

92.

K.

Kreis

and

R.

Potestio

,

J. Chem. Phys.

145

,

044104

(

2016

).

https://doi.org/10.1063/1.4959169

Google Scholar

Crossref

PubMed

93.

D.

Marx

and

M.

Parrinello

,

J. Chem. Phys.

104

,

4077

(

1996

).

https://doi.org/10.1063/1.471221

Google Scholar

Crossref

94.

M.

Tuckerman

,

D.

Marx

,

M.

Klein

, and

M.

Parrinello

,

J. Chem. Phys.

104

,

5579

(

1996

).

https://doi.org/10.1063/1.471771

Google Scholar

Crossref

95.

H. F.

Trotter

,

Proc. Am. Math. Soc.

10

,

545

(

1959

).

https://doi.org/10.1090/s0002-9939-1959-0108732-6

Google Scholar

Crossref

96.

G.

Strang

,

SIAM J. Numer. Anal.

5

,

506

(

1968

).

https://doi.org/10.1137/0705041

Google Scholar

Crossref

97.

W. R.

Rudnicki

,

B.

Lesyng

, and

S. C.

Harvey

,

Biopolymers

34

,

383

(

1994

).

https://doi.org/10.1002/bip.360340310

Google Scholar

Crossref

PubMed

98.

R.

Potestio

,

Eur. Phys. J. B

87

,

245

(

2014

).

https://doi.org/10.1140/epjb/e2014-50314-y

Google Scholar

Crossref

99.

B.

Leimkuhler

and

C.

Matthews

,

Appl. Math. Res. Express

2013

,

34

.

https://doi.org/10.1093/amrx/abs010

Crossref

100.

J.

Liu

,

D.

Li

, and

X.

Liu

,

J. Chem. Phys.

145

,

024103

(

2016

).

https://doi.org/10.1063/1.4954990

Google Scholar

Crossref

PubMed

101.

L.

Wang

,

S. D.

Fried

,

S. G.

Boxer

, and

T. E.

Markland

,

Proc. Natl. Acad. Sci. U. S. A.

111

,

18454

(

2014

).

https://doi.org/10.1073/pnas.1417923111

Google Scholar

Crossref

PubMed

102.

A.

Witt

,

S. D.

Ivanov

,

M.

Shiga

,

H.

Forbert

, and

D.

Marx

,

J. Chem. Phys.

130

,

194510

(

2009

).

https://doi.org/10.1063/1.3125009

Google Scholar

Crossref

PubMed

103.

T. J. H.

Hele

,

M. J.

Willatt

,

A.

Muolo

, and

S. C.

Althorpe

,

J. Chem. Phys.

142

,

191101

(

2015

).

https://doi.org/10.1063/1.4921234

Google Scholar

Crossref

PubMed

104.

J. D.

Halverson

,

T.

Brandes

,

O.

Lenz

,

A.

Arnold

,

S.

Bevc

,

V.

Starchenko

,

K.

Kremer

,

T.

Stuehn

, and

D.

Reith

,

Comput. Phys. Commun.

184

,

1129

(

2013

).

https://doi.org/10.1016/j.cpc.2012.12.004

Google Scholar

Crossref

105.

A. C.

Fogarty

,

R.

Potestio

, and

K.

Kremer

,

J. Chem. Phys.

142

,

195101

(

2015

).

https://doi.org/10.1063/1.4921347

Google Scholar

Crossref

PubMed

106.

K.

Kreis

,

R.

Potestio

,

K.

Kremer

, and

A. C.

Fogarty

,

J. Chem. Theory Comput.

12

,

4067

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00440

Google Scholar

Crossref

PubMed

107.

J. H.

Peters

,

R.

Klein

, and

L.

Delle Site

,

Phys. Rev. E

94

,

023309

(

2016

).

https://doi.org/10.1103/physreve.94.023309

Google Scholar

Crossref

PubMed

108.

F.

Ercolessi

and

J. B.

Adams

,

Europhys. Lett.

26

,

583

(

1994

).

https://doi.org/10.1209/0295-5075/26/8/005

Google Scholar

Crossref

109.

S.

Izvekov

and

G. A.

Voth

,

J. Chem. Phys.

123

,

134105

(

2005

).

https://doi.org/10.1063/1.2038787

Google Scholar

Crossref

PubMed

110.

S.

Izvekov

and

G. A.

Voth

,

J. Phys. Chem. B

109

,

6573

(

2005

).

https://doi.org/10.1021/jp0456685

Google Scholar

Crossref

PubMed

111.

D.

Reith

,

M.

Pütz

, and

F.

Müller-Plathe

,

J. Comput. Chem.

24

,

1624

(

2003

).

https://doi.org/10.1002/jcc.10307

Google Scholar

Crossref

PubMed

112.

J. D.

Weeks

,

D.

Chandler

, and

H. C.

Andersen

,

J. Chem. Phys.

54

,

5237

(

1971

).

https://doi.org/10.1063/1.1674820

Google Scholar

Crossref

113.

H. J. C.

Berendsen

,

J. R.

Grigera

, and

T. P.

Straatsma

,

J. Phys. Chem.

91

,

6269

(

1987

).

https://doi.org/10.1021/j100308a038

Google Scholar

Crossref

114.

Y.

Wu

,

H. L.

Tepper

, and

G. A.

Voth

,

J. Chem. Phys.

124

,

024503

(

2006

).

https://doi.org/10.1063/1.2136877

Google Scholar

Crossref

PubMed

115.

A. C.

Fogarty

and

D.

Laage

,

J. Phys. Chem. B

118

,

7715

(

2014

).

https://doi.org/10.1021/jp409805p

Google Scholar

Crossref

PubMed

116.

M.

Salvalaglio

,

C.

Perego

,

F.

Giberti

,

M.

Mazzotti

, and

M.

Parrinello

,

Proc. Natl. Acad. Sci. U. S. A.

112

,

E6

(

2015

).

https://doi.org/10.1073/pnas.1421192111

Google Scholar

Crossref

PubMed

117.

K.

Toukan

and

A.

Rahman

,

Phys. Rev. B

31

,

2643

(

1985

).

https://doi.org/10.1103/physrevb.31.2643

Google Scholar

Crossref

118.

J. G.

Kirkwood

,

J. Chem. Phys.

3

,

300

(

1935

).

https://doi.org/10.1063/1.1749657

Google Scholar

Crossref

119.

M.

Heidari

,

R.

Cortes-Huerto

,

D.

Donadio

, and

R.

Potestio

,

Eur. Phys. J.: Spec. Top.

225

,

1505

(

2016

).

https://doi.org/10.1140/epjst/e2016-60151-6

Google Scholar

Crossref

120.

D.

Eisenberg

and

W.

Kauzmann

,

The Structure and Properties of Water

(

Oxford University Press

,

New York

,

1969

).

Google Scholar

121.

H. E.

Stanley

and

N.

Ostrowzky

,

Correlations and Connectivity, Geometric Aspects of Physics, Chemistry, and Biology

(

Kluwer Academic

,

Dordrecht

,

1990

).

Google Scholar

Crossref

122.

J.

Teixeira

,

J. Phys. IV

3

,

C1-163

(

1993

).

https://doi.org/10.1051/jp4:1993114

Google Scholar

Crossref

123.

A.

Luzar

and

D.

Chandler

,

Nature

379

,

55

(

1996

).

https://doi.org/10.1038/379055a0

Google Scholar

Crossref

124.

A.

Luzar

and

D.

Chandler

,

Phys. Rev. Lett.

76

,

928

(

1996

).

https://doi.org/10.1103/physrevlett.76.928

Google Scholar

Crossref

PubMed

125.

A.

Luzar

,

J. Chem. Phys.

113

,

10663

(

2000

).

https://doi.org/10.1063/1.1320826

Google Scholar

Crossref

126.

A.

Luzar

,

Chem. Phys.

258

,

267

(

2000

).

https://doi.org/10.1016/s0301-0104(00)00163-4

Google Scholar

Crossref

127.

I.

Benjamin

,

J. Phys. Chem. B

109

,

13711

(

2005

).

https://doi.org/10.1021/jp044157f

Google Scholar

Crossref

PubMed

128.

N.

Winter

,

J.

Vieceli

, and

I.

Benjamin

,

J. Phys. Chem. B

112

,

227

(

2008

).

https://doi.org/10.1021/jp0734833

Google Scholar

Crossref

PubMed

129.

B. P.

Lambeth

,

C.

Junghans

,

K.

Kremer

,

C.

Clementi

, and

L.

Delle Site

,

J. Chem. Phys.

133

,

221101

(

2010

).

https://doi.org/10.1063/1.3522773

Google Scholar

Crossref

PubMed

130.

S.

Fritsch

,

C.

Junghans

, and

K.

Kremer

,

J. Chem. Theory Comput.

8

,

398

(

2012

).

https://doi.org/10.1021/ct200706f

Google Scholar

Crossref

PubMed

131.

A.

Agarwal

,

C.

Clementi

, and

L.

Delle Site

,

Phys. Chem. Chem. Phys.

19

,

13030

(

2017

).

https://doi.org/10.1039/c7cp01629h

Google Scholar

Crossref

PubMed

132.

D.

Marx

,

M. E.

Tuckerman

, and

G. J.

Martyna

,

Comput. Phys. Commun.

118

,

166

(

1999

).

https://doi.org/10.1016/S0010-4655(99)00208-8

Google Scholar

Crossref

133.

B.

Leimkuhler

,

D. T.

Margul

, and

M. E.

Tuckerman

,

Mol. Phys.

111

,

3579

(

2013

).

https://doi.org/10.1080/00268976.2013.844369

Google Scholar

Crossref

2017

Author(s)

From classical to quantum and back: Hamiltonian adaptive resolution path integral, ring polymer, and centroid molecular dynamics

I. INTRODUCTION

II. QUANTUM-CLASSICAL PATH INTEGRALS

III. QUANTUM-CLASSICAL PATH INTEGRAL MOLECULAR DYNAMICS

A. Evaluation of the adaptive mass and the resolution function on the centroids

B. Introducing normal modes

C. Introducing momenta

1. Constant kinetic masses (CKMs)

2. Adaptive kinetic masses (AKMs)

D. Equations of motion

E. Integration

F. Multiple time-stepping

G. Langevin thermostating

IV. APPROXIMATE QUANTUM DYNAMICS

A. Quantum-classical centroid molecular dynamics

B. Quantum-classical ring polymer molecular dynamics

V. SIMULATIONS

A. Water system

B. Setups

C. Free energy corrections

VI. RESULTS

A. Structure

B. Dynamics

C. Particle fluctuations

D. Computational speedup

VII. DISCUSSION AND CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

From classical to quantum and back: Hamiltonian adaptive resolution path integral, ring polymer, and centroid molecular dynamics

I. INTRODUCTION

II. QUANTUM-CLASSICAL PATH INTEGRALS

III. QUANTUM-CLASSICAL PATH INTEGRAL MOLECULAR DYNAMICS

A. Evaluation of the adaptive mass and the resolution function on the centroids

B. Introducing normal modes

C. Introducing momenta

1. Constant kinetic masses (CKMs)

2. Adaptive kinetic masses (AKMs)

D. Equations of motion

E. Integration

F. Multiple time-stepping

G. Langevin thermostating

IV. APPROXIMATE QUANTUM DYNAMICS

A. Quantum-classical centroid molecular dynamics

B. Quantum-classical ring polymer molecular dynamics

V. SIMULATIONS

A. Water system

B. Setups

C. Free energy corrections

VI. RESULTS

A. Structure

B. Dynamics

C. Particle fluctuations

D. Computational speedup

VII. DISCUSSION AND CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only