Measure-preserving integrators for molecular dynamics in the isothermal–isobaric ensemble derived from the Liouville operator
Graphical abstract
Density distribution of n-decane at 333 K and 1 atm using NPT(iso).
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 functionis generated. Here, and T are the applied external pressure, the system volume, and the external temperature, respectively, , 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, is the system Hamiltonian, is an arbitrary reference volume, and h is Planck’s constant. The quantity is the partition function of the ensemble given byThe 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 , respectively, then Eq. (2) becomeswhere each of the coordinate integrations is restricted by the spatial domain 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 . The partition function for this case can be determined by separating isotropic and anisotropic fluctuations. To this end, we let denote a unit box matrix with so that in d dimensions. We can, therefore, write the partition function asWe now change variables from to h using , which givesNote the additional factor of 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 and momenta in a general container of volume V described by a cell matrix h. In dimensions h would describe a generally triclinic cell. Let denote the interparticle potential and 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 , 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 and . 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 . The oxygen and the hydrogen have point charges and 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 and . The stability of the algorithms is measured by the conserved quantity , where i runs over the number of configurations , and 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)
- et al.
J. Comput. Phys.
(1999) - et al.
J. Comput. Phys.
(1999) Ann. Phys.
(2005)Phys. Lett. A
(1990)- et al.
J. Comput. Phys.
(1977) J. Comput. Phys.
(1983)- et al.
Comput. Phys. Commun.
(2000) Phys. Rev. A
(1985)- et al.
J. Chem. Phys.
(1992) - et al.
Mol. Phys.
(1996)
J. Chem. Phys.
Phys. Rev. A
J. Chem. Phys.
J. Chem. Phys.
J. Chem. Phys.
J. Chem. Phys.
J. Chem. Phys.
Europhys. Lett.
J. Chem. Phys.
J. Chem. Phys.
J. Phys. A: Math. Gen.
J. Phys. Chem. B
Cited by (52)
Structure of amorphous BaTiO<inf>3</inf> by molecular dynamics simulations using a shell model
2020, Physica B: Condensed MatterCitation Excerpt :The pressure was controlled by the Parrinello-Rahman method [16]. The temperature was controlled by the massive Nose-Hoover chain method [17–19]. The externally applied pressure was set to 0 Pa throughout the simulations.
Fullerenes containing water molecules: a study of reactive molecular dynamics simulations
2023, Physical Chemistry Chemical PhysicsStatistical mechanics: Theory and molecular simulation
2023, Statistical Mechanics: Theory and Molecular Simulation