The optimizers module also provides a StateOptimizer class. This class cannot be directly used to optimize the system, but instead it can be used as a base for you to write your own optimization algorithms in Python. To do this, create a subclass and override the optimize method to do your optimization. Your optimizer does not act directly on the atom coordinates, but instead gets a ‘state’ vector with the same number of elements as there are degrees of freedom in the system. (This allows you to also optimize rigid bodies, for example, without having to worry about the specifics of their representation.)
Several utility functions are provided:
If you want to define parameters for your optimization in the same way as the other optimizers, set '_ok_keys' appropriately and then call self.get_parameter() to get their values.
from modeller.optimizers import StateOptimizer class SteepestDescent(StateOptimizer): """Very simple steepest descent optimizer, in Python""" # Add options for our optimizer _ok_keys = StateOptimizer._ok_keys + ('min_atom_shift', 'min_e_diff', 'step_size', 'max_iterations') def __init__(self, step_size=0.0001, min_atom_shift=0.01, min_e_diff=1.0, max_iterations=None, **vars): StateOptimizer.__init__(self, step_size=step_size, min_atom_shift=min_atom_shift, min_e_diff=min_e_diff, max_iterations=max_iterations, **vars) def optimize(self, atmsel, **vars): # Do normal optimization startup StateOptimizer.optimize(self, atmsel, **vars) # Get all parameters alpha = self.get_parameter('step_size') minshift = self.get_parameter('min_atom_shift') min_ediff = self.get_parameter('min_e_diff') maxit = self.get_parameter('max_iterations') # Main optimization loop state = self.get_state() (olde, dstate) = self.energy(state) while True: for i in range(len(state)): state[i] -= alpha * dstate[i] (newe, dstate) = self.energy(state) if abs(newe - olde) < min_ediff: print("Finished at step %d due to energy criterion" % self.step) break elif self.shiftmax < minshift: print("Finished at step %d due to shift criterion" % self.step) break elif maxit is not None and self.step >= maxit: print("Finished at step %d due to step criterion" % self.step) break if newe < olde: alpha *= 2 else: alpha /= 2 olde = newe self.next_step() self.finish()
from modeller import * from modeller.optimizers import actions from modeller.scripts import complete_pdb # Load our custom steepest descent optimizer from steepest_descent import SteepestDescent env = Environ() env.io.atom_files_directory = ['../atom_files'] env.libs.topology.read(file='$(LIB)/top_heav.lib') env.libs.parameters.read(file='$(LIB)/par.lib') # Read in the initial structure: code = '1fdn' mdl = complete_pdb(env, code) atmsel = Selection(mdl) # Generate the restraints: mdl.restraints.make(atmsel, restraint_type='stereo', spline_on_site=False) # Optimize with our custom optimizer: opt = SteepestDescent(max_iterations=80) opt.optimize(atmsel, actions=actions.Trace(5))