Don't center models for mmCIF output by default.

For the sake of reproducing what was done in the modeling
as well as possible, it's better to use the same model
coordinates in the RMF and the mmCIF output (previously mmCIF
output was centered). This also ensures that trajectories
or ensembles are in the same reference frame as mmCIF and so work
correctly. Provide a Transformation3D for each model
so that this behavior can be overridden. (We do this for Nup84 since
the trajectories there are multimodel PDB files, which have already
been centered, so in this case it makes sense to center the mmCIF
output to match. Relates ihmwg/IHMCIF#30.)

pyext/src/mmcif.py

@@ -975,8 +975,11 @@ def get_transformation(self, model):
         r = [float(model.stats[prefix + '_Rotation%d' % i]) for i in range(4)]
         t = [float(model.stats[prefix + '_Translation%d' % i])
              for i in range(3)]
+        # If the model coordinates are transformed, need to back transform
+        # them first
+        inv = model.transform.get_inverse()
         return IMP.algebra.Transformation3D(IMP.algebra.Rotation3D(*r),
-                                            IMP.algebra.Vector3D(*t))
+                                            IMP.algebra.Vector3D(*t)) * inv
     def get_cross_correlation(self, model):
         """Get the cross correlation coefficient between the model projection
            and the image"""
@@ -1165,6 +1168,11 @@ def __init__(self, name):
 
 class _Model(object):
     def __init__(self, prot, simo, protocol, assembly, group):
+        # Transformation from IMP coordinates into mmCIF coordinate system.
+        # Normally we pass through coordinates unchanged, but in some cases
+        # we may want to translate them (e.g. Nup84, where the deposited PDB
+        # files are all centered; we want the mmCIF files to match)
+        self.transform = IMP.algebra.get_identity_transformation_3d()
         self.group = group
         # The _Protocol which produced this model
         self.protocol = protocol
@@ -1176,6 +1184,7 @@ def __init__(self, prot, simo, protocol, assembly, group):
         o._init_dictchain(name, prot)
         (particle_infos_for_pdb,
          self.geometric_center) = o.get_particle_infos_for_pdb_writing(name)
+        self.geometric_center = IMP.algebra.Vector3D(*self.geometric_center)
         self.entity_for_chain = {}
         self.comp_for_chain = {}
         for protname, chain_id in o.dictchain[name].items():
@@ -1249,14 +1258,13 @@ def dump_atoms(self, writer):
                 for atom in model.atoms:
                     (xyz, atom_type, residue_type, chain_id, residue_index,
                      all_indexes, radius) = atom
+                    pt = model.transform * xyz
                     l.write(id=ordinal, label_atom_id=atom_type.get_string(),
                             label_comp_id=residue_type.get_string(),
                             label_asym_id=chain_id,
                             label_entity_id=model.entity_for_chain[chain_id].id,
                             label_seq_id=residue_index,
-                            Cartn_x=xyz[0] - model.geometric_center[0],
-                            Cartn_y=xyz[1] - model.geometric_center[1],
-                            Cartn_z=xyz[2] - model.geometric_center[2],
+                            Cartn_x=pt[0], Cartn_y=pt[1], Cartn_z=pt[2],
                             model_id=model.id)
                     ordinal += 1
 
@@ -1273,14 +1281,13 @@ def dump_spheres(self, writer):
                      all_indexes, radius) = sphere
                     if all_indexes is None:
                         all_indexes = (residue_index,)
+                    pt = model.transform * xyz
                     l.write(ordinal_id=ordinal,
                             entity_id=model.entity_for_chain[chain_id].id,
                             seq_id_begin = all_indexes[0],
                             seq_id_end = all_indexes[-1],
                             asym_id=chain_id,
-                            Cartn_x=xyz[0] - model.geometric_center[0],
-                            Cartn_y=xyz[1] - model.geometric_center[1],
-                            Cartn_z=xyz[2] - model.geometric_center[2],
+                            Cartn_x=pt[0], Cartn_y=pt[1], Cartn_z=pt[2],
                             object_radius=radius,
                             rmsf=model.get_rmsf(model.comp_for_chain[chain_id],
                                                 all_indexes),
@@ -2244,8 +2251,6 @@ def add_replica_exchange_analysis(self, rex):
                 # best scoring one
                 m.name = 'Best scoring model'
                 m.stats = stats
-                # Don't alter original RMF coordinates
-                m.geometric_center = [0,0,0]
                 # Add RMSF info if available
                 for c in self.all_modeled_components:
                     e.load_rmsf(m, c)

test/test_mmcif.py

@@ -777,9 +777,9 @@ def test_model_dumper_sphere(self):
 _ihm_sphere_obj_site.object_radius
 _ihm_sphere_obj_site.rmsf
 _ihm_sphere_obj_site.model_id
-1 1 1 1 A 0.000 0.000 0.000 3.068 . 1
-2 1 2 2 A 0.000 0.000 0.000 2.997 . 1
-3 1 3 4 A 0.000 0.000 0.000 3.504 . 1
+1 1 1 1 A -8.986 11.688 -5.817 3.068 . 1
+2 1 2 2 A -8.986 11.688 -5.817 2.997 . 1
+3 1 3 4 A -8.986 11.688 -5.817 3.504 . 1
 #
 """)
 
@@ -833,8 +833,8 @@ def test_model_dumper_atom(self):
 _atom_site.Cartn_z
 _atom_site.label_entity_id
 _atom_site.model_id
-1 CA MET 1 A 0.000 0.000 0.000 1 1
-2 CA GLU 2 A 0.000 0.000 0.000 1 1
+1 CA MET 1 A -8.986 11.688 -5.817 1 1
+2 CA GLU 2 A -8.986 11.688 -5.817 1 1
 #
 #
 loop_
@@ -849,7 +849,7 @@ def test_model_dumper_atom(self):
 _ihm_sphere_obj_site.object_radius
 _ihm_sphere_obj_site.rmsf
 _ihm_sphere_obj_site.model_id
-1 1 3 4 A 0.000 0.000 0.000 3.504 . 1
+1 1 3 4 A -8.986 11.688 -5.817 3.504 . 1
 #
 """)
 
@@ -908,9 +908,9 @@ def test_model_dumper_sphere_rmsf(self):
 _ihm_sphere_obj_site.object_radius
 _ihm_sphere_obj_site.rmsf
 _ihm_sphere_obj_site.model_id
-1 1 1 1 A 0.000 0.000 0.000 3.068 4.500 1
-2 1 2 2 A 0.000 0.000 0.000 2.997 3.500 1
-3 1 3 4 A 0.000 0.000 0.000 3.504 5.500 1
+1 1 1 1 A -8.986 11.688 -5.817 3.068 4.500 1
+2 1 2 2 A -8.986 11.688 -5.817 2.997 3.500 1
+3 1 3 4 A -8.986 11.688 -5.817 3.504 5.500 1
 #
 """)