# Symplectic algorithm for constant-pressure molecular dynamics using a Nosé-Poincaré thermostat

## 1 Introduction

Traditionally, molecular-dynamics simulations are performed using constant particle number $`N`$, volume $`V`$ and energy $`E`$. However, these are not usually the conditions under which experiments are done and there has been much attention to the development of simulation methods designed to sample from other, experimentally more relevant ensembles, such as constant temperature (canonical) and/or constant pressure. Some of the most popular and useful of these are those based on so-called “extended” Hamiltonians, i.e., Hamiltonians in which extra degrees of freedom have been added to the system in order to ensure that the trajectory samples from the statistical distribution corresponding to the desired thermodynamic conditions.

For a constant pressure system, for example, Andersen introduced the volume $`V`$, along with its corresponding conjugate momentum $`\pi _V`$, as extra variables. The new variables are coupled to the system in such a way as to guarantee that the trajectory (if ergodic) samples from an isobaric statistical distribution. Similarly, to generate a constant temperature distribution, Nosé introduced a new mechanical variable $`s`$ (with conjugate momentum $`\pi _s`$) that couples into the system through the particle momenta and acts to effectively rescale time in such a way as to guarantee canonically distributed configurations. These two extensions can be combined to give a Hamiltonian whose trajectories can be shown to sample from an isothermal-isobaric ensemble.

This combined Nosé-Andersen (NA) Hamiltonian is given by

$$_{NA}=V^{2/3}\frac{p_i^2}{2m_is^2}+U(V^{1/3}𝐪)+\frac{\pi _V^2}{2Q_V}+\frac{\pi _s^2}{2Q_s}+gkT\mathrm{ln}s+P_{ext}V,$$
(1)   
where $`p_i`$ is the conjugate momentum to the scaled position $`q_i=V^{1/3}r_i`$, $`P_{ext}`$ is the external pressure and $`g`$ is given by $`N_f+1`$ where $`N_f`$ is number of degrees of freedom of the original system. The quantities $`Q_V`$ and $`Q_s`$ are the masses of the Andersen “piston” and the Nosé thermostat variable, respectively.

The equations of motion for this system are

$`\dot{p_i}`$ $`=`$ $`V^{1/3}_iU(V^{1/3}𝐪)`$ (2a)   
$`\dot{q_i}`$ $`=`$ $`{\displaystyle \frac{p_i}{s^2m_iV^{2/3}}}`$ (2b)   
$`\dot{\pi _V}`$ $`=`$ $`𝒫P_{ext}`$ (2c)   
$`\dot{V}`$ $`=`$ $`\pi _v/Q_V`$ (2d)   
$`\dot{\pi _s}`$ $`=`$ $`V^{2/3}s^3{\displaystyle \frac{p_i^2}{m_i}}{\displaystyle \frac{gkT}{s}}`$ (2e)   
$`\dot{s}`$ $`=`$ $`\pi _s/Q_s,`$ (2f)   
where the instantaneous pressure $`𝒫`$ is given by

$$𝒫=\frac{2}{3V}\underset{i}{}\frac{p_i^2}{2m_iV^{2/3}s^2}\frac{1}{3V}\underset{i}{}\frac{U}{q_i}q_i$$
(3)   
There are two major drawbacks to this approach: First, because of the time rescaling, the time variable in Nosé dynamics is not “real” time, so any discretized trajectory generated by numerically integrating the Nosé equations of motion must be transformed back into “real” time, leading to the configurations that are spaced at unequal “real”-time intervals. This is inconvenient for the construction of equilibrium averages, especially of dynamical quantities. Second, the Hamiltonian is not separable (that is, the kinetic and potential terms in the Hamiltonian are not functions only of momenta and position variables, respectively), making standard Verlet/leapfrog approaches inapplicable.

By a change of variables and a time rescaling of the equations of motion, Hoover derived new equations of motion that generate the same trajectories (for the exact solution) as the original Nosé Hamiltonian, but in real time. This Nosé-Hoover dynamics has become a standard method in molecular simulation. However, the change of variables that links the Nosé Hamiltonian to the Nosé-Hoover equations of motion is a non-canonical transformation - the total energy function of the system is still conserved, but it is no longer a Hamiltonian, since the equations of motion cannot be derived from it. Although a variety of very good time-reversible methods have been put forward, the lack of Hamiltonian structure precludes the use of symplectic integration schemes, which have been shown to have superior stability over non-symplectic methods.

## 2 The Nosé-Poincaré-Andersen (NPA) Hamiltonian

Recently, Bond, Leimkuhler and Laird have developed a new formulation of Nosé constant-temperature dynamics in which a Poincaré time transformation is applied directly to the Nosé Hamiltonian, instead of applying a time transformation to the equations of motion as in Nosé-Hoover. The result of this is a method that runs in “real” time, but is also Hamiltonian in structure. In this work we combine this new thermostat with the Andersen method for constant pressure to give an algorithm for isothermal-isobaric molecular dynamics. For a system with an Andersen piston, the new Nosé-Poincaré-Andersen (NPA) Hamiltonian is given by

