IMP logo
IMP Reference Guide  2.23.0
The Integrative Modeling Platform
version 2.23.0
domino/rigid_body_excluded_volume.py

This example shows using two rigid bodies and doing excluded volume with them.

1 ## \example domino/rigid_body_excluded_volume.py
2 # This example shows using two rigid bodies and doing excluded volume with
3 # them.
4 
5 import IMP
6 import IMP.core
7 import IMP.algebra
8 import IMP.domino
9 import IMP.atom
10 import math
11 import sys
12 
13 IMP.setup_from_argv(sys.argv, "rigid body excluded volume")
14 
15 # create a rigid body per helix
16 
17 
18 def create_representation():
19  m = IMP.Model()
20  h0 = IMP.atom.read_pdb(IMP.domino.get_example_path(
21  'helix_0.pdb'), m, IMP.atom.CAlphaPDBSelector())
22  h1 = IMP.atom.read_pdb(IMP.domino.get_example_path(
23  'helix_1.pdb'), m, IMP.atom.CAlphaPDBSelector())
24  for h in [h0, h1]:
25  IMP.atom.create_rigid_body(h)
26  return (m, [h0, h1])
27 
28 
29 def create_excluded_volume(m, helices):
30  # this is the interesting function:
31  # it uses a KClosePairsPair score to generate the list of close atoms on
32  # the fly
33  all = []
34  for h in helices:
35  all.extend(IMP.atom.get_by_type(h, IMP.atom.ATOM_TYPE))
36  lsc = IMP.container.ListSingletonContainer(m, all)
37  evr = IMP.core.ExcludedVolumeRestraint(lsc, 1)
38  evr.set_maximum_score(.01)
39  return [evr]
40 
41 
42 # creating the discrete states for domino
43 def create_discrete_states(m, helices):
44  pst = IMP.domino.ParticleStatesTable()
45  rot00 = IMP.algebra.get_rotation_about_axis(
46  IMP.algebra.Vector3D(0, 1, 0), math.pi / 2.0)
47  trs = []
48  zv = IMP.algebra.Vector3D(0, 0, 0)
49  trs.append(IMP.algebra.ReferenceFrame3D(
50  IMP.algebra.Transformation3D(rot00, zv)))
51  for dx in range(0, 15):
52  tr = IMP.algebra.Vector3D(1.0 * dx, 0, 0)
53  trs.append(IMP.algebra.ReferenceFrame3D(
54  IMP.algebra.Transformation3D(rot00, tr)))
55  pstate = IMP.domino.RigidBodyStates(trs)
56  for h in helices:
57  pst.set_particle_states(
58  IMP.core.RigidMember(h).get_rigid_body(), pstate)
59  return pst
60 
61 # setting up domino (and filters)
62 
63 
64 def create_sampler(m, rs, pst):
65  s = IMP.domino.DominoSampler(m, pst)
66  s.set_restraints(rs)
67  filters = []
68  # do not allow particles with the same ParticleStates object
69  # to have the same state index
70  filters.append(IMP.domino.ExclusionSubsetFilterTable(pst))
71  rc = IMP.domino.RestraintCache(pst)
72  rc.add_restraints(rs)
73  # filter states that score worse than the cutoffs in the Model
74  filters.append(IMP.domino.RestraintScoreSubsetFilterTable(rc))
75  states = IMP.domino.BranchAndBoundAssignmentsTable(pst, filters)
76  s.set_assignments_table(states)
77  s.set_subset_filter_tables(filters)
78  return s
79 
80 
81 def display(m, helices, name):
82  m.update()
83  w = IMP.display.PymolWriter(name)
84  for i, h in enumerate(helices):
85  g = IMP.atom.HierarchyGeometry(h)
86  g.set_color(IMP.display.get_display_color(i))
87  w.add_geometry(g)
88 
89 
90 IMP.set_log_level(IMP.SILENT)
91 print("creating representation")
92 (m, helices) = create_representation()
93 
94 print("creating score function")
95 rs = create_excluded_volume(m, helices)
96 
97 print("creating discrete states")
98 pst = create_discrete_states(m, helices)
99 
100 print("creating sampler")
101 s = create_sampler(m, rs, pst)
102 m.set_log_level(IMP.SILENT)
103 IMP.set_log_level(IMP.VERBOSE)
104 print("sampling")
105 cs = s.create_sample()
106 
107 print("found ", cs.get_number_of_configurations(), "solutions")
108 score = []
109 sf = IMP.core.RestraintsScoringFunction(rs)
110 for i in range(cs.get_number_of_configurations()):
111  cs.load_configuration(i)
112  ss = sf.evaluate(False)
113  score.append(ss)
114  print("** solution number:", i, " is:", ss)
115  display(m, helices, "sol_" + str(i) + ".pym")