A Stochastic, Resonance-Free Multiple Time-Step Algorithm for

Article pubs.acs.org/JCTC

A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps Daniel T. Margul† and Mark E. Tuckerman*,†,‡,§ †

Department of Chemistry, New York University, New York, New York 10003, United States Courant Institute of Mathematical Sciences, New York University, New York, New York 10003, United States § NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China ‡

ABSTRACT: Molecular dynamics remains one of the most widely used computational tools in the theoretical molecular sciences to sample an equilibrium ensemble distribution and/or to study the dynamical properties of a system. The efficiency of a molecular dynamics calculation is limited by the size of the time step that can be employed, which is dictated by the highest frequencies in the system. However, many properties of interest are connected to low-frequency, long time-scale phenomena, requiring many small time steps to capture. This ubiquitous problem can be ameliorated by employing multiple time-step algorithms, which assign different time steps to forces acting on different time scales. In such a scheme, fast forces are evaluated more frequently than slow forces, and as the former are often computationally much cheaper to evaluate, the savings can be significant. Standard multiple time-step approaches are limited, however, by resonance phenomena, wherein motion on the fastest time scales limits the step sizes that can be chosen for the slower time scales. In atomistic models of biomolecular systems, for example, the largest time step is typically limited to around 5 fs. Previously, we introduced an isokinetic extended phase-space algorithm (Minary et al. Phys. Rev. Lett. 2004, 93, 150201) and its stochastic analog (Leimkuhler et al. Mol. Phys. 2013, 111, 3579) that eliminate resonance phenomena through a set of kinetic energy constraints. In simulations of a fixed-charge flexible model of liquid water, for example, the time step that could be assigned to the slow forces approached 100 fs. In this paper, we develop a stochastic isokinetic algorithm for multiple time-step molecular dynamics calculations using a polarizable model based on fluctuating dipoles. The scheme developed here employs two sets of induced dipole moments, specifically, those associated with short-range interactions and those associated with a full set of interactions. The scheme is demonstrated on the polarizable AMOEBA water model. As was seen with fixed-charge models, we are able to obtain large time steps exceeding 100 fs, allowing calculations to be performed 10 to 20 times faster than standard thermostated molecular dynamics.

1. INTRODUCTION Sampling the conformational equilibria of complex systems constitutes one of the primary challenges in molecular simulation. This problem directly influences the determination of molecular crystal structures, the elucidation of the conformational preferences of biomolecules, the study of equilibria in glassy systems, and the exploration of mechanisms of first-order phase transitions, to highlight just a few examples, most of which are connected to an underlying rough potential energy landscape. Numerous approaches have been developed to accelerate the sampling of the corresponding equilibrium ensemble distribution using molecular dynamics (MD) as the engine for driving the exploration of configuration space.1−12 As powerful as these approaches are, they are often designed as single time-step methods, and therefore, they are limited in their efficiency by the size of the time step, which must be chosen small enough to integrate the highest frequency motion with sufficient accuracy. However, the complex systems these methods were developed to study are characterized by forces that give rise to motion on many time scales, a problem that is best treated using multiple time© 2016 American Chemical Society

step (MTS) algorithms. In general, the efficiency of both conventional and accelerated MD approaches can be improved by combining them with MTS algorithms. MTS methods assign different time steps to different force components within a symplectic integration scheme so that the slower, more computationally expensive forces can be evaluated less frequently than the faster, cheaper forces. Although MTS methods are typically applied in force-field-based calculations, they have also been adapted for ab initio molecular dynamics, in which the interparticle forces are determined from electronic structure calculations performed “on the fly” as the MD calculation proceeds, with notable gains in efficiency.13−17 Early MTS algorithms,18−20 particularly, the reversible Reference System Propagator Algorithm (r-RESPA),21,22 were heralded as a significant breakthrough in MD methodology. However, it was subsequently discovered that a fundamental limitation existed on the largest possible time step that could be Received: February 20, 2016 Published: April 7, 2016 2170

DOI: 10.1021/acs.jctc.6b00188 J. Chem. Theory Comput. 2016, 12, 2170−2180

Article

Journal of Chemical Theory and Computation

coupled to L auxiliary thermostat velocities v(k) 1,α, k = 1,...,L, through an isokinetic constraint enforced via a Lagrange multiplier λα, and each v(k) 1,α is coupled to a stochastic variable 39 v(k) The stochastic 2,α via a stochastic Nosé−Hoover coupling. equations of motion that result are then

achieved in these approaches. Specifically, resonance phenomena23−25 numerically couple the motion on different time scales and limit the magnitude of the time step that can be chosen for the slowest force components. Several groups have made important strides in the development of resonance-free MTS schemes26−33 and have obtained significant gains in the maximum possible time step. The method we have developed is based on the isokinetic ensemble28,33−37 in an extended phasespace formulation.28 In particular, the physical velocities in a system are coupled to a set of Nosé−Hoover chain thermostats38 or stochastic Nosé−Hoover thermostats39 via an isokinetic (constant kinetic energy) constraint. When the kinetic energy is bound through the imposition of such a constraint, energy is prevented from building up in any one particular mode, which is, therefore, prevented from becoming resonant.28,33 In MD simulations of a flexible fixed-charge water model, time steps for the slowest forces in the system could be increased to approximately 100 fs with no degradation in the structural properties generated. It should be noted, however, that this scheme is designed to sample the equilibrium isokinetic distribution and is not guaranteed to reproduce dynamics properly. Recent approaches employing colored noise thermostats32 are able to preserve dynamical properties, however, the gains in time step, while impressive, are more modest. The goal of this paper is to extend the stochastic resonancefree approach of ref 33 to polarizable models formulated in terms of fluctuating induced dipole moments. Typical MTS algorithms for fixed-charge force fields involve decomposing the forces into bonded, short-range nonbonded, and long-range nonbonded components. Although the forces in polarizable models cannot be explicitly decomposed in this fashion, the concept can, nevertheless, be adapted to the polarizable case by introducing a short-ranged induced dipole moment calculation based on a spherical cutoff, which is then used to construct induced dipole forces in the short-range steps. This is then corrected in longrange force steps by performing a full calculation of the induced dipole moments and subtracting the short-range contribution. We have implemented our approach in the TINKER7 software package,40 and in MD simulations of the AMOEBA water model,41 we find that we are able to increase the time step for the expensive long-range calculations to 120 fs with no degradation in the equilibrium properties. Importantly, we obtain computational speedup factors ranging between 10 and 25, depending on the choice of simulation parameters. This paper is organized as follows: In Section II, we review the basic equations of motion and integration strategies. We then review force and energy decomposition schemes for fixed-charge models and proceed to describe how the approach is implemented for the AMOEBA polarizable model. In Section III, we show results for molecular dynamics simulations of AMOEBA water, exploring how large the time step can be made. We also provide timings and measures of efficiency. Conclusions are given in Section IV.

drα = vαdt ⎡ F (r) ⎤ − λαvα ⎥dt dvα = ⎢ α ⎣ mα ⎦ dv1,(kα) = −λαv1,(kα)dt − v2,(kα)v1,(kα)dt dv2,(kα) =

G(v1,(kα))

dt − γv2,(kα)dt + σdwα(k)

