User-defined restraint forms

To create a new restraint form, derive a new class from the base
forms.restraint_form. You should then override the following
functions: __init__, eval, vmin, rvmin, min_mean,
vheavy, rvheavy, heavy_mean, and get_range. Note that
presently you can
only derive from this base class, not from MODELLER built-in forms.

Restraint forms can act on one or more features (each of which has an
accompanying integer modality, which you can use for any purpose), and can
take any number of floating-point parameters
as input. The features and parameters are stored in self._features
and self._parameters respectively, but for convenience the base
constructor restraint_form.__init__ can set initial values for these.

The eval function is called from MODELLER with the current feature values,
their types and modalities, and the parameter vector. You should return the
objective function contribution and, if requested, the derivatives with
respect to each feature. The feature types are required by the deltaf
function, which returns the difference between the current feature value and
the mean (a simple subtraction is not sufficient, as some feature types are
periodic). Note that you must use the passed parameter vector, as the class
is not persistent, and as such the self._parameters variable (or any other
object variable you may have set) is not available to this function.

The get_range function is used to define the feature range over which
the form is clearly non-linear. It is simply passed a similar set of
parameters to eval, and should return a 2-element tuple containing the
minimum and maximum feature values. It is only necessary to define this
function if the form acts on only a single feature and you want to be able
to convert it to a cubic spline using Restraints.spline().

The other functions are used to return the minimal
and heavy restraint violations (both absolute and relative; see
Section 5.3.1) and the means. The heavy and minimal
means correspond to the global and local minima.

frommodellerimport*frommodeller.scriptsimportcomplete_pdbenv=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')classMyGauss(forms.restraint_form):"""An implementation of Modeller's harmonic/Gaussian restraint (type 3) in pure Python"""rt=0.5900991# RT at 297.15K, in kcal/moldef__init__(self,group,feature,mean,stdev):forms.restraint_form.__init__(self,group,feature,0,(mean,stdev))defeval(self,feats,iftyp,modal,param,deriv):(mean,stdev)=paramdelt=self.deltaf(feats[0],mean,iftyp[0])val=self.rt*0.5*delt**2/stdev**2ifderiv:fderv=self.rt*delt/stdev**2returnval,[fderv]else:returnvaldefvmin(self,feats,iftyp,modal,param):(mean,stdev)=paramreturnself.deltaf(feats[0],mean,iftyp[0])defrvmin(self,feats,iftyp,modal,param):(mean,stdev)=paramreturnself.deltaf(feats[0],mean,iftyp[0])/stdevdefmin_mean(self,feats,iftyp,modal,param):(mean,stdev)=paramreturn[mean]defget_range(self,iftyp,modal,param,spline_range):(mean,stdev)=paramreturn(mean-stdev*spline_range,mean+stdev*spline_range)# There is only one minimum, so the 'heavy' mean is the same as the 'min'vheavy=vminrvheavy=rvminheavy_mean=min_meanmdl=complete_pdb(env,"1fdn")sel=selection(mdl)rsr=mdl.restraintsat=mdl.atomsrsr.add(MyGauss(group=physical.bond,feature=features.distance(at['CB:1:A'],at['CA:1:A']),mean=1.5380,stdev=0.0364))sel.energy()# Restraints using user-defined forms can be converted to splines for speed.# Only one-dimensional forms that define the get_range() method can be splined.rsr.spline(MyGauss,features.distance,physical.bond,spline_dx=0.05)sel.energy()