Overload for the SuperposeFrame function, so that a trajectory can be...

Overload for the SuperposeFrame function, so that a trajectory can be superposed on any view from any Entity

modules/mol/alg/doc/molalg.rst

@@ -34,17 +34,39 @@
   as a new coord group.
   
   :param frames: The source coord group.
-  :param from: index of the first frame
-  :param to: index of the last frame plus one. If set to -1, the value is set to 
-     the number of frames in the coord group
+  :type frames: :class:`~ost.mol.CoordGroupHandle`
   :param sel: An entity view containing the selection of atoms to be used for     
     superposition. If set to an invalid view, all atoms in the coord group are 
     used.
+  :type sel: :class:`ost.mol.EntityView`
+  :param from: index of the first frame
+  :param to: index of the last frame plus one. If set to -1, the value is set to 
+     the number of frames in the coord group
   :param ref: The index of the reference frame to use for superposition. If set 
      to -1, the each frame is superposed to the previous frame.
      
   :returns: A newly created coord group containing the superposed frames.
 
+  
+.. function:: SuperposeFrames(frames, sel, ref_view, from=0, to=-1)
+
+  Same as SuperposeFrames above, but the superposition is done on a reference
+  view and not on another frame of the trajectory.
+  
+  :param frames: The source coord group.
+  :type frames: :class:`~ost.mol.CoordGroupHandle`
+  :param sel: An entity view containing the selection of atoms of the frames to be used for     
+    superposition.
+  :type sel: :class:`ost.mol.EntityView`
+  :param ref_view: The reference view on which the frames will be superposed. The number
+    of atoms in this reference view should be equal to the number of atoms in sel.
+  :type ref_view: :class:`ost.mol.EntityView`
+  :param from: index of the first frame
+  :param to: index of the last frame plus one. If set to -1, the value is set to 
+     the number of frames in the coord group     
+  
+  :returns: A newly created coord group containing the superposed frames.
+
 .. autofunction:: ParseAtomNames
 
 .. autofunction:: MatchResidueByNum

modules/mol/alg/pymod/wrap_mol_alg.cc

