Iteration

After initialization, the center-of-mass variables are integrated using the three-step Verlet method:

$$R_J(t+dt) = -R_J(t-dt)+2R_J(t)+f_J(t){(dt)^2 \over M_J}$$ (11.24)

$$V_J(t) = ({1\over 2dt})(R_J(t+dt)-R_J(t-dt))$$ (11.25)

Friction terms are embedded in the force $f_J$.

The integration of rotational variables uses a different central difference expansion. The algorithm begins with the advancement of the Euler-Cayley parameters by ${1\over 2}dt$, using the following four steps:

$$J(t) = J(t-{1\over 2}dt)+{1\over 2}N(t)$$ (11.26)

$$J^B(t) = {\bf A}(t)J(t)$$ (11.27)

$$\omega^B_{\alpha} = J^B_{\alpha}(t)/I_\alpha \hspace*{1.0in}\alpha \in (x,y,z)$$ (11.28)

$$q(t+{1\over 2}dt) = q(t)+{1\over 2}{\bf Q}(t)\omega'(t)dt$$ (11.29)

Here $J^B$ represents the angular momentum of a rigid body transformed to the body frame, and $\omega^B$ represents the transformed angular velocity; q is the four-vector ($q_0$,$q_1$,$q_2$,$q_3$). Finally, the torques are determined by

$$N_J=\sum_{i\in J} (r_i-R_J)\times f_i$$ (11.30)

Again, all frictional terms are embedded in the force $f_J$. These half-step calculations are necessary for the final, full-step update:

$$J(t+{1\over 2}dt) = J(t-{1\over 2}dt)+N(t)dt$$ (11.31)

$$J^B(t+{1\over 2}dt) = {\bf A}(t+{1\over 2}dt)J(t+{1\over 2}dt)$$ (11.32)

$$\omega^B_{\alpha}(t+{1\over 2}dt) = J^B_{\alpha}(t+{1\over 2}dt)/I_\alpha \hspace{0.5in} \alpha \in (x,y,z)$$ (11.33)

$$q(t+dt) = q(t)+{\bf Q}(t+{1\over 2}dt)\omega'(t+{1\over 2}dt)dt$$ (11.34)

The velocities of individual atoms, which are necessary for the friction evaluation, are determined using the angular velocities of the rigid bodies. The velocity due to rotation alone is determined in the body frame, then transformed into the lab frame and added to the translational velocity of the center of mass:

$$v_i(t)=V_J(t)+{\bf A_J}^{-1}\left[\begin{array}{c} z^B_{i}\o... ...\ y^B_{i}\omega^B_{x}-x^B_{i}\omega^B_{y} \end{array} \right]$$ (11.35)

Finally, the atom positions are updated for the next energy evaluation using the new center-of-mass positions and the new Euler-Cayley parameters:

$$r_i(t+dt)=R_J(t+dt)+{\bf A^{-1}_{J}}r^B_{i}$$ (11.36)

The rigid-body dynamics routine can be run in one of three modes: FREE, LANGevin, or TCOUpling. The FREE mode simply allows frictionless dynamics, the LANGevin mode adds friction- and bath-dependent random forces, and TCOUpling uses the method of Berendsen et al. (1984). The coordinates (and velocities) are affected by using a nonzero friction coefficient in the Langevin dynamics algorithm with zero random forces. The friction coefficient is computed by the program from the equation

$$b_i=b_i^I (T_o/T -1).$$ (11.37)

where $b_i^I$ is given by the FBETa atom property. Note that FBETa is by default zero; that is, FBETa has to be specified in order to use TCOUpling. The target temperature $T_o$ is specified by the TBATh parameter.

If the simulation involves light rigid groups such as water molecules, it is recommended that a time step on the order of .25 fsec be used. For heavier groups, time steps as large as 20 or 30 fsec are adequate.