Q2

(1)

In these equations, 2

G(v1,(kα)) = Q 1v1,(kα) − kBT

(2)

where dwα is a stochastic Ornstein−Uhlenbeck (OU) process, γ is the friction constant, σ = 2kBTγ /2Q 2 , and kB is Boltzmann’s constant. Here, Q1 and Q2 are mass-like parameters (having units of energy·time2) for the thermostat variables that determine the characteristic time scale of their evolution. The presence of the stochastic Nosé−Hoover coupling in the extended phase space leads to the designation of this approach as the Stochastic Isokinetic Nosé−Hoover scheme. The dN Lagrange multipliers {λα} are used to enforce dN isokinetic constraints of the form mα vα2 +

L L+1

L

∑ Q 1(v1,(kα))2 = LkBT

(3)

k=1

The Lagrange multipliers are obtained by differentiating eq 3 with respect to time, using the equations of motion (eq 1) to substitute for v̇α and v̇(k) 1,α and then solving the resulting equations for the multipliers. This procedure yields the expression λα =

vαFα(r) − mvα2 +

L L ∑k = 1 Q 1(v1,(kα))2 v2,(kα) L+1 L L ∑k = 1 Q 1(v1,(kα))2 L+1

(4)

When eq 4 is substituted back into the equations of motion, the result is a set of non-Hamiltonian equations, which can be analyzed using the methods of ref 42. As was shown in ref 33, the analysis is done in the absence of the stochastic terms, and the non-Hamiltonian equations are shown to generate the correct isokinetic ensemble partition function Ω(N , V , β) =

∫ ddNrddNvddNLv1ddNLv2 ρisok

(5)

where β = 1/kBT, and dN

L

(k) 2

ρisok = e−β /2 ∑α=1 ∑k =1 Q 2(v2,α) e−βU (r) dN

⎛

∏ δ⎜⎜mαvα2 +

2. THEORY 2.1. Equations of Motion. Consider a system of N particles in d dimensions. We define coordinates r1,...,rdN ≡ r, momenta p1,...,pdN ≡ p, and masses m1,...,mdN associated with each degree of freedom. The system is described by a Hamiltonian H(r,p) = K(p) + U(r) with kinetic energy K(p) and potential energy U(r), which generates forces F1(r),...,FdN(r), Fα(r) = −∂U(r)/∂rα, and α = 1,...,dN. For each coordinate rα, there is an associated velocity vα = ṙα. In our stochastic MTS scheme, each physical velocity vα is

α=1

⎝

L L+1

L

⎞

∑ Q 1(v1,(kα))2 − LkBT ⎟⎟ k=1

⎠

(6)

which is canonical in the coordinates r1,...,rdN. Further analysis33 shows that the isokinetic density ρisok is preserved when the OU process is included. Additionally, using the Hö rmander conditions,43 the equations of motion were shown to be ergodic. Note that the stochastic isokinetic approach can be used as an ordinary thermostat to generate a canonical configurational distribution. Finally, although formulated as a “massive” 2171


Article


where i,j index the N position vectors r1,...,rN, nL is the lattice translation vector, n is a vector of integers, and rnij = |ri − rj + nL|. The δ0ij Kronecker δ is defined to be 1 for i = j only when n = (0,0,0). Using the error function, erf(x), and its complement, erfc(x), as smooth switches and the identity erf(x) + erfc(x) = 1, eq 8 is often evaluated using Ewald summation, wherein it is written as the sum of a screened Coulomb potential in real space and a reciprocal-space term, i.e., Uelec(r) = Ureal(r) + Urecip(r), where

thermostat here, with each degree of freedom coupled its own thermostat, larger subgroups of atoms or degrees of freedom can be coupled to a single thermostat. If the forces can be decomposed into fast and slow slow components, i.e., Fα(r) = Ffast α (r) + Fα (r), then a two timestep MTS integrator can be devised for eq 1, as discussed in ref 33, that preserves the phase-space volume element,42,44 is free of resonances, and allows a large time step to be employed for the integration of the slow component of the forces. The resulting algorithm is termed the Stochastic Isokinetic Nosé-Hoover (RESPA) or SIN(R) scheme. If the forces can be decomposed further into fast, intermediate, and slow components, i.e., Fα(r) = im slow Ffast α (r) + Fα (r) + Fα (r), and then a three time-step SIN(R) scheme results, and so forth. In the Appendix, we provide a few details with regard to the three time-step MTS scheme. 2.2. AMOEBA Polarizable Force Field. We discuss our MTS algorithm in the context of the AMOEBA force field,41,45 which is available in the TINKER molecular modeling software suite. For a system with coordinates r, the potential energy is expressed as

N

Ureal(r) =

N

Urecip(r) = =

N

N

∑ ∑ ∑ (1 − δij0) n

i=1 j≥i

(7)

N

∑ ∑ ∑ (1 − δij0)qiqj n

i=1 j≥i

1 L3

∑

Uewc(r) = −

k ≠ (0,0,0)

(9)

erf(αrijn) rijn

2π −k2 /4α 2 e |S(k)|2 + Uself k2

(10)

1 2

∑

qiqjerf(α|ri − rj|) |ri − rj|

i≠j∈)

(11)

where ) is the set of all “bonded” atoms needed to compute the first five terms in eq 7, interactions to correct for Coulomb interactions that arise when the structure factor S(k) is squared to give |S(k)|2 in eq 10. As was shown by Procacci et al.,48 this correction can also be added in reciprocal space. In practice, the real-space interactions are truncated beyond some real-space cutoff rreal, and the reciprocal-space sum is truncated beyond a value k ≡ |k| > kmax. In this case, Ureal acquires a dependence on rreal, and Urecip acquires a dependence on kmax, i.e., ̃ (r; rreal) = ∑ ∑N ∑N (1 − δij0)q q Ureal i=1 j≥i n i j

erfc(αrijn) rijn

× [1 − Θ(rijn − rreal)] ̃ (r; k max ) = 1 Urecip L3

∑ k ≠ (0,0,0)

(12)

2π −k2 /4α 2 e |S(k)|2 + Uself k2

|k|≤ k max

(13)

where Θ(·) is the Heaviside step function. The term 1−Θ(rnij − rreal) ensures that only terms for which the distance rnij ≤ rreal are included. Note that the step-function truncation can also be replaced by a smooth switch. Moreover, the sum over n is only needed for small cells, where rreal = L/2. In addition, the reciprocal-space contribution is often evaluated using an approach such as the smooth particle−mesh Ewald algorithm47 in order to increase efficiency. A straightforward approach to define intermediate and slow electrostatic contributions to the forces, termed RESPA1, is

qiqj rijn

rijn

where α is a real-space convergence parameter, k = 2πm/L, m is a α vector of integers, S(k) = ∑iqieik·ri, Uself = − π ∑i (qi2), and electrical neutrality is assumed. All of the above formulas are easily generalizable to the case of cells having shapes other than simple cubic. The term Uself removes the unphysical “selfinteraction” of charges, which results when the reciprocal-space sum is written as it appears in eq 10. Note that the use of an Ewald sum thus formulated requires the addition of an Ewald correction term Uewc(r),47

