# Optimized Constant Pressure Stochastic Dynamics

## I Introduction

Molecular Dynamics (MD) simulations  are a very efficient tool to study the statistical properties of thermodynamic systems, especially at high densities where the acceptance rates of standard Monte Carlo (MC) simulations  are small. They are also very well suited to study dynamical properties.

MD simulations are most naturally performed within the microcanonical (NVE) ensemble, while the simple Metropolis MC algorithm leads to the canonical (NVT) ensemble. However, it is often desirable to study a system within a different ensemble, and hence methods have been developed to extend both MC and MD methods to practically every ensemble of thermodynamics . For MC methods, this task is relatively straightforward: One starts from the partition function or equilibrium probability distribution function for the pertinent degrees of freedom, and defines a Markov process in the space of the latter. After verifying the condition of detailed balance, and making use of ergodicity (which often can be safely assumed, and in some cases even be shown rigorously), one has established that the process will ultimately produce fluctuations within the proper equilibrium. For example, for a simulation in the isothermal–isobaric (NPT) ensemble, the system is studied in a box with periodic boundary conditions, whose size is allowed to fluctuate. In order to keep the system homogeneous, the particle coordinates are instantaneously adjusted to these box fluctuations, via simple rescaling. We shall here be concerned with the simplest case only, where the box is just a cube of size $`L`$ in each direction. By writing

$$\stackrel{}{r}_i=L\stackrel{}{s}_i$$
(1)   
one introduces reduced coordinates $`\stackrel{}{s}_i`$ in the unit cube instead of the original coordinates $`\stackrel{}{r}_i`$ of the $`N`$ particles in the box of volume $`V=L^d`$, $`d`$ denoting the spatial dimensionality. The abovementioned adjustment of the particle configuration to the box fluctuations is facilitated by updating $`L`$ and the $`\stackrel{}{s}_i`$ independently. The partition function appropriate to the NPT ensemble is then

