This section describes the optimization methods implemented in
\MODELLER. The general form of the objective function and the structure of
optimization are similar to molecular dynamics programs, such as
\CHARMM\ \cite{CHARMM98}.

\subsection{Function}

\MODELLER\ minimizes the {\it objective function} $F$ with respect to 
Cartesian coordinates of $\sim 10,000$ atoms (3D points) that form a
{\it system} (one or more molecules):
\begin{equation}
 F = F({\bf R}) = F_{symm} + \sum_i c_i ({\bf f}_i, {\bf p}_i)
 \label{eq:defobjfunc}
\end{equation}
where $F_{symm}$ is an optional symmetry term defined in 
Eq.~\ref{eq:symmetry},
$\bf R$ are Cartesian coordinates of all atoms, $c$ is a restraint
$i$, ${\bf f}$ is a geometric feature of a molecule, and ${\bf p}$ are
parameters. For a 10,000 atom system there can be on the 
order of 200,000 restraints.  The form of
$c$ is simple; it includes a quadratic function, cosine, a weighted
sum of a few Gaussian functions, Coulomb law, Lennard-Jones potential,
cubic splines, and some other simple functions.  The geometric features 
presently include a distance, an angle, a dihedral angle, a pair of dihedral
angles between two, three, four atoms and eight atoms, respectively,
the shortest distance in the set of distances (not documented
further), solvent accessibility in $\Ang^2$, and atom density 
expressed as the number of atoms around the central atom.
A pair of dihedral angles can be used to restrain such strongly
correlated features as the mainchain dihedral angles $\Phi$ and
$\Psi$. Each of the restraints also depends on a few parameters ${\bf
p}_i$ that generally vary from a restraint to a restraint.  Some
restraints can restrain {\it pseudo-atoms} such as a gravity center of
several atoms.

\MODELLER\ allows some atoms to be {\it fixed} during optimization; 
\ie, only selected atoms are allowed to be moved. Similarly, \MODELLER\ 
also allows only a subset of all restraints to be actually used in the
calculation of the objective function. Each subset is indicated by a
list of indices specifying the selected atoms or restraints.

There are two kinds of restraints, {\it static} and {\it dynamic},
that both contribute to the objective function as indicated in
Eq.~\ref{eq:defobjfunc}:
\begin{equation}
    F = F_{symm} + F_s + F_d \; .
\end{equation}
The static restraints and their parameters are pre-defined; \ie, they
are given before the call to the optimizer and are not changed during
optimization. The dynamic restraints are re-generated repeatedly
during optimization. Usually, the CPU time is spent evenly between the
two kinds of restraints, although the dynamic restraints become more
important as the size of the system increases. All dynamic restraints
are always selected and they can restrain only pairs of atoms. In all
other respects, the two kinds of restraints are the same.

The dynamic restraints are obtained from a {\it dynamic pairs list}
(the non-bonded pairs list). Each dynamic pair corresponds to at least
one restraint, which may or may not be violated. The dynamic pairs
list includes only the pairs of atoms that satisfy the following three
conditions: (1) One or both atoms in a pair are allowed to move.  (2)
The two atoms are not connected through one, two, or three chemical bonds.
(3) The two atoms are closer than a preset cutoff distance (\eg, 4
\AA). There are on the order of 5000 atom pairs in the dynamic pairs
list when only soft-sphere overlap restraints are used. Currently, the 
restraint types on the dynamic atom pairs that can be selected include
the soft-sphere overlap, Lennard-Jones, Coulombic interactions, and
\MODELLER\ non-bonded spline restraints. xx atom density?

The existence of the dynamic pairs list is justified by the fact that 
dynamic pairs are usually a small fraction of all possible
atom--atom pairs ($N \cdot (N-1)/2$, where $N$ is the number of atoms
in a system). The use of the dynamic pairs list becomes especially
beneficent as the size of the system increases.

The actual algorithm for creating the dynamic pairs list varies with
the size of the system, whether or not all atoms are allowed to
move, or whether or not the user wants to include the fixed environment
in the calculation of non-bonded restraints involving the selected
atoms. See Section~\ref{COMMAND:optimize} for more information.

