Non-bonded interactions ----------------------- Non-bonded interactions in |Gromacs| are pair-additive: .. math:: V(\mathbf{r}_1,\ldots \mathbf{r}_N) = \sum_{i The parameters :math:C^{(12)}_{ij} and :math:C^{(6)}_{ij} depend on pairs of *atom types*; consequently they are taken from a matrix of LJ-parameters. In the Verlet cut-off scheme, the potential is shifted by a constant such that it is zero at the cut-off distance. .. _fig-lj: .. figure:: plots/f-lj.* :width: 8.00000cm The Lennard-Jones interaction. The force derived from this potential is: .. math:: \mathbf{F}_i(\mathbf{r}_ij) = \left( 12~\frac{C_{ij}^{(12)}}{{r_{ij}}^{13}} - 6~\frac{C_{ij}^{(6)}}{{r_{ij}}^7} \right) {\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}} :label: eqnljforce The LJ potential may also be written in the following form: .. math:: V_{LJ}(\mathbf{r}_ij) = 4\epsilon_{ij}\left(\left(\frac{\sigma_{ij}} {{r_{ij}}}\right)^{12} - \left(\frac{\sigma_{ij}}{{r_{ij}}}\right)^{6} \right) :label: eqnsigeps In constructing the parameter matrix for the non-bonded LJ-parameters, two types of combination rules can be used within |Gromacs|, only geometric averages (type 1 in the input section of the force-field file): .. math:: \begin{array}{rcl} C_{ij}^{(6)} &=& \left({C_{ii}^{(6)} \, C_{jj}^{(6)}}\right)^{1/2} \\ C_{ij}^{(12)} &=& \left({C_{ii}^{(12)} \, C_{jj}^{(12)}}\right)^{1/2} \end{array} :label: eqncomb or, alternatively the Lorentz-Berthelot rules can be used. An arithmetic average is used to calculate :math:\sigma_{ij}, while a geometric average is used to calculate :math:\epsilon_{ij} (type 2): .. math:: \begin{array}{rcl} \sigma_{ij} &=& \frac{1}{ 2}(\sigma_{ii} + \sigma_{jj}) \\ \epsilon_{ij} &=& \left({\epsilon_{ii} \, \epsilon_{jj}}\right)^{1/2} \end{array} :label: eqnlorentzberthelot finally an geometric average for both parameters can be used (type 3): .. math:: \begin{array}{rcl} \sigma_{ij} &=& \left({\sigma_{ii} \, \sigma_{jj}}\right)^{1/2} \\ \epsilon_{ij} &=& \left({\epsilon_{ii} \, \epsilon_{jj}}\right)^{1/2} \end{array} :label: eqnnbgeometricaverage This last rule is used by the OPLS force field. Buckingham potential ~~~~~~~~~~~~~~~~~~~~ The Buckingham potential has a more flexible and realistic repulsion term than the Lennard-Jones interaction, but is also more expensive to compute. The potential form is: .. math:: V_{bh}({r_{ij}}) = A_{ij} \exp(-B_{ij} {r_{ij}}) - \frac{C_{ij}}{{r_{ij}}^6} :label: eqnnbbuckingham .. _fig-bham: .. figure:: plots/f-bham.* :width: 8.00000cm The Buckingham interaction. See also :numref:Fig. %s . The force derived from this is: .. math:: \mathbf{F}_i({r_{ij}}) = \left[ A_{ij}B_{ij}\exp(-B_{ij} {r_{ij}}) - 6\frac{C_{ij}}{{r_{ij}}^7} \right] {\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}} :label: eqnnbbuckinghamforce .. _coul: Coulomb interaction ~~~~~~~~~~~~~~~~~~~ The Coulomb interaction between two charge particles is given by: .. math:: V_c({r_{ij}}) = f \frac{q_i q_j}{{\varepsilon_r}{r_{ij}}} :label: eqnvcoul See also :numref:Fig. %s , where :math:f = \frac{1}{4\pi \varepsilon_0} = {138.935\,458} (see chapter :ref:defunits) .. _fig-coul: .. figure:: plots/vcrf.* :width: 8.00000cm The Coulomb interaction (for particles with equal signed charge) with and without reaction field. In the latter case :math:{\varepsilon_r} was 1, :math:{\varepsilon_{rf}} was 78, and :math:r_c was 0.9 nm. The dot-dashed line is the same as the dashed line, except for a constant. The force derived from this potential is: .. math:: \mathbf{F}_i(\mathbf{r}_ij) = f \frac{q_i q_j}{{\varepsilon_r}{r_{ij}}^2}{\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}} :label: eqnfcoul A plain Coulomb interaction should only be used without cut-off or when all pairs fall within the cut-off, since there is an abrupt, large change in the force at the cut-off. In case you do want to use a cut-off, the potential can be shifted by a constant to make the potential the integral of the force. With the group cut-off scheme, this shift is only applied to non-excluded pairs. With the Verlet cut-off scheme, the shift is also applied to excluded pairs and self interactions, which makes the potential equivalent to a reaction field with :math:{\varepsilon_{rf}}=1 (see below). In |Gromacs| the relative dielectric constant :math:{\varepsilon_r} may be set in the in the input for :ref:grompp . .. _coulrf: Coulomb interaction with reaction field ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The Coulomb interaction can be modified for homogeneous systems by assuming a constant dielectric environment beyond the cut-off :math:r_c with a dielectric constant of :math:{\varepsilon_{rf}}. The interaction then reads: .. math:: V_{crf} ~=~ f \frac{q_i q_j}{{\varepsilon_r}{r_{ij}}}\left[1+\frac{{\varepsilon_{rf}}-{\varepsilon_r}}{2{\varepsilon_{rf}}+{\varepsilon_r}} \,\frac{{r_{ij}}^3}{r_c^3}\right] - f\frac{q_i q_j}{{\varepsilon_r}r_c}\,\frac{3{\varepsilon_{rf}}}{2{\varepsilon_{rf}}+{\varepsilon_r}} :label: eqnvcrf in which the constant expression on the right makes the potential zero at the cut-off :math:r_c. For charged cut-off spheres this corresponds to neutralization with a homogeneous background charge. We can rewrite :eq:eqn. %s  for simplicity as .. math:: V_{crf} ~=~ f \frac{q_i q_j}{{\varepsilon_r}}\left[\frac{1}{{r_{ij}}} + k_{rf}~ {r_{ij}}^2 -c_{rf}\right] :label: eqnvcrfrewrite with .. math:: \begin{aligned} k_{rf} &=& \frac{1}{r_c^3}\,\frac{{\varepsilon_{rf}}-{\varepsilon_r}}{(2{\varepsilon_{rf}}+{\varepsilon_r})} \end{aligned} :label: eqnkrf .. math:: \begin{aligned} c_{rf} &=& \frac{1}{r_c}+k_{rf}\,r_c^2 ~=~ \frac{1}{r_c}\,\frac{3{\varepsilon_{rf}}}{(2{\varepsilon_{rf}}+{\varepsilon_r})} \end{aligned} :label: eqncrf For large :math:{\varepsilon_{rf}} the :math:k_{rf} goes to :math:r_c^{-3}/2, while for :math:{\varepsilon_{rf}} = :math:{\varepsilon_r} the correction vanishes. In :numref:Fig. %s  the modified interaction is plotted, and it is clear that the derivative with respect to :math:{r_{ij}} (= -force) goes to zero at the cut-off distance. The force derived from this potential reads: .. math:: \mathbf{F}_i(\mathbf{r}_ij) = f \frac{q_i q_j}{{\varepsilon_r}}\left[\frac{1}{{r_{ij}}^2} - 2 k_{rf}{r_{ij}}\right]{\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}} :label: eqnfcrf The reaction-field correction should also be applied to all excluded atoms pairs, including self pairs, in which case the normal Coulomb term in :eq:eqns. %s  and :eq:%s  is absent. Tironi *et al.* have introduced a generalized reaction field in which the dielectric continuum beyond the cut-off :math:r_c also has an ionic strength :math:I :ref:71 . In this case we can rewrite the constants :math:k_{rf} and :math:c_{rf} using the inverse Debye screening length :math:\kappa: .. math:: \begin{aligned} \kappa^2 &=& \frac{2 I \,F^2}{\varepsilon_0 {\varepsilon_{rf}}RT} = \frac{F^2}{\varepsilon_0 {\varepsilon_{rf}}RT}\sum_{i=1}^{K} c_i z_i^2 \\ k_{rf} &=& \frac{1}{r_c^3}\, \frac{({\varepsilon_{rf}}-{\varepsilon_r})(1 + \kappa r_c) + {\frac{1}{2}}{\varepsilon_{rf}}(\kappa r_c)^2} {(2{\varepsilon_{rf}}+ {\varepsilon_r})(1 + \kappa r_c) + {\varepsilon_{rf}}(\kappa r_c)^2} \end{aligned} :label: eqnkgrf .. math:: \begin{aligned} c_{rf} &=& \frac{1}{r_c}\, \frac{3{\varepsilon_{rf}}(1 + \kappa r_c + {\frac{1}{2}}(\kappa r_c)^2)} {(2{\varepsilon_{rf}}+{\varepsilon_r})(1 + \kappa r_c) + {\varepsilon_{rf}}(\kappa r_c)^2} \end{aligned} :label: eqncgrf where :math:F is Faraday’s constant, :math:R is the ideal gas constant, :math:T the absolute temperature, :math:c_i the molar concentration for species :math:i and :math:z_i the charge number of species :math:i where we have :math:K different species. In the limit of zero ionic strength (:math:\kappa=0) :eq:eqns. %s  and :eq:%s  reduce to the simple forms of :eq:eqns. %s  and :eq:%s  respectively. .. _modnbint: Modified non-bonded interactions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ In |Gromacs|, the non-bonded potentials can be modified by a shift function, also called a force-switch function, since it switches the force to zero at the cut-off. The purpose of this is to replace the truncated forces by forces that are continuous and have continuous derivatives at the cut-off radius. With such forces the time integration produces smaller errors. But note that for Lennard-Jones interactions these errors are usually smaller than other errors, such as integration errors at the repulsive part of the potential. For Coulomb interactions we advise against using a shifted potential and for use of a reaction field or a proper long-range method such as PME. There is *no* fundamental difference between a switch function (which multiplies the potential with a function) and a shift function (which adds a function to the force or potential) \ :ref:72 . The switch function is a special case of the shift function, which we apply to the *force function* :math:F(r), related to the electrostatic or van der Waals force acting on particle :math:i by particle :math:j as: .. math:: \mathbf{F}_i = c \, F(r_{ij}) \frac{\mathbf{r}_ij}{r_{ij}} :label: eqnswitch For pure Coulomb or Lennard-Jones interactions :math:F(r) = F_\alpha(r) = \alpha \, r^{-(\alpha+1)}. The switched force :math:F_s(r) can generally be written as: .. math:: \begin{array}{rcl} F_s(r)~=&~F_\alpha(r) & r < r_1 \\ F_s(r)~=&~F_\alpha(r)+S(r) & r_1 \le r < r_c \\ F_s(r)~=&~0 & r_c \le r \end{array} :label: eqnswitchforce When :math:r_1=0 this is a traditional shift function, otherwise it acts as a switch function. The corresponding shifted potential function then reads: .. math:: V_s(r) = \int^{\infty}_r~F_s(x)\, dx :label: eqnswitchpotential The |Gromacs| **force switch** function :math:S_F(r) should be smooth at the boundaries, therefore the following boundary conditions are imposed on the switch function: .. math:: \begin{array}{rcl} S_F(r_1) &=&0 \\ S_F'(r_1) &=&0 \\ S_F(r_c) &=&-F_\alpha(r_c) \\ S_F'(r_c) &=&-F_\alpha'(r_c) \end{array} :label: eqnswitchforcefunction A 3\ :math:^{rd} degree polynomial of the form .. math:: S_F(r) = A(r-r_1)^2 + B(r-r_1)^3 :label: eqnswitchforcepoly fulfills these requirements. The constants A and B are given by the boundary condition at :math:r_c: .. math:: \begin{array}{rcl} A &~=~& -\alpha \, \displaystyle \frac{(\alpha+4)r_c~-~(\alpha+1)r_1} {r_c^{\alpha+2}~(r_c-r_1)^2} \\ B &~=~& \alpha \, \displaystyle \frac{(\alpha+3)r_c~-~(\alpha+1)r_1}{r_c^{\alpha+2}~(r_c-r_1)^3} \end{array} :label: eqnforceswitchboundary Thus the total force function is: .. math:: F_s(r) = \frac{\alpha}{r^{\alpha+1}} + A(r-r_1)^2 + B(r-r_1)^3 :label: eqnswitchfinalforce and the potential function reads: .. math:: V_s(r) = \frac{1}{r^\alpha} - \frac{A}{3} (r-r_1)^3 - \frac{B}{4} (r-r_1)^4 - C :label: eqnswitchfinalpotential where .. math:: C = \frac{1}{r_c^\alpha} - \frac{A}{3} (r_c-r_1)^3 - \frac{B}{4} (r_c-r_1)^4 :label: eqnswitchpotentialexp The |Gromacs| **potential-switch** function :math:S_V(r) scales the potential between :math:r_1 and :math:r_c, and has similar boundary conditions, intended to produce smoothly-varying potential and forces: .. math:: \begin{array}{rcl} S_V(r_1) &=&1 \\ S_V'(r_1) &=&0 \\ S_V''(r_1) &=&0 \\ S_V(r_c) &=&0 \\ S_V'(r_c) &=&0 \\ S_V''(r_c) &=&0 \end{array} :label: eqnpotentialswitch The fifth-degree polynomial that has these properties is .. math:: S_V(r; r_1, r_c) = \frac{1 - 10(r-r_1)^3(r_c-r_1)^2 + 15(r-r_1)^4(r_c-r_1) - 6(r-r_1)}{(r_c-r_1)^5} :label: eqn5polynomal This implementation is found in several other simulation packages,\ :ref:73 \ :ref:75  but differs from that in CHARMM.\ :ref:76  Switching the potential leads to artificially large forces in the switching region, therefore it is not recommended to switch Coulomb interactions using this function,\ :ref:72  but switching Lennard-Jones interactions using this function produces acceptable results. Modified short-range interactions with Ewald summation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When Ewald summation or particle-mesh Ewald is used to calculate the long-range interactions, the short-range Coulomb potential must also be modified. Here the potential is switched to (nearly) zero at the cut-off, instead of the force. In this case the short range potential is given by: .. math:: V(r) = f \frac{\mbox{erfc}(\beta r_{ij})}{r_{ij}} q_i q_j, :label: eqnewaldsrmod where :math:\beta is a parameter that determines the relative weight between the direct space sum and the reciprocal space sum and erfc\ :math:(x) is the complementary error function. For further details on long-range electrostatics, see sec. :ref:lrelstat`.