Enhanced sampling techniques that target a set of collective variables and that use molecular dynamics as the driving engine have seen widespread application in the computational molecular sciences as a means to explore the free-energy landscapes of complex systems. The use of molecular dynamics as the fundamental driver of the sampling requires the introduction of a time step whose magnitude is limited by the fastest motions in a system. While standard multiple time-stepping methods allow larger time steps to be employed for the slower and computationally more expensive forces, the maximum achievable increase in time step is limited by resonance phenomena, which inextricably couple fast and slow motions. Recently, we introduced deterministic and stochastic resonance-free multiple time step algorithms for molecular dynamics that solve this resonance problem and allow ten- to twenty-fold gains in the large time step compared to standard multiple time step algorithms [P. Minary et al., Phys. Rev. Lett. 93, 150201 (2004); B. Leimkuhler et al., Mol. Phys. 111, 3579–3594 (2013)]. These methods are based on the imposition of isokinetic constraints that couple the physical system to Nosé-Hoover chains or Nosé-Hoover Langevin schemes. In this paper, we show how to adapt these methods for collective variable-based enhanced sampling techniques, specifically adiabatic free-energy dynamics/temperature-accelerated molecular dynamics, unified free-energy dynamics, and by extension, metadynamics, thus allowing simulations employing these methods to employ similarly very large time steps. The combination of resonance-free multiple time step integrators with free-energy-based enhanced sampling significantly improves the efficiency of conformational exploration.

Molecular dynamics (MD)-based enhanced sampling approaches have become a staple in the computational molecular sciences as a means to explore rough energy landscapes of complex systems. Because MD allows a system to be studied under thermal conditions, use of these methods allows for efficient determination and ranking of preferred configurations based on relative free energies. One class of MD-based schemes involves targeting a pre-selected set of collective variables (CVs) for enhanced sampling, and in this class, the most popular technique remains the metadynamics algorithm1 and its variants.2–6 Metadynamics is essentially a basin-filling scheme7 in which Gaussians are deposited in the CV space as the simulation proceeds in order to fill in local free-energy minima in order to promote barrier crossing and prevent the system from revisiting previously explored regions of configuration space. Another useful MD-driven scheme is the adiabatic and driven-adiabatic free-energy dynamics (AFED, d-AFED)8,9/temperature-accelerated MD (TAMD),10 which uses high temperature and adiabatic decoupling to drive the sampling of the target CVs. The TAMD/d-AFED and metadynamics methods have recently been combined to produce the unified free-energy dynamics (UFED) approach11 in which metadynamics is used as an additional bias to TAMD/d-AFED and the temperature-accelerated sliced sampling,12 which employs metadynamics, TAMD/d-AFED, and umbrella sampling on different subsets of CVs in a given simulation. All of these approaches have proved remarkably successful in providing new insights into the thermodynamics of protein folding,13,14 crystal structure prediction,15–20 first-order phase transitions,21,22 and the screening of different poses of bound inhibitor molecules in candidate drug-binding problems,23 to name just a few examples. While in the majority of CV-based enhanced sampling applications, the CVs are chosen a priori based largely on chemical and/or physical intuition, approaches that allow CVs to be learned “on the fly” as a simulation proceeds are now being introduced.24 Because the aforementioned enhanced sampling approaches are driven by MD, their efficiency can be further increased by employing multiple time step (MTS) integration algorithms, which assign a small time step to the most rapidly varying force components and larger time steps to more slowly varying force components. The most commonly used method is the reversible reference-system propagator algorithm or r-RESPA.25 In their original form, MTS algorithms allow for modest reductions in time to the solution of a given problem. However, their performance is limited by resonance phenomena, which cause motion on the fast time scales to limit the time step that can be employed for the slow force components.26–28 Several solutions to the resonance problem have since been proposed, which include the use of colored-noise thermostatting29 and novel combinations of the isokinetic ensemble30–33 with Nosé-Hoover chain thermostats34 and Nosé-Hoover Langevin thermostats.35 The isokinetic schemes,36–39 in particular, allow very large time steps to be employed and are especially useful when equilibrium sampling is required. Moreover, unlike the colored-noise thermostat approach,29a priori parameterization of memory kernels is not needed. It is important to note, however, that the colored-noise scheme largely preserves the dynamics if the time step is not chosen too large, while the isokinetic method largely sacrifices quantitative dynamical information.

In this paper, we propose a version of the TAMD/d-AFED scheme in which two isokinetic constraints, one applied to physical phase-space variables at physical temperature and the other applied in an extended phase-space corresponding to the CVs at high temperature, are used in conjunction with MTS integration allowing enhanced sampling calculations to be performed with very large time steps. A primary reason for combining TAMD/d-AFED with a resonance-free MTS approach is that, as is described in Sec. II, stiff harmonic forces on the CVs, in addition to those already present in the chosen force field, lead to new resonant terms in the equations of motion. The approach we describe in this paper eliminates all of these resonances, thereby permitting the use of the aforementioned very large time steps. We demonstrate the performance of the new scheme on fully atomistic simulations of the solvated alanine di- and tri-peptides and show that time steps as large as 100 fs can be employed. We further show that the same approach can be extended to the UFED algorithm, suggesting that metadynamics calculations could also be performed with very large time steps using the isokinetic approach.

This paper is organized as follows. In Sec. II, we review the theoretical aspects of the double isokinetic TAMD/d-AFED and UFED approaches. In Secs. III and IV, we present the aforementioned test cases. Conclusions are given in Sec. V, and further analysis is reserved for the accompanying  Appendix.

Molecular dynamics (MD)-based enhanced sampling techniques targeting a pre-selected set of collective variables (CVs) can be generally formulated as follows. We consider a set of N particles with coordinates r1, …, rNr in a volume V at temperature T interacting via a potential U(r). Suppose a set of n < N CVs, denoted θα(r), α = 1, …, n, has been chosen to distinguish the important conformational states on the free-energy landscape of the system. The free-energy surface as a function of these CVs is obtained first by defining the so-called marginal probability distribution function at temperature T,

P(s1,,sn)P(s)=D(V)dNreβU(r)α=1nδθα(r)sα,
(1)

where β = 1/kBT and D(V) is the spatial domain defined by the confining volume. The Dirac δ-functions require that the CVs θα(r) have values sα, known as the coarse-grained variables (CGVs). The Helmholtz free-energy surface A(s) is then given by A(s) = −kBT ln P(s). Equation (1) suggests that an n-dimensional grid be introduced in the CGV space and that the marginal distribution P(s) be computed at each point of this grid via a separate simulation. This, in fact, is how the well-known bluemoon ensemble40,41 approach works and is feasible if n = 1 or perhaps, as a tour de force, if n = 2. However, for larger values of n, gridding up P(s) becomes infeasible.

A viable alternative is to sample P(s) in the same spirit that one samples the Boltzmann distribution exp(−βU(r)), e.g., via MD. However, unlike the Boltzmann distribution, the marginal P(s) is not known a priori, thus raising the question of how an unknown probability distribution can be sampled if its functional form is unknown. This is what the TAMD/d-AFED method accomplishes. In order to see how this is done, we begin by replacing the δ-functions by Gaussians as

δθα(r)sα=limκαβκα2π1/2eβκαθα(r)sα2/2.
(2)

This replacement allows the exponential in the Gaussian to be incorporated into the Boltzmann distribution as an additional harmonic potential coupling the CVs to the CGVs. In this way, the physical coordinates and Gaussian centers s can be sampled simultaneously via an extended phase-space Hamiltonian of the form

