IMP logo
IMP Reference Guide  2.16.0
The Integrative Modeling Platform
restraints/stereochemistry.py
1 """@namespace IMP.pmi.restraints.stereochemistry
2 Restraints for keeping correct stereochemistry.
3 """
4 
5 from __future__ import print_function
6 import IMP
7 import IMP.core
8 import IMP.atom
9 import IMP.container
10 import IMP.isd
11 import IMP.pmi.tools
12 from operator import itemgetter
13 import math
14 
15 
16 class ConnectivityRestraint(object):
17  """Create a restraint between consecutive TempResidue objects
18  or an entire PMI Molecule object."""
19 
20  def __init__(self,
21  objects,
22  scale=1.0,
23  disorderedlength=False,
24  upperharmonic=True,
25  resolution=1,
26  label="None"):
27  """
28  @param objects - a list of hierarchies, PMI TempResidues OR a
29  single Molecule
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
43  output/stat file
44  """
45  self.label = label
46  self.weight = 1.0
47 
48  hiers = IMP.pmi.tools.input_adaptor(objects, resolution)
49  if len(hiers) > 1:
50  raise Exception("ConnectivityRestraint: only pass stuff from "
51  "one Molecule, please")
52  hiers = hiers[0]
53 
54  self.kappa = 10 # spring constant used for the harmonic restraints
55  self.m = list(hiers)[0].get_model()
56  SortedSegments = []
57  self.rs = IMP.RestraintSet(self.m, "connectivity_restraint")
58  for h in hiers:
59  try:
60  start = IMP.atom.Hierarchy(h).get_children()[0]
61  except: # noqa: E722
62  start = IMP.atom.Hierarchy(h)
63 
64  try:
65  end = IMP.atom.Hierarchy(h).get_children()[-1]
66  except: # noqa: E722
67  end = IMP.atom.Hierarchy(h)
68 
69  startres = IMP.pmi.tools.get_residue_indexes(start)[0]
70  SortedSegments.append((start, end, startres))
71  SortedSegments = sorted(SortedSegments, key=itemgetter(2))
72 
73  # connect the particles
74  self.particle_pairs = []
75  for x in range(len(SortedSegments) - 1):
76 
77  last = SortedSegments[x][1]
78  first = SortedSegments[x + 1][0]
79 
80  apply_restraint = True
81 
82  # Apply connectivity runless ALL of the following are true:
83  # - first and last both have RigidBodyMember decorators
84  # - first and last are both RigidMembers
85  # - first and last are part of the same RigidBody object
86 
87  # Check for both in a rigid body
92  # Check if the rigid body objects for each particle are
93  # the same object.
94  # if so, skip connectivity restraint
95  if IMP.core.RigidBodyMember(first).get_rigid_body() \
96  == IMP.core.RigidBodyMember(last).get_rigid_body():
97  apply_restraint = False
98 
99  if apply_restraint:
100 
101  nreslast = len(IMP.pmi.tools.get_residue_indexes(last))
102  lastresn = IMP.pmi.tools.get_residue_indexes(last)[-1]
103  nresfirst = len(IMP.pmi.tools.get_residue_indexes(first))
104  firstresn = IMP.pmi.tools.get_residue_indexes(first)[0]
105 
106  residuegap = firstresn - lastresn - 1
107  if disorderedlength and (nreslast / 2 + nresfirst / 2
108  + residuegap) > 20.0:
109  # calculate the distance between the sphere centers
110  # using Kohn PNAS 2004
111  optdist = math.sqrt(5 / 3) * 1.93 * \
112  (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
113  if upperharmonic:
114  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
115  else:
116  hu = IMP.core.Harmonic(optdist, self.kappa)
118  else: # default
119  optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
120  if upperharmonic: # default
121  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
122  else:
123  hu = IMP.core.Harmonic(optdist, self.kappa)
125 
126  pt0 = last.get_particle()
127  pt1 = first.get_particle()
128  self.particle_pairs.append((pt0, pt1))
130  self.m, dps, (pt0.get_index(), pt1.get_index()))
131 
132  print("Adding sequence connectivity restraint between",
133  pt0.get_name(), " and ", pt1.get_name(), 'of distance',
134  optdist)
135  self.rs.add_restraint(r)
136 
137  def set_label(self, label):
138  self.label = label
139 
140  def get_weight(self):
141  return self.weight
142 
143  def add_to_model(self):
144  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
145 
146  def get_restraint(self):
147  return self.rs
148 
149  def set_weight(self, weight):
150  self.weight = weight
151  self.rs.set_weight(weight)
152 
153  def get_output(self):
154  output = {}
155  score = self.weight * self.rs.unprotected_evaluate(None)
156  output["_TotalScore"] = str(score)
157  output["ConnectivityRestraint_" + self.label] = str(score)
158  return output
159 
161  """ Returns number of connectivity restraints """
162  return len(self.rs.get_restraints())
163 
165  """ Returns the list of connected particles pairs """
166  return self.particle_pairs
167 
168  def evaluate(self):
169  return self.weight * self.rs.unprotected_evaluate(None)
170 
171 
172 class ExcludedVolumeSphere(object):
173  """A class to create an excluded volume restraint for a set of particles
174  at a given resolution.
175  Can be initialized as a bipartite restraint between two sets of particles.
176  # Potential additional function: Variable resolution for each PMI object.
177  Perhaps passing selection_tuples with (PMI_object, resolution)
178  """
179 
180  def __init__(self,
181  included_objects,
182  other_objects=None,
183  resolution=1000,
184  kappa=1.0):
185  """Constructor.
186  @param included_objects Can be one of the following inputs:
187  IMP Hierarchy, PMI System/State/Molecule/TempResidue,
188  or a list/set of them
189  @param other_objects Initializes a bipartite restraint between
190  included_objects and other_objects
191  Same format as included_objects
192  @param resolution The resolution particles at which to impose the
193  restraint. By default, the coarsest particles will be chosen.
194  If a number is chosen, for each particle, the closest
195  resolution will be used (see IMP.atom.Selection).
196  @param kappa Restraint strength
197  """
198 
199  self.weight = 1.0
200  self.kappa = kappa
201  self.label = "None"
202  self.cpc = None
203  bipartite = False
204 
205  # gather IMP hierarchies from input objects
206  hierarchies = IMP.pmi.tools.input_adaptor(included_objects,
207  resolution,
208  flatten=True)
209  included_ps = []
210  if other_objects is not None:
211  bipartite = True
212  other_hierarchies = IMP.pmi.tools.input_adaptor(other_objects,
213  resolution,
214  flatten=True)
215  other_ps = []
216 
217  # perform selection
218  if hierarchies is None:
219  raise Exception("Must at least pass included objects")
220  self.mdl = hierarchies[0].get_model()
221  included_ps = [h.get_particle() for h in hierarchies]
222  if bipartite:
223  other_ps = [h.get_particle() for h in other_hierarchies]
224 
225  # setup score
226  self.rs = IMP.RestraintSet(self.mdl, 'excluded_volume')
227  ssps = IMP.core.SoftSpherePairScore(self.kappa)
229  lsa.add(IMP.get_indexes(included_ps))
230 
231  # setup close pair container
232  if not bipartite:
234  self.cpc = IMP.container.ClosePairContainer(lsa, 0.0, rbcpf, 10.0)
235  evr = IMP.container.PairsRestraint(ssps, self.cpc)
236  else:
237  other_lsa = IMP.container.ListSingletonContainer(self.mdl)
238  other_lsa.add(IMP.get_indexes(other_ps))
240  lsa,
241  other_lsa,
242  0.0,
243  10.0)
244  evr = IMP.container.PairsRestraint(ssps, self.cpc)
245 
246  self.rs.add_restraint(evr)
247 
248  def add_excluded_particle_pairs(self, excluded_particle_pairs):
249  # add pairs to be filtered when calculating the score
250  inverted = [(p1, p0) for p0, p1 in excluded_particle_pairs]
251  lpc = IMP.container.ListPairContainer(self.mdl)
252  lpc.add(IMP.get_indexes(excluded_particle_pairs))
253  lpc.add(IMP.get_indexes(inverted))
255  self.cpc.add_pair_filter(icpf)
256 
257  def set_label(self, label):
258  self.label = label
259 
260  def add_to_model(self):
261  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.rs,
262  add_to_rmf=True)
263 
264  def get_restraint(self):
265  return self.rs
266 
267  def set_weight(self, weight):
268  self.weight = weight
269  self.rs.set_weight(weight)
270 
271  def get_output(self):
272  output = {}
273  score = self.weight * self.rs.unprotected_evaluate(None)
274  output["_TotalScore"] = str(score)
275  output["ExcludedVolumeSphere_" + self.label] = str(score)
276  return output
277 
278  def evaluate(self):
279  return self.weight * self.rs.unprotected_evaluate(None)
280 
281 
282 class HelixRestraint(object):
283  """Enforce ideal Helix dihedrals and bonds for a selection
284  at resolution 0"""
285  def __init__(self,
286  hierarchy,
287  selection_tuple,
288  weight=1.0,
289  label=''):
290  """Constructor
291  @param hierarchy the root node
292  @param selection_tuple (start, stop, molname, copynum=0)
293  @param weight
294  """
295  self.mdl = hierarchy.get_model()
296  self.rs = IMP.RestraintSet(self.mdl)
297  self.weight = weight
298  self.label = label
299  start = selection_tuple[0]
300  stop = selection_tuple[1]
301  mol = selection_tuple[2]
302  copy_index = 0
303  if len(selection_tuple) > 3:
304  copy_index = selection_tuple[3]
305 
306  sel = IMP.atom.Selection(
307  hierarchy, molecule=mol, copy_index=copy_index,
308  residue_indexes=range(start, stop+1))
309  ps = sel.get_selected_particles(with_representation=False)
310  res = [IMP.atom.Residue(p) for p in ps]
311  self.rs = IMP.RestraintSet(self.mdl, self.weight)
312  self.r = IMP.atom.HelixRestraint(res)
313  self.rs.add_restraint(self.r)
314  print('Created helix %s.%i.%i-%i with %i dihedrals and %i bonds'
315  % (mol, copy_index, start, stop, self.get_number_of_bonds(),
316  self.get_number_of_dihedrals()))
317 
318  def set_label(self, label):
319  self.label = label
320 
321  def get_weight(self):
322  return self.weight
323 
324  def add_to_model(self):
325  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.rs)
326 
327  def get_restraint(self):
328  return self.rs
329 
330  def set_weight(self, weight):
331  self.weight = weight
332  self.rs.set_weight(weight)
333 
334  def get_number_of_bonds(self):
335  return self.r.get_number_of_bonds()
336 
337  def get_number_of_dihedrals(self):
338  return self.r.get_number_of_dihedrals()
339 
340  def evaluate(self):
341  return self.weight * self.rs.unprotected_evaluate(None)
342 
343  def get_output(self):
344  output = {}
345  score = self.evaluate()
346  output["_TotalScore"] = str(score)
347  output["HelixRestraint_" + self.label] = str(score)
348  return output
349 
350 
351 class ResidueBondRestraint(object):
352  """ Add bond restraint between pair of consecutive
353  residues/beads to enforce the stereochemistry.
354  """
355  def __init__(self, objects, distance=3.78, strength=10.0, jitter=None):
356  """Constructor
357  @param objects Objects to restrain
358  @param distance Resting distance for restraint
359  @param strength Bond constant
360  @param jitter Defines the +- added to the optimal distance in
361  the harmonic well restraint used to increase the tolerance
362  """
363 
364  particles = IMP.pmi.tools.input_adaptor(objects, 1, flatten=True)
365  self.m = particles[0].get_model()
366 
367  self.rs = IMP.RestraintSet(self.m, "Bonds")
368  self.weight = 1
369  self.label = "None"
370  self.pairslist = []
371 
372  if not jitter:
373  ts = IMP.core.Harmonic(distance, strength)
374  else:
376  (distance - jitter, distance + jitter), strength)
377 
378  for ps in IMP.pmi.tools.sublist_iterator(particles, 2, 2):
379  pair = []
380  if len(ps) != 2:
381  raise ValueError("wrong length of pair")
382  for p in ps:
384  raise TypeError("%s is not a residue" % p)
385  else:
386  pair.append(p)
387  print("ResidueBondRestraint: adding a restraint between %s %s"
388  % (pair[0].get_name(), pair[1].get_name()))
389  self.rs.add_restraint(
390  IMP.core.DistanceRestraint(self.m, ts, pair[0], pair[1]))
391  self.pairslist.append(IMP.ParticlePair(pair[0], pair[1]))
392  self.pairslist.append(IMP.ParticlePair(pair[1], pair[0]))
393 
394  def set_label(self, label):
395  self.label = label
396  self.rs.set_name(label)
397  for r in self.rs.get_restraints():
398  r.set_name(label)
399 
400  def add_to_model(self):
401  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
402 
403  def get_restraint(self):
404  return self.rs
405 
406  def set_weight(self, weight):
407  self.weight = weight
408  self.rs.set_weight(weight)
409 
410  def get_excluded_pairs(self):
411  return self.pairslist
412 
413  def get_output(self):
414  output = {}
415  score = self.weight * self.rs.unprotected_evaluate(None)
416  output["_TotalScore"] = str(score)
417  output["ResidueBondRestraint_" + self.label] = str(score)
418  return output
419 
420 
421 class ResidueAngleRestraint(object):
422  """Add angular restraint between triplets of consecutive
423  residues/beads to enforce the stereochemistry.
424  """
425  def __init__(self, objects, anglemin=100.0, anglemax=140.0, strength=10.0):
426 
427  particles = IMP.pmi.tools.input_adaptor(objects, 1, flatten=True)
428  self.m = particles[0].get_model()
429 
430  self.rs = IMP.RestraintSet(self.m, "Angles")
431  self.weight = 1
432  self.label = "None"
433  self.pairslist = []
434 
436  (math.pi * anglemin / 180.0,
437  math.pi * anglemax / 180.0),
438  strength)
439 
440  for ps in IMP.pmi.tools.sublist_iterator(particles, 3, 3):
441  triplet = []
442  if len(ps) != 3:
443  raise ValueError("wrong length of triplet")
444  for p in ps:
446  raise TypeError("%s is not a residue" % p)
447  else:
448  triplet.append(p)
449  print("ResidueAngleRestraint: adding a restraint between %s %s %s"
450  % (triplet[0].get_name(), triplet[1].get_name(),
451  triplet[2].get_name()))
452  self.rs.add_restraint(
453  IMP.core.AngleRestraint(triplet[0].get_model(), ts,
454  triplet[0],
455  triplet[1],
456  triplet[2]))
457  self.pairslist.append(IMP.ParticlePair(triplet[0], triplet[2]))
458  self.pairslist.append(IMP.ParticlePair(triplet[2], triplet[0]))
459 
460  def set_label(self, label):
461  self.label = label
462  self.rs.set_name(label)
463  for r in self.rs.get_restraints():
464  r.set_name(label)
465 
466  def add_to_model(self):
467  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
468 
469  def get_restraint(self):
470  return self.rs
471 
472  def set_weight(self, weight):
473  self.weight = weight
474  self.rs.set_weight(weight)
475 
476  def get_excluded_pairs(self):
477  return self.pairslist
478 
479  def get_output(self):
480  output = {}
481  score = self.weight * self.rs.unprotected_evaluate(None)
482  output["_TotalScore"] = str(score)
483  output["ResidueAngleRestraint_" + self.label] = str(score)
484  return output
485 
486 
488  """Add dihedral restraints between quadruplet of consecutive
489  residues/beads to enforce the stereochemistry.
490  Give as input a string of "C" and "T", meaning cys (0+-40)
491  or trans (180+-40) dihedral. The length of the string must be #residue-3.
492  Without the string, the dihedral will be assumed trans.
493  """
494  def __init__(self, objects, stringsequence=None, strength=10.0):
495 
496  particles = IMP.pmi.tools.input_adaptor(objects, 1, flatten=True)
497  self.m = particles[0].get_model()
498 
499  self.rs = IMP.RestraintSet(self.m, "Angles")
500  self.weight = 1
501  self.label = "None"
502  self.pairslist = []
503 
504  if stringsequence is None:
505  stringsequence = "T" * (len(particles) - 3)
506 
507  for n, ps in enumerate(
508  IMP.pmi.tools.sublist_iterator(particles, 4, 4)):
509  quadruplet = []
510  if len(ps) != 4:
511  raise ValueError("wrong length of quadruplet")
512  for p in ps:
514  raise TypeError("%s is not a residue" % p)
515  else:
516  quadruplet.append(p)
517  dihedraltype = stringsequence[n]
518  if dihedraltype == "C":
519  anglemin = -20.0
520  anglemax = 20.0
522  (math.pi * anglemin / 180.0,
523  math.pi * anglemax / 180.0),
524  strength)
525  print("ResidueDihedralRestraint: adding a CYS restraint "
526  "between %s %s %s %s"
527  % (quadruplet[0].get_name(), quadruplet[1].get_name(),
528  quadruplet[2].get_name(), quadruplet[3].get_name()))
529  if dihedraltype == "T":
530  anglemin = 180 - 70.0
531  anglemax = 180 + 70.0
533  (math.pi * anglemin / 180.0,
534  math.pi * anglemax / 180.0),
535  strength)
536  print("ResidueDihedralRestraint: adding a TRANS restraint "
537  "between %s %s %s %s"
538  % (quadruplet[0].get_name(),
539  quadruplet[1].get_name(), quadruplet[2].get_name(),
540  quadruplet[3].get_name()))
541  self.rs.add_restraint(
542  IMP.core.DihedralRestraint(self.m, ts,
543  quadruplet[0],
544  quadruplet[1],
545  quadruplet[2],
546  quadruplet[3]))
547  self.pairslist.append(
548  IMP.ParticlePair(quadruplet[0], quadruplet[3]))
549  self.pairslist.append(
550  IMP.ParticlePair(quadruplet[3], quadruplet[0]))
551 
552  def set_label(self, label):
553  self.label = label
554  self.rs.set_name(label)
555  for r in self.rs.get_restraints():
556  r.set_name(label)
557 
558  def add_to_model(self):
559  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
560 
561  def get_restraint(self):
562  return self.rs
563 
564  def set_weight(self, weight):
565  self.weight = weight
566  self.rs.set_weight(weight)
567 
568  def get_excluded_pairs(self):
569  return self.pairslist
570 
571  def get_output(self):
572  output = {}
573  score = self.weight * self.rs.unprotected_evaluate(None)
574  output["_TotalScore"] = str(score)
575  output["ResidueDihedralRestraint_" + self.label] = str(score)
576  return output
577 
578 
580  """Add harmonic restraints between all pairs
581  """
582  def __init__(self, hierarchy, selection_tuples=None, resolution=1,
583  strength=0.01, dist_cutoff=10.0, ca_only=True):
584  """Constructor
585  @param hierarchy Root hierarchy to select from
586  @param selection_tuples Selecting regions for the restraint
587  [[start,stop,molname,copy_index=0],...]
588  @param resolution Resolution for applying restraint
589  @param strength Bond strength
590  @param dist_cutoff Cutoff for making restraints
591  @param ca_only Selects only CAlphas. Only matters if resolution=0.
592  """
593 
594  particles = []
595  self.m = hierarchy.get_model()
596  for st in selection_tuples:
597  copy_index = 0
598  if len(st) > 3:
599  copy_index = st[3]
600  if not ca_only:
601  sel = IMP.atom.Selection(
602  hierarchy, molecule=st[2],
603  residue_indexes=range(st[0], st[1]+1),
604  copy_index=copy_index)
605  else:
606  sel = IMP.atom.Selection(
607  hierarchy, molecule=st[2],
608  residue_indexes=range(st[0], st[1]+1),
609  copy_index=copy_index,
610  atom_type=IMP.atom.AtomType("CA"))
611  particles += sel.get_selected_particles()
612 
613  self.weight = 1
614  self.label = "None"
615  self.pairslist = []
616 
617  # create score
619  particles, dist_cutoff, strength)
620  for r in self.rs.get_restraints():
621  a1, a2 = r.get_inputs()
622  self.pairslist.append(IMP.ParticlePair(a1, a2))
623  self.pairslist.append(IMP.ParticlePair(a2, a1))
624  print('ElasticNetwork: created', self.rs.get_number_of_restraints(),
625  'restraints')
626 
627  def set_label(self, label):
628  self.label = label
629  self.rs.set_name(label)
630  for r in self.rs.get_restraints():
631  r.set_name(label)
632 
633  def add_to_model(self):
634  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
635 
636  def get_restraint(self):
637  return self.rs
638 
639  def set_weight(self, weight):
640  self.weight = weight
641  self.rs.set_weight(weight)
642 
643  def get_excluded_pairs(self):
644  return self.pairslist
645 
646  def get_output(self):
647  output = {}
648  score = self.weight * self.rs.unprotected_evaluate(None)
649  output["_TotalScore"] = str(score)
650  output["ElasticNetworkRestraint_" + self.label] = str(score)
651  return output
652 
653 
655  """ Enable CHARMM force field """
656  def __init__(self, root, ff_temp=300.0, zone_ps=None, zone_size=10.0,
657  enable_nonbonded=True, enable_bonded=True,
658  zone_nonbonded=False):
659  """Setup the CHARMM restraint on a selection. Expecting atoms.
660  @param root The node at which to apply the restraint
661  @param ff_temp The temperature of the force field
662  @param zone_ps Create a zone around this set of particles
663  Automatically includes the entire residue (incl. backbone)
664  @param zone_size The size for looking for neighbor residues
665  @param enable_nonbonded Allow the repulsive restraint
666  @param enable_bonded Allow the bonded restraint
667  @param zone_nonbonded EXPERIMENTAL: exclude from nonbonded all
668  sidechains that aren't in zone!
669  """
670 
671  kB = (1.381 * 6.02214) / 4184.0
672 
673  self.mdl = root.get_model()
674  self.bonds_rs = IMP.RestraintSet(self.mdl, 1.0 / (kB * ff_temp),
675  'BONDED')
676  self.nonbonded_rs = IMP.RestraintSet(self.mdl, 1.0 / (kB * ff_temp),
677  'NONBONDED')
678  self.weight = 1.0
679  self.label = ""
680 
681  # setup topology and bonds etc
682  if enable_bonded:
684  topology = ff.create_topology(root)
685  topology.apply_default_patches()
686  topology.setup_hierarchy(root)
687  if zone_ps is not None:
688  limit_to_ps = IMP.pmi.topology.get_particles_within_zone(
689  root, zone_ps, zone_size, entire_residues=True,
690  exclude_backbone=False)
692  topology,
693  limit_to_ps)
694  self.ps = limit_to_ps
695  else:
696  r = IMP.atom.CHARMMStereochemistryRestraint(root, topology)
697  self.ps = IMP.core.get_leaves(root)
698  print('init bonds score', r.unprotected_evaluate(None))
699  self.bonds_rs.add_restraint(r)
700  ff.add_radii(root)
701  ff.add_well_depths(root)
702 
703  atoms = IMP.atom.get_by_type(root, IMP.atom.ATOM_TYPE)
704  # non-bonded forces
705  if enable_nonbonded:
706  if (zone_ps is not None) and zone_nonbonded:
707  print('stereochemistry: zone_nonbonded is True')
708  # atoms list should only include backbone plus zone_ps!
709  backbone_types = ['C', 'N', 'CB', 'O']
710  sel = IMP.atom.Selection(
711  root, atom_types=[IMP.atom.AtomType(n)
712  for n in backbone_types])
713  backbone_atoms = sel.get_selected_particles()
714  sel_ps = IMP.pmi.topology.get_particles_within_zone(
715  root, zone_ps, zone_size, entire_residues=True,
716  exclude_backbone=True)
717 
719  IMP.container.ListSingletonContainer(backbone_atoms),
721  4.0)
722  else:
723  cont = IMP.container.ListSingletonContainer(self.mdl, atoms)
724  self.nbl = IMP.container.ClosePairContainer(cont, 4.0)
725  if enable_bonded:
726  self.nbl.add_pair_filter(r.get_full_pair_filter())
727  pairscore = IMP.isd.RepulsiveDistancePairScore(0, 1)
728  pr = IMP.container.PairsRestraint(pairscore, self.nbl)
729  self.nonbonded_rs.add_restraint(pr)
730  print('CHARMM is set up')
731 
732  def set_label(self, label):
733  self.label = label
734  self.rs.set_name(label)
735  for r in self.rs.get_restraints():
736  r.set_name(label)
737 
738  def add_to_model(self):
739  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.bonds_rs)
740  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.nonbonded_rs)
741 
742  def get_restraint(self):
743  return self.rs
744 
745  def get_close_pair_container(self):
746  return self.nbl
747 
748  def set_weight(self, weight):
749  self.weight = weight
750  self.rs.set_weight(weight)
751 
752  def get_output(self):
753  output = {}
754  bonds_score = self.weight * self.bonds_rs.unprotected_evaluate(None)
755  nonbonded_score = \
756  self.weight * self.nonbonded_rs.unprotected_evaluate(None)
757  score = bonds_score+nonbonded_score
758  output["_TotalScore"] = str(score)
759  output["CHARMM_BONDS"] = str(bonds_score)
760  output["CHARMM_NONBONDED"] = str(nonbonded_score)
761  return output
762 
763 
764 class PseudoAtomicRestraint(object):
765  """Add bonds and improper dihedral restraints for the CBs
766  """
767  def __init__(
768  self, rnums, representation, selection_tuple, strength=10.0, kappa=1.0,
769  jitter_angle=0.0, jitter_improper=0.0):
770  '''
771  need to add:
772  ca-ca bond
773  ca-cb is a constraint, no restraint needed
774  ca-ca-ca
775  cb-ca-ca-cb
776  '''
777 
778  self.m = representation.prot.get_model()
779  self.rs = IMP.RestraintSet(self.m, "PseudoAtomic")
780  self.rset_angles = IMP.RestraintSet(self.m, "PseudoAtomic_Angles")
781  self.rset_bonds = IMP.RestraintSet(self.m, "PseudoAtomic_Bonds")
782  self.weight = 1
783  self.label = "None"
784  self.pairslist = []
785 
786  for rnum in rnums:
787  ca, cb = self.get_ca_cb(
788  IMP.pmi.tools.select_by_tuple(representation,
789  (rnum, rnum, 'chainA'),
790  resolution=0))
791  if cb is not None:
792  nter = False
793  cter = False
794  ca_prev, cb_prev = self.get_ca_cb(
795  IMP.pmi.tools.select_by_tuple(representation,
796  (rnum - 1, rnum - 1,
797  'chainA'), resolution=0))
798  ca_next, cb_next = self.get_ca_cb(
799  IMP.pmi.tools.select_by_tuple(representation,
800  (rnum + 1, rnum + 1,
801  'chainA'), resolution=0))
802  if ca_prev is None:
803  nter = True
804  else:
805  if ca_next is None:
806  cter = True
807  else:
808  if (nter and cter):
809  continue
810 
811  # adding a bond restraint between CA and CB
812  # h=IMP.core.Harmonic(6.0,kappa)
813  # dps=IMP.core.DistancePairScore(h)
814  # pr=IMP.core.PairRestraint(dps,IMP.ParticlePair(ca,cb))
815  # self.pairslist.append((ca,cb))
816  # self.rset_bonds.add_restraint(pr)
817 
818  # creating improper dihedral restraint
819  # hus=IMP.core.Harmonic(2.09,kappa)
821  2.09 +
822  jitter_angle /
823  0.5,
824  kappa)
826  2.09 -
827  jitter_angle /
828  0.5,
829  kappa)
830  if not nter:
831  # ar13=IMP.core.AngleRestraint(hus,ca_prev,ca,cb)
832  # self.rset_angles.add_restraint(ar13)
833  ar13u = IMP.core.AngleRestraint(hupp, ca_prev, ca, cb)
834  ar13l = IMP.core.AngleRestraint(hlow, ca_prev, ca, cb)
835  self.rset_angles.add_restraint(ar13u)
836  self.rset_angles.add_restraint(ar13l)
837  if not cter:
838  # ar23=IMP.core.AngleRestraint(hus,ca_next,ca,cb)
839  # self.rset_angles.add_restraint(ar23)
840  ar23u = IMP.core.AngleRestraint(hupp, ca_next, ca, cb)
841  ar23l = IMP.core.AngleRestraint(hlow, ca_next, ca, cb)
842  self.rset_angles.add_restraint(ar23u)
843  self.rset_angles.add_restraint(ar23l)
844  if not nter and not cter:
845  # hus2=IMP.core.Harmonic(0,kappa)
846  # idr=IMP.core.DihedralRestraint(hus2,ca,ca_prev,ca_next,cb)
847  # self.rset_angles.add_restraint(idr)
848 
849  hus2upp = IMP.core.HarmonicUpperBound(
850  jitter_improper,
851  kappa)
852  hus2low = IMP.core.HarmonicLowerBound(
853  -
854  jitter_improper,
855  kappa)
857  hus2upp,
858  ca,
859  ca_prev,
860  ca_next,
861  cb)
863  hus2low,
864  ca,
865  ca_prev,
866  ca_next,
867  cb)
868  self.rset_angles.add_restraint(idru)
869  self.rset_angles.add_restraint(idrl)
870  self.rs.add_restraint(self.rset_bonds)
871  self.rs.add_restraint(self.rset_angles)
872 
873  def get_ca_cb(self, atoms):
874  ca = None
875  cb = None
876  for a in atoms:
877  if IMP.atom.Atom(a).get_atom_type() == IMP.atom.AtomType("CA"):
878  ca = a.get_particle()
879  elif IMP.atom.Atom(a).get_atom_type() == IMP.atom.AtomType("CB"):
880  cb = a.get_particle()
881  return ca, cb
882 
883  def set_label(self, label):
884  self.label = label
885  self.rs.set_name(label)
886  for r in self.rs.get_restraints():
887  r.set_name(label)
888 
889  def add_to_model(self):
890  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
891 
892  def get_restraint(self):
893  return self.rs
894 
895  def set_weight(self, weight):
896  self.weight = weight
897  self.rs.set_weight(weight)
898 
899  def get_excluded_pairs(self):
900  return self.pairslist
901 
902  def get_output(self):
903  output = {}
904  score = self.weight * self.rs.unprotected_evaluate(None)
905  output["_TotalScore"] = str(score)
906  output["PseudoAtomicRestraint_" + self.label] = str(score)
907  return output
908 
909 
910 class SymmetryRestraint(object):
911  """Create harmonic restraints between the reference and (transformed)
912  clones.
913 
914  @note Wraps IMP::core::TransformedDistancePairScore with an
915  IMP::core::Harmonic
916  """
917  def __init__(self, references, clones_list, transforms,
918  label='', strength=10.0, ca_only=False):
919  """Constructor
920  @param references Can be one of the following inputs:
921  IMP Hierarchy, PMI System/State/Molecule/TempResidue,
922  or a list/set of them
923  @param clones_list List of lists of the above
924  @param transforms Transforms moving each selection to the first
925  selection
926  @param label Label for output
927  @param strength The elastic bond strength
928  @param ca_only Optionally select so only CAlpha particles are used
929  """
930 
931  refs = IMP.pmi.tools.input_adaptor(references, flatten=True)
932  self.mdl = refs[0].get_model()
933  self.rs = IMP.RestraintSet(self.mdl, "Symmetry")
934  self.weight = 1
935  self.label = label
936  if len(clones_list) != len(transforms):
937  raise Exception(
938  'Error: There should be as many clones as transforms')
939 
940  harmonic = IMP.core.Harmonic(0., strength)
941  for tmp_clones, trans in zip(clones_list, transforms):
942  clones = IMP.pmi.tools.input_adaptor(tmp_clones, flatten=True)
943  if len(clones) != len(refs):
944  raise Exception("Error: len(references)!=len(clones)")
945  pair_score = IMP.core.TransformedDistancePairScore(harmonic, trans)
946  for p0, p1 in zip(refs, clones):
947  if not ca_only or (
948  IMP.atom.Atom(p0).get_atom_type()
949  == IMP.atom.AtomType("CA") and
950  IMP.atom.Atom(p1).get_atom_type()
951  == IMP.atom.AtomType("CA")):
953  self.mdl, pair_score, [p0.get_particle_index(),
954  p1.get_particle_index()])
955  self.rs.add_restraint(r)
956  print('created symmetry network with',
957  self.rs.get_number_of_restraints(), 'restraints')
958 
959  def set_label(self, label):
960  self.label = label
961  self.rs.set_name(label)
962  for r in self.rs.get_restraints():
963  r.set_name(label)
964 
965  def add_to_model(self):
966  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.rs)
967 
968  def get_restraint(self):
969  return self.rs
970 
971  def set_weight(self, weight):
972  self.weight = weight
973  self.rs.set_weight(weight)
974 
975  def get_excluded_pairs(self):
976  return self.pairslist
977 
978  def get_output(self):
979  output = {}
980  score = self.weight * self.rs.unprotected_evaluate(None)
981  output["SymmetryRestraint_" + self.label] = str(score)
982  output["_TotalScore"] = str(score)
983  return output
984 
985 
986 class FusionRestraint(object):
987  """Creates a restraint between the termini two polypeptides, to simulate
988  the sequence connectivity."""
989  def __init__(self, nterminal, cterminal, scale=1.0, disorderedlength=False,
990  upperharmonic=True, resolution=1, label="None"):
991  """
992  @param nterminal - single PMI2 Hierarchy/molecule at the nterminal
993  @param cterminal - single PMI2 Hierarchy/molecule at the cterminal
994  @param scale Scale the maximal distance between the beads by this
995  factor when disorderedlength is False.
996  The maximal distance is calculated as
997  ((float(residuegap) + 1.0) * 3.6) * scale.
998  @param disorderedlength - This flag uses either disordered length
999  calculated for random coil peptides (True) or zero
1000  surface-to-surface distance between beads (False)
1001  as optimal distance for the sequence connectivity
1002  restraint.
1003  @param upperharmonic - This flag uses either harmonic (False)
1004  or upperharmonic (True) in the intra-pair
1005  connectivity restraint.
1006  @param resolution - The resolution to connect things at - only used
1007  if you pass PMI objects
1008  @param label - A string to identify this restraint in the
1009  output/stat file
1010  """
1011  self.label = label
1012  self.weight = 1.0
1013  ssn = IMP.pmi.tools.get_sorted_segments(nterminal)
1014  ssc = IMP.pmi.tools.get_sorted_segments(cterminal)
1015  nter_lastres = ssn[-1][1]
1016  cter_firstres = ssc[0][0]
1017  self.m = nter_lastres.get_model()
1018 
1019  self.kappa = 10 # spring constant used for the harmonic restraints
1020 
1021  optdist = (3.6) * scale
1022  if upperharmonic: # default
1023  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1024  else:
1025  hu = IMP.core.Harmonic(optdist, self.kappa)
1027 
1028  pt0 = nter_lastres.get_particle()
1029  pt1 = cter_firstres.get_particle()
1030  r = IMP.core.PairRestraint(self.m, dps,
1031  (pt0.get_index(), pt1.get_index()))
1032  self.rs = IMP.RestraintSet(self.m, "fusion_restraint")
1033  print("Adding fusion connectivity restraint between", pt0.get_name(),
1034  " and ", pt1.get_name(), 'of distance', optdist)
1035  self.rs.add_restraint(r)
1036 
1037  def set_label(self, label):
1038  self.label = label
1039 
1040  def get_weight(self):
1041  return self.weight
1042 
1043  def add_to_model(self):
1044  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1045 
1046  def get_restraint(self):
1047  return self.rs
1048 
1049  def set_weight(self, weight):
1050  self.weight = weight
1051  self.rs.set_weight(weight)
1052 
1053  def get_output(self):
1054  output = {}
1055  score = self.weight * self.rs.unprotected_evaluate(None)
1056  output["_TotalScore"] = str(score)
1057  output["FusionRestraint_" + self.label] = str(score)
1058  return output
1059 
1060  def evaluate(self):
1061  return self.weight * self.rs.unprotected_evaluate(None)
1062 
1063 
1064 class PlaneDihedralRestraint(IMP.pmi.restraints.RestraintBase):
1065 
1066  """Restrain the dihedral between planes defined by three particles.
1067 
1068  This restraint is useful for restraining the twist of a string of
1069  more or less identical rigid bodies, so long as the curvature is mild.
1070  """
1071 
1072  def __init__(self, particle_triplets, angle=0., k=1., label=None,
1073  weight=1.):
1074  """Constructor
1075  @param particle_triplets List of lists of 3 particles. Each triplet
1076  defines a plane. Dihedrals of adjacent planes
1077  in list are scored.
1078  @param angle Angle of plane dihedral in degrees
1079  @param k Strength of restraint
1080  @param label Label for output
1081  @param weight Weight of restraint
1082  @note Particles defining planes should be rigid and more or less
1083  parallel for proper behavior
1084  """
1085  model = particle_triplets[0][0].get_model()
1086  super(PlaneDihedralRestraint, self).__init__(model, label=label,
1087  weight=weight)
1088 
1089  angle = math.pi * angle / 180.
1090  ds = IMP.core.Cosine(.5 * k, 1, -angle)
1091  for i, t1 in enumerate(particle_triplets[:-1]):
1092  t2 = particle_triplets[i + 1]
1093  q1 = [t1[1], t1[0], t2[0], t2[1]]
1094  q2 = [t1[2], t1[0], t2[0], t2[2]]
1095  self.rs.add_restraint(
1096  IMP.core.DihedralRestraint(self.model, ds, *q1))
1097  self.rs.add_restraint(
1098  IMP.core.DihedralRestraint(self.model, ds, *q2))
A filter which returns true if a container containers 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)
Definition: Residue.h:158
Enforce CHARMM stereochemistry on the given Hierarchy.
A member of a rigid body, it has internal (local) coordinates.
Definition: rigid_bodies.h:616
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:617
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
A class to store an fixed array of same-typed values.
Definition: Array.h:33
Enforce ideal Helix dihedrals and bonds for a selection at resolution 0.
Miscellaneous utilities.
Definition: tools.py:1
Cosine function.
Definition: Cosine.h:22
Creates a restraint between the termini two polypeptides, to simulate the sequence connectivity...
Apply a function to the distance between two particles after transforming the second.
Restrain the dihedral between planes defined by three particles.
Dihedral restraint between four particles.
The type of an atom.
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)
Definition: rigid_bodies.h:731
Return all spatially-proximals 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.
Definition: RestraintSet.h:37
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: tools.py:948
Store a list of ParticleIndexPairs.
A well with harmonic barriers.
Definition: HarmonicWell.h:25
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.
Definition: pmi/utilities.h:26
Add bonds and improper dihedral restraints for the CBs.
Store a list of ParticleIndexes.
A decorator for a particle representing an atom.
Definition: atom/Atom.h:234
def get_particle_pairs
Returns the list of connected particles pairs.
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:91
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.
A decorator for a residue.
Definition: Residue.h:137
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.
Definition: exception.h:49
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.
Definition: PairRestraint.h:29
def get_sorted_segments
Returns sequence-sorted segments array, each containing the first particle the last particle and the ...
Definition: tools.py:1083
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
Applies a PairScore to each Pair in a list.
def get_residue_indexes
Retrieve the residue indexes for the given particle.
Definition: tools.py:569
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 ...
def sublist_iterator
Yield all sublists of length >= lmin and <= lmax.
Definition: tools.py:644
Harmonic function (symmetric about the mean)
Definition: core/Harmonic.h:24
Restraint a set of residues to use ideal helix dihedrals and bonds.