$`Z`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑V{\displaystyle d^d\stackrel{}{r}_1\mathrm{}d^d\stackrel{}{r}_N\mathrm{exp}\left(\beta U\beta PV\right)}`$ (2)   
$`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑V{\displaystyle d^d\stackrel{}{s}_1\mathrm{}d^d\stackrel{}{s}_NV^N\mathrm{exp}\left(\beta U\beta PV\right)},`$ (3)   
where $`U`$ denotes the system’s potential energy, $`P`$ is the external pressure, and $`\beta =1/(k_BT)`$. From this one immediately reads off that one has to run a standard Metropolis algorithm on the state space $`(L,\{\stackrel{}{s}_i\})`$, using an effective Hamiltonian

$$U_{eff}=U+PVNk_BT\mathrm{ln}V.$$
(4)   
The MD approach  to non–microcanonical ensembles , pioneered by Andersen , Nosé  and Hoover , is slightly more involved. Like in MC, one defines an additional dynamical variable whose fluctuations allow to keep the thermodynamically conjugate variable fixed. In our example, this variable is $`L`$, while $`P`$ is fixed. However, the dynamics is not specified via a random process, but rather by Hamiltonian equations of motion in the extended space. This requires the definition of canonically conjugate momenta of the additional variables, and in turn the introduction of artificial masses, which have no direct physical meaning but are adjusted in order to set the time scale of the new fluctuations. The analysis then proceeds by assuming ergodicity in the extended space, such that the statistics is described by a microcanonical ensemble in that space. If the equations of motion have been chosen properly, then the equilibrium distribution function of the algorithm, which is obtained by integrating out the artificial momenta, must coincide with that of the desired ensemble. For constant–pressure simulations, often the picture of a “piston” is employed. In the last few years many contributions on extending constant–pressure methods have been made, treating cases of isotropic boxes as well as non–cubic cells . We will here however be only concerned with the case of an isotropic box, as in the original Andersen method  which produces an isobaric–isenthalpic (NPH) ensemble.

Stochastic dynamics (SD) in its classical form  can be viewed as a simulation method which is somewhere between MC and MD, sharing with the former the stochastic element (which also ensures the ergodicity of the method) and the generation of a canonical (NVT) ensemble, while being based on continuous equations of motion, and momenta, like the latter. Instead of Hamiltonian equations of motion one solves Langevin equations and adjusts the temperature via the balance between friction force and stochastic force (fluctuation–dissipation theorem). It is not surprising that this approach can be combined with the Andersen method in order to produce an NPT ensemble. It is the purpose of this paper to show that the algorithm can be derived in a very simple and straightforward manner, by exploiting the well–known equivalence of the Langevin equation with the Fokker–Planck equation , and avoiding the complicated reasoning of a recent publication . The derivation is outlined in Secs. II, III and IV, where we treat a straightforward generalization of SD to arbitrary Hamiltonian systems, and apply this to the original Andersen NPH method, which is thus modified to an NPT method. We then turn to the question of numerical implementation. Since SD reduces to standard Hamiltonian dynamics in the limit of zero friction, and practical applications are often run for rather small friction, it is useful to use an algorithm which reduces to a time–reversible symplectic integrator (TRSI) in the zero–friction limit. It is well–known that TRSIs are particularly well–suited for Hamiltonian systems , since, except for roundoff errors, they do not prefer a particular direction of time and are hence very stable. Our algorithm is derived and tested in Sec. V, where we apply the method to a simple model system of Lennard–Jones particles, and pay particular attention to the question how the parameters should be chosen for optimum performance. Finally, we conclude in Sec. VI.

## II Generalized Stochastic Dynamics

Our starting point is a set of generalized coordinates $`q_i`$ and canonically conjugate momenta $`p_i`$, such that the Hamiltonian equations of motion read

$$\dot{q}_i=\frac{}{p_i}\dot{p}_i=\frac{}{q_i},$$
(5)   
where $``$ is the Hamilton function. These equations of motion are generalized to their stochastic versions by adding friction forces with damping parameters $`\gamma _i\left(\left\{q_i\right\}\right)`$ and stochastic forces with noise strengths $`\sigma _i\left(\left\{q_i\right\}\right)`$ (note that a dependence on the generalized coordinates is permitted):

$$\dot{q}_i=\frac{}{p_i}\dot{p}_i=\frac{}{q_i}\gamma _i\frac{}{p_i}+\sigma _i\eta _i(t),$$
(6)   
where the Gaussian white noise $`\eta _i`$ satisfies the usual relations

$$\eta _i(t)=0\eta _i(t)\eta _j(t^{})=2\delta _{ij}\delta \left(tt^{}\right).$$
(7)   
Equivalently, the stochastic process is described by the Fokker–Planck equation  which describes the time evolution of the probability distribution function $`\mathrm{\Phi }`$ in the full space of stochastic variables. For an arbitrary set of variables $`x_\nu `$ (in our case both $`q_i`$ and $`p_i`$) it reads

$$\frac{\mathrm{\Phi }}{t}=i\widehat{L}_{FP}\mathrm{\Phi }=\underset{\nu }{}\frac{}{x_\nu }D_\nu ^{(1)}\mathrm{\Phi }+\underset{\mu \nu }{}\frac{^2}{x_\mu x_\nu }D_{\mu \nu }^{(2)}\mathrm{\Phi },$$
(8)   
where the right hand side of the equation defines the Fokker–Planck operator $`\widehat{L}_{FP}`$. Drift and diffusion coefficients $`D_\nu ^{(1)}`$ and $`D_{\mu \nu }^{(2)}`$ are related to the short–time behavior of the process via the Kramers–Moyal expansion :

$`D_\nu ^{(1)}`$ $`=`$ $`\underset{\mathrm{\Delta }t0}{lim}{\displaystyle \frac{1}{\mathrm{\Delta }t}}\mathrm{\Delta }x_\nu (\mathrm{\Delta }t)=\underset{\mathrm{\Delta }t0}{lim}{\displaystyle \frac{1}{\mathrm{\Delta }t}}x_\nu (t+\mathrm{\Delta }t)x_\nu (t)`$ (9)   
$`D_{\mu \nu }^{(2)}`$ $`=`$ $`\underset{\mathrm{\Delta }t0}{lim}{\displaystyle \frac{1}{2\mathrm{\Delta }t}}\mathrm{\Delta }x_\mu (\mathrm{\Delta }t)\mathrm{\Delta }x_\nu (\mathrm{\Delta }t).`$ (10)   
Straightforward evaluation of these moments for the present case yields directly the Fokker–Planck operator:

$`i\widehat{L}_{FP}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{}{q_i}}{\displaystyle \frac{}{p_i}}{\displaystyle \underset{i}{}}{\displaystyle \frac{}{p_i}}\left({\displaystyle \frac{}{q_i}}\gamma _i{\displaystyle \frac{}{p_i}}\right)+{\displaystyle \underset{i}{}}{\displaystyle \frac{^2}{p_i^2}}\sigma _i^2`$ (11)   
$`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{}{p_i}}{\displaystyle \frac{}{q_i}}+{\displaystyle \underset{i}{}}{\displaystyle \frac{}{q_i}}{\displaystyle \frac{}{p_i}}+{\displaystyle \underset{i}{}}{\displaystyle \frac{}{p_i}}\left(\gamma _i{\displaystyle \frac{}{p_i}}+\sigma _i^2{\displaystyle \frac{}{p_i}}\right).`$ (12)   
A canonical ensemble is generated if the Boltzmann distribution is the stationary distribution of the process (due to the stochastic element, the process is usually ergodic, such that only one stationary distribution exists). For the present case one finds

$$i\widehat{L}_{FP}\mathrm{exp}\left(\beta \right)=\underset{i}{}\frac{}{p_i}\left(\gamma _i\beta \sigma _i^2\right)\frac{}{p_i}\mathrm{exp}\left(\beta \right),$$
(13)   
which vanishes if

$$\sigma _i^2=k_BT\gamma _i.$$
(14)   
This relation is the generalized fluctuation–dissipation theorem which friction and noise have to satisfy in order to generate a canonical ensemble.

## III Andersen Extended System

The Andersen method  uses the box volume $`V=L^d`$ and the scaled coordinates $`\stackrel{}{s}_i`$, see Eqn. 1, as degrees of freedom. For the particle velocities one obtains

$$\dot{\stackrel{}{r}_i}=L\dot{\stackrel{}{s}_i}+\dot{L}\stackrel{}{s}_i,$$
(15)   
however, the second term is deliberately omitted in order to achieve independent fluctuations of $`L`$ and $`\stackrel{}{s}_i`$. Hence the method amounts to postulating the Lagrangian

$``$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{L^2}{2}}m_i\dot{\stackrel{}{s}_i}^2{\displaystyle \underset{i<j}{}}v_{ij}(L,\{\stackrel{}{s}_i\})+{\displaystyle \frac{Q}{2}}\dot{V}^2PV,`$ (16)   
where $`m_i`$ denotes the mass of the $`i`$th particle, $`Q`$ is the artificial piston mass or box mass, and $`v_{ij}`$ is the interaction potential between particles $`i`$ and $`j`$ (the generalization to three– and many–body forces is straightforward). The Hamiltonian is then obtained via Legendre transformation

