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.)
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'], at['C:1']), mean=1.5380, stdev=0.0364)) sel.energy()