Ah, ok, I missed that IMP.atom.create_rigid_body actually has a fix for this since 2007. That explains why in our modellings, when we use dof, the problem has not been so big, it’s quire relieving. But Riccardo, we still have problems sometimes with rigid bodies created using dof not following symmetry in some frames, that is why I am investigating. But the problem apparently must be somewhere else, as your PMI fix works (tested with an example below). I will build on the below example to try to reproduce our problem and get back to you.
And, I would keep your PMI fix rather than create_compatible_rigid_body. The create_compatible_rigid_body runs get_transformation_aligning_first_to_second so it seems to be more computationally costly, am I correct Ben?.
I could imagine things could be slightly different but here the rb rotations are not different only slightly but rather looking like the reference frames were arbitrary. Maybe the create_rigid_body function could use the PMI fix internally?
import IMP
import IMP.core
import IMP.atom
import IMP.container
import IMP.rmf
import RMF
import sys
IMP.setup_from_argv(sys.argv, "rigid bodies")
for out_id in range(20):
m = IMP.Model()
rbs = []
hs = []
for i in range(4):
mp1 = IMP.atom.read_pdb(IMP.core.get_example_path('example_protein.pdb'), m)
chains = IMP.atom.get_by_type(mp1, IMP.atom.CHAIN_TYPE)
c = IMP.atom.Hierarchy(chains[0])
com=IMP.atom.CenterOfMass.setup_particle(IMP.Particle(m),IMP.atom.get_leaves(c))
tr=IMP.algebra.Transformation3D(IMP.algebra.get_identity_rotation_3d(),comcoor)
rf = IMP.algebra.ReferenceFrame3D(tr)
rbp=IMP.Particle(m)
rb=IMP.core.RigidBody.setup_particle(rbp,rf)
for h in IMP.atom.get_leaves(c):
rb.add_member(h.get_particle())
rbs.append(rb)
hs.append(c)
for i, rb in enumerate(rbs[1:]):
# the other 3 particles are all symmetric copies of the first
IMP.core.Reference.setup_particle(rb, rbs[0])
# the symmetry operation is rotation around the z axis
tr = IMP.algebra.Transformation3D(
IMP.algebra.get_rotation_about_axis(IMP.algebra.get_basis_vector_3d(2),
2 * 3.14 / len(rbs) * (i + 1)),
IMP.algebra.Vector3D(0, 0, 0))
sm = IMP.core.TransformationSymmetry(tr)
c = IMP.core.SingletonConstraint(sm, None, m, rb)
m.add_score_state(c)
m.update()
modelrmf = RMF.create_rmf_file('model{0}.rmf'.format(out_id))
IMP.rmf.add_hierarchies(modelrmf, hs)
IMP.rmf.save_frame(modelrmf)