IMP logo
IMP Reference Guide  2.14.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 
14 
15 class ExternalBarrier(IMP.pmi.restraints.RestraintBase):
16 
17  """Restraint to keep all structures inside 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 radius Size of external barrier
23  @param hierarchies Can be one of the following inputs: IMP Hierarchy,
24  PMI System/State/Molecule/TempResidue, or a list/set of them
25  @param resolution Select which resolutions to act upon
26  @param weight Weight of restraint
27  @param center Center of the external barrier restraint
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 
61  """A simple distance restraint"""
62 
63  def __init__(self, root_hier, tuple_selection1, tuple_selection2,
64  distancemin=0, distancemax=100, resolution=1.0, kappa=1.0,
65  label=None, weight=1.):
66  """Setup distance restraint.
67  @param root_hier The hierarchy to select from
68  @param tuple_selection1 (resnum, resnum, molecule name, copy
69  number (=0))
70  @param tuple_selection2 (resnum, resnum, molecule name, copy
71  number (=0))
72  @param distancemin The minimum dist
73  @param distancemax The maximum dist
74  @param resolution For selecting particles
75  @param kappa The harmonic parameter
76  @param label A unique label to be used in outputs and
77  particle/restraint names
78  @param weight Weight of restraint
79  @note Pass the same resnum twice to each tuple_selection. Optionally
80  add a copy number (PMI2 only)
81  """
82  ts1 = IMP.core.HarmonicUpperBound(distancemax, kappa)
83  ts2 = IMP.core.HarmonicLowerBound(distancemin, kappa)
84 
85  model = root_hier.get_model()
86  copy_num1 = 0
87  if len(tuple_selection1) > 3:
88  copy_num1 = tuple_selection1[3]
89  copy_num2 = 0
90  if len(tuple_selection2) > 3:
91  copy_num2 = tuple_selection2[3]
92 
93  sel1 = IMP.atom.Selection(root_hier,
94  resolution=resolution,
95  molecule=tuple_selection1[2],
96  residue_index=tuple_selection1[0],
97  copy_index=copy_num1)
98  particles1 = sel1.get_selected_particles()
99  sel2 = IMP.atom.Selection(root_hier,
100  resolution=resolution,
101  molecule=tuple_selection2[2],
102  residue_index=tuple_selection2[0],
103  copy_index=copy_num2)
104  particles2 = sel2.get_selected_particles()
105 
106  super(DistanceRestraint, self).__init__(model, label=label,
107  weight=weight)
108  print(self.name)
109 
110  print("Created distance restraint between "
111  "%s and %s" % (particles1[0].get_name(),
112  particles2[0].get_name()))
113 
114  if len(particles1) > 1 or len(particles2) > 1:
115  raise ValueError("more than one particle selected")
116 
117  self.rs.add_restraint(
118  IMP.core.DistanceRestraint(self.model, ts1,
119  particles1[0],
120  particles2[0]))
121  self.rs.add_restraint(
122  IMP.core.DistanceRestraint(self.model, ts2,
123  particles1[0],
124  particles2[0]))
125 
126 
128  '''
129  PMI2 python restraint. Restrains particles within a
130  Cylinder aligned along the z-axis and
131  centered in x,y=0,0
132  Optionally, one can restrain the cylindrical angle
133  '''
134  import math
135  def __init__(self, m, objects, resolution, radius,mintheta=None,
136  maxtheta=None,repulsive=False,label='None'):
137  '''
138  @param objects PMI2 objects
139  @param resolution the resolution you want the restraint to be applied
140  @param radius the radius of the cylinder
141  @param mintheta minimum cylindrical angle in degrees
142  @param maxtheta maximum cylindrical angle in degrees
143  '''
144  IMP.Restraint.__init__(self, m, "CylinderRestraint %1%")
145  self.radius=radius
146  self.softness = 3.0
147  self.softness_angle = 0.5
148  self.plateau = 1e-10
149  self.weight = 1.0
150  self.m=m
151  self.mintheta=mintheta
152  self.maxtheta=maxtheta
153  self.repulsive=repulsive
154  hierarchies = IMP.pmi.tools.input_adaptor(objects,
155  resolution,
156  flatten=True)
157  self.particles = [h.get_particle() for h in hierarchies]
158  self.label=label
159 
160  def get_probability(self,p):
161  xyz=IMP.core.XYZ(p)
162  r=self.math.sqrt(xyz.get_x()**2+xyz.get_y()**2)
163  argvalue=(r-self.radius) / self.softness
164  if self.repulsive: 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) / (1.0 + self.math.exp(-max(argvalue1,argvalue2)))
176  return prob
177 
178  def unprotected_evaluate(self, da):
179  s=0.0
180  for p in self.particles:
181  s+=-self.math.log(1.0-self.get_probability(p))
182  if self.mintheta is not None and self.maxtheta is not None:
183  s+=-self.math.log(1.0-self.get_angle_probability(p))
184  return s
185 
186  def do_get_inputs(self):
187  return self.particles
188 
189  def add_to_model(self):
191 
192  def get_output(self):
193  output = {}
194  score = self.weight * self.unprotected_evaluate(None)
195  output["_TotalScore"] = str(score)
196  output["CylinderRestraint_" + self.label] = str(score)
197  return output
198 
199 
200 
202  '''
203  a python restraint with bistable potential
204  Authors: G. Bouvier, R. Pellarin. Pasteur Institute.
205  '''
206  import numpy as np
207  import math
208 
209  def __init__(self,m,p1,p2,dist1,dist2,sigma1,sigma2,weight1,weight2):
210  '''
211  input twp particles, the two equilibrium distances, their amplitudes, and their weights (populations)
212  '''
213  IMP.Restraint.__init__(self, m, "BiStableDistanceRestraint %1%")
214  self.dist1=dist1
215  self.dist2=dist2
216 
217  self.sigma1=sigma1
218  self.sigma2=sigma2
219 
220  self.weight1=weight1
221  self.weight2=weight2
222 
223  if self.weight1+self.weight2 != 1:
224  raise ValueError("The sum of the weights must be one")
225 
226  self.d1=IMP.core.XYZ(p1)
227  self.d2=IMP.core.XYZ(p2)
228  self.particle_list=[p1,p2]
229 
230  def gaussian(self,x, mu, sig, w):
231  return w*self.np.exp(-self.np.power(x - mu, 2.) / (2 * self.np.power(sig, 2.)))
232 
233  def unprotected_evaluate(self,da):
234  dist=IMP.core.get_distance(self.d1,self.d2)
235  prob=self.gaussian(dist,self.dist1,self.sigma1,self.weight1)+\
236  self.gaussian(dist,self.dist2,self.sigma2,self.weight2)
237  return -self.math.log(prob)
238 
239  def do_get_inputs(self):
240  return self.particle_list
241 
242 
243 class DistanceToPointRestraint(IMP.pmi.restraints.RestraintBase):
244 
245  """Restraint for anchoring a particle to a specific coordinate."""
246 
247  def __init__(self, root_hier, tuple_selection,
248  anchor_point=IMP.algebra.Vector3D(0, 0, 0),
249  radius=10.0, kappa=10.0, resolution=1.0, weight=1.0,
250  label=None):
251  """Setup distance restraint.
252  @param root_hier The hierarchy to select from
253  @param tuple_selection (resnum, resnum, molecule name,
254  copy number (=0))
255  @param anchor_point Point to which to restrain particle
256  (IMP.algebra.Vector3D object)
257  @param radius Size of the tolerance length
258  @param kappa The harmonic parameter
259  @param resolution For selecting a particle
260  @param weight Weight of restraint
261  @param label A unique label to be used in outputs and
262  particle/restraint names
263  @note Pass the same resnum twice to each tuple_selection. Optionally
264  add a copy number (PMI2 only)
265  """
266  model = root_hier.get_model()
267  copy_num1 = 0
268  if len(tuple_selection) > 3:
269  copy_num1 = tuple_selection[3]
270 
271  sel1 = IMP.atom.Selection(root_hier,
272  resolution=resolution,
273  molecule=tuple_selection[2],
274  residue_index=tuple_selection[0],
275  copy_index=copy_num1)
276  ps = sel1.get_selected_particles()
277  if len(ps) > 1:
278  raise ValueError("More than one particle selected")
279 
280  super(DistanceToPointRestraint, self).__init__(model, label=label,
281  weight=weight)
282  self.radius = radius
283 
284  ub3 = IMP.core.HarmonicUpperBound(self.radius, kappa)
285  if anchor_point is None:
286  c3 = IMP.algebra.Vector3D(0, 0, 0)
287  elif isinstance(anchor_point, IMP.algebra.Vector3D):
288  c3 = anchor_point
289  else:
290  raise TypeError("anchor_point must be an algebra.Vector3D object")
291  ss3 = IMP.core.DistanceToSingletonScore(ub3, c3)
292 
293  lsc = IMP.container.ListSingletonContainer(self.model)
294  lsc.add(ps)
295 
296  r3 = IMP.container.SingletonsRestraint(ss3, lsc)
297  self.rs.add_restraint(r3)
298 
299  print("\n%s: Created distance_to_point_restraint between "
300  "%s and %s" % (self.name, ps[0].get_name(), c3))
301 
302 class MembraneRestraint(IMP.pmi.restraints.RestraintBase):
303  def __init__(self,
304  hier,
305  objects_above=None,
306  objects_inside=None,
307  objects_below=None,
308  center = 0.0,
309  thickness=30.0,
310  softness=3.0,
311  plateau=0.0000000001,
312  resolution=1,
313  weight = 1.0,
314  label = None):
315  """ Setup Membrane restraint
316 
317  Simple sigmoid score calculated for particles above,
318 
319  @param objects_inside list or tuples of objects in membrane (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)) at the plateau
325  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  input a list of particles, the slope and theta of the sigmoid potential
330  theta is the cutoff distance for a protein-protein contact
331  """
332 
333  self.hier = hier
334  model = self.hier.get_model()
335 
336  rname = "MembraneRestraint"
337  super(MembraneRestraint, self).__init__(
338  model, name="MembraneRestraint", label=label, weight=weight)
339 
340  self.center = center
341  self.thickness = thickness
342  self.softness = softness
343  self.plateau = plateau
344  self.linear = 0.02
345  self.resolution = resolution
346 
347  # Create nuisance particle
348  p = IMP.Particle(model)
349  z_center = IMP.isd.Nuisance.setup_particle(p)
350  z_center.set_nuisance(self.center)
351 
352  # Setup restraint
353  mr = IMP.pmi.MembraneRestraint(model,
354  z_center.get_particle_index(),
355  self.thickness,
356  self.softness,
357  self.plateau,
358  self.linear)
359 
360  # Particles above
361  if objects_above:
362  for obj in objects_above:
363  if isinstance(obj, tuple):
364  self.particles_above = self._select_from_tuple(obj)
365 
366  elif isinstance(obj, str):
367  self.particles_above = self._select_from_string(obj)
368  mr.add_particles_above(self.particles_above)
369 
370  # Particles inside
371  if objects_inside:
372  for obj in objects_inside:
373  if isinstance(obj, tuple):
374  self.particles_inside = self._select_from_tuple(obj)
375 
376  elif isinstance(obj, str):
377  self.particles_inside = self._select_from_string(obj)
378  mr.add_particles_inside(self.particles_inside)
379 
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 
392  self.rs.add_restraint(mr)
393 
394  def get_particles_above(self):
395  return self.particles_above
396 
397  def get_particles_inside(self):
398  return self.particles_inside
399 
400  def get_particles_below(self):
401  return self.particles_below
402 
403  def _select_from_tuple(self, obj):
404  particles = IMP.atom.Selection(self.hier,
405  molecule = obj[2],
406  residue_indexes = range(obj[0], obj[1]+1, 1),
407  resolution = self.resolution).get_selected_particles()
408 
409  return particles
410 
411  def _select_from_string(self, obj):
412  particles = IMP.atom.Selection(self.hier,
413  molecule = obj,
414  resolution = self.resolution).get_selected_particles()
415  return particles
416 
417  def create_membrane_density(self, file_out='membrane_localization.mrc'):
418 
419  '''
420  Just for visualization purposes.
421  Writes density of membrane localization
422  '''
423 
424  offset = 5.0*self.thickness
425  apix = 3.0
426  resolution = 5.0
427 
428  # Create a density header of the requested size
430  IMP.algebra.Vector3D(-self.center - offset, -self.center - offset, -self.center - offset,),
431  IMP.algebra.Vector3D(self.center + offset, self.center + offset, self.center + offset))
432  dheader = IMP.em.create_density_header(bbox, apix)
433  dheader.set_resolution(resolution)
434  dmap = IMP.em.SampledDensityMap(dheader)
435 
436  for vox in range(dmap.get_header().get_number_of_voxels()):
437  c = dmap.get_location_by_voxel(vox)
438  if self._is_membrane(c[2])==1:
439  dmap.set_value(c[0], c[1], c[2], 1.0)
440  else:
441  dmap.set_value(c[0], c[1], c[2], 0.0)
442 
443  IMP.em.write_map(dmap, file_out)
444 
445  def _is_membrane(self, z):
446  if (z-self.center) < self.thickness/2.0 and (z-self.center) >= -self.thickness/2.0 :
447  return 1
448  else:
449  return 0
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
Distance restraint between two particles.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
Definition: XYZR.h:89
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:917
def __init__
input twp particles, the two equilibrium distances, their amplitudes, and their weights (populations)...
Class for sampling a density map from particles.
def __init__
Setup external barrier restraint.
Store a list of ParticleIndexes.
Restraint for anchoring a particle to a specific coordinate.
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
Restraint to keep all structures inside 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...
The general base class for IMP exceptions.
Definition: exception.h:49
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.
a python restraint with bistable potential Authors: G.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
virtual ModelObjectsTemp do_get_inputs() const =0
A restraint is a term in an IMP ScoringFunction.
Definition: Restraint.h:54