IMP logo
IMP Reference Guide  2.20.2
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(IMP.pmi.restraints.RestraintBase):
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 
46  hiers = IMP.pmi.tools.input_adaptor(objects, resolution)
47  if len(hiers) > 1:
48  raise Exception("ConnectivityRestraint: only pass stuff from "
49  "one Molecule, please")
50  hiers = hiers[0]
51  m = list(hiers)[0].get_model()
52  super(ConnectivityRestraint, self).__init__(m, label=label)
53 
54  self.kappa = 10 # spring constant used for the harmonic restraints
55  SortedSegments = []
56  for h in hiers:
57  try:
58  start = IMP.atom.Hierarchy(h).get_children()[0]
59  except: # noqa: E722
60  start = IMP.atom.Hierarchy(h)
61 
62  try:
63  end = IMP.atom.Hierarchy(h).get_children()[-1]
64  except: # noqa: E722
65  end = IMP.atom.Hierarchy(h)
66 
67  startres = IMP.pmi.tools.get_residue_indexes(start)[0]
68  SortedSegments.append((start, end, startres))
69  SortedSegments = sorted(SortedSegments, key=itemgetter(2))
70 
71  # connect the particles
72  self.particle_pairs = []
73  for x in range(len(SortedSegments) - 1):
74 
75  last = SortedSegments[x][1]
76  first = SortedSegments[x + 1][0]
77 
78  apply_restraint = True
79 
80  # Apply connectivity runless ALL of the following are true:
81  # - first and last both have RigidBodyMember decorators
82  # - first and last are both RigidMembers
83  # - first and last are part of the same RigidBody object
84 
85  # Check for both in a rigid body
90  # Check if the rigid body objects for each particle are
91  # the same object.
92  # if so, skip connectivity restraint
93  if IMP.core.RigidBodyMember(first).get_rigid_body() \
94  == IMP.core.RigidBodyMember(last).get_rigid_body():
95  apply_restraint = False
96 
97  if apply_restraint:
98 
99  nreslast = len(IMP.pmi.tools.get_residue_indexes(last))
100  lastresn = IMP.pmi.tools.get_residue_indexes(last)[-1]
101  nresfirst = len(IMP.pmi.tools.get_residue_indexes(first))
102  firstresn = IMP.pmi.tools.get_residue_indexes(first)[0]
103 
104  residuegap = firstresn - lastresn - 1
105  if disorderedlength and (nreslast / 2 + nresfirst / 2
106  + residuegap) > 20.0:
107  # calculate the distance between the sphere centers
108  # using Kohn PNAS 2004
109  optdist = math.sqrt(5 / 3) * 1.93 * \
110  (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
111  if upperharmonic:
112  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
113  else:
114  hu = IMP.core.Harmonic(optdist, self.kappa)
116  else: # default
117  optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
118  if upperharmonic: # default
119  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
120  else:
121  hu = IMP.core.Harmonic(optdist, self.kappa)
123 
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()))
129 
130  print("Adding sequence connectivity restraint between",
131  pt0.get_name(), " and ", pt1.get_name(), 'of distance',
132  optdist)
133  self.rs.add_restraint(r)
134 
136  """ Returns number of connectivity restraints """
137  return len(self.rs.get_restraints())
138 
140  """ Returns the list of connected particles pairs """
141  return self.particle_pairs
142 
143 
144 class ExcludedVolumeSphere(IMP.pmi.restraints.RestraintBase):
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)
150  """
151  _include_in_rmf = True
152 
153  def __init__(self,
154  included_objects,
155  other_objects=None,
156  resolution=1000,
157  kappa=1.0,
158  label=None):
159  """Constructor.
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
171  """
172 
173  self.kappa = kappa
174  self.cpc = None
175  bipartite = False
176 
177  # gather IMP hierarchies from input objects
178  hierarchies = IMP.pmi.tools.input_adaptor(included_objects,
179  resolution,
180  flatten=True)
181  included_ps = []
182  if other_objects is not None:
183  bipartite = True
184  other_hierarchies = IMP.pmi.tools.input_adaptor(other_objects,
185  resolution,
186  flatten=True)
187  other_ps = []
188 
189  # perform selection
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)
194 
195  included_ps = [h.get_particle() for h in hierarchies]
196  if bipartite:
197  other_ps = [h.get_particle() for h in other_hierarchies]
198 
199  # setup score
200  ssps = IMP.core.SoftSpherePairScore(self.kappa)
201  lsa = IMP.container.ListSingletonContainer(self.model)
202  lsa.add(IMP.get_indexes(included_ps))
203 
204  # setup close pair container
205  if not bipartite:
207  self.cpc = IMP.container.ClosePairContainer(lsa, 0.0, rbcpf, 10.0)
208  evr = IMP.container.PairsRestraint(ssps, self.cpc)
209  else:
210  other_lsa = IMP.container.ListSingletonContainer(self.model)
211  other_lsa.add(IMP.get_indexes(other_ps))
213  lsa,
214  other_lsa,
215  0.0,
216  10.0)
217  evr = IMP.container.PairsRestraint(ssps, self.cpc)
218 
219  self.rs.add_restraint(evr)
220 
221  def add_excluded_particle_pairs(self, excluded_particle_pairs):
222  # add pairs to be filtered when calculating the score
223  inverted = [(p1, p0) for p0, p1 in excluded_particle_pairs]
224  lpc = IMP.container.ListPairContainer(self.model)
225  lpc.add(IMP.get_indexes(excluded_particle_pairs))
226  lpc.add(IMP.get_indexes(inverted))
228  self.cpc.add_pair_filter(icpf)
229 
230 
231 class HelixRestraint(IMP.pmi.restraints.RestraintBase):
232  """Enforce ideal Helix dihedrals and bonds for a selection
233  at resolution 0"""
234  def __init__(self,
235  hierarchy,
236  selection_tuple,
237  weight=1.0,
238  label=None):
239  """Constructor
240  @param hierarchy the root node
241  @param selection_tuple (start, stop, molname, copynum=0)
242  @param weight
243  """
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]
249  copy_index = 0
250  if len(selection_tuple) > 3:
251  copy_index = selection_tuple[3]
252 
253  sel = IMP.atom.Selection(
254  hierarchy, molecule=mol, copy_index=copy_index,
255  residue_indexes=range(start, stop+1))
256  ps = sel.get_selected_particles(with_representation=False)
257  res = [IMP.atom.Residue(p) for p in ps]
258  self.r = IMP.atom.HelixRestraint(res)
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()))
263 
264  def get_number_of_bonds(self):
265  return self.r.get_number_of_bonds()
266 
267  def get_number_of_dihedrals(self):
268  return self.r.get_number_of_dihedrals()
269 
270 
271 class ResidueBondRestraint(IMP.pmi.restraints.RestraintBase):
272  """ Add bond restraint between pair of consecutive
273  residues/beads to enforce the stereochemistry.
274  """
275  def __init__(self, objects, distance=3.78, strength=10.0, jitter=None):
276  """Constructor
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
282  """
283 
284  particles = IMP.pmi.tools.input_adaptor(objects, 1, flatten=True)
285  m = particles[0].get_model()
286  super(ResidueBondRestraint, self).__init__(m)
287 
288  self.pairslist = []
289 
290  if not jitter:
291  ts = IMP.core.Harmonic(distance, strength)
292  else:
294  (distance - jitter, distance + jitter), strength)
295 
296  for ps in IMP.pmi.tools.sublist_iterator(particles, 2, 2):
297  pair = []
298  if len(ps) != 2:
299  raise ValueError("wrong length of pair")
300  for p in ps:
302  raise TypeError("%s is not a residue" % p)
303  else:
304  pair.append(p)
305  print("ResidueBondRestraint: adding a restraint between %s %s"
306  % (pair[0].get_name(), pair[1].get_name()))
307  self.rs.add_restraint(
308  IMP.core.DistanceRestraint(self.model, ts, pair[0], pair[1]))
309  self.pairslist.append(IMP.ParticlePair(pair[0], pair[1]))
310  self.pairslist.append(IMP.ParticlePair(pair[1], pair[0]))
311 
312  def get_excluded_pairs(self):
313  return self.pairslist
314 
315 
316 class ResidueAngleRestraint(IMP.pmi.restraints.RestraintBase):
317  """Add angular restraint between triplets of consecutive
318  residues/beads to enforce the stereochemistry.
319  """
320  def __init__(self, objects, anglemin=100.0, anglemax=140.0, strength=10.0):
321 
322  particles = IMP.pmi.tools.input_adaptor(objects, 1, flatten=True)
323  m = particles[0].get_model()
324  super(ResidueAngleRestraint, self).__init__(m)
325 
326  self.pairslist = []
327 
329  (math.pi * anglemin / 180.0,
330  math.pi * anglemax / 180.0),
331  strength)
332 
333  for ps in IMP.pmi.tools.sublist_iterator(particles, 3, 3):
334  triplet = []
335  if len(ps) != 3:
336  raise ValueError("wrong length of triplet")
337  for p in ps:
339  raise TypeError("%s is not a residue" % p)
340  else:
341  triplet.append(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(
346  IMP.core.AngleRestraint(triplet[0].get_model(), ts,
347  triplet[0],
348  triplet[1],
349  triplet[2]))
350  self.pairslist.append(IMP.ParticlePair(triplet[0], triplet[2]))
351  self.pairslist.append(IMP.ParticlePair(triplet[2], triplet[0]))
352 
353  def get_excluded_pairs(self):
354  return self.pairslist
355 
356 
357 class ResidueDihedralRestraint(IMP.pmi.restraints.RestraintBase):
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.
363  """
364  def __init__(self, objects, stringsequence=None, strength=10.0):
365 
366  particles = IMP.pmi.tools.input_adaptor(objects, 1, flatten=True)
367  m = particles[0].get_model()
368  super(ResidueDihedralRestraint, self).__init__(m)
369 
370  self.pairslist = []
371 
372  if stringsequence is None:
373  stringsequence = "T" * (len(particles) - 3)
374 
375  for n, ps in enumerate(
376  IMP.pmi.tools.sublist_iterator(particles, 4, 4)):
377  quadruplet = []
378  if len(ps) != 4:
379  raise ValueError("wrong length of quadruplet")
380  for p in ps:
382  raise TypeError("%s is not a residue" % p)
383  else:
384  quadruplet.append(p)
385  dihedraltype = stringsequence[n]
386  if dihedraltype == "C":
387  anglemin = -20.0
388  anglemax = 20.0
390  (math.pi * anglemin / 180.0,
391  math.pi * anglemax / 180.0),
392  strength)
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),
403  strength)
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(
410  IMP.core.DihedralRestraint(self.model, ts,
411  quadruplet[0],
412  quadruplet[1],
413  quadruplet[2],
414  quadruplet[3]))
415  self.pairslist.append(
416  IMP.ParticlePair(quadruplet[0], quadruplet[3]))
417  self.pairslist.append(
418  IMP.ParticlePair(quadruplet[3], quadruplet[0]))
419 
420 
422  """Add harmonic restraints between all pairs
423  """
424  def __init__(self, hierarchy, selection_tuples=None, resolution=1,
425  strength=0.01, dist_cutoff=10.0, ca_only=True):
426  """Constructor
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.
434  """
435 
436  particles = []
437  self.m = hierarchy.get_model()
438  for st in selection_tuples:
439  copy_index = 0
440  if len(st) > 3:
441  copy_index = st[3]
442  if not ca_only:
443  sel = IMP.atom.Selection(
444  hierarchy, molecule=st[2],
445  residue_indexes=range(st[0], st[1]+1),
446  copy_index=copy_index)
447  else:
448  sel = IMP.atom.Selection(
449  hierarchy, molecule=st[2],
450  residue_indexes=range(st[0], st[1]+1),
451  copy_index=copy_index,
452  atom_type=IMP.atom.AtomType("CA"))
453  particles += sel.get_selected_particles()
454 
455  self.weight = 1
456  self.label = "None"
457  self.pairslist = []
458 
459  # create score
461  particles, dist_cutoff, strength)
462  for r in self.rs.get_restraints():
463  a1, a2 = r.get_inputs()
464  self.pairslist.append(IMP.ParticlePair(a1, a2))
465  self.pairslist.append(IMP.ParticlePair(a2, a1))
466  print('ElasticNetwork: created', self.rs.get_number_of_restraints(),
467  'restraints')
468 
469  def set_label(self, label):
470  self.label = label
471  self.rs.set_name(label)
472  for r in self.rs.get_restraints():
473  r.set_name(label)
474 
475  def add_to_model(self):
476  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
477 
478  def get_restraint(self):
479  return self.rs
480 
481  def set_weight(self, weight):
482  self.weight = weight
483  self.rs.set_weight(weight)
484 
485  def get_excluded_pairs(self):
486  return self.pairslist
487 
488  def get_output(self):
489  output = {}
490  score = self.weight * self.rs.unprotected_evaluate(None)
491  output["_TotalScore"] = str(score)
492  output["ElasticNetworkRestraint_" + self.label] = str(score)
493  return output
494 
495 
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!
511  """
512 
513  kB = (1.381 * 6.02214) / 4184.0
514 
515  self.mdl = root.get_model()
516  self.bonds_rs = IMP.RestraintSet(self.mdl, 1.0 / (kB * ff_temp),
517  'BONDED')
518  self.nonbonded_rs = IMP.RestraintSet(self.mdl, 1.0 / (kB * ff_temp),
519  'NONBONDED')
520  self.weight = 1.0
521  self.label = ""
522 
523  # setup topology and bonds etc
524  if enable_bonded:
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)
534  topology,
535  limit_to_ps)
536  self.ps = limit_to_ps
537  else:
538  r = IMP.atom.CHARMMStereochemistryRestraint(root, topology)
539  self.ps = IMP.core.get_leaves(root)
540  print('init bonds score', r.unprotected_evaluate(None))
541  self.bonds_rs.add_restraint(r)
542  ff.add_radii(root)
543  ff.add_well_depths(root)
544 
545  atoms = IMP.atom.get_by_type(root, IMP.atom.ATOM_TYPE)
546  # non-bonded forces
547  if enable_nonbonded:
548  if (zone_ps is not None) and zone_nonbonded:
549  print('stereochemistry: zone_nonbonded is True')
550  # atoms list should only include backbone plus zone_ps!
551  backbone_types = ['C', 'N', 'CB', 'O']
552  sel = IMP.atom.Selection(
553  root, atom_types=[IMP.atom.AtomType(n)
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)
559 
561  IMP.container.ListSingletonContainer(backbone_atoms),
563  4.0)
564  else:
565  cont = IMP.container.ListSingletonContainer(self.mdl, atoms)
566  self.nbl = IMP.container.ClosePairContainer(cont, 4.0)
567  if enable_bonded:
568  self.nbl.add_pair_filter(r.get_full_pair_filter())
569  pairscore = IMP.isd.RepulsiveDistancePairScore(0, 1)
570  pr = IMP.container.PairsRestraint(pairscore, self.nbl)
571  self.nonbonded_rs.add_restraint(pr)
572  print('CHARMM is set up')
573 
574  def set_label(self, label):
575  self.label = label
576  self.rs.set_name(label)
577  for r in self.rs.get_restraints():
578  r.set_name(label)
579 
580  def add_to_model(self):
581  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.bonds_rs)
582  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.nonbonded_rs)
583 
584  def get_restraint(self):
585  return self.rs
586 
587  def get_close_pair_container(self):
588  return self.nbl
589 
590  def set_weight(self, weight):
591  self.weight = weight
592  self.rs.set_weight(weight)
593 
594  def get_output(self):
595  output = {}
596  bonds_score = self.weight * self.bonds_rs.unprotected_evaluate(None)
597  nonbonded_score = \
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)
603  return output
604 
605 
606 class PseudoAtomicRestraint(object):
607  """Add bonds and improper dihedral restraints for the CBs
608  """
609  def __init__(
610  self, rnums, representation, selection_tuple, strength=10.0, kappa=1.0,
611  jitter_angle=0.0, jitter_improper=0.0):
612  '''
613  need to add:
614  ca-ca bond
615  ca-cb is a constraint, no restraint needed
616  ca-ca-ca
617  cb-ca-ca-cb
618  '''
619 
620  self.m = representation.prot.get_model()
621  self.rs = IMP.RestraintSet(self.m, "PseudoAtomic")
622  self.rset_angles = IMP.RestraintSet(self.m, "PseudoAtomic_Angles")
623  self.rset_bonds = IMP.RestraintSet(self.m, "PseudoAtomic_Bonds")
624  self.weight = 1
625  self.label = "None"
626  self.pairslist = []
627 
628  for rnum in rnums:
629  ca, cb = self.get_ca_cb(
630  IMP.pmi.tools.select_by_tuple(representation,
631  (rnum, rnum, 'chainA'),
632  resolution=0))
633  if cb is not None:
634  nter = False
635  cter = False
636  ca_prev, cb_prev = self.get_ca_cb(
637  IMP.pmi.tools.select_by_tuple(representation,
638  (rnum - 1, rnum - 1,
639  'chainA'), resolution=0))
640  ca_next, cb_next = self.get_ca_cb(
641  IMP.pmi.tools.select_by_tuple(representation,
642  (rnum + 1, rnum + 1,
643  'chainA'), resolution=0))
644  if ca_prev is None:
645  nter = True
646  else:
647  if ca_next is None:
648  cter = True
649  else:
650  if (nter and cter):
651  continue
652 
653  # adding a bond restraint between CA and CB
654  # h=IMP.core.Harmonic(6.0,kappa)
655  # dps=IMP.core.DistancePairScore(h)
656  # pr=IMP.core.PairRestraint(dps,IMP.ParticlePair(ca,cb))
657  # self.pairslist.append((ca,cb))
658  # self.rset_bonds.add_restraint(pr)
659 
660  # creating improper dihedral restraint
661  # hus=IMP.core.Harmonic(2.09,kappa)
663  2.09 +
664  jitter_angle /
665  0.5,
666  kappa)
668  2.09 -
669  jitter_angle /
670  0.5,
671  kappa)
672  if not nter:
673  # ar13=IMP.core.AngleRestraint(hus,ca_prev,ca,cb)
674  # self.rset_angles.add_restraint(ar13)
675  ar13u = IMP.core.AngleRestraint(hupp, ca_prev, ca, cb)
676  ar13l = IMP.core.AngleRestraint(hlow, ca_prev, ca, cb)
677  self.rset_angles.add_restraint(ar13u)
678  self.rset_angles.add_restraint(ar13l)
679  if not cter:
680  # ar23=IMP.core.AngleRestraint(hus,ca_next,ca,cb)
681  # self.rset_angles.add_restraint(ar23)
682  ar23u = IMP.core.AngleRestraint(hupp, ca_next, ca, cb)
683  ar23l = IMP.core.AngleRestraint(hlow, ca_next, ca, cb)
684  self.rset_angles.add_restraint(ar23u)
685  self.rset_angles.add_restraint(ar23l)
686  if not nter and not cter:
687  # hus2=IMP.core.Harmonic(0,kappa)
688  # idr=IMP.core.DihedralRestraint(hus2,ca,ca_prev,ca_next,cb)
689  # self.rset_angles.add_restraint(idr)
690 
691  hus2upp = IMP.core.HarmonicUpperBound(
692  jitter_improper,
693  kappa)
694  hus2low = IMP.core.HarmonicLowerBound(
695  -
696  jitter_improper,
697  kappa)
699  hus2upp,
700  ca,
701  ca_prev,
702  ca_next,
703  cb)
705  hus2low,
706  ca,
707  ca_prev,
708  ca_next,
709  cb)
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)
714 
715  def get_ca_cb(self, atoms):
716  ca = None
717  cb = None
718  for a in atoms:
719  if IMP.atom.Atom(a).get_atom_type() == IMP.atom.AtomType("CA"):
720  ca = a.get_particle()
721  elif IMP.atom.Atom(a).get_atom_type() == IMP.atom.AtomType("CB"):
722  cb = a.get_particle()
723  return ca, cb
724 
725  def set_label(self, label):
726  self.label = label
727  self.rs.set_name(label)
728  for r in self.rs.get_restraints():
729  r.set_name(label)
730 
731  def add_to_model(self):
732  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
733 
734  def get_restraint(self):
735  return self.rs
736 
737  def set_weight(self, weight):
738  self.weight = weight
739  self.rs.set_weight(weight)
740 
741  def get_excluded_pairs(self):
742  return self.pairslist
743 
744  def get_output(self):
745  output = {}
746  score = self.weight * self.rs.unprotected_evaluate(None)
747  output["_TotalScore"] = str(score)
748  output["PseudoAtomicRestraint_" + self.label] = str(score)
749  return output
750 
751 
752 class SymmetryRestraint(object):
753  """Create harmonic restraints between the reference and (transformed)
754  clones.
755 
756  @note Wraps IMP::core::TransformedDistancePairScore with an
757  IMP::core::Harmonic
758  """
759  def __init__(self, references, clones_list, transforms,
760  label='', strength=10.0, ca_only=False):
761  """Constructor
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
767  selection
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
771  """
772 
773  refs = IMP.pmi.tools.input_adaptor(references, flatten=True)
774  self.mdl = refs[0].get_model()
775  self.rs = IMP.RestraintSet(self.mdl, "Symmetry")
776  self.weight = 1
777  self.label = label
778  if len(clones_list) != len(transforms):
779  raise Exception(
780  'Error: There should be as many clones as transforms')
781 
782  harmonic = IMP.core.Harmonic(0., strength)
783  for tmp_clones, trans in zip(clones_list, transforms):
784  clones = IMP.pmi.tools.input_adaptor(tmp_clones, flatten=True)
785  if len(clones) != len(refs):
786  raise Exception("Error: len(references)!=len(clones)")
787  pair_score = IMP.core.TransformedDistancePairScore(harmonic, trans)
788  for p0, p1 in zip(refs, clones):
789  if not ca_only or (
790  IMP.atom.Atom(p0).get_atom_type()
791  == IMP.atom.AtomType("CA") and
792  IMP.atom.Atom(p1).get_atom_type()
793  == IMP.atom.AtomType("CA")):
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')
800 
801  def set_label(self, label):
802  self.label = label
803  self.rs.set_name(label)
804  for r in self.rs.get_restraints():
805  r.set_name(label)
806 
807  def add_to_model(self):
808  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.rs)
809 
810  def get_restraint(self):
811  return self.rs
812 
813  def set_weight(self, weight):
814  self.weight = weight
815  self.rs.set_weight(weight)
816 
817  def get_excluded_pairs(self):
818  return self.pairslist
819 
820  def get_output(self):
821  output = {}
822  score = self.weight * self.rs.unprotected_evaluate(None)
823  output["SymmetryRestraint_" + self.label] = str(score)
824  output["_TotalScore"] = str(score)
825  return output
826 
827 
828 class FusionRestraint(object):
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"):
833  """
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
844  restraint.
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
851  output/stat file
852  """
853  self.label = label
854  self.weight = 1.0
855  ssn = IMP.pmi.tools.get_sorted_segments(nterminal)
856  ssc = IMP.pmi.tools.get_sorted_segments(cterminal)
857  nter_lastres = ssn[-1][1]
858  cter_firstres = ssc[0][0]
859  self.m = nter_lastres.get_model()
860 
861  self.kappa = 10 # spring constant used for the harmonic restraints
862 
863  optdist = (3.6) * scale
864  if upperharmonic: # default
865  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
866  else:
867  hu = IMP.core.Harmonic(optdist, self.kappa)
869 
870  pt0 = nter_lastres.get_particle()
871  pt1 = cter_firstres.get_particle()
872  r = IMP.core.PairRestraint(self.m, dps,
873  (pt0.get_index(), pt1.get_index()))
874  self.rs = IMP.RestraintSet(self.m, "fusion_restraint")
875  print("Adding fusion connectivity restraint between", pt0.get_name(),
876  " and ", pt1.get_name(), 'of distance', optdist)
877  self.rs.add_restraint(r)
878 
879  def set_label(self, label):
880  self.label = label
881 
882  def get_weight(self):
883  return self.weight
884 
885  def add_to_model(self):
886  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
887 
888  def get_restraint(self):
889  return self.rs
890 
891  def set_weight(self, weight):
892  self.weight = weight
893  self.rs.set_weight(weight)
894 
895  def get_output(self):
896  output = {}
897  score = self.weight * self.rs.unprotected_evaluate(None)
898  output["_TotalScore"] = str(score)
899  output["FusionRestraint_" + self.label] = str(score)
900  return output
901 
902  def evaluate(self):
903  return self.weight * self.rs.unprotected_evaluate(None)
904 
905 
906 class PlaneDihedralRestraint(IMP.pmi.restraints.RestraintBase):
907 
908  """Restrain the dihedral between planes defined by three particles.
909 
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.
912  """
913 
914  def __init__(self, particle_triplets, angle=0., k=1., label=None,
915  weight=1.):
916  """Constructor
917  @param particle_triplets List of lists of 3 particles. Each triplet
918  defines a plane. Dihedrals of adjacent planes
919  in list are scored.
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
926  """
927  model = particle_triplets[0][0].get_model()
928  super(PlaneDihedralRestraint, self).__init__(model, label=label,
929  weight=weight)
930 
931  angle = math.pi * angle / 180.
932  ds = IMP.core.Cosine(.5 * k, 1, -angle)
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(
938  IMP.core.DihedralRestraint(self.model, ds, *q1))
939  self.rs.add_restraint(
940  IMP.core.DihedralRestraint(self.model, ds, *q2))
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)
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:634
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:635
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
A class to store a fixed array of same-typed values.
Definition: Array.h:35
Enforce ideal Helix dihedrals and bonds for a selection at resolution 0.
Miscellaneous utilities.
Definition: tools.py:1
Cosine function.
Definition: Cosine.h:25
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:749
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.
Definition: RestraintSet.h:41
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: tools.py:888
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:238
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:92
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.
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:48
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:31
def get_sorted_segments
Returns sequence-sorted segments array, each containing the first particle the last particle and the ...
Definition: tools.py:1023
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:70
Applies a PairScore to each Pair in a list.
def get_residue_indexes
Retrieve the residue indexes for the given particle.
Definition: tools.py:512
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:587
Harmonic function (symmetric about the mean)
Definition: core/Harmonic.h:27
Restraint a set of residues to use ideal helix dihedrals and bonds.