$$=\underset{i}{}\stackrel{}{\pi }_i\dot{\stackrel{}{s}_i}+\mathrm{\Pi }_V\dot{V}=\underset{i}{}\frac{1}{2L^2m_i}\stackrel{}{\pi }_i^2+\underset{i<j}{}v_{ij}+\frac{1}{2Q}\mathrm{\Pi }_V^2+PV,$$
(17)   
where we have used the canonically conjugate momenta

$$\mathrm{\Pi }_V=\frac{}{\dot{V}}=Q\dot{V},\stackrel{}{\pi }_i=\frac{}{\dot{\stackrel{}{s}_i}}=m_iL^2\dot{\stackrel{}{s}_i},$$
(18)   
such that the Hamiltonian equations of motion read

$`\dot{\stackrel{}{s}_i}={\displaystyle \frac{1}{L^2m_i}}\stackrel{}{\pi }_i`$ $`\dot{\stackrel{}{\pi }_i}=L\stackrel{}{f}_i`$ (19)   
$`\dot{V}={\displaystyle \frac{1}{Q}}\mathrm{\Pi }_V`$ $`\dot{\mathrm{\Pi }}_V=𝒫P,`$ (20)   
where $`\stackrel{}{f}_i`$ is the force acting on the $`i`$th particle, and $`𝒫`$ abbreviates the “instantaneous” pressure

$$𝒫=\frac{L}{dV}\underset{i<j}{}\stackrel{}{f}_{ij}\stackrel{}{s}_{ij}+\frac{1}{dL^2V}\underset{i}{}\frac{1}{m_i}\stackrel{}{\pi }_i^2,$$
(21)   
$`\stackrel{}{f}_{ij}`$ being the force acting between particle $`i`$ and $`j`$, while $`\stackrel{}{s}_{ij}=\stackrel{}{s}_i\stackrel{}{s}_j`$.

Obviously, $``$ is a constant of motion. Apart from the kinetic energy of the piston, and the deviation of the simulated molecular kinetic energy from the true kinetic energy, which are both small corrections , $``$ is just the enthalpy. For this reason, the method produces the NPH ensemble.

## IV Langevin Equation for the NPT Ensemble

The idea of Feller et al.  was to replace the canonical equations of motion (Eqn. 19) by a Langevin stochastic process. It was designed to avoid oscillations of the box volume (“ringing” of the box). In Ref.  an infinite set of harmonic oscillators coupled to the box piston was used to prove the correctness of the approach. However, since the original Andersen method is based on a Hamiltonian system, the results which have been derived in Sec. II can be directly applied. Therefore the stochastic equations of motion read

