1 Non-bonded interactions
2 -----------------------
4 Non-bonded interactions in |Gromacs| are pair-additive:
6 .. math:: V(\mathbf{r}_1,\ldots \mathbf{r}_N) = \sum_{i<j}V_{ij}(\mathbf{r}_ij);
8 .. math:: \mathbf{F}_i = -\sum_j \frac{dV_{ij}(r_{ij})}{dr_{ij}} \frac{\mathbf{r}_ij}{r_{ij}}
10 Since the potential only depends on the scalar distance, interactions
11 will be centro-symmetric, i.e. the vectorial partial force on particle
12 :math:`i` from the pairwise interaction :math:`V_{ij}(r_{ij})` has the
13 opposite direction of the partial force on particle :math:`j`. For
14 efficiency reasons, interactions are calculated by loops over
15 interactions and updating both partial forces rather than summing one
16 complete nonbonded force at a time. The non-bonded interactions contain
17 a repulsion term, a dispersion term, and a Coulomb term. The repulsion
18 and dispersion term are combined in either the Lennard-Jones (or 6-12
19 interaction), or the Buckingham (or exp-6 potential). In addition,
20 (partially) charged atoms act through the Coulomb term.
24 The Lennard-Jones interaction
25 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
27 The Lennard-Jones potential :math:`V_{LJ}` between two atoms equals:
31 V_{LJ}({r_{ij}}) = \frac{C_{ij}^{(12)}}{{r_{ij}}^{12}} -
32 \frac{C_{ij}^{(6)}}{{r_{ij}}^6}
34 See also :numref:`Fig. %s <fig-lj>` The parameters :math:`C^{(12)}_{ij}` and
35 :math:`C^{(6)}_{ij}` depend on pairs of *atom types*; consequently they
36 are taken from a matrix of LJ-parameters. In the Verlet cut-off scheme,
37 the potential is shifted by a constant such that it is zero at the
42 .. figure:: plots/f-lj.*
45 The Lennard-Jones interaction.
47 The force derived from this potential is:
49 .. math:: \mathbf{F}_i(\mathbf{r}_ij) = \left( 12~\frac{C_{ij}^{(12)}}{{r_{ij}}^{13}} -
50 6~\frac{C_{ij}^{(6)}}{{r_{ij}}^7} \right) {\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}}
52 The LJ potential may also be written in the following form:
54 .. math:: V_{LJ}(\mathbf{r}_ij) = 4\epsilon_{ij}\left(\left(\frac{\sigma_{ij}} {{r_{ij}}}\right)^{12}
55 - \left(\frac{\sigma_{ij}}{{r_{ij}}}\right)^{6} \right)
58 In constructing the parameter matrix for the non-bonded LJ-parameters,
59 two types of combination rules can be used within |Gromacs|, only
60 geometric averages (type 1 in the input section of the force-field
63 .. math:: \begin{array}{rcl}
64 C_{ij}^{(6)} &=& \left({C_{ii}^{(6)} \, C_{jj}^{(6)}}\right)^{1/2} \\
65 C_{ij}^{(12)} &=& \left({C_{ii}^{(12)} \, C_{jj}^{(12)}}\right)^{1/2}
69 or, alternatively the Lorentz-Berthelot rules can be used. An
70 arithmetic average is used to calculate :math:`\sigma_{ij}`, while a
71 geometric average is used to calculate :math:`\epsilon_{ij}` (type 2):
73 .. math:: \begin{array}{rcl}
74 \sigma_{ij} &=& \frac{1}{ 2}(\sigma_{ii} + \sigma_{jj}) \\
75 \epsilon_{ij} &=& \left({\epsilon_{ii} \, \epsilon_{jj}}\right)^{1/2}
77 :label: eqnlorentzberthelot
79 finally an geometric average for both parameters can be used (type 3):
81 .. math:: \begin{array}{rcl}
82 \sigma_{ij} &=& \left({\sigma_{ii} \, \sigma_{jj}}\right)^{1/2} \\
83 \epsilon_{ij} &=& \left({\epsilon_{ii} \, \epsilon_{jj}}\right)^{1/2}
86 This last rule is used by the OPLS force field.
91 The Buckingham potential has a more flexible and realistic repulsion
92 term than the Lennard-Jones interaction, but is also more expensive to
93 compute. The potential form is:
97 V_{bh}({r_{ij}}) = A_{ij} \exp(-B_{ij} {r_{ij}}) -
98 \frac{C_{ij}}{{r_{ij}}^6}
102 .. figure:: plots/f-bham.*
105 The Buckingham interaction.
107 See also :numref:`Fig. %s <fig-bham>`. The force derived from this is:
111 \mathbf{F}_i({r_{ij}}) = \left[ A_{ij}B_{ij}\exp(-B_{ij} {r_{ij}}) -
112 6\frac{C_{ij}}{{r_{ij}}^7} \right] {\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}}
119 The Coulomb interaction between two charge particles is given by:
121 .. math:: V_c({r_{ij}}) = f \frac{q_i q_j}{{\varepsilon_r}{r_{ij}}}
124 See also :numref:`Fig. %s <fig-coul>`, where
125 :math:`f = \frac{1}{4\pi \varepsilon_0} = {138.935\,458}` (see chapter :ref:`defunits`)
129 .. figure:: plots/vcrf.*
132 The Coulomb interaction (for particles with equal signed charge) with
133 and without reaction field. In the latter case
134 :math:`{\varepsilon_r}` was 1, :math:`{\varepsilon_{rf}}` was 78, and
135 :math:`r_c` was 0.9 nm. The dot-dashed line is the same as the dashed
136 line, except for a constant.
138 The force derived from this potential is:
140 .. math:: \mathbf{F}_i(\mathbf{r}_ij) = f \frac{q_i q_j}{{\varepsilon_r}{r_{ij}}^2}{\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}}
142 A plain Coulomb interaction should only be used without cut-off or when
143 all pairs fall within the cut-off, since there is an abrupt, large
144 change in the force at the cut-off. In case you do want to use a
145 cut-off, the potential can be shifted by a constant to make the
146 potential the integral of the force. With the group cut-off scheme, this
147 shift is only applied to non-excluded pairs. With the Verlet cut-off
148 scheme, the shift is also applied to excluded pairs and self
149 interactions, which makes the potential equivalent to a reaction field
150 with :math:`{\varepsilon_{rf}}=1` (see below).
152 In |Gromacs| the relative dielectric constant :math:`{\varepsilon_r}` may
153 be set in the in the input for :ref:`grompp <gmx grompp>`.
157 Coulomb interaction with reaction field
158 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
160 The Coulomb interaction can be modified for homogeneous systems by
161 assuming a constant dielectric environment beyond the cut-off
162 :math:`r_c` with a dielectric constant of :math:`{\varepsilon_{rf}}`.
163 The interaction then reads:
165 .. math:: V_{crf} ~=~
166 f \frac{q_i q_j}{{\varepsilon_r}{r_{ij}}}\left[1+\frac{{\varepsilon_{rf}}-{\varepsilon_r}}{2{\varepsilon_{rf}}+{\varepsilon_r}}
167 \,\frac{{r_{ij}}^3}{r_c^3}\right]
168 - f\frac{q_i q_j}{{\varepsilon_r}r_c}\,\frac{3{\varepsilon_{rf}}}{2{\varepsilon_{rf}}+{\varepsilon_r}}
171 in which the constant expression on the right makes the potential zero
172 at the cut-off :math:`r_c`. For charged cut-off spheres this corresponds
173 to neutralization with a homogeneous background charge. We can rewrite
174 :eq:`eqn. %s <eqnvcrf>` for simplicity as
176 .. math:: V_{crf} ~=~ f \frac{q_i q_j}{{\varepsilon_r}}\left[\frac{1}{{r_{ij}}} + k_{rf}~ {r_{ij}}^2 -c_{rf}\right]
180 .. math:: \begin{aligned}
181 k_{rf} &=& \frac{1}{r_c^3}\,\frac{{\varepsilon_{rf}}-{\varepsilon_r}}{(2{\varepsilon_{rf}}+{\varepsilon_r})}
185 .. math:: \begin{aligned}
186 c_{rf} &=& \frac{1}{r_c}+k_{rf}\,r_c^2 ~=~ \frac{1}{r_c}\,\frac{3{\varepsilon_{rf}}}{(2{\varepsilon_{rf}}+{\varepsilon_r})}
190 For large :math:`{\varepsilon_{rf}}` the :math:`k_{rf}` goes to
191 :math:`r_c^{-3}/2`, while for :math:`{\varepsilon_{rf}}` =
192 :math:`{\varepsilon_r}` the correction vanishes. In :numref:`Fig. %s <fig-coul>` the
193 modified interaction is plotted, and it is clear that the derivative
194 with respect to :math:`{r_{ij}}` (= -force) goes to zero at the cut-off
195 distance. The force derived from this potential reads:
197 .. 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}}}}
200 The reaction-field correction should also be applied to all excluded
201 atoms pairs, including self pairs, in which case the normal Coulomb term
202 in :eq:`eqns. %s <eqnvcrf>` and :eq:`%s <eqnfcrf>` is absent.
204 Tironi *et al.* have introduced a generalized reaction field in which
205 the dielectric continuum beyond the cut-off :math:`r_c` also has an
206 ionic strength :math:`I` :ref:`71 <refTironi95>`. In this case we can
207 rewrite the constants :math:`k_{rf}` and :math:`c_{rf}` using the
208 inverse Debye screening length :math:`\kappa`:
210 .. math:: \begin{aligned}
212 \frac{2 I \,F^2}{\varepsilon_0 {\varepsilon_{rf}}RT}
213 = \frac{F^2}{\varepsilon_0 {\varepsilon_{rf}}RT}\sum_{i=1}^{K} c_i z_i^2 \\
214 k_{rf} &=& \frac{1}{r_c^3}\,
215 \frac{({\varepsilon_{rf}}-{\varepsilon_r})(1 + \kappa r_c) + {\frac{1}{2}}{\varepsilon_{rf}}(\kappa r_c)^2}
216 {(2{\varepsilon_{rf}}+ {\varepsilon_r})(1 + \kappa r_c) + {\varepsilon_{rf}}(\kappa r_c)^2}
220 .. math:: \begin{aligned}
221 c_{rf} &=& \frac{1}{r_c}\,
222 \frac{3{\varepsilon_{rf}}(1 + \kappa r_c + {\frac{1}{2}}(\kappa r_c)^2)}
223 {(2{\varepsilon_{rf}}+{\varepsilon_r})(1 + \kappa r_c) + {\varepsilon_{rf}}(\kappa r_c)^2}
227 where :math:`F` is Faraday’s constant, :math:`R` is the ideal gas
228 constant, :math:`T` the absolute temperature, :math:`c_i` the molar
229 concentration for species :math:`i` and :math:`z_i` the charge number of
230 species :math:`i` where we have :math:`K` different species. In the
231 limit of zero ionic strength (:math:`\kappa=0`) :eq:`eqns. %s <eqnkgrf>` and
232 :eq:`%s <eqncgrf>` reduce to the simple forms of :eq:`eqns. %s <eqnkrf>` and :eq:`%s <eqncrf>`
237 Modified non-bonded interactions
238 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
240 In |Gromacs|, the non-bonded potentials can be modified by a shift
241 function, also called a force-switch function, since it switches the
242 force to zero at the cut-off. The purpose of this is to replace the
243 truncated forces by forces that are continuous and have continuous
244 derivatives at the cut-off radius. With such forces the time integration
245 produces smaller errors. But note that for Lennard-Jones interactions
246 these errors are usually smaller than other errors, such as integration
247 errors at the repulsive part of the potential. For Coulomb interactions
248 we advise against using a shifted potential and for use of a reaction
249 field or a proper long-range method such as PME.
251 There is *no* fundamental difference between a switch function (which
252 multiplies the potential with a function) and a shift function (which
253 adds a function to the force or potential) \ :ref:`72 <refSpoel2006a>`. The
254 switch function is a special case of the shift function, which we apply
255 to the *force function* :math:`F(r)`, related to the electrostatic or
256 van der Waals force acting on particle :math:`i` by particle :math:`j`
259 .. math:: \mathbf{F}_i = c \, F(r_{ij}) \frac{\mathbf{r}_ij}{r_{ij}}
261 For pure Coulomb or Lennard-Jones interactions
262 :math:`F(r) = F_\alpha(r) = \alpha \, r^{-(\alpha+1)}`. The switched
263 force :math:`F_s(r)` can generally be written as:
268 F_s(r)~=&~F_\alpha(r) & r < r_1 \\
269 F_s(r)~=&~F_\alpha(r)+S(r) & r_1 \le r < r_c \\
270 F_s(r)~=&~0 & r_c \le r
273 When :math:`r_1=0` this is a traditional shift function, otherwise it
274 acts as a switch function. The corresponding shifted potential function
277 .. math:: V_s(r) = \int^{\infty}_r~F_s(x)\, dx
279 The |Gromacs| **force switch** function :math:`S_F(r)` should be smooth at
280 the boundaries, therefore the following boundary conditions are imposed
281 on the switch function:
288 S_F(r_c) &=&-F_\alpha(r_c) \\
289 S_F'(r_c) &=&-F_\alpha'(r_c)
292 A 3\ :math:`^{rd}` degree polynomial of the form
294 .. math:: S_F(r) = A(r-r_1)^2 + B(r-r_1)^3
296 fulfills these requirements. The constants A and B are given by the
297 boundary condition at :math:`r_c`:
302 A &~=~& -\alpha \, \displaystyle
303 \frac{(\alpha+4)r_c~-~(\alpha+1)r_1} {r_c^{\alpha+2}~(r_c-r_1)^2} \\
304 B &~=~& \alpha \, \displaystyle
305 \frac{(\alpha+3)r_c~-~(\alpha+1)r_1}{r_c^{\alpha+2}~(r_c-r_1)^3}
308 Thus the total force function is:
310 .. math:: F_s(r) = \frac{\alpha}{r^{\alpha+1}} + A(r-r_1)^2 + B(r-r_1)^3
312 and the potential function reads:
314 .. math:: V_s(r) = \frac{1}{r^\alpha} - \frac{A}{3} (r-r_1)^3 - \frac{B}{4} (r-r_1)^4 - C
318 .. math:: C = \frac{1}{r_c^\alpha} - \frac{A}{3} (r_c-r_1)^3 - \frac{B}{4} (r_c-r_1)^4
320 The |Gromacs| **potential-switch** function :math:`S_V(r)` scales the
321 potential between :math:`r_1` and :math:`r_c`, and has similar boundary
322 conditions, intended to produce smoothly-varying potential and forces:
335 The fifth-degree polynomial that has these properties is
337 .. 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}
339 This implementation is found in several other simulation
340 packages,\ :ref:`73 <refOhmine1988>`\ :ref:`75 <refGuenot1993>` but
341 differs from that in CHARMM.\ :ref:`76 <refSteinbach1994>` Switching the
342 potential leads to artificially large forces in the switching region,
343 therefore it is not recommended to switch Coulomb interactions using
344 this function,\ :ref:`72 <refSpoel2006a>` but switching Lennard-Jones
345 interactions using this function produces acceptable results.
347 Modified short-range interactions with Ewald summation
348 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
350 When Ewald summation or particle-mesh Ewald is used to calculate the
351 long-range interactions, the short-range Coulomb potential must also be
352 modified. Here the potential is switched to (nearly) zero at the
353 cut-off, instead of the force. In this case the short range potential is
356 .. math:: V(r) = f \frac{\mbox{erfc}(\beta r_{ij})}{r_{ij}} q_i q_j,
358 where :math:`\beta` is a parameter that determines the relative weight
359 between the direct space sum and the reciprocal space sum and
360 erfc\ :math:`(x)` is the complementary error function. For further
361 details on long-range electrostatics, see sec. :ref:`lrelstat`.