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