Introduction¶

As a default the DL_POLY_4integration algorithm`s are based on the Velocity :index:`Verlet (VV) scheme, which is both simple and time reversible ⁴. It generates trajectories in the microcanonical (NVE) ensemble in which the total energy (kinetic plus potential) is conserved. If this property drifts or fluctuates excessively in the course of a simulation it indicates that the timestep is too large or the potential cutoffs too small (relative r.m.s. fluctuations in the total energy of \(10^{-5}\) are typical with this algorithm).

The VV algorithm has two stages (VV1 and VV2). At the first stage it requires values of position (\(\underline{r}\)), velocity (\(\underline{v}\)) and force (\(\underline{f}\)) at time \(t\). The first stage is to advance the velocities to \(t+(1/2)\Delta t\) by integration of the force and then to advance the positions to a full step \(t+\Delta t\) using the new half-step velocities:

VV1:

\[\underline{v}(t + {1 \over 2} \Delta t) \leftarrow \underline{v}(t) + {\Delta t \over 2} \; {\underline{f}(t) \over m}~~,\]

where \(m\) is the mass of a site and \(\Delta t\) is the timestep

\[\underline{r}(t + \Delta t) \leftarrow \underline{r}(t) + \Delta t \; \underline{v}(t + {1 \over 2}\Delta t)\]
FF: Between the first and the second stage a recalculation of the force at time \(t+\Delta t\) is required since the positions have changed

\[\underline{f}(t + \Delta t) \leftarrow \underline{f}(t)\]
VV2: In the second stage the half-step velocities are advanced to to a full step using the new force

\[\underline{v}(t + \Delta t) \leftarrow \underline{v}(t + {1 \over 2} \Delta t) + {\Delta t \over 2} \; {\underline{f}(t + \Delta t) \over m}\]

The instantaneous kinetic energy, for example, can then be obtained from the atomic velocities as

(111)¶\[E_{kin}(t) = {1 \over 2} \sum_{1}^{\cal N} m_{i} v_{i}^{2}(t)~~, \label{ekin}\]

and assuming the system has no net momentum the instantaneous temperature is

(112)¶\[{\cal T}(t) = \frac{2}{k_{B} f} E_{kin}(t)~~, \label{tinst}\]

where \(i\) labels particles (that can be free atoms or rigid bodies), \({\cal N}\) the number of particles (free atoms and rigid bodies) in the system, \(k_{B}\) the Boltzmann’s constant and \(f\) the number of degrees of freedom in the system.

(113)¶\[f = 3{\cal N} - 3{\cal N}_{frozen} - 3{\cal N}_{shells} - {\cal N}_{constraints} - 3 - p~~. \label{freedom}\]

Here \({\cal N}_{frozen}\) indicates the number of frozen atoms in the system, \({\cal N}_{shells}\) number of core-shell units and \({\cal N}_{constraints}\) number of bond and PMF constraints. Three degrees of freedom are subtracted for the centre of mass zero net momentum (which we impose) and \(p\) is zero for periodic or three for non-periodic systems, where it accounts for fixing angular momentum about origin (which we impose).

In the case of rigid bodies (see Section Rigid Bodies and Rotational Integration Algorithms) the first part of equation (113)

\[f\prime = 3{\cal N} - 3{\cal N}_{frozen}\]

splits into

\[f\prime = \left( 3{\cal N}^{FP} - 3{\cal N}^{FP}_{frozen} \right) + \left( 3{\cal N}^{RB \mathbf{ (tra)}} - 3{\cal N}^{RB \mathbf{ (tra)}}_{frozen} \right) + \left( 3{\cal N}^{RB \mathbf{ (rot)}} - 3{\cal N}^{RB \mathbf{ (rot)}}_{frozen} \right)\]

\[f\prime = f^{FP} + f^{RB \mathbf{ (tra)}} + f^{RB \mathbf{ (rot)}}~~.\]

