Article pubs.acs.org/JCTC
Continuous Tempering Molecular Dynamics: A Deterministic Approach to Simulated Tempering Nicolas Lenner‡ and Gerald Mathias* Lehrstuhl für BioMolekulare Optik, Ludwig−Maximilians Universität München, Oettingenstrasse 67, 80538 München, Germany S Supporting Information *
ABSTRACT: Continuous tempering molecular dynamics (CTMD) generalizes simulated tempering (ST) to a continuous temperature space. Opposed to ST the CTMD equations of motion are fully deterministic and feature a conserved quantity that can be used to validate the simulation. Three variants of CTMD are discussed and compared by means of a simple test system. The implementation features of the most stable and simplest variant CTMD-< in the program package IPHIGENIE are described. Two applications − alanine dipeptide (Ac-Ala-NHMe) in explicit water and octa-alanine (Ac-(Ala)8NHMe) simulated in a dielectric continuum − demonstrate the functionality of CTMD-< . Furthermore, they serve to evaluate its sampling efficiency. Here, CTMD-< outperforms ST by 35% and replica exchange even by 75%.
1. INTRODUCTION Sampling free energy landscapes of complex systems by plain molecular dynamics (MD) simulations is often hampered by large barriers separating regions of interest that are crossed too seldom or not at all within the accessible simulation time. Relevant examples for this so-called sampling problem span the whole field of MD simulations from chemical reactions, the conformational dynamics and folding of proteins, to first-order phase transitions, just to name a few.1 Numerous methods have been proposed to enhance the sampling efficiency of MD simulations.2−23 A large class of such methods relies on a priori knowledge of the system to guide the selection of relevant reaction coordinates and their description by collective variables (CV). Sampling is then accelerated within the subspace spanned by these CV but also the free energy landscape can only be determined in this a priori chosen space. For one-dimensional free energies depending on a single CV or on a predefined reaction path methods like umbrella sampling,2,3 thermodynamic integration4,11 and λ-dynamics5,6 have emerged as prominent techniques. Higher-dimensional problems including several CV can be tackled by methods flattening the free energy landscape by a bias potential like Wang−Landau Sampling7 and filling techniques such as conformational flooding8 and metadynamics.9,10 Furthermore, an accelerated sampling of CV can be achieved by a stochastic dynamics in temperature accelerated MD12 or an extended phase space approach with adiabatic dynamics (AFED).13−15 A second large class of enhanced sampling methods are generalized ensemble methods such as temperature replica exchange(RE)-MD,16,17 simulated tempering (ST),18,19 and their solute tempering modifications.20,21 Their common idea is to simulate on a discrete set of so-called rungs l with temperatures Tl that span a ladder T0 < T1 < ... < TL−1 of size L. At high temperatures the system’s degrees of freedom uniformly © XXXX American Chemical Society
have more kinetic energy available and, thus, can surmount enthalpic barriers more easily. In RE as many replicas of the system are simulated as there are rungs on the ladder, and these rungs are occupied exclusively by a single replica. To ensure canonical sampling a proper thermostat such as Bussi,24 Langevin,25 or Nosé−Hoover Chain (NHC) thermostat26−29 has to be employed. Stochastic exchanges between replicas on different rungs are attempted periodically with an exchange probability given by the Metropolis criterion30 that preserves the canonical ensemble at each temperature Tl. An approach systematically similar to RE is ST.18,19 Here, replicas move independently between the rungs of the temperature ladder. The single replica31 ST partition function L−1
ZST =
∑ Z(Tl ) exp wl
(1)
l=0
of the generalized ensemble spanned by L rungs of temperatures Tl is the sum of the partition functions Z(Tl) of rungs l with a weighting factor exp(wl) depending on dimensionless weight parameters wl. For a system with potential energy U = U(r1, ..., rN) these weights wl modify the ST Metropolis criterion for the transition probability Γkl = min{1, exp[(βk − βl )U − (wk − wl)]}
(2)
from rung k to rung l with inverse temperatures βl = (kBTl)−1. Because the difference of the weights shifts the transition probability Γkl it equally steers the overall occurrence Pl of rung l due to the detailed balance condition PkΓkl = PlΓlk met by the Metropolis criterion (2). A uniform sampling of all rungs on the given temperature ladder is achieved by choosing Received: August 6, 2015
A
DOI: 10.1021/acs.jctc.5b00751 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation weights wl = βlF(Tl) according to the free energies F(Tl) at each temperature.18 In order for eq 2 to be valid, the momenta pi →
Tl p Tk i
considerations such as the construction of the necessary background potentials and the generation of ensembles for discrete temperature values by an on-the-fly reweighting scheme. First results of a small system of particles in an external onedimensional potential serve to demonstrate the practical usability of CTMD and help to pick the numerically most stable variant CTMD-< , which has been implemented in the molecular mechanics (MM)-MD program IPHIGENIE.39 It is employed to compute the free energy landscape of the alanine dipeptide molecule (Ac-Ala-NHMe) in water. This first minimal reference application serves to compare the results of ST and CTMD-< . A second and more challenging CTMD application determines the melting curves of the α-helical peptide octa-alanine (Ac-(Ala)8-NHMe) in a dielectric continuum described by the Hamiltonian dielectric solvent (HADES) approach.40,41 Here, we compare the sampling speeds of ST, RE, and CTMD. For the latter our interest particularly rests on the dependence of the sampling speed on the mass of the tempering particle.
(3)
of all atoms i and thereby the kinetic energy of the system have to be adapted to the new ensemble temperature Tl upon a successful exchange from rung k to l by simple scaling. If a NHC thermostat is used in the simulation its momenta pηj have to be scaled as well by the same factor. The sampling efficiency of RE and ST simulations depends on the a priori choice of parameters like the number of rungs, the rung spacing, the period of exchange tries, and the exchange algorithm.32,33 Because the ST dynamics depends only on a single replica as opposed to RE, where all multiple replicas are coupled, ST generally requires less temperature rungs than RE.21,34 Setting up ST and RE simulations favorably is not always straightforward, especially when first-order phase transitions or steep ascents within the free energy F(T) occur. Various estimates and schemes for proper initial choices as well as on-the-fly optimizations of the parameters have been suggested.33,35−37 Replacing temperature rungs by a single dynamic continuous tempering variable would make the discrete tempering parameters that specify rungs and exchanges between rungs obsolete and should be advantageous. For this purpose, however, one has to devise an alternative scheme of how the system travels through the now continuous temperature space while retaining the correct canonical distributions for each temperature. Recently, Zhang, Ma, and co-workers proposed a diffusive stochastic process in temperature space22,23,38 governed by a Langevin equation of the form d(1/β) ∂ ln ϕ(β) 2 + = U − Ũ (β) − ξ dt ∂β β
2. THEORY For a continuous version of ST we start from the canonical partition function Z (T ) =
∫
⎡ N ⎞⎤ pi2 1 ⎛⎜ ⎢ ⎟⎟⎥ d pd r exp − + U ( r , ..., r ) ∑ 1 N ⎢⎣ kBT ⎜⎝ i = 1 2mi ⎠⎥⎦ N
N
(5)
for N atoms i with masses mi at positions ri and with conjugated momenta pi of one ST rung at temperature T, which was introduced in the partition function (1) of the ST generalized ensemble. For CTMD this sum has to be converted to an integral over the continuous temperature T, which then becomes a dynamic variable of its own right. Because we chose an extended system approach, this new tempering variable T is complemented by a conjugated momentum pT and a mass RT. Note that the stochastic process (4) does not involve a momentum variable. In the resulting partition function
(4)
as a first continuous tempering alternative to ST. Here, the time evolution of the continuous tempering variable β = 1(kBT)−1 is driven by the difference of the system’s potential energy U to the average potential energy Ũ (β), a force arising from a weighting function ϕ(β) and a standard Gaussian noise ξ. In this article we pursue an extended system approach to continuous ST which provides an alternative to the stochastic approach (4). We augment the simulation system by a massive tempering “particle” together with a corresponding momentum whose dynamics are driven by deterministic equations of motions (EOM). This particle is coupled to the system’s degrees of freedom, and it moves in an external background potential W that takes over the role of the ST weights wi (eq 1). We call this approach continuous tempering molecular dynamics (CTMD). The CTMD equations are time reversible and yield a conserved energy that can be used to check the stability of the implementation and of the simulation setup. In the following Theory section we proceed by introducing the basic theoretical concepts and construct three variants of CTMD CTMD-; , CTMD-7 , and CTMD-< including their respective EOM. Conserved quantities help us to demonstrate that the CTMD EOM sample canonical momentum and position distributions of the system. Test implementations of the CTMD variants are described, along with practical
ZCT =
∫
2 ⎡ ⎞⎤ 1 ⎛⎜ pT ⎟⎟⎥ dT dpT Z(T ) exp⎢ − + W ( T ) ⎢⎣ kBT ⎜⎝ 2RT ⎠⎥⎦
(6)
designated for CTMD the summation over the rungs in ZST (1) has become the integral over dT and dpT, whereas the newly introduced background potential W(T) replaces the ST weights and has to be specified later. The next step is to devise a dynamics that will sample the canonical distributions of atomic positions and momenta for all temperatures T within a given interval and, thus, actually yields the desired ZCT (6). A prerequisite is a suitable thermostatting algorithm that yields a canonical distribution of atoms and their momenta for a given ensemble temperature T0. Here, the NHC thermostat serves as a blueprint for our construction of the CTMD EOM, which are non-Hamiltonian but feature a conserved quantity H. We provide a short recap in Section S7 of the Supporting Information (SI), which focuses on the technical aspects of the NHC construction that are utilized analogously for CTMD. B
DOI: 10.1021/acs.jctc.5b00751 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation The NHC EOM29 pi
ri̇
=
pi̇
=−
ηj̇ =
,
mi
1
pη̇
=∑ i=1
=
j′
pη̇
M
=
pη
j
Qj
pη ∂U − 1 pi Q1 ∂ri
N
pη̇
describe the tempering particle together with a conjugated momentum pρ, a mass Rρ of dimension [Rρ] = energy × time2, and a background potential W; (ρ). The variable ρ yields the temperature T(ρ) = T0 exp(−2ρ). The detailed construction of the EOM
pi2
− fkBT0 −
mi
pη2
j ′− 1
Q j ′− 1
− kBT0 −
pη2 M−1 Q M−1
pη
2
Q2
pη
pη
j ′+ 1
Q j ′+ 1
1
N
HNHC =
∑ i=1
2mi
pi̇
=−
pη̇
=∑
pη
1
,
mi
j′
− kBT0
Q j ′− 1
M
∑ j=1
pη2 j
⎛ + kBT0⎜⎜fη1 + 2Q j ⎝
M
⎞
k=2
⎠
∑ ηk⎟⎟
pη2
M−1
⎛ pη pρ ⎞ j ′+ 1 ⎟p − kBT (ρ) − ⎜⎜ + Rρ ⎟⎠ ηj ′ ⎝ Q j ′+ 1
pη̇
=
pρ̇
⎡N p2 ⎢ =2⎢∑ i + ⎢⎣ i = 1 2mi
M
Rρ
⎛ pη pρ ⎞ ⎟⎟p − fkBT (ρ) − ⎜⎜ 2 + mi Rρ ⎠ η1 ⎝ Q2
j ′− 1
(7)
pρ
ρ̇ =
,
pi2
pη2 =
j
Qj
⎛ pη pρ ⎞ ∂U ⎟⎟p − ⎜⎜ 1 + ∂ri Rρ ⎠ i ⎝ Q1
i=1
j′
pη̇
+ U (r1, ..., rN ) +
ηj̇ =
=
N
with i = 1, ... ,N, j = 1, ..., M and j′ = 2, ..., M−1 conserve the energy29,42,43 pi2
pη
pi
ri̇
Q M−1
pρ
− kBT (ρ) − M
∑ j=1
Rρ
pη , M
pη2
⎛ + kBT (ρ)⎜⎜fη1 + 2Q j ⎝ j
⎞⎤ ∂W (ρ) η ∑ k⎟⎟⎥⎥ − ; ∂ρ ⎠ k=2 ⎥⎦ M
(9)
(8)
with j = 1, ..., M and j′ = 2, ..., M−1 is carried out in Section S8.1 of the SI and is cross-referencing the equivalent technical steps of the NHC construction in Section S7. The velocities ṙi of the atoms and η̇j of the thermostatting particles on the first line are the same in the canonical NHC scheme 7. The forces acting on the pi and the pηj on the r.h.s. of lines 2−5 of eq 9 are extended by additional coupling terms ∼pρ/Rρ, similar to the coupling to the following NHC momenta pηj+1 of the chain. Furthermore, the temperature T in the NHC equation is now depending on ρ. The force acting on the ρ particle appearing on the r.h.s. of the last line of (9) depends on the kinetic energy of the system and of the NHC, the potential energy of the NHC, and on the background potential W; , which is arbitrary at this point and will be specified later. The EOM (9) conserve the quantity
where the variables ηj and momenta pηj represent the M heat− bath−chain particles of mass Qj that couple to the f dynamic degrees of freedom of the atoms. The energy conservation is easily shown by verifying dHNHC/dt = 0 with the help of the EOM (7). Despite this conserved quantity HNHC the EOM (7) are non-Hamiltonian and, thus, cannot be derived from (8) by Hamilton’s formalism. Specifying EOM and a conserved quantity, however, suffices to determine the sampling properties of the associated dynamics. The necessary mathematical framework43 as well as its application to the NHC equations43 are thoroughly discussed in a recent textbook by M. E. Tuckerman.44 Similarly, the EOM and conserved quantities H, which will be provided below for the three CTMD variants, are sufficient to specify its properties. The more technical details of the construction are left to the SI, Sections S8, S9, and S10. The η-representation of the NHC (7and 8) are the most convenient form of the NHC EOM.29 These variables originate from the transformation ηj = ln sj of the fundamental NHC scaling variables sj.29 Here, the start of the chain s1 is the variable that scales the atomic momenta pi if the system is cooler or warmer than the target temperature T0, as in the original Nosé algorithm.27 2.1. CTMD-; : Temperature Dependent Momenta. In order include new dynamical variables T and pT in the NHC EOM (7) one has to complement the conserved quantity (8) with the kinetic and potential energy terms appearing in the exponential of (6). Furthermore, in the EOM the momenta pi and pηj of the atoms and the NHC, respectively, have to couple to T, since their distributions have to depend on the temperature. Using the NHC equations as a blueprint the construction of the CTMD EOM is straightforward upon introducing a tempering variable γ = T0/T with respect to a reference temperature T0. Here, γ has a close correspondence to the NHC scaling variable s1 and is much more convenient to use than T itself. Furthermore, in analogy to the transformation η1 = ln s1 used for the NHC, γ is transformed to the variable ρ = ln γ, which is ultimately used to
N
H; =
∑ i=1
+
pi2 2mi pρ2
2Rρ
M
+U+
∑ j=1
pη2
⎛ + kBT (ρ)⎜⎜fη1 + 2Q j ⎝ j
M
⎞
∑ ηk⎟⎟ k=2
⎠
+ W; (ρ) (10)
which is easily verified by showing dH; /dt = 0 with the help of the time derivatives in (9). Having constructed the EOM (9), we need to determine the probability distribution of atomic positions and momenta which is sampled by the corresponding dynamics in the simulation. In Section S8.2 of the SI we follow the construction scheme for non-Hamiltonian EOM43 to obtain the partition function Z ; of the EOM (9), which yields the unnormalized distribution 2 ⎛ ⎞ N pi ⎜ ∑i = 1 2mi + U (r1, ..., rN ) ⎟ ϱ; (p, r | ρ) = exp⎜ − ⎟ kBT (ρ) ⎜ ⎟ ⎝ ⎠
(11)
for a fixed temperature T = T(ρ) by proper integration. According to eq 11 the coordinates and momenta sampled by the C
DOI: 10.1021/acs.jctc.5b00751 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation The unnormalized distribution
EOM (9) are canonically distributed, both with respect to the chosen temperature T(ρ). We call this first approach ‘temperature dependent momenta’ and use the label CTMD-; for it. Thus, EOM (9) yield a partition function similar to ZCT (eq 6) discussed at the beginning of Section 2, but on the downside the EOM (9) are rather complicated, particularly because of the force on the atomic momenta ṗi depends both on the NHC thermostat and on the tempering variable ρ. 2.2. CTMD-7 : Temperature Invariant Momenta. The strong dependence of the atomic momenta on ρ and η1 of CTMD-; can be relaxed by expressing the scaling due to the tempering variable explicitly, i.e. by the substitution p → λ p, λ = T/T0, for all atomic and NHC momentum variables in (8), which is analogous to the ST momentum scaling (3). Correspondingly, we call this second approach “momentum scaling” and label it CTMD-7 . With the new tempering variable λ for the tempering particle and its conjugated momentum pλ and mass Rλ the EOM λ pi
λpη
λ̇ =
pi̇
=−
pη̇
⎛ N p2 pη ⎞ =λ⎜⎜∑ i − fkBT0 − 2 pη ⎟⎟ Q2 1⎠ ⎝ i = 1 mi
pη̇
⎛ p2 ⎞ pη j ′+ 1 ⎜ ηj ′− 1 ⎟ =λ ⎜ − kBT0 − pη ⎟ ⎜ Q j ′− 1 Q j ′+ 1 j ′ ⎟ ⎝ ⎠
1
j′
pη̇
M
pλ̇
Qj
,
(14)
for a given temperature T(λ) = λT0 resulting from the EOM (12) differs from the distribution (11) of the CTMD-; EOM (9). Now, only the atomic coordinates ri are distributed according to the scaled temperature T(λ), whereas the atomic momenta are kept at the reference temperature T0. For all practical purposes this is, however, not a significant drawback, since correctly T(λ)distributed atomic momenta are obtained by the transformation pi → λ pi . 2.3. CTMD-< : Scaling the Potential. Although the 2-fold dependence of the momenta on the tempering as well as the thermostatting variables is avoided by using temperature invariant momenta, the EOM (12) are strongly modified with respect to the plain NHC EOM (7) with a λ scaling factor employed to almost every expression. A third variant of CTMD equations is obtained when the factor 1/λ within the Boltzmann factor exp(−H7 /(kBλT0)) is multiplied upon the atomic and NHC kinetic and potential energy terms of H7 , which, as a result, cancels most of the λ-factors and leaves us with a scaled potential U/λ. Thus, we call this third variant CTMD-< . Complementing these scaled terms anew with the kinetic energy of a λ particle and a background potential W< yields
pλ
=
mi
ηj̇ =
j
ri̇
,
⎡ ∑ p2 /(2m ) ⎤ ⎡ U (r1, ..., rN ) ⎤ i i ⎥ exp⎢ − ϱ7 (p, r | λ) = exp⎢ − ⎥ kBT0 kBT (λ) ⎦ ⎢⎣ ⎥⎦ ⎣
Rλ
pη ∂U − 1 λ pi ∂ri Q1
N
H< =
i=1
⎛ p2 ⎞ η =λ⎜ M − 1 − kBT0⎟ ⎜Q ⎟ ⎝ M−1 ⎠ ⎡N p2 ⎢ = − ⎢∑ i + ⎢⎣ i = 1 2mi
M
∑ j=1
pη2
⎛ j + kBT0⎜⎜fη1 + 2Q j ⎝
⎞⎤ ∂W (λ) ∑ ηk⎟⎟⎥⎥ − 7 ∂λ ⎠ k=2 ⎥⎦ M
H7 = λ ∑ i=1
+
2mi
pλ2 2Rλ
M
+ U + λ∑
+ W7(λ)
j=1
pη2
⎛ + λkBT0⎜⎜fη1 + 2Q j ⎝ j
M
pi
ri̇
=
⎛ + kBT0⎜⎜fη1 + 2Q j ⎝
∑
j
j=1
+ W
