4 This module can be used to enforce the rotation of a group of atoms, as
5 *e.g.* a protein subunit. There are a variety of rotation potentials,
6 among them complex ones that allow flexible adaptations of both the
7 rotated subunit as well as the local rotation axis during the
8 simulation. An example application can be found in ref.
9 :ref:`145 <refKutzner2011>`.
13 .. figure:: plots/rotation.*
16 Comparison of fixed and flexible axis rotation. A:
17 Rotating the sketched shape inside the white tubular cavity can
18 create artifacts when a fixed rotation axis (dashed) is used. More
19 realistically, the shape would revolve like a flexible pipe-cleaner
20 (dotted) inside the bearing (gray). B: Fixed rotation
21 around an axis :math:`\mathbf{v}` with a pivot point
22 specified by the vector :math:`\mathbf{u}`.
23 C: Subdividing the rotating fragment into slabs with
24 separate rotation axes (:math:`\uparrow`) and pivot points
25 (:math:`\bullet`) for each slab allows for flexibility. The distance
26 between two slabs with indices :math:`n` and :math:`n+1` is
29 .. _fig-equipotential:
31 .. figure:: plots/equipotential.*
34 Selection of different rotation potentials and definition of
35 notation. All four potentials :math:`V` (color coded) are shown for a
36 single atom at position :math:`\mathbf{x}_j(t)`.
37 A: Isotropic potential :math:`V^\mathrm{iso}`,
38 B: radial motion potential :math:`V^\mathrm{rm}` and
39 flexible potential :math:`V^\mathrm{flex}`, C–D: radial
40 motion2 potential :math:`V^\mathrm{rm2}` and flexible2 potential
41 :math:`V^\mathrm{flex2}` for :math:`\epsilon'\mathrm{ = }0\mathrm{ nm}^2`
42 (C) and :math:`\epsilon'\mathrm{ = }0.01\mathrm{nm}^2`
43 (D). The rotation axis is perpendicular to the plane
44 and marked by :math:`\otimes`. The light gray contours indicate
45 Boltzmann factors :math:`e^{-V/(k_B T)}` in the
46 :math:`\mathbf{x}_j`-plane for :math:`T=300`\ K and
47 :math:`k\mathrm{ = }200\mathrm{kJ}/(\mathrm{mol }\cdot\mathrm{nm}^2)`. The green
48 arrow shows the direction of the force
49 :math:`\mathbf{F}_{\!j}` acting on atom :math:`j`; the
50 blue dashed line indicates the motion of the reference position.
55 Stationary Axis with an Isotropic Potential
56 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
58 In the fixed axis approach (see :numref:`Fig. %s B <fig-rotation>`),
59 torque on a group of :math:`N` atoms with positions
60 :math:`\mathbf{x}_i` (denoted “rotation group”) is applied
61 by rotating a reference set of atomic positions – usually their initial
62 positions :math:`\mathbf{y}_i^0` – at a constant angular
63 velocity :math:`\omega` around an axis defined by a direction vector
64 :math:`\hat{\mathbf{v}}` and a pivot point
65 :math:`\mathbf{u}`. To that aim, each atom with
66 position :math:`\mathbf{x}_i` is attracted by a “virtual
67 spring” potential to its moving reference position
68 :math:`\mathbf{y}_i = \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{u})`,
69 where :math:`\mathbf{\Omega}(t)` is a matrix that describes the rotation
70 around the axis. In the simplest case, the “springs” are described by a
73 .. math:: V^\mathrm{iso} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[ \mathbf{\Omega}(t)
74 (\mathbf{y}_i^0 - \mathbf{u}) - (\mathbf{x}_i - \mathbf{u}) \right]^2
77 with optional mass-weighted prefactors :math:`w_i = N \, m_i/M` with
78 total mass :math:`M = \sum_{i=1}^N m_i`. The rotation matrix
79 :math:`\mathbf{\Omega}(t)` is
86 \cos\omega t + v_x^2{\,\xi\,}& v_x v_y{\,\xi\,}- v_z\sin\omega t & v_x v_z{\,\xi\,}+ v_y\sin\omega t\\
87 v_x v_y{\,\xi\,}+ v_z\sin\omega t & \cos\omega t + v_y^2{\,\xi\,}& v_y v_z{\,\xi\,}- v_x\sin\omega t\\
88 v_x v_z{\,\xi\,}- v_y\sin\omega t & v_y v_z{\,\xi\,}+ v_x\sin\omega t & \cos\omega t + v_z^2{\,\xi\,}\\
92 where :math:`v_x`, :math:`v_y`, and :math:`v_z` are the components of
93 the normalized rotation vector :math:`\hat{\mathbf{v}}`,
94 and :math:`{\,\xi\,}:= 1-\cos(\omega t)`. As illustrated in
95 :numref:`Fig. %s A <fig-equipotential>` for a single atom :math:`j`,
96 the rotation matrix :math:`\mathbf{\Omega}(t)` operates on the initial
98 :math:`\mathbf{y}_j^0 = \mathbf{x}_j(t_0)`
99 of atom :math:`j` at :math:`t=t_0`. At a later time :math:`t`, the
100 reference position has rotated away from its initial place (along the
101 blue dashed line), resulting in the force
103 .. math:: \mathbf{F}_{\!j}^\mathrm{iso}
104 = -\nabla_{\!j} \, V^\mathrm{iso}
106 \mathbf{\Omega}(t) (\mathbf{y}_j^0 - \mathbf{u}) - (\mathbf{x}_j - \mathbf{u} ) \right]
107 :label: eqnforcefixed
109 which is directed towards the reference position.
111 Pivot-Free Isotropic Potential
112 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
114 Instead of a fixed pivot vector :math:`\mathbf{u}` this
115 potential uses the center of mass :math:`\mathbf{x}_c` of
116 the rotation group as pivot for the rotation axis,
118 .. math:: \mathbf{x}_c = \frac{1}{M} \sum_{i=1}^N m_i \mathbf{x}_i
120 \mathbf{y}_c^0 = \frac{1}{M} \sum_{i=1}^N m_i \mathbf{y}_i^0 \ ,
123 which yields the “pivot-free” isotropic potential
125 .. math:: V^\mathrm{iso-pf} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[ \mathbf{\Omega}(t)
126 (\mathbf{y}_i^0 - \mathbf{y}_c^0) - (\mathbf{x}_i - \mathbf{x}_c) \right]^2 ,
131 .. math:: \mathbf{F}_{\!j}^\mathrm{iso-pf} = k \, w_j
133 \mathbf{\Omega}(t) ( \mathbf{y}_j^0 - \mathbf{y}_c^0)
134 - ( \mathbf{x}_j - \mathbf{x}_c )
136 :label: eqnforceisopf
138 Without mass-weighting, the pivot :math:`\mathbf{x}_c` is
139 the geometrical center of the group.
141 Parallel Motion Potential Variant
142 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
144 The forces generated by the isotropic potentials
145 (eqns. :eq:`%s <eqnpotiso>` and :eq:`%s <eqnpotisopf>`) also contain components parallel to the
146 rotation axis and thereby restrain motions along the axis of either the
147 whole rotation group (in case of :math:`V^\mathrm{iso}`) or within the
148 rotation group, in case of
150 .. math:: V^\mathrm{iso-pf}
153 unrestrained motion along the axis is preferred, we have implemented a
154 “parallel motion” variant by eliminating all components parallel to the
155 rotation axis for the potential. This is achieved by projecting the
156 distance vectors between reference and actual positions
158 .. math:: \mathbf{r}_i = \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{u}) - (\mathbf{x}_i - \mathbf{u})
160 onto the plane perpendicular to the rotation vector,
162 .. math:: \mathbf{r}_i^\perp := \mathbf{r}_i - (\mathbf{r}_i \cdot \hat{\mathbf{v}})\hat{\mathbf{v}}
167 .. math:: \begin{aligned}
169 V^\mathrm{pm} &=& \frac{k}{2} \sum_{i=1}^{N} w_i ( \mathbf{r}_i^\perp )^2 \\
170 &=& \frac{k}{2} \sum_{i=1}^{N} w_i
173 (\mathbf{y}_i^0 - \mathbf{u}) - (\mathbf{x}_i - \mathbf{u}) \right. \nonumber \\
174 && \left. - \left\lbrace
175 \left[ \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{u}) - (\mathbf{x}_i - \mathbf{u}) \right] \cdot\hat{\mathbf{v}}
176 \right\rbrace\hat{\mathbf{v}} \right\rbrace^2
182 .. math:: \mathbf{F}_{\!j}^\mathrm{pm} = k \, w_j \, \mathbf{r}_j^\perp
185 Pivot-Free Parallel Motion Potential
186 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
188 Replacing in eqn. :eq:`%s <eqnpotpm>` the fixed pivot
189 :math:`\mathbf{u}` by the center of mass
190 :math:`\mathbf{x_c}` yields the pivot-free variant of the
191 parallel motion potential. With
193 .. math:: \mathbf{s}_i = \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{y}_c^0) - (\mathbf{x}_i - \mathbf{x}_c)
195 the respective potential and forces are
197 .. math:: \begin{aligned}
198 V^\mathrm{pm-pf} &=& \frac{k}{2} \sum_{i=1}^{N} w_i ( \mathbf{s}_i^\perp )^2 \end{aligned}
201 .. math:: \begin{aligned}
202 \mathbf{F}_{\!j}^\mathrm{pm-pf} &=& k \, w_j \, \mathbf{s}_j^\perp
206 Radial Motion Potential
207 ^^^^^^^^^^^^^^^^^^^^^^^
209 In the above variants, the minimum of the rotation potential is either a
210 single point at the reference position
211 :math:`\mathbf{y}_i` (for the isotropic potentials) or a
212 single line through :math:`\mathbf{y}_i` parallel to the
213 rotation axis (for the parallel motion potentials). As a result, radial
214 forces restrict radial motions of the atoms. The two subsequent types of
215 rotation potentials, :math:`V^\mathrm{rm}` and :math:`V^\mathrm{rm2}`, drastically
216 reduce or even eliminate this effect. The first variant, :math:`V^\mathrm{rm}`
217 (:numref:`Fig. %s B <fig-equipotential>`), eliminates all force
218 components parallel to the vector connecting the reference atom and the
221 .. math:: V^\mathrm{rm} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[
223 \cdot(\mathbf{x}_i - \mathbf{u}) \right]^2 ,
231 \frac{\hat{\mathbf{v}}\times \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{u})} {\| \hat{\mathbf{v}}\times \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{u})\|} \ .
233 This variant depends only on the distance
234 :math:`\mathbf{p}_i \cdot (\mathbf{x}_i -
235 \mathbf{u})` of atom :math:`i` from the plane spanned by
236 :math:`\hat{\mathbf{v}}` and
237 :math:`\mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{u})`.
238 The resulting force is
240 .. math:: \mathbf{F}_{\!j}^\mathrm{rm} =
241 -k \, w_j \left[ \mathbf{p}_j\cdot(\mathbf{x}_j - \mathbf{u}) \right] \,\mathbf{p}_j \, .
242 :label: eqnpotrmforce
244 Pivot-Free Radial Motion Potential
245 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
247 Proceeding similar to the pivot-free isotropic potential yields a
248 pivot-free version of the above potential. With
253 \frac{\hat{\mathbf{v}}\times \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{y}_c^0)} {\| \hat{\mathbf{v}}\times \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{y}_c^0)\|} \, ,
255 the potential and force for the pivot-free variant of the radial motion
258 .. math:: \begin{aligned}
259 V^\mathrm{rm-pf} & = & \frac{k}{2} \sum_{i=1}^{N} w_i \left[
261 \cdot(\mathbf{x}_i - \mathbf{x}_c)
262 \right]^2 \, , \end{aligned}
265 .. math:: \begin{aligned}
266 \mathbf{F}_{\!j}^\mathrm{rm-pf} & = &
267 -k \, w_j \left[ \mathbf{q}_j\cdot(\mathbf{x}_j - \mathbf{x}_c) \right] \,\mathbf{q}_j
268 + k \frac{m_j}{M} \sum_{i=1}^{N} w_i \left[
269 \mathbf{q}_i\cdot(\mathbf{x}_i - \mathbf{x}_c) \right]\,\mathbf{q}_i \, .
271 :label: eqnpotrmpfforce
273 Radial Motion 2 Alternative Potential
274 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
276 As seen in :numref:`Fig. %s B <fig-equipotential>`, the force
277 resulting from :math:`V^\mathrm{rm}` still contains a small, second-order
278 radial component. In most cases, this perturbation is tolerable; if not,
279 the following alternative, :math:`V^\mathrm{rm2}`, fully eliminates the
280 radial contribution to the force, as depicted in
281 :numref:`Fig. %s C <fig-equipotential>`,
283 .. math:: V^\mathrm{rm2} =
284 \frac{k}{2} \sum_{i=1}^{N} w_i\,
285 \frac{\left[ (\hat{\mathbf{v}} \times ( \mathbf{x}_i - \mathbf{u} ))
286 \cdot \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{u}) \right]^2}
287 {\| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{u}) \|^2 +
291 where a small parameter :math:`\epsilon'` has been introduced to avoid
292 singularities. For :math:`\epsilon'\mathrm{ = }0\mathrm{nm}^2`, the
293 equipotential planes are spanned by :math:`\mathbf{x}_i -
294 \mathbf{u}` and :math:`\hat{\mathbf{v}}`,
295 yielding a force perpendicular to
296 :math:`\mathbf{x}_i - \mathbf{u}`, thus not
297 contracting or expanding structural parts that moved away from or toward
300 Choosing a small positive :math:`\epsilon'` (*e.g.*,
301 :math:`\epsilon'\mathrm{ = }0.01\mathrm{nm}^2`,
302 :numref:`Fig. %s D <fig-equipotential>`) in the denominator of
303 eqn. :eq:`%s <eqnpotrm2>` yields a well-defined potential and
304 continuous forces also close to the rotation axis, which is not the case
305 for :math:`\epsilon'\mathrm{ = }0\mathrm{nm}^2`
306 (:numref:`Fig. %s C <fig-equipotential>`). With
311 \mathbf{r}_i & := & \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{u})\\
312 \mathbf{s}_i & := & \frac{\hat{\mathbf{v}} \times (\mathbf{x}_i -
313 \mathbf{u} ) }{ \| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{u})
314 \| } \equiv \; \Psi_{i} \;\; {\hat{\mathbf{v}} \times
315 (\mathbf{x}_i-\mathbf{u} ) }\\
316 \Psi_i^{*} & := & \frac{1}{ \| \hat{\mathbf{v}} \times
317 (\mathbf{x}_i-\mathbf{u}) \|^2 + \epsilon'}\end{aligned}
319 the force on atom :math:`j` reads
321 .. math:: \mathbf{F}_{\!j}^\mathrm{rm2} =
324 (\mathbf{s}_j\cdot\mathbf{r}_{\!j})\;
325 \left[ \frac{\Psi_{\!j}^* }{\Psi_{\!j} } \mathbf{r}_{\!j}
326 - \frac{\Psi_{\!j}^{ * 2}}{\Psi_{\!j}^3}
327 (\mathbf{s}_j\cdot\mathbf{r}_{\!j})\mathbf{s}_j \right]
328 \right\rbrace \times \hat{\mathbf{v}} .
329 :label: eqnpotrm2force
331 Pivot-Free Radial Motion 2 Potential
332 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
334 The pivot-free variant of the above potential is
336 .. math:: V{^\mathrm{rm2-pf}}=
337 \frac{k}{2} \sum_{i=1}^{N} w_i\,
338 \frac{\left[ (\hat{\mathbf{v}} \times ( \mathbf{x}_i - \mathbf{x}_c ))
339 \cdot \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c) \right]^2}
340 {\| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c) \|^2 +
349 \mathbf{r}_i & := & \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c)\\
350 \mathbf{s}_i & := & \frac{\hat{\mathbf{v}} \times (\mathbf{x}_i -
351 \mathbf{x}_c ) }{ \| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c)
352 \| } \equiv \; \Psi_{i} \;\; {\hat{\mathbf{v}} \times
353 (\mathbf{x}_i-\mathbf{x}_c ) }\\ \Psi_i^{*} & := & \frac{1}{ \| \hat{\mathbf{v}} \times
354 (\mathbf{x}_i-\mathbf{x}_c) \|^2 + \epsilon'}\end{aligned}
356 the force on atom :math:`j` reads
358 .. math:: \begin{aligned}
360 \mathbf{F}_{\!j}{^\mathrm{rm2-pf}}& = &
363 (\mathbf{s}_j\cdot\mathbf{r}_{\!j})\;
364 \left[ \frac{\Psi_{\!j}^* }{\Psi_{\!j} } \mathbf{r}_{\!j}
365 - \frac{\Psi_{\!j}^{ * 2}}{\Psi_{\!j}^3}
366 (\mathbf{s}_j\cdot\mathbf{r}_{\!j})\mathbf{s}_j \right]
367 \right\rbrace \times \hat{\mathbf{v}}\\
369 + k\;\frac{m_j}{M} \left\lbrace \sum_{i=1}^{N}
370 w_i\;(\mathbf{s}_i\cdot\mathbf{r}_i) \;
371 \left[ \frac{\Psi_i^* }{\Psi_i } \mathbf{r}_i
372 - \frac{\Psi_i^{ * 2}}{\Psi_i^3} (\mathbf{s}_i\cdot\mathbf{r}_i )\;
373 \mathbf{s}_i \right] \right\rbrace \times \hat{\mathbf{v}} \, .
375 :label: eqnpotrm2pfforce
377 Flexible Axis Rotation
378 ~~~~~~~~~~~~~~~~~~~~~~
380 As sketched in :numref:`Fig. %s <fig-rotation>` A–B, the rigid body
381 behavior of the fixed axis rotation scheme is a drawback for many
382 applications. In particular, deformations of the rotation group are
383 suppressed when the equilibrium atom positions directly depend on the
384 reference positions. To avoid this limitation,
385 eqns. :eq:`%s <eqnpotrmpf>` and :eq:`%s <eqnpotrm2pf>`
386 will now be generalized towards a “flexible axis” as sketched in
387 :numref:`Fig. %s C <fig-rotation>`. This will be achieved by
388 subdividing the rotation group into a set of equidistant slabs
389 perpendicular to the rotation vector, and by applying a separate
390 rotation potential to each of these slabs.
391 :numref:`Fig. %s C <fig-rotation>` shows the midplanes of the slabs
392 as dotted straight lines and the centers as thick black dots.
394 To avoid discontinuities in the potential and in the forces, we define
395 “soft slabs” by weighing the contributions of each slab :math:`n` to the
396 total potential function :math:`V^\mathrm{flex}` by a Gaussian function
398 .. math:: g_n(\mathbf{x}_i) = \Gamma \ \mbox{exp} \left(
399 -\frac{\beta_n^2(\mathbf{x}_i)}{2\sigma^2} \right) ,
402 centered at the midplane of the :math:`n`\ th slab. Here :math:`\sigma`
403 is the width of the Gaussian function, :math:`\Delta x` the distance
404 between adjacent slabs, and
406 .. math:: \beta_n(\mathbf{x}_i) := \mathbf{x}_i \cdot \hat{\mathbf{v}} - n \, \Delta x \, .
410 .. figure:: plots/gaussians.*
413 Gaussian functions :math:`g_n` centered at :math:`n \, \Delta x` for
414 a slab distance :math:`\Delta x = 1.5` nm and :math:`n \geq -2`.
415 Gaussian function :math:`g_0` is highlighted in bold; the dashed line
416 depicts the sum of the shown Gaussian functions.
418 A most convenient choice is :math:`\sigma = 0.7 \Delta x` and
422 1/\Gamma = \sum_{n \in Z}
424 \left(-\frac{(n - \frac{1}{4})^2}{2\cdot 0.7^2}\right)
427 which yields a nearly constant sum, essentially independent of
428 :math:`\mathbf{x}_i` (dashed line in
429 :numref:`Fig. %s <fig-gaussian>`), *i.e.*,
431 .. math:: \sum_{n \in Z} g_n(\mathbf{x}_i) = 1 + \epsilon(\mathbf{x}_i) \, ,
435 :math:`| \epsilon(\mathbf{x}_i) | < 1.3\cdot 10^{-4}`.
436 This choice also implies that the individual contributions to the force
437 from the slabs add up to unity such that no further normalization is
440 To each slab center :math:`\mathbf{x}_c^n`, all atoms
441 contribute by their Gaussian-weighted (optionally also mass-weighted)
443 :math:`g_n(\mathbf{x}_i) \, \mathbf{x}_i`.
444 The instantaneous slab centers :math:`\mathbf{x}_c^n` are
445 calculated from the current positions
446 :math:`\mathbf{x}_i`,
448 .. math:: \mathbf{x}_c^n =
449 \frac{\sum_{i=1}^N g_n(\mathbf{x}_i) \, m_i \, \mathbf{x}_i}
450 {\sum_{i=1}^N g_n(\mathbf{x}_i) \, m_i} \, ,\\
453 while the reference centers :math:`\mathbf{y}_c^n` are
454 calculated from the reference positions
455 :math:`\mathbf{y}_i^0`,
457 .. math:: \mathbf{y}_c^n =
458 \frac{\sum_{i=1}^N g_n(\mathbf{y}_i^0) \, m_i \, \mathbf{y}_i^0}
459 {\sum_{i=1}^N g_n(\mathbf{y}_i^0) \, m_i} \, .
462 Due to the rapid decay of :math:`g_n`, each slab will essentially
463 involve contributions from atoms located within :math:`\approx
464 3\Delta x` from the slab center only.
466 Flexible Axis Potential
467 ^^^^^^^^^^^^^^^^^^^^^^^
469 We consider two flexible axis variants. For the first variant, the slab
470 segmentation procedure with Gaussian weighting is applied to the radial
472 (eqn. :eq:`%s <eqnpotrmpf>` / :numref:`Fig. %s B <fig-equipotential>`),
473 yielding as the contribution of slab :math:`n`
476 \frac{k}{2} \sum_{i=1}^{N} w_i \, g_n(\mathbf{x}_i)
480 (\mathbf{x}_i - \mathbf{x}_c^n)
484 and a total potential function
486 .. math:: V^\mathrm{flex} = \sum_n V^n \, .
489 Note that the global center of mass :math:`\mathbf{x}_c`
490 used in eqn. :eq:`%s <eqnpotrmpf>` is now replaced by
491 :math:`\mathbf{x}_c^n`, the center of mass of the slab.
497 \mathbf{q}_i^n & := & \frac{\hat{\mathbf{v}} \times
498 \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) }{ \| \hat{\mathbf{v}}
499 \times \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) \| } \\
500 b_i^n & := & \mathbf{q}_i^n \cdot (\mathbf{x}_i - \mathbf{x}_c^n) \, ,\end{aligned}
502 the resulting force on atom :math:`j` reads
504 .. math:: \begin{aligned}
505 \nonumber\hspace{-15mm}
506 \mathbf{F}_{\!j}^\mathrm{flex} &=&
507 - \, k \, w_j \sum_n g_n(\mathbf{x}_j) \, b_j^n \left\lbrace \mathbf{q}_j^n -
508 b_j^n \frac{\beta_n(\mathbf{x}_j)}{2\sigma^2} \hat{\mathbf{v}} \right\rbrace \\ & &
509 + \, k \, m_j \sum_n \frac{g_n(\mathbf{x}_j)}{\sum_h g_n(\mathbf{x}_h)}
510 \sum_{i=1}^{N} w_i \, g_n(\mathbf{x}_i) \, b_i^n \left\lbrace
511 \mathbf{q}_i^n -\frac{\beta_n(\mathbf{x}_j)}{\sigma^2}
512 \left[ \mathbf{q}_i^n \cdot (\mathbf{x}_j - \mathbf{x}_c^n )\right]
513 \hat{\mathbf{v}} \right\rbrace .
515 :label: eqnpotflexforce
517 Note that for :math:`V^\mathrm{flex}`, as defined, the slabs are fixed in
518 space and so are the reference centers
519 :math:`\mathbf{y}_c^n`. If during the simulation the
520 rotation group moves too far in :math:`\mathbf{v}`
521 direction, it may enter a region where – due to the lack of nearby
522 reference positions – no reference slab centers are defined, rendering
523 the potential evaluation impossible. We therefore have included a
524 slightly modified version of this potential that avoids this problem by
525 attaching the midplane of slab :math:`n=0` to the center of mass of the
526 rotation group, yielding slabs that move with the rotation group. This
527 is achieved by subtracting the center of mass
528 :math:`\mathbf{x}_c` of the group from the positions,
530 .. math:: \tilde{\mathbf{x}}_i = \mathbf{x}_i - \mathbf{x}_c \, , \mbox{\ \ \ and \ \ }
531 \tilde{\mathbf{y}}_i^0 = \mathbf{y}_i^0 - \mathbf{y}_c^0 \, ,
536 .. math:: \begin{aligned}
538 & = & \frac{k}{2} \sum_n \sum_{i=1}^{N} w_i \, g_n(\tilde{\mathbf{x}}_i)
539 \left[ \frac{\hat{\mathbf{v}} \times \mathbf{\Omega}(t)(\tilde{\mathbf{y}}_i^0
540 - \tilde{\mathbf{y}}_c^n) }{ \| \hat{\mathbf{v}} \times
541 \mathbf{\Omega}(t)(\tilde{\mathbf{y}}_i^0 -
542 \tilde{\mathbf{y}}_c^n) \| }
544 (\tilde{\mathbf{x}}_i - \tilde{\mathbf{x}}_c^n)
549 To simplify the force derivation, and for efficiency reasons, we here
550 assume :math:`\mathbf{x}_c` to be constant, and thus
551 :math:`\partial \mathbf{x}_c / \partial x =
552 \partial \mathbf{x}_c / \partial y = \partial \mathbf{x}_c / \partial z = 0`.
553 The resulting force error is small (of order :math:`O(1/N)` or
554 :math:`O(m_j/M)` if mass-weighting is applied) and can therefore be
555 tolerated. With this assumption, the forces
558 \mathbf{F}^\mathrm{flex-t}
560 have the same form as
561 eqn. :eq:`%s <eqnpotflexforce>`.
563 Flexible Axis 2 Alternative Potential
564 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
566 In this second variant, slab segmentation is applied to
567 :math:`V^\mathrm{rm2}` (eqn. :eq:`%s <eqnpotrm2pf>`), resulting in
568 a flexible axis potential without radial force contributions
569 (:numref:`Fig. %s C <fig-equipotential>`),
571 .. math:: V{^\mathrm{flex2}}=
572 \frac{k}{2} \sum_{i=1}^{N} \sum_n w_i\,g_n(\mathbf{x}_i)
573 \frac{\left[ (\hat{\mathbf{v}} \times ( \mathbf{x}_i - \mathbf{x}_c^n ))
574 \cdot \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) \right]^2}
575 {\| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c^n) \|^2 +
581 .. math:: \begin{aligned}
582 \mathbf{r}_i^n & := & \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n)\\
583 \mathbf{s}_i^n & := & \frac{\hat{\mathbf{v}} \times (\mathbf{x}_i -
584 \mathbf{x}_c^n ) }{ \| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c^n)
585 \| } \equiv \; \psi_{i} \;\; {\hat{\mathbf{v}} \times (\mathbf{x}_i-\mathbf{x}_c^n ) }\\
586 \psi_i^{*} & := & \frac{1}{ \| \hat{\mathbf{v}} \times (\mathbf{x}_i-\mathbf{x}_c^n) \|^2 + \epsilon'}\\
587 W_j^n & := & \frac{g_n(\mathbf{x}_j)\,m_j}{\sum_h g_n(\mathbf{x}_h)\,m_h}\\
589 \sum_{i=1}^{N} w_i\;g_n(\mathbf{x}_i)
590 \; (\mathbf{s}_i^n\cdot\mathbf{r}_i^n)
591 \left[ \frac{\psi_i^* }{\psi_i } \mathbf{r}_i^n
592 - \frac{\psi_i^{ * 2}}{\psi_i^3} (\mathbf{s}_i^n\cdot\mathbf{r}_i^n )\;
593 \mathbf{s}_i^n \right]
597 the force on atom :math:`j` reads
599 .. math:: \begin{aligned}
601 \mathbf{F}_{\!j}{^\mathrm{flex2}}& = &
603 \left\lbrace \sum_n w_j\;g_n(\mathbf{x}_j)\;
604 (\mathbf{s}_j^n\cdot\mathbf{r}_{\!j}^n)\;
605 \left[ \frac{\psi_j^* }{\psi_j } \mathbf{r}_{\!j}^n
606 - \frac{\psi_j^{ * 2}}{\psi_j^3} (\mathbf{s}_j^n\cdot\mathbf{r}_{\!j}^n)\;
607 \mathbf{s}_{\!j}^n \right] \right\rbrace \times \hat{\mathbf{v}} \\
610 + k \left\lbrace \sum_n W_{\!j}^n \, \mathbf{S}^n \right\rbrace \times
612 - k \left\lbrace \sum_n W_{\!j}^n \; \frac{\beta_n(\mathbf{x}_j)}{\sigma^2} \frac{1}{\psi_j}\;\;
614 \mathbf{S}^n \right\rbrace \hat{\mathbf{v}}\\
616 + \frac{k}{2} \left\lbrace \sum_n w_j\;g_n(\mathbf{x}_j)
617 \frac{\beta_n(\mathbf{x}_j)}{\sigma^2}
618 \frac{\psi_j^*}{\psi_j^2}( \mathbf{s}_j^n \cdot \mathbf{r}_{\!j}^n )^2 \right\rbrace
621 :label: eqnpotflex2force
623 Applying transformation :eq:`%s <eqntrafo>` yields a
624 “translation-tolerant” version of the flexible2 potential,
625 :math:`V{^\mathrm{flex2 - t}}`. Again, assuming that
626 :math:`\partial \mathbf{x}_c / \partial x`,
627 :math:`\partial \mathbf{x}_c /
628 \partial y`, :math:`\partial \mathbf{x}_c / \partial z`
629 are small, the resulting equations for :math:`V{^\mathrm{flex2 - t}}`
630 and :math:`\mathbf{F}{^\mathrm{flex2 - t}}` are similar
631 to those of :math:`V^\mathrm{flex2}` and
632 :math:`\mathbf{F}^\mathrm{flex2}`.
637 To apply enforced rotation, the particles :math:`i` that are to be
638 subjected to one of the rotation potentials are defined via index groups
639 ``rot-group0``, ``rot-group1``, etc., in the
640 :ref:`mdp` input file. The reference positions
641 :math:`\mathbf{y}_i^0` are read from a special
642 :ref:`trr` file provided to :ref:`grompp <gmx grompp>`. If no such
643 file is found, :math:`\mathbf{x}_i(t=0)` are used as
644 reference positions and written to :ref:`trr` such that they
645 can be used for subsequent setups. All parameters of the potentials such
646 as :math:`k`, :math:`\epsilon'`, etc.
647 (:numref:`Table %s <tab-vars>`) are provided as :ref:`mdp`
648 parameters; ``rot-type`` selects the type of the potential.
649 The option ``rot-massw`` allows to choose whether or not to
650 use mass-weighted averaging. For the flexible potentials, a cutoff value
651 :math:`g_n^\mathrm{min}` (typically :math:`g_n^\mathrm{min}=0.001`)
652 makes sure that only significant contributions to :math:`V` and
653 :math:`\mathbf{F}` are evaluated, *i.e.* terms with
654 :math:`g_n(\mathbf{x}) < g_n^\mathrm{min}` are omitted.
655 :numref:`Table %s <tab-quantities>` summarizes observables that are
656 written to additional output files and which are described below.
658 .. |ROTISO| replace:: V\ :math:`^\mathrm{iso}`
659 .. |ROTISOPF| replace:: V\ :math:`^\mathrm{iso-pf}`
660 .. |ROTPM| replace:: V\ :math:`^\mathrm{pm}`
661 .. |ROTPMPF| replace:: V\ :math:`^\mathrm{pm-pf}`
662 .. |ROTRM| replace:: V\ :math:`^\mathrm{rm}`
663 .. |ROTRMPF| replace:: V\ :math:`^\mathrm{rm-pf}`
664 .. |ROTRM2| replace:: V\ :math:`^\mathrm{rm2}`
665 .. |ROTRM2PF| replace:: V\ :math:`^\mathrm{rm2-pf}`
666 .. |ROTFL| replace:: V\ :math:`^\mathrm{flex}`
667 .. |ROTFLT| replace:: V\ :math:`^\mathrm{flex-t}`
668 .. |ROTFL2| replace:: V\ :math:`^\mathrm{flex2}`
669 .. |ROTFLT2| replace:: V\ :math:`^\mathrm{flex2-t}`
670 .. |KUNIT| replace:: :math:`\frac{\mathrm{kJ}}{\mathrm{mol} \cdot \mathrm{nm}^2}`
671 .. |BFX| replace:: **x**
672 .. |KMA| replace:: :math:`k`
673 .. |VECV| replace:: :math:`\hat{\mathbf{v}}`
674 .. |VECU| replace:: :math:`\mathbf{u}`
675 .. |OMEG| replace:: :math:`\omega`
676 .. |EPSP| replace:: :math:`{\epsilon}'`
677 .. |DELX| replace:: :math:`{\Delta}x`
678 .. |GMIN| replace:: :math:`g_n^\mathrm{min}`
679 .. |CIPS| replace:: :math:`^\circ`\ /ps
680 .. |NM2| replace:: nm\ :math:`^2`
681 .. |REF1| replace:: \ :eq:`eqnpotiso`
682 .. |REF2| replace:: \ :eq:`eqnpotisopf`
683 .. |REF3| replace:: \ :eq:`eqnpotpm`
684 .. |REF4| replace:: \ :eq:`eqnpotpmpf`
685 .. |REF5| replace:: \ :eq:`eqnpotrm`
686 .. |REF6| replace:: \ :eq:`eqnpotrmpf`
687 .. |REF7| replace:: \ :eq:`eqnpotrm2`
688 .. |REF8| replace:: \ :eq:`eqnpotrm2pf`
689 .. |REF9| replace:: \ :eq:`eqnpotflex`
690 .. |REF10| replace:: \ :eq:`eqnpotflext`
691 .. |REF11| replace:: \ :eq:`eqnpotflex2`
695 .. table:: Parameters used by the various rotation potentials.
696 |BFX| indicate which parameter is actually used for a given potential
700 +------------------------------------------+---------+--------+--------+--------+--------+-----------+-----------+
701 | parameter | |KMA| | |VECV| | |VECU| | |OMEG| | |EPSP| | |DELX| | |GMIN| |
702 +------------------------------------------+---------+--------+--------+--------+--------+-----------+-----------+
703 | :ref:`mdp` input variable name | k | vec | pivot | rate | eps | slab-dist | min-gauss |
704 +------------------------------------------+---------+--------+--------+--------+--------+-----------+-----------+
705 | unit | |KUNIT| | ``-`` | nm | |CIPS| | |NM2| | nm | ``-`` |
706 +================================+=========+=========+========+========+========+========+===========+===========+
707 | fixed axis potentials: | eqn. |
708 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
709 | isotropic | |ROTISO| | |REF1| | |BFX| | |BFX| | |BFX| | |BFX| | ``-`` | ``-`` | ``-`` |
710 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
711 | --- pivot-free | |ROTISOPF| | |REF2| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | ``-`` | ``-`` |
712 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
713 | parallel motion | |ROTPM| | |REF3| | |BFX| | |BFX| | |BFX| | |BFX| | ``-`` | ``-`` | ``-`` |
714 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
715 | --- pivot-free | |ROTPMPF| | |REF4| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | ``-`` | ``-`` |
716 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
717 | radial motion | |ROTRM| | |REF5| | |BFX| | |BFX| | |BFX| | |BFX| | ``-`` | ``-`` | ``-`` |
718 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
719 | --- pivot-free | |ROTRMPF| | |REF6| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | ``-`` | ``-`` |
720 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
721 | radial motion 2 | |ROTRM2| | |REF7| | |BFX| | |BFX| | |BFX| | |BFX| | |BFX| | ``-`` | ``-`` |
722 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
723 | --- pivot-free | |ROTRM2PF| | |REF8| | |BFX| | |BFX| | ``-`` | |BFX| | |BFX| | ``-`` | ``-`` |
724 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
725 | flexible axis potentials: | eqn. |
726 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
727 | flexible | |ROTFL| | |REF9| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | |BFX| | |BFX| |
728 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
729 | --- transl. tol | |ROTFLT| | |REF10| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | |BFX| | |BFX| |
730 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
731 | flexible 2 | |ROTFL2| | |REF11| | |BFX| | |BFX| | ``-`` | |BFX| | |BFX| | |BFX| | |BFX| |
732 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
733 | --- transl. tol | |ROTFLT2| | ``-`` | |BFX| | |BFX| | ``-`` | |BFX| | |BFX| | |BFX| | |BFX| |
734 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
739 .. |VT| replace:: :math:`V(t)`
740 .. |THET| replace:: :math:`\theta_{ \mathrm{ref}}(t)`
741 .. |THETAV| replace:: :math:`\theta_{\mathrm{av}}(t)`
742 .. |THETFIT| replace:: :math:`\theta_{\mathrm{fit}}(t)`, :math:`\theta_{\mathrm{fit}}(t,n)`
743 .. |YVEC| replace:: :math:`\mathbf{y}_{0}(n)`, :math:`\mathbf{x}_{0}(t,n)`
744 .. |TAUT| replace:: :math:`\tau(t)`
745 .. |TAUTN| replace:: :math:`\tau(t,n)`
746 .. |REFT| replace:: :numref:`see Table %s <tab-vars>`
747 .. |REFEQ| replace:: :math:`\theta_{\mathrm{ref}}(t)=\omega t`
748 .. |REF12| replace:: \ :eq:`eqnavangle`
749 .. |REF13| replace:: \ :eq:`eqnrmsdfit`
750 .. |REF14| replace:: \ :eq:`eqndefx0`\ ,\ :eq:`eqndefy0`
751 .. |REF15| replace:: \ :eq:`eqntorque`
755 .. table:: Quantities recorded in output files during enforced rotation simulations.
756 All slab-wise data is written every ``nstsout`` steps, other rotation data every ``nstrout`` steps.
760 +------------+---------+------------+--------------------+-------+----------+
761 | quantity | unit | equation | output file | fixed | flexible |
762 +============+=========+============+====================+=======+==========+
763 | |VT| | kJ/mol | |REFT| | ``rotation`` | |BFX| | |BFX| |
764 +------------+---------+------------+--------------------+-------+----------+
765 | |THET| | degrees | |REFEQ| | ``rotation`` | |BFX| | |BFX| |
766 +------------+---------+------------+--------------------+-------+----------+
767 | |THETAV| | degrees | |REF12| | ``rotation`` | |BFX| | ``-`` |
768 +------------+---------+------------+--------------------+-------+----------+
769 | |THETFIT| | degrees | |REF13| | ``rotangles`` | ``-`` | |BFX| |
770 +------------+---------+------------+--------------------+-------+----------+
771 | |YVEC| | nm | |REF14| | ``rotslabs`` | ``-`` | |BFX| |
772 +------------+---------+------------+--------------------+-------+----------+
773 | |TAUT| | kJ/mol | |REF15| | ``rotation`` | |BFX| | ``-`` |
774 +------------+---------+------------+--------------------+-------+----------+
775 | |TAUTN| | kJ/mol | |REF15| | ``rottorque`` | ``-`` | |BFX| |
776 +------------+---------+------------+--------------------+-------+----------+
781 Angle of Rotation Groups: Fixed Axis
782 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
784 For fixed axis rotation, the average angle :math:`\theta_\mathrm{av}(t)`
785 of the group relative to the reference group is determined via the
786 distance-weighted angular deviation of all rotation group atoms from
787 their reference positions,
789 .. math:: \theta_\mathrm{av} = \left. \sum_{i=1}^{N} r_i \ \theta_i \right/ \sum_{i=1}^N r_i \ .
792 Here, :math:`r_i` is the distance of the reference position to the
793 rotation axis, and the difference angles :math:`\theta_i` are determined
794 from the atomic positions, projected onto a plane perpendicular to the
795 rotation axis through pivot point :math:`\mathbf{u}` (see
796 eqn. :eq:`%s <eqnproject>` for the definition of
802 \frac{(\mathbf{y}_i-\mathbf{u})^\perp \cdot (\mathbf{x}_i-\mathbf{u})^\perp}
803 { \| (\mathbf{y}_i-\mathbf{u})^\perp \cdot (\mathbf{x}_i-\mathbf{u})^\perp
806 The sign of :math:`\theta_\mathrm{av}` is chosen such that
807 :math:`\theta_\mathrm{av} > 0` if the actual structure rotates ahead of
810 Angle of Rotation Groups: Flexible Axis
811 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
813 For flexible axis rotation, two outputs are provided, the angle of the
814 entire rotation group, and separate angles for the segments in the
815 slabs. The angle of the entire rotation group is determined by an RMSD
816 fit of :math:`\mathbf{x}_i` to the reference positions
817 :math:`\mathbf{y}_i^0` at :math:`t=0`, yielding
818 :math:`\theta_\mathrm{fit}` as the angle by which the reference has to
819 be rotated around :math:`\hat{\mathbf{v}}` for the optimal
822 .. math:: \mathrm{RMSD} \big( \mathbf{x}_i,\ \mathbf{\Omega}(\theta_\mathrm{fit})
823 \mathbf{y}_i^0 \big) \stackrel{!}{=} \mathrm{min} \, .
826 To determine the local angle for each slab :math:`n`, both reference
827 and actual positions are weighted with the Gaussian function of slab
828 :math:`n`, and :math:`\theta_\mathrm{fit}(t,n)` is calculated as in
829 eqn. :eq:`%s <eqnrmsdfit>` from the Gaussian-weighted
832 For all angles, the :ref:`mdp` input option
833 ``rot-fit-method`` controls whether a normal RMSD fit is
834 performed or whether for the fit each position
835 :math:`\mathbf{x}_i` is put at the same distance to the
836 rotation axis as its reference counterpart
837 :math:`\mathbf{y}_i^0`. In the latter case, the RMSD
838 measures only angular differences, not radial ones.
840 Angle Determination by Searching the Energy Minimum
841 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
843 Alternatively, for ``rot-fit-method = potential``, the angle
844 of the rotation group is determined as the angle for which the rotation
845 potential energy is minimal. Therefore, the used rotation potential is
846 additionally evaluated for a set of angles around the current reference
847 angle. In this case, the ``rotangles.log`` output file
848 contains the values of the rotation potential at the chosen set of
849 angles, while ``rotation.xvg`` lists the angle with minimal
855 The torque :math:`\mathbf{\tau}(t)` exerted by the
856 rotation potential is calculated for fixed axis rotation via
858 .. math:: \mathbf{\tau}(t) = \sum_{i=1}^{N} \mathbf{r}_i(t) \times \mathbf{f}_{\!i}^\perp(t) ,
861 where :math:`\mathbf{r}_i(t)` is the distance vector from
862 the rotation axis to :math:`\mathbf{x}_i(t)` and
863 :math:`\mathbf{f}_{\!i}^\perp(t)` is the force component
864 perpendicular to :math:`\mathbf{r}_i(t)` and
865 :math:`\hat{\mathbf{v}}`. For flexible axis rotation,
866 torques :math:`\mathbf{\tau}_{\!n}` are calculated for
867 each slab using the local rotation axis of the slab and the
868 Gaussian-weighted positions.