This section describes commands for creating, reading and
writing optimization schedule, and for calculating and optimizing
the objective function. For technical background, see
Section~\ref{SECTION:energydesc}.


\Command      {MAKE\_SCHEDULE}{create optimization schedule}

\Options{
\OptLine{LIBRARY\_SCHEDULE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{SCHEDULE\_SCALE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
}

\Requirements {MODEL}

\Output       {N\_SCHEDULE}

\Description  {This command constructs an optimization schedule for
the variable target function method for the current MODEL. 

The template for construction of the schedule is the
\K{LIBRARY\_SCHEDULE}-th entry in library file {\tt \$SCHED\_\-LIB}. 

The usual schedule for the variable target function part of
optimization in comparative modeling is as follows. The residue range
(\C{PICK\_RESTRAINTS} and Section~\ref{SECTION:span})
is increased with increasingly larger steps until the protein length
is reached. The scaling of homology-derived and bonded stereochemical
restraints increases from a small value to 1 in the initial few steps
to allow for imperfect starting geometries, especially those that
result from \C{RANDOMIZE\_XYZ} and long insertions or deletions. The
soft-sphere overlap restraints are slowly introduced only in the last
four steps of the variable target function method to save CPU time and
increase the radius of convergence. In comparative modeling by the
\Z{model} script in the default mode, the variable target function 
method is usually followed by simulated annealing with molecular 
dynamics. In this last stage, all homology-derived and 
stereochemical restraints are generally used with the scaling factors
of 1. There are a number of variables defined in the \Z{modlib/\_\_defs.top}
script that can be used to influence the thoroughness of both the variable
target function and molecular dynamics parts of the optimization
(Chapter~\ref{CHAPTERSCRIPTS}).

The scaling factors for all physical restraint groups, in all schedule
steps, are multiplied by the corresponding scalar in \K{SCHEDULE\_SCALE}
(1 by default). This is useful when template-derived fold restraints
have to be weakened relative to some external restraints, so that the
fold can actually reflect these external restraints, even when they
are quite different from the template-derived restraints.

This command is an alternative to the \C{READ\_SCHEDULE} command.  

Use the \C{WRITE\_SCHEDULE} command to find out what the calculated
schedule is. The schedule file written by the {\tt `model'} routine
has an extension {\tt .sch}.}


\Example{../examples/commands/make_schedule.top}




\Command      {READ\_SCHEDULE}{read optimization schedule}

\Options{
\OptLine{FILE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DIRECTORY}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{SCHEDULE\_SCALE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
}

\Output       {N\_SCHEDULE}

\Description  {This command reads a text file that contains an 
optimization schedule for the variable target function method. 

Each line in the file contains in free format the parameters for a
single step of the variable target function method. These parameters
are: step index (not used by the program), optimization method, maximal 
difference in residue indices of atoms restrained by the selected 
restraints (\C{PICK\_RESTRAINTS} and Section~\ref{SECTION:span}),
and the scaling factors for all types of restraints. The smaller the 
scaling factor, the weaker the corresponding restraint.

See \C{MAKE\_SCHEDULE} for explanation of \K{SCHEDULE\_SCALE}.

This command also sets the \TOP\ variable \K{N\_SCHEDULE} to the total
number of the variable target function steps that were read in.
}

\SeeExample{\C{MAKE\_SCHEDULE}}



\Command      {WRITE\_SCHEDULE}{write optimization schedule}

\Options{
\OptLine{FILE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{OUTPUT\_DIRECTORY}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
}

\Description  {This command writes out the schedule for the variable 
target function method. This schedule file can then be read by the
\C{READ\_SCHEDULE} command.
}

\SeeExample{\C{MAKE\_SCHEDULE}}

\beautypb

\Command      {ENERGY}{evaluate MODEL given restraints}

\Options{
\OptLine{VIOL\_REPORT\_CUT}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{VIOL\_REPORT\_CUT2}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{OUTPUT}{TYPE}{VALUES}{DEFAULT}{\Z{SHORT} \OR\ \Z{LONG} \OR\
         \Z{VERY\_LONG} \OR\ \Z{GRADIENT} \OR\ \Z{SYMMETRY} \OR\
         \Z{ENERGY\_PROFILE} \OR\ \Z{VIOLATIONS\_PROFILE} }
\OptLine{NORMALIZE\_PROFILE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{SMOOTHING\_WINDOW}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{ASGL\_OUTPUT}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{SCHEDULE\_STEP}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{TOPOLOGY\_MODEL}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{RADII\_FACTOR}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{SPHERE\_STDV}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_SPHERE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_LENNARD}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_COULOMB}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_MODELLER}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_ACCESS}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{EXCL\_LOCAL}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{LENNARD\_JONES\_SWITCH}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{COULOMB\_SWITCH}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{RELATIVE\_DIELECTRIC}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{CONTACT\_SHELL}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{UPDATE\_DYNAMIC}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{NLOGN\_USE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{COVALENT\_CYS}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{RESIDUE\_SPAN\_RANGE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
}

\Output {\K{MOLPDF}}

\Requirements {MODEL \AND\ restraints}

\Description{The main purpose of this command is to compare spatial features 
of the current MODEL with the selected restraints in order to
determine the violations of the molecular pdf. It lists variable
amounts of information about the values of the basis, feature, and
molecular pdf's for the current MODEL. All arguments that affect the 
value of the molecular pdf are also relevant for the \C{ENERGY} command.
  
Most of the output goes to the \logf\ file. The output of the \C{ENERGY}
command has to be examined carefully, at least at the end of the
optimization, when the final model is produced.  
Additional output files, for the \ASGL\ plotting program are created if 
\K{ASGL\_OUTPUT} = \V{on} (undocumented). 
% \ASGL\ is not part of the \MODELLER\ distribution, but you can obtain it 
% at URL \htmladdnormallink{{\tt http://guitar/rockefeller.edu/}}
% {http://guitar.rockefeller.edu/}. The \MODELLER\ data files for \ASGL\ 
% include {\tt viol.asgl} with average residue violations for each
% physical restraint type, which can be subsequently used with the
% \ASGL\ script {\tt viol} to plot the average feature violations. There 
% are five other \ASGL\ files that may be produced and can be used 
% subsequently with the \ASGL\ script {\tt chipred}: {\tt chi1pred1.asgl, 
% chi1pred2.asgl, chi2pred1.asgl, chi2pred2.asgl} and {\tt seq.asgl}. 
% The {\tt chipred} script plots the most likely value for $\chi_1$ and 
% $\chi_2$ sidechain dihedral angles and compares them with the actually 
% observed values for each residue of the current MODEL. 
% This option may not yet work as described.
  
\K{OUTPUT} selects various kinds of output information:
\begin{itemize}

\item \Z{LONG} writes restraint violations one per line to the \logf\ file.

\item \Z{VERY\_LONG} writes the most detailed examination of the selected
basis and feature pdf's to the \logf\ file, using several lines of
output for each restraint.

\item \Z{GRADIENT} writes the `force' gradients for the
currently selected restraints to the isotropic temperature
factors for each atom of the current MODEL.  

\item \Z{SYMMETRY} writes a comparison of equivalent distances
involved in the definition of the symmetry enforcing term
to the \logf\ file.
\end{itemize}
  
\K{VIOL\_REPORT\_CUT} is a vector with one real number for each 
physical restraint type. A restraint is reported when its `heavy 
relative violation' is larger than the corresponding cutoff.
The heavy relative violation is calculated by finding
the global minimum of a feature according to the restraint, taking
the difference between the actual feature in the model and this
global minimum, and then normalizing the difference by the standard 
deviation of the global minimum. The `minimal violation' of a restraint
is defined as the difference from the local minimum closest to 
the value of the feature in the model (with the exception of the
spline restraints; see next paragraph).

\K{VIOL\_REPORT\_CUT2} is similar to \K{VIOL\_REPORT\_CUT}, except
that it contains cutoffs for restraint `energies', not heavy
relative violations.
   
The meaning of various other reported properties of the violated restraints 
is briefly described in the \logf\ file. Note that for multi-modal 
restraints that are described by 
cubic splines (by default, all multimodal homology-derived restraints),
only one optimal value is defined, not the local and global minimum
as for the multi-modal Gaussian restraints. As a result, the
minimal violations and heaviest violations are the same.
For interpreting the seriousness of violations, use the following
rule of thumb: There should be at most a few small violations 
(\eg, 4 standard deviations) for all monomodal restraints. In comparative
modeling, the monomodal restraints include the stereochemical restraints 
and distance restraints when only one homologous structure is used. 
For the multimodal restraints, there are usually many violations reported
because the heaviest violations are used in deciding whether or not to
report a violation. In comparative modeling, the multimodal restraints
include the $\chi_i$ restraints, ($\Phi$, $\Psi$) binormal restraints
and distance restraints when more than one template is used. See also
Section~\ref{SECTION:faq}, Question~\ref{FAQ:log}.

For profiles:

This command calculates residue energies or heavy relative violations, 
depending on \K{OUTPUT}, for all physical restraint types (there are 
{\tt NPHYCNS} of them). Relative heavy violations (Table~\ref{tab:restrtypes}) 
are used because only {\it relative} violations of different features are 
comparable.
In both cases, the residue sum is the sum over all restraints that have
at least one atom in a given residue. The contribution of each restraint 
is counted exactly once for each residue, without any weighting. Restraints 
spanning more than one residue contribute equally to all of them. Thus, 
the sum of residue energies is generally larger than
molecular pdf. The command also calculates the sum of the {\tt NPHYCNS}
contributions for each residue and writes all {\tt NPHYCNS+1} columns to 
a file suitable for plotting by \ASGL. 

If \K{NORMALIZE\_PROFILE} is \V{on} the profile for each
residue is normalized by the number of terms applying to each residue.

All the curves are smoothed by the running window
averaging method if \K{SMOOTHING\_WINDOW} is larger than 0: The window
is centered on residue $i$ and extends for (\K{SMOOTHING\_WINDOW}/2) - 1
residues on each side. Thus, \K{SMOOTHING\_WINDOW} has to be an even 
number (or it is made such by the program automatically).
The only exceptions are the 
two terminii, where a smaller number of residues are available for
smoothing. The relative weight of residue $j$ when calculating the
smoothed value at residue $i$ is $(\K{SMOOTHING\_WINDOW}/2 - |j-i|)$.

The energy or the violations profile is written to the fourth column
of the MODEL atomic records (atomic isotropic temperature factors for
X-ray structures). Note that all the atoms in one residue get the
same number. This output is useful for exploring the violations on 
a graphics terminal. 


See description of \C{OPTIMIZE} for the other variables.
}


\Example{../examples/commands/energy.top}






\Command      {OPTIMIZE}{optimize MODEL given restraints \label{COMMAND:optimize}}

\Options{
\OptLine{OPTIMIZATION\_METHOD}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{SCHEDULE\_STEP}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{TOPOLOGY\_MODEL}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{RADII\_FACTOR}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{SPHERE\_STDV}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_SPHERE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_LENNARD}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_COULOMB}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_MODELLER}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DYNAMIC\_ACCESS}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{EXCL\_LOCAL}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{LENNARD\_JONES\_SWITCH}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{COULOMB\_SWITCH}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{RELATIVE\_DIELECTRIC}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{NONBONDED\_SEL\_ATOMS}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{RESIDUE\_SPAN\_RANGE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{COVALENT\_CYS}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{CONTACT\_SHELL}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{UPDATE\_DYNAMIC}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{NLOGN\_USE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{TRACE\_OUTPUT}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{MAX\_ITERATIONS}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{OUTPUT}{TYPE}{VALUES}{DEFAULT}{\Z{NO\_REPORT}}
\bull\ For conjugate gradients: \\
\OptLine{MIN\_ATOM\_SHIFTS}{TYPE}{VALUES}{DEFAULT}{\Z{NO\_REPORT} \OR\ \Z{REPORT}}
% \bull\ For Monte Carlo simulated annealing: &  do not use. \\[0.5cm]
% For Monte Carlo simulated annealing:    & \\
%    SA\_STEP       = \Real{1}            & \\
%    SA\_MVFRACT    = \Real{1}            & \\
%    SA\_TFACTR     = \Real{1}            & \\
%    SA\_T0         = \Real{1}            & \\
%    SA\_TMIN       = \Real{1}            & \\
%    SA\_MOVSPERTMP = \Integer{1}         & \\
%    SA\_MAXTMPCYC  = \Integer{1}         & \\
%    RAND\_SEED     = \Integer{1}         & \\[0.5cm]
\bull\ For molecular dynamics:            \\
\OptLine{MD\_TIME\_STEP}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{INIT\_VELOCITIES}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{TEMPERATURE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{EQUILIBRATE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{MD\_RETURN}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{CAP\_ATOM\_SHIFT}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{RAND\_SEED}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{STOP\_ON\_ERROR}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
}

\Output {\K{MOLPDF}, \K{MODELLER\_STATUS}}

\Requirements{MODEL \AND\ restraints}

\Description{This command performs a number of optimizing iterations using
a selected optimization method (\ref{SECTION:energydesc}). 
One call to \C{OPTIMIZE} corresponds to a single step of the variable
target function method. The whole variable target function method is
implemented by a \TOP\ script. The molecular pdf is optimized with
respect to the selected coordinates of the current MODEL; the
optimized coordinates are returned as the current MODEL. 

Some output may be generated during optimization; for example, a value
of the molecular pdf, average and maximal atomic shifts are
written to the current tracing file every \K{TRACE\_OUTPUT} iterations
of the optimizer if \K{TRACE\_OUTPUT} is larger than 0 (see the
\C{SWITCH\_TRACE} command).  

In addition, a summary of the optimization results is written to the 
\logf\ file after optimization, unless \K{OUTPUT} contains string 
\Z{NO\_REPORT}.

\K{OPTIMIZATION\_METHOD} = 1 selects a conjugate gradients optimization
method. \K{OPTIMIZATION\_METHOD} = 3 selects a molecular dynamics
optimization at a fixed temperature.
The conjugate gradients optimizer is a modified version of the Beale restart
conjugate gradients method \cite{ShaPhu80,ACM500}.
The molecular dynamics routine is the most basic version of the
iterative solver of the Newton's equations of motion. The integrator
uses the Verlet algorithm \cite{Ver67}. All atomic masses are set to that 
of carbon 12. A brief description of the algorithms is given in
Section~\ref{SECTION:energydesc}.

\K{SCHEDULE\_STEP} is the variable target function step. It selects some
of the optimization parameters; it refers to the line in the schedule
file which specifies (1) the optimization method (1=Conjugate Gradients, 
3=Molecular Dynamics); (2) maximal number of residues that the restraints 
are allowed to span (Section~\ref{SECTION:span}); (3) the individual scaling 
factors for all the physical restraint types. \K{OPTIMIZATION\_METHOD} 
overrides the schedule specification if it is within a defined range.
  
\K{CONTACT\_SHELL} defines the maximal distance between atoms that
flags a non-bonded atom pair. Such pairs are stored in the list of
non-bonded atom pairs. Only those non-bonded pairs that are sufficiently
close to each other will result in an actual non-boned restraint.
If undefined ($-999$), the default value is the maximum of the
three possibilities: twice the radius of the largest atom multiplied
by \K{RADII\_FACTOR} (in the case of the all non-hydrogen atoms model, this
is 3.2 \AA); \K{LENNARD\_JONES\_SWITCH[2]}; or \K{COULOMB\_SWITCH}[2].
Only those values of the three possibilities are compared that have the
corresponding \K{DYNAMIC\_SPHERE}, \K{DYNAMIC\_LENNARD}, or
\K{DYNAMIC\_COULOMB} set to \V{on}.  The best value for
\K{CONTACT\_SHELL} must be found in combination with
\K{UPDATE\_DYNAMIC} (see also below). Good values are 4\AA\ for
\K{CONTACT\_SHELL} and 0.39\AA\ for \K{UPDATE\_DYNAMIC} when no
Lennard-Jones and Coulomb terms are used; if
\K{CONTACT\_SHELL} is larger, there would be many pairs in the non-bonded
pairs list which would slow down the evaluation of the molecular pdf. If it
is too small, however, the increased frequency of the pair list recalculation
may slow down the optimization. It is useful in some simulations
to be able to set \K{CONTACT\_SHELL} to something large (\eg, 8\AA)
and \K{UPDATE\_DYNAMIC} to \V{999999.9}, so that the pairs list is
prepared only at the beginning. However, you have to make sure that
the potential energy is not invisibly pumped into the system by making
contacts that are not on the list of non-bonded pairs (see below).
  
\K{UPDATE\_DYNAMIC} sets the cumulative maximal atomic shift
that triggers recalculation of the list of atom--atom non-bonded
pairs. It should be set in combination with \K{CONTACT\_SHELL}.  For
soft-sphere overlap, to be absolutely sure that no unaccounted
contacts occur, \K{UPDATE\_DYNAMIC} has to be equal to 
(\K{CONTACT\_SHELL} -- {\tt maximal\_overlap\_distance}) / 2. {\tt
Maximal\_overlap\_distance} is equal to the diameter of the largest
atom in the model; it is 3.2 \AA\ in the case of the all non-hydrogen atoms
model. This distance is the \K{CONTACT\_SHELL} value if a default is
requested. Factor 2 comes from the fact that the moves
of both atoms can reduce the distance between them.
\K{DYNAMIC\_SPHERE} has to be set to \V{on} for the automatic 
generation of the soft-sphere overlap restraints. Another necessary
condition is that the scaled standard deviation of the
soft-sphere overlap restraints is greater than zero. It is simpler not 
to pre-calculate any soft-sphere overlap restraints and to use the 
dynamically generated restraints alone, although this may be slower.
  
Similarly, \K{DYNAMIC\_LENNARD}, \K{DYNAMIC\_COULOMB},
\K{DYNAMIC\_MODELLER} and \K{DYNAMIC\_ACCESS} determine whether the 
dynamic Lennard-Jones
terms, electrostatic interactions, \MODELLER\ non-bonded spline
restraints and \MODELLER\ atomic density restraints 
are calculated during optimization. Currently, the first derivatives
of the atom density restraints are set to 0. {\small\tt SHELL here xx.}

\K{EXCL\_LOCAL}[4] specifies whether or not the atoms in a chemical bond, 
chemical angle, dihedral/improper angle, and in the excluded pairs list 
are considered in the construction of the non-bonded atom pairs list.
This is especially useful when simplified protein representations are
used; \eg, when non-bonded restraints need to be used on \Calpha${}_i$ --
\Calpha${}_{i+2}$ terms.
  
The initial atom radii (before scaling by \K{RADII\_FACTOR}) 
depend on \K{TOPOLOGY\_MODEL} which selects a
column of radii for the specified topology model from the {\tt
\$RADII\_LIB} library file.
  
\K{RADII\_FACTOR} is the scaling factor for the atom radii as read
from the library file. The scaled radii are used only for the calculation 
of violations of the soft-sphere overlap restraints. 
  
\K{LENNARD\_JONES\_SWITCH} is a real vector of two elements. It
specifies $r_{min}$ and $r_{max}$ for the Lennard-Jones interaction
(Eq.~\ref{eq:lennard}).
The potential is smoothed down to zero between these two distances.

\K{COULOMB\_SWITCH} is a real vector of two elements. It
specifies $r_{min}$ and $r_{max}$ for the electrostatic interaction
(Eq.~\ref{eq:coulomb}).
The potential is smoothed down to zero between these two distances.

\K{RESIDUE\_SPAN\_RANGE} determines what atom pairs can possibly 
occur in the dynamic non-bonded atom pairs list (see \C{MAKE\_RESTRAINTS}).
\K{RESIDUE\_SPAN\_SIGN} is ignored in \C{OPTIMIZE}. The effect 
of \K{RESIDUE\_SPAN\_RANGE} is
modulated by \K{COVALENT\_CYS}. If \K{COVALENT\_CYS} is 
\V{on}, the disulfide bridges are taken into account when calculating 
the residue index difference between two atoms (\ie, disulfides make some 
atom pairs closer in sequence). \K{COVALENT\_CYS} = \V{on} is slow 
and only has an effect when certain statistical non-bonded potentials
are used (\ie, \K{DYNAMIC\_MODELLER} is \V{on} and the non-bonded library has
been derived considering the disulfide effect). Thus, it should generally
be set to \V{off}. The dynamic restraints include soft-sphere overlap, 
Lennard-Jones, electrostatic restraints, and general spline
restraints. The first three types of restraints can also be
generated as static restraints by \C{MAKE\_RESTRAINTS}. 

% For the dynamic restraints, only the restraints with {\it all} atoms 
% in the selected
% set of atoms are included. This means that the selected set has to be
% large enough to include atoms on the boundary of the selected set that
% will not change much as a result of optimization. This is needed so that 
% new, unaccounted for, violations are not introduced as a result of overlap
% with neighboring atoms not in the selected set.  If such a situation
% happens in the hot spot refinement, hot spot selection and refinement
% should be repeated.
  
The automatically generated dynamic restraints are always deleted after a 
command that calculates them is finished (\C{OPTIMIZE}, \C{ENERGY}, % \C{PMF},
\C{PICK\_HOT\_ATOMS}); you have to use
\C{MAKE\_RESTRAINTS} to calculate equivalent static restraints if you 
want to write the `dynamic' restraints to a file.
  
\K{MIN\_ATOM\_SHIFT} is a convergence criterion for the conjugate
gradients optimization. When the maximal atomic shift is less than the
specified value, the optimization is finished regardless of the number
of optimization cycles or function value and its change.
  
\K{MAX\_ITERATIONS} is used to prevent a waste of CPU time in the 
conjugate gradients optimization. When that many cycles are done, the
optimization is finished regardless of the maximal atomic shift.

Before calculating dynamic non-bonded restraints, \MODELLER\ determines 
which of the several routines is most appropriate and efficient for
calculating the non-bonded atom pairs list. 
The user can influence this selection by specifying two variables:
\K{NONBONDED\_SEL\_ATOMS}, which has an effect when only a subset of all
atoms is selected by the \C{PICK\_ATOMS} or \C{PICK\_HOT\_ATOMS} commands 
(set 1), and \K{NLOGN\_USE}, which has an effect when all atoms are
selected. If \K{NONBONDED\_SEL\_ATOMS} is \V{2} (default), the 
non-bonded pairs will contain only selected atoms (set 1).
This means that the optimized atoms will not ``feel'' the rest of the 
protein through the non-bonded terms at all. If \K{NONBONDED\_SEL\_ATOMS} is
\V{1}, only one of the atoms in the non-bonded pair has to be a selected atom.
This means that the selected region feels the rest of the system through 
the non-bonded terms, at the expense of longer CPU times.
When all atoms are selected, \K{NONBONDED\_SEL\_ATOMS} of course has no
effect. However, in that case, \K{NLOGN\_USE} is used to select either
a straightforward ${\cal O}(n^2)$ search or a cell-based algorithm 
which has $n \log n $ dependency of CPU time \versus\ size $n$. The latter
algorithm is used when the maximal difference in residue indices of the 
atoms in the current dynamic restraints is larger than \K{NLOGN\_USE} or 
when the box size for this algorithm would have to be larger than 8\Ang.
  
The molecular dynamics optimizer pretends that the natural logarithm
of the molecular pdf is energy in kcal/mole. \K{MD\_TIME\_STEP} is the time
step in femtoseconds. \K{TEMPERATURE} is the temperature of the system
in degrees Kelvin. \K{MAX\_ITERATIONS} determines the number of MD steps.
If \K{MD\_RETURN} is \Z{FINAL} the last structure is returned as
the MODEL. If \K{MD\_RETURN} is \Z{MINIMAL} then the structure with
the lowest value of the objective function on the whole trajectory is 
returned as the MODEL.
Rescaling of velocities is done every \K{EQUILIBRATION} steps to match
the specified temperature. Atomic shifts along one axis are limited by
\K{CAP\_ATOM\_SHIFT}. This value should be smaller than
\K{UPDATE\_DYNAMIC}. If \K{INIT\_VELOCITIES} = \V{on}, the velocity
arrays are initialized, otherwise they are not. In that case,
the final velocities from the previous run are used as the initial
velocities for the current run.
  
\K{RAND\_SEED} is the seed for the random number generator. It has
to be between $-2$ and $-50000$. Its value is changed after
the return from the optimization routine.
  
\K{MOLPDF} contains the value of the objective function at the end of
optimization.
  
\K{MODELLER\_STATUS} is set to 1 if optimization is aborted
because dynamic restraints could not be calculated as a result of a
system being too large. If \K{MODELLER\_STATUS} is equal or greater
than \K{STOP\_ON\_ERROR} the execution is stopped. Otherwise the execution
returns back to the \TOP\ routine, exiting all optimization routines
immediately. The execution then continues as if nothing happened.  It
is up to the calling \TOP\ routine to ensure that sensible action is
taken; \eg, skipping the rest of modeling for the model that resulted
in an impossible function evaluation. This option is useful when
calculating several independent models and you do not want one bad
model to abort the whole calculation. A probable reason for an interrupted
optimization is that it was far from convergence by the time the
calculation of dynamic restraints was first requested. Two possible
solutions are: (1) optimize more thoroughly (\ie\ slowly) and (2) use
a different contact pairs routine (\K{SET NLOGN\_USE} = \V{9999}).
\K{MODELLER\_STATUS} can be used in the \TOP\ routine to exit from an 
optimization of a hopeless model and to continue with another model
from a different initial conformation.}


\Example{../examples/commands/optimize.top}





\Command      {SWITCH\_TRACE}{open new optimization trace file}

\Options{
\OptLine{FILE}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DIRECTORY}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{TRACE\_OUTPUT}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
}

\Description  {This command specifies the file for the subsequent optimization 
tracing output. It is useful for separating tracing output for
different models constructed in a single run of \MODELLER. The tracing
output is only produced if \K{TRACE\_OUTPUT} is larger than 0.
The tracing file includes the iteration number, number of
function evaluations, function value, average and maximal atomic
shifts, the size of the gradient vector, kinetic energy (for molecular
dynamics `optimization' only), temperature (MD only) and total energy. 
This is written out in every \K{TRACE\_OUTPUT}-th cycle of whatever
optimization method is used, starting with the state just before the
optimization (iteration 0).

When using the {\tt model} script for comparative modeling, there is 
one {\tt .D} file for each {\tt .B} file with a model. The {\tt .D} files 
contain
information about the progress of optimization, from the beginning to
the end.  The most important column is column 3, which contains the value
of the objective function, which is being optimized, as a function of
the iteration step (every 10 steps, by default). Thus, the best model,
according to \MODELLER, is the one that has the lowest number in the
third column of the last line of its {\tt .D} file. This value is also
written out in the \V{REMARK} record of the PDB file containing the model
and in the \logf\ file.

\vspace{5mm}

\begin{table}[bt]
\begin{center}
\begin{tabular}{ll} \hline
Column & Description \\ \hline
1 & iteration number within one step of the variable target
             function method \\
2 & number of function evaluations within one step of VTFM \\
3 & objective function value \\
4 & average atomic shift \\
5 & maximal atomic shift \\
6 & proportional to the gradient \\
7 & kinetic energy \\
8 & temperature for molecular dynamics optimization \\
9 & total energy (kinetic and potential; potential = objective function
    value) \\ \hline
\end{tabular}
\end{center}
\newcommand{\opttrace}{{\em Columns in an optimization trace file.}}
\caption[\opttrace]
{\latex{\opttrace}}
\label{tab:opttrace}
\underfigtab
\end{table}
}

\SeeExample{\C{OPTIMIZE}}



\Command      {DEBUG\_FUNCTION}{test code self-consistency}

\Options{
\OptLine{DEBUG\_FUNCTION\_CUTOFF}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
\OptLine{DETAILED\_DEBUGGING}{TYPE}{VALUES}{DEFAULT}{DESCRIPTION} 
all the \C{ENERGY} options
}

\Description{This command checks the self-consistency of the code
for the objective function and its derivatives by calculating and
comparing numeric and analytical derivatives. All the parameters
influencing the evaluation of the molecular pdf are also relevant (see
\C{ENERGY}). The derivative is reported if both the absolute difference and
the fractional difference between the two kinds of evaluations are
larger than \K{DEBUG\_FUNCTION\_CUTOFF}[1] and
\K{DEBUG\_FUNCTION\_CUTOFF}[2], respectively.

When \K{DETAILED\_DEBUGGING} is \V{on}, the analytic and numeric 
derivatives of each restraint with respect to atomic positions are
also compared for the atoms `violated' by the whole molecular pdf. 
The absolute cutoff for writing out the discrepancies is
scaled by \K{DEBUG\_FUNCTION\_CUTOFF}[3]; the relative cutoff
remains the same as before.

When \MODELLER\ is compiled in double precision, this test reports a
smaller number of discrepancies.}


\Example{../examples/commands/debug_function.top}

\vfill
