User-defined feature types

To create a new feature type, derive a new class from the base features.feature. You should then set the numatoms member to the number of atoms your feature acts on, and also override the following functions: eval, deriv, and is_angle. You can also derive your class from any of the built-in MODELLER features (e.g., features.angle) if you desire.

The eval function is called from MODELLER with a model object and the indices of the atoms defining the feature. Your function should return the value of the feature. The deriv function is similar, but is also passed the current feature value; you should return the derivatives of the feature with respect to x, y and z of each defining atom. The is_angle function should return True if your feature is an angle, in which case MODELLER will automatically deal with periodicity for you, and convert any feature values to degrees for the user. (Your eval and deriv functions should, however, return angle values in radians.)

Example: examples/python/user_feat.py

from modeller import *
from modeller.scripts import complete_pdb

env = environ()

env.io.atom_files_directory = ['../atom_files']
log.verbose()
env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')

class MyDist(features.feature):
    """An implementation of Modeller's distance feature (type 1) in
       pure Python. For improved performance, see cuser_feat.py, which
       implements the feature in C."""

    numatoms = 2

    def eval(self, mdl, atom_indices):
        (a1, a2) = self.indices_to_atoms(mdl, atom_indices)
        dist = ((a1.x - a2.x) ** 2 + (a1.y - a2.y) ** 2
                + (a1.z - a2.z) ** 2) ** 0.5
        return dist

    def deriv(self, mdl, atom_indices, feat):
        (a1, a2) = self.indices_to_atoms(mdl, atom_indices)
        dx1 = (a1.x - a2.x) / feat
        dy1 = (a1.y - a2.y) / feat
        dz1 = (a1.z - a2.z) / feat
        dx2 = -dx1
        dy2 = -dy1
        dz2 = -dz1
        return ((dx1, dx2), (dy1, dy2), (dz1, dz2))

    def is_angle(self):
        return False


mdl = complete_pdb(env, "1fdn")
sel = selection(mdl)
rsr = mdl.restraints
at = mdl.atoms
rsr.add(forms.gaussian(group=physical.bond,
                       feature=MyDist(at['CA:1:A'], at['C:1:A']),
                       mean=1.5380, stdev=0.0364))
sel.energy()