Stochastic NVT Thermostats: Andersen, Langevin, and Dissipative Particle Dynamics

The Andersen Scheme. Perhaps the simplest thermostat which does correctly sample the NVT ensemble is due to Andersen [11]. Here, at each step, some prescribed number of particles is selected, and their momenta (actually, their velocities) are drawn from a Gaussian distribution at the prescribed temperature:

$$\mathscr{P}(p) = \left(\frac{\beta}{2\pi m}\right)^{3/2} \exp\left[-\beta p^2/\left(2m\right)\right]$$ (179)

This is intended to mimic collisions with bath particles at a specified $T$. The strength of the coupling to the heat bath is specified by a collision frequency, $\nu$. For each particle, a random variate is selected between 0 and 1. If this variate is less than $\nu\Delta t$, then that particle's momenta are reset.

The code mdlj_and.c implements the Andersen thermostat with the velocity Verlet algorithm for the Lennard-Jones fluid. I was curious about the importance placed on the fact that a Gaussian distribution of velocities is established during an MD simulation using the Andersen thermostat. F&S seem to indicate that the truth of this fact means that we are sampling the canonical ensemble. But they do not show data for the Berendsen thermostat, for which we should not sample NVT. The figure below is a comparison of the measured velocity distribution for a 20,000-step run at $T$ = 1.0 and $\rho = 0.8442$ using the Berendsen thermostat with $\tau = 0.1$ and a similar run with the Andersen thermostat with $\nu = 0.01$. Apparently, the Berendsen thermostat also reproduces the correct distribution.

Distributions of velocity computed from MD simulations using either the Berendsen ($\tau = 0.1$) or Andersen ($\nu = 0.01$) thermostats. $T$ = 1.0, $\rho = 0.8442$, $N$ = 256.

Although they may not be on the right track with regard to the velocity distributions, F&S do make a very good point with regard to measuring transport properties from NVT MD simulations. The Andersen thermostat destroys momentum transport because of the random velocities; hence, there is no continuity of momentum in an Andersen LJ fluid, and therefore no proper $\mathscr{D}$ or viscosity. The data they cite in Fig. 6.3 clearly shows that $\mathscr{D}$, if measured from an Andersen MD run, is incorrect.

The Langevin thermostat. In the ``Langevin'' thermostat, at each time step all particles receive a random force and have their velocities lowered using a constant friction. [12] The average magnitude of the random forces and the friction are related in a particular way, which guarantees that the ``fluctuation-dissipation'' theorem is obeyed, thereby guaranteeing NVT statistics.

In this formalism, the particle-$i$ equation of motion is modified:

$$m\ddot{\bf r}_i = -\nabla_i U - m\Gamma\dot{\bf r}_i + {\bf W}_i\left(t\right)$$ (180)

Here, $\Gamma$ is a friction coefficient with units of $\tau^{-1}$, and ${\bf W}_i$ is a random force that is uncorrelated in time and across particles, with a mean given by

$$\left<{\bf W}_i\left(t\right),{\bf W}_j\left(t^\prime\right)\right> = \delta_{ij}\delta\left(t-t^\prime\right)6k_BmT\Gamma$$ (181)

The code mdlj_lan.c implements the Langevin thermostat. The two major elements are a force initialization at each time step that adds in the random forces, ${\bf W}$, and a slight modification to the update equations in the integrator to include the effect of $\Gamma$. Note that the initialization of forces in the force routine has been removed.

One advantage of the Langevin thermostat (and to a limited extent, the Andersen thermostat and other stochastic-based thermostats) is that we can get away with a larger time step than in NVE simulations. At a density of $rho$ = 0.8442 and a mean temperature $T$ = 1.0, an NVE simulation is unstable for time-steps above about $\Delta t$ = 0.004. We can, however, run a Langevin dynamics simulation with a friction $\Gamma$ = 1.0 stably with a time-step as large as $\Delta t$ = 0.01 or even higher. This has proven invaluable in simulations of more complicated systems that simple liquids, namely linear polymers, which have very long relaxation times. MD with the Langevin thermostat is the method of choice for equilibrating samples of liquids of long bead-spring polymer chains.

Of course, the drawback of most stochastic thermostats (one exception is discussed next) is that momentum transfer is destroyed. So again, it is unadvisable to use Langeving or Andersen thermostats for runs in which you wish to compute diffusion coefficients. I echo the recommendation of F&S: when possible, use NVE to compute properties, and use thermostats for equilibration only.