Here FP stands for a free particle, i.e. a particle not participating in the constitution of a rigid body, and RB for a rigid body. In general a rigid body has 3 translational (\(\mathbf{ tra}\)) degrees of freedom, corresponding to its centre of mass being allowed to move in the 3 general direction of space, and 3 rotational (\(\mathbf{ rot}\)), corresponding to the RB being allowed to rotate around the 3 general axis in space. It is not far removed to see that for a not fully frozen rigid body one must assign 0 translational degrees of freedom but depending on the “frozenness” of the RB one may assign 1 rotational degrees of freedom when all the frozen sites are in line (i.e. rotation around one axis only) or 3 when just one site is frozen.

The routine nve_0_vv implement the Verlet algorithm in velocity verlet for free particles and calculate the instantaneous temperature. Whereas the routines nve_1_vv implements the same for systems also containing rigid bodies. The conserved quantity is the total energy of the system

\[{\cal H}_{\rm NVE} = U + E_{kin}~~,\]

where \(U\) is the potential energy of the system and \(E_{kin}\) the kinetic energy at time \(t\).

The full selection of integration algorithms within DL_POLY_4is as follows:

`nve_0_vv`	Constant E algorithm
`nve_1_vv`	The same as the above but also incorporating RB integration
`dpd_thermostat`	Constant T algorithm (DPD:cite:shardlow-03a)
`nvt_e0_vv`	Constant \(E_{kin}\) algorithm (Evans:cite:evans-84a)
`nvt_e1_vv`	The same as the above but also incorporating RB integration
`nvt_l0_vv`	Constant T algorithm (Langevin:cite:adelman-76a)
`nvt_l1_vv`	The same as the above but also incorporating RB integration
`nvt_l2_vv`	Constant T algorithm (inhomogeneous Langevin ²⁵)
`nvt_a0_vv`	Constant T algorithm (Andersen ⁶)
`nvt_a1_vv`	The same as the above but also incorporating RB integration
`nvt_b0_vv`	Constant T algorithm (Berendsen ⁹)
`nvt_b1_vv`	The same as the above but also incorporating RB integration
`nvt_h0_vv`	Constant T algorithm (Hoover ⁴³)
`nvt_h1_vv`	The same as the above but also incorporating RB integration
`nvt_g0_vv`	Constant T algorithm (GST ⁵³)
`nvt_g1_vv`	The same as the above but also incorporating RB integration
`npt_l0_vv`	Constant T,P algorithm (Langevin ⁷⁵)
`npt_l1_vv`	The same as the above but also incorporating RB integration
`npt_b0_vv`	Constant T,P algorithm (Berendsen ⁹)
`npt_b1_vv`	The same as the above but also incorporating RB integration
`npt_h0_vv`	Constant T,P algorithm (Hoover ⁴³)
`npt_h1_vv`	The same as the above but also incorporating RB integration
`npt_m0_vv`	Constant T,P algorithm (Martyna-Tuckerman-Klein ⁶⁰)
`npt_m1_vv`	The same as the above but also incorporating RB integration
`npt_l0_vv`	Constant T,\(\underline{\underline{\mathbf{\sigma}}}\) algorithm (Langevin:cite:quigley-04a)
`npt_l1_vv`	The same as the above but also incorporating RB integration
`nst_b0_vv`	Constant T,\(\underline{\underline{\mathbf{\sigma}}}\) algorithm (Berendsen ⁹)
`nst_b1_vv`	The same as the above but also incorporating RB integration
`nst_h0_vv`	Constant T,\(\underline{\underline{\mathbf{\sigma}}}\) algorithm (Hoover ⁴³)
`nst_h1_vv`	The same as the above but also incorporating RBs integration
`nst_m0_vv`	Constant T,\(\underline{\underline{\mathbf{\sigma}}}\) algorithm (Martyna-Tuckerman-Klein ⁶⁰)
`nst_m0_vv`	The same as the above but also incorporating RB integration

It is worth noting that the last four ensembles are also optionally available in an extended from to constant normal pressure and constant surface area, NP\(_{n}\)AT, or constant surface tension, NP\(_{n}\gamma\)T ⁴⁴.