% \footnote{There is some code for Monte
% Carlo optimization but it is currently not compiled into the
% executable.  The Monte Carlo optimizer implements a scheme that allows
% a faster calculation of the objective function upon a change in the
% position of a single atom. Note that this situation is different from
% the conjugate gradients and molecular dynamics optimizations where the
% objective function is re-evaluated only after {\it all} atoms have
% moved. For Monte Carlo optimization when $N_s = 1$, another
% specialized search over all other atoms is used. Should be changed to
% a search over all other non-fixed atoms only, because both atoms in a
% dynamic pair have to be non-fixed.} 

The hash-function algorithm is used to determine whether or not two atoms 
are a dynamic atom pair. This algorithm is about 20 times slower than
a lookup table but it requires much less memory and still spends a
negligible fraction of the total CPU time. A hash-function table is
prepared only once before the start of the optimization and any other operation
involving an evaluation of the objective function (\eg, \C{OPTIMIZE},
\C{ENERGY}, and \C{PICK\_HOT\_ATOMS}).

The dynamic pairs list is not necessarily re-generated each time the
objective function is evaluated, although the contribution of the
restraint to the objective function is calculated in each call to the
objective function routine with the current values of the Cartesian
coordinates. The dynamic pairs list is
re-generated only when maximal atomic shifts accumulate to a value
larger than a preset cutoff. This cutoff is chosen such that there
cannot be a violation of a restraint without having its atom pair on
the dynamic pairs list. The dynamic pairs list is recalculated in
$\sim 20\%$ and $\sim 2\%$ of the objective function calls at the
beginning and the end of optimization, respectively.

Each evaluation of the objective function or of its first derivatives
with respect to the Cartesian coordinates involves the following
steps:

\begin{enumerate}
\item Calculate non-fixed pseudo-atoms from the current atomic positions \\ 
      (routine \V{objfunc:pseudo}).

\item Update the dynamic pairs list, if necessary
      (routine \V{objfunc:upddyn}).

\item Calculate the violations of selected restraints and all other 
quantities that are shared between the calculations of the objective
function and its derivatives (routine \V{objfunc:getviol}).

\item Sum the contributions of all violated restraints to the
objective function and the derivatives (routine \V{objfunc:\-getviol}).
\end{enumerate}


\subsection{Optimizers}

\MODELLER\ currently implements a Beale restart conjugate gradients 
algorithm \cite{ShaPhu80,ACM500}
and a molecular dynamics procedure with the Verlet integrator \cite{Ver67}.
The conjugate gradients optimizer is usually
used in combination with the variable target function method 
\cite{BraGo85} which is
implemented with the \TOP\ script (Section~\ref{SECTION:flowchart}). 
The molecular
dynamics procedure can be used in a simulated annealing protocol that
is also implemented with the \TOP\ script.



\subsubsection{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:
\begin{eqnarray}
\bar{v}_x & = & 0 \\
\sigma_x  & = & \sqrt{\frac{k_B T}{m}} = 0.000263143 \sqrt{T}
\end{eqnarray}
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 \cite{Ver67}:
\begin{eqnarray}
v_x(i+1) & = & v_x(i) + \frac{\partial F}{\partial x} A \\
x(i+1)   & = & x(i) + v_x(i+1) \Delta t \\
A        & = & c \frac{\Delta t}{m} = 4.1868 \cdot 10^{-7} \frac{\Delta t}{m}
\end{eqnarray}
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 \K{CAP\_ATOM\_SHIFT}. The
velocities can be equilibrated every \K{EQUILIBRATE} steps to
stabilize temperature.  This is achieved by scaling the velocities
with a factor $f$:
\begin{eqnarray}
f       & = & \sqrt{T/E_{kin}} \\
E_{kin} & = & \frac{m}{2} \sum_i^{N_{atoms}} (v_x^2+v_y^2+v_z^2)
\end{eqnarray}
where $E_{kin}$ is the current kinetic energy of the system.
