(372e) A Hamiltonian-Based Algorithm for Rigorous Molecular Dynamics Simulation in the Nve, Nvt, Npt, and Nph Ensembles

When performing molecular dynamics simulations one attempts to model accurately a large system using a relatively small (at this time, generally less than a million) number of molecules.  It is of course the hope that the simulation results will directly correspond to experiments on a similar, but macroscopic, system.  One way to verify the validity of molecular-level simulations is to show from a theoretical point of view that the trajectories, i.e. the positions and momenta of the particles as a function of time, generated by a particular algorithm, rigorously correspond to the appropriate statistical mechanical ensemble.

For molecular dynamics (MD) simulations in the microcanonical (NVE) ensemble, a well-accepted and straightforward algorithm exists.  In the canonical (NVT) ensemble, a rigorous MD algorithm known as the Nosé-Hoover thermostat exists, which introduces an extended system containing a thermostat [1, 2].  While it is well known that the Nosé-Hoover thermostat rigorously generates trajectories in the NVT ensemble in the absence of external forces, there is still much discussion as to rigor of the algorithm in the presence of external forces.  In the isobaric-isothermal ensemble (NpT), the story is not so clear.  There are numerous NpT algorithms in the literature with varying degrees of rigor [2-10].

One of the main differences in the history of the development of NVT and NpT algorithms is that the derivation of the well-accepted and rigorous NVT algorithm (the Nosé-Hoover thermostat) started with a Hamiltonian of Nosé.  Relying on the canonical (i.e. symplectic) nature of the Hamiltonian, the equations of motion are unambiguously determined.  Nosé proved this algorithm rigorously generated trajectories in the NVT ensemble.  Hoover then added his contribution by providing a non-canonical transformation from the mathematical frame of reference to a more physical frame of reference in which analysis of the MD results is straightforward. 

In 1984, Nosé provided an NpT Hamiltonian, but did not express the equations of motion in terms of the physical variables [3]. However, for whatever reason, from there the development of NpT algorithms took a very different path.  The modern NpT algorithms do not possess a Hamiltonian.  The starting point of their derivation is simply a presumed form of the equations of motion.  Recently, Tuckerman et al. introduced a methodical procedure based in the statistical mechanics of non-Hamiltonian systems that allows one to determine the rigor of various published NpT algorithms [11, 12].  They found that the NpT algorithm of Martyna et al. [5] satisfied their criteria for rigor, while the NpT algorithms of Hoover [2] and that of Melchionna et al. [4] were not rigorous.  Curiously they did not test Nosé's NpT algorithm [3], which predates all of the other algorithms.

This paper has two purposes.  The first purpose is to demonstrate that following the same methodical procedure that was used by Nosé and Hoover to develop a rigorous algorithm for MD simulation in the NVT ensemble, one can also generate a rigorous algorithm for MD simulation in the NpT ensemble.  This process starts with a Hamiltonian, namely Nosé's 1984 NpT Hamiltonian [3], and includes a non-canonical transformation from the mathematical system to the physically meaningful system.  We prove that the algorithm is rigorous in the NpT ensemble.  To our knowledge, our NpT algorithm is different than any other NpT algorithm published to date.

The second purpose of this paper is to address the problem associated with simulation in the NVT, NpT, and NpH ensembles in the presence of external forces (or when the total linear momentum is non-zero).  We present a procedure for a systematic Hamiltonian-based development of rigorous MD algorithms in the presence or absence of external forces.  We then use this procedure to develop rigorous algorithms for the NVE, NVT, NpT, and NpH ensembles.  The resulting algorithms are new.

Results and Discussion

Our family of algorithms for MD simulation in the NVE, NVT, NpT, and NpH ensembles is completely general.  Each is valid in the absence or presence of external forces.  Each is valid regardless of the conservation of total linear momentum.  Our family of algorithms is also completely rigorous.  Following Nosé, using the statistical mechanics of Hamiltonian systems, our algorithms rigorously generate trajectories in the appropriate ensemble.  We also show that, following Tuckerman, and using the statistical mechanics of non-Hamiltonian systems, our algorithms rigorously generate trajectories in the appropriate ensemble. 

Our family of algorithms is completely self-consistent.  By self-consistent we mean that the NpT algorithm is the most general algorithm.  The other algorithms are subsets of the NpT case; we have an NVT algorithm when the barostat is turned off, an NpH algorithm when the thermostat is turned off, and an NVE algorithm when both the thermostat and barostat are inactive.

Our family of algorithms is consistent with previously existing rigorous algorithms.  For example, in the absence of external forces and when the total momentum is initialized to zero, our NVT algorithm reduces to the Nosé-Hoover thermostat.  Thus, we identify the two limiting constraints in order for the Nosé-Hoover thermostat to be valid.  However, our generalized algorithms are rigorous regardless of whether these constraints are satisfied.  As another example, our NpT algorithm reduces to the same equations of motion as Nosé's NpT algorithm, again under a set of constraints.

In short, using a Hamiltonian-based approach, we have derived rigorous algorithms for MD simulation in the NVE, NVT, NpT, and NpH in the presence or absence of external forces.  Moreover, we have provided a methodical procedure to derive these equations, which can also be used to derive nonequilibrium MD algorithms.

Finally, because our NpT algorithm is different than the NpT algorithm of Martyna et al., we discuss the validity of the two algorithms.  Our NpT algorithm satisfies both the Hamiltonian and non-Hamiltonian criteria for rigor.  The NpT algorithm of Martyna et al. satisfies only the non-Hamiltonian criteria for rigor.  We discuss the uniqueness of extended algorithms.

References

1.  Nosé, S., A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys., 1984. 52(2): p. 255-268.

2.  Hoover, W.G., Canonical dynamics:  Equilibrium phase-space distributions. Phys. Rev. A, 1985. 31(3): p. 1695-1697.

3.  Nosé, S., A Unified Formulation of the Constant Temperature Molecular-Dynamics Methods. J. Chem. Phys., 1984. 81: p. 511-519.

4.  Melchionna, S., G. Ciccotti, and B.L. Holian, Hoover NPT dynamics for systems varying in size and shape. Mol. Phys., 1993. 78(3): p. 533-544.

5.  Martyna, G.J., D.J. Tobias, and M.L. Klein, Constant pressure molecular dynamics algorithms. J. Chem. Phys., 1994. 101(5): p. 4177-4189.

6.  Andersen, H.C., Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys., 1980. 72(4): p. 2384-2393.

7.  Kalibaeva, G., M. Ferrario, and G. Ciccotti, Constant pressure-constant temperature molecular dynamics:  a correct NPT ensemble using the molecular virial. Mol. Phys., 2003. 101(6): p. 765-778.

8.  Melchionna, S., Constrained systems and statistical distribution. Phys. Rev. E, 2000. 61(6): p. 6165-6170.

9.  Melchionna, S., Erratum: Constrained systems and statistical distribution. Phys. Rev. E, 2000. 62(4): p. 5864.

10.  Parrinello, M. and A. Rahman, Crystal Structure and Pair Potentials:  A Molecular Dynamics Study. Phys. Rev. Lett., 1980. 45(14): p. 1196-1199.

11.  Tuckerman, M.E., C.J. Mundy, and G.J. Martyna, On the classical statistical mechanics of non-Hamiltonian systems. Europhys. Lett., 1999. 45(2): p. 149-155.

12.  Tuckerman, M.E., Y. Liu, G. Ciccotti, and G.J. Martyna, Non-Hamiltonian molecular dynamics: Generalizing Hamiltonian phase space principles to non-Hamiltonian systems. J. Chem. Phys., 2001. 115(4): p. 1678-1702.