Non-bonded interactions

Non-bonded interactions in GROMACS are pair-additive:

(1)V(r1,rN)=i<jVij(rij);
(2)Fi=jdVij(rij)drijrijrij

Since the potential only depends on the scalar distance, interactions will be centro-symmetric, i.e. the vectorial partial force on particle i from the pairwise interaction Vij(rij) has the opposite direction of the partial force on particle j. For efficiency reasons, interactions are calculated by loops over interactions and updating both partial forces rather than summing one complete nonbonded force at a time. The non-bonded interactions contain a repulsion term, a dispersion term, and a Coulomb term. The repulsion and dispersion term are combined in either the Lennard-Jones (or 6-12 interaction), or the Buckingham (or exp-6 potential). In addition, (partially) charged atoms act through the Coulomb term.

The Lennard-Jones interaction

The Lennard-Jones potential VLJ between two atoms equals:

(3)VLJ(rij)=C(12)ijrij12C(6)ijrij6

See also Fig. 16 The parameters C(12)ij and C(6)ij depend on pairs of atom types; consequently they are taken from a matrix of LJ-parameters. In the Verlet cut-off scheme, the potential is shifted by a constant such that it is zero at the cut-off distance.

../../_images/f-lj.png

Fig. 16 The Lennard-Jones interaction.

The force derived from this potential is:

(4)Fi(rij)=(12 C(12)ijrij136 C(6)ijrij7)rijrij

The LJ potential may also be written in the following form:

(5)VLJ(rij)=4ϵij((σijrij)12(σijrij)6)

In constructing the parameter matrix for the non-bonded LJ-parameters, two types of combination rules can be used within GROMACS, only geometric averages (type 1 in the input section of the force-field file):

(6)C(6)ij=(C(6)iiC(6)jj)1/2C(12)ij=(C(12)iiC(12)jj)1/2

or, alternatively the Lorentz-Berthelot rules can be used. An arithmetic average is used to calculate σij, while a geometric average is used to calculate ϵij (type 2):

(7)σij=12(σii+σjj)ϵij=(ϵiiϵjj)1/2

finally an geometric average for both parameters can be used (type 3):

(8)σij=(σiiσjj)1/2ϵij=(ϵiiϵjj)1/2

This last rule is used by the OPLS force field.

Buckingham potential

The Buckingham potential has a more flexible and realistic repulsion term than the Lennard-Jones interaction, but is also more expensive to compute. The potential form is:

(9)Vbh(rij)=Aijexp(Bijrij)Cijrij6
../../_images/f-bham.png

Fig. 17 The Buckingham interaction.

See also Fig. 17. The force derived from this is:

(10)Fi(rij)=[AijBijexp(Bijrij)6Cijrij7]rijrij

Coulomb interaction

The Coulomb interaction between two charge particles is given by:

(11)Vc(rij)=fqiqjεrrij

See also Fig. 18, where f=14πε0=138.935458 (see chapter Definitions and Units)

../../_images/vcrf.png

Fig. 18 The Coulomb interaction (for particles with equal signed charge) with and without reaction field. In the latter case εr was 1, εrf was 78, and rc was 0.9 nm. The dot-dashed line is the same as the dashed line, except for a constant.

The force derived from this potential is:

(12)Fi(rij)=fqiqjεrrij2rijrij

A plain Coulomb interaction should only be used without cut-off or when all pairs fall within the cut-off, since there is an abrupt, large change in the force at the cut-off. In case you do want to use a cut-off, the potential can be shifted by a constant to make the potential the integral of the force. With the group cut-off scheme, this shift is only applied to non-excluded pairs. With the Verlet cut-off scheme, the shift is also applied to excluded pairs and self interactions, which makes the potential equivalent to a reaction field with εrf=1 (see below).

In GROMACS the relative dielectric constant εr may be set in the in the input for grompp.

Coulomb interaction with reaction field

The Coulomb interaction can be modified for homogeneous systems by assuming a constant dielectric environment beyond the cut-off rc with a dielectric constant of εrf. The interaction then reads:

(13)Vcrf = fqiqjεrrij[1+εrfεr2εrf+εrrij3r3c]fqiqjεrrc3εrf2εrf+εr

