1 """@namespace IMP.pmi.restraints.stereochemistry
2 Restraints for keeping correct stereochemistry.
5 from __future__
import print_function
12 from operator
import itemgetter
17 """Create a restraint between consecutive TempResidue objects
18 or an entire PMI Molecule object."""
23 disorderedlength=
False,
28 @param objects - a list of hierarchies, PMI TempResidues OR a
30 @param scale Scale the maximal distance between the beads by this
31 factor when disorderedlength is False. The maximal distance
32 is calculated as ((float(residuegap) + 1.0) * 3.6) * scale.
33 @param disorderedlength - This flag uses either disordered length
34 calculated for random coil peptides (True) or zero
35 surface-to-surface distance between beads (False)
36 as optimal distance for the sequence connectivity restraint.
37 @param upperharmonic - This flag uses either harmonic (False)
38 or upperharmonic (True) in the intra-pair
39 connectivity restraint.
40 @param resolution - The resolution to connect things at - only used
41 if you pass PMI objects
42 @param label - A string to identify this restraint in the
48 raise Exception(
"ConnectivityRestraint: only pass stuff from "
49 "one Molecule, please")
51 m = list(hiers)[0].get_model()
52 super(ConnectivityRestraint, self).
__init__(m, label=label)
68 SortedSegments.append((start, end, startres))
69 SortedSegments = sorted(SortedSegments, key=itemgetter(2))
72 self.particle_pairs = []
73 for x
in range(len(SortedSegments) - 1):
75 last = SortedSegments[x][1]
76 first = SortedSegments[x + 1][0]
78 apply_restraint =
True
95 apply_restraint =
False
104 residuegap = firstresn - lastresn - 1
105 if disorderedlength
and (nreslast / 2 + nresfirst / 2
106 + residuegap) > 20.0:
109 optdist = math.sqrt(5 / 3) * 1.93 * \
110 (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
117 optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
124 pt0 = last.get_particle()
125 pt1 = first.get_particle()
126 self.particle_pairs.append((pt0, pt1))
128 self.model, dps, (pt0.get_index(), pt1.get_index()))
130 print(
"Adding sequence connectivity restraint between",
131 pt0.get_name(),
" and ", pt1.get_name(),
'of distance',
133 self.rs.add_restraint(r)
136 """ Returns number of connectivity restraints """
137 return len(self.rs.get_restraints())
140 """ Returns the list of connected particles pairs """
141 return self.particle_pairs
145 """A class to create an excluded volume restraint for a set of particles
146 at a given resolution.
147 Can be initialized as a bipartite restraint between two sets of particles.
148 # Potential additional function: Variable resolution for each PMI object.
149 Perhaps passing selection_tuples with (PMI_object, resolution)
151 _include_in_rmf =
True
160 @param included_objects Can be one of the following inputs:
161 IMP Hierarchy, PMI System/State/Molecule/TempResidue,
162 or a list/set of them
163 @param other_objects Initializes a bipartite restraint between
164 included_objects and other_objects
165 Same format as included_objects
166 @param resolution The resolution particles at which to impose the
167 restraint. By default, the coarsest particles will be chosen.
168 If a number is chosen, for each particle, the closest
169 resolution will be used (see IMP.atom.Selection).
170 @param kappa Restraint strength
182 if other_objects
is not None:
190 if hierarchies
is None:
191 raise Exception(
"Must at least pass included objects")
192 mdl = hierarchies[0].get_model()
193 super(ExcludedVolumeSphere, self).
__init__(mdl, label=label)
195 included_ps = [h.get_particle()
for h
in hierarchies]
197 other_ps = [h.get_particle()
for h
in other_hierarchies]
219 self.rs.add_restraint(evr)
221 def add_excluded_particle_pairs(self, excluded_particle_pairs):
223 inverted = [(p1, p0)
for p0, p1
in excluded_particle_pairs]
228 self.cpc.add_pair_filter(icpf)
232 """Enforce ideal Helix dihedrals and bonds for a selection
240 @param hierarchy the root node
241 @param selection_tuple (start, stop, molname, copynum=0)
244 m = hierarchy.get_model()
245 super(HelixRestraint, self).
__init__(m, weight=weight)
246 start = selection_tuple[0]
247 stop = selection_tuple[1]
248 mol = selection_tuple[2]
250 if len(selection_tuple) > 3:
251 copy_index = selection_tuple[3]
254 hierarchy, molecule=mol, copy_index=copy_index,
255 residue_indexes=range(start, stop+1))
256 ps = sel.get_selected_particles(with_representation=
False)
259 self.rs.add_restraint(self.r)
260 print(
'Created helix %s.%i.%i-%i with %i dihedrals and %i bonds'
261 % (mol, copy_index, start, stop, self.get_number_of_bonds(),
262 self.get_number_of_dihedrals()))
264 def get_number_of_bonds(self):
265 return self.r.get_number_of_bonds()
267 def get_number_of_dihedrals(self):
268 return self.r.get_number_of_dihedrals()
272 """ Add bond restraint between pair of consecutive
273 residues/beads to enforce the stereochemistry.
275 def __init__(self, objects, distance=3.78, strength=10.0, jitter=None):
277 @param objects Objects to restrain
278 @param distance Resting distance for restraint
279 @param strength Bond constant
280 @param jitter Defines the +- added to the optimal distance in
281 the harmonic well restraint used to increase the tolerance
285 m = particles[0].get_model()
286 super(ResidueBondRestraint, self).
__init__(m)
294 (distance - jitter, distance + jitter), strength)
299 raise ValueError(
"wrong length of pair")
302 raise TypeError(
"%s is not a residue" % p)
305 print(
"ResidueBondRestraint: adding a restraint between %s %s"
306 % (pair[0].get_name(), pair[1].get_name()))
307 self.rs.add_restraint(
312 def get_excluded_pairs(self):
313 return self.pairslist
317 """Add angular restraint between triplets of consecutive
318 residues/beads to enforce the stereochemistry.
320 def __init__(self, objects, anglemin=100.0, anglemax=140.0, strength=10.0):
323 m = particles[0].get_model()
324 super(ResidueAngleRestraint, self).__init__(m)
329 (math.pi * anglemin / 180.0,
330 math.pi * anglemax / 180.0),
336 raise ValueError(
"wrong length of triplet")
339 raise TypeError(
"%s is not a residue" % p)
342 print(
"ResidueAngleRestraint: adding a restraint between %s %s %s"
343 % (triplet[0].get_name(), triplet[1].get_name(),
344 triplet[2].get_name()))
345 self.rs.add_restraint(
353 def get_excluded_pairs(self):
354 return self.pairslist
358 """Add dihedral restraints between quadruplet of consecutive
359 residues/beads to enforce the stereochemistry.
360 Give as input a string of "C" and "T", meaning cys (0+-40)
361 or trans (180+-40) dihedral. The length of the string must be #residue-3.
362 Without the string, the dihedral will be assumed trans.
364 def __init__(self, objects, stringsequence=None, strength=10.0):
367 m = particles[0].get_model()
368 super(ResidueDihedralRestraint, self).__init__(m)
372 if stringsequence
is None:
373 stringsequence =
"T" * (len(particles) - 3)
375 for n, ps
in enumerate(
379 raise ValueError(
"wrong length of quadruplet")
382 raise TypeError(
"%s is not a residue" % p)
385 dihedraltype = stringsequence[n]
386 if dihedraltype ==
"C":
390 (math.pi * anglemin / 180.0,
391 math.pi * anglemax / 180.0),
393 print(
"ResidueDihedralRestraint: adding a CYS restraint "
394 "between %s %s %s %s"
395 % (quadruplet[0].get_name(), quadruplet[1].get_name(),
396 quadruplet[2].get_name(), quadruplet[3].get_name()))
397 if dihedraltype ==
"T":
398 anglemin = 180 - 70.0
399 anglemax = 180 + 70.0
401 (math.pi * anglemin / 180.0,
402 math.pi * anglemax / 180.0),
404 print(
"ResidueDihedralRestraint: adding a TRANS restraint "
405 "between %s %s %s %s"
406 % (quadruplet[0].get_name(),
407 quadruplet[1].get_name(), quadruplet[2].get_name(),
408 quadruplet[3].get_name()))
409 self.rs.add_restraint(
415 self.pairslist.append(
417 self.pairslist.append(
422 """Add harmonic restraints between all pairs
424 def __init__(self, hierarchy, selection_tuples=None, resolution=1,
425 strength=0.01, dist_cutoff=10.0, ca_only=
True):
427 @param hierarchy Root hierarchy to select from
428 @param selection_tuples Selecting regions for the restraint
429 [[start,stop,molname,copy_index=0],...]
430 @param resolution Resolution for applying restraint
431 @param strength Bond strength
432 @param dist_cutoff Cutoff for making restraints
433 @param ca_only Selects only CAlphas. Only matters if resolution=0.
437 self.m = hierarchy.get_model()
438 for st
in selection_tuples:
444 hierarchy, molecule=st[2],
445 residue_indexes=range(st[0], st[1]+1),
446 copy_index=copy_index)
449 hierarchy, molecule=st[2],
450 residue_indexes=range(st[0], st[1]+1),
451 copy_index=copy_index,
453 particles += sel.get_selected_particles()
461 particles, dist_cutoff, strength)
462 for r
in self.rs.get_restraints():
463 a1, a2 = r.get_inputs()
466 print(
'ElasticNetwork: created', self.rs.get_number_of_restraints(),
469 def set_label(self, label):
471 self.rs.set_name(label)
472 for r
in self.rs.get_restraints():
475 def add_to_model(self):
478 def get_restraint(self):
481 def set_weight(self, weight):
483 self.rs.set_weight(weight)
485 def get_excluded_pairs(self):
486 return self.pairslist
488 def get_output(self):
490 score = self.weight * self.rs.unprotected_evaluate(
None)
491 output[
"_TotalScore"] = str(score)
492 output[
"ElasticNetworkRestraint_" + self.label] = str(score)
497 """ Enable CHARMM force field """
498 def __init__(self, root, ff_temp=300.0, zone_ps=None, zone_size=10.0,
499 enable_nonbonded=
True, enable_bonded=
True,
500 zone_nonbonded=
False):
501 """Setup the CHARMM restraint on a selection. Expecting atoms.
502 @param root The node at which to apply the restraint
503 @param ff_temp The temperature of the force field
504 @param zone_ps Create a zone around this set of particles
505 Automatically includes the entire residue (incl. backbone)
506 @param zone_size The size for looking for neighbor residues
507 @param enable_nonbonded Allow the repulsive restraint
508 @param enable_bonded Allow the bonded restraint
509 @param zone_nonbonded EXPERIMENTAL: exclude from nonbonded all
510 sidechains that aren't in zone!
513 kB = (1.381 * 6.02214) / 4184.0
515 self.mdl = root.get_model()
526 topology = ff.create_topology(root)
527 topology.apply_default_patches()
528 topology.setup_hierarchy(root)
529 if zone_ps
is not None:
530 limit_to_ps = IMP.pmi.topology.get_particles_within_zone(
531 root, zone_ps, zone_size, entire_residues=
True,
532 exclude_backbone=
False)
536 self.ps = limit_to_ps
540 print(
'init bonds score', r.unprotected_evaluate(
None))
541 self.bonds_rs.add_restraint(r)
543 ff.add_well_depths(root)
545 atoms = IMP.atom.get_by_type(root, IMP.atom.ATOM_TYPE)
548 if (zone_ps
is not None)
and zone_nonbonded:
549 print(
'stereochemistry: zone_nonbonded is True')
551 backbone_types = [
'C',
'N',
'CB',
'O']
554 for n
in backbone_types])
555 backbone_atoms = sel.get_selected_particles()
556 sel_ps = IMP.pmi.topology.get_particles_within_zone(
557 root, zone_ps, zone_size, entire_residues=
True,
558 exclude_backbone=
True)
568 self.nbl.add_pair_filter(r.get_full_pair_filter())
571 self.nonbonded_rs.add_restraint(pr)
572 print(
'CHARMM is set up')
574 def set_label(self, label):
576 self.rs.set_name(label)
577 for r
in self.rs.get_restraints():
580 def add_to_model(self):
584 def get_restraint(self):
587 def get_close_pair_container(self):
590 def set_weight(self, weight):
592 self.rs.set_weight(weight)
594 def get_output(self):
596 bonds_score = self.weight * self.bonds_rs.unprotected_evaluate(
None)
598 self.weight * self.nonbonded_rs.unprotected_evaluate(
None)
599 score = bonds_score+nonbonded_score
600 output[
"_TotalScore"] = str(score)
601 output[
"CHARMM_BONDS"] = str(bonds_score)
602 output[
"CHARMM_NONBONDED"] = str(nonbonded_score)
607 """Add bonds and improper dihedral restraints for the CBs
610 self, rnums, representation, selection_tuple, strength=10.0, kappa=1.0,
611 jitter_angle=0.0, jitter_improper=0.0):
615 ca-cb is a constraint, no restraint needed
620 self.m = representation.prot.get_model()
629 ca, cb = self.get_ca_cb(
630 IMP.pmi.tools.select_by_tuple(representation,
631 (rnum, rnum,
'chainA'),
636 ca_prev, cb_prev = self.get_ca_cb(
637 IMP.pmi.tools.select_by_tuple(representation,
639 'chainA'), resolution=0))
640 ca_next, cb_next = self.get_ca_cb(
641 IMP.pmi.tools.select_by_tuple(representation,
643 'chainA'), resolution=0))
677 self.rset_angles.add_restraint(ar13u)
678 self.rset_angles.add_restraint(ar13l)
684 self.rset_angles.add_restraint(ar23u)
685 self.rset_angles.add_restraint(ar23l)
686 if not nter
and not cter:
710 self.rset_angles.add_restraint(idru)
711 self.rset_angles.add_restraint(idrl)
712 self.rs.add_restraint(self.rset_bonds)
713 self.rs.add_restraint(self.rset_angles)
715 def get_ca_cb(self, atoms):
720 ca = a.get_particle()
722 cb = a.get_particle()
725 def set_label(self, label):
727 self.rs.set_name(label)
728 for r
in self.rs.get_restraints():
731 def add_to_model(self):
734 def get_restraint(self):
737 def set_weight(self, weight):
739 self.rs.set_weight(weight)
741 def get_excluded_pairs(self):
742 return self.pairslist
744 def get_output(self):
746 score = self.weight * self.rs.unprotected_evaluate(
None)
747 output[
"_TotalScore"] = str(score)
748 output[
"PseudoAtomicRestraint_" + self.label] = str(score)
753 """Create harmonic restraints between the reference and (transformed)
756 @note Wraps IMP::core::TransformedDistancePairScore with an
759 def __init__(self, references, clones_list, transforms,
760 label=
'', strength=10.0, ca_only=
False):
762 @param references Can be one of the following inputs:
763 IMP Hierarchy, PMI System/State/Molecule/TempResidue,
764 or a list/set of them
765 @param clones_list List of lists of the above
766 @param transforms Transforms moving each selection to the first
768 @param label Label for output
769 @param strength The elastic bond strength
770 @param ca_only Optionally select so only CAlpha particles are used
774 self.mdl = refs[0].get_model()
778 if len(clones_list) != len(transforms):
780 'Error: There should be as many clones as transforms')
783 for tmp_clones, trans
in zip(clones_list, transforms):
785 if len(clones) != len(refs):
786 raise Exception(
"Error: len(references)!=len(clones)")
788 for p0, p1
in zip(refs, clones):
795 self.mdl, pair_score, [p0.get_particle_index(),
796 p1.get_particle_index()])
797 self.rs.add_restraint(r)
798 print(
'created symmetry network with',
799 self.rs.get_number_of_restraints(),
'restraints')
801 def set_label(self, label):
803 self.rs.set_name(label)
804 for r
in self.rs.get_restraints():
807 def add_to_model(self):
810 def get_restraint(self):
813 def set_weight(self, weight):
815 self.rs.set_weight(weight)
817 def get_excluded_pairs(self):
818 return self.pairslist
820 def get_output(self):
822 score = self.weight * self.rs.unprotected_evaluate(
None)
823 output[
"SymmetryRestraint_" + self.label] = str(score)
824 output[
"_TotalScore"] = str(score)
829 """Creates a restraint between the termini two polypeptides, to simulate
830 the sequence connectivity."""
831 def __init__(self, nterminal, cterminal, scale=1.0, disorderedlength=False,
832 upperharmonic=
True, resolution=1, label=
"None"):
834 @param nterminal - single PMI2 Hierarchy/molecule at the nterminal
835 @param cterminal - single PMI2 Hierarchy/molecule at the cterminal
836 @param scale Scale the maximal distance between the beads by this
837 factor when disorderedlength is False.
838 The maximal distance is calculated as
839 ((float(residuegap) + 1.0) * 3.6) * scale.
840 @param disorderedlength - This flag uses either disordered length
841 calculated for random coil peptides (True) or zero
842 surface-to-surface distance between beads (False)
843 as optimal distance for the sequence connectivity
845 @param upperharmonic - This flag uses either harmonic (False)
846 or upperharmonic (True) in the intra-pair
847 connectivity restraint.
848 @param resolution - The resolution to connect things at - only used
849 if you pass PMI objects
850 @param label - A string to identify this restraint in the
857 nter_lastres = ssn[-1][1]
858 cter_firstres = ssc[0][0]
859 self.m = nter_lastres.get_model()
863 optdist = (3.6) * scale
870 pt0 = nter_lastres.get_particle()
871 pt1 = cter_firstres.get_particle()
873 (pt0.get_index(), pt1.get_index()))
875 print(
"Adding fusion connectivity restraint between", pt0.get_name(),
876 " and ", pt1.get_name(),
'of distance', optdist)
877 self.rs.add_restraint(r)
879 def set_label(self, label):
882 def get_weight(self):
885 def add_to_model(self):
888 def get_restraint(self):
891 def set_weight(self, weight):
893 self.rs.set_weight(weight)
895 def get_output(self):
897 score = self.weight * self.rs.unprotected_evaluate(
None)
898 output[
"_TotalScore"] = str(score)
899 output[
"FusionRestraint_" + self.label] = str(score)
903 return self.weight * self.rs.unprotected_evaluate(
None)
908 """Restrain the dihedral between planes defined by three particles.
910 This restraint is useful for restraining the twist of a string of
911 more or less identical rigid bodies, so long as the curvature is mild.
914 def __init__(self, particle_triplets, angle=0., k=1., label=None,
917 @param particle_triplets List of lists of 3 particles. Each triplet
918 defines a plane. Dihedrals of adjacent planes
920 @param angle Angle of plane dihedral in degrees
921 @param k Strength of restraint
922 @param label Label for output
923 @param weight Weight of restraint
924 @note Particles defining planes should be rigid and more or less
925 parallel for proper behavior
927 model = particle_triplets[0][0].get_model()
928 super(PlaneDihedralRestraint, self).
__init__(model, label=label,
931 angle = math.pi * angle / 180.
933 for i, t1
in enumerate(particle_triplets[:-1]):
934 t2 = particle_triplets[i + 1]
935 q1 = [t1[1], t1[0], t2[0], t2[1]]
936 q2 = [t1[2], t1[0], t2[0], t2[2]]
937 self.rs.add_restraint(
939 self.rs.add_restraint(
A filter which returns true if a container contains the Pair.
Add dihedral restraints between quadruplet of consecutive residues/beads to enforce the stereochemist...
CHARMMParameters * get_heavy_atom_CHARMM_parameters()
Lower bound harmonic function (non-zero when feature < mean)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Enforce CHARMM stereochemistry on the given Hierarchy.
A member of a rigid body, it has internal (local) coordinates.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
Enable CHARMM force field.
A class to store a fixed array of same-typed values.
Enforce ideal Helix dihedrals and bonds for a selection at resolution 0.
Creates a restraint between the termini two polypeptides, to simulate the sequence connectivity...
Add harmonic restraints between all pairs.
Restrain the dihedral between planes defined by three particles.
Dihedral restraint between four particles.
A score on the distance between the surfaces of two spheres.
Return all close unordered pairs of particles taken from the SingletonContainer.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Return all spatially-proximal pairs of particles (a,b) from the two SingletonContainers A and B...
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
Distance restraint between two particles.
Object used to hold a set of restraints.
Store a list of ParticleIndexPairs.
A well with harmonic barriers.
Angle restraint between three particles.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
Get the indexes from a list of particle pairs.
Add bond restraint between pair of consecutive residues/beads to enforce the stereochemistry.
The standard decorator for manipulating molecular structures.
RestraintSet * create_elastic_network(const Particles &ps, Float dist_cutoff, Float strength)
Create an elastic network restraint set.
Add bonds and improper dihedral restraints for the CBs.
Store a list of ParticleIndexes.
A decorator for a particle representing an atom.
def get_particle_pairs
Returns the list of connected particles pairs.
def __init__
Setup the CHARMM restraint on a selection.
Create harmonic restraints between the reference and (transformed) clones.
Add angular restraint between triplets of consecutive residues/beads to enforce the stereochemistry...
def get_num_restraints
Returns number of connectivity restraints.
Score a pair of particles based on the distance between their centers.
A decorator for a residue.
Basic functionality that is expected to be used by a wide variety of IMP users.
Create a restraint between consecutive TempResidue objects or an entire PMI Molecule object...
A class to create an excluded volume restraint for a set of particles at a given resolution.
The general base class for IMP exceptions.
def __init__
need to add: ca-ca bond ca-cb is a constraint, no restraint needed ca-ca-ca cb-ca-ca-cb ...
Applies a PairScore to a Pair.
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Applies a PairScore to each Pair in a list.
A repulsive potential on the distance between two atoms.
Perform more efficient close pair finding when rigid bodies are involved.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
Harmonic function (symmetric about the mean)
Restraint a set of residues to use ideal helix dihedrals and bonds.