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