IMP logo
IMP Reference Guide  2.9.0
The Integrative Modeling Platform
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.algebra
9 import IMP.atom
10 import IMP.container
11 import IMP.isd
12 import itertools
13 import IMP.pmi.tools
15 from operator import itemgetter
16 from math import pi,log,sqrt
17 import sys
18 
19 
20 
21 class ConnectivityRestraint(object):
22  """ This class creates a restraint between consecutive TempResidue objects OR an entire
23  PMI MOlecule object. """
24 
25  def __init__(self,
26  objects,
27  scale=1.0,
28  disorderedlength=False,
29  upperharmonic=True,
30  resolution=1,
31  label="None"):
32  """
33  @param objects - a list of hierarchies, PMI TempResidues OR a single Molecule
34  @param scale Scale the maximal distance between the beads by this factor when disorderedlength is False.
35  The maximal distance is calculated as ((float(residuegap) + 1.0) * 3.6) * scale.
36  @param disorderedlength - This flag uses either disordered length
37  calculated for random coil peptides (True) or zero
38  surface-to-surface distance between beads (False)
39  as optimal distance for the sequence connectivity
40  restraint.
41  @param upperharmonic - This flag uses either harmonic (False)
42  or upperharmonic (True) in the intra-pair
43  connectivity restraint.
44  @param resolution - The resolution to connect things at - only used if you pass PMI objects
45  @param label - A string to identify this restraint in the output/stat file
46  """
47  self.label = label
48  self.weight = 1.0
49 
50  hiers = IMP.pmi.tools.input_adaptor(objects,resolution)
51  if len(hiers)>1:
52  raise Exception("ConnectivityRestraint: only pass stuff from one Molecule, please")
53  hiers = hiers[0]
54 
55  self.kappa = 10 # spring constant used for the harmonic restraints
56  self.m = list(hiers)[0].get_model()
57  SortedSegments = []
58  self.rs = IMP.RestraintSet(self.m, "connectivity_restraint")
59  for h in hiers:
60  try:
61  start = IMP.atom.Hierarchy(h).get_children()[0]
62  except:
63  start = IMP.atom.Hierarchy(h)
64 
65  try:
66  end = IMP.atom.Hierarchy(h).get_children()[-1]
67  except:
68  end = IMP.atom.Hierarchy(h)
69 
70  startres = IMP.pmi.tools.get_residue_indexes(start)[0]
71  endres = IMP.pmi.tools.get_residue_indexes(end)[-1]
72  SortedSegments.append((start, end, startres))
73  SortedSegments = sorted(SortedSegments, key=itemgetter(2))
74 
75  # connect the particles
76  self.particle_pairs=[]
77  for x in range(len(SortedSegments) - 1):
78 
79  last = SortedSegments[x][1]
80  first = SortedSegments[x + 1][0]
81 
82  apply_restraint = True
83 
84  # Apply connectivity runless ALL of the following are true:
85  # - first and last both have RigidBodyMember decorators
86  # - first and last are both RigidMembers
87  # - first and last are part of the same RigidBody object
88 
89  # Check for both in a rigid body
92 
93  #Check if the rigid body objects for each particle are the same object.
94  # if so, skip connectivity restraint
95  if IMP.core.RigidBodyMember(first).get_rigid_body() == IMP.core.RigidBodyMember(last).get_rigid_body():
96  apply_restraint = False
97 
98  if apply_restraint:
99 
100  nreslast = len(IMP.pmi.tools.get_residue_indexes(last))
101  lastresn = IMP.pmi.tools.get_residue_indexes(last)[-1]
102  nresfirst = len(IMP.pmi.tools.get_residue_indexes(first))
103  firstresn = IMP.pmi.tools.get_residue_indexes(first)[0]
104 
105  residuegap = firstresn - lastresn - 1
106  if disorderedlength and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
107  # calculate the distance between the sphere centers using Kohn
108  # PNAS 2004
109  optdist = sqrt(5 / 3) * 1.93 * \
110  (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
111  # optdist2=sqrt(5/3)*1.93*((nreslast)**0.6+(nresfirst)**0.6)/2
112  if upperharmonic:
113  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
114  else:
115  hu = IMP.core.Harmonic(optdist, self.kappa)
117  else: # default
118  optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
119  if upperharmonic: # default
120  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
121  else:
122  hu = IMP.core.Harmonic(optdist, self.kappa)
124 
125  pt0 = last.get_particle()
126  pt1 = first.get_particle()
127  self.particle_pairs.append((pt0,pt1))
128  r = IMP.core.PairRestraint(self.m, dps, (pt0.get_index(), pt1.get_index()))
129 
130  print("Adding sequence connectivity restraint between", pt0.get_name(), " and ", pt1.get_name(), 'of distance', optdist)
131  self.rs.add_restraint(r)
132 
133  def set_label(self, label):
134  self.label = label
135 
136  def get_weight(self):
137  return self.weight
138 
139  def add_to_model(self):
140  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
141 
142  def get_restraint(self):
143  return self.rs
144 
145  def set_weight(self, weight):
146  self.weight = weight
147  self.rs.set_weight(weight)
148 
149  def get_output(self):
150  self.m.update()
151  output = {}
152  score = self.weight * self.rs.unprotected_evaluate(None)
153  output["_TotalScore"] = str(score)
154  output["ConnectivityRestraint_" + self.label] = str(score)
155  return output
156 
158  """ Returns number of connectivity restraints """
159  return len(self.rs.get_restraints())
160 
162  """ Returns the list of connected particles pairs """
163  return self.particle_pairs
164 
165  def evaluate(self):
166  return self.weight * self.rs.unprotected_evaluate(None)
167 
168 class ExcludedVolumeSphere(object):
169  """A class to create an excluded volume restraint for a set of particles at a given resolution.
170  Can be initialized as a bipartite restraint between two sets of particles.
171  # Potential additional function: Variable resolution for each PMI object. Perhaps passing selection_tuples
172  with (PMI_object, resolution)
173  """
174 
175  def __init__(self,
176  representation=None,
177  included_objects=None,
178  other_objects=None,
179  resolution=1000,
180  kappa=1.0):
181  """Constructor.
182  @param representation DEPRECATED - just pass objects
183  @param included_objects Can be one of the following inputs:
184  IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set of them
185  @param other_objects Initializes a bipartite restraint between included_objects and other_objects
186  Same format as included_objects
187  @param resolution The resolution particles at which to impose the restraint.
188  By default, the coarsest particles will be chosen.
189  If a number is chosen, for each particle, the closest
190  resolution will be used (see IMP.atom.Selection).
191  @param kappa Restraint strength
192  """
193 
194  self.weight = 1.0
195  self.kappa = kappa
196  self.label = "None"
197  self.cpc = None
198  bipartite = False
199 
200  # gather IMP hierarchies from input objects
201  hierarchies = IMP.pmi.tools.input_adaptor(included_objects,
202  resolution,
203  flatten=True)
204  included_ps = []
205  if other_objects is not None:
206  bipartite = True
207  other_hierarchies = IMP.pmi.tools.input_adaptor(other_objects,
208  resolution,
209  flatten=True)
210  other_ps = []
211 
212 
213  # perform selection
214  if representation is None:
215  if hierarchies is None:
216  raise Exception("Must at least pass included objects")
217  self.mdl = hierarchies[0].get_model()
218  included_ps = [h.get_particle() for h in hierarchies]
219  if bipartite:
220  other_ps = [h.get_particle() for h in other_hierarchies]
221  elif isinstance(representation, IMP.pmi.representation.Representation):
222  self.mdl = representation.m
223  included_ps = IMP.pmi.tools.select(
224  representation,
225  resolution=resolution,
226  hierarchies=hierarchies)
227  if bipartite:
228  other_ps = IMP.pmi.tools.select(
229  representation,
230  resolution=resolution,
231  hierarchies=other_hierarchies)
232  else:
233  raise Exception("Must pass Representation or included_objects")
234 
235  # setup score
236  self.rs = IMP.RestraintSet(self.mdl, 'excluded_volume')
237  ssps = IMP.core.SoftSpherePairScore(self.kappa)
239  lsa.add(IMP.get_indexes(included_ps))
240 
241  # setup close pair container
242  if not bipartite:
244  self.cpc = IMP.container.ClosePairContainer(lsa, 0.0, rbcpf, 10.0)
245  evr = IMP.container.PairsRestraint(ssps, self.cpc)
246  else:
247  other_lsa = IMP.container.ListSingletonContainer(self.mdl)
248  other_lsa.add(IMP.get_indexes(other_ps))
250  lsa,
251  other_lsa,
252  0.0,
253  10.0)
254  evr = IMP.container.PairsRestraint(ssps, self.cpc)
255 
256  self.rs.add_restraint(evr)
257 
258  def add_excluded_particle_pairs(self, excluded_particle_pairs):
259  # add pairs to be filtered when calculating the score
260  inverted=[(p1,p0)for p0,p1 in excluded_particle_pairs]
261  lpc = IMP.container.ListPairContainer(self.mdl)
262  lpc.add(IMP.get_indexes(excluded_particle_pairs))
263  lpc.add(IMP.get_indexes(inverted))
265  self.cpc.add_pair_filter(icpf)
266 
267  def set_label(self, label):
268  self.label = label
269 
270  def add_to_model(self):
271  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.rs)
272 
273  def get_restraint(self):
274  return self.rs
275 
276  def set_weight(self, weight):
277  self.weight = weight
278  self.rs.set_weight(weight)
279 
280  def get_output(self):
281  self.mdl.update()
282  output = {}
283  score = self.weight * self.rs.unprotected_evaluate(None)
284  output["_TotalScore"] = str(score)
285  output["ExcludedVolumeSphere_" + self.label] = str(score)
286  return output
287 
288  def evaluate(self):
289  return self.weight * self.rs.unprotected_evaluate(None)
290 
291 class HelixRestraint(object):
292  """Enforce ideal Helix dihedrals and bonds for a selection at resolution 0"""
293  def __init__(self,
294  hierarchy,
295  selection_tuple,
296  weight=1.0,
297  label=''):
298  """Constructor
299  @param hierarchy the root node
300  @param selection_tuple (start, stop, molname, copynum=0)
301  @param weight
302  """
303  self.mdl = hierarchy.get_model()
304  self.rs = IMP.RestraintSet(self.mdl)
305  self.weight = weight
306  self.label = label
307  start = selection_tuple[0]
308  stop = selection_tuple[1]
309  mol = selection_tuple[2]
310  copy_index = 0
311  if len(selection_tuple)>3:
312  copy_index = selection_tuple[3]
313 
314  sel = IMP.atom.Selection(hierarchy,molecule=mol,copy_index=copy_index,
315  residue_indexes=range(start,stop+1))
316  ps = sel.get_selected_particles(with_representation=False)
317  res = [IMP.atom.Residue(p) for p in ps]
318  self.rs = IMP.RestraintSet(self.mdl,self.weight)
319  self.r = IMP.atom.HelixRestraint(res)
320  self.rs.add_restraint(self.r)
321  print('Created helix %s.%i.%i-%i with %i dihedrals and %i bonds'%(
322  mol,copy_index,start,stop,
323  self.get_number_of_bonds(),self.get_number_of_dihedrals()))
324  def set_label(self, label):
325  self.label = label
326 
327  def get_weight(self):
328  return self.weight
329 
330  def add_to_model(self):
331  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.rs)
332 
333  def get_restraint(self):
334  return self.rs
335 
336  def set_weight(self, weight):
337  self.weight = weight
338  self.rs.set_weight(weight)
339 
340  def get_number_of_bonds(self):
341  return self.r.get_number_of_bonds()
342 
343  def get_number_of_dihedrals(self):
344  return self.r.get_number_of_dihedrals()
345 
346  def evaluate(self):
347  return self.weight * self.rs.unprotected_evaluate(None)
348 
349  def get_output(self):
350  self.mdl.update()
351  output = {}
352  score = self.evaluate()
353  output["_TotalScore"] = str(score)
354  output["HelixRestraint_" + self.label] = str(score)
355  return output
356 
357 class ResidueBondRestraint(object):
358  """ Add bond restraint between pair of consecutive
359  residues/beads to enforce the stereochemistry.
360  """
361  def __init__(self,
362  representation=None,
363  selection_tuple=None,
364  objects=None,
365  distance=3.78,
366  strength=10.0,
367  jitter=None):
368  """Constructor
369  @param representation (PMI1)
370  @param selection_tuple Requested selection (PMI1)
371  @param objects (PMI2)
372  @param distance Resting distance for restraint
373  @param strength Bond constant
374  @param jitter Defines the +- added to the optimal distance in the harmonic well restraint
375  used to increase the tolerance
376  """
377 
378  if representation is not None and selection_tuple is not None:
379  self.m = representation.prot.get_model()
380  particles = IMP.pmi.tools.select_by_tuple(
381  representation,
382  selection_tuple,
383  resolution=1)
384 
385  elif objects is not None:
386  particles = IMP.pmi.tools.input_adaptor(objects,1,flatten=True)
387  self.m = particles[0].get_model()
388 
389  self.rs = IMP.RestraintSet(self.m, "Bonds")
390  self.weight = 1
391  self.label = "None"
392  self.pairslist = []
393 
394 
395 
396  if not jitter:
397  ts = IMP.core.Harmonic(distance, strength)
398  else:
400  (distance - jitter, distance + jitter), strength)
401 
402  for ps in IMP.pmi.tools.sublist_iterator(particles, 2, 2):
403  pair = []
404  if len(ps) != 2:
405  raise ValueError("wrong length of pair")
406  for p in ps:
408  raise TypeError("not a residue")
409  else:
410  pair.append(p)
411  print("ResidueBondRestraint: adding a restraint between %s %s" % (pair[0].get_name(), pair[1].get_name()))
412  self.rs.add_restraint(
413  IMP.core.DistanceRestraint(self.m, ts, pair[0], pair[1]))
414  self.pairslist.append(IMP.ParticlePair(pair[0], pair[1]))
415  self.pairslist.append(IMP.ParticlePair(pair[1], pair[0]))
416 
417  def set_label(self, label):
418  self.label = label
419  self.rs.set_name(label)
420  for r in self.rs.get_restraints():
421  r.set_name(label)
422 
423  def add_to_model(self):
424  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
425 
426  def get_restraint(self):
427  return self.rs
428 
429  def set_weight(self, weight):
430  self.weight = weight
431  self.rs.set_weight(weight)
432 
433  def get_excluded_pairs(self):
434  return self.pairslist
435 
436  def get_output(self):
437  self.m.update()
438  output = {}
439  score = self.weight * self.rs.unprotected_evaluate(None)
440  output["_TotalScore"] = str(score)
441  output["ResidueBondRestraint_" + self.label] = str(score)
442  return output
443 
444 
445 class ResidueAngleRestraint(object):
446  """Add angular restraint between triplets of consecutive
447  residues/beads to enforce the stereochemistry.
448  """
449  def __init__(self,
450  representation=None,
451  selection_tuple=None,
452  objects=None,
453  anglemin=100.0,
454  anglemax=140.0,
455  strength=10.0):
456 
457  if representation is not None and selection_tuple is not None:
458  self.m = representation.prot.get_model()
459  particles = IMP.pmi.tools.select_by_tuple(
460  representation,
461  selection_tuple,
462  resolution=1)
463 
464  elif objects is not None:
465  particles = IMP.pmi.tools.input_adaptor(objects,1,flatten=True)
466  self.m = particles[0].get_model()
467 
468  self.rs = IMP.RestraintSet(self.m, "Angles")
469  self.weight = 1
470  self.label = "None"
471  self.pairslist = []
472 
474  (pi * anglemin / 180.0,
475  pi * anglemax / 180.0),
476  strength)
477 
478  for ps in IMP.pmi.tools.sublist_iterator(particles, 3, 3):
479  triplet = []
480  if len(ps) != 3:
481  raise ValueError("wrong length of triplet")
482  for p in ps:
484  raise TypeError("not a residue")
485  else:
486  triplet.append(p)
487  print("ResidueAngleRestraint: adding a restraint between %s %s %s" % (triplet[0].get_name(), triplet[1].get_name(), triplet[2].get_name()))
488  self.rs.add_restraint(
489  IMP.core.AngleRestraint(triplet[0].get_model(), ts,
490  triplet[0],
491  triplet[1],
492  triplet[2]))
493  self.pairslist.append(IMP.ParticlePair(triplet[0], triplet[2]))
494  self.pairslist.append(IMP.ParticlePair(triplet[2], triplet[0]))
495 
496  def set_label(self, label):
497  self.label = label
498  self.rs.set_name(label)
499  for r in self.rs.get_restraints():
500  r.set_name(label)
501 
502  def add_to_model(self):
503  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
504 
505  def get_restraint(self):
506  return self.rs
507 
508  def set_weight(self, weight):
509  self.weight = weight
510  self.rs.set_weight(weight)
511 
512  def get_excluded_pairs(self):
513  return self.pairslist
514 
515  def get_output(self):
516  self.m.update()
517  output = {}
518  score = self.weight * self.rs.unprotected_evaluate(None)
519  output["_TotalScore"] = str(score)
520  output["ResidueAngleRestraint_" + self.label] = str(score)
521  return output
522 
523 
525  """Add dihedral restraints between quadruplet of consecutive
526  residues/beads to enforce the stereochemistry.
527  Give as input a string of "C" and "T", meaning cys (0+-40) or trans (180+-40)
528  dihedral. The length of the string must be \#residue-3.
529  Without the string, the dihedral will be assumed trans.
530  """
531  def __init__(
532  self,
533  representation=None,
534  selection_tuple=None,
535  objects=None,
536  stringsequence=None,
537  strength=10.0):
538 
539  if representation is not None and selection_tuple is not None:
540  self.m = representation.prot.get_model()
541  particles = IMP.pmi.tools.select_by_tuple(
542  representation,
543  selection_tuple,
544  resolution=1)
545 
546  elif objects is not None:
547  particles = IMP.pmi.tools.input_adaptor(objects,1,flatten=True)
548  self.m = particles[0].get_model()
549 
550  self.rs = IMP.RestraintSet(self.m, "Angles")
551  self.weight = 1
552  self.label = "None"
553  self.pairslist = []
554 
555  if stringsequence is None:
556  stringsequence = "T" * (len(particles) - 3)
557 
558  for n, ps in enumerate(IMP.pmi.tools.sublist_iterator(particles, 4, 4)):
559  quadruplet = []
560  if len(ps) != 4:
561  raise ValueError("wrong length of quadruplet")
562  for p in ps:
564  raise TypeError("not a residue")
565  else:
566  quadruplet.append(p)
567  dihedraltype = stringsequence[n]
568  if dihedraltype == "C":
569  anglemin = -20.0
570  anglemax = 20.0
572  (pi * anglemin / 180.0,
573  pi * anglemax / 180.0),
574  strength)
575  print("ResidueDihedralRestraint: adding a CYS restraint between %s %s %s %s" % (quadruplet[0].get_name(), quadruplet[1].get_name(),
576  quadruplet[2].get_name(), quadruplet[3].get_name()))
577  if dihedraltype == "T":
578  anglemin = 180 - 70.0
579  anglemax = 180 + 70.0
581  (pi * anglemin / 180.0,
582  pi * anglemax / 180.0),
583  strength)
584  print("ResidueDihedralRestraint: adding a TRANS restraint between %s %s %s %s" % (quadruplet[0].get_name(), quadruplet[1].get_name(),
585  quadruplet[2].get_name(), quadruplet[3].get_name()))
586  self.rs.add_restraint(
587  IMP.core.DihedralRestraint(self.m,ts,
588  quadruplet[0],
589  quadruplet[1],
590  quadruplet[2],
591  quadruplet[3]))
592  self.pairslist.append(
593  IMP.ParticlePair(quadruplet[0], quadruplet[3]))
594  self.pairslist.append(
595  IMP.ParticlePair(quadruplet[3], quadruplet[0]))
596 
597  def set_label(self, label):
598  self.label = label
599  self.rs.set_name(label)
600  for r in self.rs.get_restraints():
601  r.set_name(label)
602 
603  def add_to_model(self):
604  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
605 
606  def get_restraint(self):
607  return self.rs
608 
609  def set_weight(self, weight):
610  self.weight = weight
611  self.rs.set_weight(weight)
612 
613  def get_excluded_pairs(self):
614  return self.pairslist
615 
616  def get_output(self):
617  self.m.update()
618  output = {}
619  score = self.weight * self.rs.unprotected_evaluate(None)
620  output["_TotalScore"] = str(score)
621  output["ResidueDihedralRestraint_" + self.label] = str(score)
622  return output
623 
624 class SecondaryStructure(object):
625  """Experimental, requires isd_emxl for now"""
626  def __init__(self,
627  representation,
628  selection_tuple,
629  ssstring,
630  mixture=False,
631  nativeness=1.0,
632  kt_caff=0.1):
633 
634  if no_isd_emxl:
635  raise ValueError("IMP.isd_emxl is needed")
636 
637  # check that the secondary structure string
638  # is compatible with the ssstring
639 
640  self.particles = IMP.pmi.tools.select_by_tuple(
641  representation,
642  selection_tuple,
643  resolution=1)
644  self.m = representation.prot.get_model()
645  self.dihe_dict = {}
646  self.ang_dict = {}
647  self.do_mix = {}
648  self.anglfilename = IMP.isd_emxl.get_data_path("CAAngleRestraint.dat")
649  self.dihefilename = IMP.isd_emxl.get_data_path(
650  "CADihedralRestraint.dat")
651  self.nativeness = nativeness
652  self.kt_caff = kt_caff
653  self.anglrs = IMP.RestraintSet(self.m, "Angles")
654  self.dihers = IMP.RestraintSet(self.m, "Dihedrals")
655  self.bondrs = IMP.RestraintSet(self.m, "Bonds")
656  self.label = "None"
657 
658  if len(self.particles) != len(ssstring):
659  print(len(self.particles), len(ssstring))
660  print("SecondaryStructure: residue range and SS string incompatible")
661  self.ssstring = ssstring
662 
663  (bondrslist, anglrslist, diherslist,
664  pairslist) = self.get_CA_force_field()
665  self.pairslist = pairslist
666 
667  # print anglrslist, diherslist, bondrslist, self.particles
668  self.anglrs.add_restraints(anglrslist)
669  self.dihers.add_restraints(diherslist)
670  self.bondrs.add_restraints(bondrslist)
671 
672  def set_label(self, label):
673  self.label = label
674 
675  def add_to_model(self):
676  IMP.pmi.tools.add_restraint_to_model(self.m, self.anglrs)
677  IMP.pmi.tools.add_restraint_to_model(self.m, self.dihers)
678  IMP.pmi.tools.add_restraint_to_model(self.m, self.bondrs)
679 
680  def get_CA_force_field(self):
681  bondrslist = []
682  anglrslist = []
683  diherslist = []
684  pairslist = []
685  # add bonds
686  for res in range(0, len(self.particles) - 1):
687 
688  ps = self.particles[res:res + 2]
689  pairslist.append(IMP.ParticlePair(ps[0], ps[1]))
690  pairslist.append(IMP.ParticlePair(ps[1], ps[0]))
691  br = self.get_distance_restraint(ps[0], ps[1], 3.78, 416.0)
692  br.set_name('Bond_restraint')
693  bondrslist.append(br)
694  # add dihedrals
695  for res in range(0, len(self.particles) - 4):
696 
697  # if res not in dihe_dict: continue
698  # get the appropriate parameters
699  # get the particles
700  ps = self.particles[res:res + 5]
701  [phi0,
702  phi1,
703  score_dih] = self.read_potential_dihedral(
704  self.ssstring[res:res + 4],
705  True)
706  pairslist.append(IMP.ParticlePair(ps[0], ps[3]))
707  pairslist.append(IMP.ParticlePair(ps[3], ps[0]))
708  pairslist.append(IMP.ParticlePair(ps[1], ps[4]))
709  pairslist.append(IMP.ParticlePair(ps[4], ps[1]))
710  dr = IMP.isd_emxl.CADihedralRestraint(
711  ps[0],
712  ps[1],
713  ps[2],
714  ps[3],
715  ps[4],
716  phi0,
717  phi1,
718  score_dih)
719  dr.set_name('Dihedral restraint')
720  diherslist.append(dr)
721  # add angles
722  for res in range(0, len(self.particles) - 2):
723  ps = self.particles[res:res + 3]
724  [psi, score_ang] = self.read_potential_angle(
725  self.ssstring[res:res + 2], True)
726  pairslist.append(IMP.ParticlePair(ps[0], ps[2]))
727  pairslist.append(IMP.ParticlePair(ps[2], ps[0]))
728  dr = IMP.isd_emxl.CAAngleRestraint(
729  ps[0],
730  ps[1],
731  ps[2],
732  psi,
733  score_ang)
734  dr.set_name('Angle restraint')
735  anglrslist.append(dr)
736  return (bondrslist, anglrslist, diherslist, pairslist)
737 
738  def read_potential_dihedral(self, string, mix=False):
739  # read potentials for dihedral
740  score_dih = []
741  phi0 = []
742  phi1 = []
743  for i in range(0, 36):
744  phi0.append(i * 10.0 / 180.0 * pi)
745  phi1.append(i * 10.0 / 180.0 * pi)
746  for j in range(0, 36):
747  score_dih.append(0.0)
748  # open file
749  if not mix:
750  f = open(self.dihefilename, 'r')
751  for line in f.readlines():
752  riga = (line.strip()).split()
753  if (len(riga) == 4 and riga[0] == string):
754  ii = int(float(riga[1]) / 10.0)
755  jj = int(float(riga[2]) / 10.0)
756  score_dih[ii * len(phi0) + jj] = - \
757  self.kt_caff * self.log(float(riga[3]))
758  f.close()
759  if mix:
760  # mix random coil and native secondary structure
761  counts = []
762  for i in range(0, 36):
763  for j in range(0, 36):
764  counts.append(1.0)
765  f = open(self.dihefilename, 'r')
766  for line in f.readlines():
767  riga = (line.strip()).split()
768  if (len(riga) == 4 and riga[0] == string):
769  ii = int(float(riga[1]) / 10.0)
770  jj = int(float(riga[2]) / 10.0)
771  counts[ii * len(phi0) + jj] += self.nativeness * \
772  float(riga[3])
773  if (len(riga) == 4 and riga[0] == "-----"):
774  ii = int(float(riga[1]) / 10.0)
775  jj = int(float(riga[2]) / 10.0)
776  counts[ii * len(phi0) + jj] += (1.0 - self.nativeness) * \
777  float(riga[3])
778  f.close()
779  for i in range(len(counts)):
780  score_dih[i] = -self.kt_caff * self.log(counts[i])
781  return [phi0, phi1, score_dih]
782 
783  def read_potential_angle(self, string, mix=False):
784  # read potentials for angles
785  score_ang = []
786  psi = []
787  for i in range(0, 180):
788  psi.append(i / 180.0 * pi)
789  score_ang.append(0.0)
790  # read file
791  if not mix:
792  f = open(self.anglfilename, 'r')
793  for line in f.readlines():
794  riga = (line.strip()).split()
795  if (len(riga) == 3 and riga[0] == string):
796  ii = int(riga[1])
797  score_ang[ii] = -self.kt_caff * self.log(float(riga[2]))
798  f.close()
799  if mix:
800  # mix random coil and native secondary structure
801  counts = []
802  for i in range(0, 180):
803  counts.append(1.0)
804 
805  f = open(self.anglfilename, 'r')
806  for line in f.readlines():
807  riga = (line.strip()).split()
808  if (len(riga) == 3 and riga[0] == string):
809  ii = int(riga[1])
810  counts[ii] += self.nativeness * float(riga[2])
811  if (len(riga) == 3 and riga[0] == "---"):
812  ii = int(riga[1])
813  counts[ii] += (1.0 - self.nativeness) * float(riga[2])
814  f.close()
815  for i in range(0, 180):
816  score_ang[i] = -self.kt_caff * self.log(counts[i])
817  return [psi, score_ang]
818 
819  def get_excluded_pairs(self):
820  return self.pairslist
821 
822  def get_restraint(self):
823  tmprs = IMP.RestraintSet(self.m, 'tmp')
824  tmprs.add_restraint(self.anglrs)
825  tmprs.add_restraint(self.dihers)
826  tmprs.add_restraint(self.bondrs)
827  return tmprs
828 
829  def get_distance_restraint(self, p0, p1, d0, kappa):
830  h = IMP.core.Harmonic(d0, kappa)
832  pr = IMP.core.PairRestraint(self.m, dps, IMP.ParticlePair(p0, p1))
833  return pr
834 
835  def get_output(self):
836  output = {}
837  self.m.update()
838  score_angle = self.anglrs.unprotected_evaluate(None)
839  score_dihers = self.dihers.unprotected_evaluate(None)
840  score_bondrs = self.bondrs.unprotected_evaluate(None)
841  output["_TotalScore"] = str(score_angle + score_dihers + score_bondrs)
842 
843  output["SecondaryStructure_Angles_" + self.label] = str(score_angle)
844  output["SecondaryStructure_Dihedrals_" +
845  self.label] = str(score_dihers)
846  output["SecondaryStructure_Bonds_" + self.label] = str(score_bondrs)
847  return output
848 
849 
851  """Add harmonic restraints between all pairs
852  """
853  def __init__(self,representation=None,
854  selection_tuples=None,
855  resolution=1,
856  strength=0.01,
857  dist_cutoff=10.0,
858  ca_only=True,
859  hierarchy=None):
860  """Constructor
861  @param representation Representation object
862  @param selection_tuples Selecting regions for the restraint [[start,stop,molname,copy_index=0],...]
863  @param resolution Resolution for applying restraint
864  @param strength Bond strength
865  @param dist_cutoff Cutoff for making restraints
866  @param ca_only Selects only CAlphas. Only matters if resolution=0.
867  @param hierarchy Root hierarchy to select from, use this instead of representation
868  """
869 
870  particles = []
871  if representation is None and hierarchy is not None:
872  self.m = hierarchy.get_model()
873  for st in selection_tuples:
874  copy_index=0
875  if len(st)>3:
876  copy_index=st[3]
877  if not ca_only:
878  sel = IMP.atom.Selection(hierarchy,molecule=st[2],residue_indexes=range(st[0],st[1]+1),
879  copy_index=copy_index)
880  else:
881  sel = IMP.atom.Selection(hierarchy,molecule=st[2],residue_indexes=range(st[0],st[1]+1),
882  copy_index=copy_index,
883  atom_type=IMP.atom.AtomType("CA"))
884  particles+=sel.get_selected_particles()
885  elif representation is not None and type(representation)==IMP.pmi.representation.Representation:
886  self.m = representation.mdl
887  for st in selection_tuples:
888  print('selecting with',st)
889  for p in IMP.pmi.tools.select_by_tuple(representation,st,resolution=resolution):
890  if (resolution==0 and ca_only and IMP.atom.Atom(p).get_atom_type()!=IMP.atom.AtomType("CA")):
891  continue
892  else:
893  particles.append(p.get_particle())
894  else:
895  raise Exception("must pass representation or hierarchy")
896 
897  self.weight = 1
898  self.label = "None"
899  self.pairslist = []
900 
901  # create score
902  self.rs = IMP.pmi.create_elastic_network(particles,dist_cutoff,strength)
903  for r in self.rs.get_restraints():
904  a1,a2 = r.get_inputs()
905  self.pairslist.append(IMP.ParticlePair(a1,a2))
906  self.pairslist.append(IMP.ParticlePair(a2,a1))
907  print('ElasticNetwork: created',self.rs.get_number_of_restraints(),'restraints')
908 
909  def set_label(self, label):
910  self.label = label
911  self.rs.set_name(label)
912  for r in self.rs.get_restraints():
913  r.set_name(label)
914 
915  def add_to_model(self):
916  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
917 
918  def get_restraint(self):
919  return self.rs
920 
921  def set_weight(self, weight):
922  self.weight = weight
923  self.rs.set_weight(weight)
924 
925  def get_excluded_pairs(self):
926  return self.pairslist
927 
928  def get_output(self):
929  self.m.update()
930  output = {}
931  score = self.weight * self.rs.unprotected_evaluate(None)
932  output["_TotalScore"] = str(score)
933  output["ElasticNetworkRestraint_" + self.label] = str(score)
934  return output
935 
936 
938  """ Enable CHARMM force field """
939  def __init__(self,
940  root=None,
941  ff_temp=300.0,
942  zone_ps=None,
943  zone_size=10.0,
944  enable_nonbonded=True,
945  enable_bonded=True,
946  zone_nonbonded=False,
947  representation=None):
948  """Setup the CHARMM restraint on a selection. Expecting atoms.
949  @param root The node at which to apply the restraint
950  @param ff_temp The temperature of the force field
951  @param zone_ps Create a zone around this set of particles
952  Automatically includes the entire residue (incl. backbone)
953  @param zone_size The size for looking for neighbor residues
954  @param enable_nonbonded Allow the repulsive restraint
955  @param enable_bonded Allow the bonded restraint
956  @param zone_nonbonded EXPERIMENTAL: exclude from nonbonded all sidechains that aren't in zone!
957  @param representation Legacy representation object
958  """
959 
960  kB = (1.381 * 6.02214) / 4184.0
961  if representation is not None:
962  root = representation.prot
963 
964  self.mdl = root.get_model()
965  self.bonds_rs = IMP.RestraintSet(self.mdl, 1.0 / (kB * ff_temp), 'BONDED')
966  self.nonbonded_rs = IMP.RestraintSet(self.mdl, 1.0 / (kB * ff_temp), 'NONBONDED')
967  self.weight = 1.0
968  self.label = ""
969 
970  ### setup topology and bonds etc
971  if enable_bonded:
973  topology = ff.create_topology(root)
974  topology.apply_default_patches()
975  topology.setup_hierarchy(root)
976  if zone_ps is not None:
977  limit_to_ps=IMP.pmi.topology.get_particles_within_zone(
978  root,
979  zone_ps,
980  zone_size,
981  entire_residues=True,
982  exclude_backbone=False)
984  topology,
985  limit_to_ps)
986  self.ps = limit_to_ps
987  else:
988  r = IMP.atom.CHARMMStereochemistryRestraint(root, topology)
989  self.ps = IMP.core.get_leaves(root)
990  print('init bonds score',r.unprotected_evaluate(None))
991  self.bonds_rs.add_restraint(r)
992  ff.add_radii(root)
993  ff.add_well_depths(root)
994 
995  atoms = IMP.atom.get_by_type(root,IMP.atom.ATOM_TYPE)
996  ### non-bonded forces
997  if enable_nonbonded:
998  if (zone_ps is not None) and zone_nonbonded:
999  print('stereochemistry: zone_nonbonded is True')
1000  # atoms list should only include backbone plus zone_ps!
1001  backbone_types=['C','N','CB','O']
1002  sel = IMP.atom.Selection(root,atom_types=[IMP.atom.AtomType(n)
1003  for n in backbone_types])
1004  backbone_atoms = sel.get_selected_particles()
1005  #x = list(set(backbone_atoms)|set(limit_to_ps))
1006  sel_ps=IMP.pmi.topology.get_particles_within_zone(
1007  root,
1008  zone_ps,
1009  zone_size,
1010  entire_residues=True,
1011  exclude_backbone=True)
1012 
1014  IMP.container.ListSingletonContainer(backbone_atoms),
1016  4.0)
1017  else:
1018  cont = IMP.container.ListSingletonContainer(self.mdl,atoms)
1019  self.nbl = IMP.container.ClosePairContainer(cont, 4.0)
1020  if enable_bonded:
1021  self.nbl.add_pair_filter(r.get_full_pair_filter())
1022  pairscore = IMP.isd.RepulsiveDistancePairScore(0,1)
1023  pr=IMP.container.PairsRestraint(pairscore, self.nbl)
1024  self.nonbonded_rs.add_restraint(pr)
1025  print('CHARMM is set up')
1026 
1027  def set_label(self, label):
1028  self.label = label
1029  self.rs.set_name(label)
1030  for r in self.rs.get_restraints():
1031  r.set_name(label)
1032 
1033  def add_to_model(self):
1034  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.bonds_rs)
1035  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.nonbonded_rs)
1036 
1037  def get_restraint(self):
1038  return self.rs
1039 
1040  def get_close_pair_container(self):
1041  return self.nbl
1042 
1043  def set_weight(self, weight):
1044  self.weight = weight
1045  self.rs.set_weight(weight)
1046 
1047  def get_output(self):
1048  self.mdl.update()
1049  output = {}
1050  bonds_score = self.weight * self.bonds_rs.unprotected_evaluate(None)
1051  nonbonded_score = self.weight * self.nonbonded_rs.unprotected_evaluate(None)
1052  score=bonds_score+nonbonded_score
1053  output["_TotalScore"] = str(score)
1054  output["CHARMM_BONDS"] = str(bonds_score)
1055  output["CHARMM_NONBONDED"] = str(nonbonded_score)
1056  return output
1057 
1058 
1060  """Add bonds and improper dihedral restraints for the CBs
1061  """
1062  def __init__(
1063  self, rnums, representation, selection_tuple, strength=10.0, kappa=1.0,
1064  jitter_angle=0.0, jitter_improper=0.0):
1065  '''
1066  need to add:
1067  ca-ca bond
1068  ca-cb is a constraint, no restraint needed
1069  ca-ca-ca
1070  cb-ca-ca-cb
1071  '''
1072 
1073  self.m = representation.prot.get_model()
1074  self.rs = IMP.RestraintSet(self.m, "PseudoAtomic")
1075  self.rset_angles = IMP.RestraintSet(self.m, "PseudoAtomic_Angles")
1076  self.rset_bonds = IMP.RestraintSet(self.m, "PseudoAtomic_Bonds")
1077  self.weight = 1
1078  self.label = "None"
1079  self.pairslist = []
1080 
1081  # residues=IMP.pmi.tools.select_by_tuple(representation,selection_tuple,resolution=1)
1082  for rnum in rnums:
1083  ca, cb = self.get_ca_cb(
1084  IMP.pmi.tools.select_by_tuple(representation,
1085  (rnum, rnum, 'chainA'), resolution=0))
1086  if not cb is None:
1087  nter = False
1088  cter = False
1089  ca_prev, cb_prev = self.get_ca_cb(
1090  IMP.pmi.tools.select_by_tuple(representation,
1091  (rnum - 1, rnum - 1, 'chainA'), resolution=0))
1092  ca_next, cb_next = self.get_ca_cb(
1093  IMP.pmi.tools.select_by_tuple(representation,
1094  (rnum + 1, rnum + 1, 'chainA'), resolution=0))
1095  if ca_prev is None:
1096  nter = True
1097  else:
1098  if ca_next is None:
1099  cter = True
1100  else:
1101  if (nter and cter):
1102  continue
1103 
1104  # adding a bond restraint between CA and CB
1105  # h=IMP.core.Harmonic(6.0,kappa)
1106  # dps=IMP.core.DistancePairScore(h)
1107  # pr=IMP.core.PairRestraint(dps,IMP.ParticlePair(ca,cb))
1108  # self.pairslist.append((ca,cb))
1109  # self.rset_bonds.add_restraint(pr)
1110 
1111  # creating improper dihedral restraint
1112  # hus=IMP.core.Harmonic(2.09,kappa)
1114  2.09 +
1115  jitter_angle /
1116  0.5,
1117  kappa)
1119  2.09 -
1120  jitter_angle /
1121  0.5,
1122  kappa)
1123  if not nter:
1124  # ar13=IMP.core.AngleRestraint(hus,ca_prev,ca,cb)
1125  # self.rset_angles.add_restraint(ar13)
1126  ar13u = IMP.core.AngleRestraint(hupp, ca_prev, ca, cb)
1127  ar13l = IMP.core.AngleRestraint(hlow, ca_prev, ca, cb)
1128  self.rset_angles.add_restraint(ar13u)
1129  self.rset_angles.add_restraint(ar13l)
1130  if not cter:
1131  # ar23=IMP.core.AngleRestraint(hus,ca_next,ca,cb)
1132  # self.rset_angles.add_restraint(ar23)
1133  ar23u = IMP.core.AngleRestraint(hupp, ca_next, ca, cb)
1134  ar23l = IMP.core.AngleRestraint(hlow, ca_next, ca, cb)
1135  self.rset_angles.add_restraint(ar23u)
1136  self.rset_angles.add_restraint(ar23l)
1137  if not nter and not cter:
1138  # hus2=IMP.core.Harmonic(0,kappa)
1139  # idr=IMP.core.DihedralRestraint(hus2,ca,ca_prev,ca_next,cb)
1140  # self.rset_angles.add_restraint(idr)
1141 
1142  hus2upp = IMP.core.HarmonicUpperBound(
1143  jitter_improper,
1144  kappa)
1145  hus2low = IMP.core.HarmonicLowerBound(
1146  -
1147  jitter_improper,
1148  kappa)
1150  hus2upp,
1151  ca,
1152  ca_prev,
1153  ca_next,
1154  cb)
1156  hus2low,
1157  ca,
1158  ca_prev,
1159  ca_next,
1160  cb)
1161  self.rset_angles.add_restraint(idru)
1162  self.rset_angles.add_restraint(idrl)
1163  self.rs.add_restraint(self.rset_bonds)
1164  self.rs.add_restraint(self.rset_angles)
1165 
1166  def get_ca_cb(self, atoms):
1167  ca = None
1168  cb = None
1169  for a in atoms:
1170  if IMP.atom.Atom(a).get_atom_type() == IMP.atom.AtomType("CA"):
1171  ca = a.get_particle()
1172  elif IMP.atom.Atom(a).get_atom_type() == IMP.atom.AtomType("CB"):
1173  cb = a.get_particle()
1174  return ca, cb
1175 
1176  def set_label(self, label):
1177  self.label = label
1178  self.rs.set_name(label)
1179  for r in self.rs.get_restraints():
1180  r.set_name(label)
1181 
1182  def add_to_model(self):
1183  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1184 
1185  def get_restraint(self):
1186  return self.rs
1187 
1188  def set_weight(self, weight):
1189  self.weight = weight
1190  self.rs.set_weight(weight)
1191 
1192  def get_excluded_pairs(self):
1193  return self.pairslist
1194 
1195  def get_output(self):
1196  self.m.update()
1197  output = {}
1198  score = self.weight * self.rs.unprotected_evaluate(None)
1199  output["_TotalScore"] = str(score)
1200  output["PseudoAtomicRestraint_" + self.label] = str(score)
1201  return output
1202 
1203 class SymmetryRestraint(object):
1204  """Create harmonic restraints between the reference and (transformed) clones.
1205  \note Wraps IMP::core::TransformedDistancePairScore with an IMP::core::Harmonic
1206  """
1207  def __init__(self,
1208  references,
1209  clones_list,
1210  transforms,
1211  label='',
1212  strength=10.0,
1213  ca_only=False):
1214  """Constructor
1215  @param references Can be one of the following inputs:
1216  IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set of them
1217  @param clones_list List of lists of the above
1218  @param transforms Transforms moving each selection to the first selection
1219  @param label Label for output
1220  @param strength The elastic bond strength
1221  @param ca_only Optionally select so only CAlpha particles are used
1222  """
1223 
1224  refs = IMP.pmi.tools.input_adaptor(references,flatten=True)
1225  self.mdl = refs[0].get_model()
1226  self.rs = IMP.RestraintSet(self.mdl, "Symmetry")
1227  self.weight = 1
1228  self.label = label
1229  if len(clones_list)!=len(transforms):
1230  raise Exception('Error: There should be as many clones as transforms')
1231 
1232  harmonic = IMP.core.Harmonic(0.,strength)
1233  for tmp_clones,trans in zip(clones_list,transforms):
1234  clones = IMP.pmi.tools.input_adaptor(tmp_clones,flatten=True)
1235  if len(clones)!=len(refs):
1236  raise Exception("Error: len(references)!=len(clones)")
1237  pair_score = IMP.core.TransformedDistancePairScore(harmonic,trans)
1238  for p0,p1 in zip(refs,clones):
1239  if not ca_only or (IMP.atom.Atom(p0).get_atom_type()==IMP.atom.AtomType("CA")
1240  and IMP.atom.Atom(p1).get_atom_type()==IMP.atom.AtomType("CA")):
1241  r = IMP.core.PairRestraint(self.mdl, pair_score, [p0.get_particle_index(),
1242  p1.get_particle_index()])
1243  self.rs.add_restraint(r)
1244  print('created symmetry network with',self.rs.get_number_of_restraints(),'restraints')
1245 
1246  def set_label(self, label):
1247  self.label = label
1248  self.rs.set_name(label)
1249  for r in self.rs.get_restraints():
1250  r.set_name(label)
1251 
1252  def add_to_model(self):
1253  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.rs)
1254 
1255  def get_restraint(self):
1256  return self.rs
1257 
1258  def set_weight(self, weight):
1259  self.weight = weight
1260  self.rs.set_weight(weight)
1261 
1262  def get_excluded_pairs(self):
1263  return self.pairslist
1264 
1265  def get_output(self):
1266  self.mdl.update()
1267  output = {}
1268  score = self.weight * self.rs.unprotected_evaluate(None)
1269  output["SymmetryRestraint_" + self.label] = str(score)
1270  output["_TotalScore"] = str(score)
1271  return output
1272 
1273 
1274 class FusionRestraint(object):
1275  """ This class creates a restraint between the termini two polypeptides, to simulate the sequence connectivity. """
1276  def __init__(self,
1277  nterminal,
1278  cterminal,
1279  scale=1.0,
1280  disorderedlength=False,
1281  upperharmonic=True,
1282  resolution=1,
1283  label="None"):
1284  """
1285  @param nterminal - single PMI2 Hierarchy/molecule at the nterminal
1286  @param cterminal - single PMI2 Hierarchy/molecule at the cterminal
1287  @param scale Scale the maximal distance between the beads by this factor when disorderedlength is False.
1288  The maximal distance is calculated as ((float(residuegap) + 1.0) * 3.6) * scale.
1289  @param disorderedlength - This flag uses either disordered length
1290  calculated for random coil peptides (True) or zero
1291  surface-to-surface distance between beads (False)
1292  as optimal distance for the sequence connectivity
1293  restraint.
1294  @param upperharmonic - This flag uses either harmonic (False)
1295  or upperharmonic (True) in the intra-pair
1296  connectivity restraint.
1297  @param resolution - The resolution to connect things at - only used if you pass PMI objects
1298  @param label - A string to identify this restraint in the output/stat file
1299  """
1300  self.label = label
1301  self.weight = 1.0
1302  ssn=IMP.pmi.tools.get_sorted_segments(nterminal)
1303  ssc=IMP.pmi.tools.get_sorted_segments(cterminal)
1304  nter_lastres=ssn[-1][1]
1305  cter_firstres=ssc[0][0]
1306  self.m = nter_lastres.get_model()
1307 
1308  self.kappa = 10 # spring constant used for the harmonic restraints
1309 
1310  optdist = (3.6) * scale
1311  if upperharmonic: # default
1312  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1313  else:
1314  hu = IMP.core.Harmonic(optdist, self.kappa)
1316 
1317  pt0 = nter_lastres.get_particle()
1318  pt1 = cter_firstres.get_particle()
1319  r = IMP.core.PairRestraint(self.m, dps, (pt0.get_index(), pt1.get_index()))
1320  self.rs = IMP.RestraintSet(self.m, "fusion_restraint")
1321  print("Adding fusion connectivity restraint between", pt0.get_name(), " and ", pt1.get_name(), 'of distance', optdist)
1322  self.rs.add_restraint(r)
1323 
1324  def set_label(self, label):
1325  self.label = label
1326 
1327  def get_weight(self):
1328  return self.weight
1329 
1330  def add_to_model(self):
1331  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1332 
1333  def get_restraint(self):
1334  return self.rs
1335 
1336  def set_weight(self, weight):
1337  self.weight = weight
1338  self.rs.set_weight(weight)
1339 
1340  def get_output(self):
1341  self.m.update()
1342  output = {}
1343  score = self.weight * self.rs.unprotected_evaluate(None)
1344  output["_TotalScore"] = str(score)
1345  output["FusionRestraint_" + self.label] = str(score)
1346  return output
1347 
1348  def evaluate(self):
1349  return self.weight * self.rs.unprotected_evaluate(None)
1350 
1351 
1352 
1353 
1355 
1356  """Restrain the dihedral between planes defined by three particles.
1357 
1358  This restraint is useful for restraining the twist of a string of
1359  more or less identical rigid bodies, so long as the curvature is mild.
1360  """
1361 
1362  def __init__(self, particle_triplets, angle=0., k=1., label=None,
1363  weight=1.):
1364  """Constructor
1365  @param particle_triplets List of lists of 3 particles. Each triplet
1366  defines a plane. Dihedrals of adjacent planes
1367  in list are scored.
1368  @param angle Angle of plane dihedral in degrees
1369  @param k Strength of restraint
1370  @param label Label for output
1371  @param weight Weight of restraint
1372  \note Particles defining planes should be rigid and more or less
1373  parallel for proper behavior
1374  """
1375  m = particle_triplets[0][0].get_model()
1376  super(PlaneDihedralRestraint, self).__init__(m, label=label,
1377  weight=weight)
1378 
1379  angle = pi * angle / 180.
1380  ds = IMP.core.Cosine(.5 * k, 1., -angle)
1381  for i, t1 in enumerate(particle_triplets[:-1]):
1382  t2 = particle_triplets[i + 1]
1383  q1 = [t1[1], t1[0], t2[0], t2[1]]
1384  q2 = [t1[2], t1[0], t2[0], t2[2]]
1385  self.rs.add_restraint(IMP.core.DihedralRestraint(self.m, ds, *q1))
1386  self.rs.add_restraint(IMP.core.DihedralRestraint(self.m, 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:155
Enforce CHARMM stereochemistry on the given Hierarchy.
A member of a rigid body, it has internal (local) coordinates.
Definition: rigid_bodies.h:521
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:522
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
This class creates a restraint between the termini two polypeptides, to simulate the sequence connect...
Add harmonic restraints between all pairs.
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:632
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.
Representation of the system.
Object used to hold a set of restraints.
Definition: RestraintSet.h:36
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: tools.py:1618
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)
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: 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:50
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.
A decorator for a residue.
Definition: Residue.h:134
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
This class creates a restraint between consecutive TempResidue objects OR an entire PMI MOlecule obje...
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 select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:724
Experimental, requires isd_emxl for now.
Set up the representation of all proteins and nucleic acid macromolecules.
def get_sorted_segments
Returns sequence-sorted segments array, each containing the first particle the last particle and the ...
Definition: tools.py:1739
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:1061
A repulsive potential on the distance between two atoms.
Perform more efficient close pair finding when rigid bodies are involved.
Base class for PMI restraints, which wrap IMP.Restraint(s).
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:1165
Harmonic function (symmetric about the mean)
Definition: core/Harmonic.h:24
Restraint a set of residues to use ideal helix dihedrals and bonds.