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). Valid keywords are:
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).
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).
# 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')