2 % This file is part of the GROMACS molecular simulation package.
4 % Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
5 % Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 % and including many others, as listed in the AUTHORS file in the
7 % top-level source directory and at http://www.gromacs.org.
9 % GROMACS is free software; you can redistribute it and/or
10 % modify it under the terms of the GNU Lesser General Public License
11 % as published by the Free Software Foundation; either version 2.1
12 % of the License, or (at your option) any later version.
14 % GROMACS is distributed in the hope that it will be useful,
15 % but WITHOUT ANY WARRANTY; without even the implied warranty of
16 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 % Lesser General Public License for more details.
19 % You should have received a copy of the GNU Lesser General Public
20 % License along with GROMACS; if not, see
21 % http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 % Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 % If you want to redistribute modifications to GROMACS, please
25 % consider that scientific software is very special. Version
26 % control is crucial - bugs must be traceable. We will be happy to
27 % consider code for inclusion in the official distribution, but
28 % derived work must not be called official GROMACS. Details are found
29 % in the README & COPYING files - if they are missing, get the
30 % official version at http://www.gromacs.org.
32 % To help us fund GROMACS development, we humbly ask that you cite
33 % the research papers on the package. Check out http://www.gromacs.org.
35 \chapter{Special Topics}
39 \section{Free energy implementation}
40 \label{sec:dgimplement}
41 For free energy calculations, there are two things that must be
42 specified; the end states, and the pathway connecting the end states.
43 The end states can be specified in two ways. The most straightforward
44 is through the specification of end states in the topology file. Most
45 potential forms support both an $A$ state and a $B$ state. Whenever both
46 states are specified, then the $A$ state corresponds to the initial free
47 energy state, and the $B$ state corresponds to the final state.
49 In some cases, the end state can also be defined in some cases without
50 altering the topology, solely through the {\tt .mdp} file, through the use
51 of the {\tt couple-moltype},{\tt couple-lambda0}, {\tt couple-lambda1}, and
52 {\tt couple-intramol} mdp keywords. Any molecule type selected in
53 {\tt couple-moltype} will automatically have a $B$ state implicitly
54 constructed (and the $A$ state redefined) according to the {\tt couple-lambda}
55 keywords. {\tt couple-lambda0} and {\tt couple-lambda1} define the non-bonded
56 parameters that are present in the $A$ state ({\tt couple-lambda0})
57 and the $B$ state ({\tt couple-lambda1}). The choices are 'q','vdw', and
58 'vdw-q'; these indicate the Coulombic, van der Waals, or both parameters
59 that are turned on in the respective state.
61 Once the end states are defined, then the path between the end states
62 has to be defined. This path is defined solely in the .mdp file.
63 Starting in 4.6, $\lambda$ is a vector of components, with Coulombic,
64 van der Waals, bonded, restraint, and mass components all able to be
65 adjusted independently. This makes it possible to turn off the
66 Coulombic term linearly, and then the van der Waals using soft core,
67 all in the same simulation. This is especially useful for replica
68 exchange or expanded ensemble simulations, where it is important to
69 sample all the way from interacting to non-interacting states in the
70 same simulation to improve sampling.
72 {\tt fep-lambdas} is the default array of $\lambda$ values ranging
73 from 0 to 1. All of the other lambda arrays use the values in this
74 array if they are not specified. The previous behavior, where the
75 pathway is controlled by a single $\lambda$ variable, can be preserved
76 by using only {\tt fep-lambdas} to define the pathway.
78 For example, if you wanted to first to change the Coulombic terms,
79 then the van der Waals terms, changing bonded at the same time rate as
80 the van der Waals, but changing the restraints throughout the first
81 two-thirds of the simulation, then you could use this $\lambda$ vector:
84 coul-lambdas = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
85 vdw-lambdas = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
86 bonded-lambdas = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
87 restraint-lambdas = 0.0 0.0 0.1 0.2 0.3 0.5 0.7 1.0 1.0 1.0
90 This is also equivalent to:
93 fep-lambdas = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
94 coul-lambdas = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
95 restraint-lambdas = 0.0 0.0 0.1 0.2 0.3 0.5 0.7 1.0 1.0 1.0
97 The {\tt fep-lambda array}, in this case, is being used as the default to
98 fill in the bonded and van der Waals $\lambda$ arrays. Usually, it's best to fill
99 in all arrays explicitly, just to make sure things are properly
102 If you want to turn on only restraints going from $A$ to $B$, then it would be:
104 restraint-lambdas = 0.0 0.1 0.2 0.4 0.6 1.0
106 and all of the other components of the $\lambda$ vector would be left in the $A$ state.
108 To compute free energies with a vector $\lambda$ using
109 thermodynamic integration, then the TI equation becomes vector equation:
111 \Delta F = \int \langle \nabla H \rangle \cdot d\vec{\lambda}
113 or for finite differences:
115 \Delta F \approx \int \sum \langle \nabla H \rangle \cdot \Delta\lambda
118 The external {\tt pymbar} script downloaded from https://SimTK.org/home/pymbar can
119 compute this integral automatically from the {\gromacs} dhdl.xvg output.
121 \section{Potential of mean force}
123 A potential of mean force (PMF) is a potential that is obtained
124 by integrating the mean force from an ensemble of configurations.
125 In {\gromacs}, there are several different methods to calculate the mean force.
126 Each method has its limitations, which are listed below.
128 \item{\bf pull code:} between the centers of mass of molecules or groups of molecules.
129 \item{\bf free-energy code with harmonic bonds or constraints:} between single atoms.
130 \item{\bf free-energy code with position restraints:} changing the conformation of a relatively immobile group of atoms.
131 \item{\bf pull code in limited cases:} between groups of atoms that are
132 part of a larger molecule for which the bonds are constrained with
133 SHAKE or LINCS. If the pull group if relatively large,
134 the pull code can be used.
136 The pull and free-energy code a described in more detail
137 in the following two sections.
139 \subsubsection{Entropic effects}
140 When a distance between two atoms or the centers of mass of two groups
141 is constrained or restrained, there will be a purely entropic contribution
142 to the PMF due to the rotation of the two groups~\cite{RMNeumann1980a}.
143 For a system of two non-interacting masses the potential of mean force is:
145 V_{pmf}(r) = -(n_c - 1) k_B T \log(r)
147 where $n_c$ is the number of dimensions in which the constraint works
148 (i.e. $n_c=3$ for a normal constraint and $n_c=1$ when only
149 the $z$-direction is constrained).
150 Whether one needs to correct for this contribution depends on what
151 the PMF should represent. When one wants to pull a substrate
152 into a protein, this entropic term indeed contributes to the work to
153 get the substrate into the protein. But when calculating a PMF
154 between two solutes in a solvent, for the purpose of simulating
155 without solvent, the entropic contribution should be removed.
156 {\bf Note} that this term can be significant; when at 300K the distance is halved,
157 the contribution is 3.5 kJ~mol$^{-1}$.
159 \section{Non-equilibrium pulling}
160 When the distance between two groups is changed continuously,
161 work is applied to the system, which means that the system is no longer
162 in equilibrium. Although in the limit of very slow pulling
163 the system is again in equilibrium, for many systems this limit
164 is not reachable within reasonable computational time.
165 However, one can use the Jarzynski relation~\cite{Jarzynski1997a}
166 to obtain the equilibrium free-energy difference $\Delta G$
167 between two distances from many non-equilibrium simulations:
169 \Delta G_{AB} = -k_BT \log \left\langle e^{-\beta W_{AB}} \right\rangle_A
172 where $W_{AB}$ is the work performed to force the system along one path
173 from state A to B, the angular bracket denotes averaging over
174 a canonical ensemble of the initial state A and $\beta=1/k_B T$.
177 \section{The pull code}
178 \index{center-of-mass pulling}
180 The pull code applies forces or constraints between the centers
181 of mass of one or more pairs of groups of atoms.
182 Each pull reaction coordinate is called a ``coordinate'' and it operates
183 on two pull groups. A pull group can be part of one or more pull
184 coordinates. Furthermore, a coordinate can also operate on a single group
185 and an absolute reference position in space.
186 The distance between a pair of groups can be determined
187 in 1, 2 or 3 dimensions, or can be along a user-defined vector.
188 The reference distance can be constant or can change linearly with time.
189 Normally all atoms are weighted by their mass, but an additional
190 weighting factor can also be used.
192 \centerline{\includegraphics[width=6cm,angle=270]{plots/pull}}
193 \caption{Schematic picture of pulling a lipid out of a lipid bilayer
194 with umbrella pulling. $V_{rup}$ is the velocity at which the spring is
195 retracted, $Z_{link}$ is the atom to which the spring is attached and
196 $Z_{spring}$ is the location of the spring.}
200 Three different types of calculation are supported,
201 and in all cases the reference distance can be constant
202 or linearly changing with time.
204 \item{\textbf{Umbrella pulling}\swapindexquiet{umbrella}{pulling}}
205 A harmonic potential is applied between
206 the centers of mass of two groups.
207 Thus, the force is proportional to the displacement.
208 \item{\textbf{Constraint pulling\swapindexquiet{constraint}{pulling}}}
209 The distance between the centers of mass of two groups is constrained.
210 The constraint force can be written to a file.
211 This method uses the SHAKE algorithm but only needs 1 iteration to be
212 exact if only two groups are constrained.
213 \item{\textbf{Constant force pulling}}
214 A constant force is applied between the centers of mass of two groups.
215 Thus, the potential is linear.
216 In this case there is no reference distance of pull rate.
217 \item{\textbf{Flat bottom pulling}}
218 Like umbrella pulling, but the potential and force are zero for
219 negative deviations. This is useful for restraining the distance
220 between e.g. two molecules to a certain maximum distance.
223 \subsubsection{Definition of the center of mass}
225 In {\gromacs}, there are three ways to define the center of mass of a group.
226 The standard way is a ``plain'' center of mass, possibly with additional
227 weighting factors. With periodic boundary conditions it is no longer
228 possible to uniquely define the center of mass of a group of atoms.
229 Therefore, a reference atom is used. For determining the center of mass,
230 for all other atoms in the group, the closest periodic image to the reference
231 atom is used. This uniquely defines the center of mass.
232 By default, the middle (determined by the order in the topology) atom
233 is used as a reference atom, but the user can also select any other atom
234 if it would be closer to center of the group.
236 For a layered system, for instance a lipid bilayer, it may be of interest
237 to calculate the PMF of a lipid as function of its distance
238 from the whole bilayer. The whole bilayer can be taken as reference
239 group in that case, but it might also be of interest to define the
240 reaction coordinate for the PMF more locally. The {\tt .mdp} option
241 {\tt pull-coord?-geometry = cylinder} does not
242 use all the atoms of the reference group, but instead dynamically only those
243 within a cylinder with radius {\tt pull-cylinder-r} around the pull vector going
244 through the pull group. This only
245 works for distances defined in one dimension, and the cylinder is
246 oriented with its long axis along this one dimension. To avoid jumps in
247 the pull force, contributions of atoms are weighted as a function of distance
248 (in addition to the mass weighting):
250 w(r < r_\mathrm{cyl}) & = &
251 1-2 \left(\frac{r}{r_\mathrm{cyl}}\right)^2 + \left(\frac{r}{r_\mathrm{cyl}}\right)^4 \\
252 w(r \geq r_\mathrm{cyl}) & = & 0
254 Note that the radial dependence on the weight causes a radial force on
255 both cylinder group and the other pull group. This is an undesirable,
256 but unavoidable effect. To minimize this effect, the cylinder radius should
257 be chosen sufficiently large. The effective mass is 0.47 times that of
258 a cylinder with uniform weights and equal to the mass of uniform cylinder
259 of 0.79 times the radius.
262 \centerline{\includegraphics[width=6cm]{plots/pullref}}
263 \caption{Comparison of a plain center of mass reference group versus a cylinder
264 reference group applied to interface systems. C is the reference group.
265 The circles represent the center of mass of two groups plus the reference group,
266 $d_c$ is the reference distance.}
270 For a group of molecules in a periodic system, a plain reference group
271 might not be well-defined. An example is a water slab that is connected
272 periodically in $x$ and $y$, but has two liquid-vapor interfaces along $z$.
273 In such a setup, water molecules can evaporate from the liquid and they
274 will move through the vapor, through the periodic boundary, to the other
275 interface. Such a system is inherently periodic and there is no proper way
276 of defining a ``plain'' center of mass along $z$. A proper solution is to using
277 a cosine shaped weighting profile for all atoms in the reference group.
278 The profile is a cosine with a single period in the unit cell. Its phase
279 is optimized to give the maximum sum of weights, including mass weighting.
280 This provides a unique and continuous reference position that is nearly
281 identical to the plain center of mass position in case all atoms are all
282 within a half of the unit-cell length. See ref \cite{Engin2010a} for details.
284 When relative weights $w_i$ are used during the calculations, either
285 by supplying weights in the input or due to cylinder geometry
286 or due to cosine weighting,
287 the weights need to be scaled to conserve momentum:
290 \left. \sum_{j=1}^N w_j \, m_j \right/ \sum_{j=1}^N w_j^2 \, m_j
292 where $m_j$ is the mass of atom $j$ of the group.
293 The mass of the group, required for calculating the constraint force, is:
295 M = \sum_{i=1}^N w'_i \, m_i
297 The definition of the weighted center of mass is:
299 \ve{r}_{com} = \left. \sum_{i=1}^N w'_i \, m_i \, \ve{r}_i \right/ M
301 From the centers of mass the AFM, constraint, or umbrella force $\ve{F}_{\!com}$
302 on each group can be calculated.
303 The force on the center of mass of a group is redistributed to the atoms
306 \ve{F}_{\!i} = \frac{w'_i \, m_i}{M} \, \ve{F}_{\!com}
309 \subsubsection{Limitations}
310 There is one theoretical limitation:
311 strictly speaking, constraint forces can only be calculated between
312 groups that are not connected by constraints to the rest of the system.
313 If a group contains part of a molecule of which the bond lengths
314 are constrained, the pull constraint and LINCS or SHAKE bond constraint
315 algorithms should be iterated simultaneously. This is not done in {\gromacs}.
316 This means that for simulations with {\tt constraints = all-bonds}
317 in the {\tt .mdp} file pulling is, strictly speaking,
318 limited to whole molecules or groups of molecules.
319 In some cases this limitation can be avoided by using the free energy code,
321 In practice, the errors caused by not iterating the two constraint
322 algorithms can be negligible when the pull group consists of a large
323 amount of atoms and/or the pull force is small.
324 In such cases, the constraint correction displacement of the pull group
325 is small compared to the bond lengths.
329 \section{\normindex{Enforced Rotation}}
330 \index{rotational pulling|see{enforced rotation}}
331 \index{pulling, rotational|see{enforced rotation}}
334 \mathchardef\mhyphen="2D
335 \newcommand{\rotiso }{^\mathrm{iso}}
336 \newcommand{\rotisopf }{^\mathrm{iso\mhyphen pf}}
337 \newcommand{\rotpm }{^\mathrm{pm}}
338 \newcommand{\rotpmpf }{^\mathrm{pm\mhyphen pf}}
339 \newcommand{\rotrm }{^\mathrm{rm}}
340 \newcommand{\rotrmpf }{^\mathrm{rm\mhyphen pf}}
341 \newcommand{\rotrmtwo }{^\mathrm{rm2}}
342 \newcommand{\rotrmtwopf }{^\mathrm{rm2\mhyphen pf}}
343 \newcommand{\rotflex }{^\mathrm{flex}}
344 \newcommand{\rotflext }{^\mathrm{flex\mhyphen t}}
345 \newcommand{\rotflextwo }{^\mathrm{flex2}}
346 \newcommand{\rotflextwot}{^\mathrm{flex2\mhyphen t}}
348 This module can be used to enforce the rotation of a group of atoms, as {\eg}
349 a protein subunit. There are a variety of rotation potentials, among them
350 complex ones that allow flexible adaptations of both the rotated subunit as
351 well as the local rotation axis during the simulation. An example application
352 can be found in ref. \cite{Kutzner2011}.
355 \centerline{\includegraphics[width=13cm]{plots/rotation.pdf}}
356 \caption[Fixed and flexible axis rotation]{Comparison of fixed and flexible axis
358 {\sf A:} Rotating the sketched shape inside the white tubular cavity can create
359 artifacts when a fixed rotation axis (dashed) is used. More realistically, the
360 shape would revolve like a flexible pipe-cleaner (dotted) inside the bearing (gray).
361 {\sf B:} Fixed rotation around an axis \ve{v} with a pivot point
362 specified by the vector \ve{u}.
363 {\sf C:} Subdividing the rotating fragment into slabs with separate rotation
364 axes ($\uparrow$) and pivot points ($\bullet$) for each slab allows for
365 flexibility. The distance between two slabs with indices $n$ and $n+1$ is $\Delta x$.}
370 \centerline{\includegraphics[width=13cm]{plots/equipotential.pdf}}
371 \caption{Selection of different rotation potentials and definition of notation.
372 All four potentials $V$ (color coded) are shown for a single atom at position
374 {\sf A:} Isotropic potential $V\rotiso$,
375 {\sf B:} radial motion potential $V\rotrm$ and flexible potential
377 {\sf C--D:} radial motion\,2 potential $V\rotrmtwo$ and
378 flexible\,2 potential $V\rotflextwo$ for $\epsilon' = 0$\,nm$^2$ {\sf (C)}
379 and $\epsilon' = 0.01$\,nm$^2$ {\sf (D)}. The rotation axis is perpendicular to
380 the plane and marked by $\otimes$. The light gray contours indicate Boltzmann factors
381 $e^{-V/(k_B T)}$ in the $\ve{x}_j$-plane for $T=300$\,K and
382 $k=200$\,kJ/(mol$\cdot$nm$^2$). The green arrow shows the direction of the
383 force $\ve{F}_{\!j}$ acting on atom $j$; the blue dashed line indicates the
384 motion of the reference position.}
385 \label{fig:equipotential}
388 \subsection{Fixed Axis Rotation}
389 \subsubsection{Stationary Axis with an Isotropic Potential}
390 In the fixed axis approach (see \figref{rotation}B), torque on a group of $N$
391 atoms with positions $\ve{x}_i$ (denoted ``rotation group'') is applied by
392 rotating a reference set of atomic positions -- usually their initial positions
393 $\ve{y}_i^0$ -- at a constant angular velocity $\omega$ around an axis
394 defined by a direction vector $\hat{\ve{v}}$ and a pivot point \ve{u}.
395 To that aim, each atom with position $\ve{x}_i$ is attracted by a
396 ``virtual spring'' potential to its moving reference position
397 $\ve{y}_i = \mathbf{\Omega}(t) (\ve{y}_i^0 - \ve{u})$,
398 where $\mathbf{\Omega}(t)$ is a matrix that describes the rotation around the
399 axis. In the simplest case, the ``springs'' are described by a harmonic
402 V\rotiso = \frac{k}{2} \sum_{i=1}^{N} w_i \left[ \mathbf{\Omega}(t)
403 (\ve{y}_i^0 - \ve{u}) - (\ve{x}_i - \ve{u}) \right]^2 ,
406 with optional mass-weighted prefactors $w_i = N \, m_i/M$ with total mass
407 $M = \sum_{i=1}^N m_i$.
408 The rotation matrix $\mathbf{\Omega}(t)$ is
409 \newcommand{\omcost}{\,\xi\,} % abbreviation
414 \cos\omega t + v_x^2\omcost & v_x v_y\omcost - v_z\sin\omega t & v_x v_z\omcost + v_y\sin\omega t\\
415 v_x v_y\omcost + v_z\sin\omega t & \cos\omega t + v_y^2\omcost & v_y v_z\omcost - v_x\sin\omega t\\
416 v_x v_z\omcost - v_y\sin\omega t & v_y v_z\omcost + v_x\sin\omega t & \cos\omega t + v_z^2\omcost \\
420 where $v_x$, $v_y$, and $v_z$ are the components of the normalized rotation vector
421 $\hat{\ve{v}}$, and $\omcost := 1-\cos(\omega t)$. As illustrated in
422 \figref{equipotential}A for a single atom $j$, the
423 rotation matrix $\mathbf{\Omega}(t)$ operates on the initial reference positions
424 $\ve{y}_j^0 = \ve{x}_j(t_0)$ of atom $j$ at $t=t_0$. At a later
425 time $t$, the reference position has rotated away from its initial place
426 (along the blue dashed line), resulting in the force
429 = -\nabla_{\!j} \, V\rotiso
431 \mathbf{\Omega}(t) (\ve{y}_j^0 - \ve{u}) - (\ve{x}_j - \ve{u} ) \right] ,
432 \label{eqn:force_fixed}
434 which is directed towards the reference position.
437 \subsubsection{Pivot-Free Isotropic Potential}
438 Instead of a fixed pivot vector \ve{u} this potential uses the center of
439 mass $\ve{x}_c$ of the rotation group as pivot for the rotation axis,
441 \ve{x}_c = \frac{1}{M} \sum_{i=1}^N m_i \ve{x}_i
443 \mbox{\hspace{4ex}and\hspace{4ex}}
444 \ve{y}_c^0 = \frac{1}{M} \sum_{i=1}^N m_i \ve{y}_i^0 \ ,
446 which yields the ``pivot-free'' isotropic potential
448 V\rotisopf = \frac{k}{2} \sum_{i=1}^{N} w_i \left[ \mathbf{\Omega}(t)
449 (\ve{y}_i^0 - \ve{y}_c^0) - (\ve{x}_i - \ve{x}_c) \right]^2 ,
454 \mathbf{F}_{\!j}\rotisopf = k \, w_j
456 \mathbf{\Omega}(t) ( \ve{y}_j^0 - \ve{y}_c^0)
457 - ( \ve{x}_j - \ve{x}_c )
459 \label{eqn:force_isopf}
461 Without mass-weighting, the pivot $\ve{x}_c$ is the geometrical center of
465 \subsubsection{Parallel Motion Potential Variant}
466 The forces generated by the isotropic potentials
467 (\eqnsref{potiso}{potisopf}) also contain components parallel
468 to the rotation axis and thereby restrain motions along the axis of either the
469 whole rotation group (in case of $V\rotiso$) or within
470 the rotation group (in case of $V\rotisopf$). For cases where
471 unrestrained motion along the axis is preferred, we have implemented a
472 ``parallel motion'' variant by eliminating all components parallel to the
473 rotation axis for the potential. This is achieved by projecting the distance
474 vectors between reference and actual positions
476 \ve{r}_i = \mathbf{\Omega}(t) (\ve{y}_i^0 - \ve{u}) - (\ve{x}_i - \ve{u})
478 onto the plane perpendicular to the rotation vector,
481 \ve{r}_i^\perp := \ve{r}_i - (\ve{r}_i \cdot \hat{\ve{v}})\hat{\ve{v}} \ ,
486 V\rotpm &=& \frac{k}{2} \sum_{i=1}^{N} w_i ( \ve{r}_i^\perp )^2 \\
487 &=& \frac{k}{2} \sum_{i=1}^{N} w_i
490 (\ve{y}_i^0 - \ve{u}) - (\ve{x}_i - \ve{u}) \right. \nonumber \\
491 && \left. - \left\lbrace
492 \left[ \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{u}) - (\ve{x}_i - \ve{u}) \right] \cdot\hat{\ve{v}}
493 \right\rbrace\hat{\ve{v}} \right\rbrace^2 ,
498 \ve{F}_{\!j}\rotpm = k \, w_j \, \ve{r}_j^\perp .
502 \subsubsection{Pivot-Free Parallel Motion Potential}
503 Replacing in \eqnref{potpm} the fixed pivot \ve{u} by the center
504 of mass $\ve{x_c}$ yields the pivot-free variant of the parallel motion
507 \ve{s}_i = \mathbf{\Omega}(t) (\ve{y}_i^0 - \ve{y}_c^0) - (\ve{x}_i - \ve{x}_c)
509 the respective potential and forces are
511 V\rotpmpf &=& \frac{k}{2} \sum_{i=1}^{N} w_i ( \ve{s}_i^\perp )^2 \ , \\
513 \ve{F}_{\!j}\rotpmpf &=& k \, w_j \, \ve{s}_j^\perp .
514 \label{eqn:force_pmpf}
517 \subsubsection{Radial Motion Potential}
518 In the above variants, the minimum of the rotation potential is either a single
519 point at the reference position $\ve{y}_i$ (for the isotropic potentials) or a
520 single line through $\ve{y}_i$ parallel to the rotation axis (for the
521 parallel motion potentials). As a result, radial forces restrict radial motions
522 of the atoms. The two subsequent types of rotation potentials, $V\rotrm$
523 and $V\rotrmtwo$, drastically reduce or even eliminate this effect. The first
524 variant, $V\rotrm$ (\figref{equipotential}B), eliminates all force
525 components parallel to the vector connecting the reference atom and the
528 V\rotrm = \frac{k}{2} \sum_{i=1}^{N} w_i \left[
530 \cdot(\ve{x}_i - \ve{u}) \right]^2 ,
536 \frac{\hat{\ve{v}}\times \mathbf{\Omega}(t) (\ve{y}_i^0 - \ve{u})} {\| \hat{\ve{v}}\times \mathbf{\Omega}(t) (\ve{y}_i^0 - \ve{u})\|} \ .
538 This variant depends only on the distance $\ve{p}_i \cdot (\ve{x}_i -
539 \ve{u})$ of atom $i$ from the plane spanned by $\hat{\ve{v}}$ and
540 $\mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{u})$. The resulting force is
542 \mathbf{F}_{\!j}\rotrm =
543 -k \, w_j \left[ \ve{p}_j\cdot(\ve{x}_j - \ve{u}) \right] \,\ve{p}_j \, .
544 \label{eqn:potrm_force}
547 \subsubsection{Pivot-Free Radial Motion Potential}
548 Proceeding similar to the pivot-free isotropic potential yields a pivot-free
549 version of the above potential. With
552 \frac{\hat{\ve{v}}\times \mathbf{\Omega}(t) (\ve{y}_i^0 - \ve{y}_c^0)} {\| \hat{\ve{v}}\times \mathbf{\Omega}(t) (\ve{y}_i^0 - \ve{y}_c^0)\|} \, ,
554 the potential and force for the pivot-free variant of the radial motion potential read
556 V\rotrmpf & = & \frac{k}{2} \sum_{i=1}^{N} w_i \left[
558 \cdot(\ve{x}_i - \ve{x}_c)
561 \mathbf{F}_{\!j}\rotrmpf & = &
562 -k \, w_j \left[ \ve{q}_j\cdot(\ve{x}_j - \ve{x}_c) \right] \,\ve{q}_j
563 + k \frac{m_j}{M} \sum_{i=1}^{N} w_i \left[
564 \ve{q}_i\cdot(\ve{x}_i - \ve{x}_c) \right]\,\ve{q}_i \, .
565 \label{eqn:potrmpf_force}
568 \subsubsection{Radial Motion 2 Alternative Potential}
569 As seen in \figref{equipotential}B, the force resulting from
570 $V\rotrm$ still contains a small, second-order radial component. In most
571 cases, this perturbation is tolerable; if not, the following
572 alternative, $V\rotrmtwo$, fully eliminates the radial contribution to the
573 force, as depicted in \figref{equipotential}C,
576 \frac{k}{2} \sum_{i=1}^{N} w_i\,
577 \frac{\left[ (\hat{\ve{v}} \times ( \ve{x}_i - \ve{u} ))
578 \cdot \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{u}) \right]^2}
579 {\| \hat{\ve{v}} \times (\ve{x}_i - \ve{u}) \|^2 +
583 where a small parameter $\epsilon'$ has been introduced to avoid singularities.
584 For $\epsilon'=0$\,nm$^2$, the equipotential planes are spanned by $\ve{x}_i -
585 \ve{u}$ and $\hat{\ve{v}}$, yielding a force
586 perpendicular to $\ve{x}_i - \ve{u}$, thus not contracting or
587 expanding structural parts that moved away from or toward the rotation axis.
589 Choosing a small positive $\epsilon'$ ({\eg},
590 $\epsilon'=0.01$\,nm$^2$, \figref{equipotential}D) in the denominator of
591 \eqnref{potrm2} yields a well-defined potential and continuous forces also
592 close to the rotation axis, which is not the case for $\epsilon'=0$\,nm$^2$
593 (\figref{equipotential}C). With
595 \ve{r}_i & := & \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{u})\\
596 \ve{s}_i & := & \frac{\hat{\ve{v}} \times (\ve{x}_i -
597 \ve{u} ) }{ \| \hat{\ve{v}} \times (\ve{x}_i - \ve{u})
598 \| } \equiv \; \Psi_{i} \;\; {\hat{\ve{v}} \times
599 (\ve{x}_i-\ve{u} ) }\\
600 \Psi_i^{*} & := & \frac{1}{ \| \hat{\ve{v}} \times
601 (\ve{x}_i-\ve{u}) \|^2 + \epsilon'}
603 the force on atom $j$ reads
605 \ve{F}_{\!j}\rotrmtwo =
608 (\ve{s}_j\cdot\ve{r}_{\!j})\;
609 \left[ \frac{\Psi_{\!j}^* }{\Psi_{\!j} } \ve{r}_{\!j}
610 - \frac{\Psi_{\!j}^{*2}}{\Psi_{\!j}^3}
611 (\ve{s}_j\cdot\ve{r}_{\!j})\ve{s}_j \right]
612 \right\rbrace \times \hat{\ve{v}} .
613 \label{eqn:potrm2_force}
616 \subsubsection{Pivot-Free Radial Motion 2 Potential}
617 The pivot-free variant of the above potential is
620 \frac{k}{2} \sum_{i=1}^{N} w_i\,
621 \frac{\left[ (\hat{\ve{v}} \times ( \ve{x}_i - \ve{x}_c ))
622 \cdot \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{y}_c) \right]^2}
623 {\| \hat{\ve{v}} \times (\ve{x}_i - \ve{x}_c) \|^2 +
629 \ve{r}_i & := & \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{y}_c)\\
630 \ve{s}_i & := & \frac{\hat{\ve{v}} \times (\ve{x}_i -
631 \ve{x}_c ) }{ \| \hat{\ve{v}} \times (\ve{x}_i - \ve{x}_c)
632 \| } \equiv \; \Psi_{i} \;\; {\hat{\ve{v}} \times
633 (\ve{x}_i-\ve{x}_c ) }\\ \Psi_i^{*} & := & \frac{1}{ \| \hat{\ve{v}} \times
634 (\ve{x}_i-\ve{x}_c) \|^2 + \epsilon'}
636 the force on atom $j$ reads
639 \ve{F}_{\!j}\rotrmtwopf & = &
642 (\ve{s}_j\cdot\ve{r}_{\!j})\;
643 \left[ \frac{\Psi_{\!j}^* }{\Psi_{\!j} } \ve{r}_{\!j}
644 - \frac{\Psi_{\!j}^{*2}}{\Psi_{\!j}^3}
645 (\ve{s}_j\cdot\ve{r}_{\!j})\ve{s}_j \right]
646 \right\rbrace \times \hat{\ve{v}}\\
648 + k\;\frac{m_j}{M} \left\lbrace \sum_{i=1}^{N}
649 w_i\;(\ve{s}_i\cdot\ve{r}_i) \;
650 \left[ \frac{\Psi_i^* }{\Psi_i } \ve{r}_i
651 - \frac{\Psi_i^{*2}}{\Psi_i^3} (\ve{s}_i\cdot\ve{r}_i )\;
652 \ve{s}_i \right] \right\rbrace \times \hat{\ve{v}} \, .
653 \label{eqn:potrm2pf_force}
656 \subsection{Flexible Axis Rotation}
657 As sketched in \figref{rotation}A--B, the rigid body behavior of
658 the fixed axis rotation scheme is a drawback for many applications. In
659 particular, deformations of the rotation group are suppressed when the
660 equilibrium atom positions directly depend on the reference positions.
661 To avoid this limitation, \eqnsref{potrmpf}{potrm2pf}
662 will now be generalized towards a ``flexible axis'' as sketched in
663 \figref{rotation}C. This will be achieved by subdividing the
664 rotation group into a set of equidistant slabs perpendicular to
665 the rotation vector, and by applying a separate rotation potential to each
666 of these slabs. \figref{rotation}C shows the midplanes of the slabs
667 as dotted straight lines and the centers as thick black dots.
669 To avoid discontinuities in the potential and in the forces, we define
670 ``soft slabs'' by weighing the contributions of each
671 slab $n$ to the total potential function $V\rotflex$ by a Gaussian
675 g_n(\ve{x}_i) = \Gamma \ \mbox{exp} \left(
676 -\frac{\beta_n^2(\ve{x}_i)}{2\sigma^2} \right) ,
678 centered at the midplane of the $n$th slab. Here $\sigma$ is the width
679 of the Gaussian function, $\Delta x$ the distance between adjacent slabs, and
681 \beta_n(\ve{x}_i) := \ve{x}_i \cdot \hat{\ve{v}} - n \, \Delta x \, .
685 \centerline{\includegraphics[width=6.5cm]{plots/gaussians.pdf}}
686 \caption{Gaussian functions $g_n$ centered at $n \, \Delta x$ for a slab
687 distance $\Delta x = 1.5$ nm and $n \geq -2$. Gaussian function $g_0$ is
688 highlighted in bold; the dashed line depicts the sum of the shown Gaussian
690 \label{fig:gaussians}
693 A most convenient choice is $\sigma = 0.7 \Delta x$ and
695 1/\Gamma = \sum_{n \in Z}
697 \left(-\frac{(n - \frac{1}{4})^2}{2\cdot 0.7^2}\right)
700 which yields a nearly constant sum, essentially independent of $\ve{x}_i$
701 (dashed line in \figref{gaussians}), {\ie},
703 \sum_{n \in Z} g_n(\ve{x}_i) = 1 + \epsilon(\ve{x}_i) \, ,
706 with $ | \epsilon(\ve{x}_i) | < 1.3\cdot 10^{-4}$. This choice also
707 implies that the individual contributions to the force from the slabs add up to
708 unity such that no further normalization is required.
710 To each slab center $\ve{x}_c^n$, all atoms contribute by their
711 Gaussian-weighted (optionally also mass-weighted) position vectors
712 $g_n(\ve{x}_i) \, \ve{x}_i$. The
713 instantaneous slab centers $\ve{x}_c^n$ are calculated from the
714 current positions $\ve{x}_i$,
718 \frac{\sum_{i=1}^N g_n(\ve{x}_i) \, m_i \, \ve{x}_i}
719 {\sum_{i=1}^N g_n(\ve{x}_i) \, m_i} \, ,\\
721 while the reference centers $\ve{y}_c^n$ are calculated from the reference
722 positions $\ve{y}_i^0$,
726 \frac{\sum_{i=1}^N g_n(\ve{y}_i^0) \, m_i \, \ve{y}_i^0}
727 {\sum_{i=1}^N g_n(\ve{y}_i^0) \, m_i} \, .
729 Due to the rapid decay of $g_n$, each slab
730 will essentially involve contributions from atoms located within $\approx
731 3\Delta x$ from the slab center only.
733 \subsubsection{Flexible Axis Potential}
734 We consider two flexible axis variants. For the first variant,
735 the slab segmentation procedure with Gaussian weighting is applied to the radial
736 motion potential (\eqnref{potrmpf}\,/\,\figref{equipotential}B),
737 yielding as the contribution of slab $n$
740 \frac{k}{2} \sum_{i=1}^{N} w_i \, g_n(\ve{x}_i)
744 (\ve{x}_i - \ve{x}_c^n)
748 and a total potential function
750 V\rotflex = \sum_n V^n \, .
753 Note that the global center of mass $\ve{x}_c$ used in
754 \eqnref{potrmpf} is now replaced by $\ve{x}_c^n$, the center of mass of
757 \ve{q}_i^n & := & \frac{\hat{\ve{v}} \times
758 \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{y}_c^n) }{ \| \hat{\ve{v}}
759 \times \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{y}_c^n) \| } \\
760 b_i^n & := & \ve{q}_i^n \cdot (\ve{x}_i - \ve{x}_c^n) \, ,
762 the resulting force on atom $j$ reads
764 \nonumber\hspace{-15mm}
765 \ve{F}_{\!j}\rotflex &=&
766 - \, k \, w_j \sum_n g_n(\ve{x}_j) \, b_j^n \left\lbrace \ve{q}_j^n -
767 b_j^n \frac{\beta_n(\ve{x}_j)}{2\sigma^2} \hat{\ve{v}} \right\rbrace \\ & &
768 + \, k \, m_j \sum_n \frac{g_n(\ve{x}_j)}{\sum_h g_n(\ve{x}_h)}
769 \sum_{i=1}^{N} w_i \, g_n(\ve{x}_i) \, b_i^n \left\lbrace
770 \ve{q}_i^n -\frac{\beta_n(\ve{x}_j)}{\sigma^2}
771 \left[ \ve{q}_i^n \cdot (\ve{x}_j - \ve{x}_c^n )\right]
772 \hat{\ve{v}} \right\rbrace .
773 \label{eqn:potflex_force}
776 Note that for $V\rotflex$, as defined, the slabs are fixed in space and so
777 are the reference centers $\ve{y}_c^n$. If during the simulation the
778 rotation group moves too far in $\ve{v}$ direction, it may enter a
779 region where -- due to the lack of nearby reference positions -- no reference
780 slab centers are defined, rendering the potential evaluation impossible.
781 We therefore have included a slightly modified version of this potential that
782 avoids this problem by attaching the midplane of slab $n=0$ to the center of mass
783 of the rotation group, yielding slabs that move with the rotation group.
784 This is achieved by subtracting the center of mass $\ve{x}_c$ of the
785 group from the positions,
787 \tilde{\ve{x}}_i = \ve{x}_i - \ve{x}_c \, , \mbox{\ \ \ and \ \ }
788 \tilde{\ve{y}}_i^0 = \ve{y}_i^0 - \ve{y}_c^0 \, ,
794 & = & \frac{k}{2} \sum_n \sum_{i=1}^{N} w_i \, g_n(\tilde{\ve{x}}_i)
795 \left[ \frac{\hat{\ve{v}} \times \mathbf{\Omega}(t)(\tilde{\ve{y}}_i^0
796 - \tilde{\ve{y}}_c^n) }{ \| \hat{\ve{v}} \times
797 \mathbf{\Omega}(t)(\tilde{\ve{y}}_i^0 -
798 \tilde{\ve{y}}_c^n) \| }
800 (\tilde{\ve{x}}_i - \tilde{\ve{x}}_c^n)
804 To simplify the force derivation, and for efficiency reasons, we here assume
805 $\ve{x}_c$ to be constant, and thus $\partial \ve{x}_c / \partial x =
806 \partial \ve{x}_c / \partial y = \partial \ve{x}_c / \partial z = 0$. The
807 resulting force error is small (of order $O(1/N)$ or $O(m_j/M)$ if
808 mass-weighting is applied) and can therefore be tolerated. With this assumption,
809 the forces $\ve{F}\rotflext$ have the same form as
810 \eqnref{potflex_force}.
812 \subsubsection{Flexible Axis 2 Alternative Potential}
813 In this second variant, slab segmentation is applied to $V\rotrmtwo$
814 (\eqnref{potrm2pf}), resulting in a flexible axis potential without radial
815 force contributions (\figref{equipotential}C),
818 \frac{k}{2} \sum_{i=1}^{N} \sum_n w_i\,g_n(\ve{x}_i)
819 \frac{\left[ (\hat{\ve{v}} \times ( \ve{x}_i - \ve{x}_c^n ))
820 \cdot \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{y}_c^n) \right]^2}
821 {\| \hat{\ve{v}} \times (\ve{x}_i - \ve{x}_c^n) \|^2 +
827 \ve{r}_i^n & := & \mathbf{\Omega}(t)(\ve{y}_i^0 - \ve{y}_c^n)\\
828 \ve{s}_i^n & := & \frac{\hat{\ve{v}} \times (\ve{x}_i -
829 \ve{x}_c^n ) }{ \| \hat{\ve{v}} \times (\ve{x}_i - \ve{x}_c^n)
830 \| } \equiv \; \psi_{i} \;\; {\hat{\ve{v}} \times (\ve{x}_i-\ve{x}_c^n ) }\\
831 \psi_i^{*} & := & \frac{1}{ \| \hat{\ve{v}} \times (\ve{x}_i-\ve{x}_c^n) \|^2 + \epsilon'}\\
832 W_j^n & := & \frac{g_n(\ve{x}_j)\,m_j}{\sum_h g_n(\ve{x}_h)\,m_h}\\
834 \sum_{i=1}^{N} w_i\;g_n(\ve{x}_i)
835 \; (\ve{s}_i^n\cdot\ve{r}_i^n)
836 \left[ \frac{\psi_i^* }{\psi_i } \ve{r}_i^n
837 - \frac{\psi_i^{*2}}{\psi_i^3} (\ve{s}_i^n\cdot\ve{r}_i^n )\;
838 \ve{s}_i^n \right] \label{eqn:Sn}
840 the force on atom $j$ reads
843 \ve{F}_{\!j}\rotflextwo & = &
845 \left\lbrace \sum_n w_j\;g_n(\ve{x}_j)\;
846 (\ve{s}_j^n\cdot\ve{r}_{\!j}^n)\;
847 \left[ \frac{\psi_j^* }{\psi_j } \ve{r}_{\!j}^n
848 - \frac{\psi_j^{*2}}{\psi_j^3} (\ve{s}_j^n\cdot\ve{r}_{\!j}^n)\;
849 \ve{s}_{\!j}^n \right] \right\rbrace \times \hat{\ve{v}} \\
852 + k \left\lbrace \sum_n W_{\!j}^n \, \ve{S}^n \right\rbrace \times
854 - k \left\lbrace \sum_n W_{\!j}^n \; \frac{\beta_n(\ve{x}_j)}{\sigma^2} \frac{1}{\psi_j}\;\;
856 \ve{S}^n \right\rbrace \hat{\ve{v}}\\
858 + \frac{k}{2} \left\lbrace \sum_n w_j\;g_n(\ve{x}_j)
859 \frac{\beta_n(\ve{x}_j)}{\sigma^2}
860 \frac{\psi_j^*}{\psi_j^2}( \ve{s}_j^n \cdot \ve{r}_{\!j}^n )^2 \right\rbrace
862 \label{eqn:potflex2_force}
865 Applying transformation (\ref{eqn:trafo}) yields a ``translation-tolerant''
866 version of the flexible\,2 potential, $V\rotflextwot$. Again,
867 assuming that $\partial \ve{x}_c / \partial x$, $\partial \ve{x}_c /
868 \partial y$, $\partial \ve{x}_c / \partial z$ are small, the
869 resulting equations for $V\rotflextwot$ and $\ve{F}\rotflextwot$ are
870 similar to those of $V\rotflextwo$ and $\ve{F}\rotflextwo$.
873 To apply enforced rotation, the particles $i$ that are to
874 be subjected to one of the rotation potentials are defined via index groups
875 {\tt rot-group0}, {\tt rot-group1}, etc., in the {\tt .mdp} input file.
876 The reference positions $\ve{y}_i^0$ are
877 read from a special {\tt .trr} file provided to {\tt grompp}. If no such file is found,
878 $\ve{x}_i(t=0)$ are used as reference positions and written to {\tt .trr} such
879 that they can be used for subsequent setups. All parameters of the potentials
880 such as $k$, $\epsilon'$, etc. (\tabref{vars}) are provided as {\tt .mdp}
881 parameters; {\tt rot-type} selects the type of the potential.
882 The option {\tt rot-massw} allows to choose whether or not to use
883 mass-weighted averaging.
884 For the flexible potentials, a cutoff value $g_n^\mathrm{min}$
885 (typically $g_n^\mathrm{min}=0.001$) makes shure that only
886 significant contributions to $V$ and \ve{F} are evaluated, {\ie} terms with
887 $g_n(\ve{x}) < g_n^\mathrm{min}$ are omitted.
888 \tabref{quantities} summarizes observables that are written
889 to additional output files and which are described below.
893 \caption{Parameters used by the various rotation potentials.
894 {\sf x}'s indicate which parameter is actually used for a given potential.}
897 \newcommand{\kunit}{$\frac{\mathrm{kJ}}{\mathrm{mol} \cdot \mathrm{nm}^2}$}
898 \newcommand{\smtt}[1]{{\hspace{-0.5ex}\small #1\hspace{-0.5ex}}}
901 \begin{tabular}{l>{$}l<{$}rccccccc}
903 parameter & & & $k$ & $\hat{\ve{v}}$ & $\ve{u}$ & $\omega$ & $\epsilon'$ & $\Delta x$ & $g_n^\mathrm{min}$ \\
904 \multicolumn{3}{l}{{\tt .mdp} input variable name} & \smtt{k} & \smtt{vec} & \smtt{pivot} & \smtt{rate} & \smtt{eps} & \smtt{slab-dist} & \smtt{min-gauss} \\
905 unit & & & \kunit & - & nm & $^\circ$/ps & nm$^2$ & nm & \,-\, \\[1mm]
906 \hline \multicolumn{2}{l}{fixed axis potentials:} & eqn.\\
907 isotropic & V\rotiso & (\ref{eqn:potiso}) & {\sf x} & {\sf x} & {\sf x} & {\sf x} & - & - & - \\
908 --- pivot-free & V\rotisopf & (\ref{eqn:potisopf}) & {\sf x} & {\sf x} & - & {\sf x} & - & - & - \\
909 parallel motion & V\rotpm & (\ref{eqn:potpm}) & {\sf x} & {\sf x} & {\sf x} & {\sf x} & - & - & - \\
910 --- pivot-free & V\rotpmpf & (\ref{eqn:potpmpf}) & {\sf x} & {\sf x} & - & {\sf x} & - & - & - \\
911 radial motion & V\rotrm & (\ref{eqn:potrm}) & {\sf x} & {\sf x} & {\sf x} & {\sf x} & - & - & - \\
912 --- pivot-free & V\rotrmpf & (\ref{eqn:potrmpf}) & {\sf x} & {\sf x} & - & {\sf x} & - & - & - \\
913 radial motion\,2 & V\rotrmtwo & (\ref{eqn:potrm2}) & {\sf x} & {\sf x} & {\sf x} & {\sf x} & {\sf x} & - & - \\
914 --- pivot-free & V\rotrmtwopf & (\ref{eqn:potrm2pf}) & {\sf x} & {\sf x} & - & {\sf x} & {\sf x} & - & - \\ \hline
915 \multicolumn{2}{l}{flexible axis potentials:} & eqn.\\
916 flexible & V\rotflex & (\ref{eqn:potflex}) & {\sf x} & {\sf x} & - & {\sf x} & - & {\sf x} & {\sf x} \\
917 --- transl. tol. & V\rotflext & (\ref{eqn:potflext}) & {\sf x} & {\sf x} & - & {\sf x} & - & {\sf x} & {\sf x} \\
918 flexible\,2 & V\rotflextwo & (\ref{eqn:potflex2}) & {\sf x} & {\sf x} & - & {\sf x} & {\sf x} & {\sf x} & {\sf x} \\
919 --- transl. tol. & V\rotflextwot & - & {\sf x} & {\sf x} & - & {\sf x} & {\sf x} & {\sf x} & {\sf x} \\
926 \caption{Quantities recorded in output files during enforced rotation simulations.
927 All slab-wise data is written every {\tt nstsout} steps, other rotation data every {\tt nstrout} steps.}
928 \label{tab:quantities}
930 \begin{tabular}{llllcc}
932 quantity & unit & equation & output file & fixed & flexible\\ \hline
933 $V(t)$ & kJ/mol & see \ref{tab:vars} & {\tt rotation} & {\sf x} & {\sf x} \\
934 $\theta_\mathrm{ref}(t)$ & degrees & $\theta_\mathrm{ref}(t)=\omega t$ & {\tt rotation} & {\sf x} & {\sf x} \\
935 $\theta_\mathrm{av}(t)$ & degrees & (\ref{eqn:avangle}) & {\tt rotation} & {\sf x} & - \\
936 $\theta_\mathrm{fit}(t)$, $\theta_\mathrm{fit}(t,n)$ & degrees & (\ref{eqn:rmsdfit}) & {\tt rotangles} & - & {\sf x} \\
937 $\ve{y}_0(n)$, $\ve{x}_0(t,n)$ & nm & (\ref{eqn:defx0}, \ref{eqn:defy0})& {\tt rotslabs} & - & {\sf x} \\
938 $\tau(t)$ & kJ/mol & (\ref{eqn:torque}) & {\tt rotation} & {\sf x} & - \\
939 $\tau(t,n)$ & kJ/mol & (\ref{eqn:torque}) & {\tt rottorque} & - & {\sf x} \\ \hline
945 \subsubsection*{Angle of Rotation Groups: Fixed Axis}
946 For fixed axis rotation, the average angle $\theta_\mathrm{av}(t)$ of the
947 group relative to the reference group is determined via the distance-weighted
948 angular deviation of all rotation group atoms from their reference positions,
950 \theta_\mathrm{av} = \left. \sum_{i=1}^{N} r_i \ \theta_i \right/ \sum_{i=1}^N r_i \ .
953 Here, $r_i$ is the distance of the reference position to the rotation axis, and
954 the difference angles $\theta_i$ are determined from the atomic positions,
955 projected onto a plane perpendicular to the rotation axis through pivot point
956 $\ve{u}$ (see \eqnref{project} for the definition of $\perp$),
959 \frac{(\ve{y}_i-\ve{u})^\perp \cdot (\ve{x}_i-\ve{u})^\perp}
960 { \| (\ve{y}_i-\ve{u})^\perp \cdot (\ve{x}_i-\ve{u})^\perp
964 The sign of $\theta_\mathrm{av}$ is chosen such that
965 $\theta_\mathrm{av} > 0$ if the actual structure rotates ahead of the reference.
967 \subsubsection*{Angle of Rotation Groups: Flexible Axis}
968 For flexible axis rotation, two outputs are provided, the angle of the
969 entire rotation group, and separate angles for the segments in the slabs.
970 The angle of the entire rotation group is determined by an RMSD fit
972 to the reference positions $\ve{y}_i^0$ at $t=0$, yielding $\theta_\mathrm{fit}$
973 as the angle by which the reference has to be rotated around $\hat{\ve{v}}$
976 \mathrm{RMSD} \big( \ve{x}_i,\ \mathbf{\Omega}(\theta_\mathrm{fit})
977 \ve{y}_i^0 \big) \stackrel{!}{=} \mathrm{min} \, .
980 To determine the local angle for each slab $n$, both reference and actual
981 positions are weighted with the Gaussian function of slab $n$, and
982 $\theta_\mathrm{fit}(t,n)$ is calculated as in \eqnref{rmsdfit}) from the
983 Gaussian-weighted positions.
985 For all angles, the {\tt .mdp} input option {\tt rot-fit-method} controls
986 whether a normal RMSD fit is performed or whether for the fit each
987 position $\ve{x}_i$ is put at the same distance to the rotation axis as its
988 reference counterpart $\ve{y}_i^0$. In the latter case, the RMSD
989 measures only angular differences, not radial ones.
992 \subsubsection*{Angle Determination by Searching the Energy Minimum}
993 Alternatively, for {\tt rot-fit-method = potential}, the angle of the rotation
994 group is determined as the angle for which the rotation potential energy is minimal.
995 Therefore, the used rotation potential is additionally evaluated for a set of angles
996 around the current reference angle. In this case, the {\tt rotangles.log} output file
997 contains the values of the rotation potential at the chosen set of angles, while
998 {\tt rotation.xvg} lists the angle with minimal potential energy.
1001 \subsubsection*{Torque}
1003 The torque $\ve{\tau}(t)$ exerted by the rotation potential is calculated for fixed
1006 \ve{\tau}(t) = \sum_{i=1}^{N} \ve{r}_i(t) \times \ve{f}_{\!i}^\perp(t) ,
1009 where $\ve{r}_i(t)$ is the distance vector from the rotation axis to
1010 $\ve{x}_i(t)$ and $\ve{f}_{\!i}^\perp(t)$ is the force component
1011 perpendicular to $\ve{r}_i(t)$ and $\hat{\ve{v}}$. For flexible axis
1012 rotation, torques $\ve{\tau}_{\!n}$ are calculated for each slab using the
1013 local rotation axis of the slab and the Gaussian-weighted positions.
1016 \section{\normindex{Computational Electrophysiology}}
1019 The Computational Electrophysiology (CompEL) protocol \cite{Kutzner2011b} allows the simulation of
1020 ion flux through membrane channels, driven by transmembrane potentials or ion
1021 concentration gradients. Just as in real cells, CompEL establishes transmembrane
1022 potentials by sustaining a small imbalance of charges $\Delta q$ across the membrane,
1023 which gives rise to a potential difference $\Delta U$ according to the membrane capacitance:
1025 \Delta U = \Delta q / C_{membrane}
1027 The transmembrane electric field and concentration gradients are controlled by
1028 {\tt .mdp} options, which allow the user to set reference counts for the ions on either side
1029 of the membrane. If a difference between the actual and the reference numbers persists
1030 over a certain time span, specified by the user, a number of ion/water pairs are
1031 exchanged between the compartments until the reference numbers are restored.
1032 Alongside the calculation of channel conductance and ion selectivity, CompEL simulations also
1033 enable determination of the channel reversal potential, an important
1034 characteristic obtained in electrophysiology experiments.
1036 In a CompEL setup, the simulation system is divided into two compartments {\bf A} and {\bf B}
1037 with independent ion concentrations. This is best achieved by using double bilayer systems with
1038 a copy (or copies) of the channel/pore of interest in each bilayer (\figref{compelsetup} A, B).
1039 If the channel axes point in the same direction, channel flux is observed
1040 simultaneously at positive and negative potentials in this way, which is for instance
1041 important for studying channel rectification.
1044 \centerline{\includegraphics[width=13.5cm]{plots/compelsetup.pdf}}
1045 \caption{Typical double-membrane setup for CompEL simulations (A, B). Plot (C) shows
1046 the potential difference $\Delta U$ resulting
1047 from the selected charge imbalance $\Delta q_{ref}$ between the compartments.}
1048 \label{fig:compelsetup}
1051 The potential difference $\Delta U$ across the membrane is easily calculated with the
1052 {\tt gmx potential} utility. By this, the potential drop along $z$ or the
1053 pore axis is exactly known in each time interval of the simulation (\figref{compelsetup} C).
1054 Type and number of ions $n_i$ of charge $q_i$, traversing the channel in the simulation,
1055 are written to the {\tt swapions.xvg} output file, from which the average channel
1056 conductance $G$ in each interval $\Delta t$ is determined by:
1058 G = \frac{\sum_{i} n_{i}q_{i}}{\Delta t \, \Delta U} \, .
1060 The ion selectivity is calculated as the number flux ratio of different species.
1061 Best results are obtained by averaging these values over several overlapping time intervals.
1063 The calculation of reversal potentials is best achieved using a small set of simulations in which a given
1064 transmembrane concentration gradient is complemented with small ion imbalances of varying magnitude. For
1065 example, if one compartment contains 1\,M salt and the other 0.1\,M, and given charge neutrality otherwise,
1066 a set of simulations with $\Delta q = 0\,e$, $\Delta q = 2\,e$, $\Delta q = 4\,e$ could
1067 be used. Fitting a straight line through the current-voltage relationship of all obtained
1068 $I$-$U$ pairs near zero current will then yield $U_{rev}$.
1071 The following {\tt .mdp} options control the CompEL protocol:
1074 swapcoords = Z ; Swap positions: no, X, Y, Z
1075 swap-frequency = 100 ; Swap attempt frequency
1077 Choose {\tt Z} if your membrane is in the $xy$-plane (\figref{compelsetup} A, B).
1078 Ions will be exchanged between compartments depending on their $z$-positions alone.
1079 {\tt swap-frequency} determines how often a swap attempt will be made.
1080 This step requires that the positions of the ions, solvent, and swap groups are
1081 communicated between the parallel processes, so if chosen too small it can decrease the simulation
1085 split-group0 = channel0 ; Defines compartment boundary
1086 split-group1 = channel1 ; Defines other compartment boundary
1087 massw-split0 = no ; use mass-weighted center?
1090 {\tt split-group0} and {\tt split-group1} are two index groups that define the boundaries
1091 between the two compartments, which are usually the centers of the channels.
1092 If {\tt massw-split0} or {\tt massw-split1} are set to {\tt yes}, the center of mass
1093 of each index group is used as boundary, here in $z$-direction. Otherwise, the
1094 geometrical centers will be used ($\times$ in \figref{compelsetup} A). If, such as here, a membrane
1095 channel is selected as split group, the center of the channel will define the dividing
1096 plane between the compartments (dashed horizontal line in the figure). All index groups
1097 must be defined in the index file.
1100 swap-group = NA+_CL- ; Ions to be included in exchange
1101 solvent-group = SOL ; Group name of solvent molecules
1102 cyl0-r = 5.0 ; Split cylinder 0: pore radius (nm)
1103 cyl0-up = 0.75 ; Split cylinder 0 upper extension (nm)
1104 cyl0-down = 0.75 ; Split cylinder 0 lower extension (nm)
1105 cyl1-r = 5.0 ; same for other channel
1108 coupl-steps = 10 ; Average over these many swap steps
1109 threshold = 1 ; Do not swap if < threshold
1111 {\tt swap-group} identifies the index group of ions that
1112 should be involved in the flux and exchange cycles, {\tt solvent-group} defines the solvent
1113 group with which they are swapped. The cylinder options only influence the counting of
1114 ions, i.e., ions will be counted as having traveled through either channel 0 or channel 1
1115 according to the definition of (channel) cylinder radius, upper and lower extension,
1116 relative to the location of the respective split group. This will not affect the actual
1117 flux or exchange, but will provide you with the ion permeation numbers across
1118 each of the channels. Note that an ion can only be counted as passing through a particular
1119 channel if it is detected \emph{within} the defined split cylinder in a swap step.
1120 If {\tt swap-frequency} is chosen too high, a particular ion may be detected in compartment {\bf A}
1121 in one swap step, and in compartment {\bf B} in the following swap step, so it will be unclear
1122 through which of the channels it has passed.
1124 {\tt coupl-steps} sets the number of swap attempt steps. A discrepancy between
1125 actual and reference ion numbers in each compartment must persist over this many attempts
1126 before an actual exchange takes place. If {\tt coupl-steps} is set to 1, then the momentary ion distribution determines
1127 whether ions are exchanged. {\tt coupl-steps} \textgreater\ 1 will use the time-average
1128 of ion distributions over the selected number of attempt steps instead. This can be useful, for example,
1129 when ions diffuse near compartment boundaries, which would lead to numerous unproductive
1130 ion exchanges. A {\tt threshold} of 1 means that a swap is performed if the average ion
1131 count in a compartment differs by at least 1 from the requested values. Higher thresholds
1132 will lead to toleration of larger differences. Ions are exchanged until the requested
1133 number $\pm$ the threshold is reached.
1137 anionsA = -1 ; Reference count of anions in A
1138 cationsA = -1 ; ... of cations in A
1139 anionsB = -1 ; ... of anions in B
1140 cationsB = -1 ; ... of cations in B
1142 These options set the requested number of anions and cations for each of the two compartments.
1143 A number of {\tt -1} means fix the numbers found in time step 0. Note that these numbers
1144 need to add up to the total number of ions in the swap group.
1146 Note that a double-layered system for CompEL simulations can be easily prepared by
1147 duplicating an existing membrane/channel MD system in the direction of the membrane
1148 normal (typically $z$) with {\tt gmx editconf -translate 0 0 <l_z>}, where {\tt l_z}
1149 is the box length in that direction. If you have already defined index groups for
1150 the channel for the single-layered system, {\tt gmx make_ndx -n index.ndx -twin} will
1151 provide you with the groups for the double-layered system.
1153 To suppress large fluctuations of the membranes along the swap direction,
1154 it may be useful to apply a harmonic potential (acting only in the swap dimension)
1155 between each of the two channel and/or bilayer centers using umbrella pulling
1156 (see section~\ref{sec:pull}).
1158 \subsection*{Multimeric channels}
1159 If a split group consists of more than one molecule, the correct PBC image of all molecules
1160 with respect to each other has to be chosen such that the channel center can be correctly
1161 determined. \gromacs\ assumes that the starting structure in the {\tt .tpr}
1162 file has the correct PBC representation. Set the following environment variable
1163 to check whether that is the case:
1165 \item {\tt GMX_COMPELDUMP}: output the starting structure after it has been made whole to
1170 \section{Calculating a PMF using the free-energy code}
1172 \index{potentials of mean force}
1173 \index{free energy calculations}
1174 The free-energy coupling-parameter approach (see~\secref{fecalc})
1175 provides several ways to calculate potentials of mean force.
1176 A potential of mean force between two atoms can be calculated
1177 by connecting them with a harmonic potential or a constraint.
1178 For this purpose there are special potentials that avoid the generation of
1179 extra exclusions, see~\secref{excl}.
1180 When the position of the minimum or the constraint length is 1 nm more
1181 in state B than in state A, the restraint or constraint force is given
1182 by $\partial H/\partial \lambda$.
1183 The distance between the atoms can be changed as a function of $\lambda$
1184 and time by setting {\tt delta-lambda} in the {\tt .mdp} file.
1185 The results should be identical (although not numerically
1186 due to the different implementations) to the results of the pull code
1187 with umbrella sampling and constraint pulling.
1188 Unlike the pull code, the free energy code can also handle atoms that
1189 are connected by constraints.
1191 Potentials of mean force can also be calculated using position restraints.
1192 With position restraints, atoms can be linked to a position in space
1193 with a harmonic potential (see \ssecref{positionrestraint}).
1194 These positions can be made a function of the coupling parameter $\lambda$.
1195 The positions for the A and the B states are supplied to {\tt grompp} with
1196 the {\tt -r} and {\tt -rb} options, respectively.
1197 One could use this approach to do \normindex{targeted MD};
1198 note that we do not encourage the use of targeted MD for proteins.
1199 A protein can be forced from one conformation to another by using
1200 these conformations as position restraint coordinates for state A and B.
1201 One can then slowly change $\lambda$ from 0 to 1.
1202 The main drawback of this approach is that the conformational freedom
1203 of the protein is severely limited by the position restraints,
1204 independent of the change from state A to B.
1205 Also, the protein is forced from state A to B in an almost straight line,
1206 whereas the real pathway might be very different.
1207 An example of a more fruitful application is a solid system or a liquid
1208 confined between walls where one wants to measure the force required
1209 to change the separation between the boundaries or walls.
1210 Because the boundaries (or walls) already need to be fixed,
1211 the position restraints do not limit the system in its sampling.
1213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1214 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1216 \newcommand{\amine}{\sf -NH$_2$}
1217 \newcommand{\amines}{\sf -NH-}
1218 \newcommand{\aminep}{\sf -NH$_3^+$}
1219 \section{Removing fastest \swapindex{degrees of}{freedom}}
1221 The maximum time step in MD simulations is limited by the smallest
1222 oscillation period that can be found in the simulated
1223 system. Bond-stretching vibrations are in their quantum-mechanical
1224 ground state and are therefore better represented by a constraint
1225 instead of a harmonic potential.
1227 For the remaining degrees of freedom, the shortest oscillation period
1228 (as measured from a simulation) is 13~fs for bond-angle vibrations
1229 involving hydrogen atoms. Taking as a guideline that with a Verlet
1230 (leap-frog) integration scheme a minimum of 5 numerical integration
1231 steps should be performed per period of a harmonic oscillation in
1232 order to integrate it with reasonable accuracy, the maximum time step
1233 will be about 3~fs. Disregarding these very fast oscillations of
1234 period 13~fs, the next shortest periods are around 20~fs, which will
1235 allow a maximum time step of about 4~fs.
1237 Removing the bond-angle degrees of freedom from hydrogen atoms can
1238 best be done by defining them as \normindex{virtual interaction sites}
1239 instead of normal atoms. Whereas a normal atom is connected to the molecule
1240 with bonds, angles and dihedrals, a virtual site's position is calculated
1241 from the position of three nearby heavy atoms in a predefined manner
1242 (see also \secref{virtual_sites}). For the hydrogens in water and in
1243 hydroxyl, sulfhydryl, or amine groups, no degrees of freedom can be
1244 removed, because rotational freedom should be preserved. The only
1245 other option available to slow down these motions is to increase the
1246 mass of the hydrogen atoms at the expense of the mass of the connected
1247 heavy atom. This will increase the moment of inertia of the water
1248 molecules and the hydroxyl, sulfhydryl, or amine groups, without
1249 affecting the equilibrium properties of the system and without
1250 affecting the dynamical properties too much. These constructions will
1251 shortly be described in \secref{vsitehydro} and have previously
1252 been described in full detail~\cite{feenstra99}.
1254 Using both virtual sites and \swapindex{modified}{mass}es, the next
1255 bottleneck is likely to be formed by the improper dihedrals (which are
1256 used to preserve planarity or chirality of molecular groups) and the
1257 peptide dihedrals. The peptide dihedral cannot be changed without
1258 affecting the physical behavior of the protein. The improper dihedrals
1259 that preserve planarity mostly deal with aromatic residues. Bonds,
1260 angles, and dihedrals in these residues can also be replaced with
1261 somewhat elaborate virtual site constructions.
1263 All modifications described in this section can be performed using the
1264 {\gromacs} topology building tool {\tt \normindex{pdb2gmx}}. Separate
1265 options exist to increase hydrogen masses, virtualize all hydrogen atoms,
1266 or also virtualize all aromatic residues. {\bf Note} that when all hydrogen
1267 atoms are virtualized, those inside the aromatic residues will be
1268 virtualized as well, {\ie} hydrogens in the aromatic residues are treated
1269 differently depending on the treatment of the aromatic residues.
1271 Parameters for the virtual site constructions for the hydrogen atoms are
1272 inferred from the force-field parameters ({\em vis}. bond lengths and
1273 angles) directly by {\tt \normindex{grompp}} while processing the
1274 topology file. The constructions for the aromatic residues are based
1275 on the bond lengths and angles for the geometry as described in the
1276 force fields, but these parameters are hard-coded into {\tt
1277 \normindex{pdb2gmx}} due to the complex nature of the construction
1278 needed for a whole aromatic group.
1280 \subsection{Hydrogen bond-angle vibrations}
1281 \label{sec:vsitehydro}
1282 \subsubsection{Construction of virtual sites} %%%%%%%%%%%%%%%%%%%%%%%%%
1284 \centerline{\includegraphics[width=11cm]{plots/dumtypes}}
1285 \caption[Virtual site constructions for hydrogen atoms.]{The different
1286 types of virtual site constructions used for hydrogen atoms. The atoms
1287 used in the construction of the virtual site(s) are depicted as black
1288 circles, virtual sites as gray ones. Hydrogens are smaller than heavy
1289 atoms. {\sf A}: fixed bond angle, note that here the hydrogen is not a
1290 virtual site; {\sf B}: in the plane of three atoms, with fixed distance;
1291 {\sf C}: in the plane of three atoms, with fixed angle and distance;
1292 {\sf D}: construction for amine groups ({\amine} or {\aminep}), see
1294 \label{fig:vsitehydro}
1297 The goal of defining hydrogen atoms as virtual sites is to remove all
1298 high-frequency degrees of freedom from them. In some cases, not all
1299 degrees of freedom of a hydrogen atom should be removed, {\eg} in the
1300 case of hydroxyl or amine groups the rotational freedom of the
1301 hydrogen atom(s) should be preserved. Care should be taken that no
1302 unwanted correlations are introduced by the construction of virtual
1303 sites, {\eg} bond-angle vibration between the constructing atoms could
1304 translate into hydrogen bond-length vibration. Additionally, since
1305 virtual sites are by definition massless, in order to preserve total
1306 system mass, the mass of each hydrogen atom that is treated as virtual
1307 site should be added to the bonded heavy atom.
1309 Taking into account these considerations, the hydrogen atoms in a
1310 protein naturally fall into several categories, each requiring a
1311 different approach (see also \figref{vsitehydro}).
1315 \item{\em hydroxyl ({\sf -OH}) or sulfhydryl ({\sf -SH})
1316 hydrogen:\/} The only internal degree of freedom in a hydroxyl group
1317 that can be constrained is the bending of the {\sf C-O-H} angle. This
1318 angle is fixed by defining an additional bond of appropriate length,
1319 see \figref{vsitehydro}A. Doing so removes the high-frequency angle bending,
1320 but leaves the dihedral rotational freedom. The same goes for a
1321 sulfhydryl group. {\bf Note} that in these cases the hydrogen is not treated
1324 \item{\em single amine or amide ({\amines}) and aromatic hydrogens
1325 ({\sf -CH-}):\/} The position of these hydrogens cannot be constructed
1326 from a linear combination of bond vectors, because of the flexibility
1327 of the angle between the heavy atoms. Instead, the hydrogen atom is
1328 positioned at a fixed distance from the bonded heavy atom on a line
1329 going through the bonded heavy atom and a point on the line through
1330 both second bonded atoms, see \figref{vsitehydro}B.
1332 \item{\em planar amine ({\amine}) hydrogens:\/} The method used for
1333 the single amide hydrogen is not well suited for planar amine groups,
1334 because no suitable two heavy atoms can be found to define the
1335 direction of the hydrogen atoms. Instead, the hydrogen is constructed
1336 at a fixed distance from the nitrogen atom, with a fixed angle to the
1337 carbon atom, in the plane defined by one of the other heavy atoms, see
1338 \figref{vsitehydro}C.
1340 \item{\em amine group (umbrella {\amine} or {\aminep}) hydrogens:\/}
1341 Amine hydrogens with rotational freedom cannot be constructed as virtual
1342 sites from the heavy atoms they are connected to, since this would
1343 result in loss of the rotational freedom of the amine group. To
1344 preserve the rotational freedom while removing the hydrogen bond-angle
1345 degrees of freedom, two ``dummy masses'' are constructed with the same
1346 total mass, moment of inertia (for rotation around the {\sf C-N} bond)
1347 and center of mass as the amine group. These dummy masses have no
1348 interaction with any other atom, except for the fact that they are
1349 connected to the carbon and to each other, resulting in a rigid
1350 triangle. From these three particles, the positions of the nitrogen and
1351 hydrogen atoms are constructed as linear combinations of the two
1352 carbon-mass vectors and their outer product, resulting in an amine
1353 group with rotational freedom intact, but without other internal
1354 degrees of freedom. See \figref{vsitehydro}D.
1359 \centerline{\includegraphics[width=15cm]{plots/dumaro}}
1360 \caption[Virtual site constructions for aromatic residues.]{The
1361 different types of virtual site constructions used for aromatic
1362 residues. The atoms used in the construction of the virtual site(s) are
1363 depicted as black circles, virtual sites as gray ones. Hydrogens are
1364 smaller than heavy atoms. {\sf A}: phenylalanine; {\sf B}: tyrosine
1365 (note that the hydroxyl hydrogen is {\em not} a virtual site); {\sf C}:
1366 tryptophan; {\sf D}: histidine.}
1367 \label{fig:vistearo}
1370 \subsection{Out-of-plane vibrations in aromatic groups}
1371 \label{sec:vsitearo}
1372 The planar arrangements in the side chains of the aromatic residues
1373 lends itself perfectly to a virtual-site construction, giving a
1374 perfectly planar group without the inherently unstable constraints
1375 that are necessary to keep normal atoms in a plane. The basic approach
1376 is to define three atoms or dummy masses with constraints between them
1377 to fix the geometry and create the rest of the atoms as simple virtual
1378 sites type (see \secref{virtual_sites}) from these three. Each of
1379 the aromatic residues require a different approach:
1383 \item{\em Phenylalanine:\/} {\sf C}$_\gamma$, {\sf C}$_{{\epsilon}1}$,
1384 and {\sf C}$_{{\epsilon}2}$ are kept as normal atoms, but with each a
1385 mass of one third the total mass of the phenyl group. See
1386 \figref{vsitehydro}A.
1388 \item{\em Tyrosine:\/} The ring is treated identically to the
1389 phenylalanine ring. Additionally, constraints are defined between {\sf
1390 C}$_{{\epsilon}1}$, {\sf C}$_{{\epsilon}2}$, and {\sf O}$_{\eta}$.
1391 The original improper dihedral angles will keep both triangles (one
1392 for the ring and one with {\sf O}$_{\eta}$) in a plane, but due to the
1393 larger moments of inertia this construction will be much more
1394 stable. The bond-angle in the hydroxyl group will be constrained by a
1395 constraint between {\sf C}$_\gamma$ and {\sf H}$_{\eta}$. {\bf Note} that
1396 the hydrogen is not treated as a virtual site. See
1397 \figref{vsitehydro}B.
1399 \item{\em Tryptophan:\/} {\sf C}$_\beta$ is kept as a normal atom
1400 and two dummy masses are created at the center of mass of each of the
1401 rings, each with a mass equal to the total mass of the respective ring
1402 ({\sf C}$_{{\delta}2}$ and {\sf C}$_{{\epsilon}2}$ are each
1403 counted half for each ring). This keeps the overall center of mass and
1404 the moment of inertia almost (but not quite) equal to what it was. See
1405 \figref{vsitehydro}C.
1407 \item{\em Histidine:\/} {\sf C}$_\gamma$, {\sf C}$_{{\epsilon}1}$
1408 and {\sf N}$_{{\epsilon}2}$ are kept as normal atoms, but with masses
1409 redistributed such that the center of mass of the ring is
1410 preserved. See \figref{vsitehydro}D.
1414 \section{Viscosity calculation\index{viscosity}}
1416 The shear viscosity is a property of liquids that can be determined easily
1417 by experiment. It is useful for parameterizing a force field
1418 because it is a kinetic property, while most other properties
1419 which are used for parameterization are thermodynamic.
1420 The viscosity is also an important property, since it influences
1421 the rates of conformational changes of molecules solvated in the liquid.
1423 The viscosity can be calculated from an equilibrium simulation using
1424 an Einstein relation:
1426 \eta = \frac{1}{2}\frac{V}{k_B T} \lim_{t \rightarrow \infty}
1427 \frac{\mbox{d}}{\mbox{d} t} \left\langle
1428 \left( \int_{t_0}^{{t_0}+t} P_{xz}(t') \mbox{d} t' \right)^2
1431 This can be done with {\tt g_energy}.
1432 This method converges very slowly~\cite{Hess2002a}, and as such
1433 a nanosecond simulation might not
1434 be long enough for an accurate determination of the viscosity.
1435 The result is very dependent on the treatment of the electrostatics.
1436 Using a (short) cut-off results in large noise on the off-diagonal
1437 pressure elements, which can increase the calculated viscosity by an order
1440 {\gromacs} also has a non-equilibrium method for determining
1441 the viscosity~\cite{Hess2002a}.
1442 This makes use of the fact that energy, which is fed into system by
1443 external forces, is dissipated through viscous friction. The generated heat
1444 is removed by coupling to a heat bath. For a Newtonian liquid adding a
1445 small force will result in a velocity gradient according to the following
1448 a_x(z) + \frac{\eta}{\rho} \frac{\partial^2 v_x(z)}{\partial z^2} = 0
1450 Here we have applied an acceleration $a_x(z)$ in the $x$-direction, which
1451 is a function of the $z$-coordinate.
1452 In {\gromacs} the acceleration profile is:
1454 a_x(z) = A \cos\left(\frac{2\pi z}{l_z}\right)
1456 where $l_z$ is the height of the box. The generated velocity profile is:
1458 v_x(z) = V \cos\left(\frac{2\pi z}{l_z}\right)
1461 V = A \frac{\rho}{\eta}\left(\frac{l_z}{2\pi}\right)^2
1463 The viscosity can be calculated from $A$ and $V$:
1466 \eta = \frac{A}{V}\rho \left(\frac{l_z}{2\pi}\right)^2
1469 In the simulation $V$ is defined as:
1471 V = \frac{\displaystyle \sum_{i=1}^N m_i v_{i,x} 2 \cos\left(\frac{2\pi z}{l_z}\right)}
1472 {\displaystyle \sum_{i=1}^N m_i}
1474 The generated velocity profile is not coupled to the heat bath. Moreover,
1475 the velocity profile is excluded from the kinetic energy.
1476 One would like $V$ to be as large as possible to get good statistics.
1477 However, the shear rate should not be so high that the system gets too far
1478 from equilibrium. The maximum shear rate occurs where the cosine is zero,
1481 \mbox{sh}_{\max} = \max_z \left| \frac{\partial v_x(z)}{\partial z} \right|
1482 = A \frac{\rho}{\eta} \frac{l_z}{2\pi}
1484 For a simulation with: $\eta=10^{-3}$ [kg\,m$^{-1}$\,s$^{-1}$],
1485 $\rho=10^3$\,[kg\,m$^{-3}$] and $l_z=2\pi$\,[nm],
1486 $\mbox{sh}_{\max}=1$\,[ps\,nm$^{-1}$] $A$.
1487 This shear rate should be smaller than one over the longest
1488 correlation time in the system. For most liquids, this will be the rotation
1489 correlation time, which is around 10 ps. In this case, $A$ should
1490 be smaller than 0.1\,[nm\,ps$^{-2}$].
1491 When the shear rate is too high, the observed viscosity will be too low.
1492 Because $V$ is proportional to the square of the box height,
1493 the optimal box is elongated in the $z$-direction.
1494 In general, a simulation length of 100 ps is enough to obtain an
1495 accurate value for the viscosity.
1497 The heat generated by the viscous friction is removed by coupling to a heat
1498 bath. Because this coupling is not instantaneous the real temperature of the
1499 liquid will be slightly lower than the observed temperature.
1500 Berendsen derived this temperature shift~\cite{Berendsen91}, which can
1501 be written in terms of the shear rate as:
1503 T_s = \frac{\eta\,\tau}{2 \rho\,C_v} \mbox{sh}_{\max}^2
1505 where $\tau$ is the coupling time for the Berendsen thermostat and
1506 $C_v$ is the heat capacity. Using the values of the example above,
1507 $\tau=10^{-13}$ [s] and $C_v=2 \cdot 10^3$\,[J kg$^{-1}$\,K$^{-1}$], we
1508 get: $T_s=25$\,[K\,ps$^{-2}$]\,sh$_{\max}^2$. When we want the shear
1509 rate to be smaller than $1/10$\,[ps$^{-1}$], $T_s$ is smaller than
1510 0.25\,[K], which is negligible.
1512 {\bf Note} that the system has to build up the velocity profile when starting
1513 from an equilibrium state. This build-up time is of the order of the
1514 correlation time of the liquid.
1516 Two quantities are written to the energy file, along with their averages
1517 and fluctuations: $V$ and $1/\eta$, as obtained from (\ref{visc}).
1519 \section{Tabulated interaction functions\index{tabulated interaction functions}}
1520 \subsection{Cubic splines for potentials}
1521 \label{subsec:cubicspline}
1522 In some of the inner loops of {\gromacs}, look-up tables are used
1523 for computation of potential and forces.
1524 The tables are interpolated using a cubic
1526 There are separate tables for electrostatic, dispersion, and repulsion
1528 but for the sake of caching performance these have been combined
1529 into a single array.
1530 The cubic spline interpolation for $x_i \leq x < x_{i+1}$ looks like this:
1532 V_s(x) = A_0 + A_1 \,\epsilon + A_2 \,\epsilon^2 + A_3 \,\epsilon^3
1535 where the table spacing $h$ and fraction $\epsilon$ are given by:
1537 h &=& x_{i+1} - x_i \\
1538 \epsilon&=& (x - x_i)/h
1540 so that $0 \le \epsilon < 1$.
1541 From this, we can calculate the derivative in order to determine the forces:
1544 -\frac{{\rm d}V_s(x)}{{\rm d}\epsilon}\frac{{\rm d}\epsilon}{{\rm d}x} ~=~
1545 -(A_1 + 2 A_2 \,\epsilon + 3 A_3 \,\epsilon^2)/h
1547 The four coefficients are determined from the four conditions
1548 that $V_s$ and $-V_s'$ at both ends of each interval should match
1549 the exact potential $V$ and force $-V'$.
1550 This results in the following errors for each interval:
1552 |V_s - V |_{max} &=& V'''' \frac{h^4}{384} + O(h^5) \\
1553 |V_s' - V' |_{max} &=& V'''' \frac{h^3}{72\sqrt{3}} + O(h^4) \\
1554 |V_s''- V''|_{max} &=& V'''' \frac{h^2}{12} + O(h^3)
1556 V and V' are continuous, while V'' is the first discontinuous
1558 The number of points per nanometer is 500 and 2000
1559 for mixed- and double-precision versions of {\gromacs}, respectively.
1560 This means that the errors in the potential and force will usually
1561 be smaller than the mixed precision accuracy.
1563 {\gromacs} stores $A_0$, $A_1$, $A_2$ and $A_3$.
1564 The force routines get a table with these four parameters and
1565 a scaling factor $s$ that is equal to the number of points per nm.
1566 ({\bf Note} that $h$ is $s^{-1}$).
1567 The algorithm goes a little something like this:
1569 \item Calculate distance vector (\ve{r}$_{ij}$) and distance r$_{ij}$
1570 \item Multiply r$_{ij}$ by $s$ and truncate to an integer value $n_0$
1571 to get a table index
1572 \item Calculate fractional component ($\epsilon$ = $s$r$_{ij} - n_0$)
1574 \item Do the interpolation to calculate the potential $V$ and the scalar force $f$
1575 \item Calculate the vector force \ve{F} by multiplying $f$ with \ve{r}$_{ij}$
1578 {\bf Note} that table look-up is significantly {\em
1579 slower} than computation of the most simple Lennard-Jones and Coulomb
1580 interaction. However, it is much faster than the shifted Coulomb
1581 function used in conjunction with the PPPM method. Finally, it is much
1582 easier to modify a table for the potential (and get a graphical
1583 representation of it) than to modify the inner loops of the MD
1586 \subsection{User-specified potential functions}
1587 \label{subsec:userpot}
1588 You can also use your own \index{potential function}s
1589 without editing the {\gromacs} code.
1590 The potential function should be according to the following equation
1592 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})
1594 where $f$, $g$, and $h$ are user defined functions. {\bf Note} that if $g(r)$ represents a
1595 normal dispersion interaction, $g(r)$ should be $<$ 0. C$_6$, C$_{12}$
1596 and the charges are read from the topology. Also note that combination
1597 rules are only supported for Lennard-Jones and Buckingham, and that
1598 your tables should match the parameters in the binary topology.
1600 When you add the following lines in your {\tt .mdp} file:
1611 {\tt mdrun} will read a single non-bonded table file,
1612 or multiple when {\tt energygrp-table} is set (see below).
1613 The name of the file(s) can be set with the {\tt mdrun} option {\tt -table}.
1614 The table file should contain seven columns of table look-up data in the
1615 order: $x$, $f(x)$, $-f'(x)$, $g(x)$, $-g'(x)$, $h(x)$, $-h'(x)$.
1616 The $x$ should run from 0 to $r_c+1$ (the value of {\tt table_extension} can be
1617 changed in the {\tt .mdp} file).
1618 You can choose the spacing you like; for the standard tables {\gromacs}
1619 uses a spacing of 0.002 and 0.0005 nm when you run in mixed
1620 and double precision, respectively. In this
1621 context, $r_c$ denotes the maximum of the two cut-offs {\tt rvdw} and
1622 {\tt rcoulomb} (see above). These variables need not be the same (and
1623 need not be 1.0 either). Some functions used for potentials contain a
1624 singularity at $x = 0$, but since atoms are normally not closer to each
1625 other than 0.1 nm, the function value at $x = 0$ is not important.
1627 possible to combine a standard Coulomb with a modified LJ potential
1628 (or vice versa). One then specifies {\eg} {\tt coulombtype = Cut-off} or
1629 {\tt coulombtype = PME}, combined with {\tt vdwtype = User}. The table file must
1630 always contain the 7 columns however, and meaningful data (i.e. not
1631 zeroes) must be entered in all columns. A number of pre-built table
1632 files can be found in the {\tt GMXLIB} directory for 6-8, 6-9, 6-10, 6-11, and 6-12
1633 Lennard-Jones potentials combined with a normal Coulomb.
1635 If you want to have different functional forms between different
1636 groups of atoms, this can be set through energy groups.
1637 Different tables can be used for non-bonded interactions between
1638 different energy groups pairs through the {\tt .mdp} option {\tt energygrp-table}
1639 (see details in the User Guide).
1640 Atoms that should interact with a different potential should
1641 be put into different energy groups.
1642 Between group pairs which are not listed in {\tt energygrp-table},
1643 the normal user tables will be used. This makes it easy to use
1644 a different functional form between a few types of atoms.
1646 \section{Mixed Quantum-Classical simulation techniques}
1648 In a molecular mechanics (MM) force field, the influence of electrons
1649 is expressed by empirical parameters that are assigned on the basis of
1650 experimental data, or on the basis of results from high-level quantum
1651 chemistry calculations. These are valid for the ground state of a
1652 given covalent structure, and the MM approximation is usually
1653 sufficiently accurate for ground-state processes in which the overall
1654 connectivity between the atoms in the system remains
1655 unchanged. However, for processes in which the connectivity does
1656 change, such as chemical reactions, or processes that involve multiple
1657 electronic states, such as photochemical conversions, electrons can no
1658 longer be ignored, and a quantum mechanical description is required
1659 for at least those parts of the system in which the reaction takes
1662 One approach to the simulation of chemical reactions in solution, or
1663 in enzymes, is to use a combination of quantum mechanics (QM) and
1664 molecular mechanics (MM). The reacting parts of the system are treated
1665 quantum mechanically, with the remainder being modeled using the
1666 force field. The current version of {\gromacs} provides interfaces to
1667 several popular Quantum Chemistry packages (MOPAC~\cite{mopac},
1668 GAMESS-UK~\cite{gamess-uk}, Gaussian~\cite{g03} and CPMD~\cite{Car85a}).
1670 {\gromacs} interactions between the two subsystems are
1671 either handled as described by Field {\em et al.}~\cite{Field90a} or
1672 within the ONIOM approach by Morokuma and coworkers~\cite{Maseras96a,
1675 \subsection{Overview}
1677 Two approaches for describing the interactions between the QM and MM
1678 subsystems are supported in this version:
1681 \item{\textbf{Electronic Embedding}} The electrostatic interactions
1682 between the electrons of the QM region and the MM atoms and between
1683 the QM nuclei and the MM atoms are included in the Hamiltonian for
1684 the QM subsystem: \beq H^{QM/MM} =
1685 H^{QM}_e-\sum_i^n\sum_J^M\frac{e^2Q_J}{4\pi\epsilon_0r_{iJ}}+\sum_A^N\sum_J^M\frac{e^2Z_AQ_J}{e\pi\epsilon_0R_{AJ}},
1686 \eeq where $n$ and $N$ are the number of electrons and nuclei in the
1687 QM region, respectively, and $M$ is the number of charged MM
1688 atoms. The first term on the right hand side is the original
1689 electronic Hamiltonian of an isolated QM system. The first of the
1690 double sums is the total electrostatic interaction between the QM
1691 electrons and the MM atoms. The total electrostatic interaction of the
1692 QM nuclei with the MM atoms is given by the second double sum. Bonded
1693 interactions between QM and MM atoms are described at the MM level by
1694 the appropriate force-field terms. Chemical bonds that connect the two
1695 subsystems are capped by a hydrogen atom to complete the valence of
1696 the QM region. The force on this atom, which is present in the QM
1697 region only, is distributed over the two atoms of the bond. The cap
1698 atom is usually referred to as a link atom.
1700 \item{\textbf{ONIOM}} In the ONIOM approach, the energy and gradients
1701 are first evaluated for the isolated QM subsystem at the desired level
1702 of {\it{ab initio}} theory. Subsequently, the energy and gradients of
1703 the total system, including the QM region, are computed using the
1704 molecular mechanics force field and added to the energy and gradients
1705 calculated for the isolated QM subsystem. Finally, in order to correct
1706 for counting the interactions inside the QM region twice, a molecular
1707 mechanics calculation is performed on the isolated QM subsystem and
1708 the energy and gradients are subtracted. This leads to the following
1709 expression for the total QM/MM energy (and gradients likewise): \beq
1710 E_{tot} = E_{I}^{QM}
1711 +E_{I+II}^{MM}-E_{I}^{MM}, \eeq where the
1712 subscripts I and II refer to the QM and MM subsystems,
1713 respectively. The superscripts indicate at what level of theory the
1714 energies are computed. The ONIOM scheme has the
1715 advantage that it is not restricted to a two-layer QM/MM description,
1716 but can easily handle more than two layers, with each layer described
1717 at a different level of theory.
1722 To make use of the QM/MM functionality in {\gromacs}, one needs to:
1725 \item introduce link atoms at the QM/MM boundary, if needed;
1726 \item specify which atoms are to be treated at a QM level;
1727 \item specify the QM level, basis set, type of QM/MM interface and so on.
1730 \subsubsection{Adding link atoms}
1732 At the bond that connects the QM and MM subsystems, a link atoms is
1733 introduced. In {\gromacs} the link atom has special atomtype, called
1734 LA. This atomtype is treated as a hydrogen atom in the QM calculation,
1735 and as a virtual site in the force-field calculation. The link atoms, if
1736 any, are part of the system, but have no interaction with any other
1737 atom, except that the QM force working on it is distributed over the
1738 two atoms of the bond. In the topology, the link atom (LA), therefore,
1739 is defined as a virtual site atom:
1744 LA QMatom MMatom 1 0.65
1747 See~\secref{vsitetop} for more details on how virtual sites are
1748 treated. The link atom is replaced at every step of the simulation.
1750 In addition, the bond itself is replaced by a constraint:
1755 QMatom MMatom 2 0.153
1758 {\bf Note} that, because in our system the QM/MM bond is a carbon-carbon
1759 bond (0.153 nm), we use a constraint length of 0.153 nm, and dummy
1760 position of 0.65. The latter is the ratio between the ideal C-H
1761 bond length and the ideal C-C bond length. With this ratio, the link
1762 atom is always 0.1 nm away from the {\tt QMatom}, consistent with the
1763 carbon-hydrogen bond length. If the QM and MM subsystems are connected
1764 by a different kind of bond, a different constraint and a different
1765 dummy position, appropriate for that bond type, are required.
1767 \subsubsection{Specifying the QM atoms}
1769 Atoms that should be treated at a QM level of theory, including the
1770 link atoms, are added to the index file. In addition, the chemical
1771 bonds between the atoms in the QM region are to be defined as
1772 connect bonds (bond type 5) in the topology file:
1781 \subsubsection{Specifying the QM/MM simulation parameters}
1783 In the {\tt .mdp} file, the following parameters control a QM/MM simulation.
1787 \item[\tt QMMM = no]\mbox{}\\ If this is set to {\tt yes}, a QM/MM
1788 simulation is requested. Several groups of atoms can be described at
1789 different QM levels separately. These are specified in the QMMM-grps
1790 field separated by spaces. The level of {\it{ab initio}} theory at which the
1791 groups are described is specified by {\tt QMmethod} and {\tt QMbasis}
1792 Fields. Describing the groups at different levels of theory is only
1793 possible with the ONIOM QM/MM scheme, specified by {\tt QMMMscheme}.
1795 \item[\tt QMMM-grps =]\mbox{}\\groups to be described at the QM level
1797 \item[\tt QMMMscheme = normal]\mbox{}\\Options are {\tt normal} and
1798 {\tt ONIOM}. This selects the QM/MM interface. {\tt normal} implies
1799 that the QM subsystem is electronically embedded in the MM
1800 subsystem. There can only be one {\tt QMMM-grps} that is modeled at
1801 the {\tt QMmethod} and {\tt QMbasis} level of {\it{ ab initio}}
1802 theory. The rest of the system is described at the MM level. The QM
1803 and MM subsystems interact as follows: MM point charges are included
1804 in the QM one-electron Hamiltonian and all Lennard-Jones interactions
1805 are described at the MM level. If {\tt ONIOM} is selected, the
1806 interaction between the subsystem is described using the ONIOM method
1807 by Morokuma and co-workers. There can be more than one QMMM-grps each
1808 modeled at a different level of QM theory (QMmethod and QMbasis).
1810 \item[\tt QMmethod = ]\mbox{}\\Method used to compute the energy
1811 and gradients on the QM atoms. Available methods are AM1, PM3, RHF,
1812 UHF, DFT, B3LYP, MP2, CASSCF, MMVB and CPMD. For CASSCF, the number of
1813 electrons and orbitals included in the active space is specified by
1814 {\tt CASelectrons} and {\tt CASorbitals}. For CPMD, the plane-wave
1815 cut-off is specified by the {\tt planewavecutoff} keyword.
1817 \item[\tt QMbasis = ]\mbox{}\\Gaussian basis set used to expand the
1818 electronic wave-function. Only Gaussian basis sets are currently
1819 available, i.e. STO-3G, 3-21G, 3-21G*, 3-21+G*, 6-21G, 6-31G, 6-31G*,
1820 6-31+G*, and 6-311G. For CPMD, which uses plane wave expansion rather
1821 than atom-centered basis functions, the {\tt planewavecutoff} keyword
1822 controls the plane wave expansion.
1824 \item[\tt QMcharge = ]\mbox{}\\The total charge in {\it{e}} of the {\tt
1825 QMMM-grps}. In case there are more than one {\tt QMMM-grps}, the total
1826 charge of each ONIOM layer needs to be specified separately.
1828 \item[\tt QMmult = ]\mbox{}\\The multiplicity of the {\tt
1829 QMMM-grps}. In case there are more than one {\tt QMMM-grps}, the
1830 multiplicity of each ONIOM layer needs to be specified separately.
1832 \item[\tt CASorbitals = ]\mbox{}\\The number of orbitals to be
1833 included in the active space when doing a CASSCF computation.
1835 \item[\tt CASelectrons = ]\mbox{}\\The number of electrons to be
1836 included in the active space when doing a CASSCF computation.
1838 \item[\tt SH = no]\mbox{}\\If this is set to yes, a QM/MM MD
1839 simulation on the excited state-potential energy surface and enforce a
1840 diabatic hop to the ground-state when the system hits the conical
1841 intersection hyperline in the course the simulation. This option only
1842 works in combination with the CASSCF method.
1848 The energies and gradients computed in the QM calculation are added to
1849 those computed by {\gromacs}. In the {\tt .edr} file there is a section
1850 for the total QM energy.
1852 \subsection{Future developments}
1854 Several features are currently under development to increase the
1855 accuracy of the QM/MM interface. One useful feature is the use of
1856 delocalized MM charges in the QM computations. The most important
1857 benefit of using such smeared-out charges is that the Coulombic
1858 potential has a finite value at interatomic distances. In the point
1859 charge representation, the partially-charged MM atoms close to the QM
1860 region tend to ``over-polarize'' the QM system, which leads to artifacts
1863 What is needed as well is a transition state optimizer.
1865 \section{\normindex{Adaptive Resolution Scheme}}
1866 \newcommand{\adress}{AdResS}
1867 The adaptive resolution scheme~\cite{Praprotnik2005,Praprotnik2008} (\seeindex{\adress}{Adaptive Resolution Scheme}) couples two systems with different resolutions by a force interpolation scheme.
1868 In contrast to the mixed Quantum-Classical simulation techniques of the previous section, the number of high resolution particles is not fixed, but can vary over the simulation time.
1870 Below we discuss {\adress} for a double resolution (atomistic and coarse grained) representation of the same system. See \figref{adress} for illustration.
1871 The details of implementation described in this section were published in~\cite{Junghans2010,Fritsch2012}.
1875 \includegraphics[width=0.5\textwidth]{plots/adress}
1876 \caption{A schematic illustration of the {\adress} method for water.}
1880 Every molecule needs a well-defined mapping point (usually the center of mass)
1881 but any other linear combination of particle coordinates is also sufficient.
1882 In the topology the mapping point is defined by a virtual site. The forces in the coarse-grained region are functions of the mapping point positions only.
1883 In this implementation molecules are modeled by charge groups or sets of charge groups, which actually allows one to have multiple mapping points per molecule. This can be useful for bigger molecules like polymers. In that case one has to also extend the AdResS description to bonded interactions~\cite{Praprotnik2011}, which will be implemented into \gromacs in one of the future versions.
1885 The force between two molecules is given by~\cite{Praprotnik2005}
1886 \footnote{Note that the equation obeys Newton's third law, which is not the case for other interpolation schemes~\cite{DelleSite2007}.}:
1888 \vec{F}_{\alpha\beta}=w_\alpha w_\beta \vec{F}^\mathrm{ex,mol}_{\alpha\beta} + \left[1-w_\alpha w_\beta\right] \vec{F}^\mathrm{cg,mol}_{\alpha\beta}~,
1889 \label{eqn:interpolation}
1891 where $\alpha$ and $\beta$ label the two molecules and $w_\alpha$, $w_\beta$ are the adaptive weights of the two molecules.
1893 The first part, which represents the explicit interaction of the molecules, can be written as:
1895 \vec{F}^\mathrm{ex,mol}_{\alpha\beta}=\sum_{i\in\alpha}\sum_{j\in\beta} \vec{F}^\mathrm{ex}_{ij}~,
1897 where $\vec{F}^\mathrm{ex}_{ij}$ is the force between the $i$th atom in $\alpha$th molecule and the $j$th atom in the $\beta$th molecule, which is given by an explicit force field.
1898 The second part of \eqnref{interpolation} comes from the coarse-grained interaction of the molecules.
1899 In \gromacs a slightly extended case is implemented:
1901 \vec{F}_{\alpha\beta}=\sum_{i\in\alpha}\sum_{j\in\beta} w_i w_j \vec{F}^\mathrm{ex}_{ij} + \left[1-w_\alpha w_\beta\right] \vec{F}^\mathrm{cg,mol}_{\alpha\beta}~,
1902 \label{eqn:interpolation2}
1904 where $w_i$ and $w_j$ are atom-wise weights, which are determined by the {\tt adress-site} option. For {\tt adress-site} being the center of mass, atom $i$ has the weight of the center of mass of its \emph{charge group}.
1905 The weight $w_\alpha$ of molecule $\alpha$ is determined by the position of coarse-grained particle, which is constructed as a virtual site from the atomistic particles as specified in the topology.
1906 This extension allows one to perform all kind of AdResS variations, but the common case can be recovered by using a center of mass virtual site in the topology, {\tt adress-site=COM} and putting all atoms (except the virtual site representing the coarse-grained interaction) of a molecule into one charge group.
1907 For big molecules, it is sometimes useful to use an atom-based weight, which can be either be achieved by setting {\tt adress-site=atomperatom} or putting every atom into a separate charge group (the center of mass of a charge group with one atom is the atom itself).
1909 The coarse-grained force field $\vec{F}^\mathrm{cg}$ is usually derived from the atomistic system by structure-based coarse-graining (see \secref{cg-forcefields}). To specify which atoms belong to a coarse-grained representation, energy groups are used.
1910 Each coarse-grained interaction has to be associated with a specific energy group, which is why the virtual sites representing the coarse-grained interactions also have to be in different charge groups. The energy groups which are treated as coarse-grained interactions are then listed in {\tt adress_cg_grp_names}.
1911 The most important element of this interpolation (see \eqnref{interpolation} and \eqnref{interpolation2}) is the adaptive weighting function (for illustration see \figref{adress}):
1914 \left\{\begin{array}{c@{\;:\;}l}
1915 1&\mathrm{atomistic/explicit\;region}\\
1916 0<w<1&\mathrm{hybrid\;region}\\
1917 0&\mathrm{coarse-grained\;region}
1918 \end{array}\right.~,
1919 \label{equ:weighting}
1921 which has a value between 0 and 1.
1922 This definition of $w$ gives a purely explicit force in the explicit region and a purely coarse-grained force in the coarse-grained region,
1923 so essentially \eqnref{interpolation} only the hybrid region has mixed interactions which would not appear in a standard simulation.
1924 In {\gromacs}, a $\cos^2$-like function is implemented as a weighting function:
1927 \left\{\begin{array}{c@{\;:\;}r@{x}l}
1928 0&&>d_\mathrm{ex}+d_\mathrm{hy}\\
1929 \cos^2\left(\frac{\pi}{2d_\mathrm{hy}}(x-d_\mathrm{ex})\right)&d_\mathrm{ex}+d_\mathrm{hy}>&>d_\mathrm{ex}\\
1931 \end{array}\right.~,
1934 where $d_\mathrm{ex}$ and $d_\mathrm{hy}$ are the sizes of the explicit and the hybrid region, respectively.
1935 Depending on the physical interest of the research, other functions could be implemented as long as the following boundary conditions are fulfilled:
1936 The function is 1) continuous, 2) monotonic and 3) has zero derivatives at the boundaries.
1937 Spherical and one-dimensional splitting of the simulation box has been implemented ({\tt adress-type} option)
1938 and depending on this, the distance $x$ to the center of the explicit region is calculated as follows:
1942 \begin{array}{c@{\;:\;}l}
1943 |(\vec{R}_\alpha-\vec{R}_\mathrm{ct})\cdot\hat{e}|&\mathrm{splitting\;in\;}\hat{e}\mathrm{\;direction}\\
1944 |\vec{R}_\alpha-\vec{R}_\mathrm{ct}|&\mathrm{spherical\;splitting}
1948 where $\vec{R}_\mathrm{ct}$ is the center of the explicit zone (defined by {\tt adress-reference-coords} option). $\vec{R}_\alpha$ is the mapping point of the $\alpha$th molecule. For the center of mass mapping, it is given by:
1950 R_\alpha=\frac{\sum_{i\in\alpha}m_i r_i}{\sum_{i\in\alpha}m_i}
1953 Note that the value of the weighting function depends exclusively on the mapping of the molecule.
1955 The interpolation of forces (see \eqnref{interpolation2}) can produce inhomogeneities in the density and affect the structure of the system in the hybrid region.
1957 One way of reducing the density inhomogeneities is by the application of the so-called thermodynamic force (TF)~\cite{Poblete2010}.
1958 Such a force consists of a space-dependent external field applied in the hybrid region on the coarse-grained site of each molecule. It can be specified for each of the species of the system.
1959 The TF compensates the pressure profile~\cite{Fritsch2012b} that emerges under a homogeneous density profile. Therefore, it can correct the local density inhomogeneities in the hybrid region and it also allows the coupling of atomistic and coarse-grained representations which by construction have different pressures at the target density.
1960 The field can be determined by an iterative procedure, which is described in detail in the \href{http://code.google.com/p/votca/downloads/list?&q=manual}{manual} of the \normindex{VOTCA package}~\cite{ruehle2009}. Setting the {\tt adress-interface-correction} to {\tt thermoforce} enables the TF correction and\newline{\tt adress-tf-grp-names} defines the energy groups to act on.
1962 \subsection{Example: Adaptive resolution simulation of water}\label{subsec:adressexample}
1963 In this section the set up of an adaptive resolution simulation coupling atomistic SPC ~\cite{Berendsen81} water to its coarse-grained representation will be explained (as used in \cite{Fritsch2012b}).
1964 The following steps are required to setup the simulation:
1966 \item Perform a reference all-atom simulation
1967 \item Create a coarse-grained representation and save it as tabulated interaction function
1968 \item Create a hybrid topology for the SPC water
1969 \item Modify the atomistic coordinate file to include the coarse grained representation
1970 \item Define the geometry of the adaptive simulation in the grompp input file
1971 \item Create an index file
1973 The coarse-grained representation of the interaction is stored as tabulated interaction function see \ssecref{userpot}. The convention is to use the $C^{(12)}$ columns with the $C^{(12)}$- coefficient set to 1. All other columns should be zero. The VOTCA manual has detailed instructions and a tutorial for SPC water on how to coarse-grain the interaction using various techniques.
1974 Here we named the coarse grained interaction CG, so the corresponding tabulated file is {\tt table_CG_CG.xvg}. To create the topology one can start from the atomistic topology file (e.g. share/gromacs/top/oplsaa.ff/spc.itp), we are assuming rigid water here. In the VOTCA tutorial the file is named {\tt hybrid_spc.itp}.
1975 The only difference to the atomistic topology is the addition of a coarse-grained virtual site:
1983 ; nr type resnr residue atom cgnr charge mass
1984 1 opls_116 1 SOL OW 1 -0.82
1985 2 opls_117 1 SOL HW1 1 0.41
1986 3 opls_117 1 SOL HW2 1 0.41
2000 ; Site from funct a d
2001 4 1 2 3 1 0.05595E+00 0.05595E+00
2003 The virtual site type 3 with the specified coefficients places the virtual site in the center of mass of the molecule (for larger molecules virtual_sitesn has to be used).
2004 We now need to include our modified water model in the topology file and define the type {\tt CG}. In {\tt topol.top}:
2007 #include "ffoplsaa.itp"
2010 ;name mass charge ptype sigma epsilon
2011 CG 0.00000 0.0000 V 1 0.25
2013 #include "hybrid_spc.itp"
2020 The $\sigma$ and $\epsilon$ values correspond to $C_6=1$ and $C_{12}=1$ and thus the table file should contain the coarse-grained interaction in either the $C_6$ or $C_{12}$ column. In the example the OPLS force field is used where $\sigma$ and $\epsilon$ are specified.
2021 Note that for force fields which define atomtypes directly in terms of $C_6$ and $C_{12}$, one can simply set $C_6=0$ and $C_{12}=1$. See section \ssecref{userpot} for more details on tabulated interactions. Since now the water molecule has a virtual site the coordinate file also needs to include that.
2024 adaptive water coordinates
2026 1SOL OW 1 0.283 0.886 0.647
2027 1SOL HW1 2 0.359 0.884 0.711
2028 1SOL HW2 3 0.308 0.938 0.566
2029 1SOL CG 4 0.289 0.889 0.646
2030 1SOL OW 5 1.848 0.918 0.082
2031 1SOL HW1 6 1.760 0.930 0.129
2032 1SOL HW2 7 1.921 0.912 0.150
2033 1SOL CG 8 1.847 0.918 0.088
2036 This file can be created manually or using the VOTCA tool {\tt csg_map } with the {\tt --hybrid} option.\\
2037 In the grompp input file the AdResS feature needs to be enabled and the geometry defined.
2041 ; AdResS relevant options
2043 energygrp_table = CG CG
2045 ; Method for doing Van der Waals
2049 adress_type = xsplit
2050 adress_ex_width = 1.5
2051 adress_hy_width = 1.5
2052 adress_interface_correction = off
2053 adress_reference_coords = 8 0 0
2054 adress_cg_grp_names = CG
2057 Here we are defining an energy group {\tt CG} which consists of the coarse-grained virtual site.
2058 As discussed above, the coarse-grained interaction is usually tabulated. This requires the {\tt vdw-type} parameter to be set to {\tt user}. In the case where multi-component systems are coarse-grained, an energy group has to be defined for each component. Note that all the energy groups defining coarse-grained representations have to be listed again in {\tt adress_cg_grp_names} to distinguish them from regular energy groups.\\
2059 The index file has to be updated to have a group CG which includes all the coarse-grained virtual sites. This can be done easily using the {\tt make_ndx} tool of gromacs.
2061 \section{Using VMD plug-ins for trajectory file I/O}
2062 \index{VMD plug-ins}\index{trajectory file}{\gromacs} tools are able
2063 to use the plug-ins found in an existing installation of
2064 \href{http://www.ks.uiuc.edu/Research/vmd}{VMD} in order to read and
2065 write trajectory files in formats that are not native to
2066 {\gromacs}. You will be able to supply an AMBER DCD-format trajectory
2067 filename directly to {\gromacs} tools, for example.
2069 This requires a VMD installation not older than version 1.8, that your
2070 system provides the dlopen function so that programs can determine at
2071 run time what plug-ins exist, and that you build shared libraries when
2072 building {\gromacs}. CMake will find the vmd executable in your path, and
2073 from it, or the environment variable {\tt VMDDIR} at configuration or
2074 run time, locate the plug-ins. Alternatively, the {\tt VMD_PLUGIN_PATH}
2075 can be used at run time to specify a path where these plug-ins can be
2076 found. Note that these plug-ins are in a binary format, and that format
2077 must match the architecture of the machine attempting to use them.
2080 \section{\normindex{Interactive Molecular Dynamics}}
2081 {\gromacs} supports the interactive molecular dynamics (IMD) protocol as implemented
2082 by \href{http://www.ks.uiuc.edu/Research/vmd}{VMD} to control a running simulation
2083 in NAMD. IMD allows to monitor a running {\gromacs} simulation from a VMD client.
2084 In addition, the user can interact with the simulation by pulling on atoms, residues
2085 or fragments with a mouse or a force-feedback device. Additional information about
2086 the {\gromacs} implementation and an exemplary {\gromacs} IMD system can be found
2087 \href{http://www.mpibpc.mpg.de/grubmueller/interactivemd}{on this homepage}.
2089 \subsection{Simulation input preparation}
2090 The {\gromacs} implementation allows transmission and interaction with a part of the
2091 running simulation only, e.g.\ in cases where no water molecules should be transmitted
2092 or pulled. The group is specified via the {\tt .mdp} option {\tt IMD-group}. When
2093 {\tt IMD-group} is empty, the IMD protocol is disabled and cannot be enabled via the
2094 switches in {\tt mdrun}. To interact with the entire system, {\tt IMD-group} can
2095 be set to {\tt System}. When using {\tt grompp}, a {\tt .gro} file
2096 to be used as VMD input is written out ({\tt -imd} switch of {\tt grompp}).
2098 \subsection{Starting the simulation}
2099 Communication between VMD and {\gromacs} is achieved via TCP sockets and thus enables
2100 controlling an {\tt mdrun} running locally or on a remote cluster. The port for the
2101 connection can be specified with the {\tt -imdport} switch of {\tt mdrun}, 8888 is
2102 the default. If a port number of 0 or smaller is provided, {\gromacs} automatically
2103 assigns a free port to use with IMD.
2105 Every $N$ steps, the {\tt mdrun} client receives the applied forces from VMD and sends the new
2106 positions to the client. VMD permits increasing or decreasing the communication frequency
2108 By default, the simulation starts and runs even if no IMD client is connected. This
2109 behavior is changed by the {\tt -imdwait} switch of {\tt mdrun}. After startup and
2110 whenever the client has disconnected, the integration stops until reconnection of
2112 When the {\tt -imdterm} switch is used, the simulation can be terminated by pressing
2113 the stop button in VMD. This is disabled by default.
2114 Finally, to allow interacting with the simulation (i.e. pulling from VMD) the {\tt -imdpull}
2115 switch has to be used.
2116 Therefore, a simulation can only be monitored but not influenced from the VMD client
2117 when none of {\tt -imdwait}, {\tt -imdterm} or {\tt -imdpull} are set. However, since
2118 the IMD protocol requires no authentication, it is not advisable to run simulations on
2119 a host directly reachable from an insecure environment. Secure shell forwarding of TCP
2120 can be used to connect to running simulations not directly reachable from the interacting host.
2121 Note that the IMD command line switches of {\tt mdrun} are hidden by default and show
2122 up in the help text only with {\tt gmx mdrun -h -hidden}.
2124 \subsection{Connecting from VMD}
2125 In VMD, first the structure corresponding to the IMD group has to be loaded ({\it File
2126 $\rightarrow$ New Molecule}). Then the IMD connection window has to be used
2127 ({\it Extensions $\rightarrow$ Simulation $\rightarrow$ IMD Connect (NAMD)}). In the IMD
2128 connection window, hostname and port have to be specified and followed by pressing
2129 {\it Connect}. {\it Detach Sim} allows disconnecting without terminating the simulation, while
2130 {\it Stop Sim} ends the simulation on the next neighbor searching step (if allowed by
2133 The timestep transfer rate allows adjusting the communication frequency between simulation
2134 and IMD client. Setting the keep rate loads every $N^\mathrm{th}$ frame into VMD instead
2135 of discarding them when a new one is received. The displayed energies are in SI units
2136 in contrast to energies displayed from NAMD simulations.
2139 % LocalWords: PMF pmf kJ mol Jarzynski BT bilayer rup mdp AFM fepmf fecalc rb
2140 % LocalWords: posre grompp fs Verlet dihedrals hydrogens hydroxyl sulfhydryl
2141 % LocalWords: vsitehydro planarity chirality pdb gmx virtualize virtualized xz
2142 % LocalWords: vis massless tryptophan histidine phenyl parameterizing ij PPPM
2143 % LocalWords: parameterization Berendsen rlist coulombtype rcoulomb vdwtype LJ
2144 % LocalWords: rvdw energygrp mdrun pre GMXLIB MOPAC GAMESS CPMD ONIOM
2145 % LocalWords: Morokuma iJ AQ AJ initio atomtype QMatom MMatom QMMM grps et al
2146 % LocalWords: QMmethod QMbasis QMMMscheme RHF DFT LYP CASSCF MMVB CASelectrons
2147 % LocalWords: CASorbitals planewavecutoff STO QMcharge QMmult diabatic edr
2148 % LocalWords: hyperline delocalized Coulombic indices nm ccc th ps
2149 % LocalWords: GTX CPUs GHz md sd bd vv avek tcoupl andersen tc OPLSAA GROMOS
2150 % LocalWords: OBC obc CCMA tol pbc xyz barostat pcoupl acc gpu PLUGIN Cmake GX
2151 % LocalWords: MSVC gcc installpath cmake DGMX DCMAKE functionalities GPGPU GTS
2152 % LocalWords: OpenCL DeviceID gromacs gpus html GeForce Quadro FX Plex CX GF
2153 % LocalWords: VMD DCD GROningen MAchine BIOSON Groningen der Spoel Drunen Comp
2154 % LocalWords: Phys Comm moltype intramol vdw Waals fep multivector pf
2155 % LocalWords: pymbar dhdl xvg LINCS Entropic entropic solutes ref com iso pm
2156 % LocalWords: rm prefactors equipotential potiso potisopf potpm trr
2157 % LocalWords: potrm potrmpf midplanes midplane gaussians potflex vars massw av
2158 % LocalWords: shure observables rccccccc vec eps dist min eqn transl nstsout
2159 % LocalWords: nstrout rotangles rotslabs rottorque RMSD rmsdfit excl NH amine
2160 % LocalWords: positionrestraint es SH phenylalanine solvated sh nanometer QM
2161 % LocalWords: Lennard Buckingham UK Svensson ab vsitetop co UHF MP interatomic
2162 % LocalWords: AdResS adress cg grp coords TF VOTCA thermoforce tf SPC userpot
2163 % LocalWords: itp sitesn atomtypes ff csg ndx Tesla CHARMM AA gauss
2164 % LocalWords: CMAP nocmap fators Monte performant lib LD DIR llllcc
2165 % LocalWords: CMake homepage DEV overclocking GT dlopen vmd VMDDIR
2166 % LocalWords: versa PME atomperatom graining forcefields hy spc OPLS
2167 % LocalWords: topol multi