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.

## I. INTRODUCTION

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 algorithm^{1} and its variants.^{2–6} Metadynamics is essentially a basin-filling scheme^{7} 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) approach^{11} 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 thermostatting^{29} and novel combinations of the isokinetic ensemble^{30–33} with Nosé-Hoover chain thermostats^{34} 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,^{29} *a 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.

## II. THEORY

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 **r**_{1}, …, **r**_{N} ≡ **r** 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*,

where *β* = 1/*k*_{B}*T* 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**) = −*k*_{B}*T* 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 ensemble^{40,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

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

where the momenta **p**_{1}, …, **p**_{N} ≡ **p** are conjugate to the physical coordinates, *m*_{1}, …, *m*_{N} 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 proved

^{8,9}that the marginal distribution is properly sampled. Introducing the adiabatic decoupling

^{8,9,10,42,66}also allows the sampling to be accelerated via the introduction of a high temperature

*T*

_{s}≫

*T*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

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

In Refs. 11 and 43, however, it was shown that use of free-energy gradient information, obtained by averages of the harmonic force $\kappa \alpha \theta \alpha (r)\u2212s\alpha $, 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

where **F**_{i}(**r**) = −*∂U*/*∂***r**_{i} 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 phenomena^{26–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 algorithm^{37–39} and its deterministic Nosé-Hoover chain antecedent^{36} 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 *q*_{k} denote any one of the 3*N* physical Cartesian degrees of freedom, and let $vk=q\u0307k$ 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 *ν*_{α} = *π*_{α}/*μ*_{α} = $s\u0307\alpha $ 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

Here, *F*_{k}(**r**) is the *k*^{th} force component, *G*(*y*) = *Q*_{1}*y*^{2} − *k*_{B}*T*, $G(y)=Q1y2\u2212kBTs$, *Q*_{1} and *Q*_{2} 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, $\sigma k=2\gamma k/\beta Q2$, and $\Sigma \alpha =2\Gamma \alpha /\beta sQ2$, where, *β*_{s} = 1/*k*_{B}*T*_{s}. The stochastic processes d$wk$ and d*W*_{α} 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:

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

where $F\alpha (r,s\alpha )=\kappa \alpha (\theta \alpha (r)\u2212s\alpha )$. 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 *U*_{bias}(**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 *t*_{l} < *t*, where *t* is the current time of the simulation, i.e.,

where ∥⋅∥ is the usual *L*^{2} 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 −*∂U*_{bias}/*∂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.

## III. BENCHMARK SIMULATIONS

### A. General simulation protocols for solvated peptides

All the simulations were performed using the PINY_MD package^{51} 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}/rad^{2}, and the harmonic coupling constants *κ*_{Φ,Ψ} were both set to 2.78 × 10^{3} (kcal/mol)/rad^{2}. The temperature *T* was set as 300 K for the physical system. The extended system was maintained at temperature *T*_{s} = 1000 K for UFED and at *T*_{s} = 1500 K for d-AFED. In all the simulations, the values of $Q1$ and $Q2$ were determined by the relation $Q1,2$ = *k*_{B}*T*_{s}*τ*^{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 *n*_{sy} = 3, *n*_{res} = 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 thermostat^{34} 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 model^{57} 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 *L*^{1} error of the FESs was calculated in the following way:

where *N*_{b} is the number of bins used to represent the free-energy surface. **s**_{i} ≡ ((*s*_{1,i}, …, *s*_{n,i})) is a point in the extended space. The free energy *A*(**s**_{i}) was compared to the benchmark free energy *A*_{conv}(**s**_{i}) at that point.

### B. Alanine dipeptide (aqueous phase)

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).

### C. Alanine tripeptide (aqueous phase)

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}

## IV. RESULTS

### A. One-dimensional quartic double well coupled to a harmonic oscillator

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

where *D*_{0} = 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 *m*_{y} = 1 and *T*_{y} = 1. In the AFED simulations, the temperature for *x* was *T*_{x} = 10. The mass of *x* was *m*_{x} ≫ *m*_{y}. 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 *n*_{sy} = 3, the RESPA factor *n*_{res} = 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 *m*_{x} ≫ *m*_{y}, the slow and fast forces on *x* and *y* are, respectively,

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

were updated every inner time step *δt*. The probability distribution of *P*(*x*) was compared with the analytical probability distribution *P*_{an}(*x*) using the *L*^{1} error measure

where *N*_{b} was the number of bins used to generate the distribution. Figure 1 shows that the *L*^{1} error decreases with increasing *m*_{x}. 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, *m*_{x} ≥ 300. Since further increase of *m*_{x} has little influence of the *L*^{1} error, we chose *m*_{x} = 300 in the proceeding studies.

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.

The sensitivity of the error to the values of thermostat parameters *Q*_{1}, *Q*_{2}, and *γ* was also studied (Fig. 3). For a wide range of values of these parameters, only when the values of *γ* and *Q*_{2} are either too small or too large is the error seen to become large. For most values of these parameters, the *L*^{1} error is small and roughly constant.

### B. Alanine dipeptide

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 *L*^{1} error in the two-dimensional free-energy surface. In particular, note that this *L*^{1} 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 *L*^{1} 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 *L*^{1} 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.

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.

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 *t*_{bench} per time step *δt* = 0.5 fs of the double isokinetic d-AFED simulation without RESPA. The time ratio *r* = *t*_{bench}/*t*_{res}, where *t*_{res} 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.

### C. Alanine tripeptide

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 *L*^{1} 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.

Figure 9 shows the *L*^{1} 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.

## V. CONCLUSION

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 approach^{2,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.

## ACKNOWLEDGMENTS

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

### APPENDIX: ANALYSIS OF THE DOUBLE ISOKINETIC METHOD

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 *m*_{x} and *m*_{y}, respectively, where *m*_{x} ≫ *m*_{y}, with associated temperatures *T*_{x} and *T*_{y}, where *T*_{x} ≫ *T*_{y}. 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

Here, *F*_{x} = −*∂U*/*∂x*, *F*_{y} = −*∂U*/*∂y*, and we have taken the thermostat masses all to be 1. Then *G*(*u*, *T*) = *u*^{2} − *k*_{B}*T*. The parameters *λ* and *μ* are Lagrange multipliers for enforcing the constraints

where *β*_{x} = 1/*k*_{B}*T*_{x} and *β*_{y} = 1/*k*_{B}*T*_{y}. The Lagrange multipliers are given by

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

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\u0307=\xi (x)$, where *ξ*(x) is a vector field, the compressibility factor is given by

Once *κ*(x) is determined, we seek to express it as a total time derivative *κ*(x) = d$w$(x)/d*t* for some function $w$(x), in which case it follows that the metric factor is $g(x)=exp(\u2212w(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

Because the two temperatures *T*_{x} and *T*_{y} 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

which can be written as

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

and the local partition function generated is

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

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

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

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.

## REFERENCES

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 *T*_{e} ≪ *T* 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 (*T*_{s} ≫ *T*) 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.