Optimizers

MODELLER currently implements a Beale restart conjugate gradients algorithm [Shanno & Phua, 1980,Shanno & Phua, 1982] and a molecular dynamics procedure with the leap-frog Verlet integrator [Verlet, 1967]. The conjugate gradients optimizer is usually used in combination with the variable target function method [Braun & Gõ, 1985] which is implemented with the AutoModel class (Section A.4). The molecular dynamics procedure can be used in a simulated annealing protocol that is also implemented with the AutoModel class.

Molecular dynamics

Force in MODELLER is obtained by equating the objective function $F$ with internal energy in kcal/mole. The atomic masses are all set to that of C¹² (MODELLER unit is kg/mole). The initial velocities at a given temperature are obtained from a Gaussian random number generator with a mean and standard deviation of:

$$\bar{v}_x =$$ 0 (A.7)

$$\sigma_x = \sqrt{\frac{k_B T}{m_i}}$$ (A.8)

where $k_B$ is the Boltzmann constant, $m_i$ is the mass of one C¹² atom, and the velocity is expressed in angstroms/femtosecond.

The Newtonian equations of motion are integrated by the leap-frog Verlet algorithm [Verlet, 1967]:

$$\dot{\vec{r}}_i\left(t+\frac{\delta t}{2}\right) = \dot{\vec{r}}_i\left(t-\frac{\delta t}{2}\right) - \frac{\partial F}{\partial \vec{r}_i(t)}\frac{\delta t}{m_i}$$ (A.9)

$$\vec{r}_i\left(t+\delta t\right) = \vec{r}_i(t) + \dot{\vec{r}}_i\left(t+\frac{\delta t}{2}\right)\delta t$$ (A.10)

where $\vec{r}_i$ is the position of atom $i$. In addition, velocity is capped at a maximum value, before calculating the shift, such that the maximal shift along one axis can only be cap_atom_shift. The velocities can be equilibrated every equilibrate steps to stabilize temperature. This is achieved by scaling the velocities with a factor $f$:

$$f = \sqrt{T/T_{kin}}$$ (A.11)

$$T_{kin} = \frac{E_{kin}}{\frac{1}{2} k_{b} N_{f}}$$ (A.12)

$$E_{kin} = \frac{1}{2} \sum_i^{N_{atoms}} m_i\dot{\vec{r}}_i^2$$ (A.13)

where $k_B$ is the Boltzmann constant, $N_f$ the number of degrees of freedom, $E_{kin}$ the current kinetic energy and $T_{kin}$ the current kinetic temperature.

Langevin dynamics

Langevin dynamics (LD) are implemented as in [Loncharich et al., 1992]. The equations of motion (Equation A.9) are modified as follows:

$$\dot{\vec{r}}_i\left(t+\frac{\delta t}{2}\right) = \dot{\vec{r}}_... ...\vec{r}_i(t)}\right) \frac{\delta t}{m_i} \frac{1}{1+\frac{1}{2}\gamma\delta t}$$ (A.14)

where $\gamma$ is a friction factor (in $fs^{-1}$) and $\vec{R}_i$ a random force, chosen to have zero mean and standard deviation

$$\sigma(\vec{R}_i) = \sqrt{\frac{2\gamma m_i k_{B} T}{\delta t}}$$ (A.15)

Self-guided MD and LD

MODELLER also implements the self-guided MD [Wu & Wang, 1999] and LD [Wu & Brooks, 2003] methods. For self-guided MD, the equations of motion (Equation A.9) are modified as follows:

$$\vec{g}_i(t) = \left(1 - \frac{\delta t}{t_l}\right) \vec{g}_i(t - \delta t) + \... ...\lambda \vec{g}_i(t-\delta t) - \frac{\partial F}{\partial \vec{r}_i(t)}\right)$$ (A.16)

$$\dot{\vec{r}}_i\left(t+\frac{\delta t}{2}\right) = \dot{\vec{r}}_i\left(t-\frac{\delta t}{2}\right) + \left(\lambda \vec{g}_i(t) - \frac{\partial F}{\partial \vec{r}_i(t)}\right)\frac{\delta t}{m_i}$$ (A.17)

