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()