1Enforced Rotation 2----------------- 3 4This 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, 6among them complex ones that allow flexible adaptations of both the 7rotated subunit as well as the local rotation axis during the 8simulation. An example application can be found in ref. 9:ref:`145 <refKutzner2011>`. 10 11.. _fig-rotation: 12 13.. figure:: plots/rotation.* 14 :width: 13.00000cm 15 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 27 :math:`\Delta x`. 28 29.. _fig-equipotential: 30 31.. figure:: plots/equipotential.* 32 :width: 13.00000cm 33 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. 51 52Fixed Axis Rotation 53^^^^^^^^^^^^^^^^^^^ 54 55Stationary Axis with an Isotropic Potential 56~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 57 58In the fixed axis approach (see :numref:`Fig. %s B <fig-rotation>`), 59torque on a group of :math:`N` atoms with positions 60:math:`\mathbf{x}_i` (denoted “rotation group”) is applied 61by rotating a reference set of atomic positions – usually their initial 62positions :math:`\mathbf{y}_i^0` – at a constant angular 63velocity :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 66position :math:`\mathbf{x}_i` is attracted by a “virtual 67spring” potential to its moving reference position 68:math:`\mathbf{y}_i = \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{u})`, 69where :math:`\mathbf{\Omega}(t)` is a matrix that describes the rotation 70around the axis. In the simplest case, the “springs” are described by a 71harmonic potential, 72 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 75 :label: eqnpotiso 76 77with optional mass-weighted prefactors :math:`w_i = N \, m_i/M` with 78total mass :math:`M = \sum_{i=1}^N m_i`. The rotation matrix 79:math:`\mathbf{\Omega}(t)` is 80 81.. math:: \mathbf{\Omega}(t) = 82 \left( 83 \begin{array}{ccc} 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\,}\\ 87 \end{array} 88 \right) 89 :label: eqnrotmat 90 91where :math:`v_x`, :math:`v_y`, and :math:`v_z` are the components of 92the normalized rotation vector :math:`\hat{\mathbf{v}}`, 93and :math:`{\,\xi\,}:= 1-\cos(\omega t)`. As illustrated in 94:numref:`Fig. %s A <fig-equipotential>` for a single atom :math:`j`, 95the rotation matrix :math:`\mathbf{\Omega}(t)` operates on the initial 96reference positions 97:math:`\mathbf{y}_j^0 = \mathbf{x}_j(t_0)` 98of atom :math:`j` at :math:`t=t_0`. At a later time :math:`t`, the 99reference position has rotated away from its initial place (along the 100blue dashed line), resulting in the force 101 102.. math:: \mathbf{F}_{\!j}^\mathrm{iso} 103 = -\nabla_{\!j} \, V^\mathrm{iso} 104 = k \, w_j \left[ 105 \mathbf{\Omega}(t) (\mathbf{y}_j^0 - \mathbf{u}) - (\mathbf{x}_j - \mathbf{u} ) \right] 106 :label: eqnforcefixed 107 108which is directed towards the reference position. 109 110Pivot-Free Isotropic Potential 111^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 112 113Instead of a fixed pivot vector :math:`\mathbf{u}` this 114potential uses the center of mass :math:`\mathbf{x}_c` of 115the rotation group as pivot for the rotation axis, 116 117.. math:: \mathbf{x}_c = \frac{1}{M} \sum_{i=1}^N m_i \mathbf{x}_i 118 \mbox{and} 119 \mathbf{y}_c^0 = \frac{1}{M} \sum_{i=1}^N m_i \mathbf{y}_i^0 \ , 120 :label: eqncom 121 122which yields the “pivot-free” isotropic potential 123 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 , 126 :label: eqnpotisopf 127 128with forces 129 130.. math:: \mathbf{F}_{\!j}^\mathrm{iso-pf} = k \, w_j 131 \left[ 132 \mathbf{\Omega}(t) ( \mathbf{y}_j^0 - \mathbf{y}_c^0) 133 - ( \mathbf{x}_j - \mathbf{x}_c ) 134 \right] . 135 :label: eqnforceisopf 136 137Without mass-weighting, the pivot :math:`\mathbf{x}_c` is 138the geometrical center of the group. 139 140Parallel Motion Potential Variant 141^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 142 143The forces generated by the isotropic potentials 144(eqns. :eq:`%s <eqnpotiso>` and :eq:`%s <eqnpotisopf>`) also contain components parallel to the 145rotation axis and thereby restrain motions along the axis of either the 146whole rotation group (in case of :math:`V^\mathrm{iso}`) or within the 147rotation group, in case of :math:`V^\mathrm{iso-pf}`. 148 149For cases where unrestrained motion along the axis is preferred, we have implemented a 150“parallel motion” variant by eliminating all components parallel to the 151rotation axis for the potential. This is achieved by projecting the 152distance vectors between reference and actual positions 153 154.. math:: \mathbf{r}_i = \mathbf{\Omega}(t) (\mathbf{y}_i^0 - \mathbf{u}) - (\mathbf{x}_i - \mathbf{u}) 155 :label: eqnrotdistvectors 156 157onto the plane perpendicular to the rotation vector, 158 159.. math:: \mathbf{r}_i^\perp := \mathbf{r}_i - (\mathbf{r}_i \cdot \hat{\mathbf{v}})\hat{\mathbf{v}} 160 :label: eqnproject 161 162yielding 163 164.. math:: \begin{aligned} 165 \nonumber 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 168 \left\lbrace 169 \mathbf{\Omega}(t) 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 174 \end{aligned} 175 :label: eqnpotpm 176 177and similarly 178 179.. math:: \mathbf{F}_{\!j}^\mathrm{pm} = k \, w_j \, \mathbf{r}_j^\perp 180 :label: eqnforcepm 181 182Pivot-Free Parallel Motion Potential 183^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 184 185Replacing 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 188parallel motion potential. With 189 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 192 193the respective potential and forces are 194 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} 197 :label: eqnpotpmpf 198 199.. math:: \begin{aligned} 200 \mathbf{F}_{\!j}^\mathrm{pm-pf} &=& k \, w_j \, \mathbf{s}_j^\perp 201 \end{aligned} 202 :label: eqnforcepmpf 203 204Radial Motion Potential 205^^^^^^^^^^^^^^^^^^^^^^^ 206 207In the above variants, the minimum of the rotation potential is either a 208single point at the reference position 209:math:`\mathbf{y}_i` (for the isotropic potentials) or a 210single line through :math:`\mathbf{y}_i` parallel to the 211rotation axis (for the parallel motion potentials). As a result, radial 212forces restrict radial motions of the atoms. The two subsequent types of 213rotation potentials, :math:`V^\mathrm{rm}` and :math:`V^\mathrm{rm2}`, drastically 214reduce or even eliminate this effect. The first variant, :math:`V^\mathrm{rm}` 215(:numref:`Fig. %s B <fig-equipotential>`), eliminates all force 216components parallel to the vector connecting the reference atom and the 217rotation axis, 218 219.. math:: V^\mathrm{rm} = \frac{k}{2} \sum_{i=1}^{N} w_i \left[ 220 \mathbf{p}_i 221 \cdot(\mathbf{x}_i - \mathbf{u}) \right]^2 , 222 :label: eqnpotrm 223 224with 225 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 229 230This 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})`. 235The resulting force is 236 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 240 241Pivot-Free Radial Motion Potential 242^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 243 244Proceeding similar to the pivot-free isotropic potential yields a 245pivot-free version of the above potential. With 246 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 250 251the potential and force for the pivot-free variant of the radial motion 252potential read 253 254.. math:: \begin{aligned} 255 V^\mathrm{rm-pf} & = & \frac{k}{2} \sum_{i=1}^{N} w_i \left[ 256 \mathbf{q}_i 257 \cdot(\mathbf{x}_i - \mathbf{x}_c) 258 \right]^2 \, , \end{aligned} 259 :label: eqnpotrmpf 260 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 \, . 266 \end{aligned} 267 :label: eqnpotrmpfforce 268 269Radial Motion 2 Alternative Potential 270^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 271 272As seen in :numref:`Fig. %s B <fig-equipotential>`, the force 273resulting from :math:`V^\mathrm{rm}` still contains a small, second-order 274radial component. In most cases, this perturbation is tolerable; if not, 275the following alternative, :math:`V^\mathrm{rm2}`, fully eliminates the 276radial contribution to the force, as depicted in 277:numref:`Fig. %s C <fig-equipotential>`, 278 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 + 284 \epsilon'} \, , 285 :label: eqnpotrm2 286 287where a small parameter :math:`\epsilon'` has been introduced to avoid 288singularities. For :math:`\epsilon'\mathrm{ = }0\mathrm{nm}^2`, the 289equipotential planes are spanned by :math:`\mathbf{x}_i - 290\mathbf{u}` and :math:`\hat{\mathbf{v}}`, 291yielding a force perpendicular to 292:math:`\mathbf{x}_i - \mathbf{u}`, thus not 293contracting or expanding structural parts that moved away from or toward 294the rotation axis. 295 296Choosing 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 299eqn. :eq:`%s <eqnpotrm2>` yields a well-defined potential and 300continuous forces also close to the rotation axis, which is not the case 301for :math:`\epsilon'\mathrm{ = }0\mathrm{nm}^2` 302(:numref:`Fig. %s C <fig-equipotential>`). With 303 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 313 314the force on atom :math:`j` reads 315 316.. math:: \mathbf{F}_{\!j}^\mathrm{rm2} = 317 - k\; 318 \left\lbrace w_j\; 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 325 326Pivot-Free Radial Motion 2 Potential 327^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 328 329The pivot-free variant of the above potential is 330 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 + 336 \epsilon'} \, . 337 :label: eqnpotrm2pf 338 339With 340 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 349 350the force on atom :math:`j` reads 351 352.. math:: \begin{aligned} 353 \nonumber 354 \mathbf{F}_{\!j}{^\mathrm{rm2-pf}}& = & 355 - k\; 356 \left\lbrace w_j\; 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}}\\ 362 & & 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}} \, . 368 \end{aligned} 369 :label: eqnpotrm2pfforce 370 371Flexible Axis Rotation 372~~~~~~~~~~~~~~~~~~~~~~ 373 374As sketched in :numref:`Fig. %s <fig-rotation>` A–B, the rigid body 375behavior of the fixed axis rotation scheme is a drawback for many 376applications. In particular, deformations of the rotation group are 377suppressed when the equilibrium atom positions directly depend on the 378reference positions. To avoid this limitation, 379eqns. :eq:`%s <eqnpotrmpf>` and :eq:`%s <eqnpotrm2pf>` 380will now be generalized towards a “flexible axis” as sketched in 381:numref:`Fig. %s C <fig-rotation>`. This will be achieved by 382subdividing the rotation group into a set of equidistant slabs 383perpendicular to the rotation vector, and by applying a separate 384rotation potential to each of these slabs. 385:numref:`Fig. %s C <fig-rotation>` shows the midplanes of the slabs 386as dotted straight lines and the centers as thick black dots. 387 388To 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 390total potential function :math:`V^\mathrm{flex}` by a Gaussian function 391 392.. math:: g_n(\mathbf{x}_i) = \Gamma \ \mbox{exp} \left( 393 -\frac{\beta_n^2(\mathbf{x}_i)}{2\sigma^2} \right) , 394 :label: eqngaussian 395 396centered at the midplane of the :math:`n`\ th slab. Here :math:`\sigma` 397is the width of the Gaussian function, :math:`\Delta x` the distance 398between adjacent slabs, and 399 400.. math:: \beta_n(\mathbf{x}_i) := \mathbf{x}_i \cdot \hat{\mathbf{v}} - n \, \Delta x \, . 401 :label: eqngaussianpart2 402 403.. _fig-gaussian: 404 405.. figure:: plots/gaussians.* 406 :width: 6.50000cm 407 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. 412 413A most convenient choice is :math:`\sigma = 0.7 \Delta x` and 414 415.. math:: 1/\Gamma = \sum_{n \in Z} 416 \mbox{exp} 417 \left(-\frac{(n - \frac{1}{4})^2}{2\cdot 0.7^2}\right) 418 \approx 1.75464 \, , 419 :label: eqngaussianpart3 420 421which yields a nearly constant sum, essentially independent of 422:math:`\mathbf{x}_i` (dashed line in 423:numref:`Fig. %s <fig-gaussian>`), *i.e.*, 424 425.. math:: \sum_{n \in Z} g_n(\mathbf{x}_i) = 1 + \epsilon(\mathbf{x}_i) \, , 426 :label: eqnnormal 427 428with 429:math:`| \epsilon(\mathbf{x}_i) | < 1.3\cdot 10^{-4}`. 430This choice also implies that the individual contributions to the force 431from the slabs add up to unity such that no further normalization is 432required. 433 434To each slab center :math:`\mathbf{x}_c^n`, all atoms 435contribute by their Gaussian-weighted (optionally also mass-weighted) 436position vectors 437:math:`g_n(\mathbf{x}_i) \, \mathbf{x}_i`. 438The instantaneous slab centers :math:`\mathbf{x}_c^n` are 439calculated from the current positions 440:math:`\mathbf{x}_i`, 441 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} \, ,\\ 445 :label: eqndefx0 446 447while the reference centers :math:`\mathbf{y}_c^n` are 448calculated from the reference positions 449:math:`\mathbf{y}_i^0`, 450 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} \, . 454 :label: eqndefy0 455 456Due to the rapid decay of :math:`g_n`, each slab will essentially 457involve contributions from atoms located within :math:`\approx 4583\Delta x` from the slab center only. 459 460Flexible Axis Potential 461^^^^^^^^^^^^^^^^^^^^^^^ 462 463We consider two flexible axis variants. For the first variant, the slab 464segmentation procedure with Gaussian weighting is applied to the radial 465motion potential 466(eqn. :eq:`%s <eqnpotrmpf>` / :numref:`Fig. %s B <fig-equipotential>`), 467yielding as the contribution of slab :math:`n` 468 469.. math:: V^n = 470 \frac{k}{2} \sum_{i=1}^{N} w_i \, g_n(\mathbf{x}_i) 471 \left[ 472 \mathbf{q}_i^n 473 \cdot 474 (\mathbf{x}_i - \mathbf{x}_c^n) 475 \right]^2 , 476 :label: eqnflexpot 477 478and a total potential function 479 480.. math:: V^\mathrm{flex} = \sum_n V^n \, . 481 :label: eqnpotflex 482 483Note that the global center of mass :math:`\mathbf{x}_c` 484used in eqn. :eq:`%s <eqnpotrmpf>` is now replaced by 485:math:`\mathbf{x}_c^n`, the center of mass of the slab. 486With 487 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 494 495the resulting force on atom :math:`j` reads 496 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 . 507 \end{aligned} 508 :label: eqnpotflexforce 509 510Note that for :math:`V^\mathrm{flex}`, as defined, the slabs are fixed in 511space and so are the reference centers 512:math:`\mathbf{y}_c^n`. If during the simulation the 513rotation group moves too far in :math:`\mathbf{v}` 514direction, it may enter a region where – due to the lack of nearby 515reference positions – no reference slab centers are defined, rendering 516the potential evaluation impossible. We therefore have included a 517slightly modified version of this potential that avoids this problem by 518attaching the midplane of slab :math:`n=0` to the center of mass of the 519rotation group, yielding slabs that move with the rotation group. This 520is achieved by subtracting the center of mass 521:math:`\mathbf{x}_c` of the group from the positions, 522 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 \, , 525 :label: eqntrafo 526 527such that 528 529.. math:: \begin{aligned} 530 V^\mathrm{flex-t} 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) \| } 536 \cdot 537 (\tilde{\mathbf{x}}_i - \tilde{\mathbf{x}}_c^n) 538 \right]^2 . 539 \end{aligned} 540 :label: eqnpotflext 541 542To simplify the force derivation, and for efficiency reasons, we here 543assume :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`. 546The 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 548tolerated. With this assumption, the forces :math:`\mathbf{F}^\mathrm{flex-t}` 549have the same form as eqn. :eq:`%s <eqnpotflexforce>`. 550 551Flexible Axis 2 Alternative Potential 552^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 553 554In this second variant, slab segmentation is applied to 555:math:`V^\mathrm{rm2}` (eqn. :eq:`%s <eqnpotrm2pf>`), resulting in 556a flexible axis potential without radial force contributions 557(:numref:`Fig. %s C <fig-equipotential>`), 558 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 + 564 \epsilon'} 565 :label: eqnpotflex2 566 567With 568 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}\\ 576 \mathbf{S}^n & := & 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] 582 \end{aligned} 583 :label: eqnSn 584 585the force on atom :math:`j` reads 586 587.. math:: \begin{aligned} 588 \nonumber 589 \mathbf{F}_{\!j}{^\mathrm{flex2}}& = & 590 - k\; 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}} \\ 596 \nonumber 597 & & 598 + k \left\lbrace \sum_n W_{\!j}^n \, \mathbf{S}^n \right\rbrace \times 599 \hat{\mathbf{v}} 600 - k \left\lbrace \sum_n W_{\!j}^n \; \frac{\beta_n(\mathbf{x}_j)}{\sigma^2} \frac{1}{\psi_j}\;\; 601 \mathbf{s}_j^n \cdot 602 \mathbf{S}^n \right\rbrace \hat{\mathbf{v}}\\ 603 & & 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 607 \hat{\mathbf{v}} . 608 \end{aligned} 609 :label: eqnpotflex2force 610 611Applying 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` 617are small, the resulting equations for :math:`V{^\mathrm{flex2 - t}}` 618and :math:`\mathbf{F}{^\mathrm{flex2 - t}}` are similar 619to those of :math:`V^\mathrm{flex2}` and 620:math:`\mathbf{F}^\mathrm{flex2}`. 621 622Usage 623~~~~~ 624 625To apply enforced rotation, the particles :math:`i` that are to be 626subjected 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 631file is found, :math:`\mathbf{x}_i(t=0)` are used as 632reference positions and written to :ref:`trr` such that they 633can be used for subsequent setups. All parameters of the potentials such 634as :math:`k`, :math:`\epsilon'`, etc. 635(:numref:`Table %s <tab-vars>`) are provided as :ref:`mdp` 636parameters; ``rot-type`` selects the type of the potential. 637The option ``rot-massw`` allows to choose whether or not to 638use mass-weighted averaging. For the flexible potentials, a cutoff value 639:math:`g_n^\mathrm{min}` (typically :math:`g_n^\mathrm{min}=0.001`) 640makes 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 644written to additional output files and which are described below. 645 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` 680 681.. _tab-vars: 682 683.. table:: Parameters used by the various rotation potentials. 684 |BFX| indicate which parameter is actually used for a given potential 685 :widths: auto 686 :align: center 687 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 +-------------------+------------+---------+---------+--------+--------+--------+--------+-----------+-----------+ 723 724 725| 726 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` 740 741.. _tab-quantities: 742 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. 745 :widths: auto 746 :align: center 747 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 +------------+---------+------------+--------------------+-------+----------+ 765 766 767 768 769Angle of Rotation Groups: Fixed Axis 770^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 771 772For fixed axis rotation, the average angle :math:`\theta_\mathrm{av}(t)` 773of the group relative to the reference group is determined via the 774distance-weighted angular deviation of all rotation group atoms from 775their reference positions, 776 777.. math:: \theta_\mathrm{av} = \left. \sum_{i=1}^{N} r_i \ \theta_i \right/ \sum_{i=1}^N r_i \ . 778 :label: eqnavangle 779 780Here, :math:`r_i` is the distance of the reference position to the 781rotation axis, and the difference angles :math:`\theta_i` are determined 782from the atomic positions, projected onto a plane perpendicular to the 783rotation axis through pivot point :math:`\mathbf{u}` (see 784eqn. :eq:`%s <eqnproject>` for the definition of 785:math:`\perp`), 786 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 790 \| } \ . 791 :label: eqnavanglepart2 792 793The sign of :math:`\theta_\mathrm{av}` is chosen such that 794:math:`\theta_\mathrm{av} > 0` if the actual structure rotates ahead of 795the reference. 796 797Angle of Rotation Groups: Flexible Axis 798^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 799 800For flexible axis rotation, two outputs are provided, the angle of the 801entire rotation group, and separate angles for the segments in the 802slabs. The angle of the entire rotation group is determined by an RMSD 803fit 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 806be rotated around :math:`\hat{\mathbf{v}}` for the optimal 807fit, 808 809.. math:: \mathrm{RMSD} \big( \mathbf{x}_i,\ \mathbf{\Omega}(\theta_\mathrm{fit}) 810 \mathbf{y}_i^0 \big) \stackrel{!}{=} \mathrm{min} \, . 811 :label: eqnrmsdfit 812 813To determine the local angle for each slab :math:`n`, both reference 814and actual positions are weighted with the Gaussian function of slab 815:math:`n`, and :math:`\theta_\mathrm{fit}(t,n)` is calculated as in 816eqn. :eq:`%s <eqnrmsdfit>` from the Gaussian-weighted 817positions. 818 819For all angles, the :ref:`mdp` input option 820``rot-fit-method`` controls whether a normal RMSD fit is 821performed or whether for the fit each position 822:math:`\mathbf{x}_i` is put at the same distance to the 823rotation axis as its reference counterpart 824:math:`\mathbf{y}_i^0`. In the latter case, the RMSD 825measures only angular differences, not radial ones. 826 827Angle Determination by Searching the Energy Minimum 828^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 829 830Alternatively, for ``rot-fit-method = potential``, the angle 831of the rotation group is determined as the angle for which the rotation 832potential energy is minimal. Therefore, the used rotation potential is 833additionally evaluated for a set of angles around the current reference 834angle. In this case, the ``rotangles.log`` output file 835contains the values of the rotation potential at the chosen set of 836angles, while ``rotation.xvg`` lists the angle with minimal 837potential energy. 838 839Torque 840^^^^^^ 841 842The torque :math:`\mathbf{\tau}(t)` exerted by the 843rotation potential is calculated for fixed axis rotation via 844 845.. math:: \mathbf{\tau}(t) = \sum_{i=1}^{N} \mathbf{r}_i(t) \times \mathbf{f}_{\!i}^\perp(t) , 846 :label: eqntorque 847 848where :math:`\mathbf{r}_i(t)` is the distance vector from 849the rotation axis to :math:`\mathbf{x}_i(t)` and 850:math:`\mathbf{f}_{\!i}^\perp(t)` is the force component 851perpendicular to :math:`\mathbf{r}_i(t)` and 852:math:`\hat{\mathbf{v}}`. For flexible axis rotation, 853torques :math:`\mathbf{\tau}_{\!n}` are calculated for 854each slab using the local rotation axis of the slab and the 855Gaussian-weighted positions. 856