Elsevier

Chemical Physics

Volume 370, Issues 1–3, 12 May 2010, Pages 294-305
Chemical Physics

Measure-preserving integrators for molecular dynamics in the isothermal–isobaric ensemble derived from the Liouville operator

https://doi.org/10.1016/j.chemphys.2010.02.014Get rights and content

Abstract

The Liouville operator approach is employed to derive a new measure-preserving geometric integrator for molecular dynamics simulations in the isothermal–isobaric (NPT) ensemble. Recently, we introduced such a scheme for NPT simulations with isotropic cell fluctuations in the absence of holonomic constraints [M.E. Tuckerman et al., J. Phys. A 39 (2006) 5629]. Here, we extend this approach to include both fully flexible cell fluctuations and holonomic constraints via a new and simpler formulation of the ROLL algorithm of Martyna et al. [Martyna et al., Mol. Phys. 87 (1996) 1117]. The new algorithm improves on earlier schemes in that it possesses a simpler mathematical structure and rigorously preserves the phase space metric. The new algorithm is illustrated on two example systems, ice and liquid n-decane.

Graphical abstract

Density distribution of n-decane at 333 K and 1 atm using NPT(iso).

  1. Download : Download full-size image

Introduction

Over the last several decades, molecular dynamics (MD) has become one of the most important and commonly used approaches for studying condensed phase systems. MD calculations generally serve two often complementary purposes. First, MD can be employed as a means of generating a collection of classical microscopic configurations in a particular equilibrium ensemble. This application shows that MD is intimately connected with statistical mechanics and can serve as a computational tool for solving statistical mechanical problems. Second, an MD simulation can be used to study dynamical properties such as time correlation functions. Formally, a time correlation function should be computed first by sampling an ensemble of initial conditions from an equilibrium distribution and then launching a trajectory from each one. These points underscore the importance of having efficient and rigorous techniques capable of generating equilibrium distributions.

Starting from the most fundamental equilibrium ensemble, the microcanonical or NVE ensemble, which is generated simply by integrating Hamilton’s equations of motion, the most popular methods for generating a canonical (NVT) ensemble build on the equations for the microcanonical ensemble by introducing a mechanism for generating the proper energy fluctuations and ensuring that the average kinetic energy is consistent with the desired target temperature T. Such methods include Langevin dynamics, Nosé–Hoover dynamics [1], the Nosé-Poincaré algorithm [2], and the Nosé–Hoover chain scheme [3]. The last three of these fall into the class of techniques known as “extended phase-space” methods, in which the physical phase space is supplemented by additional dynamical variables that couple to the physical degrees of freedom in a manner that produces the desired kinetic energy fluctuations via a carefully designed set of equations of motion. These equations can be either Hamiltonian or non-Hamiltonian by construction. Devising robust numerical integrators for such approaches has proved to be relatively straightforward even for systems subject to holonomic constraints [4], and consequently, the canonical ensemble is often the ensemble of choice for many studies and serves as a common starting point for the development of enhanced configurational sampling methods. However, it is often necessary and/or desirable to study a system under conditions of constant pressure P and temperature T. Together with constant particle number N, these conditions define the isothermal–isobaric or NPT ensemble. The NPT ensemble is often used for purposes of equilibrating a system, allowing it to adjust to an appropriate density, for computing equilibrium properties under isobaric conditions, such as the Gibbs free energy or average density, and for studying structural phase transitions.

Among the equilibrium ensembles, the NPT ensemble is somewhat more difficult to generate as both the correct kinetic energy and instantaneous pressure fluctuations must be generated as prescribed by the ensemble distribution function, which means that both the total energy and the volume must be allowed to vary. If the volume is restricted to isotropic fluctuations only, then the ensemble distribution functionf(x,V)=1N!V0h3Nexp[-β(H(x)+PV)]Δ(N,P,T)is generated. Here, P,V and T are the applied external pressure, the system volume, and the external temperature, respectively, β=1/kT, x is the phase space vector, i.e. the vector whose components are the momenta and coordinates of all of the particles in the system, H(x) is the system Hamiltonian, V0 is an arbitrary reference volume, and h is Planck’s constant. The quantity Δ(N,P,T) is the partition function of the ensemble given byΔ(N,P,T)=1N!V0h3N0dVdxexp[-β(H(x)+PV)]=1V00dVe-βPVQ(N,V,T).The 1/N! combinatorial prefactor in Eqs. (1), (2) correspond to one-component systems but can be easily generalized for multi-component systems. If the Cartesian momenta and coordinates of the N particles are denoted p1,,pN,r1,,rNp,r, respectively, then Eq. (2) becomesΔ(N,P,T)=1N!V0h3N0dVD(V)dNrdNp[-β(H(p,r)+PV)]where each of the coordinate integrations is restricted by the spatial domain D(V) defined by the volume at each value of the volume integration.