$`\dot{\stackrel{}{\pi }_i}`$ $`=`$ $`L\stackrel{}{f}_i{\displaystyle \frac{\gamma _i}{L^2m_i}}\stackrel{}{\pi }_i+\sqrt{k_BT\gamma _i}\stackrel{}{\eta }_i(t)`$ (22)   
$`\dot{\mathrm{\Pi }}_V`$ $`=`$ $`𝒫P{\displaystyle \frac{\gamma _V}{Q}}\mathrm{\Pi }_V+\sqrt{k_BT\gamma _V}\stackrel{}{\eta }_V(t)`$ (23)   
$`\eta _i^\alpha =\eta _V`$ $`=`$ $`0`$ (24)   
$`\eta _i^\alpha (t)\eta _j^\beta (t^{})`$ $`=`$ $`2\delta _{ij}\delta _{\alpha \beta }\delta (tt^{})`$ (25)   
$`\eta _V(t)\eta _V(t^{})`$ $`=`$ $`2\delta (tt^{})`$ (26)   
$`\eta _i^\alpha (t)\eta _V(t^{})`$ $`=`$ $`0,`$ (27)   
while the equations of motion for $`\stackrel{}{s}_i`$ and $`V`$ remain unchanged. Here $`\alpha `$ and $`\beta `$ denote Cartesian directions, and the fluctuation–dissipation relation has already been taken into account.

This method generates the correct NPT ensemble, as is seem from writing down the partition function which naturally arises from the algorithm (cf. Sec. II):

$$Z=𝑑\mathrm{\Pi }_Vd^d\stackrel{}{\pi }_1\mathrm{}d^d\stackrel{}{\pi }_N𝑑Vd^d\stackrel{}{s}_1\mathrm{}d^d\stackrel{}{s}_N\mathrm{exp}\left(\beta \right),$$
(28)   
where one obviously has to use the Hamiltonian of the Andersen method, see Eqn. 17. Integrating out the momenta, one sees directly that this is, apart from unimportant prefactors, identical to Eqn. 2, i. e. the correct partition function of the NPT ensemble.

Note that there is still considerable freedom in the choice of the friction parameters. Feller et al. chose $`\gamma _i=0`$ and $`\gamma _V=\text{const.}`$. This is tantamount to coupling only the piston degree of freedom to the heat bath. Since this degree of freedom is tightly coupled to all the others, the method produces the same NPT ensemble as the more general case $`\gamma _i0`$. However, we view this latter case, which also includes a direct coupling of the particles to the heat bath, as more advantageous, not for fundamental reasons, but rather for practical ones: As in standard SD, every degree of freedom is thermostatted individually. For this reason, local instabilities, arising from discretization errors, are efficiently corrected for, without spreading throughout the system. Loosely spoken, a particle which happens to be too “hot” will be “cooled down” by its local friction, while in the opposite case it will be “heated up” by the noise. For this reason, the SD version with $`\gamma _i0`$ allows for a slightly larger time step than the pure Andersen method, while the Feller et al. version ($`\gamma _i=0`$, $`\gamma _V0`$) is not more stable, involving only global thermostatting. We have chosen

$$\gamma _i=\gamma _0L^2,$$
(29)   
while choosing a constant value for $`\gamma _V`$. Then the Langevin equation for $`\stackrel{}{\pi }_i`$ is rewritten as

$$\dot{\stackrel{}{\pi }_i}=L\left(\stackrel{}{f}_i\frac{\gamma _0}{m_i}\frac{\stackrel{}{\pi }_i}{L}+\sqrt{k_BT\gamma _0}\stackrel{}{\eta }_i(t)\right);$$
(30)   
this fits naturally to standard SD, where the friction force is $`\gamma _0\stackrel{}{p}_i/m_i=\gamma _0\stackrel{}{\pi }_i/(Lm_i)`$. Further details on the implementation will be given in Sec. V; there we will also discuss the choice of the parameters $`\gamma _0`$ and $`\gamma _V`$, as well as of the piston mass $`Q`$, which are all irrelevant for the statistical properties of the system, but of great importance for the equilibration properties of the algorithm.

## V Implementation

### A Symplectic Integrator

Symplectic time–reversible integrators are known to be extremely useful for molecular dynamics simulations of Hamiltonian systems . This is so because, except for roundoff errors, they conserve the phase–space volume exactly (note the intimate relation to entropy production), and do not mark any particular direction of time (note that a global drift in the algorithm breaks this time–reversal symmetry). Hence they are very stable and allow for a large time step. The most common example (actually, the lowest–order scheme) is the well–known Verlet algorithm in its various formulations .

A particularly transparent way to construct these algorithms  is based on the Liouville operator $`i\widehat{L}`$, which is just the Fokker–Planck operator in the special case of vanishing friction (see Eqn. 11). Noticing that the formal solution of Eqn. 8 is just

$$\mathrm{\Phi }(t)=\mathrm{exp}\left(i\widehat{L}t\right)\mathrm{\Phi }(0),$$
(31)   
one focuses directly on the operator $`\mathrm{exp}\left(i\widehat{L}t\right)`$, and decomposes $`\widehat{L}`$ into a sum of simpler operators,

$$i\widehat{L}=\underset{k=1}{\overset{n}{}}i\widehat{L}_k,$$
(32)   
where the imaginary unit is extracted for convenience; $`\widehat{L}`$ as well as each of the $`\widehat{L}_k`$ is self–adjoint. Now, the exact time development within a time step $`\mathrm{\Delta }t`$ is replaced by an approximate one by using the factorization