H(r,p,s,π)=i=1Npi22mi+α=1nπα22μα+U(r)+α=1n12καθα(r)sα2,
(3)

where the momenta p1, …, pNp are conjugate to the physical coordinates, m1, …, mN are the particle masses, π1, …, πnπ are the momenta conjugate to the CGVs, and μα are a set of fictitious masses for the CGVs. Note that the spring constants κ1, …, κnκ are approximated as finite in Eq. (3). Simple use of Eq. (3) does not guarantee that the unknown marginal P(s) will be sampled. However, if we also affect an adiabatic decoupling between the CGVs and the physical variables by choosing the fictitious masses μα sufficiently large, then it can be proved8,9 that the marginal distribution is properly sampled. Introducing the adiabatic decoupling8,9,10,42,66 also allows the sampling to be accelerated via the introduction of a high temperature TsT for the CGVs, allowing for facile barrier crossing. In Refs. 8, 9, and 42, it was shown that the combination of high-temperature and adiabatic decoupling generates a distribution Pκ,adb(s) related to the true marginal P(s) via

limκPκ,adb(s)=P(s)T/Ts,
(4)

which we term the adiabatic marginal distribution. Note that a similar relation also emerges in the context of the so-called well-tempered metadynamics method.2,3 Thus, the free energy at temperature T can, in principle, be obtained via

A(s)=limκkBTslnPκ,adb(s).
(5)

In Refs. 11 and 43, however, it was shown that use of free-energy gradient information, obtained by averages of the harmonic force καθα(r)sα, leads to greater accuracy and faster convergence of the free-energy surface A(s).

The adiabatic marginal distribution can be generated using standard MD equations of motion coupled to two separate heat baths

mir̈i=Fi(r)α=1nκαθα(r)sαθαri+Bath(T),μαs̈α=καθα(r)sα+Bath(Ts),
(6)

where Fi(r) = −∂U/ri is the force from the potential.44 The bath coupling can be stochastic or deterministic, and in previous applications of this scheme, we have made use of the robust generalized Gaussian moment thermostat (GGMT).45 

The efficiency of these equations of motion in generating the adiabatic marginal distribution can be increased when they are integrated using multiple time-stepping algorithms.25,46,59 However, as noted in the Introduction, resonance phenomena26–28 place a fundamental limit on the magnitude of the largest time step that can be employed. For typical atomistic force fields, this limit is around 5-10 fs.47 The stochastic isokinetic Nosé-Hoover algorithm37–39 and its deterministic Nosé-Hoover chain antecedent36 circumvent resonances by preventing energy from building up in any mode in a system through the imposition of isokinetic constraints on the dynamics. Adapting this method for TAMD/d-AFED requires imposition of different isokinetic constraints on the physical and extended phase-spaces. Let qk denote any one of the 3N physical Cartesian degrees of freedom, and let vk=q̇k denote its velocity. Each of these physical variables is coupled to L additional pairs of Nosé-Hoover Langevin variables denoted v1,k,j and v2,k,j, where j = 1, …, L. Although for most properties of interest, it is found that L = 1 is likely sufficient,37 we, nevertheless, present the formalism for arbitrary L, as particular cases could arise in which employing L > 1 could be useful (see Ref. 38 for one such example). Similarly, let να = πα/μα = ṡα denote the velocity for sα. Each of these is coupled to L pairs of Nosé-Hoover Langevin variables denoted ν1,α,j and ν1,α,j, where j = 1, …, L. The “double” stochastic isokinetic equations of motion for TAMD/d-AFED then take the forms

dqk=vkdt,dvk=Fk(r)mkdtλkvkdt,dv1,k,j=λkv1,k,jdtv2,k,jv1,k,jdt,dv2,k,j=G(v1,j,k)Q2dtγkv2,k,jdt+σkdwk,dsα=ναdt,dνα=καμαθα(r)sαdtηαναdt,dν1,α,j=ηαν1,α,jdtν2,α,jν1,α,jdt,dν2,α,j=G(ν1,α,j)Q2dtΓαν2,α,jdt+ΣαdWα.
(7)

Here, Fk(r) is the kth force component, G(y) = Q1y2kBT, G(y)=Q1y2kBTs, Q1 and Q2 and Q1 and Q2 are time scale parameters for the two sets of Nosé-Hoover Langevin thermostats on the physical and extended phase-space variables, respectively, σk=2γk/βQ2, and Σα=2Γα/βsQ2, where, βs = 1/kBTs. The stochastic processes dwk and dWα are taken to be Gaussian random processes for the physical and extended phase-space variables, respectively. Finally, the parameters γk and Γα are friction constants and λk and ηα are Lagrange multipliers for enforcing the following constraint conditions:

mkvk2+LL+1j=1LQ1v1,k,j2=LkBT,μανα2+LL+1j=1LQ1ν1,α,j2=LkBTs.
(8)

Exact analytical expressions for these multipliers can be obtained using Gauss’ principle of least constraint,30,31,33,48 yielding

λk=vkFk(r)LL+1j=1LQ1v1,k,j2v2,k,jmkvk2+LL+1j=1LQ1v1,k,j2,ηα=ναFα(r,sα)LL+1j=1LQ1ν1,α,j2ν2,α,jμανα2+LL+1j=1LQ1ν1,α,j2,
(9)

where Fα(r,sα)=κα(θα(r)sα). In the  Appendix, we provide an analysis of these equations and show that they reproduce the correct free-energy surface A(s) at the physical temperature T. In Ref. 37, we devised MTS integrators for these equations of motion. These integrators involve three time steps, the smallest being assigned to bonded interactions, the next to short-range non-bonded interactions, and the largest or “outer” time step to long-range non-bonded interactions. When electrostatic interactions are treated via Ewald summation, two schemes exist for dividing up the sum. In one of these, known as RESPA1, both the real-space and reciprocal-space sums are divided into short-range and long-range contributions, the latter according to the magnitude of the corresponding reciprocal-space vectors. In the second scheme, known as RESPA2, the entire reciprocal-space sum is assigned to the outer time step, which introduces some short-range character at this largest time step.37,49 In order to correct for this, a short-range correction is subtracted out in real space. Although the cancellation is not perfect, with a correct choice of parameters, the error in this cancellation can be made negligible. The scheme of Procacci et al.50 could presumably also be employed to improve the accuracy of the correction. Finally, in Refs. 38 and 39, we show how polarization forces in polarizable models of the AMOEBA type can be divided based on short-range and long-range self-consistent polarization calculations. For all MTS schemes, we term the algorithm the “stochastic isokinetic Nosé-Hoover RESPA” or SIN(R) algorithm. When the Nosé-Hoover variables are coupled to the shortest time step, the scheme is referred to as a XI-RESPA (eXtended-Inner-RESPA) scheme, if they are coupled to the intermediate time step, it is referred to as XM-RESPA (eXtended-Middle-RESPA), and if they are coupled to the outer time step, we refer to the method as XO-RESPA (eXtended-Outer-RESPA), following the nomenclature in Ref. 46. In the unified free-energy dynamics (UFED) approach,11 TAMD/d-AFED is supplemented with a bias potential Ubias(s, t) derived from the metadynamics scheme.1 This bias takes the form of a sum of Gaussian functions of height h and width σ added to the CGVs at specific times tl < t, where t is the current time of the simulation, i.e.,

Ubias(s,t)=hl,tl<tess(tl)2/2σ2,
(10)