The first five terms are, respectively, the bond, angle, Urey− Bradley, improper torsion, and normal torsion terms, which comprise the bonded or intramolecular energies. The remaining terms are, respectively, the van der Waals, permanent electrostatic and induced dipolar electrostatic energies, which comprise the nonbonded terms. In most MTS implementations, the fastest forces Ffast α (r) are those arising from the intramolecular terms, i.e., intra Ffast α (r) = Fα (r), and these are associated with the smallest time step, typically around 0.5 fs. In fact, it is sometimes useful to subdivide these further so that normal and improper torsions, which are somewhat more costly to evaluate than bond, bend, and UB terms, are integrated with a larger time step. The definition of the intermediate and slow forces involves the nonbonded contributions to the potential and is the focus of the remaining part of our discussion. In AMOEBA, the van der Waals interaction is a buffered 14-7 Lennard-Jones potential, and the separation of these van der Waals forces into short- and longrange components, which contribute to the intermediate and slow forces, respectively, is described in ref 33. Briefly, this requires expressing the van der Waals forces Fα,vdw(r) as a sum res long res res Fshort α,vdW(r;rvdw) + Fα,vdW(r;rvdw), where rvdw is a short-range cutoff dividing the short-range from the long-range contributions. In the rest of this section, we will first review force decompositions for fixed-charge models and then discuss the technical details of decomposing the permanent electrostatic and induced polarind ization energies, Uperm elec (r) and Uelec(r), respectively, in the AMOEBA model into short- and long-range contributions for use with the SIN(R) MTS algorithm. 2.3. Multiple Time Stepping. 2.3.1. Review of FixedCharge Models. Let us first review the problem of dividing the electrostatic potential of a fixed-charge model into short- and long-range contributions. This can be done in several ways, as has been described in the literature.33,46 For a system in a cubic, periodic box of side length L, the electrostatic potential is Uelec(r) =

erfc(αrijn)

i=1 j≥i

n

U (r) = Ubond(r) + Uangle(r) + UUB(r) + Uimp(r) perm ind + Utors(r) + UvdW(r) + Uelec (r) + Uelec (r)

N

∑ ∑ ∑ (1 − δij0)qiqj

(8) 2172


Article


formulation, termed RESPA2 in refs 46 and 33, formulates the slow contribution to the electrostatic potential as

based on a division of the real-space term in eq 12 into short- and long-range components. That is, we define corresponding slow im potentials Uim elec(r) and Uelec (r), such that Uelec(r) = Uelec(r) + slow Uelec (r), where

slow ̃ long(r) + Urecip ̃ (r; k max ) Uelec (r) = Ureal

im ̃ short(r) Uelec (r) = Ureal slow Uelec (r)

=

−

̃ long(r; rreal) Ureal

̃ (r; k max ) + Urecip

i

(14)

̃ short(r) = Ureal ̃ long(r; rreal) Ureal

erfc(αrij0) 0 res qiqj S(rij ; relec , rij0 i=1 j>i N

=

N

∑ ∑ ∑ (1 − n

δij0)qiqj

λ)

̃ short(r) = Ureal

erfc(αrijn)

i=1 j≥i

The function S(·) is a smooth, differentiable switch ensuring that real-space interactions are only included in the short-range potential if two particles are within a short-range cutoff distance rres elec ≤ rreal, and λ is a healing length that determines the range of the switch. An example of such a switch is

(16)

where (17)

N

∑∑

qiqj rij0

res S(rij0 − relec )

(20)

μ p μj = αj Etot j = αj(Ej + Ej )

and u=

res S(rij0 − relec )

In the RESPA2 decomposition, the short-range component of the reciprocal-space energy is effectively canceled in real-space. However, this cancellation between real and reciprocal spaces can be imperfect, and therefore, a scheme such as that suggested by Procacci et al.48 could also be adapted for use within the RESPA2 algorithm. 2.3.2. The AMOEBA Model. In fixed-charge models, implementation of an MTS scheme is straightforward because the electrostatic potential can be decomposed into a sum of long short- and long-range components, i.e., (Uelec = Ushort elec + Uelec ). A short- and long-range division of the potential cannot be performed directly in the AMOEBA polarization scheme due to the need to solve a set of self-consistent equations for the induced dipole moments. Therefore, a different approach is required. The self-consistent problem arises from the fact that the complete set of dipole moments is determined from the local electric field, i.e.,

(15)

g (r ; rc , λ) = 1 + u3(15u − 6u 2 − 10)

N

i=1 j>i

rijn

