Re: [IMP-users] getting the MCCGsampler to work (Josh Bullock)

15 Jul 2014


      hi list, sorry, just wanted to clarify a little more my ( lack of )
understanding re particles ...
so when i create a rigid body, i'm decorating a selection of particles -
but it's the particles that hold all the information. and what i want to
add to the root hierarchy are the decorated particles, not the decorator
itself, correct ? i think that this is why i was getting no mass before.
so when i load a pdb into IMP, a particle is created for each atom ? I've
been loading 4 pdb structures into IMP with a total number of 308 amino
acids, but i seem to be creating at least 10,515 particles ( shown by
printing IMP.kernel.Particle.get_index(IMP.Particle(m)) ) and i can't make
that add up - conceptually i'm a bit lost ...
and finally is there a way to access each set of decorated particles
iteratively ? i think my problems with creating the RigidBodyMover are
because i'm trying to add the entire model at once ( as shown ) - and there
are some particles in there which i haven't been able to decorate.
rb_mover = IMP.core.RigidBodyMover(m,
IMP.kernel.Particle.get_index(IMP.Particle(m)), maxTranslation, maxRotation
)
eagerly anticipating enlightenment,
josh
p.s. full script for reference
https://gist.github.com/mysticvision/abb8d42e24f329b5388b
On 14 July 2014 17:09, imp-users-request@salilab.org wrote:
> Send IMP-users mailing list submissions to
>         imp-users@salilab.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
>         https://salilab.org/mailman/listinfo/imp-users
> or, via email, send a message with subject or body 'help' to
>         imp-users-request@salilab.org
>
> You can reach the person managing the list at
>         imp-users-owner@salilab.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of IMP-users digest..."
>
>
> Today's Topics:
>
>    1. Re: getting the MCCGsampler to work (Josh Bullock)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Mon, 14 Jul 2014 17:09:08 +0100
> From: Josh Bullock jma.bullock@gmail.com
> To: imp-users@salilab.org
> Subject: Re: [IMP-users] getting the MCCGsampler to work
> Message-ID:
>         <
> CAHh_40-ROvuH3Koy+cyBvjpi127Y67VBLk3pAo-FX0HEFYJpWg@mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
>
> ah i see, yes i had no movers or optimizers ...
>
> so now i'm trying to set up a RigidBodyMover but i'm struggling to create
> the rigid bodies. There seem to be two different ways of doing this, either
> IMP.atom.create_rigid_body or IMP.atom.setup_as_rigid_body.
>
> If I use create_rigid_body then when i later try to add restraints i get
> told "leaf 'A' ... does not have mass"
>
> If I use setup_as_rigid_body i can add the restraints but i when i try to
> add it to the RigidBodyMover i get told "Rigid body passed to
> RigidBodyMover must be set to be optimized" even
> though rb.set_coordinates_are_optimized(True)
>
> what is the function of choice here ?
>
> thanks !
>
> josh
>
> ------------------------------
>
> def create_representation():
>     m= IMP.Model()
>     all=IMP.atom.Hierarchy.setup_particle(IMP.Particle(m))
>     all.set_name("the universe")
>     # create a protein, represented as a set of connected balls of
> appropriate
>     # radii and number, chose by the resolution parameter and the number of
>     # amino acids.
>
>     def create_protein_from_pdbs(name, files):
>
>         def create_from_pdb(file):
>             sls=IMP.SetLogState(IMP.NONE)
>             datadir = os.getcwd()
>             print datadir
>     t=IMP.atom.read_pdb( datadir+'/' + file, m,
>                                  IMP.atom.ATOMPDBSelector())
>             del sls
>             #IMP.atom.show_molecular_hierarchy(t)
>             c=IMP.atom.Chain(IMP.atom.get_by_type(t,
> IMP.atom.CHAIN_TYPE)[0])
>
>             # there is no reason to use all atoms, just approximate the pdb
> shape instead
>             s=IMP.atom.create_simplified_along_backbone(c,1)
>
>             # make the simplified structure rigid
>             rb=IMP.atom.create_rigid_body(s)
>             #rb=IMP.atom.setup_as_rigid_body(s)
>             #rb=IMP.atom.create_rigid_body(c)
>             rb.set_coordinates_are_optimized(True)
>             print rb.get_coordinates_are_optimized()
>
>             return rb
>
>         h= create_from_pdb(files[0])
>         h.set_name(name)
>         all.add_child(h)
>
>     create_protein_from_pdbs("A", [Firstpdb])
>
> (m,all)= create_representation()
>
>
> On 11 July 2014 19:42, imp-users-request@salilab.org wrote:
>
> > Send IMP-users mailing list submissions to
> >         imp-users@salilab.org
> >
> > To subscribe or unsubscribe via the World Wide Web, visit
> >         https://salilab.org/mailman/listinfo/imp-users
> > or, via email, send a message with subject or body 'help' to
> >         imp-users-request@salilab.org
> >
> > You can reach the person managing the list at
> >         imp-users-owner@salilab.org
> >
> > When replying, please edit your Subject line so it is more specific
> > than "Re: Contents of IMP-users digest..."
> >
> >
> > Today's Topics:
> >
> >    1. getting the MCCGsampler to work (Josh Bullock)
> >    2. Re: getting the MCCGsampler to work (Ben Webb)
> >    3. Re: getting the MCCGsampler to work (Ben Webb)
> >
> >
> > ----------------------------------------------------------------------
> >
> > Message: 1
> > Date: Fri, 11 Jul 2014 13:12:01 +0100
> > From: Josh Bullock jma.bullock@gmail.com
> > To: imp-users@salilab.org
> > Subject: [IMP-users] getting the MCCGsampler to work
> > Message-ID:
> >         <
> > CAHh_40_B9Ciag3LMdDMgFLJf2MCx2-tGPp3tmMke9MQriHwMJA@mail.gmail.com>
> > Content-Type: text/plain; charset="utf-8"
> >
> > Hello again,
> >
> > So I have the same overall problem as before - creating an ensemble of a
> > 4-subunit complex using MSconnectivity restraints. Having visualised the
> > output (via RMF - thanks Barak :?) , it's clear that no matter how many
> > steps of CG or MC I put, the models do not change from their initial
> random
> > placement. I know that the restraints are present because the models are
> > evaluated and scored appropriately.
> >
> > So I saw on an old (2011) nup84 example that MCCG can't handle rigid
> > bodies, is this still the case ? If so, should I switch to the DOMINO
> > sampler ? or it does work and likely there's an error in my code ...
> >
> > Many thanks !
> >
> > Josh
> >
> > -------------------------------------------------
> >
> > import IMP
> > import IMP.atom
> > import IMP.rmf
> > import inspect
> > import IMP.container
> > import IMP.display
> > import IMP.statistics
> > #import IMP.example
> > import sys, math, os, optparse
> > import RMF
> >
> > from optparse import OptionParser
> >
> >
> > # Convert the arguments into strings and number
> > Firstpdb = str(sys.argv[1])
> > Secondpdb = str(sys.argv[2])
> > Thirdpdb = str(sys.argv[3])
> > Fourthpdb = str(sys.argv[4])
> > models = float(sys.argv[5])
> >
> > #*****************************************
> >
> > # the spring constant to use, it doesnt really matter
> > k=100
> > # the target resolution for the representation, this is used to specify
> how
> > detailed
> > # the representation used should be
> > resolution=300
> > # the box to perform everything
> > bb=IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(0,0,0),
> >                              IMP.algebra.Vector3D(100, 100, 100))
> >
> >
> > # this function creates the molecular hierarchies for the various
> involved
> > proteins
> > def create_representation():
> >     m= IMP.Model()
> >     all=IMP.atom.Hierarchy.setup_particle(IMP.Particle(m))
> >     all.set_name("the universe")
> >     # create a protein, represented as a set of connected balls of
> > appropriate
> >     # radii and number, chose by the resolution parameter and the number
> of
> >     # amino acids.
> >
> >     def create_protein_from_pdbs(name, files):
> >
> >         def create_from_pdb(file):
> >             sls=IMP.SetLogState(IMP.NONE)
> >             datadir = os.getcwd()
> >             print datadir
> >     t=IMP.atom.read_pdb( datadir+'/' + file, m,
> >                                  IMP.atom.ATOMPDBSelector())
> >             del sls
> >             #IMP.atom.show_molecular_hierarchy(t)
> >             c=IMP.atom.Chain(IMP.atom.get_by_type(t,
> > IMP.atom.CHAIN_TYPE)[0])
> >             if c.get_number_of_children()==0:
> >                 IMP.atom.show_molecular_hierarchy(t)
> >             # there is no reason to use all atoms, just approximate the
> pdb
> > shape instead
> >             s=IMP.atom.create_simplified_along_backbone(c,
> >                                                         1)
> >             #IMP.atom.destroy(t)
> >             # make the simplified structure rigid
> >             rb=IMP.atom.create_rigid_body(s)
> >             rb=IMP.atom.create_rigid_body(c)
> >             rb.set_coordinates_are_optimized(True)
> >             return s       # <------- swapping c with s will give a
> coarse
> > grain representation - much faster !
> > #            return c
> >
> >         h= create_from_pdb(files[0])
> >         h.set_name(name)
> >         all.add_child(h)
> >
> >     create_protein_from_pdbs("A", [Firstpdb])
> >     create_protein_from_pdbs("B", [Secondpdb])
> >     create_protein_from_pdbs("C", [Thirdpdb])
> >     create_protein_from_pdbs("D", [Fourthpdb])
> >     #create_protein_from_pdbs("C", ["rpt3_imp.pdb"])
> >     return (m, all)
> >
> > # create the needed restraints and add them to the model
> >
> > def create_restraints(m, all):
> >     def add_connectivity_restraint(s):
> >
> >         tr= IMP.core.TableRefiner()
> >         rps=[]
> >         for sc in s:
> >             ps= sc.get_selected_particles()
> >             rps.append(ps[0])
> >             tr.add_particle(ps[0], ps)
> >
> >         # duplicate the IMP.atom.create_connectivity_restraint
> > functionality
> >
> >         score=
> >
> >
> IMP.core.KClosePairsPairScore(IMP.core.HarmonicSphereDistancePairScore(0,1),tr)
> >
> >         #ub = IMP.core.HarmonicUpperBound(1.0, 0.1)
> >         #ss = IMP.core.DistancePairScore(ub)
> >
> >         r= IMP.core.MSConnectivityRestraint(m,score)
> >
> >         iA = r.add_type([rps[0]])
> >         iB = r.add_type([rps[1]])
> >         iC = r.add_type([rps[2]])
> >         iD = r.add_type([rps[3]])
> >         #n1 = r.add_composite([iA, iB])
> >         n1 = r.add_composite([iA, iB, iC, iD])
> >         n2 = r.add_composite([iA, iB],n1)
> >         n3 = r.add_composite([iB, iD],n1)
> >         n4 = r.add_composite([iA, iB, iC],n1)
> >         n5 = r.add_composite([iB, iC, iD],n1)
> >
> >         m.add_restraint(r)
> >
> >     evr=IMP.atom.create_excluded_volume_restraint([all])
> >     m.add_restraint(evr)
> >     # a Selection allows for natural specification of what the restraints
> > act on
> >     S= IMP.atom.Selection
> >     sA=S(hierarchy=all, molecule="A")
> >     sB=S(hierarchy=all, molecule="B")
> >     sC=S(hierarchy=all, molecule="C")
> >     sD=S(hierarchy=all, molecule="D")
> >     add_connectivity_restraint([sA, sB, sC, sD])
> >
> >     nbl = IMP.container.ClosePairContainer([all], 0, 2)
> >     h = IMP.core.HarmonicLowerBound(0, 1)
> >     sd = IMP.core.SphereDistancePairScore(h)
> >     # use the lower bound on the inter-sphere distance to push the
> spheres
> > apart
> >     nbr = IMP.container.PairsRestraint(sd, nbl)
> >     m.add_restraint(nbr)
> >
> >
> >    # r1 = IMP.core.ExcludedVolumeRestraint(all)
> >    # m.add_restraint(r1)
> >
> >
> > # find acceptable conformations of the model
> > def get_conformations(m):
> >     sampler= IMP.core.MCCGSampler(m)
> >     sampler.set_bounding_box(bb)
> >     # magic numbers, experiment with them and make them large enough for
> > things to work
> >     sampler.set_number_of_conjugate_gradient_steps(200)
> >     sampler.set_number_of_monte_carlo_steps(40)
> >     sampler.set_number_of_attempts(models)
> >     # We don't care to see the output from the sampler
> >     #sampler.set_log_level(IMP.SILENT)
> >     # return the IMP.ConfigurationSet storing all the found
> configurations
> > that
> >     # meet the various restraint maximum scores.
> >     cs= sampler.create_sample()
> >     return cs
> >
> >
> > # cluster the conformations and write them to a file
> > def analyze_conformations(cs, all, gs):
> >     # we want to cluster the configurations to make them easier to
> > understand
> >     # in this case, the clustering is pretty meaningless
> >     embed= IMP.statistics.ConfigurationSetXYZEmbedding(cs,
> >
> >  IMP.container.ListSingletonContainer(IMP.atom.get_leaves(all)), True)
> >     cluster= IMP.statistics.create_lloyds_kmeans(embed, 10, 10000)
> >     # dump each cluster center to a file so it can be viewed.
> >     for i in range(cluster.get_number_of_clusters()):
> >         center= cluster.get_cluster_center(i)
> >         cs.load_configuration(i)
> >         h = IMP.atom.Hierarchy.get_children(all)
> >         #tfn = IMP.create_temporary_file_name("cluster%d"%i, ".rmf")
> >         huh = "./models/CLUSTER%d"%i
> >         huh = huh +".rmf"
> >         #print "file is", tfn
> >         print "file is", huh
> >         rh = RMF.create_rmf_file(huh)
> >
> >
> >         IMP.rmf.add_hierarchies(rh, h)
> >
> >         # add the current configuration to the file as frame 0
> >         IMP.rmf.save_frame(rh)
> >
> >         #for g in gs:
> >          #   rh.add_geometry(g)
> >
> >
> >
> >
> #******************************************************************************************
> > # now do the actual work
> >
> > (m,all)= create_representation()
> > #IMP.atom.show_molecular_hierarchy(all)
> > create_restraints(m, all)
> >
> > # in order to display the results, we need something that maps the
> > particles onto
> > # geometric objets. The IMP.display.Geometry objects do this mapping.
> > # IMP.display.XYZRGeometry map an IMP.core.XYZR particle onto a sphere
> > gs=[]
> > for i in range(all.get_number_of_children()):
> >     color= IMP.display.get_display_color(i)
> >     n= all.get_child(i)
> >     name= n.get_name()
> >     g= IMP.atom.HierarchyGeometry(n)
> >     g.set_color(color)
> >     gs.append(g)
> >
> > cs= get_conformations(m)
> >
> > print "found", cs.get_number_of_configurations(), "solutions"
> >
> > ListScores = []
> > for i in range(0, cs.get_number_of_configurations()):
> >         cs.load_configuration(i)
> >         # print the configuration
> >         print "solution number: ",i,"scored :", m.evaluate(False)
> >         ListScores.append(m.evaluate(False))
> >
> > # for each of the configuration, dump it to a file to view in pymol
> > for i in range(0, cs.get_number_of_configurations()):
> >     cs.load_configuration(i)
> >     h = IMP.atom.Hierarchy.get_children(all)
> >     #tfn = IMP.create_temporary_file_name("josh%d"%i, ".rmf")
> >     #print "file is", tfn
> >     huh = "./models/IMP%d"%i
> >     huh = huh +".rmf"
> >     print "file is", huh
> >     rh = RMF.create_rmf_file(huh)
> >
> >     # add the hierarchy to the file
> >     IMP.rmf.add_hierarchies(rh, h)
> >
> >     # add the current configuration to the file as frame 0
> >     IMP.rmf.save_frame(rh)
> >
> >     #for g in gs:
> >      #   w.add_geometry(g)
> >
> > analyze_conformations(cs, all, gs)
> >