4 Exclusions and 1-4 Interactions.
5 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
7 Atoms within a molecule that are close by in the chain, *i.e.* atoms
8 that are covalently bonded, or linked by one or two atoms are called
9 *first neighbors, second neighbors* and *third neighbors*, respectively
10 (see :numref:`Fig. %s <fig-chain>`). Since the interactions of atom **i** with atoms
11 **i+1** and **i+2** are mainly quantum mechanical, they can not be
12 modeled by a Lennard-Jones potential. Instead it is assumed that these
13 interactions are adequately modeled by a harmonic bond term or
14 constraint (**i, i+1**) and a harmonic angle term (**i, i+2**). The
15 first and second neighbors (atoms **i+1** and **i+2**) are therefore
16 *excluded* from the Lennard-Jones interaction list of atom **i**; atoms
17 **i+1** and **i+2** are called *exclusions* of atom **i**.
21 .. figure:: plots/chain.*
24 Atoms along an alkane chain.
26 For third neighbors, the normal Lennard-Jones repulsion is sometimes
27 still too strong, which means that when applied to a molecule, the
28 molecule would deform or break due to the internal strain. This is
29 especially the case for carbon-carbon interactions in a
30 *cis*-conformation (*e.g.* *cis*-butane). Therefore, for some of these
31 interactions, the Lennard-Jones repulsion has been reduced in the GROMOS
32 force field, which is implemented by keeping a separate list of 1-4 and
33 normal Lennard-Jones parameters. In other force fields, such as
34 OPLS \ :ref:`103 <refJorgensen88>`, the standard Lennard-Jones
35 parameters are reduced by a factor of two, but in that case also the
36 dispersion (r\ :math:`^{-6}`) and the Coulomb interaction are scaled.
37 |Gromacs| can use either of these methods.
42 In principle, the force calculation in MD is an :math:`O(N^2)` problem.
43 Therefore, we apply a cut-off for non-bonded force (NBF) calculations;
44 only the particles within a certain distance of each other are
45 interacting. This reduces the cost to :math:`O(N)` (typically
46 :math:`100N` to :math:`200N`) of the NBF. It also introduces an error,
47 which is, in most cases, acceptable, except when applying the cut-off
48 implies the creation of charges, in which case you should consider using
49 the lattice sum methods provided by |Gromacs|.
51 Consider a water molecule interacting with another atom. If we would
52 apply a plain cut-off on an atom-atom basis we might include the
53 atom-oxygen interaction (with a charge of :math:`-0.82`) without the
54 compensating charge of the protons, and as a result, induce a large
55 dipole moment over the system. Therefore, we have to keep groups of
56 atoms with total charge 0 together. These groups are called *charge
57 groups*. Note that with a proper treatment of long-range electrostatics
58 (e.g. particle-mesh Ewald (sec. :ref:`pme`), keeping charge groups
59 together is not required.
61 Treatment of Cut-offs in the group scheme
62 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
64 |Gromacs| is quite flexible in treating cut-offs, which implies there can
65 be quite a number of parameters to set. These parameters are set in the
66 input file for grompp. There are two sort of parameters that affect the
67 cut-off interactions; you can select which type of interaction to use in
68 each case, and which cut-offs should be used in the neighbor searching.
70 For both Coulomb and van der Waals interactions there are interaction
71 type selectors (termed vdwtype and coulombtype) and two parameters, for
72 a total of six non-bonded interaction parameters. See the User Guide for
73 a complete description of these parameters.
75 In the group cut-off scheme, all of the interaction functions in
76 :numref:`Table %s <tab-funcparm>` require that neighbor searching be done with a
77 radius at least as large as the :math:`r_c` specified for the functional
78 form, because of the use of charge groups. The extra radius is typically
79 of the order of 0.25 nm (roughly the largest distance between two atoms
80 in a charge group plus the distance a charge group can diffuse within
81 neighbor list updates).
83 .. |CPCOP| replace:: :math:`r_c`, :math:`{\varepsilon}_{r}`
84 .. |CRFP| replace:: :math:`r_c`, :math:`{\varepsilon}_{rf}`
85 .. |CSHFP| replace:: :math:`r_1`, :math:`r_c`, :math:`{\varepsilon}_{r}`
86 .. |CSWFP| replace:: :math:`r_1`, :math:`r_c`, :math:`{\varepsilon}_{r}`
87 .. |VPCOP| replace:: :math:`r_c`
88 .. |VSHFP| replace:: :math:`r_1`, :math:`r_c`
89 .. |VSWFP| replace:: :math:`r_1`, :math:`r_c`
93 .. table:: Parameters for the different functional forms of the
94 non-bonded interactions.
96 +----------------------------+------------+
98 +=========+==================+============+
99 | Coulomb | Plain cut-off | |CPCOP| |
100 | +------------------+------------+
101 | | Reaction field | |CRFP| |
102 | +------------------+------------+
103 | | Shift function | |CSHFP| |
104 | +------------------+------------+
105 | | Switch function | |CSWFP| |
106 +---------+------------------+------------+
107 | VdW | Plain cut-off | |VPCOP| |
108 | +------------------+------------+
109 | | Shift function | |VSHFP| |
110 | +------------------+------------+
111 | | Switch function | |VSWFP| |
112 +---------+------------------+------------+
116 Virtual interaction sites
117 -------------------------
119 Virtual interaction sites (called dummy atoms in
120 |Gromacs| versions before 3.3) can be used in |Gromacs| in a number of ways.
121 We write the position of the virtual site :math:`\mathbf{r}_s` as a function
122 of the positions of other particles
123 :math:`\mathbf{r}`\ :math:`_i`: :math:`\mathbf{r}_s =
124 f(\mathbf{r}_1..\mathbf{r}_n)`. The virtual site, which may carry charge or be
125 involved in other interactions, can now be used in the force
126 calculation. The force acting on the virtual site must be redistributed
127 over the particles with mass in a consistent way. A good way to do this
128 can be found in ref. \ :ref:`104 <refBerendsen84b>`. We can write the
131 .. math:: V = V(\mathbf{r}_s,\mathbf{r}_1,\ldots,\mathbf{r}_n) = V^*(\mathbf{r}_1,\ldots,\mathbf{r}_n)
134 The force on the particle :math:`i` is then:
136 .. math:: \mathbf{F}_i = -\frac{\partial V^*}{\partial \mathbf{r}_i}
137 = -\frac{\partial V}{\partial \mathbf{r}_i} -
138 \frac{\partial V}{\partial \mathbf{r}_s}
139 \frac{\partial \mathbf{r}_s}{\partial \mathbf{r}_i}
140 = \mathbf{F}_i^{direct} + \mathbf{F}_i
141 :label: eqnvsiteforce
143 The first term is the normal force. The second term is the force on
144 particle :math:`i` due to the virtual site, which can be written in
147 .. math:: \mathbf{F}_i = \left[\begin{array}{ccc}
148 {\displaystyle\frac{\partial x_s}{\partial x_i}} & {\displaystyle\frac{\partial y_s}{\partial x_i}} & {\displaystyle\frac{\partial z_s}{\partial x_i}} \\[1ex]
149 {\displaystyle\frac{\partial x_s}{\partial y_i}} & {\displaystyle\frac{\partial y_s}{\partial y_i}} & {\displaystyle\frac{\partial z_s}{\partial y_i}} \\[1ex]
150 {\displaystyle\frac{\partial x_s}{\partial z_i}} & {\displaystyle\frac{\partial y_s}{\partial z_i}} & {\displaystyle\frac{\partial z_s}{\partial z_i}} \end{array}\right]\mathbf{F}_{s}
153 where :math:`\mathbf{F}_{s}` is the force on the virtual site and
154 :math:`x_s`, :math:`y_s` and :math:`z_s` are the coordinates of the
155 virtual site. In this way, the total force and the total torque are
156 conserved \ :ref:`104 <refBerendsen84b>`.
158 The computation of the virial (:eq:`eqn. %s <eqnXi>`) is non-trivial when
159 virtual sites are used. Since the virial involves a summation over all
160 the atoms (rather than virtual sites), the forces must be redistributed
161 from the virtual sites to the atoms (using :eq:`eqn. %s <eqnfvsite>`) *before*
162 computation of the virial. In some special cases where the forces on the
163 atoms can be written as a linear combination of the forces on the
164 virtual sites (types 2 and 3 below) there is no difference between
165 computing the virial before and after the redistribution of forces.
166 However, in the general case redistribution should be done first.
170 .. figure:: plots/dummies.*
173 The seven different types of virtual site construction. The
174 constructing atoms are shown as black circles, the virtual sites in
177 There are six ways to construct virtual sites from surrounding atoms in
178 |Gromacs|, which we classify by the number of constructing atoms. **Note**
179 that all site types mentioned can be constructed from types 3fd
180 (normalized, in-plane) and 3out (non-normalized, out of plane). However,
181 the amount of computation involved increases sharply along this list, so
182 we strongly recommended using the first adequate virtual site type that
183 will be sufficient for a certain purpose. :numref:`Fig. %s <fig-vsites>` depicts 6 of
184 the available virtual site constructions. The conceptually simplest
185 construction types are linear combinations:
187 .. math:: \mathbf{r}_s = \sum_{i=1}^N w_i \, \mathbf{r}_i
188 :label: eqnvsitelincomb
190 The force is then redistributed using the same weights:
192 .. math:: \mathbf{F}_i = w_i \, \mathbf{F}_{s}
193 :label: eqnvsitelincombforce
195 The types of virtual sites supported in |Gromacs| are given in the list
196 below. Constructing atoms in virtual sites can be virtual sites
197 themselves, but only if they are higher in the list, i.e. virtual sites
198 can be constructed from “particles” that are simpler virtual sites. The
199 virtual site velocities are reported, but not used in the integration
200 of the virtual site positions.
205 - This allows giving an atom multiple atom types and
206 with that also assigned multiple, different bonded interactions. This
207 can especially be of use in free-energy calculations.
209 - The coordinates of the virtual site equal that of the constructing atom:
211 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i
214 - The force is moved to the constructing atom:
216 .. math:: \mathbf{F}_i ~=~ \mathbf{F}_{s}
217 :label: eqnvsite1force
219 - The velocity of the virtual site equals that of the constructing atom:
221 .. math:: \mathbf{v}_s ~=~ \mathbf{v}_i
224 As a linear combination of two atoms (:numref:`Fig. %s <fig-vsites>` 2)
225 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
227 - The weights are calculated as
229 .. math:: w_i = 1 - a ~,~~ w_j = a
230 :label: eqnvsitelin2atom
232 - In this case the virtual site is on the line through atoms :math:`i`
235 - The velocity of the virtual site is a linear combination of the
236 velocities of the constructing atoms
238 On the line through two atoms, with a fixed distance (:numref:`Fig. %s <fig-vsites>` 2fd)
239 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
241 - The position is calculated as:
243 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \frac{ \mathbf{r}_{ij} }
244 { | \mathbf{r}_{ij} | }
245 :label: eqnvsite2fdatom
247 - In this case the virtual site is on the line through the other two
248 particles at a distance of :math:`|a|` from :math:`i`. The force on
249 particles :math:`i` and :math:`j` due to the force on the virtual site
252 .. math:: \begin{array}{lcr}
253 \mathbf{F}_i &=& \displaystyle \mathbf{F}_{s} - \gamma ( \mathbf{F}_{is} - \mathbf{p} ) \\[1ex]
254 \mathbf{F}_j &=& \displaystyle \gamma (\mathbf{F}_{s} - \mathbf{p}) \\[1ex]
258 \displaystyle \gamma = \frac{a}{ | \mathbf{r}_{ij} | } \\[2ex]
259 \displaystyle \mathbf{p} = \frac{ \mathbf{r}_{is} \cdot \mathbf{F}_{s} }
260 { \mathbf{r}_{is} \cdot \mathbf{r}_{is} } \mathbf{r}_{is}
262 :label: eqnvsite2fdforce
264 - The velocity is calculated as:
266 .. math:: \mathbf{v}_{s} = \mathbf{v}_{i} + \frac{a}{|\mathbf{r}_{ij}|}
267 \left(\mathbf{v}_{ij} - \mathbf{r}_{ij}
268 \frac{\mathbf{v}_{ij}\cdot\mathbf{r}_{ij}}
269 {|\mathbf{r}_{ij}|^2}\right)
270 :label: eqnvsite2fdatomvel
272 As a linear combination of three atoms (:numref:`Fig. %s <fig-vsites>` 3)
273 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
275 - The weights are calculated as:
277 .. math:: w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b
278 :label: eqnvsitelin3atom
280 - In this case the virtual site is in the plane of the other three
283 In the plane of three atoms, with a fixed distance (:numref:`Fig. %s <fig-vsites>` 3fd)
284 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
286 - The position is calculated as:
288 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + b \frac{ \mathbf{r}_{ijk} } { | \mathbf{r}_{ijk} | }
290 \mathbf{r}_{ijk} ~=~ (1 - a) \mathbf{r}_{ij} + a \mathbf{r}_{jk}
291 :label: eqnvsiteplane3atom
293 - In this case the virtual site is in the plane of the other three
294 particles at a distance of :math:`|b|` from :math:`i`. The force on
295 particles :math:`i`, :math:`j` and :math:`k` due to the force on the
296 virtual site can be computed as:
298 .. math:: \begin{array}{lcr}
299 \mathbf{F}_i &=& \displaystyle \mathbf{F}_{s} - \gamma ( \mathbf{F}_{is} - \mathbf{p} ) \\[1ex]
300 \mathbf{F}_j &=& \displaystyle (1-a)\gamma (\mathbf{F}_{s} - \mathbf{p}) \\[1ex]
301 \mathbf{F}_k &=& \displaystyle a \gamma (\mathbf{F}_{s} - \mathbf{p}) \\
305 \displaystyle \gamma = \frac{b}{ | \mathbf{r}_{ij} + a \mathbf{r}_{jk} | } \\[2ex]
306 \displaystyle \mathbf{p} = \frac{ \mathbf{r}_{is} \cdot \mathbf{F}_{s} }
307 { \mathbf{r}_{is} \cdot \mathbf{r}_{is} } \mathbf{r}_{is}
309 :label: eqnvsiteplane3atomforce
311 - The velocity is calculated as:
313 .. math:: \mathbf{v}_{s} ~=~ \mathbf{v}_{i} +
314 \frac{b}{|\mathbf{r}_{ijk}|}
315 \left(\dot{\mathbf{r}}_{ijk} -
316 \mathbf{r}_{ijk}\frac{\dot{\mathbf{r}}_{ijk}\cdot\mathbf{r}_{ijk}}
317 {|\mathbf{r}_{ijk}|^2}\right)
318 :label: eqnvsiteplane3atomvel
320 In the plane of three atoms, with a fixed angle and distance (:numref:`Fig. %s <fig-vsites>` 3fad)
321 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
323 - The position is calculated as:
325 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i +
326 d \cos \theta \frac{\mathbf{r}_{ij}}{ | \mathbf{r}_{ij} | } +
327 d \sin \theta \frac{\mathbf{r}_\perp}{ | \mathbf{r}_\perp | }
329 \mathbf{r}_\perp ~=~ \mathbf{r}_{jk} -
330 \frac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} }
331 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} }
333 :label: eqnvsite2fadF
335 - In this case the virtual site is in the plane of the other three
336 particles at a distance of :math:`|d|` from :math:`i` at an angle of
337 :math:`\alpha` with :math:`\mathbf{r}_{ij}`. Atom
338 :math:`k` defines the plane and the direction of the angle. **Note**
339 that in this case :math:`b` and :math:`\alpha` must be specified,
340 instead of :math:`a` and :math:`b` (see also sec. :ref:`vsitetop`).
341 The force on particles :math:`i`, :math:`j` and :math:`k` due to the
342 force on the virtual site can be computed as (with
343 :math:`\mathbf{r}_\perp` as defined in
344 :eq:`eqn. %s <eqnvsite2fadF>`):
346 .. math:: \begin{array}{c}
347 \begin{array}{lclllll}
348 \mathbf{F}_i &=& \mathbf{F}_{s} &-&
349 \dfrac{d \cos \theta}{ | \mathbf{r}_{ij} | } \mathbf{F}_1 &+&
350 \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \left(
351 \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} }
352 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{F}_2 +
353 \mathbf{F}_3 \right) \\[3ex]
355 \dfrac{d \cos \theta}{ | \mathbf{r}_{ij} | } \mathbf{F}_1 &-&
356 \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \left(
358 \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} }
359 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{F}_2 +
360 \mathbf{F}_3 \right) \\[3ex]
361 \mathbf{F}_k &=& && &&
362 \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \mathbf{F}_2 \\[3ex]
365 \mathbf{F}_1 = \mathbf{F}_{s} -
366 \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{F}_{s} }
367 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{r}_{ij}
369 \mathbf{F}_2 = \mathbf{F}_1 -
370 \dfrac{ \mathbf{r}_\perp \cdot \mathbf{F}_{s} }
371 { \mathbf{r}_\perp \cdot \mathbf{r}_\perp } \mathbf{r}_\perp
373 \mathbf{F}_3 = \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{F}_{s} }
374 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{r}_\perp
376 :label: eqnvsite2fadFforce
378 - The velocity is calculated as:
380 .. math:: \mathbf{v}_{s} &= \mathbf{v}_{i} + d\cos\theta\ \frac{\delta}{\delta t}\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} +
381 d\sin\theta\ \frac{\delta}{\delta t}\frac{\mathbf{r}_{\perp}}{|\mathbf{r}_{\perp}|} \\
383 \frac{\delta}{\delta t}\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} &=
384 \frac{1}{|\mathbf{r}_{ij}|}\left(\mathbf{v}_{ij} - \mathbf{r}_{ij}
385 \frac{\mathbf{v}_{ij}\cdot\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|^2}\right)\\
386 \frac{\delta}{\delta t}\frac{\mathbf{r}_{\perp}}{|\mathbf{r}_{\perp}|} &=
387 \frac{1}{|\mathbf{r}_{\perp}|}
388 \left(\dot{\mathbf{r}}_{\perp} - \mathbf{r}_{\perp}\frac{\dot{\mathbf{r}}_{\perp}\cdot\mathbf{r}_{\perp}}{|\mathbf{r}_{\perp}|^2}\right)\\
389 \dot{\mathbf{r}}_\perp &= \mathbf{v}_{jk} - \mathbf{r}_{ij}
390 \frac{|\mathbf{r}_{ij}|^2(\mathbf{v}_{ij}\cdot\mathbf{r}_{jk} + \mathbf{r}_{ij}\cdot\mathbf{v}_{jk}) -
391 (\mathbf{r}_{ij}\cdot\mathbf{r}_{jk})(2\mathbf{r}_{ij}\cdot\mathbf{v}_{ij})} {|\mathbf{r}_{ij}|^4} -
392 \frac{\mathbf{r}_{ij}\cdot\mathbf{r}_{jk}}{|\mathbf{r}_{ij}|^2}\ \mathbf{v}_{ij}
393 :label: eqnvsite2fadvel
395 As a non-linear combination of three atoms, out of plane (:numref:`Fig. %s <fig-vsites>` 3out)
396 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
398 - The position is calculated as:
400 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \mathbf{r}_{ij} + b \mathbf{r}_{ik} +
401 c (\mathbf{r}_{ij} \times \mathbf{r}_{ik})
402 :label: eqnvsitenonlin3atom
404 - This enables the construction of virtual sites out of the plane of
405 the other atoms. The force on particles :math:`i,j` and :math:`k` due
406 to the force on the virtual site can be computed as:
408 .. math:: \begin{array}{lcl}
409 \mathbf{F}_j &=& \left[\begin{array}{ccc}
410 a & -c\,z_{ik} & c\,y_{ik} \\[0.5ex]
411 c\,z_{ik} & a & -c\,x_{ik} \\[0.5ex]
412 -c\,y_{ik} & c\,x_{ik} & a
413 \end{array}\right]\mathbf{F}_{s} \\
414 \mathbf{F}_k &=& \left[\begin{array}{ccc}
415 b & c\,z_{ij} & -c\,y_{ij} \\[0.5ex]
416 -c\,z_{ij} & b & c\,x_{ij} \\[0.5ex]
417 c\,y_{ij} & -c\,x_{ij} & b
418 \end{array}\right]\mathbf{F}_{s} \\
419 \mathbf{F}_i &=& \mathbf{F}_{s} - \mathbf{F}_j - \mathbf{F}_k
421 :label: eqnvsitenonlin3atomforce
423 - The velocity is calculated as:
425 .. math:: \mathbf{v}_{s} ~=~ \mathbf{v}_{i} + \frac{c}{|\mathbf{r}_{m}|}\left(\dot{\mathbf{r}}_{m} -
426 \mathbf{r}_{m} \frac{\dot{\mathbf{r}}_{m}\cdot\mathbf{r}_{m}}{|\mathbf{r}_{m}|^2}\right)
427 :label: eqnvsitenonlin3atomvel
429 From four atoms, with a fixed distance, see separate :numref:`Fig. %s <fig-vsite4fdn>`
430 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
431 - This construction is a bit complex,
432 in particular since the previous type (4fd) could be unstable which
433 forced us to introduce a more elaborate construction:
437 .. figure:: plots/vsite-4fdn.*
440 The new 4fdn virtual site construction, which is stable even when
441 all constructing atoms are in the same plane.
443 - The position is calculated as
445 .. math:: \begin{aligned}
446 \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
447 \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
448 \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\
449 \mathbf{r}_s &=& \mathbf{r}_i + c \frac{\mathbf{r}_m}{ | \mathbf{r}_m | }
453 - The velocity is calculated as:
455 .. math:: \mathbf{v}_{s} = \mathbf{v}_{i} + \frac{c}{|\mathbf{r}_{m}|}\left(\dot{\mathbf{r}}_{m} - \mathbf{r}_{m} \frac{\dot{\mathbf{r}}_{m}\cdot\mathbf{r}_{m}}{|\mathbf{r}_{m}|^2}\right)\\
457 \dot{\mathbf{r}}_{m} &= \dot{\mathbf{r}}_{ja} \times \mathbf{r}_{jb} + \mathbf{r}_{ja} \times \dot{\mathbf{r}}_{jb}
460 - In this case the virtual site is at a distance of :math:`|c|` from
461 :math:`i`, while :math:`a` and :math:`b` are parameters. **Note**
462 that the vectors :math:`\mathbf{r}_{ik}` and :math:`\mathbf{r}_{ij}`
463 are not normalized to save floating-point operations. The force on
464 particles :math:`i`, :math:`j`, :math:`k` and :math:`l` due to the
465 force on the virtual site are computed through chain rule derivatives
466 of the construction expression. This is exact and conserves energy,
467 but it does lead to relatively lengthy expressions that we do not
468 include here (over 200 floating-point operations). The interested
469 reader can look at the source code in ``vsite.c``. Fortunately, this
470 vsite type is normally only used for chiral centers such as
471 :math:`C_{\alpha}` atoms in proteins.
473 The new 4fdn construct is identified with a ‘type’ value of 2 in the
474 topology. The earlier 4fd type is still supported internally (‘type’
475 value 1), but it should not be used for new simulations. All current
476 |Gromacs| tools will automatically generate type 4fdn instead.
478 A linear combination of :math:`N` atoms with relative weights :math:`a_i`
479 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
481 - The weight for atom :math:`i` is:
483 .. math:: w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1}
484 :label: eqnvsiterelweight
486 - There are three options for setting the weights:
488 - center of geometry: equal weights
490 - center of mass: :math:`a_i` is the mass of atom :math:`i`; when in
491 free-energy simulations the mass of the atom is changed, only the
492 mass of the A-state is used for the weight
494 - center of weights: :math:`a_i` is defined by the user