@@ -38,8 +38,12 @@ Real (*ldt_a)(const mol::EntityView&, const mol::EntityView& ref, Real, Real, co
 Real (*ldt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistTest;
 mol::EntityView (*fc_a)(const mol::EntityView&, Real,bool)=&mol::alg::FilterClashes;
 mol::EntityView (*fc_b)(const mol::EntityHandle&, Real, bool)=&mol::alg::FilterClashes;
+mol::CoordGroupHandle (*superpose_frames1)(mol::CoordGroupHandle&, mol::EntityView&, int, int, int)=&mol::alg::SuperposeFrames;
+mol::CoordGroupHandle (*superpose_frames2)(mol::CoordGroupHandle&,  mol::EntityView&, mol::EntityView&, int, int)=&mol::alg::SuperposeFrames;
 }
 
+
+
 BOOST_PYTHON_MODULE(_ost_mol_alg)
 {
   export_svdSuperPose();
@@ -52,7 +56,10 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
   def("LocalDistTest", ldt_b, (arg("ref_index")=0, arg("mdl_index")=1));
   def("FilterClashes", fc_a, (arg("ent"), arg("tolerance")=0.1, arg("always_remove_bb")=false));
   def("FilterClashes", fc_b, (arg("ent"), arg("tolerance")=0.1, arg("always_remove_bb")=false));
-  def("SuperposeFrames", &ost::mol::alg::SuperposeFrames, 
+  def("SuperposeFrames", superpose_frames1, 
       (arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0, 
        arg("end")=-1, arg("ref")=-1));
+  def("SuperposeFrames", superpose_frames2, 
+  (arg("source"), arg("sel")=ost::mol::EntityView(), arg("ref_view")=ost::mol::EntityView(),arg("begin")=0, 
+   arg("end")=-1));
 }

modules/mol/alg/src/superpose_frames.cc

@@ -75,9 +75,59 @@ inline EMatX3 col_sub(const EMatX3& m, const ECVec3& s)
   return r;
 }
 
-CoordGroupHandle SuperposeFrames(CoordGroupHandle cg, EntityView sel,
+  
+  void AddSuperposedFrame(CoordGroupHandle& superposed, EMatX3& ref_mat,EMatX3& ref_centered,ECVec3& ref_center,CoordFramePtr frame,std::vector<unsigned long>& indices)
+{
+  // This function superposes and then adds a CoordFrame (frame) to a CoordGroup (superposed).
+  // ref_mat, ref_centered and ref_center contain respectively the positions, centered positions and
+  // vector to center of the reference points for superposition.
+  // indices is a vector of the indices of the atoms to be superposed on the reference positions
+  EMat3X frame_centered=EMat3X::Zero(3, indices.size());
+  EMat3X frame_mat=EMat3X::Zero(3, indices.size());
+  ERVec3 frame_center;
+  for (size_t j=0; j<indices.size(); ++j) {
+    frame_mat.col(j)=gvec_to_rvec((*frame)[indices[j]]);
+  }
+  std::vector<geom::Vec3> frame_data=*frame;
+  frame_center=frame_mat.rowwise().sum()/frame_mat.cols();
+  frame_centered=row_sub(frame_mat, frame_center);
+  //single value decomposition
+  Eigen::SVD<EMat3> svd(frame_centered*ref_centered);
+  EMat3 matrixVT=svd.matrixV().transpose();
+  //determine rotation
+  Real detv=matrixVT.determinant();
+  Real dett=svd.matrixU().determinant();
+  Real det=detv*dett;
+  EMat3 e_rot;
+  if (det<0) {
+    EMat3 tmat=EMat3::Identity();
+    tmat(2,2)=-1;
+    e_rot=(svd.matrixU()*tmat)*matrixVT;
+  }else{
+    e_rot=svd.matrixU()*matrixVT;
+  }
+  //  prepare rmsd calculation
+  geom::Vec3 shift=rvec_to_gvec(ref_center);
+  geom::Vec3 com_vec=-cvec_to_gvec(frame_center);
+  geom::Mat3 rot=emat_to_gmat(e_rot);
+  geom::Mat4 mat4_com, mat4_rot, mat4_shift;
+  mat4_rot.PasteRotation(rot);
+  mat4_shift.PasteTranslation(shift);
+  mat4_com.PasteTranslation(com_vec);
+  geom::Mat4 tf=geom::Mat4(mat4_shift*mat4_rot*mat4_com);
+  for (std::vector<geom::Vec3>::iterator c=frame_data.begin(), 
+       e2=frame_data.end(); c!=e2; ++c) {
+    *c=geom::Vec3(tf*geom::Vec4(*c));
+  }
+  superposed.AddFrame(frame_data);
+}  
+
+
+CoordGroupHandle SuperposeFrames(CoordGroupHandle& cg, EntityView& sel,
                                  int begin, int end, int ref)
 {
+  //This function superposes the frames of a CoordGroup (cg) with indices between begin and end
+  //onto the frame with index ref. The superposition is done on a selection of atoms given by the EntityView sel.
   int real_end=end==-1 ? cg.GetFrameCount() : end;
   CoordFramePtr ref_frame;
   std::vector<unsigned long> indices;
@@ -97,10 +147,7 @@ CoordGroupHandle SuperposeFrames(CoordGroupHandle cg, EntityView sel,
   CoordGroupHandle superposed=CreateCoordGroup(cg.GetEntity().GetAtomList());
   EMatX3 ref_mat=EMatX3::Zero(indices.size(), 3);
   EMatX3 ref_centered=EMatX3::Zero(indices.size(), 3);
-  EMat3X frame_centered=EMat3X::Zero(3, indices.size());
-  EMat3X frame_mat=EMat3X::Zero(3, indices.size());
   ECVec3 ref_center;
-  ERVec3 frame_center;
   if (ref!=-1) {
     ref_frame=cg.GetFrame(ref);
     for (size_t i=0; i<indices.size(); ++i) {
@@ -125,46 +172,59 @@ CoordGroupHandle SuperposeFrames(CoordGroupHandle cg, EntityView sel,
       ref_centered=col_sub(ref_mat, ref_center);
     }
     CoordFramePtr frame=cg.GetFrame(i);
-    for (size_t j=0; j<indices.size(); ++j) {
-      frame_mat.col(j)=gvec_to_rvec((*frame)[indices[j]]);
-    }
-    std::vector<geom::Vec3> frame_data=*frame;
-    frame_center=frame_mat.rowwise().sum()/frame_mat.cols();
-    frame_centered=row_sub(frame_mat, frame_center);
-    //single value decomposition
-    Eigen::SVD<EMat3> svd(frame_centered*ref_centered);
-    EMat3 matrixVT=svd.matrixV().transpose();
-    //determine rotation
-    Real detv=matrixVT.determinant();
-    Real dett=svd.matrixU().determinant();
-    Real det=detv*dett;
-    EMat3 e_rot;
-    if (det<0) {
-      EMat3 tmat=EMat3::Identity();
-      tmat(2,2)=-1;
-      e_rot=(svd.matrixU()*tmat)*matrixVT;
-    }else{
-      e_rot=svd.matrixU()*matrixVT;
-    }
-    //  prepare rmsd calculation
-    geom::Vec3 shift=rvec_to_gvec(ref_center);
-    geom::Vec3 com_vec=-cvec_to_gvec(frame_center);
-    geom::Mat3 rot=emat_to_gmat(e_rot);
-    geom::Mat4 mat4_com, mat4_rot, mat4_shift;
-    mat4_rot.PasteRotation(rot);
-    mat4_shift.PasteTranslation(shift);
-    mat4_com.PasteTranslation(com_vec);
-    geom::Mat4 tf=geom::Mat4(mat4_shift*mat4_rot*mat4_com);
-    for (std::vector<geom::Vec3>::iterator c=frame_data.begin(), 
-         e2=frame_data.end(); c!=e2; ++c) {
-      *c=geom::Vec3(tf*geom::Vec4(*c));
-    }
-    superposed.AddFrame(frame_data);
+    AddSuperposedFrame(superposed,ref_mat,ref_centered,ref_center,frame,indices);    
     if (ref==-1) {
       ref_frame=superposed.GetFrame(superposed.GetFrameCount()-1);
     }
   }
   return superposed;
 }
-
+  
+  
+CoordGroupHandle SuperposeFrames(CoordGroupHandle& cg, EntityView& sel,
+                                  EntityView& ref_view, int begin, int end)
+{
+  //This function superposes the frames of a CoordGroup (cg) with indices between begin and end,
+  //using a selection of atoms (sel), onto an EntityView (ref_view).
+  if (!ref_view.IsValid()){
+    throw std::runtime_error("Invalid reference view");
+  }
+  int real_end=end==-1 ? cg.GetFrameCount() : end;
+  CoordFramePtr ref_frame;
+  std::vector<unsigned long> indices;
+  if (!sel.IsValid()) {
+    indices.reserve(cg.GetAtomCount());
+    for (size_t i=0;i<cg.GetAtomCount(); ++i) {
+      indices.push_back(i);
+    }
+  } else {
+    AtomViewList atoms=sel.GetAtomList();
+    indices.reserve(atoms.size());
+    for (AtomViewList::const_iterator i=atoms.begin(), 
+         e=atoms.end(); i!=e; ++i) {
+      indices.push_back(i->GetIndex());
+    }
+  }
+  if (int(indices.size())!=ref_view.GetAtomCount()){
+    throw std::runtime_error("atom counts of the two views are not equal");
+  }
+  CoordGroupHandle superposed=CreateCoordGroup(cg.GetEntity().GetAtomList());
+  EMatX3 ref_mat=EMatX3::Zero(indices.size(), 3);
+  EMatX3 ref_centered=EMatX3::Zero(indices.size(), 3);
+  ECVec3 ref_center;
+  AtomViewList atoms=ref_view.GetAtomList();
+  int i=0;
+  for (AtomViewList::const_iterator a=atoms.begin(), 
+       e=atoms.end(); a!=e; a++, i++) {
+    ref_mat.row(i)=gvec_to_cvec((*a).GetPos());
+  }
+  ref_center=ref_mat.colwise().sum()/ref_mat.rows();
+  ref_centered=col_sub(ref_mat, ref_center);
+  for (int i=begin; i<real_end; ++i) {
+    CoordFramePtr frame=cg.GetFrame(i);
+    AddSuperposedFrame(superposed,ref_mat,ref_centered,ref_center,frame,indices);    
+  }
+  return superposed;
+}
+  
 }}}

modules/mol/alg/src/superpose_frames.hh

@@ -24,11 +24,15 @@
 #include <ost/mol/alg/module_config.hh>
 namespace ost { namespace mol { namespace alg {
 
-/// \brief returns a superposed version of coord group
-CoordGroupHandle DLLEXPORT_OST_MOL_ALG SuperposeFrames(CoordGroupHandle cg, 
-                                                       EntityView sel,
+/// \brief returns a superposed version of coord group, superposed on a reference frame
+CoordGroupHandle DLLEXPORT_OST_MOL_ALG SuperposeFrames(CoordGroupHandle& cg, 
+                                                       EntityView& sel,
                                                        int begin=0, int end=-1, 
                                                        int ref=-1);
+/// \brief returns a superposed version of coord group, superposed on a reference view
+CoordGroupHandle DLLEXPORT_OST_MOL_ALG SuperposeFrames(CoordGroupHandle& cg, 
+                                                       EntityView& sel, EntityView& ref_view,
+                                                       int begin=0, int end=-1);
 }}}
 
 #endif