1Methods 2------- 3 4Exclusions and 1-4 Interactions. 5~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 6 7Atoms within a molecule that are close by in the chain, *i.e.* atoms 8that 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 12modeled by a Lennard-Jones potential. Instead it is assumed that these 13interactions are adequately modeled by a harmonic bond term or 14constraint (**i, i+1**) and a harmonic angle term (**i, i+2**). The 15first 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**. 18 19.. _fig-chain: 20 21.. figure:: plots/chain.* 22 :width: 8.00000cm 23 24 Atoms along an alkane chain. 25 26For third neighbors, the normal Lennard-Jones repulsion is sometimes 27still too strong, which means that when applied to a molecule, the 28molecule would deform or break due to the internal strain. This is 29especially the case for carbon-carbon interactions in a 30*cis*-conformation (*e.g.* *cis*-butane). Therefore, for some of these 31interactions, the Lennard-Jones repulsion has been reduced in the GROMOS 32force field, which is implemented by keeping a separate list of 1-4 and 33normal Lennard-Jones parameters. In other force fields, such as 34OPLS \ :ref:`103 <refJorgensen88>`, the standard Lennard-Jones 35parameters are reduced by a factor of two, but in that case also the 36dispersion (r\ :math:`^{-6}`) and the Coulomb interaction are scaled. 37|Gromacs| can use either of these methods. 38 39Charge Groups 40~~~~~~~~~~~~~ 41 42In principle, the force calculation in MD is an :math:`O(N^2)` problem. 43Therefore, we apply a cut-off for non-bonded force (NBF) calculations; 44only the particles within a certain distance of each other are 45interacting. This reduces the cost to :math:`O(N)` (typically 46:math:`100N` to :math:`200N`) of the NBF. It also introduces an error, 47which is, in most cases, acceptable, except when applying the cut-off 48implies the creation of charges, in which case you should consider using 49the lattice sum methods provided by |Gromacs|. 50 51Consider a water molecule interacting with another atom. If we would 52apply a plain cut-off on an atom-atom basis we might include the 53atom-oxygen interaction (with a charge of :math:`-0.82`) without the 54compensating charge of the protons, and as a result, induce a large 55dipole moment over the system. Therefore, we have to keep groups of 56atoms with total charge 0 together. These groups are called *charge 57groups*. Note that with a proper treatment of long-range electrostatics 58(e.g. particle-mesh Ewald (sec. :ref:`pme`), keeping charge groups 59together is not required. 60 61Treatment of Cut-offs in the group scheme 62~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 63 64|Gromacs| is quite flexible in treating cut-offs, which implies there can 65be quite a number of parameters to set. These parameters are set in the 66input file for grompp. There are two sort of parameters that affect the 67cut-off interactions; you can select which type of interaction to use in 68each case, and which cut-offs should be used in the neighbor searching. 69 70For both Coulomb and van der Waals interactions there are interaction 71type selectors (termed vdwtype and coulombtype) and two parameters, for 72a total of six non-bonded interaction parameters. See the User Guide for 73a complete description of these parameters. 74 75In 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 77radius at least as large as the :math:`r_c` specified for the functional 78form, because of the use of charge groups. The extra radius is typically 79of the order of 0.25 nm (roughly the largest distance between two atoms 80in a charge group plus the distance a charge group can diffuse within 81neighbor list updates). 82 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` 90 91.. _tab-funcparm: 92 93.. table:: Parameters for the different functional forms of the 94 non-bonded interactions. 95 96 +----------------------------+------------+ 97 | Type | Parameters | 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 +---------+------------------+------------+ 113 114.. _virtualsites: 115 116Virtual interaction sites 117------------------------- 118 119Virtual interaction sites (called dummy atoms in 120|Gromacs| versions before 3.3) can be used in |Gromacs| in a number of ways. 121We write the position of the virtual site :math:`\mathbf{r}_s` as a function 122of the positions of other particles 123:math:`\mathbf{r}`\ :math:`_i`: :math:`\mathbf{r}_s = 124f(\mathbf{r}_1..\mathbf{r}_n)`. The virtual site, which may carry charge or be 125involved in other interactions, can now be used in the force 126calculation. The force acting on the virtual site must be redistributed 127over the particles with mass in a consistent way. A good way to do this 128can be found in ref. \ :ref:`104 <refBerendsen84b>`. We can write the 129potential energy as: 130 131.. math:: V = V(\mathbf{r}_s,\mathbf{r}_1,\ldots,\mathbf{r}_n) = V^*(\mathbf{r}_1,\ldots,\mathbf{r}_n) 132 :label: eqnvsiteepot 133 134The force on the particle :math:`i` is then: 135 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 142 143The first term is the normal force. The second term is the force on 144particle :math:`i` due to the virtual site, which can be written in 145tensor notation: 146 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} 151 :label: eqnfvsite 152 153where :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 155virtual site. In this way, the total force and the total torque are 156conserved \ :ref:`104 <refBerendsen84b>`. 157 158The computation of the virial (:eq:`eqn. %s <eqnXi>`) is non-trivial when 159virtual sites are used. Since the virial involves a summation over all 160the atoms (rather than virtual sites), the forces must be redistributed 161from the virtual sites to the atoms (using :eq:`eqn. %s <eqnfvsite>`) *before* 162computation of the virial. In some special cases where the forces on the 163atoms can be written as a linear combination of the forces on the 164virtual sites (types 2 and 3 below) there is no difference between 165computing the virial before and after the redistribution of forces. 166However, in the general case redistribution should be done first. 167 168.. _fig-vsites: 169 170.. figure:: plots/dummies.* 171 :width: 15.00000cm 172 173 The seven different types of virtual site construction. The 174 constructing atoms are shown as black circles, the virtual sites in 175 gray. 176 177There are six ways to construct virtual sites from surrounding atoms in 178|Gromacs|, which we classify by the number of constructing atoms. **Note** 179that all site types mentioned can be constructed from types 3fd 180(normalized, in-plane) and 3out (non-normalized, out of plane). However, 181the amount of computation involved increases sharply along this list, so 182we strongly recommended using the first adequate virtual site type that 183will be sufficient for a certain purpose. :numref:`Fig. %s <fig-vsites>` depicts 6 of 184the available virtual site constructions. The conceptually simplest 185construction types are linear combinations: 186 187.. math:: \mathbf{r}_s = \sum_{i=1}^N w_i \, \mathbf{r}_i 188 :label: eqnvsitelincomb 189 190The force is then redistributed using the same weights: 191 192.. math:: \mathbf{F}_i = w_i \, \mathbf{F}_{s} 193 :label: eqnvsitelincombforce 194 195The types of virtual sites supported in |Gromacs| are given in the list 196below. Constructing atoms in virtual sites can be virtual sites 197themselves, but only if they are higher in the list, i.e. virtual sites 198can be constructed from “particles” that are simpler virtual sites. 199 200- On top of an atom. This allows giving an atom multiple atom types and 201 with that also assigned multiple, different bonded interactions. This 202 can espically be of use in free-energy calculations. 203 204- The coordinates of the virtual site equal that of the constructing atom: 205 206 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i 207 :label: eqnvsite1 208 209- The force is moved to the constructing atom: 210 211 .. math:: \mathbf{F}_i ~=~ \mathbf{F}_{s} 212 :label: eqnvsite1force 213 214- As a linear combination of two atoms 215 (:numref:`Fig. %s <fig-vsites>` 2): 216 217 .. math:: w_i = 1 - a ~,~~ w_j = a 218 :label: eqnvsitelin2atom 219 220- In this case the virtual site is on the line through atoms :math:`i` 221 and :math:`j`. 222 223- On the line through two atoms, with a fixed distance 224 (:numref:`Fig. %s <fig-vsites>` 2fd): 225 226 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \frac{ \mathbf{r}_{ij} } 227 { | \mathbf{r}_{ij} | } 228 :label: eqnvsite2fdatom 229 230- In this case the virtual site is on the line through the other two 231 particles at a distance of :math:`|a|` from :math:`i`. The force on 232 particles :math:`i` and :math:`j` due to the force on the virtual site 233 can be computed as: 234 235 .. math:: \begin{array}{lcr} 236 \mathbf{F}_i &=& \displaystyle \mathbf{F}_{s} - \gamma ( \mathbf{F}_{is} - \mathbf{p} ) \\[1ex] 237 \mathbf{F}_j &=& \displaystyle \gamma (\mathbf{F}_{s} - \mathbf{p}) \\[1ex] 238 \end{array} 239 ~\mbox{ where }~ 240 \begin{array}{c} 241 \displaystyle \gamma = \frac{a}{ | \mathbf{r}_{ij} | } \\[2ex] 242 \displaystyle \mathbf{p} = \frac{ \mathbf{r}_{is} \cdot \mathbf{F}_{s} } 243 { \mathbf{r}_{is} \cdot \mathbf{r}_{is} } \mathbf{r}_{is} 244 \end{array} 245 :label: eqnvsite2fdforce 246 247- As a linear combination of three atoms 248 (:numref:`Fig. %s <fig-vsites>` 3): 249 250 .. math:: w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b 251 :label: eqnvsitelin3atom 252 253- In this case the virtual site is in the plane of the other three 254 particles. 255 256- In the plane of three atoms, with a fixed distance 257 (:numref:`Fig. %s <fig-vsites>` 3fd): 258 259 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + b \frac{ (1 - a) \mathbf{r}_{ij} + a \mathbf{r}_{jk} } 260 { | (1 - a) \mathbf{r}_{ij} + a \mathbf{r}_{jk} | } 261 :label: eqnvsiteplane3atom 262 263- In this case the virtual site is in the plane of the other three 264 particles at a distance of :math:`|b|` from :math:`i`. The force on 265 particles :math:`i`, :math:`j` and :math:`k` due to the force on the 266 virtual site can be computed as: 267 268 .. math:: \begin{array}{lcr} 269 \mathbf{F}_i &=& \displaystyle \mathbf{F}_{s} - \gamma ( \mathbf{F}_{is} - \mathbf{p} ) \\[1ex] 270 \mathbf{F}_j &=& \displaystyle (1-a)\gamma (\mathbf{F}_{s} - \mathbf{p}) \\[1ex] 271 \mathbf{F}_k &=& \displaystyle a \gamma (\mathbf{F}_{s} - \mathbf{p}) \\ 272 \end{array} 273 ~\mbox{ where }~ 274 \begin{array}{c} 275 \displaystyle \gamma = \frac{b}{ | \mathbf{r}_{ij} + a \mathbf{r}_{jk} | } \\[2ex] 276 \displaystyle \mathbf{p} = \frac{ \mathbf{r}_{is} \cdot \mathbf{F}_{s} } 277 { \mathbf{r}_{is} \cdot \mathbf{r}_{is} } \mathbf{r}_{is} 278 \end{array} 279 :label: eqnvsiteplane3atomforce 280 281- In the plane of three atoms, with a fixed angle and 282 distance (:numref:`Fig. %s <fig-vsites>` 3fad): 283 284 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + 285 d \cos \theta \frac{\mathbf{r}_{ij}}{ | \mathbf{r}_{ij} | } + 286 d \sin \theta \frac{\mathbf{r}_\perp}{ | \mathbf{r}_\perp | } 287 ~\mbox{ where }~ 288 \mathbf{r}_\perp ~=~ \mathbf{r}_{jk} - 289 \frac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} } 290 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } 291 \mathbf{r}_{ij} 292 :label: eqnvsite2fadF 293 294- In this case the virtual site is in the plane of the other three 295 particles at a distance of :math:`|d|` from :math:`i` at an angle of 296 :math:`\alpha` with :math:`\mathbf{r}_{ij}`. Atom 297 :math:`k` defines the plane and the direction of the angle. **Note** 298 that in this case :math:`b` and :math:`\alpha` must be specified, 299 instead of :math:`a` and :math:`b` (see also sec. :ref:`vsitetop`). 300 The force on particles :math:`i`, :math:`j` and :math:`k` due to the 301 force on the virtual site can be computed as (with 302 :math:`\mathbf{r}_\perp` as defined in 303 :eq:`eqn. %s <eqnvsite2fadF>`): 304 305 .. math:: \begin{array}{c} 306 \begin{array}{lclllll} 307 \mathbf{F}_i &=& \mathbf{F}_{s} &-& 308 \dfrac{d \cos \theta}{ | \mathbf{r}_{ij} | } \mathbf{F}_1 &+& 309 \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \left( 310 \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} } 311 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{F}_2 + 312 \mathbf{F}_3 \right) \\[3ex] 313 \mathbf{F}_j &=& && 314 \dfrac{d \cos \theta}{ | \mathbf{r}_{ij} | } \mathbf{F}_1 &-& 315 \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \left( 316 \mathbf{F}_2 + 317 \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} } 318 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{F}_2 + 319 \mathbf{F}_3 \right) \\[3ex] 320 \mathbf{F}_k &=& && && 321 \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \mathbf{F}_2 \\[3ex] 322 \end{array} \\[5ex] 323 ~\mbox{where }~ 324 \mathbf{F}_1 = \mathbf{F}_{s} - 325 \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{F}_{s} } 326 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{r}_{ij} 327 ~\mbox{, }~ 328 \mathbf{F}_2 = \mathbf{F}_1 - 329 \dfrac{ \mathbf{r}_\perp \cdot \mathbf{F}_{s} } 330 { \mathbf{r}_\perp \cdot \mathbf{r}_\perp } \mathbf{r}_\perp 331 ~\mbox{and }~ 332 \mathbf{F}_3 = \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{F}_{s} } 333 { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{r}_\perp 334 \end{array} 335 :label: eqnvsite2fadFforce 336 337- As a non-linear combination of three atoms, out of 338 plane (:numref:`Fig. %s <fig-vsites>` 3out): 339 340 .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \mathbf{r}_{ij} + b \mathbf{r}_{ik} + 341 c (\mathbf{r}_{ij} \times \mathbf{r}_{ik}) 342 :label: eqnvsitenonlin3atom 343 344- This enables the construction of virtual sites out of the plane of 345 the other atoms. The force on particles :math:`i,j` and :math:`k` due 346 to the force on the virtual site can be computed as: 347 348 .. math:: \begin{array}{lcl} 349 \mathbf{F}_j &=& \left[\begin{array}{ccc} 350 a & -c\,z_{ik} & c\,y_{ik} \\[0.5ex] 351 c\,z_{ik} & a & -c\,x_{ik} \\[0.5ex] 352 -c\,y_{ik} & c\,x_{ik} & a 353 \end{array}\right]\mathbf{F}_{s} \\ 354 \mathbf{F}_k &=& \left[\begin{array}{ccc} 355 b & c\,z_{ij} & -c\,y_{ij} \\[0.5ex] 356 -c\,z_{ij} & b & c\,x_{ij} \\[0.5ex] 357 c\,y_{ij} & -c\,x_{ij} & b 358 \end{array}\right]\mathbf{F}_{s} \\ 359 \mathbf{F}_i &=& \mathbf{F}_{s} - \mathbf{F}_j - \mathbf{F}_k 360 \end{array} 361 :label: eqnvsitenonlin3atomforce 362 363- From four atoms, with a fixed distance, see 364 separate :numref:`Fig. %s <fig-vsite4fdn>`. This construction is a bit complex, 365 in particular since the previous type (4fd) could be unstable which 366 forced us to introduce a more elaborate construction: 367 368.. _fig-vsite4fdn: 369 370.. figure:: plots/vsite-4fdn.* 371 :width: 5.00000cm 372 373 The new 4fdn virtual site construction, which is stable even when 374 all constructing atoms are in the same plane. 375 376- 377 .. math:: \begin{aligned} 378 \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\ 379 \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\ 380 \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\ 381 \mathbf{x}_s &=& \mathbf{x}_i + c \frac{\mathbf{r}_m}{ | \mathbf{r}_m | } 382 \end{aligned} 383 :label: eqnvsite 384 385- In this case the virtual site is at a distance of :math:`|c|` from 386 :math:`i`, while :math:`a` and :math:`b` are parameters. **Note** 387 that the vectors :math:`\mathbf{r}_{ik}` and :math:`\mathbf{r}_{ij}` 388 are not normalized to save floating-point operations. The force on 389 particles :math:`i`, :math:`j`, :math:`k` and :math:`l` due to the 390 force on the virtual site are computed through chain rule derivatives 391 of the construction expression. This is exact and conserves energy, 392 but it does lead to relatively lengthy expressions that we do not 393 include here (over 200 floating-point operations). The interested 394 reader can look at the source code in ``vsite.c``. Fortunately, this 395 vsite type is normally only used for chiral centers such as 396 :math:`C_{\alpha}` atoms in proteins. 397 398 The new 4fdn construct is identified with a ‘type’ value of 2 in the 399 topology. The earlier 4fd type is still supported internally (‘type’ 400 value 1), but it should not be used for new simulations. All current 401 |Gromacs| tools will automatically generate type 4fdn instead. 402 403- A linear combination of :math:`N` atoms with relative 404 weights :math:`a_i`. The weight for atom :math:`i` is: 405 406 .. math:: w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1} 407 :label: eqnvsiterelweight 408 409- There are three options for setting the weights: 410 411 - center of geometry: equal weights 412 413 - center of mass: :math:`a_i` is the mass of atom :math:`i`; when in 414 free-energy simulations the mass of the atom is changed, only the 415 mass of the A-state is used for the weight 416 417 - center of weights: :math:`a_i` is defined by the user 418 419