Ensembles¶

The aim of the time evolution in MD simulations is to find the phase trajectory

\[\{ \boldsymbol{r}_i(t_1), ~\boldsymbol{v}_{i}(t_1)\}_{i=1}^N,~ \{ \boldsymbol{r}_i(t_2), ~\boldsymbol{v}_{i}(t_2)\}_{i=1}^N,~ \cdots\]

starting from the initial phase point

\[\{ \boldsymbol{r}_i(t_0), ~\boldsymbol{v}_{i}(t_0)\}_{i=1}^N.\]

The time interval between two time points \(\Delta t=t_1-t_0=t_2-t_1=\cdots\) is called the time step.

The algorithm for integrating by one step depends on the ensemble type and other external conditions. There are many ensembles used in MD simulations, but we only consider the following 3 in the current version:

NVE ensemble¶

The NVE ensmeble is also called the micro-canonical ensemble. We use the velocity-Verlet integration method with the following equations:

\[\begin{split}\boldsymbol{v}_i(t_{m+1}) \approx \boldsymbol{v}_i(t_{m}) + \frac{\boldsymbol{F}_i(t_m)+\boldsymbol{F}_i(t_{m+1})}{2m_i}\Delta t \\ \boldsymbol{r}_i(t_{m+1}) \approx \boldsymbol{r}_i(t_{m}) + \boldsymbol{v}_i(t_m) \Delta t + \frac{1}{2} \frac{\boldsymbol{F}_i(t_m)}{m_i} (\Delta t)^2.\end{split}\]

Here, \(\boldsymbol{v}_i(t_{m})\) is the velocity vector of particle \(i\) at time \(t_{m}\). \(\boldsymbol{r}_i(t_{m})\) is the position vector of particle \(i\) at time \(t_{m}\). \(\boldsymbol{F}_i(t_{m})\) is the force vector of particle \(i\) at time \(t_{m}\). \(m_i\) is the mass of particle \(i\) and \(\Delta t\) is the time step for integration.

NVT ensemble¶

The NVT ensemble is also called the canonical ensemble. We have implemented several thermostats for the NVT ensemble.

Berendsen thermostat¶

The velocities are scaled in the Berendsen thermostat [Berendsen1984] as follows:

\[\boldsymbol{v}_i^{\text{scaled}} = \boldsymbol{v}_i \sqrt{1 + \frac{\Delta t}{\tau_T} \left(\frac{T_0}{T} - 1\right)}.\]

Here, \(\tau_T\) is the coupling time (relaxation time) parameter, \(T_0\) is the target temperature, and \(T\) is the instant temperature calculated from the current velocities. We require that \(\tau_T/\Delta t \geq 1\).

When \(\tau_T/\Delta t = 1\), the above formula reduces to the simple velocity-scaling formula:

\[\boldsymbol{v}_i^{\text{scaled}} = \boldsymbol{v}_i \sqrt{\frac{T_0}{T}}.\]

A larger value of \(\tau_T/\Delta t\) represents a weaker coupling between the system and the thermostat. We recommend a value of \(\tau_T/\Delta t \approx 100\). The above re-scaling is applied at each time step after the velocity-Verlet integration. This thermostat is usually used for reaching equilibrium and is not recommended for sampling the canonical ensemble.

Nose-Hoover chain thermostat¶

In the Nose-Hoover chain (NHC) method [Tuckerman2010], the equations of motion for the particles in the thermostatted region are (those for the thermostat variables are not presented):

\[\begin{split}\frac{d \boldsymbol{r}_i}{dt} = \frac{\boldsymbol{p}_i}{m_i} \\ \frac{d \boldsymbol{p}_i}{dt} = \boldsymbol{F}_i - \frac{\pi_0}{Q_0} \boldsymbol{p}_i.\end{split}\]

