docs/reference-manual/functions/long-range-vdw.rst

   1 Long Range Van der Waals interactions
   2 -------------------------------------
   3
   4 Dispersion correction
   5 ~~~~~~~~~~~~~~~~~~~~~
   6
   7 In this section, we derive long-range corrections due to the use of a
   8 cut-off for Lennard-Jones or Buckingham interactions. We assume that the
   9 cut-off is so long that the repulsion term can safely be neglected, and
  10 therefore only the dispersion term is taken into account. Due to the
  11 nature of the dispersion interaction (we are truncating a potential
  12 proportional to :math:`-r^{-6}`), energy and pressure corrections are
  13 both negative. While the energy correction is usually small, it may be
  14 important for free energy calculations where differences between two
  15 different Hamiltonians are considered. In contrast, the pressure
  16 correction is very large and can not be neglected under any
  17 circumstances where a correct pressure is required, especially for any
  18 NPT simulations. Although it is, in principle, possible to parameterize
  19 a force field such that the pressure is close to the desired
  20 experimental value without correction, such a method makes the
  21 parameterization dependent on the cut-off and is therefore undesirable.
  22
  23 .. _ecorr:
  24
  25 Energy
  26 ^^^^^^
  27
  28 The long-range contribution of the dispersion interaction to the virial
  29 can be derived analytically, if we assume a homogeneous system beyond
  30 the cut-off distance :math:`r_c`. The dispersion energy between two
  31 particles is written as:
  32
  33 .. math:: V({r_{ij}}) ~=~- C_6\,{r_{ij}}^{-6}
  34           :label: eqnlrljEdisp
  35
  36 and the corresponding force is:
  37
  38 .. math:: \mathbf{F}_{ij} ~=~- 6\,C_6\,r_{ij}^{-8}\mathbf{r}_{ij}
  39           :label: eqnlrljFdisp
  40
  41 In a periodic system it is not easy to calculate the full potentials,
  42 so usually a cut-off is applied, which can be abrupt or smooth. We will
  43 call the potential and force with cut-off :math:`V_c` and
  44 :math:`\mathbf{F}_c`. The long-range contribution to the
  45 dispersion energy in a system with :math:`N` particles and particle
  46 density :math:`\rho` = :math:`N/V` is:
  47
  48 .. math:: V_{lr} ~=~ {\frac{1}{2}}N \rho\int_0^{\infty}   4\pi r^2 g(r) \left( V(r) -V_c(r) \right) {{{\rm d}r}}
  49           :label: eqnenercorr
  50
  51 We will integrate this for the shift function, which is the most
  52 general form of van der Waals interaction available in |Gromacs|. The
  53 shift function has a constant difference :math:`S` from 0 to :math:`r_1`
  54 and is 0 beyond the cut-off distance :math:`r_c`. We can integrate
  55 :eq:`eqn. %s <eqnenercorr>`, assuming that the density in the sphere within
  56 :math:`r_1` is equal to the global density and the radial distribution
  57 function :math:`g(r)` is 1 beyond :math:`r_1`:
  58
  59 .. math:: \begin{aligned}
  60           V_{lr}  &=& {\frac{1}{2}}N \left(
  61             \rho\int_0^{r_1}  4\pi r^2 g(r) \, C_6 \,S\,{{{\rm d}r}}
  62           + \rho\int_{r_1}^{r_c}  4\pi r^2 \left( V(r) -V_c(r) \right) {{{\rm d}r}}
  63           + \rho\int_{r_c}^{\infty}  4\pi r^2 V(r) \, {{{\rm d}r}}
  64           \right) \\
  65           & = & {\frac{1}{2}}N \left(\left(\frac{4}{3}\pi \rho r_1^{3} - 1\right) C_6 \,S
  66           + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {{{\rm d}r}}
  67           -\frac{4}{3} \pi N \rho\, C_6\,r_c^{-3}
  68           \right)\end{aligned}
  69           :label: eqnlrljshift
  70
  71 where the term :math:`-1` corrects for the self-interaction. For a
  72 plain cut-off we only need to assume that :math:`g(r)` is 1 beyond
  73 :math:`r_c` and the correction reduces to \ :ref:`108 <refAllen87>`:
  74
  75 .. math:: \begin{aligned}
  76           V_{lr} & = & -\frac{2}{3} \pi N \rho\, C_6\,r_c^{-3}\end{aligned}
  77           :label: eqnlrljcorrreduced
  78
  79 If we consider, for example, a box of pure water, simulated with a
  80 cut-off of 0.9 nm and a density of 1 g cm\ :math:`^{-3}` this correction
  81 is :math:`-0.75` kJ mol\ :math:`^{-1}` per molecule.
  82
  83 For a homogeneous mixture we need to define an *average dispersion constant*:
  84
  85 .. math:: {\left< C_6 \right>}= \frac{2}{N(N-1)}\sum_i^N\sum_{j>i}^N C_6(i,j)\\
  86           :label: eqnavcsix
  87
  88 In |Gromacs|, excluded pairs of atoms do not contribute to the average.
  89
  90 In the case of inhomogeneous simulation systems, *e.g.* a system with a
  91 lipid interface, the energy correction can be applied if
  92 :math:`{\left< C_6 \right>}` for both components is comparable.
  93
  94 .. _virial:
  95
  96 Virial and pressure
  97 ^^^^^^^^^^^^^^^^^^^
  98
  99 The scalar virial of the system due to the dispersion interaction
 100 between two particles :math:`i` and :math:`j` is given by:
 101
 102 .. math:: \Xi~=~-{\frac{1}{2}} \mathbf{r}_{ij} \cdot \mathbf{F}_{ij} ~=~ 3\,C_6\,r_{ij}^{-6}
 103           :label: eqnlrljdispvirial
 104
 105 The pressure is given by:
 106
 107 .. math:: P~=~\frac{2}{3\,V}\left(E_{kin} - \Xi\right)
 108           :label: eqnlrljpressure
 109
 110 The long-range correction to the virial is given by:
 111
 112 .. math:: \Xi_{lr} ~=~ {\frac{1}{2}}N \rho \int_0^{\infty} 4\pi r^2 g(r) (\Xi -\Xi_c) \,{{\rm d}r}
 113           :label: eqnlrljcorrvirial
 114
 115 We can again integrate the long-range contribution to the virial
 116 assuming :math:`g(r)` is 1 beyond :math:`r_1`:
 117
 118 .. math:: \begin{aligned}
 119           \Xi_{lr}&=&   {\frac{1}{2}}N \rho \left(
 120               \int_{r_1}^{r_c}  4 \pi r^2 (\Xi -\Xi_c)  \,{{\rm d}r}+ \int_{r_c}^{\infty} 4 \pi r^2 3\,C_6\,{r_{ij}}^{-6}\,  {{\rm d}r}\right)  \nonumber\\
 121                   &=&     {\frac{1}{2}}N \rho \left(
 122               \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \, {{\rm d}r}+ 4 \pi C_6 \, r_c^{-3} \right)\end{aligned}
 123           :label: eqnlrljvirialcontrib
 124
 125 For a plain cut-off the correction to the pressure
 126 is \ :ref:`108 <refAllen87>`:
 127
 128 .. math:: P_{lr}~=~-\frac{4}{3} \pi C_6\, \rho^2 r_c^{-3}
 129           :label: eqnlrljpressurecorr
 130
 131 Using the same example of a water box, the correction to the virial is
 132 0.75 kJ mol\ :math:`^{-1}` per molecule, the corresponding correction to
 133 the pressure for SPC water is approximately :math:`-280` bar.
 134
 135 For homogeneous mixtures, we can again use the average dispersion
 136 constant :math:`{\left< C_6 \right>}` (:eq:`eqn. %s <eqnavcsix>`):
 137
 138 .. math:: P_{lr}~=~-\frac{4}{3} \pi {\left< C_6 \right>}\rho^2 r_c^{-3}
 139           :label: eqnpcorr
 140
 141 For inhomogeneous systems, :eq:`eqn. %s <eqnpcorr>` can be applied under the
 142 same restriction as holds for the energy (see sec. :ref:`ecorr`).
 143
 144 Lennard-Jones PME
 145 ~~~~~~~~~~~~~~~~~
 146
 147 In order to treat systems, using Lennard-Jones potentials, that are
 148 non-homogeneous outside of the cut-off distance, we can instead use the
 149 Particle-mesh Ewald method as discussed for electrostatics above. In
 150 this case the modified Ewald equations become
 151
 152 .. math:: \begin{aligned}
 153           V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex]
 154           V_{\mathrm{dir}} &=& -\frac{1}{2} \sum_{i,j}^{N}
 155           \sum_{n_x}\sum_{n_y}
 156           \sum_{n_{z}*} \frac{C^{ij}_6 g(\beta {r}_{ij,{\bf n}})}{{r_{ij,{\bf n}}}^6}
 157           \end{aligned}
 158           :label: eqnljpmerealspace
 159
 160 .. math:: \begin{aligned}
 161           V_{\mathrm{rec}} &=& \frac{{\pi}^{\frac{3}{2}} \beta^{3}}{2V} \sum_{m_x}\sum_{m_y}\sum_{m_{z}*}
 162           f(\pi | {\mathbf m} | /\beta) \times \sum_{i,j}^{N} C^{ij}_6 {\mathrm{exp}}\left[-2\pi i {\bf m}\cdot({\bf r_i}-{\bf r_j})\right] \\[0.5ex]
 163           V_{0} &=& -\frac{\beta^{6}}{12}\sum_{i}^{N} C^{ii}_6\end{aligned}
 164           :label: eqnljpmerealspace2
 165
 166 where :math:`{\bf m}=(m_x,m_y,m_z)`, :math:`\beta` is the parameter
 167 determining the weight between direct and reciprocal space, and
 168 :math:`{C^{ij}_6}` is the combined dispersion parameter for particle
 169 :math:`i` and :math:`j`. The star indicates that terms with
 170 :math:`i = j` should be omitted when :math:`((n_x,n_y,n_z)=(0,0,0))`,
 171 and :math:`{\bf r}_{ij,{\bf n}}` is the real distance between the
 172 particles. Following the derivation by
 173 Essmann \ :ref:`15 <refEssmann95>`, the functions :math:`f` and :math:`g`
 174 introduced above are defined as
 175
 176 .. math:: \begin{aligned}
 177           f(x)&=&1/3\left[(1-2x^2){\mathrm{exp}}(-x^2) + 2{x^3}\sqrt{\pi}\,{\mathrm{erfc}}(x) \right] \\
 178           g(x)&=&{\mathrm{exp}}(-x^2)(1+x^2+\frac{x^4}{2}).\end{aligned}
 179           :label: eqnljpmerealdistance
 180
 181 The above methodology works fine as long as the dispersion parameters
 182 can be combined geometrically (:eq:`eqn. %s <eqncomb>`) in the same way as the
 183 charges for electrostatics
 184
 185 .. math:: C^{ij}_{6,\mathrm{geom}} = \left(C^{ii}_6 \, C^{jj}_6\right)^{1/2}
 186           :label: eqnljpmegeom
 187
 188 For Lorentz-Berthelot combination rules (:eq:`eqn. %s <eqnlorentzberthelot>`),
 189 the reciprocal part of this sum has to be calculated seven times due to
 190 the splitting of the dispersion parameter according to
 191
 192 .. math:: C^{ij}_{6,\mathrm{L-B}} = (\sigma_i+\sigma_j)^6=\sum_{n=0}^{6} P_{n}\sigma_{i}^{n}\sigma_{j}^{(6-n)},
 193           :label: eqnljpmelorenztberthelot
 194
 195 for :math:`P_{n}` the Pascal triangle coefficients. This introduces a
 196 non-negligible cost to the reciprocal part, requiring seven separate
 197 FFTs, and therefore this has been the limiting factor in previous
 198 attempts to implement LJ-PME. A solution to this problem is to use
 199 geometrical combination rules in order to calculate an approximate
 200 interaction parameter for the reciprocal part of the potential, yielding
 201 a total interaction of
 202
 203 .. math:: \begin{aligned}
 204           V(r<r_c) & = & \underbrace{C^{\mathrm{dir}}_6 g(\beta r) r^{-6}}_{\mathrm{Direct \  space}} + \underbrace{C^\mathrm{recip}_{6,\mathrm{geom}} [1 - g(\beta r)] r^{-6}}_{\mathrm{Reciprocal \  space}} \nonumber \\
 205           &=& C^\mathrm{recip}_{6,\mathrm{geom}}r^{-6} + \left(C^{\mathrm{dir}}_6-C^\mathrm{recip}_{6,\mathrm{geom}}\right)g(\beta r)r^{-6} \\
 206           V(r>r_c) & = & \underbrace{C^\mathrm{recip}_{6,\mathrm{geom}} [1 - g(\beta r)] r^{-6}}_{\mathrm{Reciprocal \  space}}.\end{aligned}
 207           :label: eqnpmearith
 208
 209 This will preserve a well-defined Hamiltonian and significantly
 210 increase the performance of the simulations. The approximation does
 211 introduce some errors, but since the difference is located in the
 212 interactions calculated in reciprocal space, the effect will be very
 213 small compared to the total interaction energy. In a simulation of a
 214 lipid bilayer, using a cut-off of 1.0 nm, the relative error in total
 215 dispersion energy was below 0.5%. A more thorough discussion of this can
 216 be found in :ref:`109 <refWennberg13>`.
 217
 218 In |Gromacs| we now perform the proper calculation of this interaction by
 219 subtracting, from the direct-space interactions, the contribution made
 220 by the approximate potential that is used in the reciprocal part
 221
 222 .. math:: V_\mathrm{dir} = C^{\mathrm{dir}}_6 r^{-6} - C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}.
 223           :label: eqnljpmedirectspace
 224
 225 This potential will reduce to the expression in
 226 :eq:`eqn. %s <eqnljpmerealspace>` when
 227 :math:`C^{\mathrm{dir}}_6 = C^\mathrm{recip}_6`, and the total
 228 interaction is given by
 229
 230 .. math:: \begin{aligned}
 231           \nonumber V(r<r_c) &=& \underbrace{C^{\mathrm{dir}}_6 r^{-6} - C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}}_{\mathrm{Direct \  space}} + \underbrace{C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}}_{\mathrm{Reciprocal \  space}} \\
 232           &=&C^{\mathrm{dir}}_6 r^{-6}
 233           \end{aligned}
 234           :label: eqnljpmecorr2
 235
 236 .. math:: \begin{aligned}
 237           V(r>r_c) &=& C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}.\end{aligned}
 238           :label: eqnljpmecorr3
 239
 240 For the case when :math:`C^{\mathrm{dir}}_6 \neq C^\mathrm{recip}_6`
 241 this will retain an unmodified LJ force up to the cut-off, and the error
 242 is an order of magnitude smaller than in simulations where the
 243 direct-space interactions do not account for the approximation used in
 244 reciprocal space. When using a VdW interaction modifier of
 245 potential-shift, the constant
 246
 247 .. math:: \left(-C^{\mathrm{dir}}_6 + C^\mathrm{recip}_6 [1 - g(\beta r_c)]\right) r_c^{-6}
 248           :label: eqnljpmeconstant
 249
 250 is added to :eq:`eqn. %s <eqnljpmecorr2>` in order to ensure that the potential
 251 is continuous at the cutoff. Note that, in the same way as
 252 :eq:`eqn. %s <eqnljpmedirectspace>`, this degenerates into the expected
 253 :math:`-C_6g(\beta r_c)r^{-6}_c` when :math:`C^{\mathrm{dir}}_6 =
 254 C^\mathrm{recip}_6`. In addition to this, a long-range dispersion
 255 correction can be applied to correct for the approximation using a
 256 combination rule in reciprocal space. This correction assumes, as for
 257 the cut-off LJ potential, a uniform particle distribution. But since the
 258 error of the combination rule approximation is very small this
 259 long-range correction is not necessary in most cases. Also note that
 260 this homogenous correction does not correct the surface tension, which
 261 is an inhomogeneous property.
 262
 263 Using LJ-PME
 264 ^^^^^^^^^^^^
 265
 266 As an example for using Particle-mesh Ewald summation for Lennard-Jones
 267 interactions in |Gromacs|, specify the following lines in your :ref:`mdp` file:
 268
 269 ::
 270
 271     vdwtype          = PME
 272     rvdw             = 0.9
 273     vdw-modifier     = Potential-Shift
 274     rlist            = 0.9
 275     rcoulomb         = 0.9
 276     fourierspacing   = 0.12
 277     pme-order        = 4
 278     ewald-rtol-lj    = 0.001
 279     lj-pme-comb-rule = geometric
 280
 281 The same Fourier grid and interpolation order are used if both LJ-PME
 282 and electrostatic PME are active, so the settings for
 283 ``fourierspacing`` and ``pme-order`` are common
 284 to both. ``ewald-rtol-lj`` controls the splitting between
 285 direct and reciprocal space in the same way as
 286 ``ewald-rtol``. In addition to this, the combination rule to
 287 be used in reciprocal space is determined by
 288 ``lj-pme-comb-rule``. If the current force field uses
 289 Lorentz-Berthelot combination rules, it is possible to set
 290 ``lj-pme-comb-rule = geometric`` in order to gain a
 291 significant increase in performance for a small loss in accuracy. The
 292 details of this approximation can be found in the section above.
 293
 294 Note that the use of a complete long-range dispersion correction means
 295 that as with Coulomb PME, ``rvdw`` is now a free parameter
 296 in the method, rather than being necessarily restricted by the
 297 force-field parameterization scheme. Thus it is now possible to optimize
 298 the cutoff, spacing, order and tolerance terms for accuracy and best
 299 performance.
 300
 301 Naturally, the use of LJ-PME rather than LJ cut-off adds computation and
 302 communication done for the reciprocal-space part, so for best
 303 performance in balancing the load of parallel simulations using PME-only
 304 ranks, more such ranks should be used. It may be possible to improve
 305 upon the automatic load-balancing used by :ref:`mdrun <gmx mdrun>`.