where ∥⋅∥ is the usual L2 norm. This bias has the effect of filling in basins as they are visited during the TAMD/d-AFED run, thereby further accelerating barrier crossing. The bias serves only this purpose and is not ultimately used in the construction of the FES. In order to implement this bias, one only needs to add a force term −∂Ubias/∂sα into the equation of motion for sα in Eqs. (6) or να in Eqs. (7). In actual applications, evaluating the sum in Eq. (10) for a long trajectory can become computationally expensive. In order to circumvent this, the bias is mapped onto a regular grid, and an interpolation algorithm is used to calculate the bias potential and its gradient at arbitrary values of the CGVs.

All the simulations were performed using the PINY_MD package51 on fully solvated alanine di- and tripeptide systems. We have chosen the CHARMM22 force field,52 as we have benchmark data using this force field for both of these peptides.11 The CVs are taken to be the two backbone Ramachandran dihedral angles Φ and Ψ for the dipeptide and the four corresponding Ramachandran angles for the tripeptide. The mass-like parameters μΦ,Ψ were assigned values 168.0 amu Å2/rad2, and the harmonic coupling constants κΦ,Ψ were both set to 2.78 × 103 (kcal/mol)/rad2. The temperature T was set as 300 K for the physical system. The extended system was maintained at temperature Ts = 1000 K for UFED and at Ts = 1500 K for d-AFED. In all the simulations, the values of Q1 and Q2 were determined by the relation Q1,2 = kBTsτ2, where τ = 80 fs. In double isokinetic simulations, the SIN(R) thermostat was used for both the physical and the extended systems with γ = 0.03 ps−1. The Nosé-Hoover Langevin part of the algorithm is integrated with a Suzuki-Yoshida scheme with nsy = 3, nres = 2,46,53–55 and L is taken to be 1. The benchmark results were the converged free energy surfaces (FESs) obtained using the XO-RESPA UFED method with a time step of δt = 0.2 fs for the intramolecular interactions and an intermolecular time step of Δt = 1.0 fs. A Nosé-Hoover chain thermostat34 was used for the physical system, and a generalized Gaussian moment thermostat (GGMT)45,56 was used for the extended system for the benchmark simulations. The flexible SPC/E water model57 was used for explicit solvation. Electrostatic interactions were treated using the smooth particle-mesh Ewald (SPME)58 method.

For simulations using the XI-, XM-, and XO-RESPA algorithms,46,60 time steps of 0.5 fs, 1.0 fs, and 5.0 fs were used for the intramolecular interactions and harmonic coupling between the CVs and CGVs, the torsion angle interactions, and the short-range interactions, respectively. With this choice, we then varied the largest or “outermost” time step in order to determine how large this time step could be chosen for the expensive long-range interactions.

The L1 error of the FESs was calculated in the following way:

ζ=1Nbi=1Nb|A(si)Aconv(si)|,
(11)

where Nb is the number of bins used to represent the free-energy surface. si ≡ ((s1,i, …, sn,i)) is a point in the extended space. The free energy A(si) was compared to the benchmark free energy Aconv(si) at that point.

An alanine dipeptide solvated in a cubic box of 216 water molecules was equilibrated consecutively at constant volume for 100 ps, at a constant pressure of 1 atm for 100 ps, and at the obtained equilibrium volume for 100 ps. The real space truncation was 9.25 Å, corresponding to half the box length obtained in the NPT run. In the XI-RESPA1 and XM-RESPA1 simulations, the short-range cutoff was 6.0 Å. A healing length of 3.0 Å was employed in order to smoothly switch off the short-range interactions.

In the XI-RESPA2 and XM-RESPA2 simulations, the short-range cutoff was chosen to be 6.8 Å. A healing length of 3.0 Å was chosen, meaning that the switch begins at 3.8 Å. The system was equilibrated for 100 ps. The Gaussian height in the UFED simulations was taken to be 0.2 kcal/mol. The Gaussian width and deposition rate were 0.2 rad and 0.01 fs−1. A basis-set expansion method was used to construct the FESs.61,62 In particular, the FESs were expanded in a basis of 256 plane waves on a 16 × 16 Fourier grid (see Fig. 2 of Ref. 11 for the Ramachandran free-energy surface).

An alanine tripeptide was solvated in a cubic box of 512 water molecules. NVT, NPT, and NVT simulations, each 200 ps in length, were performed in succession in order to equilibrate the system, as in the case of the alanine dipeptide. A real space truncation of 12.0 Å was employed. In the XI-RESPA1 and XM-RESPA1 simulations, the short-range cutoff was 6.0 Å and the healing length was 3.0 Å. In the XI-RESPA2 and XM-RESPA2 simulations, the short-range cutoff was chosen to be 8.0 Å. The healing length was 3.0 Å, meaning that the switch begins at 3.0 Å for RESPA1 and 5.0 Å for RESPA2. In the UFED simulations, the Gaussian height and width were 2.0 kcal/mol and 1.25 rad, respectively. The deposition rate was 0.01 fs−1. The FESs were constructed with a sparse binning scheme.11 

As a simple example, we chose a two-variable system with a potential of the form

V(x,y)=D0(x2a2)2+12κy2+λxy,
(12)

where D0 = 5, a = 1, κ = 1, and λ = 2.878. In these simulations, the earlier adiabatic free energy dynamics (AFED) scheme was used.8 The mass and temperature of y were my = 1 and Ty = 1. In the AFED simulations, the temperature for x was Tx = 10. The mass of x was mxmy. Each degree of freedom (x and y) was coupled to a separate thermostat. In double isokinetic simulations, the SIN(R) thermostat was used and was compared to standard AFED with a GGMT thermostat. The thermostat parameters were the Suzuki-Yoshida factor nsy = 3, the RESPA factor nres = 2, and L = 1. For all simulations, the total number of steps simulated was the same, i.e.,1, 000, 000/Δt, which clearly changes depending on the choice of Δt. In simulations using XI-RESPA, since mxmy, the slow and fast forces on x and y are, respectively,

fx,slow(x)=4D0(x2a2)x,fy,slow(x)=λx.
(13)

These were updated at every outer time step Δt. The fast forces

fx,fast(y)=λy,fy,fast(y)=κy
(14)

were updated every inner time step δt. The probability distribution of P(x) was compared with the analytical probability distribution Pan(x) using the L1 error measure

ζ=1Nbi=1Nb|P(xi)Pan(xi)|,
(15)

where Nb was the number of bins used to generate the distribution. Figure 1 shows that the L1 error decreases with increasing mx. This result is consistent with previous AFED/d-AFED studies.8,9 The adiabaticity condition is only met when the mass of x is sufficiently large, for example, mx ≥ 300. Since further increase of mx has little influence of the L1 error, we chose mx = 300 in the proceeding studies.

FIG. 1.

L1 error of the probability distribution P(x) for a double well coupled to a harmonic oscillator in Eq. (12). The figure shows the convergence of the probability distribution as a function of mx with Tx = 10. In all cases, Ty = 1, my = 1, and the time step Δt = 0.000 25. The analytical probability distribution Pan(x) is the benchmark.

FIG. 1.

L1 error of the probability distribution P(x) for a double well coupled to a harmonic oscillator in Eq. (12). The figure shows the convergence of the probability distribution as a function of mx with Tx = 10. In all cases, Ty = 1, my = 1, and the time step Δt = 0.000 25. The analytical probability distribution Pan(x) is the benchmark.

Close modal

The comparison of the AFED method with and without XI-RESPA is shown in Fig. 2. The inner time step is δt = 0.000 25, and the outer time step Δt is varied. We see that the error is nearly constant in the simulations with RESPA. The SIN(R) scheme allows the use of very large time steps for the slow force with minimal residual error.

