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