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