res , λ)][1 − Θ(rijn − rreal)] × [1 − S(rij0 ; relec

⎧ 1 r ≤ rc − λ ⎪ S(r ; rc , λ) = ⎨ g (r ; rc , λ) rc − λ ≤ r ≤ rc ⎪ 0 r ≥ rc , ⎩

j>i

rij0

and corrects for this by adding the last term in this expression back into the short-range potential, using the fact that erf(αr0ij) + erfc(αr0ij) = 1, which gives a short-range potential of the form

N

∑∑

erf(αrij0)

(19)

and N

∑ ∑ (1 −

δij0)qiqj

1 (r − rc + λ) λ

Nμ

= αj(Epj +

(18)

∑

Tjk ·μk ) (21)

k = 1, k ≠ j

Note that this switch can also be applied directly to the forces. The remaining real-space interactions, as well as all reciprocalspace interactions, are included in the long-range potential. It is assumed that rres elec is sufficiently small compared to the size of the box that only interactions in the primary cell need to be included in the definition of the short-range real-space potential. Depending on implementation details, the balance between the number of real- and reciprocal-space terms can be adjusted by the choice of α and kmax. An extreme example, for which there is precedence, is to choose α sufficiently large that Uslow elec (r) = Ũ recip(r;kmax) with no residual real-space contribution.49 The reciprocal-space term in the long-range potential clearly includes rapidly varying short-range terms from the larger kvectors. Fortunately, these terms tend to be small due to the Gaussian damping in eq 13. Nevertheless, this is not a perfect decomposition scheme, and an alternative formulation of RESPA1, described in ref 33, replaces Ũ recip(r;kmax) in the longrange term with Ũ recip(r;kres), where kres < kmax is a long-range reciprocal-space cutoff. The incomplete reciprocal-space calculation in the long-range potential is corrected for by adding Ũ recip(r;kmax) − Ũ recip(r;kres) in the short-range potential. However, this scheme requires specialized reciprocal-space sums to achieve efficiency, something very few implementations and packages incorporate.50 Thus, a second alternative

In this equation, the total dipole at the jth polarizable site μj depends on its polarizability αj and two contributions to the total electric field experienced at that site that comprise Etot j . The “static” electric field Epj results from permanent monopole, dipole, and quadrupole electrostatics at site j from all other sites k ≠ j.41 The “induced” electric field depends on the set of mutual induced dipoles and the dipole tensor41 Tjk = λ5,jk

3rjrk r jk5

− λ3,jk

1 r jk3

(22)

where 1 is the 3 × 3 identity matrix, and λ3,jk and λ5,jk are defined as λ3,jk = 1 − exp( −αvjk3 ) λ5,jk = 1 − (1 + αvjk3 )exp( −αvjk3 )

(23)

In the above definitions, vjk = rjk/(αjαk)1/6, where αj and αk are the atomic polarizabilities of j and k, and α with no subscript is the Ewald convergence parameter. Therefore, the induced dipolar energy term Uind elec(r) becomes 2173


Article

Journal of Chemical Theory and Computation ind Uelec (r) = U0 −

∑ μj ·Etot j

rec Uelec (r) =

j

=

∑ j

μj2 2αj

Nμ

∑ μj ·Epj −

−

j=1

1 2

Nμ

μj ·Tjk ·μk

self Uelec =−

k=1

∑

(μj ·∇j Tjk ·μk )

k = 1, k ≠ j

(25)

In previous work, Masella proposed that one could exploit the rapid decay of the polarization forces with increasing interatomic distances within a multiple time-step scheme for polarizable models.51 Due to the nonadditivity of these forces, a simple short-range and long-range division of the forces is not possible. The scheme suggested by Masella introduces a set of short-range induced dipoles that satisfy the self-consistent equation

k = 1, rjk ≤ rc

2

(μ ) =∑ S j

Nμ

2αj

j=1

(27)

With this in mind, the electrostatic potential energy for a periodic system with charges q1,...,qN, dipole vectors μ1,...,μN, and quadrupole tensors Q1,...,QN takes the form perm ind Uelec(r) = Uelec (r) + Uelec (r)

(29)

52 rec self The terms Udir elec(r), Uelec(r), and Uelec are defined as

dir (r) = Uelec

⎧ ⎪

∑ ∑ ∑ ⎨(qi + μi ·∇i + Q i: Ωi) n

i

j≥i

⎪

⎩

⎛ erfc(αrijn) ⎞⎫ ⎪ ⎟⎟⎬ × (qj + μj ·∇j + Q j: Ωj)⎜⎜ n ⎪ rij ⎝ ⎠⎭ × [1 − Θ(rijn − rreal)]

1 2

⎧ ⎪ ⎨(qi + μi ·∇i + Q i: Ωi) ⎪ i≠j∈) ⎩

∑

erf(α|ri − rj|) ⎫ ⎪ ⎬ |ri − rj| ⎪ ⎭

(33)

(34)

(35)

(28)

ind dir rec self = Uelec (r) + Uelec (r) + Uelec (r) + Uelec

(32)

In eq 35, the short-range potential is subtracted from the full potential to generate the long-range contribution. The shortrange potential calculation during a long-range step is necessary, although if the outer time step is long enough, the additional cost of this step will be negligible. As was done in the fixed-charge case, a smooth switch such as that of eq 16 should be used either on the potential or directly on the forces in order to avoid a harsh transition between the short- and long-range regions. 2.3.3. TINKER Implementation. The scheme described in this section has been implemented in the TINKER package within a three time-step framework (see Appendix): A time step dt is assigned to fast forces, which include the bonds, bends, torsions, improper torsions, and UB terms. A time step δt > dt is assigned to the intermediate forces, which include all short-range electrostatic and van der Waals interactions and the Ewald exclusion correction. Finally, a time step Δt > δt > dt is assigned to all long-range interactions. The short-range electrostatic force and energy calculation in TINKER proceeds as follows: First, the real-space portion of the field due to the permanent multipoles is calculated for multipole pairs that are within a distance of rres elec, and the dipole tensor T is constructed. A polynomial predictor is then employed for a first guess of the induced dipoles, which initiates a self-consistent equation solver54,55 for the induced dipoles. The real space portion of the Ewald energy due to induced dipoles and static multipoles for multipole pairs within the distance rres elec is then calculated, and the switching function in eq 16 is applied, in the present study, to the forces. Because the reciprocal-space and long-range real-space interactions are excluded from the intermediate level forces, we expect considerable computational savings if a large outer time step Δt can be used. At the outer step,

In this scheme, short-range induced dipoles are computed more frequently using the aforementioned scheme; when long-range forces are needed, they are obtained by performing a full induced dipole moment calculation and obtaining the long-range forces by subtraction, i.e., ind ind Find L , j = F j − FS , j

⎝

slow ind ind dir rec self Uelec (r) = Uelec (r) − Ushort (r) + Ulong (r) + Uelec (r) + Uelec

∑ μjS ·ESj ,tot j=1

i

⎞ 2α 2 2 4α 4 |μi | + ∥Q i∥2F ⎟ 3 5 ⎠

im ind dir Uelec (r) = Ushort (r) + Ushort (r) + Uewc(r)

Nμ

−

⎫ ⎬ ⎭

We can now define intermediate and slow electrostatic potential energy contributions for use in an MTS scheme, which here, we will refer to as the “RESPA1-pol” algorithm:

(26)

where the static and induced electric fields experienced at polarizable site j are generated from dipoles within a cutoff distance rc of site j. Additionally, the short-range static field ES,p j is generated for multipole contributions within the cutoff rc. The solution to eq 26 has a short-range polarization energy given by ind Uelec, S(r)

⎛

∑ ⎜qi2 +

×(qj + μj ·∇j + Q j: Ωj)

Tjk ·μkS )

∑

α π

Uewc(r) = −

Nμ