$$e^{i\widehat{L}\mathrm{\Delta }t}e^{i\widehat{L}_1\mathrm{\Delta }t/2}e^{i\widehat{L}_2\mathrm{\Delta }t/2}\mathrm{}e^{i\widehat{L}_{n1}\mathrm{\Delta }t/2}e^{i\widehat{L}_n\mathrm{\Delta }t}e^{i\widehat{L}_{n1}\mathrm{\Delta }t/2}\mathrm{}e^{i\widehat{L}_2\mathrm{\Delta }t/2}e^{i\widehat{L}_1\mathrm{\Delta }t/2}.$$
(33)   
This scheme is automatically phase–space conserving, since each of the operators is unitary, and time–reversible symmetric, since the inverse operator is just the original operator, evaluated for $`\mathrm{\Delta }t`$. Now, if each of the operators $`\widehat{L}_k`$ is simple enough, the action of $`e^{i\widehat{L}_kt}`$ on a phase space point can be calculated trivially, such that the algorithm is a succession of simple updating steps.

Specifically, for the Andersen equations of motion (see Sec. III) we choose the operators

$`i\widehat{L}_1`$ $`=`$ $`{\displaystyle \underset{i}{}}L\stackrel{}{f}_i{\displaystyle \frac{}{\stackrel{}{\pi }_i}}`$ (34)   
$`i\widehat{L}_2`$ $`=`$ $`\left(𝒫P\right){\displaystyle \frac{}{\mathrm{\Pi }_V}}`$ (35)   
$`i\widehat{L}_3`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Pi }_V}{Q}}{\displaystyle \frac{}{V}}`$ (36)   
$`i\widehat{L}_4`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{\stackrel{}{\pi }_i}{L^2m_i}}{\displaystyle \frac{}{\stackrel{}{s}_i}},`$ (37)   
resulting in the following updating scheme:

1. $`\stackrel{}{\pi }_i(t)\stackrel{}{\pi }_i(t+\mathrm{\Delta }t/2)=\stackrel{}{\pi }_i(t)+L(t)\stackrel{}{f}_i(t)\mathrm{\Delta }t/2`$
2. $`\mathrm{\Pi }_V(t)\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)=\mathrm{\Pi }_V(t)+\left(𝒫P\right)\mathrm{\Delta }t/2`$

(note that for the evaluation of $`𝒫`$, one has to take the old positions $`\stackrel{}{s}_i(t)`$ and the old box size $`L(t)`$, but already the updated momenta $`\stackrel{}{\pi }_i(t+\mathrm{\Delta }t/2)`$)
3. $`V(t)V(t+\mathrm{\Delta }t/2)=V(t)+Q^1\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t/2`$
4. $`\stackrel{}{s}_i(t)\stackrel{}{s}_i(t+\mathrm{\Delta }t)=\stackrel{}{s}_i(t)+\frac{\stackrel{}{\pi }_i(t+\mathrm{\Delta }t/2)}{L^2(t+\mathrm{\Delta }t/2)m_i}\mathrm{\Delta }t`$
5. $`V(t+\mathrm{\Delta }t/2)V(t+\mathrm{\Delta }t)=V(t+\mathrm{\Delta }t/2)+Q^1\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t/2`$
6. $`\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)\mathrm{\Pi }_V(t+\mathrm{\Delta }t)=\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)+\left(𝒫P\right)\mathrm{\Delta }t/2`$

(for evaluation of $`𝒫`$, one uses $`\stackrel{}{s}_i(t+\mathrm{\Delta }t)`$, $`L(t+\mathrm{\Delta }t)`$, and $`\stackrel{}{\pi }_i(t+\mathrm{\Delta }t/2)`$)
7. $`\stackrel{}{\pi }_i(t+\mathrm{\Delta }t/2)\stackrel{}{\pi }_i(t+\mathrm{\Delta }t)=\stackrel{}{\pi }_i(t+\mathrm{\Delta }t/2)+L(t+\mathrm{\Delta }t)\stackrel{}{f}_i(t+\mathrm{\Delta }t)\mathrm{\Delta }t/2`$ .

It is often convenient to formulate the algorithm in terms of the conventional variables

$$\stackrel{}{r}_i(t)=L(t)\stackrel{}{s}_i(t)\stackrel{}{p}_i(t)=L(t)^1\stackrel{}{\pi }_i(t),$$
(38)   
which are however not canonically conjugate with respect to each other. The pressure, in terms of these variables, is written as

$$𝒫=\frac{1}{dV}\underset{i<j}{}\stackrel{}{f}_{ij}\stackrel{}{r}_{ij}+\frac{1}{dV}\underset{i}{}\frac{1}{m_i}\stackrel{}{p}_i^{\mathrm{\hspace{0.17em}2}},$$
(39)   
and the updating scheme, which now involves various rescaling steps, proceeds as follows:

