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
81 .. math:: \mathbf{\Omega}(t) =
84 \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\\
85 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\\
86 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\,}\\
91 where :math:`v_x`, :math:`v_y`, and :math:`v_z` are the components of
92 the normalized rotation vector :math:`\hat{\mathbf{v}}`,
93 and :math:`{\,\xi\,}:= 1-\cos(\omega t)`. As illustrated in
94 :numref:`Fig. %s A <fig-equipotential>` for a single atom :math:`j`,
95 the rotation matrix :math:`\mathbf{\Omega}(t)` operates on the initial
97 :math:`\mathbf{y}_j^0 = \mathbf{x}_j(t_0)`
98 of atom :math:`j` at :math:`t=t_0`. At a later time :math:`t`, the
99 reference position has rotated away from its initial place (along the
100 blue dashed line), resulting in the force
102 .. math:: \mathbf{F}_{\!j}^\mathrm{iso}
103 = -\nabla_{\!j} \, V^\mathrm{iso}
105 \mathbf{\Omega}(t) (\mathbf{y}_j^0 - \mathbf{u}) - (\mathbf{x}_j - \mathbf{u} ) \right]
106 :label: eqnforcefixed
108 which is directed towards the reference position.
110 Pivot-Free Isotropic Potential
111 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
113 Instead of a fixed pivot vector :math:`\mathbf{u}` this
114 potential uses the center of mass :math:`\mathbf{x}_c` of
115 the rotation group as pivot for the rotation axis,
117 .. math:: \mathbf{x}_c = \frac{1}{M} \sum_{i=1}^N m_i \mathbf{x}_i
119 \mathbf{y}_c^0 = \frac{1}{M} \sum_{i=1}^N m_i \mathbf{y}_i^0 \ ,
122 which yields the “pivot-free” isotropic potential
124 .. math:: V^\mathrm{iso-pf} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[ \mathbf{\Omega}(t)
125 (\mathbf{y}_i^0 - \mathbf{y}_c^0) - (\mathbf{x}_i - \mathbf{x}_c) \right]^2 ,
130 .. math:: \mathbf{F}_{\!j}^\mathrm{iso-pf} = k \, w_j
132 \mathbf{\Omega}(t) ( \mathbf{y}_j^0 - \mathbf{y}_c^0)
133 - ( \mathbf{x}_j - \mathbf{x}_c )
135 :label: eqnforceisopf
137 Without mass-weighting, the pivot :math:`\mathbf{x}_c` is
138 the geometrical center of the group.
140 Parallel Motion Potential Variant
141 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
143 The forces generated by the isotropic potentials
144 (eqns. :eq:`%s <eqnpotiso>` and :eq:`%s <eqnpotisopf>`) also contain components parallel to the
145 rotation axis and thereby restrain motions along the axis of either the
146 whole rotation group (in case of :math:`V^\mathrm{iso}`) or within the
147 rotation group, in case of :math:`V^\mathrm{iso-pf}`.
149 For cases where unrestrained motion along the axis is preferred, we have implemented a
150 “parallel motion” variant by eliminating all components parallel to the
151 rotation axis for the potential. This is achieved by projecting the
152 distance vectors between reference and actual positions
154 .. math:: \mathbf{r}_i = \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{u}) - (\mathbf{x}_i - \mathbf{u})
155 :label: eqnrotdistvectors
157 onto the plane perpendicular to the rotation vector,
159 .. math:: \mathbf{r}_i^\perp := \mathbf{r}_i - (\mathbf{r}_i \cdot \hat{\mathbf{v}})\hat{\mathbf{v}}
164 .. math:: \begin{aligned}
166 V^\mathrm{pm} &=& \frac{k}{2} \sum_{i=1}^{N} w_i ( \mathbf{r}_i^\perp )^2 \\
167 &=& \frac{k}{2} \sum_{i=1}^{N} w_i
170 (\mathbf{y}_i^0 - \mathbf{u}) - (\mathbf{x}_i - \mathbf{u}) \right. \nonumber \\
171 && \left. - \left\lbrace
172 \left[ \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{u}) - (\mathbf{x}_i - \mathbf{u}) \right] \cdot\hat{\mathbf{v}}
173 \right\rbrace\hat{\mathbf{v}} \right\rbrace^2
179 .. math:: \mathbf{F}_{\!j}^\mathrm{pm} = k \, w_j \, \mathbf{r}_j^\perp
182 Pivot-Free Parallel Motion Potential
183 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
185 Replacing in eqn. :eq:`%s <eqnpotpm>` the fixed pivot
186 :math:`\mathbf{u}` by the center of mass
187 :math:`\mathbf{x_c}` yields the pivot-free variant of the
188 parallel motion potential. With
190 .. math:: \mathbf{s}_i = \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{y}_c^0) - (\mathbf{x}_i - \mathbf{x}_c)
191 :label: eqnparrallelpotential
193 the respective potential and forces are
195 .. math:: \begin{aligned}
196 V^\mathrm{pm-pf} &=& \frac{k}{2} \sum_{i=1}^{N} w_i ( \mathbf{s}_i^\perp )^2 \end{aligned}
199 .. math:: \begin{aligned}
200 \mathbf{F}_{\!j}^\mathrm{pm-pf} &=& k \, w_j \, \mathbf{s}_j^\perp
204 Radial Motion Potential
205 ^^^^^^^^^^^^^^^^^^^^^^^
207 In the above variants, the minimum of the rotation potential is either a
208 single point at the reference position
209 :math:`\mathbf{y}_i` (for the isotropic potentials) or a
210 single line through :math:`\mathbf{y}_i` parallel to the
211 rotation axis (for the parallel motion potentials). As a result, radial
212 forces restrict radial motions of the atoms. The two subsequent types of
213 rotation potentials, :math:`V^\mathrm{rm}` and :math:`V^\mathrm{rm2}`, drastically
214 reduce or even eliminate this effect. The first variant, :math:`V^\mathrm{rm}`
215 (:numref:`Fig. %s B <fig-equipotential>`), eliminates all force
216 components parallel to the vector connecting the reference atom and the
219 .. math:: V^\mathrm{rm} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[
221 \cdot(\mathbf{x}_i - \mathbf{u}) \right]^2 ,
226 .. math:: \mathbf{p}_i :=
227 \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})\|} \ .
228 :label: eqnpotrmpart2
230 This variant depends only on the distance
231 :math:`\mathbf{p}_i \cdot (\mathbf{x}_i -
232 \mathbf{u})` of atom :math:`i` from the plane spanned by
233 :math:`\hat{\mathbf{v}}` and
234 :math:`\mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{u})`.
235 The resulting force is
237 .. math:: \mathbf{F}_{\!j}^\mathrm{rm} =
238 -k \, w_j \left[ \mathbf{p}_j\cdot(\mathbf{x}_j - \mathbf{u}) \right] \,\mathbf{p}_j \, .
239 :label: eqnpotrmforce
241 Pivot-Free Radial Motion Potential
242 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
244 Proceeding similar to the pivot-free isotropic potential yields a
245 pivot-free version of the above potential. With
247 .. math:: \mathbf{q}_i :=
248 \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)\|} \, ,
249 :label: eqnpotrmpfpart1
251 the potential and force for the pivot-free variant of the radial motion
254 .. math:: \begin{aligned}
255 V^\mathrm{rm-pf} & = & \frac{k}{2} \sum_{i=1}^{N} w_i \left[
257 \cdot(\mathbf{x}_i - \mathbf{x}_c)
258 \right]^2 \, , \end{aligned}
261 .. math:: \begin{aligned}
262 \mathbf{F}_{\!j}^\mathrm{rm-pf} & = &
263 -k \, w_j \left[ \mathbf{q}_j\cdot(\mathbf{x}_j - \mathbf{x}_c) \right] \,\mathbf{q}_j
264 + k \frac{m_j}{M} \sum_{i=1}^{N} w_i \left[
265 \mathbf{q}_i\cdot(\mathbf{x}_i - \mathbf{x}_c) \right]\,\mathbf{q}_i \, .
267 :label: eqnpotrmpfforce
269 Radial Motion 2 Alternative Potential
270 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
272 As seen in :numref:`Fig. %s B <fig-equipotential>`, the force
273 resulting from :math:`V^\mathrm{rm}` still contains a small, second-order
274 radial component. In most cases, this perturbation is tolerable; if not,
275 the following alternative, :math:`V^\mathrm{rm2}`, fully eliminates the
276 radial contribution to the force, as depicted in
277 :numref:`Fig. %s C <fig-equipotential>`,
279 .. math:: V^\mathrm{rm2} =
280 \frac{k}{2} \sum_{i=1}^{N} w_i\,
281 \frac{\left[ (\hat{\mathbf{v}} \times ( \mathbf{x}_i - \mathbf{u} ))
282 \cdot \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{u}) \right]^2}
283 {\| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{u}) \|^2 +
287 where a small parameter :math:`\epsilon'` has been introduced to avoid
288 singularities. For :math:`\epsilon'\mathrm{ = }0\mathrm{nm}^2`, the
289 equipotential planes are spanned by :math:`\mathbf{x}_i -
290 \mathbf{u}` and :math:`\hat{\mathbf{v}}`,
291 yielding a force perpendicular to
292 :math:`\mathbf{x}_i - \mathbf{u}`, thus not
293 contracting or expanding structural parts that moved away from or toward
296 Choosing a small positive :math:`\epsilon'` (*e.g.*,
297 :math:`\epsilon'\mathrm{ = }0.01\mathrm{nm}^2`,
298 :numref:`Fig. %s D <fig-equipotential>`) in the denominator of
299 eqn. :eq:`%s <eqnpotrm2>` yields a well-defined potential and
300 continuous forces also close to the rotation axis, which is not the case
301 for :math:`\epsilon'\mathrm{ = }0\mathrm{nm}^2`
302 (:numref:`Fig. %s C <fig-equipotential>`). With
304 .. math:: \begin{aligned}
305 \mathbf{r}_i & := & \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{u})\\
306 \mathbf{s}_i & := & \frac{\hat{\mathbf{v}} \times (\mathbf{x}_i -
307 \mathbf{u} ) }{ \| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{u})
308 \| } \equiv \; \Psi_{i} \;\; {\hat{\mathbf{v}} \times
309 (\mathbf{x}_i-\mathbf{u} ) }\\
310 \Psi_i^{*} & := & \frac{1}{ \| \hat{\mathbf{v}} \times
311 (\mathbf{x}_i-\mathbf{u}) \|^2 + \epsilon'}\end{aligned}
312 :label: eqnpotrm2forcepart1
314 the force on atom :math:`j` reads
316 .. math:: \mathbf{F}_{\!j}^\mathrm{rm2} =
319 (\mathbf{s}_j\cdot\mathbf{r}_{\!j})\;
320 \left[ \frac{\Psi_{\!j}^* }{\Psi_{\!j} } \mathbf{r}_{\!j}
321 - \frac{\Psi_{\!j}^{ * 2}}{\Psi_{\!j}^3}
322 (\mathbf{s}_j\cdot\mathbf{r}_{\!j})\mathbf{s}_j \right]
323 \right\rbrace \times \hat{\mathbf{v}} .
324 :label: eqnpotrm2force
326 Pivot-Free Radial Motion 2 Potential
327 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
329 The pivot-free variant of the above potential is
331 .. math:: V{^\mathrm{rm2-pf}}=
332 \frac{k}{2} \sum_{i=1}^{N} w_i\,
333 \frac{\left[ (\hat{\mathbf{v}} \times ( \mathbf{x}_i - \mathbf{x}_c ))
334 \cdot \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c) \right]^2}
335 {\| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c) \|^2 +
341 .. math:: \begin{aligned}
342 \mathbf{r}_i & := & \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c)\\
343 \mathbf{s}_i & := & \frac{\hat{\mathbf{v}} \times (\mathbf{x}_i -
344 \mathbf{x}_c ) }{ \| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c)
345 \| } \equiv \; \Psi_{i} \;\; {\hat{\mathbf{v}} \times
346 (\mathbf{x}_i-\mathbf{x}_c ) }\\ \Psi_i^{*} & := & \frac{1}{ \| \hat{\mathbf{v}} \times
347 (\mathbf{x}_i-\mathbf{x}_c) \|^2 + \epsilon'}\end{aligned}
348 :label: eqnpotrm2pfpart2
350 the force on atom :math:`j` reads
352 .. math:: \begin{aligned}
354 \mathbf{F}_{\!j}{^\mathrm{rm2-pf}}& = &
357 (\mathbf{s}_j\cdot\mathbf{r}_{\!j})\;
358 \left[ \frac{\Psi_{\!j}^* }{\Psi_{\!j} } \mathbf{r}_{\!j}
359 - \frac{\Psi_{\!j}^{ * 2}}{\Psi_{\!j}^3}
360 (\mathbf{s}_j\cdot\mathbf{r}_{\!j})\mathbf{s}_j \right]
361 \right\rbrace \times \hat{\mathbf{v}}\\
363 + k\;\frac{m_j}{M} \left\lbrace \sum_{i=1}^{N}
364 w_i\;(\mathbf{s}_i\cdot\mathbf{r}_i) \;
365 \left[ \frac{\Psi_i^* }{\Psi_i } \mathbf{r}_i
366 - \frac{\Psi_i^{ * 2}}{\Psi_i^3} (\mathbf{s}_i\cdot\mathbf{r}_i )\;
367 \mathbf{s}_i \right] \right\rbrace \times \hat{\mathbf{v}} \, .
369 :label: eqnpotrm2pfforce
371 Flexible Axis Rotation
372 ~~~~~~~~~~~~~~~~~~~~~~
374 As sketched in :numref:`Fig. %s <fig-rotation>` A–B, the rigid body
375 behavior of the fixed axis rotation scheme is a drawback for many
376 applications. In particular, deformations of the rotation group are
377 suppressed when the equilibrium atom positions directly depend on the
378 reference positions. To avoid this limitation,
379 eqns. :eq:`%s <eqnpotrmpf>` and :eq:`%s <eqnpotrm2pf>`
380 will now be generalized towards a “flexible axis” as sketched in
381 :numref:`Fig. %s C <fig-rotation>`. This will be achieved by
382 subdividing the rotation group into a set of equidistant slabs
383 perpendicular to the rotation vector, and by applying a separate
384 rotation potential to each of these slabs.
385 :numref:`Fig. %s C <fig-rotation>` shows the midplanes of the slabs
386 as dotted straight lines and the centers as thick black dots.
388 To avoid discontinuities in the potential and in the forces, we define
389 “soft slabs” by weighing the contributions of each slab :math:`n` to the
390 total potential function :math:`V^\mathrm{flex}` by a Gaussian function
392 .. math:: g_n(\mathbf{x}_i) = \Gamma \ \mbox{exp} \left(
393 -\frac{\beta_n^2(\mathbf{x}_i)}{2\sigma^2} \right) ,
396 centered at the midplane of the :math:`n`\ th slab. Here :math:`\sigma`
397 is the width of the Gaussian function, :math:`\Delta x` the distance
398 between adjacent slabs, and
400 .. math:: \beta_n(\mathbf{x}_i) := \mathbf{x}_i \cdot \hat{\mathbf{v}} - n \, \Delta x \, .
401 :label: eqngaussianpart2
405 .. figure:: plots/gaussians.*
408 Gaussian functions :math:`g_n` centered at :math:`n \, \Delta x` for
409 a slab distance :math:`\Delta x = 1.5` nm and :math:`n \geq -2`.
410 Gaussian function :math:`g_0` is highlighted in bold; the dashed line
411 depicts the sum of the shown Gaussian functions.
413 A most convenient choice is :math:`\sigma = 0.7 \Delta x` and
415 .. math:: 1/\Gamma = \sum_{n \in Z}
417 \left(-\frac{(n - \frac{1}{4})^2}{2\cdot 0.7^2}\right)
419 :label: eqngaussianpart3
421 which yields a nearly constant sum, essentially independent of
422 :math:`\mathbf{x}_i` (dashed line in
423 :numref:`Fig. %s <fig-gaussian>`), *i.e.*,
425 .. math:: \sum_{n \in Z} g_n(\mathbf{x}_i) = 1 + \epsilon(\mathbf{x}_i) \, ,
429 :math:`| \epsilon(\mathbf{x}_i) | < 1.3\cdot 10^{-4}`.
430 This choice also implies that the individual contributions to the force
431 from the slabs add up to unity such that no further normalization is
434 To each slab center :math:`\mathbf{x}_c^n`, all atoms
435 contribute by their Gaussian-weighted (optionally also mass-weighted)
437 :math:`g_n(\mathbf{x}_i) \, \mathbf{x}_i`.
438 The instantaneous slab centers :math:`\mathbf{x}_c^n` are
439 calculated from the current positions
440 :math:`\mathbf{x}_i`,
442 .. math:: \mathbf{x}_c^n =
443 \frac{\sum_{i=1}^N g_n(\mathbf{x}_i) \, m_i \, \mathbf{x}_i}
444 {\sum_{i=1}^N g_n(\mathbf{x}_i) \, m_i} \, ,\\
447 while the reference centers :math:`\mathbf{y}_c^n` are
448 calculated from the reference positions
449 :math:`\mathbf{y}_i^0`,
451 .. math:: \mathbf{y}_c^n =
452 \frac{\sum_{i=1}^N g_n(\mathbf{y}_i^0) \, m_i \, \mathbf{y}_i^0}
453 {\sum_{i=1}^N g_n(\mathbf{y}_i^0) \, m_i} \, .
456 Due to the rapid decay of :math:`g_n`, each slab will essentially
457 involve contributions from atoms located within :math:`\approx
458 3\Delta x` from the slab center only.
460 Flexible Axis Potential
461 ^^^^^^^^^^^^^^^^^^^^^^^
463 We consider two flexible axis variants. For the first variant, the slab
464 segmentation procedure with Gaussian weighting is applied to the radial
466 (eqn. :eq:`%s <eqnpotrmpf>` / :numref:`Fig. %s B <fig-equipotential>`),
467 yielding as the contribution of slab :math:`n`
470 \frac{k}{2} \sum_{i=1}^{N} w_i \, g_n(\mathbf{x}_i)
474 (\mathbf{x}_i - \mathbf{x}_c^n)
478 and a total potential function
480 .. math:: V^\mathrm{flex} = \sum_n V^n \, .
483 Note that the global center of mass :math:`\mathbf{x}_c`
484 used in eqn. :eq:`%s <eqnpotrmpf>` is now replaced by
485 :math:`\mathbf{x}_c^n`, the center of mass of the slab.
488 .. math:: \begin{aligned}
489 \mathbf{q}_i^n & := & \frac{\hat{\mathbf{v}} \times
490 \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) }{ \| \hat{\mathbf{v}}
491 \times \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) \| } \\
492 b_i^n & := & \mathbf{q}_i^n \cdot (\mathbf{x}_i - \mathbf{x}_c^n) \, ,\end{aligned}
493 :label: eqnflexpotpart2
495 the resulting force on atom :math:`j` reads
497 .. math:: \begin{aligned}
498 \nonumber\hspace{-15mm}
499 \mathbf{F}_{\!j}^\mathrm{flex} &=&
500 - \, k \, w_j \sum_n g_n(\mathbf{x}_j) \, b_j^n \left\lbrace \mathbf{q}_j^n -
501 b_j^n \frac{\beta_n(\mathbf{x}_j)}{2\sigma^2} \hat{\mathbf{v}} \right\rbrace \\ & &
502 + \, k \, m_j \sum_n \frac{g_n(\mathbf{x}_j)}{\sum_h g_n(\mathbf{x}_h)}
503 \sum_{i=1}^{N} w_i \, g_n(\mathbf{x}_i) \, b_i^n \left\lbrace
504 \mathbf{q}_i^n -\frac{\beta_n(\mathbf{x}_j)}{\sigma^2}
505 \left[ \mathbf{q}_i^n \cdot (\mathbf{x}_j - \mathbf{x}_c^n )\right]
506 \hat{\mathbf{v}} \right\rbrace .
508 :label: eqnpotflexforce
510 Note that for :math:`V^\mathrm{flex}`, as defined, the slabs are fixed in
511 space and so are the reference centers
512 :math:`\mathbf{y}_c^n`. If during the simulation the
513 rotation group moves too far in :math:`\mathbf{v}`
514 direction, it may enter a region where – due to the lack of nearby
515 reference positions – no reference slab centers are defined, rendering
516 the potential evaluation impossible. We therefore have included a
517 slightly modified version of this potential that avoids this problem by
518 attaching the midplane of slab :math:`n=0` to the center of mass of the
519 rotation group, yielding slabs that move with the rotation group. This
520 is achieved by subtracting the center of mass
521 :math:`\mathbf{x}_c` of the group from the positions,
523 .. math:: \tilde{\mathbf{x}}_i = \mathbf{x}_i - \mathbf{x}_c \, , \mbox{\ \ \ and \ \ }
524 \tilde{\mathbf{y}}_i^0 = \mathbf{y}_i^0 - \mathbf{y}_c^0 \, ,
529 .. math:: \begin{aligned}
531 & = & \frac{k}{2} \sum_n \sum_{i=1}^{N} w_i \, g_n(\tilde{\mathbf{x}}_i)
532 \left[ \frac{\hat{\mathbf{v}} \times \mathbf{\Omega}(t)(\tilde{\mathbf{y}}_i^0
533 - \tilde{\mathbf{y}}_c^n) }{ \| \hat{\mathbf{v}} \times
534 \mathbf{\Omega}(t)(\tilde{\mathbf{y}}_i^0 -
535 \tilde{\mathbf{y}}_c^n) \| }
537 (\tilde{\mathbf{x}}_i - \tilde{\mathbf{x}}_c^n)
542 To simplify the force derivation, and for efficiency reasons, we here
543 assume :math:`\mathbf{x}_c` to be constant, and thus
544 :math:`\partial \mathbf{x}_c / \partial x =
545 \partial \mathbf{x}_c / \partial y = \partial \mathbf{x}_c / \partial z = 0`.
546 The resulting force error is small (of order :math:`O(1/N)` or
547 :math:`O(m_j/M)` if mass-weighting is applied) and can therefore be
548 tolerated. With this assumption, the forces :math:`\mathbf{F}^\mathrm{flex-t}`
549 have the same form as eqn. :eq:`%s <eqnpotflexforce>`.
551 Flexible Axis 2 Alternative Potential
552 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
554 In this second variant, slab segmentation is applied to
555 :math:`V^\mathrm{rm2}` (eqn. :eq:`%s <eqnpotrm2pf>`), resulting in
556 a flexible axis potential without radial force contributions
557 (:numref:`Fig. %s C <fig-equipotential>`),
559 .. math:: V{^\mathrm{flex2}}=
560 \frac{k}{2} \sum_{i=1}^{N} \sum_n w_i\,g_n(\mathbf{x}_i)
561 \frac{\left[ (\hat{\mathbf{v}} \times ( \mathbf{x}_i - \mathbf{x}_c^n ))
562 \cdot \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n) \right]^2}
563 {\| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c^n) \|^2 +
569 .. math:: \begin{aligned}
570 \mathbf{r}_i^n & := & \mathbf{\Omega}(t)(\mathbf{y}_i^0 - \mathbf{y}_c^n)\\
571 \mathbf{s}_i^n & := & \frac{\hat{\mathbf{v}} \times (\mathbf{x}_i -
572 \mathbf{x}_c^n ) }{ \| \hat{\mathbf{v}} \times (\mathbf{x}_i - \mathbf{x}_c^n)
573 \| } \equiv \; \psi_{i} \;\; {\hat{\mathbf{v}} \times (\mathbf{x}_i-\mathbf{x}_c^n ) }\\
574 \psi_i^{*} & := & \frac{1}{ \| \hat{\mathbf{v}} \times (\mathbf{x}_i-\mathbf{x}_c^n) \|^2 + \epsilon'}\\
575 W_j^n & := & \frac{g_n(\mathbf{x}_j)\,m_j}{\sum_h g_n(\mathbf{x}_h)\,m_h}\\
577 \sum_{i=1}^{N} w_i\;g_n(\mathbf{x}_i)
578 \; (\mathbf{s}_i^n\cdot\mathbf{r}_i^n)
579 \left[ \frac{\psi_i^* }{\psi_i } \mathbf{r}_i^n
580 - \frac{\psi_i^{ * 2}}{\psi_i^3} (\mathbf{s}_i^n\cdot\mathbf{r}_i^n )\;
581 \mathbf{s}_i^n \right]
585 the force on atom :math:`j` reads
587 .. math:: \begin{aligned}
589 \mathbf{F}_{\!j}{^\mathrm{flex2}}& = &
591 \left\lbrace \sum_n w_j\;g_n(\mathbf{x}_j)\;
592 (\mathbf{s}_j^n\cdot\mathbf{r}_{\!j}^n)\;
593 \left[ \frac{\psi_j^* }{\psi_j } \mathbf{r}_{\!j}^n
594 - \frac{\psi_j^{ * 2}}{\psi_j^3} (\mathbf{s}_j^n\cdot\mathbf{r}_{\!j}^n)\;
595 \mathbf{s}_{\!j}^n \right] \right\rbrace \times \hat{\mathbf{v}} \\
598 + k \left\lbrace \sum_n W_{\!j}^n \, \mathbf{S}^n \right\rbrace \times
600 - k \left\lbrace \sum_n W_{\!j}^n \; \frac{\beta_n(\mathbf{x}_j)}{\sigma^2} \frac{1}{\psi_j}\;\;
602 \mathbf{S}^n \right\rbrace \hat{\mathbf{v}}\\
604 + \frac{k}{2} \left\lbrace \sum_n w_j\;g_n(\mathbf{x}_j)
605 \frac{\beta_n(\mathbf{x}_j)}{\sigma^2}
606 \frac{\psi_j^*}{\psi_j^2}( \mathbf{s}_j^n \cdot \mathbf{r}_{\!j}^n )^2 \right\rbrace
609 :label: eqnpotflex2force
611 Applying transformation :eq:`%s <eqntrafo>` yields a
612 “translation-tolerant” version of the flexible2 potential,
613 :math:`V{^\mathrm{flex2 - t}}`. Again, assuming that
614 :math:`\partial \mathbf{x}_c / \partial x`,
615 :math:`\partial \mathbf{x}_c /
616 \partial y`, :math:`\partial \mathbf{x}_c / \partial z`
617 are small, the resulting equations for :math:`V{^\mathrm{flex2 - t}}`
618 and :math:`\mathbf{F}{^\mathrm{flex2 - t}}` are similar
619 to those of :math:`V^\mathrm{flex2}` and
620 :math:`\mathbf{F}^\mathrm{flex2}`.
625 To apply enforced rotation, the particles :math:`i` that are to be
626 subjected to one of the rotation potentials are defined via index groups
627 ``rot-group0``, ``rot-group1``, etc., in the
628 :ref:`mdp` input file. The reference positions
629 :math:`\mathbf{y}_i^0` are read from a special
630 :ref:`trr` file provided to :ref:`grompp <gmx grompp>`. If no such
631 file is found, :math:`\mathbf{x}_i(t=0)` are used as
632 reference positions and written to :ref:`trr` such that they
633 can be used for subsequent setups. All parameters of the potentials such
634 as :math:`k`, :math:`\epsilon'`, etc.
635 (:numref:`Table %s <tab-vars>`) are provided as :ref:`mdp`
636 parameters; ``rot-type`` selects the type of the potential.
637 The option ``rot-massw`` allows to choose whether or not to
638 use mass-weighted averaging. For the flexible potentials, a cutoff value
639 :math:`g_n^\mathrm{min}` (typically :math:`g_n^\mathrm{min}=0.001`)
640 makes sure that only significant contributions to :math:`V` and
641 :math:`\mathbf{F}` are evaluated, *i.e.* terms with
642 :math:`g_n(\mathbf{x}) < g_n^\mathrm{min}` are omitted.
643 :numref:`Table %s <tab-quantities>` summarizes observables that are
644 written to additional output files and which are described below.
646 .. |ROTISO| replace:: V\ :math:`^\mathrm{iso}`
647 .. |ROTISOPF| replace:: V\ :math:`^\mathrm{iso-pf}`
648 .. |ROTPM| replace:: V\ :math:`^\mathrm{pm}`
649 .. |ROTPMPF| replace:: V\ :math:`^\mathrm{pm-pf}`
650 .. |ROTRM| replace:: V\ :math:`^\mathrm{rm}`
651 .. |ROTRMPF| replace:: V\ :math:`^\mathrm{rm-pf}`
652 .. |ROTRM2| replace:: V\ :math:`^\mathrm{rm2}`
653 .. |ROTRM2PF| replace:: V\ :math:`^\mathrm{rm2-pf}`
654 .. |ROTFL| replace:: V\ :math:`^\mathrm{flex}`
655 .. |ROTFLT| replace:: V\ :math:`^\mathrm{flex-t}`
656 .. |ROTFL2| replace:: V\ :math:`^\mathrm{flex2}`
657 .. |ROTFLT2| replace:: V\ :math:`^\mathrm{flex2-t}`
658 .. |KUNIT| replace:: :math:`\frac{\mathrm{kJ}}{\mathrm{mol} \cdot \mathrm{nm}^2}`
659 .. |BFX| replace:: **x**
660 .. |KMA| replace:: :math:`k`
661 .. |VECV| replace:: :math:`\hat{\mathbf{v}}`
662 .. |VECU| replace:: :math:`\mathbf{u}`
663 .. |OMEG| replace:: :math:`\omega`
664 .. |EPSP| replace:: :math:`{\epsilon}'`
665 .. |DELX| replace:: :math:`{\Delta}x`
666 .. |GMIN| replace:: :math:`g_n^\mathrm{min}`
667 .. |CIPS| replace:: :math:`^\circ`\ /ps
668 .. |NM2| replace:: nm\ :math:`^2`
669 .. |REF1| replace:: \ :eq:`eqnpotiso`
670 .. |REF2| replace:: \ :eq:`eqnpotisopf`
671 .. |REF3| replace:: \ :eq:`eqnpotpm`
672 .. |REF4| replace:: \ :eq:`eqnpotpmpf`
673 .. |REF5| replace:: \ :eq:`eqnpotrm`
674 .. |REF6| replace:: \ :eq:`eqnpotrmpf`
675 .. |REF7| replace:: \ :eq:`eqnpotrm2`
676 .. |REF8| replace:: \ :eq:`eqnpotrm2pf`
677 .. |REF9| replace:: \ :eq:`eqnpotflex`
678 .. |REF10| replace:: \ :eq:`eqnpotflext`
679 .. |REF11| replace:: \ :eq:`eqnpotflex2`
683 .. table:: Parameters used by the various rotation potentials.
684 |BFX| indicate which parameter is actually used for a given potential
688 +------------------------------------------+---------+--------+--------+--------+--------+-----------+-----------+
689 | parameter | |KMA| | |VECV| | |VECU| | |OMEG| | |EPSP| | |DELX| | |GMIN| |
690 +------------------------------------------+---------+--------+--------+--------+--------+-----------+-----------+
691 | :ref:`mdp` input variable name | k | vec | pivot | rate | eps | slab-dist | min-gauss |
692 +------------------------------------------+---------+--------+--------+--------+--------+-----------+-----------+
693 | unit | |KUNIT| | ``-`` | nm | |CIPS| | |NM2| | nm | ``-`` |
694 +================================+=========+=========+========+========+========+========+===========+===========+
695 | fixed axis potentials: | eqn. |
696 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
697 | isotropic | |ROTISO| | |REF1| | |BFX| | |BFX| | |BFX| | |BFX| | ``-`` | ``-`` | ``-`` |
698 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
699 | --- pivot-free | |ROTISOPF| | |REF2| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | ``-`` | ``-`` |
700 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
701 | parallel motion | |ROTPM| | |REF3| | |BFX| | |BFX| | |BFX| | |BFX| | ``-`` | ``-`` | ``-`` |
702 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
703 | --- pivot-free | |ROTPMPF| | |REF4| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | ``-`` | ``-`` |
704 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
705 | radial motion | |ROTRM| | |REF5| | |BFX| | |BFX| | |BFX| | |BFX| | ``-`` | ``-`` | ``-`` |
706 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
707 | --- pivot-free | |ROTRMPF| | |REF6| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | ``-`` | ``-`` |
708 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
709 | radial motion 2 | |ROTRM2| | |REF7| | |BFX| | |BFX| | |BFX| | |BFX| | |BFX| | ``-`` | ``-`` |
710 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
711 | --- pivot-free | |ROTRM2PF| | |REF8| | |BFX| | |BFX| | ``-`` | |BFX| | |BFX| | ``-`` | ``-`` |
712 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
713 | flexible axis potentials: | eqn. |
714 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
715 | flexible | |ROTFL| | |REF9| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | |BFX| | |BFX| |
716 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
717 | --- transl. tol | |ROTFLT| | |REF10| | |BFX| | |BFX| | ``-`` | |BFX| | ``-`` | |BFX| | |BFX| |
718 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
719 | flexible 2 | |ROTFL2| | |REF11| | |BFX| | |BFX| | ``-`` | |BFX| | |BFX| | |BFX| | |BFX| |
720 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
721 | --- transl. tol | |ROTFLT2| | ``-`` | |BFX| | |BFX| | ``-`` | |BFX| | |BFX| | |BFX| | |BFX| |
722 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+
727 .. |VT| replace:: :math:`V(t)`
728 .. |THET| replace:: :math:`\theta_{ \mathrm{ref}}(t)`
729 .. |THETAV| replace:: :math:`\theta_{\mathrm{av}}(t)`
730 .. |THETFIT| replace:: :math:`\theta_{\mathrm{fit}}(t)`, :math:`\theta_{\mathrm{fit}}(t,n)`
731 .. |YVEC| replace:: :math:`\mathbf{y}_{0}(n)`, :math:`\mathbf{x}_{0}(t,n)`
732 .. |TAUT| replace:: :math:`\tau(t)`
733 .. |TAUTN| replace:: :math:`\tau(t,n)`
734 .. |REFT| replace:: :numref:`see Table %s <tab-vars>`
735 .. |REFEQ| replace:: :math:`\theta_{\mathrm{ref}}(t)=\omega t`
736 .. |REF12| replace:: \ :eq:`eqnavangle`
737 .. |REF13| replace:: \ :eq:`eqnrmsdfit`
738 .. |REF14| replace:: \ :eq:`eqndefx0`\ ,\ :eq:`eqndefy0`
739 .. |REF15| replace:: \ :eq:`eqntorque`
743 .. table:: Quantities recorded in output files during enforced rotation simulations.
744 All slab-wise data is written every ``nstsout`` steps, other rotation data every ``nstrout`` steps.
748 +------------+---------+------------+--------------------+-------+----------+
749 | quantity | unit | equation | output file | fixed | flexible |
750 +============+=========+============+====================+=======+==========+
751 | |VT| | kJ/mol | |REFT| | ``rotation`` | |BFX| | |BFX| |
752 +------------+---------+------------+--------------------+-------+----------+
753 | |THET| | degrees | |REFEQ| | ``rotation`` | |BFX| | |BFX| |
754 +------------+---------+------------+--------------------+-------+----------+
755 | |THETAV| | degrees | |REF12| | ``rotation`` | |BFX| | ``-`` |
756 +------------+---------+------------+--------------------+-------+----------+
757 | |THETFIT| | degrees | |REF13| | ``rotangles`` | ``-`` | |BFX| |
758 +------------+---------+------------+--------------------+-------+----------+
759 | |YVEC| | nm | |REF14| | ``rotslabs`` | ``-`` | |BFX| |
760 +------------+---------+------------+--------------------+-------+----------+
761 | |TAUT| | kJ/mol | |REF15| | ``rotation`` | |BFX| | ``-`` |
762 +------------+---------+------------+--------------------+-------+----------+
763 | |TAUTN| | kJ/mol | |REF15| | ``rottorque`` | ``-`` | |BFX| |
764 +------------+---------+------------+--------------------+-------+----------+
769 Angle of Rotation Groups: Fixed Axis
770 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
772 For fixed axis rotation, the average angle :math:`\theta_\mathrm{av}(t)`
773 of the group relative to the reference group is determined via the
774 distance-weighted angular deviation of all rotation group atoms from
775 their reference positions,
777 .. math:: \theta_\mathrm{av} = \left. \sum_{i=1}^{N} r_i \ \theta_i \right/ \sum_{i=1}^N r_i \ .
780 Here, :math:`r_i` is the distance of the reference position to the
781 rotation axis, and the difference angles :math:`\theta_i` are determined
782 from the atomic positions, projected onto a plane perpendicular to the
783 rotation axis through pivot point :math:`\mathbf{u}` (see
784 eqn. :eq:`%s <eqnproject>` for the definition of
787 .. math:: \cos \theta_i =
788 \frac{(\mathbf{y}_i-\mathbf{u})^\perp \cdot (\mathbf{x}_i-\mathbf{u})^\perp}
789 { \| (\mathbf{y}_i-\mathbf{u})^\perp \cdot (\mathbf{x}_i-\mathbf{u})^\perp
791 :label: eqnavanglepart2
793 The sign of :math:`\theta_\mathrm{av}` is chosen such that
794 :math:`\theta_\mathrm{av} > 0` if the actual structure rotates ahead of
797 Angle of Rotation Groups: Flexible Axis
798 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
800 For flexible axis rotation, two outputs are provided, the angle of the
801 entire rotation group, and separate angles for the segments in the
802 slabs. The angle of the entire rotation group is determined by an RMSD
803 fit of :math:`\mathbf{x}_i` to the reference positions
804 :math:`\mathbf{y}_i^0` at :math:`t=0`, yielding
805 :math:`\theta_\mathrm{fit}` as the angle by which the reference has to
806 be rotated around :math:`\hat{\mathbf{v}}` for the optimal
809 .. math:: \mathrm{RMSD} \big( \mathbf{x}_i,\ \mathbf{\Omega}(\theta_\mathrm{fit})
810 \mathbf{y}_i^0 \big) \stackrel{!}{=} \mathrm{min} \, .
813 To determine the local angle for each slab :math:`n`, both reference
814 and actual positions are weighted with the Gaussian function of slab
815 :math:`n`, and :math:`\theta_\mathrm{fit}(t,n)` is calculated as in
816 eqn. :eq:`%s <eqnrmsdfit>` from the Gaussian-weighted
819 For all angles, the :ref:`mdp` input option
820 ``rot-fit-method`` controls whether a normal RMSD fit is
821 performed or whether for the fit each position
822 :math:`\mathbf{x}_i` is put at the same distance to the
823 rotation axis as its reference counterpart
824 :math:`\mathbf{y}_i^0`. In the latter case, the RMSD
825 measures only angular differences, not radial ones.
827 Angle Determination by Searching the Energy Minimum
828 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
830 Alternatively, for ``rot-fit-method = potential``, the angle
831 of the rotation group is determined as the angle for which the rotation
832 potential energy is minimal. Therefore, the used rotation potential is
833 additionally evaluated for a set of angles around the current reference
834 angle. In this case, the ``rotangles.log`` output file
835 contains the values of the rotation potential at the chosen set of
836 angles, while ``rotation.xvg`` lists the angle with minimal
842 The torque :math:`\mathbf{\tau}(t)` exerted by the
843 rotation potential is calculated for fixed axis rotation via
845 .. math:: \mathbf{\tau}(t) = \sum_{i=1}^{N} \mathbf{r}_i(t) \times \mathbf{f}_{\!i}^\perp(t) ,
848 where :math:`\mathbf{r}_i(t)` is the distance vector from
849 the rotation axis to :math:`\mathbf{x}_i(t)` and
850 :math:`\mathbf{f}_{\!i}^\perp(t)` is the force component
851 perpendicular to :math:`\mathbf{r}_i(t)` and
852 :math:`\hat{\mathbf{v}}`. For flexible axis rotation,
853 torques :math:`\mathbf{\tau}_{\!n}` are calculated for
854 each slab using the local rotation axis of the slab and the
855 Gaussian-weighted positions.