Here, \(\boldsymbol{r}_i\) is the position of atom \(i\). \(\boldsymbol{p}_i\) is the momentum of atom \(i\). \(m_i\) is the mass of atom \(i\). \(\boldsymbol{F}_i\) is the total force on atom \(i\) resulting from the potential model used. \(Q_0=N_{\rm f} k_{\rm B} T_0 \tau_T^2\) is the ‘’mass’’ of the thermostat variable directly coupled to the system and \(\pi_0\) is the corresponding ‘’momentum’’. \(N_{\rm f}\) is the degree of freedom in the thermostatted region. \(k_{\rm B}\) is Boltzmann’s constant and \(T_0\) is the target temperature. \(\tau_T\) is a time parameter, and we suggest a value of \(\tau_T/\Delta t \approx 100\), where \(\Delta t\) is the time step.

We use a fixed chain length of 4.

Langevin thermostat¶

In the Langevin method, the equations of motion for the particles in the thermostatted region are

\[\begin{split}\frac{d \boldsymbol{r}_i}{dt} = \frac{\boldsymbol{p}_i}{m_i} \\ \frac{d \boldsymbol{p}_i}{dt} = \boldsymbol{F}_i - \frac{\boldsymbol{p}_i}{\tau_T} + \boldsymbol{f}_i,\end{split}\]

Here, \(\boldsymbol{r}_i\) is the position of atom \(i\). \(\boldsymbol{p}_i\) is the momentum of atom \(i\). \(m_i\) is the mass of atom \(i\). \(\boldsymbol{F}_i\) is the total force on atom \(i\) resulted from the potential model used. \(\boldsymbol{f}_i\) is a random force with a variation determined by the fluctuation-dissipation relation to recover the canonical ensemble distribution with the target temperature. \(\tau_T\) is a time parameter, and we suggest a value of \(\tau_T/\Delta t \approx 100\), where \(\Delta t\) is the time step. We implemented the integrators proposed in [Bussi2007a] and [Leimkuhler2013].

Bussi-Donadio-Parrinello thermostat¶

The Berendsen thermostat does not generate a true NVT ensemble. As an extension of the Berendsen thermostat, the Bussi-Donadio-Parrinello (BDP) thermostat [Bussi2007b] incorporates a proper randomness into the velocity re-scaling factor and generates a true NVT ensemble. It is also called the stochastic velocity rescaling (SVR) thermostat.

In the BDP thermostat, the velocities are scaled in the following way:

\[\boldsymbol{v}_i^{\text{scaled}} = \alpha \boldsymbol{v}_i\]

where

\[\alpha^2= e^{-\Delta t/\tau_T} + \frac{T_0}{TN_f} \left( 1-e^{-\Delta t/\tau_T} \right) \left( R_1^2 + \sum_{i=2}^{N_f}R_i^2 \right) + 2e^{-\Delta t/2\tau_T} R_1 \sqrt{\frac{T_0}{TN_f} \left( 1-e^{-\Delta t/\tau_T} \right) }.\]

Here, \(\boldsymbol{v}_i\) is the velocity of atom \(i\) before the re-scaling. \(N_{\rm f}\) is the degree of freedom in the thermostatted region. \(T\) is instant temperature and \(T_0\) is the target temperature. \(\Delta t\) is the time step for integration. \(\tau_T\) is a time parameter, and we suggest a value of \(\tau_T/\Delta t \approx 100\), where \(\Delta t\) is the time step. \(\{R_i\}_{i=1}^{N_f}\) are \(N_{\rm f}\) Gaussian distributed random numbers with zero mean and unit variance.

NPT ensemble¶

The NPT ensemble is also called the isothermal-isobaric ensemble.

Berendsen barostat¶

The Berendsen barostat [Berendsen1984] is used with the Berendsen thermostat discussed above. The barostat scales the box and positions as follows:

\[\begin{split}\left( \begin{array}{ccc} a_x^{\rm scaled} & b_x^{\rm scaled} & c_x^{\rm scaled} \\ a_y^{\rm scaled} & b_y^{\rm scaled} & c_y^{\rm scaled} \\ a_z^{\rm scaled} & b_z^{\rm scaled} & c_z^{\rm scaled} \end{array} \right) = \left( \begin{array}{ccc} \mu_{xx} & \mu_{xy} & \mu_{xz} \\ \mu_{yx} & \mu_{yy} & \mu_{yz} \\ \mu_{zx} & \mu_{zy} & \mu_{zz} \\ \end{array} \right) \left( \begin{array}{ccc} a_x & b_x & c_x \\ a_y & b_y & c_y \\ a_z & b_z & c_z \end{array} \right)\end{split}\]

