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