docs/reference-manual/tabulated-interaction-functions.rst

   1 Tabulated interaction functions
   2 -------------------------------
   3
   4 .. _cubicspline:
   5
   6 Cubic splines for potentials
   7 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   8
   9 In some of the inner loops of |Gromacs|, look-up tables are used for
  10 computation of potential and forces. The tables are interpolated using a
  11 cubic spline algorithm. There are separate tables for electrostatic,
  12 dispersion, and repulsion interactions, but for the sake of caching
  13 performance these have been combined into a single array. The cubic
  14 spline interpolation for :math:`x_i \leq x < x_{i+1}` looks like this:
  15
  16 .. math::  V_s(x) = A_0 + A_1 \,\epsilon + A_2 \,\epsilon^2 + A_3 \,\epsilon^3
  17            :label: eqnspline
  18
  19 where the table spacing :math:`h` and fraction :math:`\epsilon` are
  20 given by:
  21
  22 .. math::
  23
  24    \begin{aligned}
  25    h    &=&     x_{i+1} - x_i   \\
  26    \epsilon&=&  (x - x_i)/h\end{aligned}
  27
  28 so that :math:`0 \le \epsilon < 1`. From this, we can calculate the
  29 derivative in order to determine the forces:
  30
  31 .. math::
  32
  33    -V_s'(x) ~=~
  34    -\frac{{\rm d}V_s(x)}{{\rm d}\epsilon}\frac{{\rm d}\epsilon}{{\rm d}x} ~=~
  35    -(A_1 + 2 A_2 \,\epsilon + 3 A_3 \,\epsilon^2)/h
  36
  37 The four coefficients are determined from the four conditions that
  38 :math:`V_s` and :math:`-V_s'` at both ends of each interval should match
  39 the exact potential :math:`V` and force :math:`-V'`. This results in the
  40 following errors for each interval:
  41
  42 .. math::
  43
  44    \begin{aligned}
  45    | V_s  - V  | _{max} &=& V'''' \frac{h^4}{384} + O(h^5) \\
  46    | V_s' - V' | _{max} &=& V'''' \frac{h^3}{72\sqrt{3}} + O(h^4) \\
  47    | V_s''- V''| _{max} &=& V'''' \frac{h^2}{12}  + O(h^3)\end{aligned}
  48
  49 V and V’ are continuous, while V” is the first discontinuous
  50 derivative. The number of points per nanometer is 500 and 2000 for
  51 mixed- and double-precision versions of |Gromacs|, respectively. This
  52 means that the errors in the potential and force will usually be smaller
  53 than the mixed precision accuracy.
  54
  55 |Gromacs| stores :math:`A_0`, :math:`A_1`, :math:`A_2` and :math:`A_3`.
  56 The force routines get a table with these four parameters and a scaling
  57 factor :math:`s` that is equal to the number of points per nm. (**Note**
  58 that :math:`h` is :math:`s^{-1}`). The algorithm goes a little something
  59 like this:
  60
  61 #. Calculate distance vector
  62    (:math:`\mathbf{r}_{ij}`) and distance
  63    r\ :math:`_{ij}`
  64
  65 #. Multiply r\ :math:`_{ij}` by :math:`s` and truncate to an integer
  66    value :math:`n_0` to get a table index
  67
  68 #. Calculate fractional component (:math:`\epsilon` =
  69    :math:`s`\ r\ :math:`_{ij} - n_0`) and :math:`\epsilon^2`
  70
  71 #. Do the interpolation to calculate the potential :math:`V` and the
  72    scalar force :math:`f`
  73
  74 #. Calculate the vector force :math:`\mathbf{F}` by
  75    multiplying :math:`f` with
  76    :math:`\mathbf{r}_{ij}`
  77
  78 **Note** that table look-up is significantly *slower* than computation
  79 of the most simple Lennard-Jones and Coulomb interaction. However, it is
  80 much faster than the shifted Coulomb function used in conjunction with
  81 the PPPM method. Finally, it is much easier to modify a table for the
  82 potential (and get a graphical representation of it) than to modify the
  83 inner loops of the MD program.
  84
  85 User-specified potential functions
  86 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  87
  88 You can also use your own potential functions without editing the
  89 |Gromacs| code. The potential function should be according to the
  90 following equation
  91
  92 .. math:: V(r_{ij}) ~=~ \frac{q_i q_j}{4 \pi\epsilon_0} f(r_{ij}) + C_6 \,g(r_{ij}) + C_{12} \,h(r_{ij})
  93
  94 where :math:`f`, :math:`g`, and :math:`h` are user defined functions.
  95 **Note** that if :math:`g(r)` represents a normal dispersion
  96 interaction, :math:`g(r)` should be :math:`<` 0. C\ :math:`_6`,
  97 C\ :math:`_{12}` and the charges are read from the topology. Also note
  98 that combination rules are only supported for Lennard-Jones and
  99 Buckingham, and that your tables should match the parameters in the
 100 binary topology.
 101
 102 When you add the following lines in your :ref:`mdp` file:
 103
 104 ::
 105
 106     rlist           = 1.0
 107     coulombtype     = User
 108     rcoulomb        = 1.0
 109     vdwtype         = User
 110     rvdw            = 1.0
 111
 112 :ref:`mdrun <gmx mdrun>` will read a single non-bonded table file, or
 113 multiple when ``energygrp-table`` is set (see below). The
 114 name of the file(s) can be set with the :ref:`mdrun <gmx mdrun>` option
 115 ``-table``. The table file should contain seven columns of
 116 table look-up data in the order: :math:`x`, :math:`f(x)`,
 117 :math:`-f'(x)`, :math:`g(x)`, :math:`-g'(x)`, :math:`h(x)`,
 118 :math:`-h'(x)`. The :math:`x` should run from 0 to :math:`r_c+1` (the
 119 value of ``table_extension`` can be changed in the :ref:`mdp` file). You can
 120 choose the spacing you like; for the standard tables |Gromacs| uses a
 121 spacing of 0.002 and 0.0005 nm when you run in mixed and double
 122 precision, respectively. In this context, :math:`r_c` denotes the
 123 maximum of the two cut-offs ``rvdw`` and ``rcoulomb`` (see above). These
 124 variables need not be the same (and need not be 1.0 either). Some
 125 functions used for potentials contain a singularity at :math:`x = 0`,
 126 but since atoms are normally not closer to each other than 0.1 nm, the
 127 function value at :math:`x = 0` is not important. Finally, it is also
 128 possible to combine a standard Coulomb with a modified LJ potential (or
 129 vice versa). One then specifies *e.g.* ``coulombtype = Cut-off`` or
 130 ``coulombtype = PME``, combined with ``vdwtype = User``. The table file must
 131 always contain the 7 columns however, and meaningful data (i.e. not
 132 zeroes) must be entered in all columns. A number of pre-built table
 133 files can be found in the ``GMXLIB`` directory for 6-8, 6-9, 6-10, 6-11, and
 134 6-12 Lennard-Jones potentials combined with a normal Coulomb.
 135
 136 If you want to have different functional forms between different groups
 137 of atoms, this can be set through energy groups. Different tables can be
 138 used for non-bonded interactions between different energy groups pairs
 139 through the :ref:`mdp` option ``energygrp-table`` (see details in the User Guide).
 140 Atoms that should interact with a different potential should be put into
 141 different energy groups. Between group pairs which are not listed in
 142 ``energygrp-table``, the normal user tables will be used. This makes it easy
 143 to use a different functional form between a few types of atoms.