| edat = <energy_data> | objective function parameters | |
| schedule_scale = <physical.values> | 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 | factors for physical restraint types in scaling the schedule |
| residue_span_range = <int:2> | 0 99999 | range of residues spanning the allowed distances; for MAKE_RESTRAINTS, PICK_RESTRAINTS, non-bonded dynamic pairs |
| max_iterations = <int:1> | 200 | maximal iterations in optimization |
| output = <str:1> | 'LONG' | 'NO_REPORT' | 'REPORT' |
| min_atom_shift = <float:1> | 0.010 | minimal atomic shift for the optimization convergence test |
| actions = [] | list of periodic actions |
This command creates a new Python optimizer object. Calling the object's optimize method with an atom selection then performs a number of optimizing iterations using a modified version of the Beale restart conjugate gradients method [Shanno & Phua, 1980,Shanno & Phua, 1982]. A brief description of the algorithm is given in Section A.2.
The optimization can be controlled with a number of keyword arguments, which can be specified either when the object is created, or when the optimize method is called (if the same keyword is specified in both, that for the optimize method takes precedence). If unspecified, the variables default to whatever is already set for the atom selection.
min_atom_shift is a convergence criterion for the 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 optimization. When that many cycles are done, the optimization is finished regardless of the maximal atomic shift.
A summary of the optimization results is written to the log file after optimization, unless output contains string 'NO_REPORT'.
edat is an energy_data object containing objective function parameters. See section 6.3 for more information.
The scaling factors for the physical restraint types are given by schedule_scale.
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 non-bonded atom pairs list used for dynamic restraints (see Section 5.3).
The optimize method, when called, returns molpdf, 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).
If actions is set, it should be a list of periodic actions. Each is a Python object containing an action which is carried out periodically during the optimization, after every step. For example, actions.write_structure() can be used to write out a PDB file with structure snapshots during the run, while actions.trace() writes basic information about the optimization to a trace file. If multiple actions are given, they are run in the order they are given.
# Example for: conjugate_gradients(), molecular_dynamics(), model.switch_trace()
# This will optimize stereochemistry of a given model, including
# non-bonded contacts.
from modeller import *
from modeller.scripts import complete_pdb
from modeller.optimizers import conjugate_gradients, molecular_dynamics, actions
env = environ()
env.io.atom_files_directory = '../atom_files'
env.edat.dynamic_sphere = True
env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')
code = '1fas'
mdl = complete_pdb(env, code)
mdl.write(file=code+'.ini')
# Select all atoms:
atmsel = selection(mdl)
# Generate the restraints:
mdl.restraints.make(atmsel, restraint_type='stereo', spline_on_site=False)
mdl.restraints.write(file=code+'.rsr')
mpdf = atmsel.energy()
# Create optimizer objects and set defaults for all further optimizations
cg = conjugate_gradients(output='REPORT')
md = molecular_dynamics(output='REPORT')
# Open a file to get basic stats on each optimization
trcfil = file(code+'.D00000001', 'w')
# Run CG on the all-atom selection; write stats every 5 steps
cg.optimize(atmsel, max_iterations=20, actions=actions.trace(5, trcfil))
# Run MD; write out a PDB structure (called '1fas.D9999xxxx.pdb') every
# 10 steps during the run, and write stats every 10 steps
md.optimize(atmsel, temperature=300, max_iterations=50,
actions=[actions.write_structure(10, code+'.D9999%04d.pdb'),
actions.trace(10, trcfil)])
# Finish off with some more CG, and write stats every 5 steps
cg.optimize(atmsel, max_iterations=20,
actions=[actions.trace(5, trcfil)])
mpdf = atmsel.energy()
mdl.write(file=code+'.B')