1. $`\stackrel{}{p}_i^{}=\stackrel{}{p}_i(t)+\stackrel{}{f}_i(t)\mathrm{\Delta }t/2`$
2. $`𝒫`$ is evaluated using Eqn. 39 with $`\stackrel{}{r}_i(t)`$, $`L(t)`$, and $`\stackrel{}{p}_i^{}`$; then $`\mathrm{\Pi }_V`$ is updated as before:

$`\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)=\mathrm{\Pi }_V(t)+\left(𝒫P\right)\mathrm{\Delta }t/2`$
3. $`V(t+\mathrm{\Delta }t/2)=V(t)+Q^1\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t/2`$
4. $`\stackrel{}{r}_i^{}=\stackrel{}{r}_i(t)+\frac{L^2(t)}{L^2(t+\mathrm{\Delta }t/2)}\frac{\stackrel{}{p}_i^{}}{m_i}\mathrm{\Delta }t`$
5. $`V(t+\mathrm{\Delta }t)=V(t+\mathrm{\Delta }t/2)+Q^1\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t/2`$

followed by two rescaling steps:

* $`\stackrel{}{r}_i(t+\mathrm{\Delta }t)=\frac{L(t+\mathrm{\Delta }t)}{L(t)}\stackrel{}{r}_i^{}`$
    * $`\stackrel{}{p}_i^{\prime \prime }=\frac{L(t)}{L(t+\mathrm{\Delta }t)}\stackrel{}{p}_i^{}`$
6. $`𝒫`$ is evaluated using Eqn. 39 with $`\stackrel{}{r}_i(t+\mathrm{\Delta }t)`$, $`L(t+\mathrm{\Delta }t)`$, and $`\stackrel{}{p}_i^{\prime \prime }`$; then

$`\mathrm{\Pi }_V(t+\mathrm{\Delta }t)=\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)+\left(𝒫P\right)\mathrm{\Delta }t/2`$
7. $`\stackrel{}{p}_i(t+\mathrm{\Delta }t)=\stackrel{}{p}_i^{\prime \prime }+\stackrel{}{f}_i(t+\mathrm{\Delta }t)\mathrm{\Delta }t/2`$ .

So far the algorithm has been developed for the case without friction and noise. For the case with friction and noise, we simply use the scheme given above, and introduce the following replacements:

$`\stackrel{}{f}_i\mathrm{\Delta }t/2`$ $``$ $`\stackrel{}{f}_i\mathrm{\Delta }t/2\gamma _0{\displaystyle \frac{\stackrel{}{p}_i}{m_i}}\mathrm{\Delta }t/2+\sqrt{k_BT\gamma _0\mathrm{\Delta }t}\stackrel{}{z}_i`$ (40)   
$`\left(𝒫P\right)\mathrm{\Delta }t/2`$ $``$ $`\left(𝒫P\right)\mathrm{\Delta }t/2\gamma _V{\displaystyle \frac{\mathrm{\Pi }_V}{Q}}\mathrm{\Delta }t/2+\sqrt{k_BT\gamma _V\mathrm{\Delta }t}z_V.`$ (41)   
Here $`z_i^\alpha `$ and $`z_V`$ denote uncorrelated random numbers with zero mean and unit variance; for simplicity we sample them from a uniform distribution via

$$z=\sqrt{12}\left(u\frac{1}{2}\right)$$
(42)   
where $`u`$ is uniformly distributed on the unit interval. The momenta which occur in Eqn. 40 are, for simplicity, taken as $`\stackrel{}{p}_i(t)`$ in step (1), $`\mathrm{\Pi }_V(t)`$ in step (2), $`\mathrm{\Pi }_V(t+\mathrm{\Delta }t/2)`$ in step (6), and $`\stackrel{}{p}_i^{\prime \prime }`$ in step (7).

### B Choice of Parameters

We start from the observation that a molecular system is characterized by a typical molecular frequency $`\omega _0`$, which can be viewed as the “Einstein” frequency of oscillations of an atom in its “cage” . With use of the intermolecular potential $`v(\stackrel{}{r})`$ one gets

$`\omega _0^2={\displaystyle \frac{\rho }{dm}}{\displaystyle d^d\stackrel{}{r}g(\stackrel{}{r})^2v(\stackrel{}{r})},`$ (43)   
with $`\rho `$ being the particle number density, $`m`$ the mass of the molecules, and $`g(\stackrel{}{r})`$ the pair distribution function. Alternatively, one can define a molecular time scale by the time which a sound wave needs for traveling the nearest neighbor distance . However, both frequencies coincide by order of magnitude. This frequency governs the time step $`\mathrm{\Delta }t`$ which one has to choose in order to keep the MD algorithm stable; a typical rule of thumb says $`\mathrm{\Delta }t=(1/50)(2\pi /\omega _0)`$.

