OPTIMIZE -- optimize MODEL given restraints

OPTIMIZATION_METHOD = <integer:1>	999	type of optimization method: 1 | 3
SCHEDULE_STEP = <integer:1>	1	schedule step for optimization
TOPOLOGY_MODEL = <integer:1>	3	selects topology library: 1-10
RADII_FACTOR = <real:1>	0.82	factor for van der Waals radii
SPHERE_STDV = <real:1>	0.05	standard deviation of soft-sphere repulsion
DYNAMIC_SPHERE = <logical:1>	on	whether to use dynamic soft-sphere repulsion terms
DYNAMIC_LENNARD = <logical:1>	off	whether to use dynamic Lennard-Jones energy terms
DYNAMIC_COULOMB = <logical:1>	off	whether to use dynamic Coulomb energy terms
DYNAMIC_MODELLER = <logical:1>	off	whether to use dynamic MODELLER non-bonded restraints
DYNAMIC_ACCESS = <logical:1>	off	whether to use dynamic accessibility energy terms
EXCL_LOCAL = <logical:4>	on on on on	whether to exclude bonds, angles, dihedrals, explicit excl pairs from the homology-derived distance rsrs
LENNARD_JONES_SWITCH = <real:2>	6.5 7.5	the range for Lennard-Jones interaction smoothing to 0
COULOMB_SWITCH = <real:2>	6.5 7.5	the range for Coulomb interaction smoothing to 0
RELATIVE_DIELECTRIC = <real:1>	1.0	relative dielectric constant
NONBONDED_SEL_ATOMS = <integer:1>	1	a non-bonded pair has to have at least as many selected atoms
RESIDUE_SPAN_RANGE = <integer:2>	0 99999	range of residues spanning the allowed distances; for MAKE_RESTRAINTS, PICK_RESTRAINTS, non-bonded dynamic pairs
COVALENT_CYS = <logical:1>	off	whether to consider SG-SG covalent bond similar to polypeptide chain when proximity of residues along the sequence is considered. If PATCH_SS_MODEL is done, then make it ON.
CONTACT_SHELL = <real:1>	4.0	distance cutoff for calculation of the non-bonded pairs list
UPDATE_DYNAMIC = <real:1>	0.39	when to update non-bonded pairs list
NLOGN_USE = <integer:1>	15	number of residues at which to begin using the N Log N non-bonded pairs routine
TRACE_OUTPUT = <integer:1>	0	modulus for writing information about optimization iterations: 0 for nothing
MAX_ITERATIONS = <integer:1>	200	maximal iterations in optimization
OUTPUT = <string:1>	'LONG'	'NO_REPORT' | 'REPORT'
$\bullet$ For conjugate gradients:
MIN_ATOM_SHIFT = <real:1>	0.010	minimal atomic shift for the optimization convergence test
$\bullet$ For molecular dynamics:
MD_TIME_STEP = <real:1>	4.0	time step for MD in fs
INIT_VELOCITIES = <logical:1>	on	whether to initialize velocities before MD
TEMPERATURE = <real:1>	293.0	temperature for MD simulation in K
EQUILIBRATE = <integer:1>	999999	equilibrate during MD every that many steps
MD_RETURN = <string:1>	'FINAL'	return MODEL with 'MINIMAL' energy or 'FINAL' MODEL
CAP_ATOM_SHIFT = <real:1>	0.2	limit for atomic shifts in optimization
RAND_SEED = <integer:1>	8123	random seed from -50000 to -2
STOP_ON_ERROR = <integer:1>	1	whether to stop on error

Output:

MOLPDF, MODELLER_STATUS

Requirements:

MODEL & restraints

Description:

This command performs a number of optimizing iterations using a selected optimization method (5.2). One call to 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 TRACE_OUTPUT iterations of the optimizer if TRACE_OUTPUT is larger than 0 (see the SWITCH_TRACE command).

In addition, a summary of the optimization results is written to the log file after optimization, unless OUTPUT contains string 'NO_REPORT'.

OPTIMIZATION_METHOD = 1 selects a conjugate gradients optimization method. 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 [Shanno & Phua, 1980,Shanno & Phua, 1982]. 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 [Verlet, 1967]. All atomic masses are set to that of carbon 12. A brief description of the algorithms is given in Section 5.2.

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 2.5.3); (3) the individual scaling factors for all the physical restraint types. OPTIMIZATION_METHOD overrides the schedule specification if it is within a defined range.

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 RADII_FACTOR (in the case of the all non-hydrogen atoms model, this is 3.2 Å); LENNARD_JONES_SWITCH[2]; or COULOMB_SWITCH[2]. Only those values of the three possibilities are compared that have the corresponding DYNAMIC_SPHERE, DYNAMIC_LENNARD, or DYNAMIC_COULOMB set to on. The best value for CONTACT_SHELL must be found in combination with UPDATE_DYNAMIC (see also below). Good values are 4Å for CONTACT_SHELL and 0.39Å for UPDATE_DYNAMIC when no Lennard-Jones and Coulomb terms are used; if 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 CONTACT_SHELL to something large (e.g., 8Å) and UPDATE_DYNAMIC to 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).

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 CONTACT_SHELL. For soft-sphere overlap, to be absolutely sure that no unaccounted contacts occur, UPDATE_DYNAMIC has to be equal to (CONTACT_SHELL - maximal_overlap_distance) / 2. Maximal_overlap_distance is equal to the diameter of the largest atom in the model; it is 3.2 Å in the case of the all non-hydrogen atoms model. This distance is the 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. DYNAMIC_SPHERE has to be set to 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, DYNAMIC_LENNARD, DYNAMIC_COULOMB, DYNAMIC_MODELLER and 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. SHELL here xx.

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; e.g., when non-bonded restraints need to be used on ${C}_\alpha$${}_i$ - ${C}_\alpha$${}_{i+2}$ terms.

