Long Range Van der Waals interactions#

Dispersion correction#

In this section, we derive long-range corrections due to the use of a cut-off for Lennard-Jones or Buckingham interactions. We assume that the cut-off is so long that the repulsion term can safely be neglected, and therefore only the dispersion term is taken into account. Due to the nature of the dispersion interaction (we are truncating a potential proportional to r6), energy and pressure corrections are both negative. While the energy correction is usually small, it may be important for free energy calculations where differences between two different Hamiltonians are considered. In contrast, the pressure correction is very large and can not be neglected under any circumstances where a correct pressure is required, especially for any NPT simulations. Although it is, in principle, possible to parameterize a force field such that the pressure is close to the desired experimental value without correction, such a method makes the parameterization dependent on the cut-off and is therefore undesirable.

Energy#

The long-range contribution of the dispersion interaction to the virial can be derived analytically, if we assume a homogeneous system beyond the cut-off distance rc. The dispersion energy between two particles is written as:

(303)#V(rij) = C6rij6

and the corresponding force is:

(304)#Fij = 6C6rij8rij

In a periodic system it is not easy to calculate the full potentials, so usually a cut-off is applied, which can be abrupt or smooth. We will call the potential and force with cut-off Vc and Fc. The long-range contribution to the dispersion energy in a system with N particles and particle density ρ = N/V is:

(305)#Vlr = 12Nρ04πr2g(r)(V(r)Vc(r))dr

We will integrate this for the shift function, which is the most general form of van der Waals interaction available in GROMACS. The shift function has a constant difference S from 0 to r1 and is 0 beyond the cut-off distance rc. We can integrate (305), assuming that the density in the sphere within r1 is equal to the global density and the radial distribution function g(r) is 1 beyond r1:

(306)#Vlr=12N(ρ0r14πr2g(r)C6Sdr+ρr1rc4πr2(V(r)Vc(r))dr+ρrc4πr2V(r)dr)=12N((43πρr131)C6S+ρr1rc4πr2(V(r)Vc(r))dr43πNρC6rc3)

where the term 1 corrects for the self-interaction. For a plain cut-off we only need to assume that g(r) is 1 beyond rc and the correction reduces to 108:

(307)#Vlr=23πNρC6rc3

If we consider, for example, a box of pure water, simulated with a cut-off of 0.9 nm and a density of 1 g cm3 this correction is 0.75 kJ mol1 per molecule.

For a homogeneous mixture we need to define an average dispersion constant:

(308)#C6=2N(N1)iNj>iNC6(i,j)

In GROMACS, excluded pairs of atoms do not contribute to the average.

In the case of inhomogeneous simulation systems, e.g. a system with a lipid interface, the energy correction can be applied if C6 for both components is comparable.

Virial and pressure#

The scalar virial of the system due to the dispersion interaction between two particles i and j is given by:

(309)#Ξ = 12rijFij = 3C6rij6

The pressure is given by:

(310)#P = 23V(EkinΞ)

The long-range correction to the virial is given by:

(311)#Ξlr = 12Nρ04πr2g(r)(ΞΞc)dr

We can again integrate the long-range contribution to the virial assuming g(r) is 1 beyond r1:

(312)#Ξlr=12Nρ(r1rc4πr2(ΞΞc)dr+rc4πr23C6rij6dr)=12Nρ(r1rc4πr2(ΞΞc)dr+4πC6rc3)

For a plain cut-off the correction to the pressure is 108:

(313)#Plr = 43πC6ρ2rc3

Using the same example of a water box, the correction to the virial is 0.75 kJ mol1 per molecule, the corresponding correction to the pressure for SPC water is approximately 280 bar.

For homogeneous mixtures, we can again use the average dispersion constant C6 ((308)):

(314)#Plr = 43πC6ρ2rc3

For inhomogeneous systems, (314) can be applied under the same restriction as holds for the energy (see sec. Energy).

Lennard-Jones PME#

In order to treat systems, using Lennard-Jones potentials, that are non-homogeneous outside of the cut-off distance, we can instead use the Particle-mesh Ewald method as discussed for electrostatics above. In this case the modified Ewald equations become

(315)#V=Vdir+Vrec+V0Vdir=12i,jNnxnynzC6ijg(βrij,n)rij,n6
(316)#Vrec=π32β32Vmxmymzf(π|m|/β)×i,jNC6ijexp[2πim(rirj)]V0=β612iNC6ii

where m=(mx,my,mz), β is the parameter determining the weight between direct and reciprocal space, and C6ij is the combined dispersion parameter for particle i and j. The star indicates that terms with i=j should be omitted when ((nx,ny,nz)=(0,0,0)), and rij,n is the real distance between the particles. Following the derivation by Essmann 15, the functions f and g introduced above are defined as