$$_{NPA}=[_{NA}_{NA}(t=0)]s,$$
(4)   
where $`_{NA}`$ is given in Eq. 1. As discussed in Ref. 12, the above form of the Hamiltonian (a specific case of a Poincaré time transformation) will generate the same trajectories as the original Nosé-Andersen Hamiltonian, except with time rescaled by $`s`$ (which puts the trajectories back into real time), The resulting equations of motion (except for $`\pi _s`$) for this constant pressure and temperature Nosé-Poincaré Hamiltonian are the as those given above for the Nosé-Andersen system (Eqs. 2a-2f), except that the right-hand side is multiplied by the thermostat variable $`s`$. For $`\pi _s`$ we have

$$\dot{\pi _s}=s\frac{}{s}\mathrm{\Delta }H=V^{2/3}\frac{p_i^2}{m_is^2}gkT\mathrm{\Delta }H$$
(5)   
where $`\mathrm{\Delta }_{NA}_{NA}(t=0)`$.

## 3 Integrating the NPA Equations of motion

The NPA Hamiltonian is nonseperable since the kinetic energy contains the extended “position” variables $`s`$ and $`V`$. The equations of motion for a general time-independent, non-seperable Hamiltonian can be written (for general positions $`Q`$ and conjugate momenta $`P`$)

$`\dot{Q}`$ $`=`$ $`G(P,Q)`$   
$`\dot{P}`$ $`=`$ $`F(P,Q),`$ (6)   
where $`G(P,Q)=\frac{}{P}`$ and $`F(P,Q)=\frac{}{Q}`$. (For a seperable Hamiltonian $`G`$ is only a function of $`P`$ and $`F`$ is only a function of $`Q`$.) For such a nonseperable system, standard symplectic splitting methods, such as the Verlet/leapfrog algorithm, are not directly applicable. However, symplectic methods specifically for nonseperable systems have been developed. One simple example that is second-order and time-reversible is the Generalized Leapfrog Algorithm (GLA)

$`P_{n+1/2}`$ $`=`$ $`P_n+hF(P_{n+1/2},Q_n)/2`$   
$`Q_{n+1}`$ $`=`$ $`Q_n+h[G(P_{n+1/2},Q_n)+G(P_{n+1/2},Q_{n+1})]/2`$   
$`P_{n+1}`$ $`=`$ $`P_{n+1/2}+hF(P_{n+1/2},Q_{n+1})/2,`$ (7)   
where $`h`$ is the time step and $`P_n`$ and $`Q_n`$ are the approximations to $`P(t)`$ and $`Q(t)`$ at $`t=t_n=nh`$. (This method can be obtained as the concatenation of the Symplectic Euler method,

$`P_{n+1}`$ $`=`$ $`P_n+hF(P_{n+1},Q_n)`$   
$`Q_{n+1}`$ $`=`$ $`Q_n+hG(P_{n+1},Q_n),`$ (8)   
with its adjoint. The concatenation of a integrator with its adjoint guarantees a time-reversible method.) This method is a simple example of a class of symplectic integrators for nonseperable Hamiltonians.

Applying the GLA to the NPA equations of motion gives

$`p_{i,n+1/2}`$ $`=`$ $`p_{i,n}{\displaystyle \frac{h}{2}}s_nV_n^{1/3}_iU(V^{1/3}𝐪)`$ (9a)   
$`\pi _{v,n+1/2}`$ $`=`$ $`\pi _{v,n}+{\displaystyle \frac{h}{2}}s_n[𝒫(𝐪_n,𝐩_{n+1/2},V_n,s_n)P_{ext}]`$ (9b)   
$`\pi _{s,n+1/2}`$ $`=`$ $`\pi _{s,n}+{\displaystyle \frac{h}{2}}\left({\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{p_{i,n+1/2}^2}{m_iV_n^{2/3}s_n^2}}gk_BT\right)`$ (9c)   
$`{\displaystyle \frac{h}{2}}\mathrm{\Delta }(𝐪_n,𝐩_{n+1/2},V_n,\pi _{v,n+1/2},s_n,\pi _{s,n+1/2})`$   
$`s_{n+1}`$ $`=`$ $`s_n+{\displaystyle \frac{h}{2}}(s_n+s_{n+1}){\displaystyle \frac{\pi _{s,n+1/2}}{Q_s}}`$ (9d)   
$`V_{n+1}`$ $`=`$ $`V_n+{\displaystyle \frac{h}{2}}(s_n+s_{n+1}){\displaystyle \frac{\pi _{v,n+1/2}}{Q_v}}`$ (9e)   
$`q_{i,n+1}`$ $`=`$ $`q_{i,n}+{\displaystyle \frac{h}{2}}\left({\displaystyle \frac{1}{s_nV_n^{2/3}}}+{\displaystyle \frac{1}{s_{n+1}V_{n+1}^{2/3}}}\right){\displaystyle \frac{p_{i,n+1/2}}{m_i}}`$ (9f)   
$`\pi _{s,n+1}`$ $`=`$ $`\pi _{s,n+1/2}+{\displaystyle \frac{h}{2}}\left({\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{p_{i,n+1/2}^2}{m_iV_{n+1}^{2/3}s_{n+1}^2}}gk_BT\right)`$ (9g)   
$`{\displaystyle \frac{h}{2}}\mathrm{\Delta }(𝐪_{n+1},𝐩_{n+1/2},V_{n+1},\pi _{v,n+1/2},s_{n+1},\pi _{s,n+1/2})`$   
$`\pi _{v,n+1}`$ $`=`$ $`\pi _{v,n+1/2}+{\displaystyle \frac{h}{2}}s_{n+1}[𝒫(𝐪_{n+1},𝐩_{n+1/2},V_{n+1},s_{n+1})P_{ext}]`$ (9h)   
$`p_{i,n+1}`$ $`=`$ $`p_{i,n+1/2}+{\displaystyle \frac{h}{2}}s_{n+1}V_{n+1}^{1/3}_iU(V_{n+1}^{1/3}𝐪_{𝐧+\mathrm{𝟏}})`$ (9i)   
As in the case of the constant volume Nosé-Poincare algorithm, the GLA for the NPA is explicit - this is not necessarily the case for a general nonseperable Hamiltonian. Note that Eq. 9c requires the solution of a scalar quadratic equation for $`\pi _{s,n+1/2}`$. Details of how to solve this equation without involving subtractive cancellation can be found in Ref. 12.