Similarly, the piston degree of freedom performs oscillations, if it is simulated with very weak friction in the NVT ensemble. Following Nosé , we can estimate their frequency $`\mathrm{\Omega }_0`$ quite easily. Within a linearized approximation, the isothermal compressibility

$$\kappa _T=\frac{1}{V}\frac{V}{P}=\frac{1}{Vk_BT}\left(V^2V^2\right)$$
(44)   
controls the relation between pressure fluctuations $`\delta 𝒫=𝒫P`$ and volume fluctuations $`\delta V=VV`$ via

$$\delta 𝒫=\frac{P}{V}\delta V=\frac{1}{V\kappa _T}\delta V.$$
(45)   
Therefore, one concludes from Eqn. 19 by trivial insertion

$$\frac{d^2}{dt^2}\delta V=\frac{1}{QV\kappa _T}\delta V,$$
(46)   
which is the equation of motion of a harmonic oscillator with frequency

$$\mathrm{\Omega }_0^2=\frac{1}{QV\kappa _T}.$$
(47)   
Obviously, the piston mass has to be chosen small enough such that the system can adjust its volume sufficiently fast. On the other hand, it cannot be chosen too small, since then $`\mathrm{\Omega }_0`$ becomes too large, see Eqn. 47. Clearly, one does not want $`\mathrm{\Omega }_0`$ to exceed $`\omega _0`$, since otherwise the simulation would need an unnecessarily small time step. The optimum piston mass is thus found from the resonance condition $`\mathrm{\Omega }_0=\omega _0`$ , which yields a quite different value for $`Q`$ (by a factor of $`L^2`$) than Andersen’s original suggestion  — this original criterion has turned out to be not correct. The similar frequencies of the molecular oscillator and the box volume lead to a very quick energy transfer between them, resulting in a very efficient equilibration. However, one will often choose a substantially larger value for $`Q`$ in order to separate the time scales, such that the molecular motion on short length and time scales is largely unaffected by the piston motion. Regardless of the precise choice for $`Q`$, one should note that keeping $`\mathrm{\Omega }_0`$ constant implies a scaling of $`Q`$ with the inverse system size.

When the coupling to the heat bath with friction and noise is added, the question arises how to choose the damping parameters $`\gamma _0`$ and $`\gamma _V`$. Let us hence study a (deterministic) damped harmonic oscillator,

$$m\ddot{x}+\gamma \dot{x}+m\omega _0^2x=0.$$
(48)   
Obviously, for small $`\gamma `$ the damping can practically be neglected. On the other hand, for $`\gamma m\omega _0`$, damping force and harmonic force are of the same order of magnitude (for a harmonic oscillator, we can estimate the velocity via $`\dot{x}\omega _0x`$). The exact calculation yields

$$\gamma _c=2m\omega _0$$
(49)   
for the “aperiodic limit”. This value quantifies the qualitative distinction between “weak” and “strong” damping, denoting the boundary between oscillatory behavior ($`\gamma <\gamma _c`$) and pure relaxational dynamics ($`\gamma >\gamma _c`$). Only in the weak damping case $`\gamma \gamma _c`$, the fastest time scale (which governs the time step) is given by $`\omega _0`$, while for $`\gamma \gamma _c`$ the damping term dominates, requiring a smaller time step. For this reason, $`\gamma `$ values beyond $`\gamma _c`$ are clearly undesirable. Even worse, for $`\gamma \gamma _c`$ the relaxation contains also a very slow component, whose characteristic time is, for large $`\gamma `$, given by $`\gamma /(m\omega _0^2)`$. Therefore, $`\gamma =\gamma _c`$ is clearly the optimum value for fast equilibration.

However, for the single–particle damping $`\gamma _0`$ we typically choose a value which is between one and two orders of magnitude smaller than $`\gamma _c`$. This is in accord with the philosophy of simulating the system rather close to its true Hamiltonian dynamics, such that at least on the local scales, both spatially and temporally, the dynamics can be considered as realistic. We hence use the coupling to the heat bath mainly for additional stabilization, deliberately keeping the molecular oscillations in the simulation.

On the other hand, there is no analogous argument of “physical realism” for the box motion, which is intrinsically unphysical. One is therefore clearly led to the choice $`\gamma _V=\gamma _{Vc}`$, such that the ringing is just avoided, while the box motions are still sufficiently fast. The same conclusion was obtained by Feller et al. , where however no theoretical background was provided.