FIG. 2.

L1 error of the probability distribution P(x) for a double well coupled to a harmonic oscillator in Eq. (12) for different time steps Δt. For XI-RESPA methods, the inner time step δt is fixed at 0.000 25 and Δt is the outer time step. In all cases, Tx = 10, Ty = 1, mx = 300, and my = 1. The analytical probability distribution Pan(x) is the benchmark.

FIG. 2.

L1 error of the probability distribution P(x) for a double well coupled to a harmonic oscillator in Eq. (12) for different time steps Δt. For XI-RESPA methods, the inner time step δt is fixed at 0.000 25 and Δt is the outer time step. In all cases, Tx = 10, Ty = 1, mx = 300, and my = 1. The analytical probability distribution Pan(x) is the benchmark.

Close modal

The sensitivity of the error to the values of thermostat parameters Q1, Q2, and γ was also studied (Fig. 3). For a wide range of values of these parameters, only when the values of γ and Q2 are either too small or too large is the error seen to become large. For most values of these parameters, the L1 error is small and roughly constant.

FIG. 3.

L1 error of the probability distribution P(x) for a double well coupled to a harmonic oscillator in Eq. (12) for different values of Q1, Q2, and γ. These are double isokinetic AFED results. In all cases, Δt = 0.1, Tx = 10, Ty = 1, mx = 300, and my = 1. The analytical probability distribution Pan(x) is the benchmark. Red and blue dots separate values of ζ into two regions: ζ ≤ 0.01 is designated with red dots, and ζ > 0.01 is designated with blue dots.

FIG. 3.

L1 error of the probability distribution P(x) for a double well coupled to a harmonic oscillator in Eq. (12) for different values of Q1, Q2, and γ. These are double isokinetic AFED results. In all cases, Δt = 0.1, Tx = 10, Ty = 1, mx = 300, and my = 1. The analytical probability distribution Pan(x) is the benchmark. Red and blue dots separate values of ζ into two regions: ζ ≤ 0.01 is designated with red dots, and ζ > 0.01 is designated with blue dots.

Close modal

In simulations that do not employ multiple time-stepping, we have tested the largest time step can be used within the double isokinetic d-AFED and UFED methods in the aqueous phase. The results are shown in Fig. 4, which shows the L1 error in the two-dimensional free-energy surface. In particular, note that this L1 error is computed using either the full free-energy surface or using only those bins whose free-energy values lie below a chosen maximum free-energy cutoff. In the latter comparison, two cutoff values, 7.5 kcal/mol and 5.0 kcal/mol, were chosen. In Fig. 4, we see that the L1 error begins to increase when the time step is larger than 3 fs. In the d-AFED simulation, the convergence of the high free energy area is slow, which can be seen from the fact that the L1 error changes depending on the choice of the free-energy cutoff. By comparison, the UFED method improves convergence in a shorter simulation time, consistent with what was observed in Ref. 11.

FIG. 4.

L1 error of FESs of alanine dipeptide in aqueous phase as a function of time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 10 ns, 0.2 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature for extended variables Ts in the d-AFED simulations are 20 ns and 1500 K, respectively. The benchmark is the converged FES of 15 ns XO-RESPA UFED simulation with Gaussian height h = 0.02 kcal/mol.

FIG. 4.

L1 error of FESs of alanine dipeptide in aqueous phase as a function of time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 10 ns, 0.2 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature for extended variables Ts in the d-AFED simulations are 20 ns and 1500 K, respectively. The benchmark is the converged FES of 15 ns XO-RESPA UFED simulation with Gaussian height h = 0.02 kcal/mol.

Close modal

We next apply the XI-RESPA1 and XI-RESPA2 and XM-RESPA1 and XM-RESPA2 algorithms to the double isokinetic d-AFED and UFED methods for the solvated alanine dipeptide. The results are shown in Figs. 5 and 6. The XI-RESPA and XM-RESPA double isokinetic schemes generally work well for this system. Regions having free energies lower than 5.0 kcal/mol are in good agreement with the benchmark. The FESs generated from UFED converge better than those from d-AFED, again consistent with Ref. 11. Importantly, the error in the XI-RESPA1 results remains nearly constant for different choices of the outer time step. For the XI-RESPA2 result, the error tends to increase more when the outer time step increases. The reason for this is likely that the convergence of the XI-RESPA2 and XM-RESPA2 FESs is more sensitive to the particular choice of the RESPA simulation parameters.37 However, the UFED results are much less sensitive to the details of the RESPA protocol. As a test of the dependence on the number L of stochastic Nosé-Hoover baths, we also generated XI-RESPA2 double isokinetic d-AFED free-energy surfaces for the solvated alanine dipeptide with L = 4. The errors (curves not shown), which should be compared to those in the third panel of Fig. 5, are generally reduced by approximately 0.1 kcal/mol or less, depending on the free-energy cutoff. The potential energies of the physical system and the extended system have clearly convergent averages over the trajectory as seen in Fig. 7, which indicates that the XI-RESPA double isokinetic method does not suffer from energy accumulation when a large outer time step is used.

FIG. 5.

L1 error of FESs of alanine dipeptide in aqueous phase as a function of outer time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 10 ns, 0.2 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature for extended variables Ts in the d-AFED simulations are 20 ns and 1500 K, respectively. The benchmark is the converged FES of 15 ns XO-RESPA UFED simulation with Gaussian height h = 0.02 kcal/mol.

FIG. 5.

L1 error of FESs of alanine dipeptide in aqueous phase as a function of outer time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 10 ns, 0.2 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature for extended variables Ts in the d-AFED simulations are 20 ns and 1500 K, respectively. The benchmark is the converged FES of 15 ns XO-RESPA UFED simulation with Gaussian height h = 0.02 kcal/mol.

Close modal
FIG. 6.

L1 error of FESs of alanine dipeptide in aqueous phase as a function of outer time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 10 ns, 0.2 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature for extended variables Ts in the d-AFED simulations are 20 ns and 1500 K, respectively. The benchmark is the converged FES of 15 ns XO-RESPA UFED simulation with Gaussian height h = 0.02 kcal/mol.

FIG. 6.

L1 error of FESs of alanine dipeptide in aqueous phase as a function of outer time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 10 ns, 0.2 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature for extended variables Ts in the d-AFED simulations are 20 ns and 1500 K, respectively. The benchmark is the converged FES of 15 ns XO-RESPA UFED simulation with Gaussian height h = 0.02 kcal/mol.

Close modal
FIG. 7.

The instant and average potential energy (PE) of the physical and extended system after 5 ns XI-RESPA1 double isokinetic d-AFED simulation of alanine dipeptide in aqueous phase. The outer time step Δt is 100 fs and the temperature for extended variables is 1500 K.

FIG. 7.

The instant and average potential energy (PE) of the physical and extended system after 5 ns XI-RESPA1 double isokinetic d-AFED simulation of alanine dipeptide in aqueous phase. The outer time step Δt is 100 fs and the temperature for extended variables is 1500 K.

Close modal

