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  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  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.