where $\lambda$ is the guiding factor (the same for all atoms), $t_l$ the guide time in femtoseconds, and $\vec{g}_i$ a guiding force, set to zero at the start of the simulation. (Position $\vec{r}_i$ is updated in the usual way.)

For self-guided Langevin dynamics, the guiding forces are determined as follows (terms are as defined in Equation A.14):

$$\vec{g}_i(t) = \left(1 - \frac{\delta t}{t_l}\right) \vec{g}_i(t ... ...\frac{\delta t}{t_l}\gamma m_i \dot{\vec{r}}_i\left(t-\frac{\delta t}{2}\right)$$ (A.18)

A scaling parameter χ is then determined by first making an unconstrained half step:

$$\dot{\vec{r}}_i'(t) = \dot{\vec{r}}_i\left(t - \frac{\delta t}{2}\right) + \frac{1}{2}\... ... + \vec{R}_i - \frac{\delta F}{\delta \vec{r}_i(t)}\right) \frac{\delta t}{m_i}$$ (A.19)

$$\zeta = \left(1+\frac{\gamma\delta t}{2}\right)\frac{\sum_i^N \lambda \vec{g}_i(t)\dot{\vec{r}}_i'(t)}{\sum_i^N m_i \dot{\vec{r}}_i'^2(t)}$$ (A.20)

$$\chi = \left(1 + \frac{(\gamma + \zeta)\delta t}{2}\right)^{-1}$$ (A.21)

Finally, the velocities are advanced using the scaling factor:

$$\dot{\vec{r}}_i\left(t+\frac{\delta t}{2}\right) = (2\chi - 1)\do... ...vec{R}_i - \frac{\partial F}{\partial \vec{r}_i(t)}\right) \frac{\delta t}{m_i}$$ (A.22)

Rigid bodies

Where rigid bodies are used, these are optimized separately from the other atoms in the system. This has the additional advantage of reducing the number of degrees of freedom.

Rigid molecular dynamics

The state of each rigid body is specified by the position of the center of mass, $\vec{r}_{COM}$, and an orientation quaternion, $\tilde{q}$ [Goldstein, 1980]. (The quaternion has 4 components, $q_1$ through $q_4$, of which the first three refer to the vector part, and the last to the scalar.) The translational and rotational motions of each body are separated. Each body is translated about its center of mass using the standard Verlet equations (Equation A.9) using the force:

$$\frac{\partial{F}}{\partial{\vec{r}_{COM}}} = \sum_i \frac{\partial{F}}{\partial{\vec{r_i}}}$$ (A.23)

where the sum $i$ operates over all atoms in the rigid body, and $\vec{r_i}$ is the position of atom $i$ in real space.

For the rotational motion, the orientation quaternions are again integrated using the same Verlet equations. For this, the quaternion accelerations are calculated using the following relation [Rapaport, 1997]:

$$\ddot{\tilde{q}} = \frac{1}{2} \mathsfsl{W}^T\left( \begin{ar... ...dot{\omega}'_z \\ -2\sum_m \dot{q}_m^2 \\ \end{array} \right)$$ (A.24)

where $\mathsfsl{W}$ is the orthogonal matrix

$$\mathsfsl{W} = \left(\begin{array}{rrrr} q_4 & q_3 & -q_2 & -q_1 ... ..._2 \\ q_2 & -q_1 & q_4 & -q_3 \\ q_1 & q_2 & q_3 & q_4 \\ \end{array}\right)$$ (A.25)

and $\dot{\omega}'_k$ is the first derivative of the angular velocity (in the body-fixed frame) about axis $k$ - i.e., the angular acceleration. These angular accelerations are in turn calculated from the Euler equations for rigid body rotation, such as:

$$\dot{\omega}'_x = \frac{T_x + (I_y - I_z)\omega'_y\omega'_z}{I_x}$$ (A.26)

(Similar equations exist for the $y$ and $z$ components.) The angular velocities $\vec{\omega}'$ are obtained from the quaternion velocities:

$$\left(\begin{array}{c} \omega'_x \\ \omega'_y \\ \omega'_z \\ 0 \\ \end{array}\right) = 2\mathsfsl{W}\dot{\tilde{q}}$$ (A.27)

The torque, $\vec{T}$, in the body-fixed frame, is calculated as

$$\vec{T} = \mathsfsl{A} \sum_i (\vec{r_i} - \vec{r}_{COM}) \times -\frac{\partial{F}}{\partial{\vec{r_i}}}$$ (A.28)

and $\mathsfsl{A}$ is the rotation matrix to convert from world space to body space

$$\mathsfsl{A} = 2 \left( \begin{array}{ccc} q_1^2 + q_4^2 - \frac{... ...+ q_2q_4 & q_2q_3 - q_1q_4 & q_3^2 + q_4^2 - \frac{1}{2} \\ \end{array}\right)$$ (A.29)

and finally the $x$ component of the inertia tensor, $I_x$, is given by

$$I_x = \sum_i m_i ({r'}_{i,y}^2 + {r'}_{i,z}^2)$$ (A.30)

where $\vec{r'_i}$ is the position of each atom in body space (i.e. relative to the center of mass, and unrotated), and $m_i$ is the mass of atom $i$ (taken to be the mass of one $C^{12}$ atom, as above). Similar relations exist for the $y$ and $z$ components.

The kinetic energy of each rigid body (used for temperature control) is given as a combination of translation and rotational components:

$$E_{kin}^{body} = \frac{1}{2}(\sum_i m) \dot{\vec{r}}_{COM}^{2} + \frac{1}{2}(I_x {\omega'}_x^2 + I_y {\omega'}_y^2 + I_z {\omega'}_z^2)$$ (A.31)

Initial translational and rotational velocities of each rigid body are set in the same way as for atomistic dynamics.

Rigid minimization

The state of each rigid body is specified by 6 parameters: the position of the center of mass, $\vec{r}_{COM}$, and the rotations in radians about the body-fixed axes: $\theta_x$, $\theta_y$, and $\theta_z$. The first derivative of the objective function $F$ with respect to the center of mass is obtained from Equation A.23, and those with respect to the angles from:

$$\frac{\partial{F}}{\partial{\theta_k}} = \mathsfsl{M_k}\vec{r'_i} \cdot \frac{\partial{F}}{\partial{\vec{r_i}}}$$ (A.32)

The transformation matrices $\mathsfsl{M_k}$ are given as:

$$\mathsfsl{M_x} = \left[ \begin{array}{ccc} 0 & -\sin{\theta_z}\si... ...{\theta_y}\cos{\theta_x} & -\cos{\theta_y}\sin{\theta_x} \\ \end{array}\right]$$ (A.33)

$$\mathsfsl{M_y} = \left[ \begin{array}{ccc} -\cos{\theta_z}\sin{\t... ...{\theta_y}\sin{\theta_x} & -\sin{\theta_y}\cos{\theta_x} \\ \end{array}\right]$$ (A.34)

$$\mathsfsl{M_z} = \left[ \begin{array}{ccc} -\sin{\theta_z}\cos{\t... ...\cos{\theta_z}\sin{\theta_y}\cos{\theta_x} \\ 0 & 0 & 0 \\ \end{array}\right]$$ (A.35)

The atomic positions $\vec{r_i}$ are reconstructed when necessary from the body's orientation by means of the following relation:

$$\vec{r_i} = \mathsfsl{M}\vec{r'_i} + \vec{r}_{COM}$$ (A.36)

where $\mathsfsl{M}$ is the rotation matrix

$$\mathsfsl{M} = \left[ \begin{array}{ccc} \cos{\theta_z}\cos{\thet... ...s{\theta_y}\sin{\theta_x} & \cos{\theta_y}\cos{\theta_x} \\ \end{array}\right]$$ (A.37)