The computational cost of long-range interactions is rather high for solvated systems. We have investigated the simulation time saved by using the XI-RESPA and XM-RESPA schemes within the particular implementation in the PINY_MD package.51 The benchmark is the average CPU time tbench per time step δt = 0.5 fs of the double isokinetic d-AFED simulation without RESPA. The time ratio r = tbench/tres, where tres is the average CPU time of 0.5 fs XI-RESPA simulation. The time saved by using XI-RESPA is roughly 10 for the largest time step of Δt = 100 fs. There is negligible difference in the time saving between the d-AFED and UFED methods, which means that the extra time spent on calculating forces from the bias potential is negligible compared to the total simulation time. The corresponding time saved using XM-RESPA is around 20-25, i.e., roughly a factor of 2 better than XI-RESPA. This is similar to the gains obtained in Ref. 38 within the TINKER package.63 Generally, however, gains associated with a particular implementation will depend sensitively on the details of the implementation.

We chose the alanine tripeptide in aqueous solution, a larger system with 4 CVs, to demonstrate the consistency of the performance of the XI-RESPA double isokinetic d-AFED/UFED. Figure 8 shows the L1 error of different outer time steps. Although the high free-energy area of the 4D FES converges relatively slowly, regions of the free-energy surface below 5.0 kcal/mol are correctly mapped by both the d-AFED and the UFED methods within 0.5 kcal/mol. Since regions of interest are usually of relatively low free energy, the results illustrate that the XI-RESPA double isokinetic d-AFED/UFED can be applied to multiple CVs and that the accuracy achieved using a large outer time step is comparable to the benchmark.

FIG. 8.

L1 error of FESs of alanine tripeptide in aqueous phase as a function of outer time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 50 ns, 2.0 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature of extended variables Ts in the d-AFED simulations are 120 ns, and 1500 K, respectively. The benchmark is the FES of 270 ns XO-RESPA UFED simulation with Gaussian height h = 0.2 kcal/mol.

FIG. 8.

L1 error of FESs of alanine tripeptide in aqueous phase as a function of outer time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 50 ns, 2.0 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature of extended variables Ts in the d-AFED simulations are 120 ns, and 1500 K, respectively. The benchmark is the FES of 270 ns XO-RESPA UFED simulation with Gaussian height h = 0.2 kcal/mol.

Close modal

Figure 9 shows the L1 error for different outer time steps of the XM-RESPA double isokinetic d-AFED/UFED. The CPU time ratio for the PINY_MD implementation is, again, approximately 10 for XI-RESPA and 20-25 for XM-RESPA.

FIG. 9.

L1 error of FESs of alanine tripeptide in aqueous phase as a function of outer time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 50 ns, 2.0 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature of extended variables Ts in the d-AFED simulations are 120 ns and 1500 K, respectively. The benchmark is the FES of 270 ns XO-RESPA UFED simulation with Gaussian height h = 0.2 kcal/mol.

FIG. 9.

L1 error of FESs of alanine tripeptide in aqueous phase as a function of outer time step Δt. The values of simulation time, Gaussian height h, and temperature of extended variables Ts in the UFED simulations are 50 ns, 2.0 kcal/mol, and 1000 K, respectively. The values of simulation time and temperature of extended variables Ts in the d-AFED simulations are 120 ns and 1500 K, respectively. The benchmark is the FES of 270 ns XO-RESPA UFED simulation with Gaussian height h = 0.2 kcal/mol.

Close modal

In this paper, we have shown how to combine resonance-free multiple time step integration, SIN(R), with the free-energy-based enhanced sampling methods, driven-adiabatic free-energy dynamics (d-AFED)/temperature-accelerated molecular dynamics (TAMD) and unified free-energy dynamics (UFED). Since metadynamics is a special case of UFED where Ts=T, the double isokinetic approach also, therefore, extends to metadynamics. We expect that the extension to the well-tempered metadynamics approach2,3 would be straightforward and would yield similar gains in time step. The combined scheme was applied to the test cases of the solvated alanine di- and tripeptides. Two RESPA force divisions, termed RESPA1 and RESPA2,37,49 which differ in how the reciprocal-space part of the Ewald sum is handled, were investigated. For the dipeptide, the two Ramachandran angles ϕ and ψ were targeted as collective variables (CVs), and it was found that when three time steps are employed, the time step associated with the costly long-range forces could be pushed to values as large as 100 fs for d-AFED calculations with RESPA1 and UFED calculations with RESPA1 and RESPA2 with negligible change in the Ramachandran free-energy surface. For free-energy values below a cutoff of 7.5 kcal/mol, even d-AFED calculations with 100 fs could be employed. Similar results were found for the solvated alanine tripeptide, for which a four-dimensional surface in the two pairs of Ramachandran angles was studied. For the particular implementation of the SIN(R) in conjunction with RESPA1 or RESPA2 and the d-AFED/TAMD and UFED algorithms, we obtain a speedup of roughly 20-25 when the largest time step is pushed to 100 fs. The gains in both time step and efficiency possible by combining the SIN(R) resonance-free multiple time step algorithm with the enhanced sampling methods employed here are significant, and we expect that combining SIN(R) with other enhanced conformational sampling approaches has the potential to yield similar significant gains.

This work was supported by the National Science Foundation (NSF) No. CHE-1565980.

In this appendix, we provide an analysis of Eqs. (7) and show that in the adiabatic limit, they generate the correct free-energy profile in the coarse-grained variables. In order to keep the analysis simple, we will work with two variables x and y with masses mx and my, respectively, where mxmy, with associated temperatures Tx and Ty, where TxTy. Under these conditions, we seek to show that a double isokinetic scheme will generate the potential of mean force and free energy in x. These two variables are assumed to be coupled by a potential U(x, y).

The analysis follows closely that presented in Ref. 37. The equations are analyzed in the absence of the friction and random force terms; hence, the double isokinetic equations of motion for this xy system (taking L = 1) read

ẋ=vx,v̇x=Fx(x,y)mxλvx,v̇1x=λv1xv2xv1x,v̇2x=G(v1x,Tx),ẏ=vy,v̇y=Fy(x,y)myμvy,v̇1y=μv1yv2yv1y,v̇2y=G(v1y,Ty).
(A1)

Here, Fx = −∂U/∂x, Fy = −∂U/∂y, and we have taken the thermostat masses all to be 1. Then G(u, T) = u2kBT. The parameters λ and μ are Lagrange multipliers for enforcing the constraints

Λx=mxvx2+12v1x2=βx1,Λy=myvy2+12v1y2=βy1,
(A2)

where βx = 1/kBTx and βy = 1/kBTy. The Lagrange multipliers are given by

λ=1ΛxvxFx(x,y)12v1x2v2x,μ=1ΛyvyFy(x,y)12v1y2v2y.
(A3)

Since the equations of motion are inherently non-Hamiltonian, we employ the classical non-Hamiltonian statistical mechanical analysis framework in Refs. 64 and 65. This frame generates an effective “microcanonical” partition function Ω corresponding to the equations of motion that takes the form

Ω=dxg(x)δΛx(vx,v1x)βx1δΛy(vy,v1y)βy1.
(A4)

Here x = {x, y, vx, vy, v1x, v2x, v1y, v2y} is the full phase-space vector. The metric factor g(x) is determined by the phase compressibility factor. When the equations of motion are written compactly as ẋ=ξ(x), where ξ(x) is a vector field, the compressibility factor is given by

κ(x)=xξ(x).
(A5)

Once κ(x) is determined, we seek to express it as a total time derivative κ(x) = dw(x)/dt for some function w(x), in which case it follows that the metric factor is g(x)=exp(w(x)).

The calculation of κ(x) for the isokinetic equations was worked out in Ref. 37. In Eq. (19) of Ref. 37, the compressibility is given by

κ=βxFx(x,y)vx+v2xv1x2v2x+βyFy(x,y)vy+v2yv1y2v2y.
(A6)

