The DL_POLY_2 Force Field

The force field is the set of functions needed to define the interactions in a molecular system. These may have a wide variety of analytical forms, with some basis in chemical physics, which must be parameterised to give the correct energy and forces. A huge variety of forms is possible and for this reason the DL_POLY_2 force field is designed to be adaptable. While it is not supplied with its own force field parameters, many of the functions familiar to GROMOS, [3] Dreiding [8] and AMBER [4] users have been coded in the package, as well as less familiar forms. In addition DL_POLY_2 retains the possibility of the user defining additional potentials.

In DL_POLY_2 the total configuration energy of a molecular system may be written as:

$$U(\mbox{$\underline{r}$}_{1},\mbox{$\underline{r}$}_{2},\ldots,\mbox{$\underline{r}$}_{N}) = \sum_{i_{bond}=1}^{N_{bond}} U_{bond}(i_{bond},\mbox{$\underline{r}$}_{a},\mbox{$\underline{r}$}_{b})$$

$$+\sum_{i_{angle}=1}^{N_{angle}} U_{angle}(i_{angle},\mbox{$\underline{r}$}_{a},\mbox{$\underline{r}$}_{b},\mbox{$\underline{r}$}_{c})$$

$$+\sum_{i_{dihed}=1}^{N_{dihed}} U_{dihed}(i_{dihed},\mbox{$\under... ...box{$\underline{r}$}_{b},\mbox{$\underline{r}$}_{c},\mbox{$\underline{r}$}_{d})$$

$$+\sum_{i_{inv}=1}^{N_{inv}} U_{inv}(i_{inv},\mbox{$\underline{r}$... ...box{$\underline{r}$}_{b},\mbox{$\underline{r}$}_{c},\mbox{$\underline{r}$}_{d})$$

$$+\sum_{i=1}^{N-1}\sum_{j>i}^{N} U_{pair}(i,j,\vert\mbox{$\underline{r}$}_{i}-\mbox{$\underline{r}$}_{j}\vert)$$

$$+\sum_{i=1}^{N-2}\sum_{j>i}^{N-1}\sum_{k>j}^{N} U_{3\_body}(i,j,k,\mbox{$\underline{r}$}_{i},\mbox{$\underline{r}$}_{j},\mbox{$\underline{r}$}_{k})$$

$$+\sum_{i=1}^{N-3}\sum_{j>i}^{N-2}\sum_{k>j}^{N-1}\sum_{n>k}^{N} U... ...box{$\underline{r}$}_{j},\mbox{$\underline{r}$}_{k},\mbox{$\underline{r}$}_{n})$$

$$+\sum_{i=1}^{N}U_{extn}(i,\mbox{$\underline{r}$}_{i},\mbox{$\underline{v}$}_{i})$$ (2.1)

where $U_{bond},~U_{angle},~U_{dihed},~U_{inv}$, $U_{pair}$, $U_{3\_body}$ and $U_{4\_body}$ are empirical interaction functions representing chemical bonds, valence angles, dihedral angles, inversion angles, pair-body, three-body and four-body forces respectively. The first four are regarded by DL_POLY_2 as intra-molecular interactions and the next three as inter-molecular interactions. The final term $U_{extn}$ represents an external field potential. The position vectors $\mbox{$\underline{r}$}_{a},\mbox{$\underline{r}$}_{b},\mbox{$\underline{r}$}_{c}$ and $\mbox{$\underline{r}$}_{d}$ refer to the positions of the atoms specifically involved in a given interaction. (Almost universally, it is the differences in position that determine the interaction.) The numbers $N_{bond},~N_{angle}$, $N_{dihed}$ and $N_{inv}$ refer to the total numbers of these respective interactions present in the simulated system, and the indices $i_{bond},~i_{angle},~i_{inv}$ and $i_{dihed}$ uniquely specify an individual interaction of each type. It is important to note that there is no global specification of the intramolecular interactions in DL_POLY_2 - all bonds, valence angles and dihedrals must be individually cited.

The indices $i$, $j$ (and $k$, $n$) appearing in the pair-body (and three or four-body) terms indicate the atoms involved in the interaction. There is normally a very large number of these and they are therefore specified according to atom types rather than indices. In DL_POLY_2 it is assumed that the pair-body terms arise from van der Waals and/or electrostatic (Coulombic) forces. The former are regarded as short ranged interactions and the latter as long ranged. Long range forces require special techniques to evaluate accurately (see section 2.4.) In DL_POLY_2 the three-body terms are restricted to valence angle and H-bond forms. The nonbonded, three-body and four-body interactions are globally specified according to the types of atoms involved. DL_POLY_2 also has the ability to handle metals via density dependent functions (see below). Though essentially many-body potentials they are handled in DL_POLY_2 as special forms of pair potential.

In DL_POLY_2 the intramolecular bonded terms are handled using bookkeeping arrays, which specify the atoms involved in a particular interaction and point to the appropriate arrays of parameters that define the potential. The calculation of bonded forces therefore follows the simple scheme:

1. Every atom in the simulated system is assigned a unique index number from $1$ to $N$;
2. Every intramolecular bonded term $U_{type}$ in the system has a unique index number $i_{type}$: from $1$ to $N_{type}$ where $type$ represents a bond, angle or dihedral.
3. A pointer array $key_{type}(n_{type},i_{type})$ carries the indices of the specific atoms involved in the potential term labelled $i_{type}$. The dimension $n_{type}$ will be $2,~3$ or $4$, if the term represents a bond, valence angle, dihedral/inversion.
4. The array $key_{type}(n_{type},i_{type})$ is used to identify the atoms in a bonded term and the appropriate form of interaction and thus to calculate the energy and forces.

DL_POLY_2 calculates the nonbonded pair interactions using a Verlet neighbour list [12] which is reconstructed at intervals during the simulation. This list records the indices of all `secondary' atoms within a certain radius of each `primary' atom; the radius being the cut-off radius ($r_{cut}$) normally applied to the nonbonded potential function, plus an additional increment ( $\Delta r_{cut}$). The neighbour list removes the need to scan over all atoms in the simulation at every timestep. The larger radius ( $r_{cut}+\Delta r_{cut}$) means the same list can be used for several timesteps without requiring an update. The frequency at which the list must be updated depends on the thickness of the region $\Delta r_{cut}$. DL_POLY_2 has two methods for constructing the neighbour list: the first is based on the Brode-Ahlrichs scheme [22] and is used when $r_{cut}$ is large in comparison with the simulation cell; the second uses the link-cell algorithm [23] when $r_{cut}$ is relatively small. The potential energy and forces arising from the nonbonded interactions are calculated using interpolation tables.

A complication in the construction of the Verlet neighbour list for macromolecules is the concept of excluded atoms, which arises from the need to exclude certain atom pairs from the overall list. Which atom pairs need to be excluded is dependent on the precise nature of the force field model, but as a minimum atom pairs linked via extensible bonds or constraints and atoms (grouped in pairs) linked via valence angles are probable candidates. The assumption behind this requirement is that atoms that are formally bonded in a chemical sense, should not participate in nonbonded interactions. (However this is not a universal requirement of all force fields.) The same considerations are needed in dealing with charged excluded atoms. DL_POLY_2 has several subroutines available for constructing the Verlet neighbour list, while taking care of the excluded atoms (see chapters 3 and 8 for further information.)

Three- and four-body nonbonded forces are assumed to be short ranged and therefore calculated using the link-cell algorithm [23]. They ignore the possibility of there being any excluded interactions involving the atoms concerned.

Throughout this section the description of the force field assumes the simulated system is described as an assembly of atoms. This is for convenience only and readers should understand that DL_POLY_2 does recognise molecular entities, defined either through constraint bonds or rigid bodies. In the case of rigid bodies, the atomic forces are resolved into molecular forces and torques. These matters are discussed in greater detail later in sections 2.5.1 and 2.5.6).

Note that the subroutines mentioned in the subsections of this chapter are described in greater detail in chapter 8 of the Reference Manual.