The initial atom radii (before scaling by RADII_FACTOR) depend on TOPOLOGY_MODEL which selects a column of radii for the specified topology model from the $RADII_LIB library file.

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.

LENNARD_JONES_SWITCH is a real vector of two elements. It specifies $r_{min}$ and $r_{max}$ for the Lennard-Jones interaction (Eq. 5.63). The potential is smoothed down to zero between these two distances.

COULOMB_SWITCH is a real vector of two elements. It specifies $r_{min}$ and $r_{max}$ for the electrostatic interaction (Eq. 5.60). The potential is smoothed down to zero between these two distances.

RESIDUE_SPAN_RANGE determines what atom pairs can possibly occur in the dynamic non-bonded atom pairs list (see MAKE_RESTRAINTS). RESIDUE_SPAN_SIGN is ignored in OPTIMIZE. The effect of RESIDUE_SPAN_RANGE is modulated by COVALENT_CYS. If COVALENT_CYS is on, the disulfide bridges are taken into account when calculating the residue index difference between two atoms (i.e., disulfides make some atom pairs closer in sequence). COVALENT_CYS = on is slow and only has an effect when certain statistical non-bonded potentials are used (i.e., DYNAMIC_MODELLER is on and the non-bonded library has been derived considering the disulfide effect). Thus, it should generally be set to 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 MAKE_RESTRAINTS.

The automatically generated dynamic restraints are always deleted after a command that calculates them is finished (OPTIMIZE, ENERGY, PICK_HOT_ATOMS); you have to use MAKE_RESTRAINTS to calculate equivalent static restraints if you want to write the `dynamic' restraints to a file.

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.

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: NONBONDED_SEL_ATOMS, which has an effect when only a subset of all atoms is selected by the PICK_ATOMS or PICK_HOT_ATOMS commands (set 1), and NLOGN_USE, which has an effect when all atoms are selected. If NONBONDED_SEL_ATOMS is 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 NONBONDED_SEL_ATOMS is 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, NONBONDED_SEL_ATOMS of course has no effect. However, in that case, 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 NLOGN_USE or when the box size for this algorithm would have to be larger than 8.

The molecular dynamics optimizer pretends that the natural logarithm of the molecular pdf is energy in kcal/mole. MD_TIME_STEP is the time step in femtoseconds. TEMPERATURE is the temperature of the system in degrees Kelvin. MAX_ITERATIONS determines the number of MD steps. If MD_RETURN is 'FINAL' the last structure is returned as the MODEL. If MD_RETURN is '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 EQUILIBRATION steps to match the specified temperature. Atomic shifts along one axis are limited by CAP_ATOM_SHIFT. This value should be smaller than UPDATE_DYNAMIC. If INIT_VELOCITIES = 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.

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.

MOLPDF contains the value of the objective function at the end of optimization.

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 MODELLER_STATUS is equal or greater than 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; e.g., 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 (i.e. slowly) and (2) use a different contact pairs routine (SET NLOGN_USE = 9999). 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:

# Example for: OPTIMIZE, SWITCH_TRACE

# This will optimize stereochemistry of a given model, including
# non-bonded contacts.

READ_TOPOLOGY FILE = '$(LIB)/top_heav.lib'
READ_PARAMETERS FILE = '$(LIB)/par.lib'
READ_MODEL FILE = '1fas'
SEQUENCE_TO_ALI ATOM_FILES = '1fas', ALIGN_CODES = '1fas'
SEQUENCE_TO_ALI ADD_SEQUENCE = on, ATOM_FILES = ATOM_FILES '1fas.ini', ;
                ALIGN_CODES = ALIGN_CODES '1fas-ini'
GENERATE_TOPOLOGY SEQUENCE = '1fas-ini'
TRANSFER_XYZ
BUILD_MODEL INITIALIZE_XYZ = off
WRITE_MODEL FILE = '1fas.ini'

# Generate the restraints:
MAKE_RESTRAINTS RESTRAINT_TYPE = 'stereo'
WRITE_RESTRAINTS FILE = '1fas.rsr'
ENERGY DYNAMIC_SPHERE = on
SWITCH_TRACE TRACE_OUTPUT = 1, FILE = '1fas.trc'
OPTIMIZE OPTIMIZATION_METHOD = 1, MAX_ITERATIONS = 20
OPTIMIZE OPTIMIZATION_METHOD = 3, TEMPERATURE = 300, MAX_ITERATIONS = 50
OPTIMIZE OPTIMIZATION_METHOD = 1, MAX_ITERATIONS = 20
ENERGY

WRITE_MODEL FILE = '1fas.B'