and

\[\begin{split}\left( \begin{array}{c} x^{\rm scaled}_i \\ y^{\rm scaled}_i \\ z^{\rm scaled}_i \end{array} \right) = \left( \begin{array}{ccc} \mu_{xx} & \mu_{xy} & \mu_{xz} \\ \mu_{yx} & \mu_{yy} & \mu_{yz} \\ \mu_{zx} & \mu_{zy} & \mu_{zz} \\ \end{array} \right) \left( \begin{array}{c} x_i \\ y_i \\ z_i \end{array} \right).\end{split}\]

We consider the following three pressure-controlling conditions:

Condition 1: The simulation box is orthogonal and only the hydrostatic pressure (trace of the pressure tensor) is controlled. The simulation box must be periodic in all three directions. The scaling matrix only has nonzero diagonal components and the diagonal components can be written as:

\[\mu_{xx}=\mu_{yy}=\mu_{zz}= 1-\frac{\beta_{\rm hydro} \Delta t}{3 \tau_p} (p^{\rm target}_{\rm hydro} - p^{\rm instant}_{\rm hydro}).\]
Condition 2: The simulation box is orthogonal and the three diagonal pressure components are controlled independently. The simulation box can be periodic or non-periodic in any of the three directions. Pressure is only controlled for periodic directions. The diagonal components of the scaling matrix can be written as:

\[\begin{split}\mu_{xx}= 1-\frac{\beta_{xx} \Delta t}{3 \tau_p} (p^{\rm target}_{xx} - p^{\rm instant}_{xx}) \\ \mu_{yy}= 1-\frac{\beta_{yy} \Delta t}{3 \tau_p} (p^{\rm target}_{yy} - p^{\rm instant}_{yy}) \\ \mu_{zz}= 1-\frac{\beta_{zz} \Delta t}{3 \tau_p} (p^{\rm target}_{zz} - p^{\rm instant}_{zz}).\end{split}\]
Condition 3: The simulation box is triclinic and the 6 nonequivalent pressure components are controlled independently. The simulation box must be periodic in all three directions. The scaling matrix components are:

\[\mu_{\alpha\beta}= 1-\frac{\beta_{\alpha\beta} \Delta t}{3 \tau_p} (p^{\rm target}_{\alpha\beta} - p^{\rm instant}_{\alpha\beta}).\]

The parameter \(\beta_{\alpha\beta}\) is the isothermal compressibility, which is the inverse of the elastic modulus. \(\Delta t\) is the time step and \(\tau_p\) is the pressure coupling time (relaxation time).

Stochastic cell rescaling barostat¶

The Berendsen method does not generate a true NPT ensemble. As an extension of the Berendsen method, the stochastic cell rescaling (SCR) barostat [Bernetti2020], combined with the BDP thermostat, incorporates a proper randomness into the box and position rescaling factor and generates a true NPT ensemble.

In the SCR barostat, the scaling matrix is a sum of the scaling matrix as in the Berendsen barostat and a stochastic one. The stochastic scaling matrix components are

\[\mu^{\rm stochastic}_{\alpha\beta} = \sqrt{ \frac{1}{D_{\rm couple}}} \sqrt{ \frac{\beta_{\alpha\beta} \Delta t}{3\tau_p} \frac{2k_{\rm B} T^{\rm target}}{V} } R_{\alpha\beta}.\]

Here, \(\beta_{\alpha\beta}\), \(\Delta t\), and \(\tau_p\) have the same meanings as in the Berendsen barostat. \(k_{\rm B}\) is Boltzmann’s constant. \(T^{\rm target}\) is the target temperature. \(V\) is the current volume of the system. \(R_{\alpha\beta}\) is a Gaussian random number with zero mean and unit variance. \(D_{\rm couple}\) is the number of directions that are coupled together. It is 3, 1, and 1, respectively, for condition 1, condition 2, and condition 3 as discussed above.