6 import modeller.scripts
7 import modeller.optimizers
11 if not hasattr(modeller.terms,
'EnergyTerm'):
12 modeller.terms.EnergyTerm = modeller.terms.energy_term
13 modeller.Selection = modeller.selection
16 class _TempDir(object):
17 """Make a temporary directory that is deleted when the object is."""
20 self.tmpdir = tempfile.mkdtemp()
23 shutil.rmtree(self.tmpdir, ignore_errors=
True)
27 """A Modeller restraint which evaluates an IMP scoring function.
28 This can be used to incorporate IMP Restraints into an existing
29 comparative modeling pipeline, or to use Modeller optimizers or
33 _physical_type = modeller.physical.absposition
35 def __init__(self, particles, scoring_function=None):
37 @param particles A list of the IMP atoms (as Particle objects),
38 same order as the Modeller atoms.
39 @param scoring_function An IMP::ScoringFunction object that will
40 be incorporated into the Modeller score (molpdf).
41 @note since Modeller, unlike IMP, is sensitive to the ordering
42 of atoms, it usually makes sense to create the model in
43 Modeller and then use ModelLoader to load it into IMP,
44 since that will preserve the Modeller atom ordering in IMP.
46 modeller.terms.EnergyTerm.__init__(self)
47 self._particles = particles
49 self._sf = scoring_function
51 self._sf = particles[0].get_model()
53 def eval(self, mdl, deriv, indats):
54 atoms = self.indices_to_atoms(mdl, indats)
55 _copy_modeller_coords_to_imp(atoms, self._particles)
56 if len(self._particles) == 0:
59 score = self._sf.evaluate(deriv)
61 dvx = [0.] * len(indats)
62 dvy = [0.] * len(indats)
63 dvz = [0.] * len(indats)
64 _get_imp_derivs(self._particles, dvx, dvy, dvz)
65 return (score, dvx, dvy, dvz)
71 """An IMP restraint using all defined Modeller restraints.
72 This is useful if you want to use Modeller restraints with an IMP
73 optimizer, or in combination with IMP restraints.
75 @note Currently only the coordinates of the atoms are translated
76 between Modeller and IMP; thus, a Modeller restraint which
77 uses any other attribute (e.g. charge) will not react if
78 this attribute is changed by IMP.
81 def __init__(self, model, modeller_model, particles):
83 @param model The IMP Model object.
84 @param modeller_model The Modeller model object.
85 @param particles A list of the IMP atoms (as Particle objects),
86 in the same order as the Modeller atoms.
87 @note since Modeller, unlike IMP, is sensitive to the ordering
88 of atoms, it usually makes sense to create the model in
89 Modeller and then use ModelLoader to load it into IMP,
90 since that will preserve the Modeller atom ordering in IMP.
93 if hasattr(x,
'get_particle'):
94 return x.get_particle()
97 IMP.Restraint.__init__(self, model,
"ModellerRestraints %1%")
98 self._modeller_model = modeller_model
99 self._particles = [get_particle(x)
for x
in particles]
102 atoms = self._modeller_model.atoms
103 sel = modeller.Selection(self._modeller_model)
104 _copy_imp_coords_to_modeller(self._particles, atoms)
105 energies = sel.energy()
107 _add_modeller_derivs_to_imp(atoms, self._particles, accum)
113 def do_show(self, fh):
114 fh.write(
"ModellerRestraints")
116 return self._particles
119 def _copy_imp_coords_to_modeller(particles, atoms):
120 """Copy atom coordinates from IMP to Modeller"""
124 for (num, at)
in enumerate(atoms):
125 at.x = particles[num].get_value(xkey)
126 at.y = particles[num].get_value(ykey)
127 at.z = particles[num].get_value(zkey)
130 def _copy_modeller_coords_to_imp(atoms, particles):
131 """Copy atom coordinates from Modeller to IMP"""
135 for (num, at)
in enumerate(atoms):
136 particles[num].set_value(xkey, at.x)
137 particles[num].set_value(ykey, at.y)
138 particles[num].set_value(zkey, at.z)
141 def _add_modeller_derivs_to_imp(atoms, particles, accum):
142 """Add atom derivatives from Modeller to IMP"""
143 for (num, at)
in enumerate(atoms):
145 xyz.add_to_derivative(0, at.dvx, accum)
146 xyz.add_to_derivative(1, at.dvy, accum)
147 xyz.add_to_derivative(2, at.dvz, accum)
150 def _get_imp_derivs(particles, dvx, dvy, dvz):
151 """Move atom derivatives from IMP to Modeller"""
155 for idx
in range(0, len(dvx)):
156 dvx[idx] = particles[idx].get_derivative(xkey)
157 dvy[idx] = particles[idx].get_derivative(ykey)
158 dvz[idx] = particles[idx].get_derivative(zkey)
162 def _HarmonicLowerBoundGenerator(parameters, modalities):
163 (mean, stdev) = parameters
167 def _HarmonicUpperBoundGenerator(parameters, modalities):
168 (mean, stdev) = parameters
172 def _HarmonicGenerator(parameters, modalities):
173 (mean, stdev) = parameters
177 def _CosineGenerator(parameters, modalities):
178 (phase, force_constant) = parameters
179 (periodicity,) = modalities
182 def _LinearGenerator(parameters, modalities):
183 (scale,) = parameters
186 def _SplineGenerator(parameters, modalities):
187 (open, low, high, delta, lowderiv, highderiv) = parameters[:6]
189 for v
in parameters[6:]:
197 _unary_func_generators = {
198 1: _HarmonicLowerBoundGenerator,
199 2: _HarmonicUpperBoundGenerator,
200 3: _HarmonicGenerator,
203 10: _SplineGenerator,
207 def _DistanceRestraintGenerator(form, modalities, atoms, parameters):
208 unary_func_gen = _unary_func_generators[form]
210 unary_func_gen(parameters, modalities),
213 def _AngleRestraintGenerator(form, modalities, atoms, parameters):
214 unary_func_gen = _unary_func_generators[form]
216 unary_func_gen(parameters, modalities),
217 atoms[0], atoms[1], atoms[2])
219 def _MultiBinormalGenerator(form, modalities, atoms, parameters):
220 nterms = modalities[0]
221 if len(parameters) != nterms * 6:
222 raise ValueError(
"Incorrect number of parameters (%d) for multiple "
223 "binormal restraint - expecting %d (%d terms * 6)" \
224 % (len(parameters), nterms * 6, nterms))
226 atoms[:4], atoms[4:8])
227 for i
in range(nterms):
229 t.set_weight(parameters[i])
230 t.set_means((parameters[nterms + i * 2],
231 parameters[nterms + i * 2 + 1]))
232 t.set_standard_deviations((parameters[nterms * 3 + i * 2],
233 parameters[nterms * 3 + i * 2 + 1]))
234 t.set_correlation(parameters[nterms * 5 + i])
238 def _DihedralRestraintGenerator(form, modalities, atoms, parameters):
240 return _MultiBinormalGenerator(form, modalities, atoms, parameters)
241 unary_func_gen = _unary_func_generators[form]
243 unary_func_gen(parameters, modalities),
244 atoms[0], atoms[1], atoms[2], atoms[3])
246 def _get_protein_atom_particles(protein):
247 """Given a protein particle, get the flattened list of all child atoms"""
249 for ichain
in range(protein.get_number_of_children()):
250 chain = protein.get_child(ichain)
251 for ires
in range(chain.get_number_of_children()):
252 residue = chain.get_child(ires)
253 for iatom
in range(residue.get_number_of_children()):
254 atom = residue.get_child(iatom)
255 atom_particles.append(atom.get_particle())
256 return atom_particles
258 def _load_restraints_line(line, atom_particles):
259 """Parse a single Modeller restraints file line and return the
260 corresponding IMP restraint."""
263 if typ ==
'MODELLER5':
266 raise NotImplementedError(
"Only 'R' lines currently read from " + \
267 "Modeller restraints files")
268 form = int(spl.pop(0))
269 modalities = [int(spl.pop(0))]
270 features = [int(spl.pop(0))]
273 natoms = [int(spl.pop(0))]
274 nparam = int(spl.pop(0))
275 nfeat = int(spl.pop(0))
276 for i
in range(nfeat - 1):
277 modalities.append(int(spl.pop(0)))
278 features.append(int(spl.pop(0)))
279 natoms.append(int(spl.pop(0)))
280 atoms = [int(spl.pop(0))
for x
in range(natoms[0])]
281 for i
in range(len(atoms)):
282 atoms[i] = atom_particles[atoms[i] - 1]
283 parameters = [float(spl.pop(0))
for x
in range(nparam)]
284 restraint_generators = {
285 1 : _DistanceRestraintGenerator,
286 2 : _AngleRestraintGenerator,
287 3 : _DihedralRestraintGenerator,
288 4 : _DihedralRestraintGenerator,
290 restraint_gen = restraint_generators[features[0]]
291 return restraint_gen(form, modalities, atoms, parameters)
294 def _load_entire_restraints_file(filename, protein):
295 """Yield a set of IMP restraints from a Modeller restraints file."""
296 atoms = _get_protein_atom_particles(protein)
297 with open(filename,
'r') as fh:
300 rsr = _load_restraints_line(line, atoms)
303 except Exception
as err:
304 print(
"Cannot read restraints file line:\n" + line)
308 def _copy_residue(r, model):
309 """Copy residue information from modeller to imp"""
314 p.set_name(str(
"residue "+r.num));
318 def _copy_atom(a, model):
319 """Copy atom information from modeller"""
325 if hasattr(a,
'charge'):
327 if hasattr(a,
'type'):
329 ap.set_input_index(a.index)
332 def _copy_chain(c, model):
333 """Copy chain information from modeller"""
340 def _get_forcefield(submodel):
352 """Add radii to the hierarchy using the Modeller radius library, radii.lib.
353 Each radius is scaled by the given scale (Modeller usually scales radii
354 by a factor of 0.82). submodel specifies the topology submodel, which is
355 the column in radii.lib to use."""
359 with open(filename)
as fh:
361 if line.startswith(
'#'):
continue
364 radii[spl[0]] = float(spl[submodel])
365 atoms = IMP.atom.get_by_type(hierarchy, IMP.atom.ATOM_TYPE)
370 radius = radii[ct] * scale
378 """Read a Modeller model into IMP. After creating this object, the atoms
379 in the Modeller model can be loaded into IMP using the load_atoms()
380 method, then optionally any Modeller static restraints can be read in
381 with load_static_restraints() or load_static_restraints_file().
383 This class can also be used to read Modeller alignment structures;
384 however, only load_atoms() will be useful in such a case (since
385 alignment structures don't have restraints or other information).
391 @param modeller_model The Modeller model or alignment structure
394 self._modeller_model = modeller_model
397 """Construct an IMP::atom::Hierarchy that contains the same atoms as
398 the Modeller model or alignment structure.
400 IMP atoms created from a Modeller model will be given charges and
401 CHARMM types, extracted from the model. Alignment structures don't
402 contain this information, so the IMP atoms won't either.
404 @param model The IMP::Model object in which the hierarchy will be
405 created. The highest level hierarchy node is a PROTEIN.
406 @return the newly-created root IMP::atom::Hierarchy.
411 for chain
in self._modeller_model.chains:
416 for residue
in chain.residues:
417 rp = _copy_residue(residue, model)
420 for atom
in residue.atoms:
421 ap = _copy_atom(atom, model)
424 self._atoms[atom.index] = ap
426 self._modeller_hierarchy = hpp
429 def _get_nonbonded_list(self, atoms, pair_filter, edat, distance):
434 if pair_filter
is None:
436 if edat.excl_local[0]:
437 pair_filter.set_bonds(list(self.
load_bonds()))
438 if edat.excl_local[1]:
440 if edat.excl_local[2]:
442 nbl.add_pair_filter(pair_filter)
446 """Load the Modeller bond topology into the IMP model. Each bond is
447 represented in IMP as an IMP::atom::Bond, with no defined length
448 or stiffness. These bonds are primarily useful as input to
449 IMP::atom::StereochemistryPairFilter, to exclude bond interactions
450 from the nonbonded list. Typically the contribution to the scoring
451 function from the bonds is included in the Modeller static restraints
452 (use load_static_restraints() or load_static_restraints_file() to
453 load these). If you want to regenerate the stereochemistry in IMP,
454 do not use these functions (as then stereochemistry scoring terms
455 and exclusions would be double-counted) and instead use the
456 IMP::atom::CHARMMTopology class.
458 You must call load_atoms() prior to using this function.
459 @see load_angles(), load_dihedrals(), load_impropers()
460 @return A generator listing all of the bonds.
462 if not hasattr(self,
'_modeller_hierarchy'):
463 raise ValueError(
"Call load_atoms() first.")
464 for (maa, mab)
in self._modeller_model.bonds:
465 pa = self._atoms[maa.index]
466 pb = self._atoms[mab.index]
476 IMP.atom.Bond.SINGLE).get_particle()
479 """Load the Modeller angle topology into the IMP model.
480 See load_bonds() for more details."""
481 return self._internal_load_angles(self._modeller_model.angles,
485 """Load the Modeller dihedral topology into the IMP model.
486 See load_bonds() for more details."""
487 return self._internal_load_angles(self._modeller_model.dihedrals,
491 """Load the Modeller improper topology into the IMP model.
492 See load_bonds() for more details."""
493 return self._internal_load_angles(self._modeller_model.impropers,
496 def _internal_load_angles(self, angles, angle_class):
497 if not hasattr(self,
'_modeller_hierarchy'):
498 raise ValueError(
"Call load_atoms() first.")
499 for modeller_atoms
in angles:
500 imp_particles = [self._atoms[x.index]
for x
in modeller_atoms]
502 a = angle_class.setup_particle(p,
504 yield a.get_particle()
507 """Convert a Modeller static restraints file into equivalent
508 IMP::Restraints. load_atoms() must have been called first to read
509 in the atoms that the restraints will act upon.
510 @param filename Name of the Modeller restraints file. The restraints
511 in this file are assumed to act upon the model read in by
512 load_atoms(); no checking is done to enforce this.
513 @return A Python generator of the newly-created IMP::Restraint
516 if not hasattr(self,
'_modeller_hierarchy'):
517 raise ValueError(
"Call load_atoms() first.")
518 return _load_entire_restraints_file(filename, self._modeller_hierarchy)
522 """Convert the current set of Modeller static restraints into equivalent
523 IMP::Restraints. load_atoms() must have been called first to read
524 in the atoms that the restraints will act upon.
525 @return A Python generator of the newly-created IMP::Restraint
528 class _RestraintGenerator(object):
529 """Simple generator wrapper"""
532 def __iter__(self, *args, **keys):
534 def close(self, *args, **keys):
535 return self._gen.close(*args, **keys)
537 return next(self._gen)
539 def send(self, *args, **keys):
540 return self._gen.send(*args, **keys)
541 def throw(self, *args, **keys):
542 return self._gen.throw(*args, **keys)
545 rsrfile = os.path.join(t.tmpdir,
'restraints.rsr')
546 self._modeller_model.restraints.write(file=rsrfile)
549 wrap = _RestraintGenerator(gen)
555 """Convert Modeller dynamic restraints into IMP::Restraint objects.
557 For each currently active Modeller dynamic restraint
558 (e.g. soft-sphere, electrostatics) an equivalent IMP::Restraint
560 load_atoms() must have been called first to read
561 in the atoms that the restraints will act upon.
563 If pair_filter is given, it is an IMP::PairFilter object to exclude
564 pairs from the nonbonded lists used by the dynamic restraints.
565 Otherwise, an IMP::atom::StereochemistryPairFilter object is created
566 to exclude Modeller bonds, angles and dihedrals, as specified by
567 edat.excl_local. (Note that this calls load_bonds(), load_angles()
568 and load_dihedrals(), so will create duplicate lists of bonds if
569 those methods are called manually as well.)
571 @note Currently only soft-sphere, electrostatic and Lennard-Jones
572 restraints are loaded.
573 @return A Python generator of the newly-created IMP::Restraint
576 if not hasattr(self,
'_modeller_hierarchy'):
577 raise ValueError(
"Call load_atoms() first.")
578 edat = self._modeller_model.env.edat
579 libs = self._modeller_model.env.libs
581 m = atoms[0].get_model()
584 if edat.dynamic_sphere:
587 nbl = self._get_nonbonded_list(atoms, pair_filter, edat, 0.)
590 libs.topology.submodel, edat.radii_factor)
597 if edat.dynamic_lennard
or edat.dynamic_coulomb:
599 d = max(edat.contact_shell - 3.0, 0.0)
600 nbl = self._get_nonbonded_list(atoms, pair_filter, edat, d)
601 ff = _get_forcefield(libs.topology.submodel)
602 ff.add_radii(self._modeller_hierarchy)
604 if edat.dynamic_lennard:
605 ff.add_well_depths(self._modeller_hierarchy)
607 edat.lennard_jones_switch[1])
611 if edat.dynamic_coulomb:
613 edat.coulomb_switch[1])
615 ps.set_relative_dielectric(edat.relative_dielectric)
619 __version__ =
"2.21.0"
622 '''Return the version of this module, as a string'''
626 '''Return the fully-qualified name of this module'''
627 return "IMP::modeller"
630 '''Return the full path to one of this module's data files'''
632 return IMP._get_module_data_path(
"modeller", fname)
635 '''Return the full path to one of this module's example files'''
637 return IMP._get_module_example_path(
"modeller", fname)
static CHARMMAtom setup_particle(Model *m, ParticleIndex pi, String charmm_type)
Create a decorator with the passed CHARMM type.
static Charged setup_particle(Model *m, ParticleIndex pi, Float charge)
Lower bound harmonic function (non-zero when feature < mean)
def load_static_restraints
Convert the current set of Modeller static restraints into equivalent IMP::Restraints.
def get_module_version
Return the version of this module, as a string.
def load_bonds
Load the Modeller bond topology into the IMP model.
A Modeller restraint which evaluates an IMP scoring function.
Coulomb (electrostatic) score between a pair of particles.
static Atom setup_particle(Model *m, ParticleIndex pi, Atom other)
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
static XYZR setup_particle(Model *m, ParticleIndex pi)
A decorator for a particle which has bonds.
def load_impropers
Load the Modeller improper topology into the IMP model.
virtual double unprotected_evaluate(DerivativeAccumulator *da) const
Return the unweighted score for the restraint.
def get_example_path
Return the full path to one of this module's example files.
Dihedral restraint between four particles.
def load_static_restraints_file
Convert a Modeller static restraints file into equivalent IMP::Restraints.
A score on the distance between the surfaces of two spheres.
Return all close unordered pairs of particles taken from the SingletonContainer.
static Residue setup_particle(Model *m, ParticleIndex pi, ResidueType t, int index, int insertion_code)
A single binormal term in a MultipleBinormalRestraint.
def load_dynamic_restraints
Convert Modeller dynamic restraints into IMP::Restraint objects.
Distance restraint between two particles.
static XYZ setup_particle(Model *m, ParticleIndex pi)
An IMP restraint using all defined Modeller restraints.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
CHARMM force field parameters.
Bond create_bond(Bonded a, Bonded b, Bond o)
Connect the two wrapped particles by a custom bond.
def load_dihedrals
Load the Modeller dihedral topology into the IMP model.
static Float get_k_from_standard_deviation(Float sd, Float t=297.15)
Return the k to use for a given Gaussian standard deviation.
def get_data_path
Return the full path to one of this module's data files.
Angle restraint between three particles.
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
Get the indexes from a list of particle pairs.
Lennard-Jones score between a pair of particles.
A particle that describes an angle between three particles.
The standard decorator for manipulating molecular structures.
Store a list of ParticleIndexes.
static Bonded setup_particle(Model *m, ParticleIndex pi)
def load_atoms
Construct an IMP::atom::Hierarchy that contains the same atoms as the Modeller model or alignment str...
Version and module information for Objects.
A decorator for a particle with x,y,z coordinates.
Modeller-style multiple binormal (phi/psi) restraint.
def add_soft_sphere_radii
Add radii to the hierarchy using the Modeller radius library, radii.lib.
A particle that describes a dihedral angle between four particles.
def load_angles
Load the Modeller angle topology into the IMP model.
def get_module_name
Return the fully-qualified name of this module.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
virtual VersionInfo get_version_info() const
Get information about the module and version of the object.
A filter that excludes bonds, angles and dihedrals.
Read a Modeller model into IMP.
Class to handle individual particles of a Model object.
std::string get_data_path(std::string file_name)
Return the full path to one of this module's data files.
Smooth interaction scores by switching the derivatives (force switch).
Functionality for loading, creating, manipulating and scoring atomic structures.
static Chain setup_particle(Model *m, ParticleIndex pi, std::string id)
Hierarchies get_leaves(const Selection &h)
Applies a PairScore to each Pair in a list.
virtual ModelObjectsTemp do_get_inputs() const =0
Closed cubic spline function.
A decorator for an atom that has a defined CHARMM type.
A restraint is a term in an IMP ScoringFunction.
Harmonic function (symmetric about the mean)
A decorator for a particle with x,y,z coordinates and a radius.