IMP logo
IMP Reference Guide  2.6.1
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 class ConnectivityRestraint(object):
21  """ This class creates a restraint between consecutive TempResidue objects OR an entire
22  PMI MOlecule object. """
23 
24  def __init__(self,
25  objects,
26  scale=1.0,
27  disorderedlength=False,
28  upperharmonic=True,
29  resolution=1):
30  """
31  @param objects - a list of hierarchies, PMI TempResidues OR a single Molecule
32  @param scale - Riccardo knows what this is
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
37  restraint.
38  @param upperharmonic - This flag uses either harmonic (False)
39  or upperharmonic (True) in the intra-pair
40  connectivity restraint.
41  @parm resolution The resolution to connect things at - only used if you pass PMI objects
42  """
43  self.label = "None"
44  self.weight = 1.0
45 
46  hiers = IMP.pmi.tools.input_adaptor(objects,resolution)
47  if len(hiers)>1:
48  raise Exception("ConnectivityRestraint: only pass stuff from one Molecule, please")
49  hiers = hiers[0]
50 
51  self.kappa = 10 # No idea what this is
52  self.m = list(hiers)[0].get_model()
53  SortedSegments = []
54  self.rs = IMP.RestraintSet(self.m, "connectivity_restraint")
55  for h in hiers:
56  try:
57  start = IMP.atom.Hierarchy(h).get_children()[0]
58  except:
59  start = IMP.atom.Hierarchy(h)
60 
61  try:
62  end = IMP.atom.Hierarchy(h).get_children()[-1]
63  except:
64  end = IMP.atom.Hierarchy(h)
65 
66  startres = IMP.pmi.tools.get_residue_indexes(start)[0]
67  endres = IMP.pmi.tools.get_residue_indexes(end)[-1]
68  SortedSegments.append((start, end, startres))
69  SortedSegments = sorted(SortedSegments, key=itemgetter(2))
70 
71  # connect the particles
72  for x in range(len(SortedSegments) - 1):
73 
74  last = SortedSegments[x][1]
75  first = SortedSegments[x + 1][0]
76 
77  apply_restraint = True
78 
79  # Apply connectivity runless ALL of the following are true:
80  # - first and last both have RigidBodyMember decorators
81  # - first and last are both RigidMembers
82  # - first and last are part of the same RigidBody object
83 
84  # Check for both in a rigid body
87 
88  #Check if the rigid body objects for each particle are the same object.
89  # if so, skip connectivity restraint
90  if IMP.core.RigidBodyMember(first).get_rigid_body() == IMP.core.RigidBodyMember(last).get_rigid_body():
91  apply_restraint = False
92 
93  if apply_restraint:
94 
95  nreslast = len(IMP.pmi.tools.get_residue_indexes(last))
96  lastresn = IMP.pmi.tools.get_residue_indexes(last)[-1]
97  nresfirst = len(IMP.pmi.tools.get_residue_indexes(first))
98  firstresn = IMP.pmi.tools.get_residue_indexes(first)[0]
99 
100  residuegap = firstresn - lastresn - 1
101  if disorderedlength and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
102  # calculate the distance between the sphere centers using Kohn
103  # PNAS 2004
104  optdist = sqrt(5 / 3) * 1.93 * \
105  (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
106  # optdist2=sqrt(5/3)*1.93*((nreslast)**0.6+(nresfirst)**0.6)/2
107  if upperharmonic:
108  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
109  else:
110  hu = IMP.core.Harmonic(optdist, self.kappa)
112  else: # default
113  optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
114  if upperharmonic: # default
115  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
116  else:
117  hu = IMP.core.Harmonic(optdist, self.kappa)
119 
120  pt0 = last.get_particle()
121  pt1 = first.get_particle()
122  r = IMP.core.PairRestraint(self.m, dps, (pt0.get_index(), pt1.get_index()))
123 
124  print("Adding sequence connectivity restraint between", pt0.get_name(), " and ", pt1.get_name(), 'of distance', optdist)
125  self.rs.add_restraint(r)
126 
127 
128  def set_label(self, label):
129  self.label = label
130 
131  def get_weight(self):
132  return self.weight
133 
134  def add_to_model(self):
135  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
136 
137  def get_restraint(self):
138  return self.rs
139 
140  def set_weight(self, weight):
141  self.weight = weight
142  self.rs.set_weight(weight)
143 
144  def get_output(self):
145  self.m.update()
146  output = {}
147  score = self.weight * self.rs.unprotected_evaluate(None)
148  output["_TotalScore"] = str(score)
149  output["ConnectivityRestraint_" + self.label] = str(score)
150  return output
151 
153  """ Returns number of connectivity restraints """
154  return len(self.rs.get_restraints())
155 
156  def evaluate(self):
157  return self.weight * self.rs.unprotected_evaluate(None)
158 
159 class ExcludedVolumeSphere(object):
160  """A class to create an excluded volume restraint for a set of particles at a given resolution.
161  Can be initialized as a bipartite restraint between two sets of particles.
162  # Potential additional function: Variable resolution for each PMI object. Perhaps passing selection_tuples
163  with (PMI_object, resolution)
164  """
165 
166  def __init__(self,
167  representation=None,
168  included_objects=None,
169  other_objects=None,
170  resolution=1000,
171  kappa=1.0):
172  """Constructor.
173  @param representation DEPRECATED - just pass objects
174  @param included_objects Can be one of the following inputs:
175  IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set of them
176  @param other_objects Initializes a bipartite restraint between included_objects and other_objects
177  Same format as included_objects
178  @param resolution The resolution particles at which to impose the restraint.
179  By default, the coarsest particles will be chosen.
180  If a number is chosen, for each particle, the closest
181  resolution will be used (see IMP.atom.Selection).
182  @param kappa Restraint strength
183  """
184 
185  self.weight = 1.0
186  self.kappa = kappa
187  self.label = "None"
188  self.cpc = None
189  bipartite = False
190 
191  # gather IMP hierarchies from input objects
192  hierarchies = IMP.pmi.tools.input_adaptor(included_objects,
193  resolution,
194  flatten=True)
195  included_ps = []
196  if other_objects is not None:
197  bipartite = True
198  other_hierarchies = IMP.pmi.tools.input_adaptor(other_objects,
199  resolution,
200  flatten=True)
201  other_ps = []
202 
203  # perform selection
204  if representation is None:
205  if hierarchies is None:
206  raise Exception("Must at least pass included objects")
207  self.mdl = hierarchies[0].get_model()
208  included_ps = [h.get_particle() for h in hierarchies]
209  if bipartite:
210  other_ps = [h.get_particle() for h in other_hierarchies]
211  elif type(representation) is IMP.pmi.representation.Representation or \
212  type(representation) is IMP.pmi.representation.SimplifiedModel:
213  self.mdl = representation.m
214  included_ps = IMP.pmi.tools.select(
215  representation,
216  resolution=resolution,
217  hierarchies=hierarchies)
218  if bipartite:
219  other_particles = IMP.pmi.tools.select(
220  representation,
221  resolution=resolution,
222  hierarchies=other_hierarchies)
223  else:
224  raise Exception("Must pass Representation or included_objects")
225 
226  # setup score
227  self.rs = IMP.RestraintSet(self.mdl, 'excluded_volume')
228  ssps = IMP.core.SoftSpherePairScore(self.kappa)
230  lsa.add(IMP.get_indexes(included_ps))
231 
232  # setup close pair container
233  if not bipartite:
235  self.cpc = IMP.container.ClosePairContainer(lsa, 0.0, rbcpf, 10.0)
236  evr = IMP.container.PairsRestraint(ssps, self.cpc)
237  else:
238  other_lsa = IMP.container.ListSingletonContainer(self.mdl)
239  other_lsa.add_particles(other_ps)
241  lsa,
242  other_lsa,
243  0.0,
244  10.0)
245  evr = IMP.container.PairsRestraint(ssps, self.cpc)
246 
247  self.rs.add_restraint(evr)
248 
249  def add_excluded_particle_pairs(self, excluded_particle_pairs):
250  # add pairs to be filtered when calculating the score
251  lpc = IMP.container.ListPairContainer(self.mdl)
252  lpc.add(IMP.get_indexes(excluded_particle_pairs))
254  self.cpc.add_pair_filter(icpf)
255 
256  def set_label(self, label):
257  self.label = label
258 
259  def add_to_model(self):
260  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.rs)
261 
262  def get_restraint(self):
263  return self.rs
264 
265  def set_weight(self, weight):
266  self.weight = weight
267  self.rs.set_weight(weight)
268 
269  def get_output(self):
270  self.mdl.update()
271  output = {}
272  score = self.weight * self.rs.unprotected_evaluate(None)
273  output["_TotalScore"] = str(score)
274  output["ExcludedVolumeSphere_" + self.label] = str(score)
275  return output
276 
277  def evaluate(self):
278  return self.weight * self.rs.unprotected_evaluate(None)
279 
280 
281 class ResidueBondRestraint(object):
282  """ Add bond restraint between pair of consecutive
283  residues/beads to enforce the stereochemistry.
284  """
285  def __init__(self,
286  representation,
287  selection_tuple,
288  distance=3.78,
289  strength=10.0,
290  jitter=None):
291  """Constructor
292  @param representation
293  @param selection_tuple Requested selection
294  @param distance Resting distance for restraint
295  @param strength Bond constant
296  @param jitter Defines the +- added to the optimal distance in the harmonic well restraint
297  used to increase the tolerance
298  """
299  self.m = representation.prot.get_model()
300  self.rs = IMP.RestraintSet(self.m, "Bonds")
301  self.weight = 1
302  self.label = "None"
303  self.pairslist = []
304 
305  particles = IMP.pmi.tools.select_by_tuple(
306  representation,
307  selection_tuple,
308  resolution=1)
309 
310  if not jitter:
311  ts = IMP.core.Harmonic(distance, strength)
312  else:
314  (distance - jitter, distance + jitter), strength)
315 
316  for ps in IMP.pmi.tools.sublist_iterator(particles, 2, 2):
317  pair = []
318  if len(ps) != 2:
319  raise ValueError("wrong length of pair")
320  for p in ps:
322  raise TypeError("not a residue")
323  else:
324  pair.append(p)
325  print("ResidueBondRestraint: adding a restraint between %s %s" % (pair[0].get_name(), pair[1].get_name()))
326  self.rs.add_restraint(
327  IMP.core.DistanceRestraint(self.m, ts, pair[0], pair[1]))
328  self.pairslist.append(IMP.ParticlePair(pair[0], pair[1]))
329  self.pairslist.append(IMP.ParticlePair(pair[1], pair[0]))
330 
331  def set_label(self, label):
332  self.label = label
333  self.rs.set_name(label)
334  for r in self.rs.get_restraints():
335  r.set_name(label)
336 
337  def add_to_model(self):
338  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
339 
340  def get_restraint(self):
341  return self.rs
342 
343  def set_weight(self, weight):
344  self.weight = weight
345  self.rs.set_weight(weight)
346 
347  def get_excluded_pairs(self):
348  return self.pairslist
349 
350  def get_output(self):
351  self.m.update()
352  output = {}
353  score = self.weight * self.rs.unprotected_evaluate(None)
354  output["_TotalScore"] = str(score)
355  output["ResidueBondRestraint_" + self.label] = str(score)
356  return output
357 
358 
359 class ResidueAngleRestraint(object):
360  """Add angular restraint between triplets of consecutive
361  residues/beads to enforce the stereochemistry.
362  """
363  def __init__(self,
364  representation,
365  selection_tuple,
366  anglemin=100.0,
367  anglemax=140.0,
368  strength=10.0):
369  self.m = representation.prot.get_model()
370  self.rs = IMP.RestraintSet(self.m, "Angles")
371  self.weight = 1
372  self.label = "None"
373  self.pairslist = []
374 
375  particles = IMP.pmi.tools.select_by_tuple(
376  representation,
377  selection_tuple,
378  resolution=1)
379 
381  (pi * anglemin / 180.0,
382  pi * anglemax / 180.0),
383  strength)
384 
385  for ps in IMP.pmi.tools.sublist_iterator(particles, 3, 3):
386  triplet = []
387  if len(ps) != 3:
388  raise ValueError("wrong length of triplet")
389  for p in ps:
391  raise TypeError("not a residue")
392  else:
393  triplet.append(p)
394  print("ResidueAngleRestraint: adding a restraint between %s %s %s" % (triplet[0].get_name(), triplet[1].get_name(), triplet[2].get_name()))
395  self.rs.add_restraint(
396  IMP.core.AngleRestraint(triplet[0].get_model(), ts,
397  triplet[0],
398  triplet[1],
399  triplet[2]))
400  self.pairslist.append(IMP.ParticlePair(triplet[0], triplet[2]))
401  self.pairslist.append(IMP.ParticlePair(triplet[2], triplet[0]))
402 
403  def set_label(self, label):
404  self.label = label
405  self.rs.set_name(label)
406  for r in self.rs.get_restraints():
407  r.set_name(label)
408 
409  def add_to_model(self):
410  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
411 
412  def get_restraint(self):
413  return self.rs
414 
415  def set_weight(self, weight):
416  self.weight = weight
417  self.rs.set_weight(weight)
418 
419  def get_excluded_pairs(self):
420  return self.pairslist
421 
422  def get_output(self):
423  self.m.update()
424  output = {}
425  score = self.weight * self.rs.unprotected_evaluate(None)
426  output["_TotalScore"] = str(score)
427  output["ResidueAngleRestraint_" + self.label] = str(score)
428  return output
429 
430 
432  """Add dihedral restraints between quadruplet of consecutive
433  residues/beads to enforce the stereochemistry.
434  Give as input a string of "C" and "T", meaning cys (0+-40) or trans (180+-40)
435  dihedral. The length of the string must be \#residue-3.
436  Without the string, the dihedral will be assumed trans.
437  """
438  def __init__(
439  self,
440  representation,
441  selection_tuple,
442  stringsequence=None,
443  strength=10.0):
444  self.m = representation.prot.get_model()
445  self.rs = IMP.RestraintSet(self.m, "Angles")
446  self.weight = 1
447  self.label = "None"
448  self.pairslist = []
449 
450  particles = IMP.pmi.tools.select_by_tuple(
451  representation,
452  selection_tuple,
453  resolution=1)
454 
455  if stringsequence is None:
456  stringsequence = "T" * (len(particles) - 3)
457 
458  for n, ps in enumerate(IMP.pmi.tools.sublist_iterator(particles, 4, 4)):
459  quadruplet = []
460  if len(ps) != 4:
461  raise ValueError("wrong length of quadruplet")
462  for p in ps:
464  raise TypeError("not a residue")
465  else:
466  quadruplet.append(p)
467  dihedraltype = stringsequence[n]
468  if dihedraltype == "C":
469  anglemin = -20.0
470  anglemax = 20.0
472  (pi * anglemin / 180.0,
473  pi * anglemax / 180.0),
474  strength)
475  print("ResidueDihedralRestraint: adding a CYS restraint between %s %s %s %s" % (quadruplet[0].get_name(), quadruplet[1].get_name(),
476  quadruplet[2].get_name(), quadruplet[3].get_name()))
477  if dihedraltype == "T":
478  anglemin = 180 - 70.0
479  anglemax = 180 + 70.0
481  (pi * anglemin / 180.0,
482  pi * anglemax / 180.0),
483  strength)
484  print("ResidueDihedralRestraint: adding a TRANS restraint between %s %s %s %s" % (quadruplet[0].get_name(), quadruplet[1].get_name(),
485  quadruplet[2].get_name(), quadruplet[3].get_name()))
486  self.rs.add_restraint(
488  quadruplet[0],
489  quadruplet[1],
490  quadruplet[2],
491  quadruplet[3]))
492  self.pairslist.append(
493  IMP.ParticlePair(quadruplet[0], quadruplet[3]))
494  self.pairslist.append(
495  IMP.ParticlePair(quadruplet[3], quadruplet[0]))
496 
497  def set_label(self, label):
498  self.label = label
499  self.rs.set_name(label)
500  for r in self.rs.get_restraints():
501  r.set_name(label)
502 
503  def add_to_model(self):
504  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
505 
506  def get_restraint(self):
507  return self.rs
508 
509  def set_weight(self, weight):
510  self.weight = weight
511  self.rs.set_weight(weight)
512 
513  def get_excluded_pairs(self):
514  return self.pairslist
515 
516  def get_output(self):
517  self.m.update()
518  output = {}
519  score = self.weight * self.rs.unprotected_evaluate(None)
520  output["_TotalScore"] = str(score)
521  output["ResidueDihedralRestraint_" + self.label] = str(score)
522  return output
523 
524 class SecondaryStructure(object):
525  """Experimental, requires isd_emxl for now"""
526  def __init__(self,
527  representation,
528  selection_tuple,
529  ssstring,
530  mixture=False,
531  nativeness=1.0,
532  kt_caff=0.1):
533 
534  if no_isd_emxl:
535  raise ValueError("IMP.isd_emxl is needed")
536 
537  # check that the secondary structure string
538  # is compatible with the ssstring
539 
540  self.particles = IMP.pmi.tools.select_by_tuple(
541  representation,
542  selection_tuple,
543  resolution=1)
544  self.m = representation.prot.get_model()
545  self.dihe_dict = {}
546  self.ang_dict = {}
547  self.do_mix = {}
548  self.anglfilename = IMP.isd_emxl.get_data_path("CAAngleRestraint.dat")
549  self.dihefilename = IMP.isd_emxl.get_data_path(
550  "CADihedralRestraint.dat")
551  self.nativeness = nativeness
552  self.kt_caff = kt_caff
553  self.anglrs = IMP.RestraintSet(self.m, "Angles")
554  self.dihers = IMP.RestraintSet(self.m, "Dihedrals")
555  self.bondrs = IMP.RestraintSet(self.m, "Bonds")
556  self.label = "None"
557 
558  if len(self.particles) != len(ssstring):
559  print(len(self.particles), len(ssstring))
560  print("SecondaryStructure: residue range and SS string incompatible")
561  self.ssstring = ssstring
562 
563  (bondrslist, anglrslist, diherslist,
564  pairslist) = self.get_CA_force_field()
565  self.pairslist = pairslist
566 
567  # print anglrslist, diherslist, bondrslist, self.particles
568  self.anglrs.add_restraints(anglrslist)
569  self.dihers.add_restraints(diherslist)
570  self.bondrs.add_restraints(bondrslist)
571 
572  def set_label(self, label):
573  self.label = label
574 
575  def add_to_model(self):
576  IMP.pmi.tools.add_restraint_to_model(self.m, self.anglrs)
577  IMP.pmi.tools.add_restraint_to_model(self.m, self.dihers)
578  IMP.pmi.tools.add_restraint_to_model(self.m, self.bondrs)
579 
580  def get_CA_force_field(self):
581  bondrslist = []
582  anglrslist = []
583  diherslist = []
584  pairslist = []
585  # add bonds
586  for res in range(0, len(self.particles) - 1):
587 
588  ps = self.particles[res:res + 2]
589  pairslist.append(IMP.ParticlePair(ps[0], ps[1]))
590  pairslist.append(IMP.ParticlePair(ps[1], ps[0]))
591  br = self.get_distance_restraint(ps[0], ps[1], 3.78, 416.0)
592  br.set_name('Bond_restraint')
593  bondrslist.append(br)
594  # add dihedrals
595  for res in range(0, len(self.particles) - 4):
596 
597  # if res not in dihe_dict: continue
598  # get the appropriate parameters
599  # get the particles
600  ps = self.particles[res:res + 5]
601  [phi0,
602  phi1,
603  score_dih] = self.read_potential_dihedral(
604  self.ssstring[res:res + 4],
605  True)
606  pairslist.append(IMP.ParticlePair(ps[0], ps[3]))
607  pairslist.append(IMP.ParticlePair(ps[3], ps[0]))
608  pairslist.append(IMP.ParticlePair(ps[1], ps[4]))
609  pairslist.append(IMP.ParticlePair(ps[4], ps[1]))
610  dr = IMP.isd_emxl.CADihedralRestraint(
611  ps[0],
612  ps[1],
613  ps[2],
614  ps[3],
615  ps[4],
616  phi0,
617  phi1,
618  score_dih)
619  dr.set_name('Dihedral restraint')
620  diherslist.append(dr)
621  # add angles
622  for res in range(0, len(self.particles) - 2):
623  ps = self.particles[res:res + 3]
624  [psi, score_ang] = self.read_potential_angle(
625  self.ssstring[res:res + 2], True)
626  pairslist.append(IMP.ParticlePair(ps[0], ps[2]))
627  pairslist.append(IMP.ParticlePair(ps[2], ps[0]))
628  dr = IMP.isd_emxl.CAAngleRestraint(
629  ps[0],
630  ps[1],
631  ps[2],
632  psi,
633  score_ang)
634  dr.set_name('Angle restraint')
635  anglrslist.append(dr)
636  return (bondrslist, anglrslist, diherslist, pairslist)
637 
638  def read_potential_dihedral(self, string, mix=False):
639  # read potentials for dihedral
640  score_dih = []
641  phi0 = []
642  phi1 = []
643  for i in range(0, 36):
644  phi0.append(i * 10.0 / 180.0 * pi)
645  phi1.append(i * 10.0 / 180.0 * pi)
646  for j in range(0, 36):
647  score_dih.append(0.0)
648  # open file
649  if not mix:
650  f = open(self.dihefilename, 'r')
651  for line in f.readlines():
652  riga = (line.strip()).split()
653  if (len(riga) == 4 and riga[0] == string):
654  ii = int(float(riga[1]) / 10.0)
655  jj = int(float(riga[2]) / 10.0)
656  score_dih[ii * len(phi0) + jj] = - \
657  self.kt_caff * self.log(float(riga[3]))
658  f.close()
659  if mix:
660  # mix random coil and native secondary structure
661  counts = []
662  for i in range(0, 36):
663  for j in range(0, 36):
664  counts.append(1.0)
665  f = open(self.dihefilename, 'r')
666  for line in f.readlines():
667  riga = (line.strip()).split()
668  if (len(riga) == 4 and riga[0] == string):
669  ii = int(float(riga[1]) / 10.0)
670  jj = int(float(riga[2]) / 10.0)
671  counts[ii * len(phi0) + jj] += self.nativeness * \
672  float(riga[3])
673  if (len(riga) == 4 and riga[0] == "-----"):
674  ii = int(float(riga[1]) / 10.0)
675  jj = int(float(riga[2]) / 10.0)
676  counts[ii * len(phi0) + jj] += (1.0 - self.nativeness) * \
677  float(riga[3])
678  f.close()
679  for i in range(len(counts)):
680  score_dih[i] = -self.kt_caff * self.log(counts[i])
681  return [phi0, phi1, score_dih]
682 
683  def read_potential_angle(self, string, mix=False):
684  # read potentials for angles
685  score_ang = []
686  psi = []
687  for i in range(0, 180):
688  psi.append(i / 180.0 * pi)
689  score_ang.append(0.0)
690  # read file
691  if not mix:
692  f = open(self.anglfilename, 'r')
693  for line in f.readlines():
694  riga = (line.strip()).split()
695  if (len(riga) == 3 and riga[0] == string):
696  ii = int(riga[1])
697  score_ang[ii] = -self.kt_caff * self.log(float(riga[2]))
698  f.close()
699  if mix:
700  # mix random coil and native secondary structure
701  counts = []
702  for i in range(0, 180):
703  counts.append(1.0)
704 
705  f = open(self.anglfilename, 'r')
706  for line in f.readlines():
707  riga = (line.strip()).split()
708  if (len(riga) == 3 and riga[0] == string):
709  ii = int(riga[1])
710  counts[ii] += self.nativeness * float(riga[2])
711  if (len(riga) == 3 and riga[0] == "---"):
712  ii = int(riga[1])
713  counts[ii] += (1.0 - self.nativeness) * float(riga[2])
714  f.close()
715  for i in range(0, 180):
716  score_ang[i] = -self.kt_caff * self.log(counts[i])
717  return [psi, score_ang]
718 
719  def get_excluded_pairs(self):
720  return self.pairslist
721 
722  def get_restraint(self):
723  tmprs = IMP.RestraintSet(self.m, 'tmp')
724  tmprs.add_restraint(self.anglrs)
725  tmprs.add_restraint(self.dihers)
726  tmprs.add_restraint(self.bondrs)
727  return tmprs
728 
729  def get_distance_restraint(self, p0, p1, d0, kappa):
730  h = IMP.core.Harmonic(d0, kappa)
732  pr = IMP.core.PairRestraint(self.m, dps, IMP.ParticlePair(p0, p1))
733  return pr
734 
735  def get_output(self):
736  output = {}
737  self.m.update()
738  score_angle = self.anglrs.unprotected_evaluate(None)
739  score_dihers = self.dihers.unprotected_evaluate(None)
740  score_bondrs = self.bondrs.unprotected_evaluate(None)
741  output["_TotalScore"] = str(score_angle + score_dihers + score_bondrs)
742 
743  output["SecondaryStructure_Angles_" + self.label] = str(score_angle)
744  output["SecondaryStructure_Dihedrals_" +
745  self.label] = str(score_dihers)
746  output["SecondaryStructure_Bonds_" + self.label] = str(score_bondrs)
747  return output
748 
749 
751  """Add harmonic restraints between all pairs
752  """
753  def __init__(self,representation=None,
754  selection_tuples=None,
755  resolution=1,
756  strength=0.01,
757  dist_cutoff=10.0,
758  ca_only=True,
759  hierarchy=None):
760  """Constructor
761  @param representation Representation object
762  @param selection_tuples Selecting regions for the restraint
763  @param resolution Resolution for applying restraint
764  @param strength Bond strength
765  @param dist_cutoff Cutoff for making restraints
766  @param ca_only Selects only CAlphas. Only matters if resolution=0.
767  @param hierarchy Root hierarchy to select from, use this instead of representation
768  """
769 
770  particles = []
771  if representation is None and hierarchy is not None:
772  self.m = hierarchy.get_model()
773  for st in selection_tuples:
774  if not ca_only:
775  sel = IMP.atom.Selection(hierarchy,chain=st[0],residue_indexes=range(st[1],st[2]+1))
776  else:
777  sel = IMP.atom.Selection(hierarchy,chain=st[0],residue_indexes=range(st[1],st[2]+1),
778  atom_type=IMP.atom.AtomType("CA"))
779  particles+=sel.get_selected_particles()
780  elif representation is not None and type(representation)==IMP.pmi.representation.Representation:
781  self.m = representation.mdl
782  for st in selection_tuples:
783  print('selecting with',st)
784  for p in IMP.pmi.tools.select_by_tuple(representation,st,resolution=resolution):
785  if (resolution==0 and ca_only and IMP.atom.Atom(p).get_atom_type()!=IMP.atom.AtomType("CA")):
786  continue
787  else:
788  particles.append(p.get_particle())
789  else:
790  raise Exception("must pass representation or hierarchy")
791 
792  self.weight = 1
793  self.label = "None"
794  self.pairslist = []
795 
796  # create score
797  self.rs = IMP.pmi.create_elastic_network(particles,dist_cutoff,strength)
798  for r in self.rs.get_restraints():
799  a1,a2 = r.get_inputs()
800  self.pairslist.append(IMP.ParticlePair(a1,a2))
801  self.pairslist.append(IMP.ParticlePair(a2,a1))
802  print('created',self.rs.get_number_of_restraints(),'restraints')
803 
804  def set_label(self, label):
805  self.label = label
806  self.rs.set_name(label)
807  for r in self.rs.get_restraints():
808  r.set_name(label)
809 
810  def add_to_model(self):
811  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
812 
813  def get_restraint(self):
814  return self.rs
815 
816  def set_weight(self, weight):
817  self.weight = weight
818  self.rs.set_weight(weight)
819 
820  def get_excluded_pairs(self):
821  return self.pairslist
822 
823  def get_output(self):
824  self.m.update()
825  output = {}
826  score = self.weight * self.rs.unprotected_evaluate(None)
827  output["_TotalScore"] = str(score)
828  output["ElasticNetworkRestraint_" + self.label] = str(score)
829  return output
830 
831 
833  """ Enable CHARMM force field """
834  def __init__(self,
835  root=None,
836  ff_temp=300.0,
837  zone_ps=None,
838  zone_size=10.0,
839  enable_nonbonded=True,
840  enable_bonded=True,
841  zone_nonbonded=False,
842  representation=None):
843  """Setup the CHARMM restraint on a selection. Expecting atoms.
844  @param root The node at which to apply the restraint
845  @param ff_temp The temperature of the force field
846  @param zone_ps Create a zone around this set of particles
847  Automatically includes the entire residue (incl. backbone)
848  @param zone_size The size for looking for neighbor residues
849  @param enable_nonbonded Allow the repulsive restraint
850  @param enable_bonded Allow the bonded restraint
851  @param zone_nonbonded EXPERIMENTAL: exclude from nonbonded all sidechains that aren't in zone!
852  @param representation Legacy representation object
853  """
854 
855  kB = (1.381 * 6.02214) / 4184.0
856  if representation is not None:
857  root = representation.prot
858 
859  self.mdl = root.get_model()
860  self.bonds_rs = IMP.RestraintSet(self.mdl, 1.0 / (kB * ff_temp), 'BONDED')
861  self.nonbonded_rs = IMP.RestraintSet(self.mdl, 1.0 / (kB * ff_temp), 'NONBONDED')
862  self.weight = 1.0
863  self.label = ""
864 
865  ### setup topology and bonds etc
866  if enable_bonded:
868  topology = ff.create_topology(root)
869  topology.apply_default_patches()
870  topology.setup_hierarchy(root)
871  if zone_ps is not None:
872  limit_to_ps=IMP.pmi.topology.get_particles_within_zone(
873  root,
874  zone_ps,
875  zone_size,
876  entire_residues=True,
877  exclude_backbone=False)
879  topology,
880  limit_to_ps)
881  self.ps = limit_to_ps
882  else:
883  r = IMP.atom.CHARMMStereochemistryRestraint(root, topology)
884  self.ps = IMP.core.get_leaves(root)
885  self.bonds_rs.add_restraint(r)
886  ff.add_radii(root)
887  ff.add_well_depths(root)
888 
889  atoms = IMP.atom.get_by_type(root,IMP.atom.ATOM_TYPE)
890  ### non-bonded forces
891  if enable_nonbonded:
892  if (zone_ps is not None) and zone_nonbonded:
893  print('stereochemistry: zone_nonbonded is True')
894  # atoms list should only include backbone plus zone_ps!
895  backbone_types=['C','N','CB','O']
896  sel = IMP.atom.Selection(root,atom_types=[IMP.atom.AtomType(n)
897  for n in backbone_types])
898  backbone_atoms = sel.get_selected_particles()
899  #x = list(set(backbone_atoms)|set(limit_to_ps))
900  sel_ps=IMP.pmi.topology.get_particles_within_zone(
901  root,
902  zone_ps,
903  zone_size,
904  entire_residues=True,
905  exclude_backbone=True)
906 
908  IMP.container.ListSingletonContainer(backbone_atoms),
910  4.0)
911  else:
912  cont = IMP.container.ListSingletonContainer(self.mdl,atoms)
913  self.nbl = IMP.container.ClosePairContainer(cont, 4.0)
914  if enable_bonded:
915  self.nbl.add_pair_filter(r.get_full_pair_filter())
916  pairscore = IMP.isd.RepulsiveDistancePairScore(0,1)
917  pr=IMP.container.PairsRestraint(pairscore, self.nbl)
918  self.nonbonded_rs.add_restraint(pr)
919  print('CHARMM is set up')
920 
921  def set_label(self, label):
922  self.label = label
923  self.rs.set_name(label)
924  for r in self.rs.get_restraints():
925  r.set_name(label)
926 
927  def add_to_model(self):
928  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.bonds_rs)
929  IMP.pmi.tools.add_restraint_to_model(self.mdl, self.nonbonded_rs)
930 
931  def get_restraint(self):
932  return self.rs
933 
934  def get_close_pair_container(self):
935  return self.nbl
936 
937  def set_weight(self, weight):
938  self.weight = weight
939  self.rs.set_weight(weight)
940 
941  def get_output(self):
942  self.mdl.update()
943  output = {}
944  bonds_score = self.weight * self.bonds_rs.unprotected_evaluate(None)
945  nonbonded_score = self.weight * self.nonbonded_rs.unprotected_evaluate(None)
946  score=bonds_score+nonbonded_score
947  output["_TotalScore"] = str(score)
948  output["CHARMM_BONDS"] = str(bonds_score)
949  output["CHARMM_NONBONDED"] = str(nonbonded_score)
950  return output
951 
952 
953 class PseudoAtomicRestraint(object):
954  """Add bonds and improper dihedral restraints for the CBs
955  """
956  def __init__(
957  self, rnums, representation, selection_tuple, strength=10.0, kappa=1.0,
958  jitter_angle=0.0, jitter_improper=0.0):
959  '''
960  need to add:
961  ca-ca bond
962  ca-cb is a constraint, no restraint needed
963  ca-ca-ca
964  cb-ca-ca-cb
965  '''
966 
967  self.m = representation.prot.get_model()
968  self.rs = IMP.RestraintSet(self.m, "PseudoAtomic")
969  self.rset_angles = IMP.RestraintSet(self.m, "PseudoAtomic_Angles")
970  self.rset_bonds = IMP.RestraintSet(self.m, "PseudoAtomic_Bonds")
971  self.weight = 1
972  self.label = "None"
973  self.pairslist = []
974 
975  # residues=IMP.pmi.tools.select_by_tuple(representation,selection_tuple,resolution=1)
976  for rnum in rnums:
977  ca, cb = self.get_ca_cb(
978  IMP.pmi.tools.select_by_tuple(representation,
979  (rnum, rnum, 'chainA'), resolution=0))
980  if not cb is None:
981  nter = False
982  cter = False
983  ca_prev, cb_prev = self.get_ca_cb(
984  IMP.pmi.tools.select_by_tuple(representation,
985  (rnum - 1, rnum - 1, 'chainA'), resolution=0))
986  ca_next, cb_next = self.get_ca_cb(
987  IMP.pmi.tools.select_by_tuple(representation,
988  (rnum + 1, rnum + 1, 'chainA'), resolution=0))
989  if ca_prev is None:
990  nter = True
991  else:
992  if ca_next is None:
993  cter = True
994  else:
995  if (nter and cter):
996  continue
997 
998  # adding a bond restraint between CA and CB
999  # h=IMP.core.Harmonic(6.0,kappa)
1000  # dps=IMP.core.DistancePairScore(h)
1001  # pr=IMP.core.PairRestraint(dps,IMP.ParticlePair(ca,cb))
1002  # self.pairslist.append((ca,cb))
1003  # self.rset_bonds.add_restraint(pr)
1004 
1005  # creating improper dihedral restraint
1006  # hus=IMP.core.Harmonic(2.09,kappa)
1008  2.09 +
1009  jitter_angle /
1010  0.5,
1011  kappa)
1013  2.09 -
1014  jitter_angle /
1015  0.5,
1016  kappa)
1017  if not nter:
1018  # ar13=IMP.core.AngleRestraint(hus,ca_prev,ca,cb)
1019  # self.rset_angles.add_restraint(ar13)
1020  ar13u = IMP.core.AngleRestraint(hupp, ca_prev, ca, cb)
1021  ar13l = IMP.core.AngleRestraint(hlow, ca_prev, ca, cb)
1022  self.rset_angles.add_restraint(ar13u)
1023  self.rset_angles.add_restraint(ar13l)
1024  if not cter:
1025  # ar23=IMP.core.AngleRestraint(hus,ca_next,ca,cb)
1026  # self.rset_angles.add_restraint(ar23)
1027  ar23u = IMP.core.AngleRestraint(hupp, ca_next, ca, cb)
1028  ar23l = IMP.core.AngleRestraint(hlow, ca_next, ca, cb)
1029  self.rset_angles.add_restraint(ar23u)
1030  self.rset_angles.add_restraint(ar23l)
1031  if not nter and not cter:
1032  # hus2=IMP.core.Harmonic(0,kappa)
1033  # idr=IMP.core.DihedralRestraint(hus2,ca,ca_prev,ca_next,cb)
1034  # self.rset_angles.add_restraint(idr)
1035 
1036  hus2upp = IMP.core.HarmonicUpperBound(
1037  jitter_improper,
1038  kappa)
1039  hus2low = IMP.core.HarmonicLowerBound(
1040  -
1041  jitter_improper,
1042  kappa)
1044  hus2upp,
1045  ca,
1046  ca_prev,
1047  ca_next,
1048  cb)
1050  hus2low,
1051  ca,
1052  ca_prev,
1053  ca_next,
1054  cb)
1055  self.rset_angles.add_restraint(idru)
1056  self.rset_angles.add_restraint(idrl)
1057  self.rs.add_restraint(self.rset_bonds)
1058  self.rs.add_restraint(self.rset_angles)
1059 
1060  def get_ca_cb(self, atoms):
1061  ca = None
1062  cb = None
1063  for a in atoms:
1064  if IMP.atom.Atom(a).get_atom_type() == IMP.atom.AtomType("CA"):
1065  ca = a.get_particle()
1066  elif IMP.atom.Atom(a).get_atom_type() == IMP.atom.AtomType("CB"):
1067  cb = a.get_particle()
1068  return ca, cb
1069 
1070  def set_label(self, label):
1071  self.label = label
1072  self.rs.set_name(label)
1073  for r in self.rs.get_restraints():
1074  r.set_name(label)
1075 
1076  def add_to_model(self):
1077  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1078 
1079  def get_restraint(self):
1080  return self.rs
1081 
1082  def set_weight(self, weight):
1083  self.weight = weight
1084  self.rs.set_weight(weight)
1085 
1086  def get_excluded_pairs(self):
1087  return self.pairslist
1088 
1089  def get_output(self):
1090  self.m.update()
1091  output = {}
1092  score = self.weight * self.rs.unprotected_evaluate(None)
1093  output["_TotalScore"] = str(score)
1094  output["PseudoAtomicRestraint_" + self.label] = str(score)
1095  return output
1096 
1097 class SymmetryRestraint(object):
1098  """Create harmonic restraints between the reference and (transformed) clones.
1099  \note Wraps IMP::core::TransformedDistancePairScore with an IMP::core::Harmonic
1100  """
1101  def __init__(self,
1102  references,
1103  clones_list,
1104  transforms,
1105  label='',
1106  strength=10.0):
1107  """Constructor
1108  @param references List of particles for symmetry reference
1109  @param clones_list List of lists of clone particles
1110  @param transforms Transforms moving each selection to the first selection
1111  @param label Label for output
1112  @param strength The elastic bond strength
1113  \note You will have to perform an IMP::atom::Selection to get the particles you need.
1114  Will check to make sure the numbers match.
1115  """
1116  self.mdl = root.get_model()
1117  self.rs = IMP.RestraintSet(self.mdl, "Symmetry")
1118  self.weight = 1
1119  self.label = label
1120  if len(clones_list)!=len(transforms):
1121  raise Exception('Error: There should be as many clones as transforms')
1122 
1123  harmonic = IMP.core.Harmonic(0.,strength)
1124  for clones,trans in zip(clones_list,transforms):
1125  if len(clones)!=len(references):
1126  raise Exception("Error: len(references)!=len(clones)")
1127  pair_score = IMP.core.TransformedDistancePairScore(harmonic,trans)
1128  for p0,p1 in zip(references,clones):
1129  r = IMP.core.PairRestraint(pair_score,(p0,p1))
1130  self.rs.add_restraint(r)
1131 
1132  print('created symmetry network with',self.rs.get_number_of_restraints(),'restraints')
1133 
1134  def set_label(self, label):
1135  self.label = label
1136  self.rs.set_name(label)
1137  for r in self.rs.get_restraints():
1138  r.set_name(label)
1139 
1140  def add_to_model(self):
1141  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1142 
1143  def get_restraint(self):
1144  return self.rs
1145 
1146  def set_weight(self, weight):
1147  self.weight = weight
1148  self.rs.set_weight(weight)
1149 
1150  def get_excluded_pairs(self):
1151  return self.pairslist
1152 
1153  def get_output(self):
1154  self.mdl.update()
1155  output = {}
1156  score = self.weight * self.rs.unprotected_evaluate(None)
1157  output["SymmetryRestraint_" + self.label] = str(score)
1158  output["_TotalScore"] = str(score)
1159  return output
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:358
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:359
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
Miscellaneous utilities.
Definition: tools.py:1
Add harmonic restraints between all pairs.
Apply a function to the distance between two particles after transforming the first.
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:469
Return all close ordered pairs of particles taken from the two SingletonContainers.
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:1455
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:23
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 add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:37
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.
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:702
Experimental, requires isd_emxl for now.
Set up the representation of all proteins and nucleic acid macromolecules.
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:1000
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:1096
Harmonic function (symmetric about the mean)
Definition: core/Harmonic.h:24