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
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
and the corresponding force is:
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
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
where the term
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 cm
For a homogeneous mixture we need to define an average dispersion constant:
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
Virial and pressure#
The scalar virial of the system due to the dispersion interaction
between two particles
The pressure is given by:
The long-range correction to the virial is given by:
We can again integrate the long-range contribution to the virial
assuming
For a plain cut-off the correction to the pressure is 108:
Using the same example of a water box, the correction to the virial is
0.75 kJ mol
For homogeneous mixtures, we can again use the average dispersion
constant
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
where
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
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
for
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
This potential will reduce to the expression in
(315) when
For the case when
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
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.