IMP logo
IMP Reference Guide  2.18.0
The Integrative Modeling Platform
restraints/basic.py
1 """@namespace IMP.pmi.restraints.basic
2 Some miscellaneous simple restraints.
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.pmi.tools
12 import IMP.pmi.restraints
13 import math
14 
15 
16 class ExternalBarrier(IMP.pmi.restraints.RestraintBase):
17  """Keeps all structures inside a sphere."""
18 
19  def __init__(self, hierarchies, radius=10.0, resolution=10, weight=1.0,
20  center=None, label=None):
21  """Setup external barrier restraint.
22  @param hierarchies Can be one of the following inputs: IMP Hierarchy,
23  PMI System/State/Molecule/TempResidue, or a list/set of them
24  @param radius Size of external barrier
25  @param resolution Select which resolutions to act upon
26  @param weight Weight of restraint
27  @param center Center of the external barrier
28  (IMP.algebra.Vector3D object)
29  @param label A unique label to be used in outputs and
30  particle/restraint names.
31  """
32  hiers = IMP.pmi.tools.input_adaptor(hierarchies, resolution,
33  flatten=True)
34  model = hiers[0].get_model()
35  particles = [h.get_particle() for h in hiers]
36 
37  super(ExternalBarrier, self).__init__(model, label=label,
38  weight=weight)
39  self.radius = radius
40 
41  if center is None:
42  c3 = IMP.algebra.Vector3D(0, 0, 0)
43  elif type(center) is IMP.algebra.Vector3D:
44  c3 = center
45  else:
46  raise Exception(
47  "%s: @param center must be an IMP.algebra.Vector3D object" % (
48  self.name))
49 
50  ub3 = IMP.core.HarmonicUpperBound(radius, 10.0)
52  lsc = IMP.container.ListSingletonContainer(self.model)
53 
54  lsc.add(particles)
56  self.rs.add_restraint(r3)
57 
58 
59 class DistanceRestraint(IMP.pmi.restraints.RestraintBase):
60  """A simple distance restraint"""
61 
62  def __init__(self, root_hier, tuple_selection1, tuple_selection2,
63  distancemin=0, distancemax=100, resolution=1.0, kappa=1.0,
64  label=None, weight=1.):
65  """Setup distance restraint.
66  @param root_hier The hierarchy to select from
67  @param tuple_selection1 (resnum, resnum, molecule name, copy
68  number (=0))
69  @param tuple_selection2 (resnum, resnum, molecule name, copy
70  number (=0))
71  @param distancemin The minimum dist
72  @param distancemax The maximum dist
73  @param resolution For selecting particles
74  @param kappa The harmonic parameter
75  @param label A unique label to be used in outputs and
76  particle/restraint names
77  @param weight Weight of restraint
78  @note Pass the same resnum twice to each tuple_selection. Optionally
79  add a copy number.
80  """
81  ts1 = IMP.core.HarmonicUpperBound(distancemax, kappa)
82  ts2 = IMP.core.HarmonicLowerBound(distancemin, kappa)
83 
84  model = root_hier.get_model()
85  copy_num1 = 0
86  if len(tuple_selection1) > 3:
87  copy_num1 = tuple_selection1[3]
88  copy_num2 = 0
89  if len(tuple_selection2) > 3:
90  copy_num2 = tuple_selection2[3]
91 
92  sel1 = IMP.atom.Selection(root_hier,
93  resolution=resolution,
94  molecule=tuple_selection1[2],
95  residue_index=tuple_selection1[0],
96  copy_index=copy_num1)
97  particles1 = sel1.get_selected_particles()
98  sel2 = IMP.atom.Selection(root_hier,
99  resolution=resolution,
100  molecule=tuple_selection2[2],
101  residue_index=tuple_selection2[0],
102  copy_index=copy_num2)
103  particles2 = sel2.get_selected_particles()
104 
105  super(DistanceRestraint, self).__init__(model, label=label,
106  weight=weight)
107  print(self.name)
108 
109  print("Created distance restraint between "
110  "%s and %s" % (particles1[0].get_name(),
111  particles2[0].get_name()))
112 
113  if len(particles1) > 1 or len(particles2) > 1:
114  raise ValueError("more than one particle selected")
115 
116  self.rs.add_restraint(
117  IMP.core.DistanceRestraint(self.model, ts1,
118  particles1[0],
119  particles2[0]))
120  self.rs.add_restraint(
121  IMP.core.DistanceRestraint(self.model, ts2,
122  particles1[0],
123  particles2[0]))
124 
125 
127  """Restrain particles within (or outside) a cylinder.
128  The cylinder is aligned along the z-axis and with center x=y=0.
129  Optionally, one can restrain the cylindrical angle
130  """
131  import math
132 
133  def __init__(self, m, objects, resolution, radius, mintheta=None,
134  maxtheta=None, repulsive=False, label='None'):
135  '''
136  @param objects PMI2 objects to restrain
137  @param resolution the resolution you want the restraint to be applied
138  @param radius the radius of the cylinder
139  @param mintheta minimum cylindrical angle in degrees
140  @param maxtheta maximum cylindrical angle in degrees
141  @param repulsive If True, restrain the particles to be outside
142  of the cylinder instead of inside
143  @param label A unique label to be used in outputs and
144  particle/restraint names
145  '''
146  IMP.Restraint.__init__(self, m, "CylinderRestraint %1%")
147  self.radius = radius
148  self.softness = 3.0
149  self.softness_angle = 0.5
150  self.plateau = 1e-10
151  self.weight = 1.0
152  self.m = m
153  self.mintheta = mintheta
154  self.maxtheta = maxtheta
155  self.repulsive = repulsive
156  hierarchies = IMP.pmi.tools.input_adaptor(objects,
157  resolution,
158  flatten=True)
159  self.particles = [h.get_particle() for h in hierarchies]
160  self.label = label
161 
162  def get_probability(self, p):
163  xyz = IMP.core.XYZ(p)
164  r = self.math.sqrt(xyz.get_x()**2+xyz.get_y()**2)
165  argvalue = (r-self.radius) / self.softness
166  if self.repulsive:
167  argvalue = -argvalue
168  prob = (1.0 - self.plateau) / (1.0 + self.math.exp(-argvalue))
169  return prob
170 
171  def get_angle_probability(self, p):
172  xyz = IMP.core.XYZ(p)
173  angle = self.math.atan2(xyz.get_y(), xyz.get_x())*180.0/self.math.pi
174  anglediff = (angle - self.maxtheta + 180 + 360) % 360 - 180
175  argvalue1 = anglediff / self.softness_angle
176  anglediff = (angle - self.mintheta + 180 + 360) % 360 - 180
177  argvalue2 = -anglediff / self.softness_angle
178  prob = ((1.0-self.plateau)
179  / (1.0 + self.math.exp(-max(argvalue1, argvalue2))))
180  return prob
181 
182  def unprotected_evaluate(self, da):
183  s = 0.0
184  for p in self.particles:
185  s += -self.math.log(1.0-self.get_probability(p))
186  if self.mintheta is not None and self.maxtheta is not None:
187  s += -self.math.log(1.0-self.get_angle_probability(p))
188  return s
189 
190  def do_get_inputs(self):
191  return self.particles
192 
193  def add_to_model(self):
195 
196  def get_output(self):
197  output = {}
198  score = self.weight * self.unprotected_evaluate(None)
199  output["_TotalScore"] = str(score)
200  output["CylinderRestraint_" + self.label] = str(score)
201  return output
202 
203 
205  '''Distance restraint with bistable potential
206  Authors: G. Bouvier, R. Pellarin. Pasteur Institute.
207  '''
208  import numpy as np
209  import math
210 
211  def __init__(self, m, p1, p2, dist1, dist2, sigma1, sigma2, weight1,
212  weight2):
213  '''
214  input two particles, the two equilibrium distances, their amplitudes,
215  and their weights (populations)
216  '''
217  IMP.Restraint.__init__(self, m, "BiStableDistanceRestraint %1%")
218  self.dist1 = dist1
219  self.dist2 = dist2
220 
221  self.sigma1 = sigma1
222  self.sigma2 = sigma2
223 
224  self.weight1 = weight1
225  self.weight2 = weight2
226 
227  if self.weight1+self.weight2 != 1:
228  raise ValueError("The sum of the weights must be one")
229 
230  self.d1 = IMP.core.XYZ(p1)
231  self.d2 = IMP.core.XYZ(p2)
232  self.particle_list = [p1, p2]
233 
234  def gaussian(self, x, mu, sig, w):
235  return (w*self.np.exp(-self.np.power(x - mu, 2.)
236  / (2 * self.np.power(sig, 2.))))
237 
238  def unprotected_evaluate(self, da):
239  dist = IMP.core.get_distance(self.d1, self.d2)
240  prob = self.gaussian(dist, self.dist1, self.sigma1, self.weight1) + \
241  self.gaussian(dist, self.dist2, self.sigma2, self.weight2)
242  return -self.math.log(prob)
243 
244  def do_get_inputs(self):
245  return self.particle_list
246 
247 
248 class DistanceToPointRestraint(IMP.pmi.restraints.RestraintBase):
249  """Anchor a particle to a specific coordinate."""
250 
251  def __init__(self, root_hier, tuple_selection,
252  anchor_point=IMP.algebra.Vector3D(0, 0, 0),
253  radius=10.0, kappa=10.0, resolution=1.0, weight=1.0,
254  label=None):
255  """Setup distance restraint.
256  @param root_hier The hierarchy to select from
257  @param tuple_selection (resnum, resnum, molecule name,
258  copy number (=0))
259  @param anchor_point Point to which to restrain particle
260  (IMP.algebra.Vector3D object)
261  @param radius Maximum distance the particle can move from the
262  fixed point before it is restrained
263  @param kappa Strength of the harmonic restraint for point-particle
264  distance greater than the radius
265  @param resolution For selecting a particle
266  @param weight Weight of restraint
267  @param label A unique label to be used in outputs and
268  particle/restraint names
269  @note Pass the same resnum twice to each tuple_selection. Optionally
270  add a copy number
271  """
272  model = root_hier.get_model()
273  copy_num1 = 0
274  if len(tuple_selection) > 3:
275  copy_num1 = tuple_selection[3]
276 
277  sel1 = IMP.atom.Selection(root_hier,
278  resolution=resolution,
279  molecule=tuple_selection[2],
280  residue_index=tuple_selection[0],
281  copy_index=copy_num1)
282  ps = sel1.get_selected_particles()
283  if len(ps) > 1:
284  raise ValueError("More than one particle selected")
285 
286  super(DistanceToPointRestraint, self).__init__(model, label=label,
287  weight=weight)
288  self.radius = radius
289 
290  ub3 = IMP.core.HarmonicUpperBound(self.radius, kappa)
291  if anchor_point is None:
292  c3 = IMP.algebra.Vector3D(0, 0, 0)
293  elif isinstance(anchor_point, IMP.algebra.Vector3D):
294  c3 = anchor_point
295  else:
296  raise TypeError("anchor_point must be an algebra.Vector3D object")
297  ss3 = IMP.core.DistanceToSingletonScore(ub3, c3)
298 
299  lsc = IMP.container.ListSingletonContainer(self.model)
300  lsc.add(ps)
301 
302  r3 = IMP.container.SingletonsRestraint(ss3, lsc)
303  self.rs.add_restraint(r3)
304 
305  print("\n%s: Created distance_to_point_restraint between "
306  "%s and %s" % (self.name, ps[0].get_name(), c3))
307 
308 
309 class MembraneRestraint(IMP.pmi.restraints.RestraintBase):
310  """Restrain particles to be above, below, or inside a planar membrane.
311  The membrane is defined to lie on the xy plane with a given z
312  coordinate and thickness, and particles are restrained (by their
313  z coordinates) with a simple sigmoid score.
314  """
315 
316  def __init__(self, hier, objects_above=None, objects_inside=None,
317  objects_below=None, center=0.0, thickness=30.0,
318  softness=3.0, plateau=0.0000000001, resolution=1,
319  weight=1.0, label=None):
320  """Setup the restraint.
321 
322  @param objects_inside list or tuples of objects in membrane
323  (e.g. ['p1', (10, 30,'p2')])
324  @param objects_above list or tuples of objects above membrane
325  @param objects_below list or tuples of objects below membrane
326  @param thickness Thickness of the membrane along the z-axis
327  @param softness Softness of the limiter in the sigmoid function
328  @param plateau Parameter to set the probability (=1- plateau))
329  at the plateau phase of the sigmoid
330  @param weight Weight of restraint
331  @param label A unique label to be used in outputs and
332  particle/restraint names.
333  """
334 
335  self.hier = hier
336  model = self.hier.get_model()
337 
338  super(MembraneRestraint, self).__init__(
339  model, name="MembraneRestraint", label=label, weight=weight)
340 
341  self.center = center
342  self.thickness = thickness
343  self.softness = softness
344  self.plateau = plateau
345  self.linear = 0.02
346  self.resolution = resolution
347 
348  # Create nuisance particle
349  p = IMP.Particle(model)
350  z_center = IMP.isd.Nuisance.setup_particle(p)
351  z_center.set_nuisance(self.center)
352 
353  # Setup restraint
354  mr = IMP.pmi.MembraneRestraint(model,
355  z_center.get_particle_index(),
356  self.thickness,
357  self.softness,
358  self.plateau,
359  self.linear)
360 
361  # Particles above
362  if objects_above:
363  for obj in objects_above:
364  if isinstance(obj, tuple):
365  self.particles_above = self._select_from_tuple(obj)
366 
367  elif isinstance(obj, str):
368  self.particles_above = self._select_from_string(obj)
369  mr.add_particles_above(self.particles_above)
370 
371  # Particles inside
372  if objects_inside:
373  for obj in objects_inside:
374  if isinstance(obj, tuple):
375  self.particles_inside = self._select_from_tuple(obj)
376 
377  elif isinstance(obj, str):
378  self.particles_inside = self._select_from_string(obj)
379  mr.add_particles_inside(self.particles_inside)
380 
381  # Particles below
382  if objects_below:
383  for obj in objects_below:
384  if isinstance(obj, tuple):
385  self.particles_below = self._select_from_tuple(obj)
386 
387  elif isinstance(obj, str):
388  self.particles_below = self._select_from_string(obj)
389  mr.add_particles_below(self.particles_below)
390 
391  self.rs.add_restraint(mr)
392 
393  def get_particles_above(self):
394  return self.particles_above
395 
396  def get_particles_inside(self):
397  return self.particles_inside
398 
399  def get_particles_below(self):
400  return self.particles_below
401 
402  def _select_from_tuple(self, obj):
403  particles = IMP.atom.Selection(
404  self.hier, molecule=obj[2],
405  residue_indexes=range(obj[0], obj[1]+1, 1),
406  resolution=self.resolution).get_selected_particles()
407 
408  return particles
409 
410  def _select_from_string(self, obj):
411  particles = IMP.atom.Selection(
412  self.hier, molecule=obj,
413  resolution=self.resolution).get_selected_particles()
414  return particles
415 
416  def create_membrane_density(self, file_out='membrane_localization.mrc'):
417  """Create an MRC density file to visualize the membrane."""
418  offset = 5.0 * self.thickness
419  apix = 3.0
420  resolution = 5.0
421 
422  # Create a density header of the requested size
424  IMP.algebra.Vector3D(-self.center - offset, -self.center - offset,
425  -self.center - offset),
426  IMP.algebra.Vector3D(self.center + offset, self.center + offset,
427  self.center + offset))
428  dheader = IMP.em.create_density_header(bbox, apix)
429  dheader.set_resolution(resolution)
430  dmap = IMP.em.SampledDensityMap(dheader)
431 
432  for vox in range(dmap.get_header().get_number_of_voxels()):
433  c = dmap.get_location_by_voxel(vox)
434  if self._is_membrane(c[2]) == 1:
435  dmap.set_value(c[0], c[1], c[2], 1.0)
436  else:
437  dmap.set_value(c[0], c[1], c[2], 0.0)
438 
439  IMP.em.write_map(dmap, file_out)
440 
441  def _is_membrane(self, z):
442  if ((z-self.center) < self.thickness/2.0 and
443  (z-self.center) >= -self.thickness/2.0):
444  return 1
445  else:
446  return 0
447 
448 
449 class ResidueProteinProximityRestraint(IMP.pmi.restraints.RestraintBase):
450  """Restrain residue/residues to bind to unknown location in a target"""
451 
452  def __init__(self, hier, selection, cutoff=6., sigma=3., xi=0.01,
453  resolution=1.0, weight=1.0, label=None):
454  """
455  Constructor
456  @param hier Hierarchy of the system
457  @param section Selection of residues and target;
458  syntax is (prot, r1, r2, target_prot) or
459  (prot1, r1, r2, target_prot, target_r1, target_r2)
460  @param cutoff Distance cutoff between selected segment and target
461  protein
462  @param sigma Distance variance between selected fragments
463  @param xi Slope of a distance-linear scoring function that
464  funnels the score when the particles are too
465  far away
466  @param resolution Resolution at which to apply restraint
467  @param weight Weight of the restraint
468  @param label Extra text to label the restraint so that it is
469  searchable in the output
470  """
471  self.hier = hier
472  m = self.hier.get_model()
473 
474  super(ResidueProteinProximityRestraint, self).__init__(
475  m, name="ResidueProteinProximityRestraint", label=label,
476  weight=weight)
477 
478  self.cutoff = cutoff
479  self.sigma = sigma
480  self.xi = xi
481  self.resolution = resolution
482 
483  # Check selection
484  print('selection', selection, isinstance(selection, tuple))
485  if (not isinstance(selection, tuple)
486  and not isinstance(selection, list)):
487  raise ValueError("Selection should be a tuple or list")
488  if len(selection) < 4:
489  raise ValueError(
490  "Selection should be (prot, r1, r2, target_prot) or "
491  "(prot1, r1, r2, target_prot, target_r1, target_r2)")
492 
493  # Selection
494  self.prot1 = selection[0]
495  self.r1 = int(selection[1])
496  self.r2 = int(selection[2])
497 
498  self.prot2 = selection[3]
499  if len(selection) == 6:
500  self.tr1 = int(selection[4])
501  self.tr2 = int(selection[5])
502 
503  if self.r1 == self.r2:
504  sel_resi = IMP.atom.Selection(
505  self.hier, molecule=self.prot1, residue_index=self.r1,
506  resolution=self.resolution).get_selected_particles()
507  else:
508  sel_resi = IMP.atom.Selection(
509  self.hier, molecule=self.prot1,
510  residue_indexes=range(self.r1, self.r2+1, 1),
511  resolution=self.resolution).get_selected_particles()
512 
513  if len(selection) == 4:
514  sel_target = IMP.atom.Selection(
515  self.hier, molecule=self.prot2,
516  resolution=self.resolution).get_selected_particles()
517 
518  elif len(selection) == 6:
519  sel_target = IMP.atom.Selection(
520  self.hier, molecule=self.prot2,
521  residue_indexes=range(self.tr1, self.tr2+1, 1),
522  resolution=self.resolution).get_selected_particles()
523 
524  self.included_ps = sel_resi + sel_target
525 
526  # Setup restraint
527  distance = 0.0
528  slack = cutoff*2
529 
531  m, self.cutoff, self.sigma, self.xi, True,
532  'ResidueProteinProximityRestraint')
533 
534  print('Selected fragment and target lengths:', len(sel_resi),
535  len(sel_target))
536 
537  # Setup close pair container
538  # Find close pair within included_resi and included_target
540  lsa_target.add(IMP.get_indexes(sel_target))
541 
543  lsa_resi.add(IMP.get_indexes(sel_resi))
544 
545  self.cpc = IMP.container.CloseBipartitePairContainer(lsa_resi,
546  lsa_target,
547  distance,
548  slack)
549 
550  br.add_pairs_container(self.cpc)
551 
552  br.add_contribution_particles(sel_resi, sel_target)
553 
554  # Compute interpolation paramaters
555  yi = ((cutoff**2/(2*sigma**2)
556  - math.log(1/math.sqrt(2*math.pi*sigma*sigma))+cutoff*xi/2.)
557  / (cutoff/2.))
558  interpolation_factor = -(cutoff/2.)*(xi-yi)
559  max_p = (math.exp(-((distance+slack)**2)/(2*sigma**2))
560  / math.sqrt(2*math.pi*sigma*sigma))
561  max_score = -math.log(max_p)
562 
563  # Add interpolation parameters
564  br.set_yi(yi)
565  br.set_interpolation_factor(interpolation_factor)
566  br.set_max_score(max_score)
567 
568  self.rs.add_restraint(br)
569 
570  self.restraint_sets = [self.rs] + self.restraint_sets[1:]
571 
573  """ Get particles in the close pair container """
574  return self.cpc.get_indexes()
575 
576  def get_output(self):
577  output = {}
578  score = self.weight * self.rs.unprotected_evaluate(None)
579  output["ResidueProteinProximityRestraint_score_" + self.label] \
580  = str(score)
581 
582  return output
Applies a SingletonScore to each Singleton in a list.
static Nuisance setup_particle(Model *m, ParticleIndex pi)
Definition: Nuisance.h:31
def __init__
Setup distance restraint.
Lower bound harmonic function (non-zero when feature < mean)
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
Miscellaneous utilities.
Definition: tools.py:1
virtual double unprotected_evaluate(DerivativeAccumulator *da) const
Return the unweighted score for the restraint.
Return all spatially-proximals pairs of particles (a,b) from the two SingletonContainers A and B...
Distance restraint between two particles.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
Definition: XYZR.h:89
Restrain residue/residues to bind to unknown location in a target.
Classes to handle different kinds of restraints.
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: tools.py:946
def __init__
input two particles, the two equilibrium distances, their amplitudes, and their weights (populations)...
Class for sampling a density map from particles.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
Get the indexes from a list of particle pairs.
def __init__
Setup external barrier restraint.
Store a list of ParticleIndexes.
Anchor a particle to a specific coordinate.
def create_membrane_density
Create an MRC density file to visualize the membrane.
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:91
DensityHeader create_density_header(const algebra::BoundingBoxD< 3 > &bb, float spacing)
Create a header from a bounding box in 3D.
Apply a function to the distance to a fixed point.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
Restrain particles to be above, below, or inside a planar membrane.
Keeps all structures inside a sphere.
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...
def get_container_pairs
Get particles in the close pair container.
The general base class for IMP exceptions.
Definition: exception.h:48
VectorD< 3 > Vector3D
Definition: VectorD.h:421
Class to handle individual particles of a Model object.
Definition: Particle.h:41
Functionality for loading, creating, manipulating and scoring atomic structures.
Distance restraint with bistable potential Authors: G.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:70
virtual ModelObjectsTemp do_get_inputs() const =0
Restrain particles within (or outside) a cylinder.
A restraint is a term in an IMP ScoringFunction.
Definition: Restraint.h:53