Because the two temperatures Tx and Ty are not equal, Eq. (A6) cannot be written as a total time derivative of any function. However, in the limit that x and y are fully decoupled adiabatically, we can assume that the dynamics causes the y subsystem to sample all of its local phase space in a neighborhood where x is essentially fixed and vx = v1x = v2x ≈ 0. In this limit, the compressibility can be approximated as

κβyFy(x,y)vy+v2yv1y2v2y,
(A7)

which can be written as

κβyddtU(x,y)+12ddtv2y2,
(A8)

where the value of x is held fixed. Hence, the metric factor becomes

g(x)eβyU(x,y)eβyv2y2/2
(A9)

and the local partition function generated is

Ω(x)=dvydv1ydv2ydyeβyU(x,y)eβyv2y2/2×δΛy(vy,v1y)βy1eβyA(x),
(A.10)

where the last line follows from the definition of the marginal distribution function and corresponding Helmholtz free energy.

The authors of Ref. 10 referred to the limiting dynamics of x when it is fully adiabatically decoupled from y. In this limit, Eqs. (A1) reduce to

ẋ=vx,v̇x=Fx(x)mxλvx,v̇1x=λv1xv2xv1x,v̇2x=G(v1x,Tx),
(A11)

where Fx(x)=dA(x)/dx is the gradient of the free energy. The Lagrange multiplier λ is then given by

λ=1ΛxvxFx(x)12v1x2v2x.
(A12)

The analysis of this isokinetic subsystem then shows that the distribution it will generate in x is ρadb(x, vx, v1x, v2x) ∝ exp(−βxA(x)) or, more precisely,

ρadb(x,vx,v1x,v2x)=dxeβxA(x)dv2xeβxv2x2/2×dvxdv1xdv2xδΛx(vx,v1x)βx1dxeβyA(x)Ty/TxdxdyeβyU(x,y)Ty/Tx
(A13)

which is the expected adiabatic distribution.8,9,42

The analysis illustrates the importance of the adiabatic decoupling of the two subsystems in obtaining the correct free-energy surface. It does not indicate how strictly, in practice, the adiabatic decoupling must be enforced. The examples presented in this study suggest, however, that the adiabatic decoupling condition need not be particularly stringent in order to obtain reasonably accurate results, particularly for UFED, with very large time steps.