The Dissipative Particle Dynamics thermostat. The DPD thermostat [13,14] adds pairwise random and dissipative forces to all particles, and has been shown to preserve momentum transport. Hence, it is the only stochastic thermostat so far that should even be considered for use if one wishes to compute transport properties.

The DPD thermostat is implemented by slight modification of the force routine to add in the pairwise random and dissipative forces. For the $ij$ pair, the dissipative force is defined as

$${\bf f}_{ij}^D = -\gamma\omega^D\left(r_{ij}\right)\left({\bf v}_{ij}\cdot\hat{\bf r}_{ij}\right)\hat{\bf r}_{ij}$$ (182)

Here, $\gamma$ is a friction coefficient, $\omega$ is a cut-off function for the force as a function of the scalar distance between $i$ and $j$ which simply limits the interaction range of the dissipative (and random) forces, ${\bf v}_{ij} = {\bf v}_i - {\bf v}_j$ is the relative velocity of $i$ to $j$, and $\hat{\bf r}_{ij} = {\bf r}_{ij}/r_{ij}$ is the unit vector pointing from $j$ to $i$. The random force is defined as

$${\bf f}_{ij}^R = \sigma\omega^R\left(r_{ij}\right)\zeta_{ij}\hat{\bf r}_{ij}$$ (183)

Here, $\sigma$ is the strength of the random force, $\omega^R$ is a cut-off, and $\zeta_{ij}$ is a Gaussian random number with zero mean and unit variance, and $\zeta_{ij} = \zeta_{ji}$.

The update of velocity uses these new forces:

$${\bf v}_i\left(t + \Delta t\right) = {\bf v}_i\left(t\right)... ...\Delta t}{m}{\bf f}_i^D + \frac{\sqrt{\Delta t}}{m}{\bf f}_i^R$$ (184)

where

$${\bf f}_i^D = \sum_{j\ne i} {\bf f}_{ij}^D$$ (185)

$${\bf f}_i^R = \sum_{j\ne i} {\bf f}_{ij}^R$$ (186)

The parameters $\gamma$ and $\sigma$ are linked by a fluctuation-dissipation theorem:

$$\sigma^2 = 2\gamma k_B T$$ (187)

So, in practice, one must specify either $\gamma$ or $\sigma$, and then a setpoint temperature, $T$, in order to use the DPD thermostat.

The cutoff functions are also related:

$$\omega^D\left(r_{ij}\right) = \left[\omega^R\left(r_{ij}\right)\right]^2$$ (188)

This is the only real constraint on the cutoffs; we are otherwise allowed to use any cutoff we like. The simplest uses the cutoff radius of the pair potential, $r_c$:

$$\omega\left(r\right) = \left\{\begin{tabular}{ll} 1 & \mbox{$r < r_c$}\\ 0 & \mbox{$r > r_c$} \end{tabular}\right.$$ (189)

Note that, with this choice, $\left[\omega^R\left(r_{ij}\right)\right]^2$ = $\omega^R\left(r_{ij}\right)$ = $\omega^D\left(r_{ij}\right) = \omega$.

The code mdlj_dpd.c implements the DPD thermostat in an MD simulation of the Lennard-Jones liquid. The major changes (compared to mdlj.c) are to the force routine, which now requires several more arguments, including particle velocities, and parameters for the thermostat. Inside the pair loop, the force on each particle is updated by the conservative, dissipative, and random pairwise force components. The random force is divided by $\sqrt{\Delta t}$ so that the velocity Verlet algorithm need not be altered to implement Eq. 184.

The behavior of the DPD thermostat can be assessed in a similar fashion as was the Berendsen thermostat above. Below is a lin-log plot of system temperature after the setpoint is changed from 1.0 to 2.0, for a Lennard-Jones system $N$ = 256 and $rho$ = 0.8442. Again, each curve corresponds to a different value of $\gamma$, and they increase by factors of 10. The corresponding time at which the setpoint $T$ is reached is also seen to increase by the same factor.

Instantaneous temperature, $T$, vs. time in an MD simulation of 256 particles at a density of 0.8442, with temperature controlled by the dissipative particle dynamics thermostat [13,14], for various values of the parameter $\gamma$.

If you have been paying attention, you will notice that in all implementations given here, the thermostat is applied each and every time step. In principle, one could modulate the strength of a thermostat by controlling the frequency with which it is applied. As an exercise, you can pick one of these implementations and modify it such that the thermostat is applied once per $M$ time steps, where $M$ is specified on the command-line. You can compare the performance of the thermostat in the same way we have already seen (time required to go from $T$ = 1.0 to $T$ = 2.0 upon an instantaneous change in the setpoint $T$) for a given friction or strength and for various values of $M$.