## 4 Test Simulation Results and Summary

In order to evaluate this method, simulations were performed using an embedded atom potential for aluminum. We report test simulations on a system of 256 particles with periodic boundary conditions for an aluminum melt at $`T=1000K`$ and $`P=0`$. For this model mass is measured in amu, distance in Å and energy in eV, the natural time unit of the simulation is then 10.181fs, that is, a simulation time step of 0.1 corresponds to an acutal time step of 1.0181fs.

First the stability of the method was tested. Fig. 1(a) shows the value of the NPA Hamiltonian (a conserved quantity) as a function of time in a long run. The trajectory shown here was begun after initial equilibration at 1000K for 2$`\times 10^6`$ time steps (2.03ns). In this simulation, the values of $`Q_v`$ and $`Q_s`$ (in reduced units) were 10<sup>-4</sup> and 2.5, respectively. The stability of the method is excellent, giving no noticeable drift in $`_{\text{NPA}}`$ over the course of a long trajectory. The pressure and temperature trajectories for this run are also shown in Figs. 1(b) and 1(c), respectively.

In Fig. 2(a)-(c), we show the ability of the method to regulate the pressure, temperature and density, respectively, for three sets of extended variable masses. (The values for $`Q_s`$ were larger here than that used in Figure 1 because small values of that variable lead to instabilities when the initial temperature is far from the target temperature.) The system was initialized to an fcc lattice of initial density 0.06021Å<sup>-3</sup> with the individual velocity components chosen from a Maxwell-Boltzmann distribution at 100K. The simulation was then run with the NPA with $`T=1000K`$ and $`P=0`$. In all cases, the instantaneous pressure, temperature and density evolve quickly and stabilize about their desired values.

The GLA has a global error that is second-order in the time step. To demonstrate that this also is true in our results, a series of simulations were performed using various values for the time step. The system was initialized in an identical manner to that described in the last paragraph and then run for a total time of 2.0362 ps. Figure  shows, for several combinations of extended variable masses, a log-log plot of the energy error, as estimated by the standard deviation of $`_{\text{NPA}}`$ (Eq. 3), versus the time step. One notes that here smaller values of $`Q_v`$ lead to smaller fluctuations in $`_{\text{NPA}}`$.

Finally, to demonstrate that the method yields relevant dynamical quantities, the normalized velocity autocorrelation function, $`C(t)=𝐯(t)𝐯(0)/𝐯(0)𝐯(0)`$, was calculated using our constant NPT algorithm (with $`Q_v`$ and $`Q_s`$ as in Fig. 1) and compared to the same quantity calculated using standard constant NVE molecular dynamics (with a velocity-Verlet integrator). The NVE simulations were run at an energy and density corresponding to the average energy and density for the constant NPT simulations. This comparison is shown in Figure . Both systems were first equilibrated at 1000K for 200,000 steps (203.6 ps) and run for 20,000 steps (20.36 ps) to collect averages. $`C(t)`$ for the Nosé-Poincaré-Andersen method for constant NPT molecular dynamics is indistinguishable in this figure from that of the NVE simulation.

## 5 Acknowledgements

We gratefully acknowledge Steve Bond and Ben Leimkuhler for helpful conversations and invaluable advice, as well as the Kansas Center for Advanced Scientific Computing for the use of their computer facilities. We also would like to thank the National Science Foundation (under grants CHE-9500211 and DMS-9627330) as well as the University of Kansas General Research Fund for their generous support of this research.