model.energy() -- evaluate MODEL given restraints

edat = <energy_data>		objective function parameters
viol_report_cut = <float:35>	4.5 4.5 4.5 4.5 4.5 4.5 4.5 4.5 4.5 4.5 4.5 4.5 4.5 999 999 999 999 4.5 4.5 4.5 4.5 4.5 4.5 999 6.5 4.5 4.5 4.5 4.5 4.5 999 999 999 4.5 4.5	cutoffs for reporting relative violations
viol_report_cut2 = <float:35>	2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
output = <str:1>	'LONG'	'SHORT' | 'LONG' | 'VERY_LONG' | 'GRADIENT' | 'SYMMETRY' | 'ENERGY_PROFILE' | 'VIOLATIONS_PROFILE'
normalize_profile = <bool:1>	False	whether to normalize energy/violations profiles or not, by the number of terms per residue
smoothing_window = <int:1>	3	profiles are smoothed over 2*SW + 1 residues
schedule_scale = <float:35>	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	factors for physical restraint types in scaling the schedule
file = <str:1>	'default'	partial or complete filename
asgl_output = <bool:1>	False	whether to write output for ASGL
residue_span_range = <int:2>	0 99999	range of residues spanning the allowed distances; for MAKE_RESTRAINTS, PICK_RESTRAINTS, non-bonded dynamic pairs

Output:

molpdf

Requirements:

restraints

Description:

The main purpose of this command is to compare spatial features of the current MODEL with the selected restraints in order to determine the violations of the molecular pdf. It lists variable amounts of information about the values of the basis, feature, and molecular pdf's for the current MODEL. All arguments that affect the value of the molecular pdf are also relevant for the model.energy() command.

Within this routine only, the scaling factors for the physical restraint types are obtained from the model.schedule.step step of the current schedule, multipled by schedule_scale (the original values are returned upon exit from the routine). This allows easy reporting of only a selected subset of all restraints.

Most of the output goes to the log file. The output of the model.energy() command has to be examined carefully, at least at the end of the optimization, when the final model is produced. Additional output files, for the ASGL plotting program are created if asgl_output = True (undocumented).

output selects various kinds of output information:

• 'LONG' writes restraint violations one per line to the log file.

• 'VERY_LONG' writes the most detailed examination of the selected basis and feature pdf's to the log file, using several lines of output for each restraint.

• 'GRADIENT' writes the `force' gradients for the currently selected restraints to the isotropic temperature factors for each atom of the current MODEL.

• 'SYMMETRY' writes a comparison of equivalent distances involved in the definition of the symmetry enforcing term to the log file.

viol_report_cut is a vector with one real number for each physical restraint type. A restraint is reported when its `heavy relative violation' is larger than the corresponding cutoff. The heavy relative violation is calculated by finding the global minimum of a feature according to the restraint, taking the difference between the actual feature in the model and this global minimum, and then normalizing the difference by the standard deviation of the global minimum. The `minimal violation' of a restraint is defined as the difference from the local minimum closest to the value of the feature in the model (with the exception of the spline restraints; see next paragraph).

viol_report_cut2 is similar to viol_report_cut, except that it contains cutoffs for restraint `energies', not heavy relative violations.

The meaning of various other reported properties of the violated restraints is briefly described in the log file. Note that for multi-modal restraints that are described by cubic splines (by default, all multimodal homology-derived restraints), only one optimal value is defined, not the local and global minimum as for the multi-modal Gaussian restraints. As a result, the minimal violations and heaviest violations are the same. For interpreting the seriousness of violations, use the following rule of thumb: There should be at most a few small violations (e.g., 4 standard deviations) for all monomodal restraints. In comparative modeling, the monomodal restraints include the stereochemical restraints and distance restraints when only one homologous structure is used. For the multimodal restraints, there are usually many violations reported because the heaviest violations are used in deciding whether or not to report a violation. In comparative modeling, the multimodal restraints include the $\chi_i$ restraints, ($\Phi$, $\Psi$) binormal restraints and distance restraints when more than one template is used. See also Section 1.8, Question 15.

For profiles:

This command calculates residue energies or heavy relative violations, depending on output, for all physical restraint types (there are NPHYCNS of them). Relative heavy violations (Table 3.3) are used because only relative violations of different features are comparable. In both cases, the residue sum is the sum over all restraints that have at least one atom in a given residue. The contribution of each restraint is counted exactly once for each residue, without any weighting. Restraints spanning more than one residue contribute equally to all of them. Thus, the sum of residue energies is generally larger than molecular pdf. The command also calculates the sum of the NPHYCNS contributions for each residue and writes all NPHYCNS+1 columns to a file suitable for plotting by ASGL.

If normalize_profile is True the profile for each residue is normalized by the number of terms applying to each residue.

All the curves are smoothed by the running window averaging method if smoothing_window is larger than 0: The window is centered on residue $i$ and extends for (smoothing_window/2) - 1 residues on each side. Thus, smoothing_window has to be an even number (or it is made such by the program automatically). The only exceptions are the two terminii, where a smaller number of residues are available for smoothing. The relative weight of residue $j$ when calculating the smoothed value at residue $i$ is $({\sf smoothing\_window}\index{smoothing\_window@{\sf smoothing\_window}}/2 - \vert j-i\vert)$.

The energy or the violations profile is written to the fourth column of the MODEL atomic records (atomic isotropic temperature factors for X-ray structures). Note that all the atoms in one residue get the same number. This output is useful for exploring the violations on a graphics terminal.

See description of model.optimize() for the other variables.

Example: examples/commands/energy.py

# Example for: model.energy()

# This will calculate the stereochemical energy (bonds,
# angles, dihedrals, impropers) for a given model.

env = environ()
env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')

aln = alignment(env)
code = "1fas"
mdl = model(env, file=code)
aln.append_model(mdl, atom_files=code, align_codes=code)
aln.append_model(mdl, atom_files=code+'.ini', align_codes=code+'-ini')

mdl.generate_topology(aln, sequence=code+'-ini')
# Must patch disulfides here to calculate the non-bonded
# energy properly. Also, when you use hydrogens, disulfides
# must always be patched so that sulfhydril hydrogens are
# removed from the model.
for ids in [ ('17', '39'),
             ( '3', '22'),
             ('53', '59'),
             ('41', '52') ]:
    mdl.patch(residue_type='DISU', residue_ids=ids)

mdl.transfer_xyz(aln)
mdl.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')

mdl.restraints.make(aln, restraint_type='stereo', spline_on_site=False)
molpdf = mdl.energy(edat=energy_data(dynamic_sphere=True))