Using this choice, one is left with only one time scale for the piston motion, given by $`\mathrm{\Omega }_0`$, and this is in turn adjusted according to the needs of the simulation: If one is only interested in statics, then “resonant” coupling is desirable (i. e. $`\mathrm{\Omega }_0\omega _0`$, and also $`\gamma _0\gamma _c`$), at the expense of distorting the motions even on the molecular scale. If, on the other hand, it is desired to realistically simulate the molecular oscillations, one should enforce a separation of time scales by choosing both $`\mathrm{\Omega }_0\omega _0`$ and $`\gamma _0\gamma _c`$.

### C Numerical Test

We study a system containing 100 particles interacting via a truncated Lennard–Jones potential whose attractive part is cut off:

$$U_{LJ}(r)=4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6+\frac{1}{4}\right];r<2^{1/6}\sigma .$$
(50)   
We choose Lennard–Jones units where the parameters $`\sigma `$ and $`ϵ`$ as well as the particle mass $`m`$ are set to unity.

Figure 1 illustrates the problem of the “ringing” in this particular system, by displaying the autocorrelation function of the pressure fluctuations for various simulation parameters. The simulated temperature is $`k_BT=1.0`$, while the pressure was fixed at the rather low value $`P=1.0`$. The mean volume for the 100 particles is $`V=262.7`$, corresponding to a density $`\rho =0.38`$. The compressibility at this state point is $`\kappa _T=0.3`$, such that for a box mass of $`Q=0.1`$ one finds $`\mathrm{\Omega }_00.36`$ for the ringing frequency or $`2\pi /\mathrm{\Omega }_018`$ for the oscillation time. The figure shows that the box indeed oscillates, and that the frequency of the oscillations has been estimated correctly. The left part of the figure is for pure undamped Andersen dynamics where the dependence on $`Q`$ is shown. Choosing a value of $`Q=0.0001`$ leads to a very fast relaxation of the pressure autocorrelation function. The theoretical prediction for the best box mass for the molecular frequency of $`\omega _08.5`$ in this system leads to $`Q_{opt}=1/(V\kappa _T\omega _0^2)0.00018`$. But as seen in Fig. 2, even for this value of $`Q`$ the oscillations still remain on a very short time scale. Only a value of $`\gamma _V=0.001`$, close to $`\gamma _{Vc}=2Q\mathrm{\Omega }_00.002`$ suppresses the fluctuations efficiently and the autocorrelation function resembles the autocorrelation function in the NVE ensemble. Interestingly, it is also seen that for constant volume the pressure relaxation is considerably slower if the molecular damping is turned on. Conversely, this behavior is practically absent in the constant pressure case, where actually the “best” autocorrelation function was found for $`\gamma _V=0.001`$ (as discussed), combined with some additional molecular friction $`\gamma _0=0.5`$.

The statistical accuracy of an observable is given by the ratio between simulation time and the integrated autocorrelation time $`\tau `$, i. e. the value of the time integral over the normalized autocorrelation function . From that perspective, a slow decay with many oscillations is actually not particularly harmful, since the integral value is rather small, due to cancellation. However, the result of Ref.  holds only in the asymptotic limit where the simulation time is substantially longer than the decay of the correlation function. Moreover, the numerical integration of an oscillatory function converges only slowly and is hence rather awkward. For these reasons, a simulation algorithm which avoids oscillations is clearly preferable. In order to illustrate this further, Table I lists $`\tau `$ for the autocorrelation functions shown in Fig. 2. The smallest $`\tau `$ is actually found for the pure Andersen NPH simulation. However, turning on the damping increases $`\tau `$ only by less than a factor of two, while the decay time decreases nearly by a factor of ten.

The case of a box whose mass has, for reasons of separation of time scales, been chosen substantially larger than the optimum value, is illustrated in the upper right part of Fig. 1 ($`Q=0.1`$, i. e. three orders of magnitude larger than $`Q_{opt}`$). The autocorrelation function decays most rapidly for $`\gamma _V=0.1`$, as theoretically expected ($`\gamma _{Vc}=0.072`$). Compared to undamped dynamics at the same mass, this is a considerable improvement. Nevertheless, this decay is still substantially slower than what one can obtain if also the mass is chosen optimally (Fig. 2) — this is simply the price which is being paid for achieving realistic molecular motion. For practical purposes, the decay obtained for $`Q=0.1`$, $`\gamma _V=0.1`$ is quite acceptable.

## VI Conclusion

We have discussed the algorithm of Andersen  and the generalization to stochastic piston motion by Feller et al. , generalizing it even further to also include stochastic motion of the particles. We gave a straightforward proof that the NPT ensemble is produced. The implementation by means of a symplectic algorithm is particularly stable and well–suited for MD problems. Another important point, the choice of the right simulation parameters, was studied both theoretically and numerically, and a guideline for their optimum values was given. We view this algorithm as a particularly useful realization of the constant–pressure ensemble.

## Acknowledgments

Fruitful and critical discussions with Patrick Ahlrichs, Markus Deserno, Alex Bunker and Kristian Müller–Nedebock are gratefully acknowledged.