User-defined optimizers

The optimizers module also provides a state_optimizer 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:

• 'self.get_state()': returns a state vector representing the current state of the system (x,y,z coordinates of all non-rigid atoms in the selection, and center of mass and rotation angles of all selected rigid bodies).

• 'self.energy(state)': given a state vector, returns the system energy and a similar vector of state gradients. Also updates the atom shifts self.shiftavr and self.shiftmax from the previous state.

• 'self.next_step()': updates the step counter self.step, and does any periodic actions, if defined.

• 'self.finish()': does any cleanup at the end of the optimization.

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.

Example: examples/python/steepest_descent.py

from modeller.optimizers import state_optimizer

class steepest_descent(state_optimizer):
    """Very simple steepest descent optimizer, in Python"""

    # Add options for our optimizer
    _ok_keys = state_optimizer._ok_keys + ('min_atom_shift', 'min_e_diff',
                                           'step_size', 'max_iterations')

    def optimize(self, atmsel, **vars):
        # Do normal optimization startup
        state_optimizer.optimize(self, atmsel, **vars)

        # Get all parameters
        alpha = self.get_parameter('step_size', 0.0001)
        minshift = self.get_parameter('min_atom_shift', 0.01)
        min_ediff = self.get_parameter('min_e_diff', 1.0)
        maxit = self.get_parameter('max_iterations', None)

        # Main optimization loop
        state = self.get_state()
        (olde, dstate) = self.energy(state)
        while True:
            for i in range(len(state)):
                state[i] = 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()

Example: examples/python/steepest_descent_test.py

from modeller import *
from modeller.optimizers import actions
from modeller.scripts import complete_pdb

# Load our custom steepest descent optimizer
from steepest_descent import steepest_descent

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 = '1fdx'
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 = steepest_descent(max_iterations=80)
opt.optimize(atmsel, actions=actions.trace(5))