model.optimize() -- optimize MODEL given restraints

edat = <energy_data>		objective function parameters
optimization_method = <int:1>	999	type of optimization method: 1 | 3
residue_span_range = <int:2>	0 99999	range of residues spanning the allowed distances; for MAKE_RESTRAINTS, PICK_RESTRAINTS, non-bonded dynamic pairs
trace_output = <int:1>	0	modulus for writing information about optimization iterations: 0 for nothing
max_iterations = <int:1>	200	maximal iterations in optimization
output = <str:1>	'LONG'	'NO_REPORT' | 'REPORT'
$\bullet$ For conjugate gradients:
min_atom_shift = <float:1>	0.010	minimal atomic shift for the optimization convergence test
$\bullet$ For molecular dynamics:
md_time_step = <float:1>	4.0	time step for MD in fs
init_velocities = <bool:1>	True	whether to initialize velocities before MD
temperature = <float:1>	293.0	temperature for MD simulation in K
equilibrate = <int:1>	999999	equilibrate during MD every that many steps
md_return = <str:1>	'FINAL'	return MODEL with 'MINIMAL' energy or 'FINAL' MODEL
cap_atom_shift = <float:1>	0.2	limit for atomic shifts in optimization

Output:

molpdf

Requirements:

restraints

Description:

This command performs a number of optimizing iterations using a selected optimization method (5.2). One call to model.optimize() corresponds to a single step of the variable target function method. The whole variable target function method is implemented by the automodel class. 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 model.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.

model.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 3.7.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.

edat is an energy_data object containing objective function parameters. See section 3.3 for more information.

It is useful in some simulations to be able to set energy_data.contact_shell to something large (e.g., 8Å) and energy_data.update_dynamic to 999999.9, so that the pairs list is prepared only at the beginning of the optimization. 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).

residue_span_range determines what atom pairs can possibly occur in the dynamic non-bonded atom pairs list (see model.restraints.make()). residue_span_sign is ignored in model.optimize(). 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 model.restraints.make().

The automatically generated dynamic restraints are always deleted after a command that calculates them is finished (model.optimize(), model.energy(), model.pick_hot_atoms()); you have to use model.restraints.make() 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.

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 energy_data.update_dynamic. If init_velocities = True, 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.

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

An exception is raised if optimization is aborted because dynamic restraints could not be calculated as a result of a system being too large. It is up to the calling script 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 energy_data.nlogn_use = 9999).

Example: examples/commands/optimize.py

# Example for: model.optimize(), misc.switch_trace()

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

env = environ()
env.edat.dynamic_sphere = True

env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')

code = '1fas'
aln = alignment(env)
mdl = model(env, file=code)
aln.append_model(mdl, align_codes=code, atom_files=code)
aln.append_model(mdl, align_codes=code+'-ini', atom_files=code+'.ini')
mdl.generate_topology(aln, sequence=code+'-ini')
mdl.transfer_xyz(aln)
mdl.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')
mdl.write(file=code+'.ini')

# Generate the restraints:
mdl.restraints.make(aln, restraint_type='stereo', spline_on_site=False)
mdl.restraints.write(file=code+'.rsr')

mpdf = mdl.energy()
mdl.switch_trace(file=code+'.trc')
mdl.optimize(optimization_method=1, max_iterations=20, output='REPORT',
             trace_output=1)
mdl.optimize(optimization_method=3, temperature=300, max_iterations=50,
             output='REPORT', trace_output=1)
mdl.optimize(optimization_method=1, max_iterations=20, output='REPORT',
             trace_output=1)
mpdf = mdl.energy()

mdl.write(file=code+'.B')