It is also possible to incorporate anisotropic volume fluctuations in the NPT ensemble. This is done by introducing the box matrix or cell matrix h whose columns are the three cell vectors a, b and c. The volume is determined from V=a·b×c=det(h). The partition function for this case can be determined by separating isotropic and anisotropic fluctuations. To this end, we let h0 denote a unit box matrix with det(h0)=1 so that h=V1/dh0 in d dimensions. We can, therefore, write the partition function asΔ(N,P,T)=1V00dVdh0e-βPVQ(N,V,h0,T)δ(det(h0)-1).We now change variables from h0 to h using dh0=(V-d)dh, which givesΔ(N,P,T)=1V00dVdhV-de-βPVQ(N,h,T)δ1Vdet(h)-1=1V00dVdhV1-de-βPVQ(N,h,T)δ(det(h)-V)=1V0dh[det(h)]1-de-βPdet(h)Q(N,h,T).Note the additional factor of [det(h)]1-d in the integration measure. When designing molecular dynamics algorithms for the NPT ensemble with anisotropic cell fluctuations, it is essential that both the equations of motion and numerical solver generate this factor as part of the phase-space volume element.

As with the canonical ensemble, the most commonly used algorithms for generating the NPT ensemble are based on the extended-phase-space approach [5], [1], [6], [7], [8], in which the volume is introduced as an additional dynamical variable, along with a corresponding momentum variable, in order to maintain isobaric conditions, while additional thermostat variables [9], [1], [3], [10], [11], [12] are used to control the kinetic energy fluctuations. Additional mass-like parameters are assigned to each of these additional variables, and it is these parameters that determine the time scales on which the additional variables evolve. Equations of motion are specified as a means of generating the phase-space distribution of the ensemble. Given that simulations in the NPT ensemble are typically more challenging than in the NVE or NVT ensembles, it is important to have robust numerical integrators for the equations of motion. Unlike the NVE and NVT ensemble cases, this has proved more challenging for the NPT ensemble, particularly when the system is subject to a set of holonomic constraints. If non-Hamiltonian equations of motion are employed, as is often the case, then the numerical integrator should be consistent with a non-Hamiltonian generalization of Liouville’s theorem [13], [14]. This generalization attaches a phase-space metric factor to the usual volume element that must be preserved, and a numerical integrator that achieves this is known as a “measure-preserving” algorithm. Note that measure-preserving algorithms are also important when non-Hamiltonian dynamical systems are used to generate the NVT ensemble [15]. In this paper, we briefly review and analyze equations of motion for generating the NPT ensemble and then discuss the problem of constructing a measure-preserving numerical integrator for the equations of motion that can be easily adapted for systems subject to holonomic constraints.

Of the various dynamical schemes that have been proposed for generating the NPT ensemble, we will focus our attention on the algorithm of Martyna, Tobias and Klein (MTK) [7], which has been shown to correctly reproduce the isothermal–isobaric distribution for both isotropic and anisotropic volume fluctuations [7], [14]. We point out, however, that Hamiltonian formulations of the NPT ensemble have also been presented [8], although these will not be discussed here. Algorithms for integrating the MTK equations of motion introduced previously, including a scheme known as “ROLL” for enforcing holonomic constraints [4], have proved very useful, but unfortunately, they are not strictly measure preserving. Recently, we introduced a measure-preserving algorithm for molecular dynamics simulations in the NPT ensemble for isotropic volume fluctuations in the absence of holonomic constraints [16] and applied the approach in simulations of aqueous solutions of ethanolamines [17]. A multiple time-scale version of the integrator based on the r-RESPA [18] approach was also presented [16]. Here, we extend this approach to include anisotropic volume fluctuations and holonomic constraints via a new and simpler version of the ROLL algorithm [4]. The new algorithm is generally more straightforward to implement than that of Ref. [4], especially given a working algorithm for the canonical (NVT) ensemble employing Nosé–Hoover chains [3]. The new approach will be illustrated by two representative systems, namely ice and liquid n-decane.