in which the constant expression on the right makes the potential zero at the cut-off rc. For charged cut-off spheres this corresponds to neutralization with a homogeneous background charge. We can rewrite (13) for simplicity as

(14)Vcrf = fqiqjεr[1rij+krf rij2crf]

with

(15)krf=1r3cεrfεr(2εrf+εr)
(16)crf=1rc+krfr2c = 1rc3εrf(2εrf+εr)

For large εrf the krf goes to r3c/2, while for εrf = εr the correction vanishes. In Fig. 18 the modified interaction is plotted, and it is clear that the derivative with respect to rij (= -force) goes to zero at the cut-off distance. The force derived from this potential reads:

(17)Fi(rij)=fqiqjεr[1rij22krfrij]rijrij

The reaction-field correction should also be applied to all excluded atoms pairs, including self pairs, in which case the normal Coulomb term in (13) and (17) is absent.

Modified non-bonded interactions

In GROMACS, the non-bonded potentials can be modified by a shift function, also called a force-switch function, since it switches the force to zero at the cut-off. The purpose of this is to replace the truncated forces by forces that are continuous and have continuous derivatives at the cut-off radius. With such forces the time integration produces smaller errors. But note that for Lennard-Jones interactions these errors are usually smaller than other errors, such as integration errors at the repulsive part of the potential. For Coulomb interactions we advise against using a shifted potential and for use of a reaction field or a proper long-range method such as PME.

There is no fundamental difference between a switch function (which multiplies the potential with a function) and a shift function (which adds a function to the force or potential) 72. The switch function is a special case of the shift function, which we apply to the force function F(r), related to the electrostatic or van der Waals force acting on particle i by particle j as:

(18)Fi=cF(rij)rijrij

For pure Coulomb or Lennard-Jones interactions F(r)=Fα(r)=αr(α+1). The switched force Fs(r) can generally be written as:

(19)Fs(r) = Fα(r)r<r1Fs(r) = Fα(r)+S(r)r1r<rcFs(r) = 0rcr

When r1=0 this is a traditional shift function, otherwise it acts as a switch function. The corresponding shifted potential function then reads:

(20)Vs(r)=r Fs(x)dx

The GROMACS force switch function SF(r) should be smooth at the boundaries, therefore the following boundary conditions are imposed on the switch function:

(21)SF(r1)=0SF(r1)=0SF(rc)=Fα(rc)SF(rc)=Fα(rc)

A 3rd degree polynomial of the form

(22)SF(r)=A(rr1)2+B(rr1)3

fulfills these requirements. The constants A and B are given by the boundary condition at rc:

(23)A = α(α+4)rc  (α+1)r1rα+2c (rcr1)2B = α(α+3)rc  (α+1)r1rα+2c (rcr1)3

Thus the total force function is:

(24)Fs(r)=αrα+1+A(rr1)2+B(rr1)3

and the potential function reads:

(25)Vs(r)=1rαA3(rr1)3B4(rr1)4C

where

(26)C=1rαcA3(rcr1)3B4(rcr1)4

The GROMACS potential-switch function SV(r) scales the potential between r1 and rc, and has similar boundary conditions, intended to produce smoothly-varying potential and forces:

(27)SV(r1)=1SV(r1)=0SV

The fifth-degree polynomial that has these properties is

(28)S_V(r; r_1, r_c) = \frac{1 - 10(r-r_1)^3(r_c-r_1)^2 + 15(r-r_1)^4(r_c-r_1) - 6(r-r_1)}{(r_c-r_1)^5}

This implementation is found in several other simulation packages,7375 but differs from that in CHARMM.76 Switching the potential leads to artificially large forces in the switching region, therefore it is not recommended to switch Coulomb interactions using this function,72 but switching Lennard-Jones interactions using this function produces acceptable results.

Modified short-range interactions with Ewald summation

When Ewald summation or particle-mesh Ewald is used to calculate the long-range interactions, the short-range Coulomb potential must also be modified. Here the potential is switched to (nearly) zero at the cut-off, instead of the force. In this case the short range potential is given by:

(29)V(r) = f \frac{\mbox{erfc}(\beta r_{ij})}{r_{ij}} q_i q_j,

where \beta is a parameter that determines the relative weight between the direct space sum and the reciprocal space sum and erfc(x) is the complementary error function. For further details on long-range electrostatics, see sec. Long Range Electrostatics.