Molecular Dynamics at Constant Pressure: The Berendsen Barostat

As with temperature control, there are different classes of pressure control for MD simulation. The only one we consider here is the length-scaling technique of Berendsen. It should be noted that one can also use the the extended Nosé-Hoover (extended Lagrangian) formalism of Martyna, which is mentioned in F&S; in the interest of time, we will forego a discussion of this technique.

Here we consider implementation of the Berendsen barostat [10]. Recall that the working definition of instantaneous pressure, $P$, is given by:

$$P = \rho T + {\rm vir}/V$$ (196)

where ${\rm vir}$ is the virial:

$${\rm vir} = \frac{1}{3}\sum_{i>j} {\bf f}\left(r_{ij}\right)\cdot{\bf r}_{ij}$$ (197)

and $V$ is the system volume. ${\bf f}\left(r_{ij}\right)$ is the force exerted on particle $i$ by particle $j$.

Consider a cubic system, where $V = L^3$. The Berendsen barostat uses a scale factor, $\mu$, which is a function of $P$, to scale lengths in the system:

$${\bf r}_i \rightarrow \mu{\bf r}_i$$ (198)

$$L \rightarrow \mu L$$ (199)

$\mu$ is given by

$$\mu = \left[1-\frac{\Delta t}{\tau_P}\left(P - P_0\right)\right]^{1/3}$$ (200)

Here, $\Delta t$ is the integrator time-step, $\tau _P$ is the ``rise time'' of the barostat, and $P_0$ is the setpoint pressure. Berendsen discusses the tensor-based analog for non-cubic systems [10].

The code mdlj_berp.c implements the Berendsen barostat. Below, I show results of using the Berendsen barostat to induce a pressure jump from 1.0 to 6.0 in a sample of LJ fluid, for various values of the rise time, $\tau _P$. Notice that temperature is not controlled, but rises from about 1.3 to 2.5 due to the increase of pressure.

Instantaneous pressure, $P$, vs. time (upper), and Instantaneou temperature, $T$, vs. time (lower) in an MD simulation of 256 particles at a density of 0.8442, using the Berendsen barostat [10] to impose an instantaneous pressure jump from 1.0 to 6.0. Each curve corresponds to a different value of the rise time, $\tau _P$.

Length scaling at each time step using a global scale factor, while effective in this instance, can be lead to violent oscillations of pressure in more ordered systems, and is therefore not recommended for production MD runs. However, it is common to find length scaling barostats used in the literature without reporting how effective they are, measured at least in terms of pressure and its fluctuations. But they can be useful for pre-equilibrating samples at some $P$ prior to beginning an NVE simulation during which one hopes the instantaneous pressure fluctuates about the previous setpoint. It is easy to implement both the Berendsen thermostat and barostat in the same simulation program, to allow pre-equilibration at setpoint $T$ and $P$.