(317)#f(x)=1/3[(12x2)exp(x2)+2x3πerfc(x)]g(x)=exp(x2)(1+x2+x42).

The above methodology works fine as long as the dispersion parameters can be combined geometrically ((145)) in the same way as the charges for electrostatics

(318)#C6,geomij=(C6iiC6jj)1/2

For Lorentz-Berthelot combination rules ((146)), the reciprocal part of this sum has to be calculated seven times due to the splitting of the dispersion parameter according to

(319)#C6,LBij=(σi+σj)6=n=06Pnσinσj(6n),

for Pn the Pascal triangle coefficients. This introduces a non-negligible cost to the reciprocal part, requiring seven separate FFTs, and therefore this has been the limiting factor in previous attempts to implement LJ-PME. A solution to this problem is to use geometrical combination rules in order to calculate an approximate interaction parameter for the reciprocal part of the potential, yielding a total interaction of

(320)#V(r<rc)=C6dirg(βr)r6Direct space+C6,geomrecip[1g(βr)]r6Reciprocal space=C6,geomrecipr6+(C6dirC6,geomrecip)g(βr)r6V(r>rc)=C6,geomrecip[1g(βr)]r6Reciprocal space.

This will preserve a well-defined Hamiltonian and significantly increase the performance of the simulations. The approximation does introduce some errors, but since the difference is located in the interactions calculated in reciprocal space, the effect will be very small compared to the total interaction energy. In a simulation of a lipid bilayer, using a cut-off of 1.0 nm, the relative error in total dispersion energy was below 0.5%. A more thorough discussion of this can be found in 109.

In GROMACS we now perform the proper calculation of this interaction by subtracting, from the direct-space interactions, the contribution made by the approximate potential that is used in the reciprocal part

(321)#Vdir=C6dirr6C6recip[1g(βr)]r6.

This potential will reduce to the expression in (315) when C6dir=C6recip, and the total interaction is given by

(322)#V(r<rc)=C6dirr6C6recip[1g(βr)]r6Direct space+C6recip[1g(βr)]r6Reciprocal space=C6dirr6
(323)#V(r>rc)=C6recip[1g(βr)]r6.

For the case when C6dirC6recip this will retain an unmodified LJ force up to the cut-off, and the error is an order of magnitude smaller than in simulations where the direct-space interactions do not account for the approximation used in reciprocal space. When using a VdW interaction modifier of potential-shift, the constant

(324)#(C6dir+C6recip[1g(βrc)])rc6

is added to (322) in order to ensure that the potential is continuous at the cutoff. Note that, in the same way as (321), this degenerates into the expected C6g(βrc)rc6 when C6dir=C6recip. In addition to this, a long-range dispersion correction can be applied to correct for the approximation using a combination rule in reciprocal space. This correction assumes, as for the cut-off LJ potential, a uniform particle distribution. But since the error of the combination rule approximation is very small this long-range correction is not necessary in most cases. Also note that this homogenous correction does not correct the surface tension, which is an inhomogeneous property.

Using LJ-PME#

As an example for using Particle-mesh Ewald summation for Lennard-Jones interactions in GROMACS, specify the following lines in your mdp file:

vdwtype          = PME
rvdw             = 0.9
vdw-modifier     = Potential-Shift
rlist            = 0.9
rcoulomb         = 0.9
fourierspacing   = 0.12
pme-order        = 4
ewald-rtol-lj    = 0.001
lj-pme-comb-rule = geometric

The same Fourier grid and interpolation order are used if both LJ-PME and electrostatic PME are active, so the settings for fourierspacing and pme-order are common to both. ewald-rtol-lj controls the splitting between direct and reciprocal space in the same way as ewald-rtol. In addition to this, the combination rule to be used in reciprocal space is determined by lj-pme-comb-rule. If the current force field uses Lorentz-Berthelot combination rules, it is possible to set lj-pme-comb-rule = geometric in order to gain a significant increase in performance for a small loss in accuracy. The details of this approximation can be found in the section above.

Note that the use of a complete long-range dispersion correction means that as with Coulomb PME, rvdw is now a free parameter in the method, rather than being necessarily restricted by the force-field parameterization scheme. Thus it is now possible to optimize the cutoff, spacing, order and tolerance terms for accuracy and best performance.

Naturally, the use of LJ-PME rather than LJ cut-off adds computation and communication done for the reciprocal-space part, so for best performance in balancing the load of parallel simulations using PME-only ranks, more such ranks should be used. It may be possible to improve upon the automatic load-balancing used by mdrun.