docs/reference-manual/long-range-electrostatics.rst

   1 .. _lrelstat:
   2
   3 Long Range Electrostatics
   4 -------------------------
   5
   6 Ewald summation
   7 ~~~~~~~~~~~~~~~
   8
   9 The total electrostatic energy of :math:`N` particles and their periodic
  10 images is given by
  11
  12 .. math:: V=\frac{f}{2}\sum_{n_x}\sum_{n_y}
  13           \sum_{n_{z}*} \sum_{i}^{N} \sum_{j}^{N}
  14           \frac{q_i q_j}{{\bf r}_{ij,{\bf n}}}.
  15           :label: eqntotalcoulomb
  16
  17 :math:`(n_x,n_y,n_z)={\bf n}` is the box index vector, and the star
  18 indicates that terms with :math:`i=j` should be omitted when
  19 :math:`(n_x,n_y,n_z)=(0,0,0)`. The distance :math:`{\bf r}_{ij,{\bf n}}`
  20 is the real distance between the charges and not the minimum-image. This
  21 sum is conditionally convergent, but very slow.
  22
  23 Ewald summation was first introduced as a method to calculate long-range
  24 interactions of the periodic images in crystals \ :ref:`105 <refEwald21>`. The idea
  25 is to convert the single slowly-converging sum :eq:`eqn. %s <eqntotalcoulomb>`
  26 into two quickly-converging terms and a constant term:
  27
  28 .. math::
  29
  30    \begin{aligned}
  31    V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex]
  32    V_{\mathrm{dir}} &=& \frac{f}{2} \sum_{i,j}^{N}
  33    \sum_{n_x}\sum_{n_y}
  34    \sum_{n_{z}*} q_i q_j \frac{\mbox{erfc}(\beta {r}_{ij,{\bf n}} )}{{r}_{ij,{\bf n}}} \\[0.5ex]
  35    V_{\mathrm{rec}} &=& \frac{f}{2 \pi V} \sum_{i,j}^{N} q_i q_j
  36    \sum_{m_x}\sum_{m_y}
  37    \sum_{m_{z}*} \frac{\exp{\left( -(\pi {\bf m}/\beta)^2 + 2 \pi i
  38          {\bf m} \cdot ({\bf r}_i - {\bf r}_j)\right)}}{{\bf m}^2} \\[0.5ex]
  39    V_{0} &=& -\frac{f \beta}{\sqrt{\pi}}\sum_{i}^{N} q_i^2,\end{aligned}
  40
  41 where :math:`\beta` is a parameter that determines the relative weight
  42 of the direct and reciprocal sums and :math:`{\bf m}=(m_x,m_y,m_z)`. In
  43 this way we can use a short cut-off (of the order of :math:`1`  nm) in
  44 the direct space sum and a short cut-off in the reciprocal space sum
  45 (*e.g.* 10 wave vectors in each direction). Unfortunately, the
  46 computational cost of the reciprocal part of the sum increases as
  47 :math:`N^2` (or :math:`N^{3/2}` with a slightly better algorithm) and it
  48 is therefore not realistic for use in large systems.
  49
  50 Using Ewald
  51 ^^^^^^^^^^^
  52
  53 Don’t use Ewald unless you are absolutely sure this is what you want -
  54 for almost all cases the PME method below will perform much better. If
  55 you still want to employ classical Ewald summation enter this in your
  56 :ref:`mdp` file, if the side of your box is about :math:`3`  nm:
  57
  58 ::
  59
  60     coulombtype     = Ewald
  61     rvdw            = 0.9
  62     rlist           = 0.9
  63     rcoulomb        = 0.9
  64     fourierspacing  = 0.6
  65     ewald-rtol      = 1e-5
  66
  67 The ratio of the box dimensions and the fourierspacing parameter
  68 determines the highest magnitude of wave vectors :math:`m_x,m_y,m_z` to
  69 use in each direction. With a 3-nm cubic box this example would use
  70 :math:`11` wave vectors (from :math:`-5` to :math:`5`) in each
  71 direction. The ``ewald-rtol`` parameter is the relative strength of the
  72 electrostatic interaction at the cut-off. Decreasing this gives you a
  73 more accurate direct sum, but a less accurate reciprocal sum.
  74
  75 .. _pme:
  76
  77 PME
  78 ~~~
  79
  80 Particle-mesh Ewald is a method proposed by Tom
  81 Darden \ :ref:`14 <refDarden93>` to improve the performance of the reciprocal sum.
  82 Instead of directly summing wave vectors, the charges are assigned to a
  83 grid using interpolation. The implementation in |Gromacs| uses cardinal
  84 B-spline interpolation \ :ref:`15 <refEssmann95>`, which is referred to as
  85 smooth PME (SPME). The grid is then Fourier transformed with a 3D FFT
  86 algorithm and the reciprocal energy term obtained by a single sum over
  87 the grid in k-space.
  88
  89 The potential at the grid points is calculated by inverse
  90 transformation, and by using the interpolation factors we get the forces
  91 on each atom.
  92
  93 The PME algorithm scales as :math:`N \log(N)`, and is substantially
  94 faster than ordinary Ewald summation on medium to large systems. On very
  95 small systems it might still be better to use Ewald to avoid the
  96 overhead in setting up grids and transforms. For the parallelization of
  97 PME see the section on MPMD PME (:ref:`mpmdpme`).
  98
  99 With the Verlet cut-off scheme, the PME direct space potential is
 100 shifted by a constant such that the potential is zero at the cut-off.
 101 This shift is small and since the net system charge is close to zero,
 102 the total shift is very small, unlike in the case of the Lennard-Jones
 103 potential where all shifts add up. We apply the shift anyhow, such that
 104 the potential is the exact integral of the force.
 105
 106 Using PME
 107 ^^^^^^^^^
 108
 109 As an example for using Particle-mesh Ewald summation in |Gromacs|,
 110 specify the following lines in your :ref:`mdp` file:
 111
 112 ::
 113
 114     coulombtype     = PME
 115     rvdw            = 0.9
 116     rlist           = 0.9
 117     rcoulomb        = 0.9
 118     fourierspacing  = 0.12
 119     pme-order       = 4
 120     ewald-rtol      = 1e-5
 121
 122 In this case the ``fourierspacing`` parameter determines the
 123 maximum spacing for the FFT grid (i.e. minimum number of grid points),
 124 and ``pme-order`` controls the interpolation order. Using
 125 fourth-order (cubic) interpolation and this spacing should give
 126 electrostatic energies accurate to about :math:`5\cdot10^{-3}`. Since
 127 the Lennard-Jones energies are not this accurate it might even be
 128 possible to increase this spacing slightly.
 129
 130 Pressure scaling works with PME, but be aware of the fact that
 131 anisotropic scaling can introduce artificial ordering in some systems.
 132
 133 P3M-AD
 134 ~~~~~~
 135
 136 The Particle-Particle
 137 Particle-Mesh
 138 methods of Hockney & Eastwood can also be applied in |Gromacs| for the
 139 treatment of long range electrostatic
 140 interactions \ :ref:`106 <refHockney81>`. Although the P3M method was the first efficient long-range
 141 electrostatics method for molecular simulation, the smooth PME (SPME)
 142 method has largely replaced P3M as the method of choice in atomistic
 143 simulations. One performance disadvantage of the original P3M method was
 144 that it required 3 3D-FFT back transforms to obtain the forces on the
 145 particles. But this is not required for P3M and the forces can be
 146 derived through analytical differentiation of the potential, as done in
 147 PME. The resulting method is termed P3M-AD. The only remaining
 148 difference between P3M-AD and PME is the optimization of the lattice
 149 Green influence function for error minimization that P3M uses. However,
 150 in 2012 it has been shown that the SPME influence function can be
 151 modified to obtain P3M \ :ref:`107 <refBallenegger2012>`. This means
 152 that the advantage of error minimization in P3M-AD can be used at the
 153 same computational cost and with the same code as PME, just by adding a
 154 few lines to modify the influence function. However, at optimal
 155 parameter setting the effect of error minimization in P3M-AD is less
 156 than 10%. P3M-AD does show large accuracy gains with interlaced (also
 157 known as staggered) grids, but that is not supported in |Gromacs| (yet).
 158
 159 P3M is used in |Gromacs| with exactly the same options as used with PME by
 160 selecting the electrostatics type:
 161
 162 ::
 163
 164     coulombtype     = P3M-AD
 165
 166 Optimizing Fourier transforms and PME calculations
 167 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 168
 169 It is recommended to optimize the parameters for calculation of
 170 electrostatic interaction such as PME grid dimensions and cut-off radii.
 171 This is particularly relevant to do before launching long production
 172 runs.
 173
 174 :ref:`gmx mdrun` will automatically do a lot of PME
 175 optimization, and |Gromacs| also includes a special tool,
 176 :ref:`gmx tune_pme`, which automates the process of selecting
 177 the optimal number of PME-only ranks.