1.
A.
Laio
and
M.
Parrinello
, “
Escaping free energy minima
,”
Proc. Natl. Acad. Sci. U. S. A.
99
,
12562
(
2002
).
2.
A.
Barducci
,
G.
Bussi
, and
M.
Parrinello
, “
Well-tempered metadynamics: A smoothly converging and tunable free-energy method
,”
Phys. Rev. Lett.
100
,
020603
(
2008
).
3.
M.
Bonomi
and
M.
Parrinello
, “
Enhanced sampling in the well-tempered ensemble
,”
Phys. Rev. Lett.
104
,
190601
(
2010
).
4.
D.
Branduardi
,
G.
Bussi
, and
M.
Parrinello
, “
Metadynamics with adaptive Gaussians
,”
J. Chem. Theory Comput.
8
,
2247
2254
(
2012
).
5.
S.
Piana
and
A.
Laio
, “
A bias-exchange approach to protein folding
,”
J. Phys. Chem. B
111
,
4553
(
2007
).
6.
G. A.
Tribello
,
A.
Gareth
,
M.
Ceriotti
, and
M.
Parrinello
, “
A self-learning algorithm for biased molecular dynamics
,”
Proc. Natl. Acad. Sci. U. S. A.
107
,
17509
17514
(
2010
).
7.
T.
Huber
,
A. E.
Torda
, and
W.
van Gunsteren
, “
Local elevation: A method for improving the searching properties of molecular dynamics simulation
,”
J. Comput.-Aided Mol. Des.
8
,
695
708
(
1994
).
8.
L.
Rosso
,
P.
Mináry
,
Z.
Zhu
, and
M. E.
Tuckerman
, “
On the use of the adiabatic molecular dynamics technique in the calculation of free energy profiles
,”
J. Chem. Phys.
116
,
4389
4402
(
2002
).
9.
J. B.
Abrams
and
M. E.
Tuckerman
, “
Efficient and direct generation of multidimensional free energy surfaces via adiabatic dynamics without coordinate transformations
,”
J. Phys. Chem. B
112
,
15742
15757
(
2008
).
10.
L.
Maragliano
and
E.
Vanden-Eijnden
, “
A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations
,”
Chem. Phys. Lett.
426
,
168
175
(
2006
).
11.
M.
Chen
,
M. A.
Cuendet
, and
M. E.
Tuckerman
, “
Heating and flooding: A unified approach for rapid generation of free energy surfaces
,”
J. Chem. Phys.
137
,
024102
(
2012
).
12.
S.
Awasthi
and
N. N.
Nair
, “
Exploring high dimensional free energy landscapes: Temperature accelerated sliced sampling
,”
J. Chem. Phys.
146
,
094108
(
2017
).
13.
M.
Bonomi
,
D.
Branduardi
,
F. L.
Gervasio
, and
M.
Parrinello
, “
The unfolded ensemble and folding mechanism of the C-terminal GB1 beta-hairpin
,”
J. Am. Chem. Soc.
130
,
13938
13944
(
2008
).
14.
F.
Pietrucci
and
A.
Laio
, “
A collective variable for the efficient exploration of protein beta-sheet structures: Application to SH3 and GB1
,”
J. Chem. Theory Comput.
5
,
2197
2201
(
2009
).
15.
P.
Raiteri
,
R.
Martonak
, and
M.
Parrinello
, “
Exploring polymorphism: The case of benzene
,”
Angew. Chem., Int. Ed.
44
,
3769
3773
(
2005
).
16.
L.
Filion
,
M.
Marechal
,
B.
von Oorschot
,
D.
Pelt
,
F.
Smallenburg
, and
M.
Dijkstra
, “
Efficient method for predicting crystal structures at finite temperature: Variable box shape simulations
,”
Phys. Rev. Lett.
103
,
188302
(
2009
).
17.
T. Q.
Yu
and
M. E.
Tuckerman
, “
A temperature accelerated approach for exploring polymorphism in molecular crystals
,”
Phys. Rev. Lett.
107
,
015701
(
2011
).
18.
T. Q.
Yu
,
P. Y.
Chen
,
M.
Chen
,
A.
Samanta
,
E.
Vanden-Eijnden
, and
M. E.
Tuckerman
, “
Order-parameter-aided temperature-accelerated sampling for the exploration of crystal polymorphism and solid-liquid phase transitions
,”
J. Chem. Phys.
140
,
214109
(
2014
).
19.
A. M.
Reilly
,
R. I.
Cooper
,
C. S.
Adjiman
,
S.
Bhattacharya
,
A. D.
Boese
,
J. G.
Brandenburg
,
P. J.
Bygrave
,
R.
Bylsma
,
J. E.
Campbell
,
R.
Car
,
D. H.
Case
,
R.
Chadha
,
J. C.
Cole
,
K.
Cosburn
,
H. M.
Cuppen
,
F.
Curtis
,
G. M.
Day
,
R. A.
DiStasio
,
A.
Dzyabchenko
,
B. P.
van Eijck
,
D. M.
Elking
,
J. A.
van den Ende
,
J. C.
Facelli
,
M. B.
Ferraro
,
L.
Fusti-Molnar
,
C.-A.
Gatsiou
,
T. S.
Gee
,
R.
de Gelder
,
L. M.
Ghiringhelli
,
H.
Goto
,
S.
Grimme
,
R.
Guo
,
D. W. M.
Hofmann
,
J.
Hoja
,
R. K.
Hylton
,
L.
Iuzzolino
,
W.
Jankiewicz
,
D. T.
de Jong
,
J.
Kendrick
,
N. J. J.
de Klerk
,
H.-Y.
Ko
,
L. N.
Kuleshova
,
X.
Li
,
S.
Lohani
,
F. J. J.
Leusen
,
A. M.
Lund
,
J.
Lv
,
Y.
Ma
,
N.
Marom
,
A. E.
Masunov
,
P.
McCabe
,
D. P.
McMahon
,
H.
Meekes
,
M. P.
Metz
,
A. J.
Misquitta
,
S.
Mohamed
,
B.
Monserrat
,
R. J.
Needs
,
M. A.
Neumann
,
J.
Nyman
,
S.
Obata
,
H.
Oberhofer
,
A. R.
Oganov
,
A. M.
Orendt
,
G. I.
Pagola
,
C. C.
Pantelides
,
C. J.
Pickard
,
R.
Podeszwa
,
L. S.
Price
,
S. L.
Price
,
A.
Pulido
,
M. G.
Read
,
K.
Reuter
,
E.
Schneider
,
C.
Schober
,
G. P.
Shields
,
P.
Singh
,
I. J.
Sugden
,
K.
Szalewicz
,
C. R.
Taylor
,
A.
Tkatchenko
,
M. E.
Tuckerman
,
F.
Vacarro
,
M.
Vasileiadis
,
A.
Vazquez-Mayagoitia
,
L.
Vogt
,
Y.
Wang
,
R. E.
Watson
,
G. A.
de Wijs
,
J.
Yang
,
Q.
Zhu
, and
C. R.
Groom
, “
Report on the sixth blind test of organic crystal structure prediction methods
,”
Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater.
72
,
439
459
(
2016
).
20.
E.
Schneider
,
L.
Vogt
, and
M. E.
Tuckerman
, “
Exploring polymorphism of benzene and naphthalene with free energy based enhanced molecular dynamics
,”
Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater.
72
,
542
550
(
2016
).
21.
O.
Valsson
and
M.
Parrinello
, “
Thermodynamical description of a quasi-first-order phase transition from the well-tempered ensemble
,”
J. Chem. Theory Comput.
9
,
5267
5276
(
2013
).
22.
A.
Samanta
,
M. E.
Tuckerman
,
T. Q.
Yu
, and
W.
E
, “
Microscopic mechanisms of the equilibrium melting of a solid
,”
Science
346
,
729
(
2014
).
23.
A. J.
Clark
,
R.
Tiwary
,
K.
Borrelli
,
S. L.
Feng
,
E. B.
Miller
,
R.
Abel
,
R. A.
Friesner
, and
B. J.
Berne
, “
Prediction of protein ligand binding poses via a combination of induced fit docking and metadynamics simulations
,”
J. Chem. Theory Comput.
12
,
2990
2998
(
2016
).
24.
E.
Chiavazzo
,
R.
Covino
,
R. R.
Coifman
,
C. W.
Gear
,
A. S.
Georgiou
,
G.
Hummer
, and
I. G.
Kevrekidis
, “
Intrinsic map dynamics exploration for uncharted effective free-energy landscapes
,”
Proc. Natl. Acad. Sci. U. S. A.
114
,
E494
(
2017
).
25.
M.
Tuckerman
,
G.
Martyna
, and
B.
Berne
, “
Reversible multiple time scale molecular dynamics
,”
J. Chem. Phys.
97
,
1990
(
1992
).
26.
A.
Sandu
and
T.
Schlick
, “
Masking resonance artifacts in force-splitting methods for biomolecular simulations by extrapolative langevin dynamics
,”
J. Comput. Phys.
151
,
74
(
1999
).
27.
E.
Barth
and
T.
Schlick
, “
Extrapolation versus impulse in multiple timestepping schemes. II. Linear analysis and applications to Newtonian and Langevin dynamics
,”
J. Chem. Phys.
109
,
1633
1642
(
1998
).
28.
Q.
Ma
,
J. A.
Izaguirre
, and
R. D.
Skeel
, “
Verlet-I/R-RESPA/impulse is limited by nonlinear instabilities
,”
SIAM J. Sci. Comput.
24
,
1951
1973
(
2003
).
29.
J. A.
Morrone
,
T. E.
Markland
,
M.
Ceriotti
, and
B. J.
Berne
, “
Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances
,”
J. Chem. Phys.
134
,
014103
(
2011
).
30.
D. J.
Evans
and
O. P.
Morriss
, “
Non-Newtonian molecular dynamics
,”
Comput. Phys. Rep.
1
,
297
(
1984
).
31.
D. J.
Evans
and
G. P.
Morriss
,
Statistical Mechanics of Nonequilibrium Liquids
(
Harcourt Brace Javanovich
,
London
,
1980
).
32.
P.
Minary
,
G. J.
Martyna
, and
M. E.
Tuckerman
, “
Algorithms and novel applications based on the isokinetic ensemble. I. Biophysical and path integral molecular dynamics
,”
J. Chem. Phys.
118
,
2510
(
2003
).
33.
M. E.
Tuckerman
,
Statistical Mechanics: Theory and Molecular Simulation
(
Oxford University Press
,
Oxford
,
2010
).
34.
G. J.
Martyna
,
M. L.
Klein
, and
M.
Tuckerman
, “
Nosé–Hoover chains: The canonical ensemble via continuous dynamics
,”
J. Chem. Phys.
97
,
2635
2643
(
1992
).
35.
B.
Leimkuhler
,
E.
Noorizadeh
, and
F.
Theil
, “
A gentle stochastic thermostat for molecular dynamics
,”
J. Stat. Phys.
135
,
261
277
(
2009
).
36.
P.
Minary
,
M. E.
Tuckerman
, and
G. J.
Martyna
, “
Long time molecular dynamics for enhanced conformational sampling in biomolecular systems
,”
Phys. Rev. Lett.
93
,
150201
(
2004
).
37.
B.
Leimkuhler
,
D. T.
Margul
, and
M. E.
Tuckerman
, “
Stochastic, resonance-free multiple time-step algorithm for molecular dynamics with very large time steps
,”
Mol. Phys.
111
,
3579
3594
(
2013
).
38.
D. T.
Margul
and
M. E.
Tuckerman
, “
A stochastic, resonance-free multiple time-step algorithm for polarizable models that permits very large time steps
,”
J. Chem. Theory Comput.
12
,
2170
2180
(
2016
).
39.
A.
Albaugh
,
H. A.
Boateng
,
R. T.
Bradshaw
,
O. N.
Demerdash
,
J.
Dziedzic
,
Y. Z.
Mao
,
D. T.
Margul
,
J.
Swails
,
Q.
Zeng
,
D. A.
Case
,
P.
Eastman
,
L. P.
Wang
,
J. W.
Essex
,
M.
Head-Gordon
,
V. S.
Pande
,
J. W.
Ponder
,
Y. H.
Shao
,
C. K.
Skylaris
,
I. T.
Todorov
,
M. E.
Tuckerman
, and
T.
Head-Gordon
, “
Advanced potential energy surfaces for molecular simulation
,”
J. Phys. Chem. B
120
,
9811
9832
(
2016
).
40.
E. A.
Carter
,
G.
Ciccotti
,
J. T.
Hynes
, and
R.
Kapral
, “
Constrained reaction coordinate dynamics for the simulation of rare events
,”
Chem. Phys. Lett.
156
,
472
(
1989
).
41.
M.
Sprik
and
G.
Ciccotti
, “
Free energy from constrained molecular dynamics
,”
J. Chem. Phys.
109
,
7737
(
1998
).
42.
S.
Samuelson
,
A.
Hughes
, and
G. J.
Martyna
, “
Solvent, force field, temperature and quantum effects on the folding free energy surface of blocked alanine tripeptide
,”
J. Chim. Phys.
94
,
1503
(
1997
).
43.
M. A.
Cuendet
and
M. E.
Tuckerman
, “
Free energy reconstruction from metadynamics or adiabatic free energy dynamics simulations
,”
J. Chem. Theory Comput.
10
,
2975
(
2014
).
44.