μjS = αj(ESj , p +

2

In these equations, Ωi is a second-order differential operator Ωi = ∇i∇i = ∂2/∂ri ∂ri, A:B is the Frobenius inner product of matrices A and B, and ∥·∥F is the Frobenius norm. As in the fixed-charge case, an Ewald exclusion energy Uewc(r) must be added to the bonded potential that corrects for the permanent monopole, dipole, and quadrupole interactions among bonded terms,53 where

Nμ

∑ μk ·∇j (Ekp) +

2

e−k /4α k2

(31)

where U0 = ∑j μ2j /(2αj) is the energy needed to create Nμ dipoles. Since the induced dipoles satisfy ∂Uind elec/∂μj = 0, the forces due to the induced dipole potential become Nμ

max

j = 1 k = 1, k ≠ j

(24)

Find j =

{

∑|kk≠|≤(0,0,0) ∑i ∑j (qi + μi ·∇i + Q i: Ωi) k

× (qj + μj ·∇j + Q j: Ωj)e 2πi k·(rj − ri)

Nμ

∑ ∑

2π L3

(30) 2174


Article

Journal of Chemical Theory and Computation Table 1. Computational Details of SIN(R) Simulations and NVE and NVT Benchmark Simulations name

time step

short range vdW cutoff

short range vdW healing length

short range electrostatic cutoff

short range electrostatic healing length

XO-SIN(R) XM-SIN(R) XI-SIN(R) NVE NVT

75 fs 120 fs 120 fs 0.5 fs 0.5 fs

7.0 Å 7.0 Å 7.0 Å N/A N/A

0.5 Å 0.5 Å 0.5 Å N/A N/A

5.0 Å 5.0 Å 5.0 Å N/A N/A

0.5 Å 0.5 Å 0.5 Å N/A N/A

Figure 1. Two measures of computational efficiency: On the left, the speedup factor η(rreal) as a function of the real space cutoff rreal is shown for SIN(R) simulations; η(rreal) is calculated as the CPU time required for a 600 fs NVT simulation divided by that required for SIN(R). On the right, the slowdown factor τ0/τ is shown for NVT and SIN(R) simulations with rreal = 9.0 Å. For each simulation, τ0 is the average CPU time required for 100 time steps with dipoles converged to within 0.01 D; τ is the average CPU time required for an equal number of time steps for a given dipole convergence criterion.

Figure 2. Oxygen−oxygen (left), oxygen−hydrogen (middle), and hydrogen−hydrogen (right) radial distribution functions gOO(r), gOH(r), and gHH(r) generated by 120 ps NVT and SIN(R) simulations.

For all SIN(R) simulations, the bonded interactions are integrated with a time step of dt = 0.5 fs, and the short-range nonbonded interactions are integrated with a time step δt = 3.0 fs. Full-range nonbonded real-space interactions utilize a cutoff of rreal = 10.4 Å. In the XO-SIN(R) simulations, each velocity degree of freedom is coupled to L = 4 extended phase-space stochastic Nosé−Hoover variable sets, while for XM-SIN(R) and XISIN(R) L = 1. For the induced dipole computation, the dipoles are considered converged when the root-mean-square change of the dipoles is 10−2 Debye or when 500 iterations have been performed. Finally, the Ewald convergence parameter α = 0.3617 Å−1 for all simulations using the 10.4 Å cutoff. For those simulations in which rreal is varied, α is adjusted accordingly. 3.1. Computational Efficiency. A primary motivation for the development of the SIN(R) method is its high computational efficiency relative to existing molecular dynamics (MD) integrators. In the absence of multiple time stepping, MD with polarizable electrostatics can be made faster only by varying the balance between real and reciprocal spaces. Hence, it is important to study the relative performance of NVT (or NVE) and the various SIN(R) algorithms across a range of values of rreal. In Figure 1, it is shown that for rreal ≥ 6.0 Å and a dipole convergence criterion of 10−2 Debye, every SIN(R) variant is between 10 and 25 times more efficient than the NVT benchmark. It is also shown that for an increasingly stringent dipole convergence condition, NVT simulations slow down at a greater rate than the SIN(R) simulations. Thus, a more strict

where the long-range interactions are needed, the short-range interactions are obtained and stored. Then, a full force calculation is performed in both real and reciprocal space, and the short-range forces are subtracted from the full forces to generate the long-range components.

3. RESULTS In this section, the SIN(R) integrator with the RESPA1-pol multipole electrostatics scheme is tested on a periodic cubic system of 512 H2O molecules interacting via the AMOEBA 2003 force field41 in a box of length length of 24.83 Å at a temperature of 298 K. In Table 1, computational details are provided for the SIN(R) simulations as well as for the NVE and NVT benchmarks. The table shows three versions of SIN(R): In the XO-SIN(R) algorithm, the extended phase space variables are updated at the largest time step Δt. In XI-SIN(R), it is updated at the shortest time step dt. In XM-SIN(R), a previously unreported method, it is updated at the intermediate time step δt. It is important to note that the velocity Verlet NVE algorithm in TINKER currently has no multiple time-step option. The NVT algorithm is a Nosé− Hoover chain scheme38 with global thermostatting and chain length of 4 with the option to use a two time-step r-RESPA decomposition.21,56 In the benchmark NVT runs, an inner time step of 0.25 fs is used for bonded interactions, and an outer time step is used for nonbonded interactions. In Table 1, all time steps given are the outermost step length for each method. 2175


Article


Figure 3. L1 errors ξOO(t), ξOH(t), and ξHH(t) of oxygen−oxygen, oxygen−hydrogen, and hydrogen−hydrogen radial distribution functions, respectively, generated by NVT and SIN(R) simulations.

Figure 4. Distributions of hydrogen bonding angles α (left), β (middle), and θ (right) generated by 480 ps NVT and SIN(R) simulations.

Figure 5. L1 error ξα(t), ξβ(t), and ξθ(t) of distributions of hydrogen bond angles α, β, θ, respectively, generated by NVT and SIN(R) simulations.

Figure 6. Distributions (left) of the molecular dipole moment of water generated by NVT and SIN(R) integrators and the average value ⟨P(t)⟩ of those distributions over time (right).

dipole convergence is possible within the SIN(R) scheme with only minimal reduction in computational efficiency. 3.2. Radial Distribution Functions. In Figure 2, the radial distribution functions (RDFs) of AMOEBA water are shown for NVT and SIN(R) simulations. It can be seen that there is very good agreement between the NVT benchmark and SIN(R) results for all of the SIN(R) variants, showing that simulations with outer time steps as large as 120 fs can be performed without loss of accuracy in these RDFs. In Figure 3, the L1 norms of the differences between the distributions at time t, and their longtime, converged values are shown. The L1 norm ξ of a distribution P(xi) from a final distribution Pfinal(xi) is calculated as

ξ=

1 Nb

Nb

∑ |P(xi) − Pfinal(xi)| i=1

(36)

where Nb is the number of bins used to construct the distributions. As was found for fixed-charge models,33 the large outer time step of the SIN(R) simulations does not affect the rate of convergence of the computed distribution functions. 3.3. Hydrogen Bond Angle Distribution Functions. We next analyze the hydrogen bond angle distribution functions. In a water−water hydrogen bond, let the donor oxygen be denoted OD, the donor hydrogen, HD, the acceptor oxygen, OA, and the acceptor hydrogen, HA. Using these four atoms, we define three angles: angle α is the OD−HD−OA angle, β is the OA−OD−HD angle, and θ is the HD−OA−HA angle. 2176


Article


Figure 7. Time evolution of the three contributions to the average molecular dipole moment. Contributions from static dipoles ⟨Pstatic(t)⟩ (left), induced dipoles ⟨Pinduce(t)⟩ (middle), and monopoles ⟨Pmono(t)⟩ (right) are shown.

diffusion coefficient for AMOEBA water at 298 K, namely, 1.96 × 10−5cm2/s.57 The SIN(R) simulations XO-SIN(R), XM-SIN(R), and XI-SIN(R), which are detailed in Table 1, yield selfdiffusion coefficients of 1.8 × 10−5, 1.7 × 10−5, and 1.8 × 10−5 cm2/s, at 298 K, respectively. These coefficients are not in perfect agreement with the NVE results, although they are within approximately 10% of the NVE value, showing that at least this dynamical property is only weakly perturbed by the isokinetic scheme. Surprisingly, SIN(R) with the largest outer time step seems to give a better diffusion constant than the approach of ref 32 at its largest outer time step.

Distribution functions for the hydrogen bond angles α, β, and θ were calculated for the NVT and SIN(R) simulations. These distributions are shown in Figure 4. The hydrogen bond angle distributions from the SIN(R) simulations agree well with those from the NVT benchmark simulations, showing again that the use of a large outer time step does not affect the accuracy of the results. The L1 norms of the differences between the distributions at time t and their converged values are shown in Figure 5. Again, it is observed that the convergence of these distributions is insensitive to the size of the outer time step in the SIN(R) simulations. 3.4. Molecular and Total Cell Dipoles. We next analyze the distribution of molecular dipoles, which is shown in Figure 6. This figure also shows the convergence of the average molecular dipole moment. The average of the molecular dipole moment from the NVT benchmark simulation is 2.79 D, while the values from the XO-SIN(R), XM-SINR(R), and XI-SINR(R) simulations are 2.76, 2.74, and 2.74 D, respectively. Although there is a slight shift in the distribution, as shown in Figure 6, the difference in the average is under 2%. This already small difference can be reduced. With an outer time step of 15 fs, the average molecular dipole moment from XM-SIN(R) and XI-SIN(R) simulations is 2.77 D. The origin of the small shift in the average molecular dipole moment can be analyzed by examining the different contributions to the overall molecular dipole moment: static dipoles, induced dipoles, and monopole terms of the form ∑iqiri, where i indexes the three atoms in water, qi is the charge of an atom, and ri is its position. In Figure 7, the convergence of these contributions over the simulation is shown. It is shown that the main source of disagreement is the induced dipole contribution: Compared to the NVT benchmark, the largest deviation from this value is a 4.2% underestimate from XMSIN(R). With a 15 fs outer time step, this deviation reduces to 1.1%. SIN(R) values for static dipole and monopole contributions each disagree with the NVT values by less than 1%. 3.5. Diffusion Coefficient. The question often arises as to how well SIN(R) simulations are able to reproduce dynamical properties. It is important to point out that SIN(R) is meant as a sampling algorithm and is not designed to reproduce dynamics. Nevertheless, the question is one worth exploring. Therefore, we computed self-diffusion coefficients of AMOEBA water for the NVT, NVE, and SIN(R) simulations. For this, we employed the usual Einstein relation Dself =

1 d 1 lim 6 t →∞ dt N

4. CONCLUSION We have extended our resonance-free multiple time-step algorithm, developed previously for fixed-charge models,33 to a polarizable model. In this approach, bonded interactions are treated as “fast”, short-range electrostatics and van der Waals interactions comprise an “intermediate” time scale, and longrange electrostatics and van der Waals interactions are treated as “slow.” At the intermediate level, we include a self-consistent solution of induced dipoles due only to neighboring dipoles within a short-range cutoff following a scheme introduced by Masella.51 Slow forces include not only long-range interactions in real space but also all reciprocal-space forces; hence, we term the scheme RESPA1-pol. As was observed for fixed-charge models,33 the new algorithm allows very large outer time steps (∼120 fs) to be employed, which greatly enhances the efficiency of the calculations. Depending on the choice of a long-range real-space cutoff, our approach can converge equilibrium properties in 5%− 10% of the computational time required for standard NVT molecular dynamics with our implementation in the TINKER package with negligible loss in accuracy. Surprisingly, we also find that diffusion constants are only minimally perturbed, suggesting that at least slow dynamics are largely preserved, at least qualitatively, although this issue merits further investigation. We expect that the force decomposition described for polarizable models can be applied to a wide variety of complex systems in biomolecular and materials modeling and that other types of RESPA force decompositions33,46 can be employed. The method can also be adapted for self-consistent solutions schemes that aim to improve on the conjugate gradient self-consistent field method, such as a hybrid self-consistent field/extended Lagrangian58 or a purely extended Lagrangian59 method. The method should also be applicable to other types of polarizable models such as the Drude model.60−65 Beyond classical molecular dynamics with force fields, we believe that the stochastic isokinetic multiple time-step integrator can be used for computational cost savings at higher levels of theory such as ab initio molecular dynamics66 and path-integral molecular dynamics.67

N

∑ ⟨|ri(t ) − ri(0)|2 ⟩ i=1

(37)

where i = 1,...,N indexes molecular centers of mass. From the mean square deviations, we obtain self-diffusion coefficients of 2.0 × 10−5 and 1.9 × 10−5 cm2/s for NVT and NVE simulations, respectively. These are in agreement with the previously reported 2177


Article


■

so that λα = λFα + λNα . If x denotes the full (physical and extended) set of phase-space variables, then the evolution of the system from an initial condition x0 to a phase-space vector xt at time t is given by xt = exp(iLt)x0, where exp(iLt) is the classical propagator. Of course, exp(iLt) cannot be evaluated analytically; however, for a single time step Δt, it can be approximated via a Trotter factorization, which then leads to a numerical solver for the equations of motion.21,22,33 Suppose the forces Fα(r) can be decomposed into fast, im intermediate, and slow components, i.e., Fα(r) = Ffast α (r) + Fα (r) slow + Fα (r). We wish to factorize the propagator such that time steps Δt, δt, and dt can be used for the slow, intermediate, and fast components, Δt = nδt = nmdt for integers n and m, where n,m ≥ 1. In order to generate such a factorization, we need to define corresponding contributions to the total Liouville operator, taking into account the fact that the Lagrange multipliers also depend on these three components of the force. With this in mind, and following the pattern of ref 33, we introduce the proceeding decomposition of the operator iLv as

APPENDIX

SIN(R) with Three Time Steps

In this Appendix, we derive a multiple time-step (MTS) integrator for the equations of motion (eq 1), based on the reference system propagator algorithm (RESPA) introduced by Tuckerman et al.21 The derivation begins with the Liouville operator for the equations of motion given by iL = iLr + iLv + iLN + iLOU

(A.1)

where iLr =

∑ vα α

iLv =

∂ ∂rα

⎡⎛

L ⎤ ⎞ ∂ Fα(r) ∂ − λαF vα⎟ − λαF ∑ v1,(kα) (k) ⎥ ⎠ ∂vα ∂v1, α ⎥⎦ ⎣⎝ mα k=1

∑ ⎢⎢⎜ α

iLN = −

∑ [λαN vα α

L

+

∑ v2,(kα)v1,(kα) k=1

∂ + λαN ∂vα ∂

∂v1,(kα)

L

−

∑ k=1

L

∑ v1,(kα) k=1

∂ ∂v1,(kα)

iLvfast =

G(v1,(kα)) ∂ ] Q 2 ∂v2,(kα)

α

iLvim =

(A.2)

λαF =

λαN

mvα2

=−

+

iLvslow =

α

⎣⎝ mα

L ⎤ ⎞ ∂ ∂ − λαimvα ⎟ − λαim ∑ v1,(kα) (k) ⎥ ⎠ ∂vα ∂v1, α ⎥⎦ k=1

⎡⎛ slow F (r)

∑ ⎢⎢⎜⎜ α

⎣⎝

α

mα

L ⎤ ⎞ ∂ ∂ − λαslow vα ⎟⎟ − λαslow ∑ v1,(kα) (k) ⎥ ∂v1, α ⎥⎦ ⎠ ∂vα k=1

(A.4) 2 Setting Λα = mv2α + [L/(L + 1)]∑Lk=1 Q1 (v(k) 1,α) , we have fast fast im im introduced λα = vαFα (r)/Λα, λα = vαFα (r)/Λα, and λslow = α F fast slow vαFslow + λim α (r)/Λα, so that λα = λα α + λα . With this decomposition, and following the suggestion of Matthews and Leimkuhler for the placement of the Ornstein−Uhlenbeck operator,68 the three time-step XO-RESPA factorization for SIN(R) then becomes

vαFα(r) L L ∑k = 1 Q 1(v1,(kα))2 L+1

L L ∑k = 1 Q 1(v1,(kα))2 v2,(kα) L+1 L L mvα2 + L + 1 ∑k = 1 Q 1(v1,(kα))2

⎡⎛ im F (r)

∑ ⎢⎢⎜ α

and iLOU corresponds to the Ornstein−Uhlenbeck stochastic process applied to v(k) 2,α. The Lagrange multiplier has been divided into force-dependent and thermostat contributions given, respectively, by

L ⎤ ⎞ ∂ Fαfast(r) ∂ − λαfastvα ⎟⎟ − λαfast ∑ v1,(kα) (k) ⎥ ∂v1, α ⎥⎦ ⎠ ∂vα ⎣⎝ mα k=1

⎡⎛

∑ ⎢⎢⎜⎜

(A.3)

⎛ Δt ⎞ ⎛ Δt ⎞⎧ ⎛ δt ⎞ exp(iLΔt ) = exp⎜iLN ⎟exp⎜iLvslow ⎟⎨exp⎜iLvim ⎟ ⎝ ⎝ ⎠ ⎝ ⎠ 2 2 ⎩ 2⎠ m ⎡ ⎛ ⎛ dt ⎞ ⎛ dt ⎞ ⎛ dt ⎞ dt ⎞⎤ × ⎢exp⎜iLvfast ⎟exp⎜iLr ⎟exp(iLOUdt )exp⎜iLr ⎟exp⎜iLvfast ⎟⎥ ⎝ 2⎠ ⎝ ⎣ ⎝ 2⎠ ⎝ 2⎠ 2 ⎠⎦

⎛ ⎛ δt ⎞⎫ Δt ⎞ ⎛ Δt ⎞ × exp⎜iLvim ⎟⎬ exp⎜iLvslow ⎟exp⎜iLN ⎟ ⎝ ⎠ ⎝ ⎭ 2 2 ⎠ ⎝ 2 ⎠ n

In this factorization, it can be seen that the deterministic part of the Nosé−Hoover operator is applied only when the slow forces are required, which, when large time steps are used, can be very infrequent and might not provide sufficient thermalization of the system. As a compromise, this operator can be pulled inside the RESPA loops and applied more frequently, which leads to the

(A.5)

schemes referred to as XM-RESPA and XI-RESPA in the main text. In the XM-RESPA scheme, the deterministic Nosé−Hoover operator is applied when the intermediate forces are required, which can be expressed via the following factorization:

2178


Article

Journal of Chemical Theory and Computation ⎛ δt ⎞ ⎛ Δt ⎞ ⎛ δt ⎞ exp(iLΔt ) = exp⎜iLN ⎟exp⎜iLvslow ⎟exp⎜ −iLN ⎟ ⎝ ⎠ ⎝ ⎝ ⎠ 2⎠ 2 2 ⎧ ⎛ δt ⎞ ⎛ δt ⎞ × ⎨exp⎜iLN ⎟exp⎜iLvim ⎟ ⎝ ⎠ ⎝ ⎩ 2 2⎠ m ⎡ ⎛ ⎛ dt ⎞ ⎛ dt ⎞ ⎛ dt ⎞ dt ⎞⎤ × ⎢exp⎜iLvfast ⎟exp⎜iLr ⎟exp(iLOUdt )exp⎜iLr ⎟exp⎜iLvfast ⎟⎥ ⎝ 2⎠ ⎝ ⎣ ⎝ 2⎠ ⎝ 2⎠ 2 ⎠⎦

⎛ δt ⎞ ⎛ δt ⎞⎫ × exp iLvim ⎟exp⎜iLN ⎟⎬ ⎝ 2⎠ ⎝ 2 ⎠⎭

n

⎜

⎛ Δt ⎞ ⎛ δt ⎞ δt ⎞ ⎛ × exp⎜ −iLN ⎟exp⎜iLvslow ⎟exp⎜iLN ⎟ ⎝ 2⎠ ⎝ 2 ⎠ ⎝ 2⎠

■

As discussed in ref 56, the operator exp(−iLNδt/2) is never explicitly applied. It is included in eq A.6 as a mathematical device to cancel the operator exp(iLNδt/2) in the curly braces on the first and last steps of the short-range RESPA loop so that the thermostats are applied only after the particle velocities are updated with the appropriate combination of forces. Similarly, the XI-RESPA scheme, in which the thermostats are applied when the fast forces are required, has the factorization form

⎧ ⎛ dt ⎞ ⎛ δt ⎞ ⎛ dt ⎞ × ⎨exp⎜iLN ⎟exp⎜iLvim ⎟exp⎜− iLN ⎟ ⎩ ⎝ 2⎠ ⎝ 2⎠ ⎝ 2⎠ ⎡ ⎛ dt ⎞ ⎛ dt ⎞ ⎛ dt ⎞ × ⎢exp⎜iLN ⎟exp⎜iLvfast ⎟exp⎜iLr ⎟ ⎝ ⎠ ⎝ ⎣ 2 2⎠ ⎝ 2⎠ ⎛ dt ⎞ ⎛ dt ⎞ ⎛ dt ⎞ × exp(iLOUdt )exp⎜iLr ⎟exp⎜iLvfast ⎟exp⎜iLN ⎟]m ⎝ 2⎠ ⎝ 2⎠ ⎝ 2⎠ n

⎛ dt ⎞ ⎛ Δt ⎞ ⎛ dt ⎞ × exp⎜ − iLN ⎟exp⎜iLvslow ⎟exp⎜iLN ⎟ ⎝ ⎝ ⎠ 2 2 ⎠ ⎝ 2⎠ (A.7)

The action of each of the operators appearing in these three factorization schemes can be worked out analytically, and the details are given in ref 33. As the operator structure is exactly the same as that of a standard thermostated RESPA code, as described in ref 56, all that is needed to implement the SIN(R) scheme is a modification of individual steps employing the formulas given in ref 33.

■

REFERENCES

(1) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141−151. (2) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (3) Faradjian, A. K.; Elber, R. J. Chem. Phys. 2004, 120, 10880−10889. (4) West, A. M. A.; Elber, R.; Shalloway, D. J. Chem. Phys. 2007, 126, 145104. (5) Henin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. J. Chem. Theory Comput. 2010, 6, 35−47. (6) Baumketner, A.; Shea, J. E. Theor. Chem. Acc. 2006, 116, 262−273. (7) Maragliano, L.; Vanden-Eijnden, E. J. Chem. Phys. 2008, 128, 184110. (8) Rosso, L.; Minary, P.; Zhu, Z.; Tuckerman, M. E. J. Chem. Phys. 2002, 116, 4389−4402. (9) Rosso, L.; Abrams, J. B.; Tuckerman, M. E. J. Phys. Chem. B 2005, 109, 4162−4167. (10) Abrams, J. B.; Tuckerman, M. E. J. Phys. Chem. B 2008, 112, 15742−15757. (11) Maragliano, L.; Vanden-Eijnden, E. Chem. Phys. Lett. 2006, 426, 168−175. (12) Chen, M.; Cuendet, M. A.; Tuckerman, M. E. J. Chem. Phys. 2012, 137, 1−11. (13) Steele, R. P. J. Chem. Phys. 2013, 139, 011102. (14) Tuckerman, M. E.; Parrinello, M. J. Chem. Phys. 1994, 101, 1316− 1329. (15) Luehr, N.; Markland, T. E.; Martinez, T. J. J. Chem. Phys. 2014, 140, 084116. (16) Fatehi, F.; Steele, R. P. J. Chem. Theory Comput. 2015, 11, 884− 898. (17) Steele, R. P. J. Phys. Chem. A 2015, 119, 12119−12130. (18) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1990, 93, 1287−1291. (19) Tuckerman, M. E.; Berne, B. J.; Rossi, A. J. Chem. Phys. 1991, 94, 1465−1469. (20) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1991, 94, 6811−6815. (21) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990−2001. (22) Tuckerman, M. E. Statistical Mechanics: Theory and Molecular Simulations; Oxford University Press, Cambridge, U.K., 2008. (23) Schlick, T.; Mandziuk, M.; Skeel, R. D.; Srinivas, K. J. Comput. Phys. 1998, 140, 1−29. (24) Sandu, A.; Schlick, T. J. Comput. Phys. 1999, 151, 74−113. (25) Ma, Q.; Izaguirre, J. A.; Skeel, R. D. SIAM J. Sci. Comput. 2003, 24, 1951−1973. (26) Izaguirre, J. A.; Reich, S.; Skeel, R. D. J. Chem. Phys. 1999, 110, 9853−9864. (27) Izaguirre, J. A.; Catarello, D. P.; Wozniak, J. M.; Skeel, R. D. J. Chem. Phys. 2001, 114, 2090−2098.

⎛ dt ⎞ ⎛ Δt ⎞ ⎛ dt ⎞ exp(iLΔt ) = exp⎜iLN ⎟exp⎜iLvslow ⎟exp⎜− iLN ⎟ ⎝ 2⎠ ⎝ 2 ⎠ ⎝ 2⎠

⎛ δt ⎞ ⎛ dt ⎞⎫ dt ⎞ ⎛ × exp⎜ − iLN ⎟exp⎜iLvim ⎟exp⎜iLN ⎟⎬ ⎝ 2⎠ ⎝ 2⎠ ⎝ 2 ⎠⎭

(A.6)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge funding from the National Science Foundation Grant No. CHE-1265889 for support of this work. D.T.M. is thankful to Minzhong Xu for useful discussions. 2179


Article

Journal of Chemical Theory and Computation (28) Minary, P.; Tuckerman, M. E.; Martyna, G. J. Phys. Rev. Lett. 2004, 93, 1−4. (29) Omelyan, I. P.; Kovalenko, A. J. Chem. Phys. 2011, 135, 234107. (30) Omelyan, I. P.; Kovalenko, A. J. Chem. Theory Comput. 2012, 8, 6−16. (31) Omelyan, I. P.; Kovalenko, A. Phys. Rev. E 2012, 85, 026706. (32) Morrone, J. A.; Markland, T. E.; Ceriotti, M.; Berne, B. J. J. Chem. Phys. 2011, 134, 014103. (33) Leimkuhler, B.; Margul, D. T.; Tuckerman, M. E. Mol. Phys. 2013, 111, 3579−3594. (34) Evans, D. J.; Morriss, G. P. Comput. Phys. Rep. 1984, 1, 299−343. (35) Evans, D. J.; Morriss, G. P. Statistical Mechanics of Nonequilibrium Liquids; Harcourt Brace Javanovich: London, 1980. (36) Zhang, F. J. Chem. Phys. 1997, 106, 6102−6106. (37) Minary, P.; Martyna, G. J.; Tuckerman, M. E. J. Chem. Phys. 2003, 118, 2510−2526. (38) Martyna, G. J.; Klein, M. L.; Tuckerman, M. E. J. Chem. Phys. 1992, 97, 2635−2643. (39) Leimkuhler, B.; Noorizadeh, E.; Theil, F. J. Stat. Phys. 2009, 135, 261−277. (40) Ponder, J. W. TINKER - software tools for molecular design, 2004. http://dasher.wustl.edu/ffe/downloads/guide.pdf (accessed April 2016). (41) Ren, P.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 5933−5947. (42) Tuckerman, M. E.; Liu, Y.; Ciccotti, G.; Martyna, G. J. J. Chem. Phys. 2001, 115, 1678−1702. (43) Rogers, L. C. G.; Williams, D. Diffusions, Markov Processes, and Martingales, Vol.2: Itô Calculus; John Wiley & Sons: New York, 1987. (44) Tuckerman, M. E.; Mundy, C. J.; Martyna, G. J. Europhys. Lett. 1999, 45, 149−155. (45) Ren, P.; Ponder, J. W. J. Comput. Chem. 2002, 23, 1497−1506. (46) Morrone, J. A.; Zhou, R.; Berne, B. J. J. Chem. Theory Comput. 2010, 6, 1798−1804. (47) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T. A.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8593. (48) Procacci, P.; Marchi, M.; Martyna, G. J. J. Chem. Phys. 1998, 108, 8799−8803. (49) Zhou, R.; Harder, E.; Xu, H.; Berne, B. J. J. Chem. Phys. 2001, 115, 2348−2358. (50) Tuckerman, M. E.; Yarne, D. A.; Samuelson, S. O.; Hughes, A. L.; Martyna, G. J. Comput. Phys. Commun. 2000, 128, 333−376. (51) Masella, M. Mol. Phys. 2006, 104, 415−428. (52) Smith, W. CCP5 Inform. Q. 1998, 46, 18−30. (53) Sagui, C.; Pedersen, L. G.; Darden, T. A. J. Chem. Phys. 2004, 120, 73−87. (54) Kolafa, J. J. Comput. Chem. 2004, 25, 335−342. (55) Wang, W.; Skeel, R. D. J. Chem. Phys. 2005, 123, 164107. (56) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117−1157. (57) Ren, P.; Ponder, J. W. J. Phys. Chem. B 2004, 108, 13427−13437. (58) Albaugh, A.; Demerdash, O.; Head-Gordon, T. J. Chem. Phys. 2015, 143, 174104. (59) Niklasson, A. M. Phys. Rev. Lett. 2008, 100, 123004. (60) Lamoureux, G.; Roux, B. J. Chem. Phys. 2003, 119, 3025−3039. (61) Lamoureux, G.; MacKerell, A. D.; Roux, B. J. Chem. Phys. 2003, 119, 5185−5197. (62) Whitfield, T. W.; Martyna, G. J. Chem. Phys. Lett. 2006, 424, 409− 413. (63) Jones, A. P.; Crain, J.; Cipcigan, F. S.; Sokhan, V. P.; Modani, M.; Martyna, G. J. Mol. Phys. 2013, 111, 3465−3477. (64) Lemkul, J. A.; Roux, B.; van der Spoel, D.; MacKerell, A. D. J. Comput. Chem. 2015, 36, 1473−1479. (65) Sokhan, V. P.; Jones, A. P.; Cipcigan, F. S.; Crain, J.; Martyna, G. J. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 6341−6346. (66) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press: Cambridge, U.K., 2009. (67) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. J. Chem. Phys. 1993, 99, 2796−2808.

(68) Leimkuhler, B.; Matthews, C. Appl. Math. Res. Exp. 2013, 2013, 34−56.

2180


A Stochastic, Resonance-Free Multiple Time-Step Algorithm for

Recommend Documents