Section snippets

Equations of motion

We begin by considering a classical system of N particles in d dimensions with positions r1,,rN and momenta p1,,pN in a general container of volume V described by a cell matrix h. In d=3 dimensions h would describe a generally triclinic cell. Let U(r1,,rN) denote the interparticle potential and Fi=-U/ri denote the corresponding interparticle forces. For the moment, we consider systems in which no constraints are present. The MTK [7] algorithm for generating the isothermal–isobaric

Measure-preserving integrators

It has been shown that numerical integrators for Hamiltonian systems gain special stability when they are symplectic [25]. Symplectic integrators also preserve the phase space measure and, therefore, obey Liouville’s theorem. To our knowledge, no generalization of the symplectic condition for general non-Hamiltonian systems exists at this time. However, we can still insist that a numerical solver for a non-Hamiltonian system preserve the generalized phase-space volume g(xt)dxt, which renders

Holonomic constraints in the NPT ensembles

Holonomic constraints are commonly used in MD simulations as a means of freezing out high-frequency vibrations, thereby allowing a larger time step to be employed. The combination of the SHAKE [32] and RATTLE [33] algorithms is designed to work in conjunction with the velocity-verlet method. In the NVT ensemble, these constraint methods can be implemented exactly as they are in the NVE ensemble: The SHAKE step is called immediately after the positions are updated, and the RATTLE step is called

Simulation details

As a test, the new constraint NPT algortihm is applied to study ice Ih at T=250K and Pext=1atm. The TIP4P [36] model, which is a four-site model including two hydrogens, one oxygen and one ghost atom, is chosen. The ghost atom is placed 0.15 Å away from the oxygen along the H–O–H bisector and has a negative charge of -1.04e. The oxygen and the hydrogen have point charges 0e and +0.52e respectively. Bond constraints are applied to the O–H bonds and the H–H distance in the simulation to keep the

Results and discussion

In this work, both NPT(iso) and NPT(flex) simulations have been run for ice Ih at T=250K and Pext=1atm. The stability of the algorithms is measured by the conserved quantity Econv=1Nci=1Nc|Ei-E0|/E0, where i runs over the number of configurations Nc, and Ei is the total energy of the system including the atoms, barostat and thermostat variables, at the ith step. In Fig. 1(a) and Fig. 3(a), the instantaneous and cumulative averages of the energy conservation measure are shown for the NPT(iso)

Acknowledgments

The authors wish to acknowledge Michael Shirts for helpful discussions. This work was supported by NSF CHE-0704036 and DOE FG05-08OR23334. J.A. thanks Conacyt for financial support.

References (44)

  • S.D. Bond et al.

    J. Comput. Phys.

    (1999)
  • S.D. Bond et al.

    J. Comput. Phys.

    (1999)
  • V.E. Tarasov

    Ann. Phys.

    (2005)
  • H. Yoshida

    Phys. Lett. A

    (1990)
  • J.P. Ryckaert et al.

    J. Comput. Phys.

    (1977)
  • H.C. Andersen

    J. Comput. Phys.

    (1983)
  • M.E. Tuckerman et al.

    Comput. Phys. Commun.

    (2000)
  • W.G. Hoover

    Phys. Rev. A

    (1985)
  • G.J. Martyna et al.

    J. Chem. Phys.

    (1992)
  • G.J. Martyna et al.

    Mol. Phys.

    (1996)
  • H.C. Andersen

    J. Chem. Phys.

    (1980)
  • W.G. Hoover

    Phys. Rev. A

    (1986)
  • G.J. Martyna et al.

    J. Chem. Phys.

    (1994)
  • J.B. Sturgeon et al.

    J. Chem. Phys.

    (2000)
  • S. Nosé

    J. Chem. Phys.

    (1984)
  • Y. Liu et al.

    J. Chem. Phys.

    (2000)
  • B.J. Leimkuhler et al.

    J. Chem. Phys.

    (2004)
  • M.E. Tuckerman et al.

    Europhys. Lett.

    (1999)
  • M.E. Tuckerman et al.

    J. Chem. Phys.

    (2001)
  • M.E. Tuckerman et al.

    J. Chem. Phys.

    (1999)
  • M.E. Tuckerman et al.

    J. Phys. A: Math. Gen.

    (2006)
  • R. Lopez-Rendom et al.

    J. Phys. Chem. B

    (2006)
  • Cited by (52)

    View all citing articles on Scopus
    View full text