The equations of TAMD/d-AFED reflect something like an inversion of the well-known Car-Parrinello (CP)66 algorithm for ab initio molecular dynamics. In the CP approach, fictitious dynamics of the electronic degrees of freedom are kept at a low temperature TeT and are assigned a fictitious mass μ such that they move on a time scale that is fast compared to the nuclei. In this limit, they generate a potential of mean force for the physical nuclear degrees of freedom that approximates the true Born-Oppenheimer electronic ground-state surface. In TAMD/d-AFED, the fictitious degrees of freedom are hot compared to the physical degrees of freedom (TsT) and are made to move slowly via the fictitious mass assignment. In this limit, the physical degrees of freedom generate a potential of mean force for the fictitious degrees of freedom that approximates the true Helmholtz free-energy surface as a function of the fictitious variables.

45.
Y.
Liu
and
M. E.
Tuckerman
, “
Generalized Gaussian moment thermostatting: A new continuous dynamical approach to the canonical ensemble
,”
J. Chem. Phys.
112
,
1685
(
2000
).
46.
G. J.
Martyna
,
M. E.
Tuckerman
,
D. J.
Tobias
, and
M. L.
Klein
, “
Explicit reversible integrators for extended systems dynamics
,”
Mol. Phys.
87
,
1117
1157
(
1996
).
47.
Z.
Zhang
,
X.
Liu
,
Z.
Chen
,
H.
Zheng
,
K.
Yan
, and
J.
Liu
, “
A unified thermostat scheme for efficient configurational sampling for classical/quantum canonical ensembles via molecular dynamics
,”
J. Chem. Phys.
147
,
034109
(
2017
).
48.
K. F.
Gauss
, “
Über ein allgemeines grundgesetz der mechanik
,”
Reine Angew. Math.
IV
,
232
(
1829
).
49.
J. A.
Morrone
,
R. H.
Zhou
, and
B. J.
Berne
, “
Molecular dynamics with multiple time scales: How to avoid pitfalls
,”
J. Chem. Theory Comput.
6
,
1798
1804
(
2010
).
50.
P.
Procacci
,
M.
Marchi
, and
G. J.
Martyna
, “
Electrostatic calculations and multiple time scales in molecular dynamics simulation of flexible molecular systems
,”
J. Chem. Phys.
108
,
8799
8803
(
1998
).
51.
M. E.
Tuckerman
,
D. A.
Yarne
,
S. O.
Samuelson
,
A. L.
Hughes
, and
G. J.
Martyna
, “
Exploiting multiple levels of parallelism in molecular dynamics based calculations via modern techniques and software paradigms on distributed memory computers
,”
Comput. Phys. Commun.
128
,
333
376
(
2000
).
52.
A. D.
MacKerell
,
D.
Bashford
,
M.
Bellott
,
R. L.
Dunbrack
,
J. D.
Evanseck
,
M. J.
Field
,
S.
Fischer
,
J.
Gao
,
H.
Guo
,
S.
Ha
,
D.
Joseph-McCarthy
,
L.
Kuchnir
,
K.
Kuczera
,
F. T.
Lau
,
C.
Mattos
,
S.
Michnick
,
T.
Ngo
,
D. T.
Nguyen
,
B.
Prodhom
,
W. E.
Reiher
,
B.
Roux
,
M.
Schlenkrich
,
J. C.
Smith
,
R.
Stote
,
J.
Straub
,
M.
Watanabe
,
J.
Wiórkiewicz-Kuczera
,
D.
Yin
, and
M.
Karplus
, “
All-atom empirical potential for molecular modeling and dynamics studies of proteins
,”
J. Phys. Chem. B
102
,
3586
3616
(
1998
).
53.
M.
Suzuki
, “
General theory of fractal path-integrals with applications to many-body theories and statistical physics
,”
J. Math. Phys.
32
,
400
(
1991
).
54.
M.
Suzuki
, “
General nonsymmetric higher-order decomposition of exponential operators and symplectic integrators
,”
J. Phys. Soc. Jpn.
61
,
3015
(
1992
).
55.
H.
Yoshida
, “
Construction of higher-order symplectic integrators
,”
Phys. Lett. A
150
,
262
(
1990
).
56.
L.
Rosso
and
M. E.
Tuckerman
, “
An adiabatic molecular dynamics method for the calculation of free energy profiles
,”
Mol. Simulat.
28
,
91
(
2002
).
57.
F.
Paesani
,
W.
Zhang
,
D. A.
Case
,
T. E.
Cheatham
, and
G. A.
Voth
, “
An accurate and simple quantum model for liquid water
,”
J. Chem. Phys.
125
,
184507
(
2006
).
58.
U.
Essmann
,
L.
Perera
,
M. L.
Berkowitz
,
T.
Darden
,
H.
Lee
, and
L. G.
Pedersen
, “
A smooth particle mesh Ewald method
,”
J. Chem. Phys.
103
,
8577
8593
(
1995
).
59.
D. A.
Gibson
and
E. A.
Carter
, “
Time-reversible multiple time scale ab initio molecular dynamics
,”
J. Phys. Chem.
97
,
13429
13434
(
1993
).
60.
S. O.
Samuelson
and
G. J.
Martyna
, “
Two dimensional umbrella sampling techniques for the computer simulation study of helical peptides at thermal equilibrium: The 3K(I) peptide in vacuo and solution
,”
J. Chem. Phys.
109
,
11061
(
1998
).
61.
E.
Darve
,
D.
Rodríguez-Gómez
, and
A.
Pohorille
, “
Adaptive biasing force method for scalar and vector free energy calculations
,”
J. Chem. Phys.
128
,
144120
(
2008
).
62.
L.
Maragliano
and
E.
Vanden-Eijnden
, “
Single-sweep methods for free energy calculations
,”
J. Chem. Phys.
128
,
184110
(
2008
).
63.
J. W.
Ponder
, TINKER – software tools for molecular design, 2004, http://dasher.wustl.edu/ffe/downloads/guide.pdf.
64.
M. E.
Tuckerman
,
C. J.
Mundy
, and
G. J.
Martyna
, “
On the classical statistical mechanics of non-Hamiltonian systems
,”
Europhys. Lett.
45
,
149
(
1999
).
65.
M. E.
Tuckerman
,
Y.
Liu
,
G.
Ciccotti
, and
G. J.
Martyna
, “
Non-Hamiltonian molecular dynamics: Generalizing Hamiltonian phase space principles to non-Hamiltonian systems
,”
J. Chem. Phys.
115
,
1678
(
2001
).
66.
R.
Car
and
M.
Parrinello
, “
Unified approach for molecular dynamics and density-functional theory
,”
Phys. Rev. Lett.
55
,
2471
(
1985
).