Optimizers

MODELLER currently implements a Beale restart conjugate gradients algorithm [Shanno & Phua, 1980,Shanno & Phua, 1982] and a molecular dynamics procedure with the 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 6.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$^{12}$ (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$$ (6.7)

$$\sigma_x = \sqrt{\frac{k_B T}{m}} = 0.000263143 \sqrt{T}$$ (6.8)

where $k_B$ is the Boltzmann constant, $m$ is the mass of one C$^{12}$ atom, and the velocity is expressed in angstroms/femtosecond.

The Newton's equations of motion are integrated by the Verlet algorithm [Verlet, 1967]:

$$v_x(i+1) = v_x(i) + \frac{\partial F}{\partial x} A$$ (6.9)

$$x(i+1) = x(i) + v_x(i+1) \Delta t$$ (6.10)

$$A = c \frac{\Delta t}{m} = 4.1868 \cdot 10^{-7} \frac{\Delta t}{m}$$ (6.11)

where velocities $v(i+1)$ are for $t + \Delta t / 2$ and positions $x(i+1)$ for $t + \Delta t$. Parameter $c$ is a scaling factor so that positions are expressed in angstroms, time in femtoseconds, and velocities in angstroms/femtosecond, given that the objective function is in kcal/mole and atomic mass in kg/mole. 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/E_{kin}}$$ (6.12)

$$E_{kin} = \frac{m}{2} \sum_i^{N_{atoms}} (v_x^2+v_y^2+v_z^2)$$ (6.13)

where $E_{kin}$ is the current kinetic energy of the system.