Trajectory analysis tools added

modules/mol/alg/doc/molalg.rst

@@ -82,3 +82,117 @@ The following function detects steric clashes in atomic structures. Two atoms ar
   :returns: The filtered :class:`~ost.mol.EntityView`
 
 
+Trajectory Analyses
+--------------------------------------------------------------------------------
+
+This is a set of funcitons used for basic trajectory analysis such as extracting positions, 
+distances, angles and RMSDs. The organization is such that most functions have their counterpart 
+at the individual frame level (CoordFrame) so that they can alsobe called on one frame instead 
+of the whole trajectory.
+
+All these functions have a "stride" argument that defaults to stride=1, which is used to skip frames in the anlysis.
+
+
+.. function:: AnalyzeAtomPos(traj, atom1, stride=1)
+
+  This function extracts the position of an atom from a trajectory. It returns 
+  a vector containing the position of the atom for each analyzed frame.
+  
+  :param traj: The trajectory to be analyzed.
+  :type traj: :class:`~ost.mol.CoordGroupHandle`
+  :param atom1: The :class:`~ost.mol.AtomHandle`.
+  :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+  
+  
+.. function:: AnalyzeCenterOfMassPos(traj, sele, stride=1)
+
+  This function extracts the position of the center-of-mass of a selection 
+  (:class:`~ost.mol.EntityView`) from a trajectory and returns it as a vector.
+  
+  :param traj: The trajectory to be analyzed.
+  :type traj: :class:`~ost.mol.CoordGroupHandle`
+  :param sele: The selection from which the center of mass is computed
+  :type sele: :class:`~ost.mol.EntityView`.
+  :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+  
+
+.. function:: AnalyzeDistanceBetwAtoms(traj, atom1, atom2, stride=1)
+
+  This function extracts the distance between two atoms from a trajectory 
+  and returns it as a vector.
+  
+  :param traj: The trajectory to be analyzed.
+  :type traj: :class:`~ost.mol.CoordGroupHandle`
+  :param atom1: The first :class:`~ost.mol.AtomHandle`.
+  :param atom2: The second :class:`~ost.mol.AtomHandle`.
+  :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.  
+
+
+.. function:: AnalyzeAngle(traj, atom1, atom2, atom3, stride=1)
+
+  This function extracts the angle between three atoms from a trajectory 
+  and returns it as a vector. The second atom is taken as being the central
+  atom, so that the angle is between the vectors (atom1.pos-atom2.pos) and
+  (atom3.pos-atom2.pos).
+  
+  :param traj: The trajectory to be analyzed.
+  :type traj: :class:`~ost.mol.CoordGroupHandle`
+  :param atom1: The first :class:`~ost.mol.AtomHandle`.
+  :param atom2: The second :class:`~ost.mol.AtomHandle`.
+  :param atom3: The third :class:`~ost.mol.AtomHandle`.
+  :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+.. function:: AnalyzeDihedralAngle(traj, atom1, atom2, atom3, atom4, stride=1)
+
+  This function extracts the dihedral angle between four atoms from a trajectory 
+  and returns it as a vector. The angle is between the planes containing the first 
+  three and the last three atoms.
+  
+  :param traj: The trajectory to be analyzed.
+  :type traj: :class:`~ost.mol.CoordGroupHandle`
+  :param atom1: The first :class:`~ost.mol.AtomHandle`.
+  :param atom2: The second :class:`~ost.mol.AtomHandle`.
+  :param atom3: The third :class:`~ost.mol.AtomHandle`.
+  :param atom4: The fourth :class:`~ost.mol.AtomHandle`.
+  :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+.. function:: AnalyzeDistanceBetwCenterOfMass(traj, sele1, sele2, stride=1)
+
+  This function extracts the distance between the center-of-mass of two selections 
+  (:class:`~ost.mol.EntityView`) from a trajectory and returns it as a vector.
+
+  :param traj: The trajectory to be analyzed.
+  :type traj: :class:`~ost.mol.CoordGroupHandle`
+  :param sele1: The selection from which the first center of mass is computed
+  :type sele1: :class:`~ost.mol.EntityView`.
+  :param sele2: The selection from which the second center of mass is computed
+  :type sele2: :class:`~ost.mol.EntityView`.
+  :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+.. function:: AnalyzeRMSD(traj, reference_view, sele_view)
+
+  This function extracts the rmsd between two :class:`~ost.mol.EntityView` and returns it as a vector
+  The views don't have to be from the same entity. The reference positions are taken directly
+  from the reference_view, evaluated only once. The positions from the sele_view are evaluated for
+  each frame.
+  If you want to compare to frame i of the trajectory t, first use t.CopyFrame(i) for example:
+  eh=io.LoadPDB(...),t=io.LoadCHARMMTraj(eh,...);sele=eh.Select(...);t.CopyFrame(0);mol.alg.AnalyzeRMSD(t,sele,sele)
+
+  :param traj: The trajectory to be analyzed.
+  :type traj: :class:`~ost.mol.CoordGroupHandle`
+  :param reference_view: The selection used as reference structure
+  :type reference_view: :class:`~ost.mol.EntityView`.
+  :param sele_view: The selection compared to the reference_view
+  :type sele_view: :class:`~ost.mol.EntityView`.
+  :param stride: Size of the increment of the frame's index between two consecutive frames analyzed.
+
+
+
+
+
+
+
+
+

modules/mol/alg/pymod/export_trajectory_analysis.cc

@@ -25,6 +25,7 @@ using namespace ost;
 using namespace ost::mol::alg;
 
 //"thin wrappers" for default parameters
+/*
 BOOST_PYTHON_FUNCTION_OVERLOADS(extractapos, ExtractAtomPosition, 2, 3)
 BOOST_PYTHON_FUNCTION_OVERLOADS(extractcmpos, ExtractCMPosition, 2, 3)
 BOOST_PYTHON_FUNCTION_OVERLOADS(extractdist, ExtractDistance, 3, 4)
@@ -32,14 +33,14 @@ BOOST_PYTHON_FUNCTION_OVERLOADS(extractang, ExtractAngle, 4, 5)
 BOOST_PYTHON_FUNCTION_OVERLOADS(extractdih, ExtractDihedral, 5, 6)
 BOOST_PYTHON_FUNCTION_OVERLOADS(extractcmdist, ExtractCMDistance, 3, 4)
 BOOST_PYTHON_FUNCTION_OVERLOADS(extractrmsd, ExtractRMSD, 3, 4)
-
+*/
 void export_TrajectoryAnalysis()
 {
-  def("ExtractAtomPosition",ExtractAtomPosition,extractapos());
-  def("ExtractCMPosition",&ExtractCMPosition,extractcmpos());
-  def("ExtractDistance",&ExtractDistance,extractdist());
-  def("ExtractAngle",&ExtractAngle,extractang());
-  def("ExtractDihedral",&ExtractDihedral,extractdih());
-  def("ExtractCMDistance",&ExtractCMDistance,extractcmdist());
-  def("ExtractRMSD",&ExtractRMSD,extractrmsd());
+  def("AnalyzeAtomPos",&AnalyzeAtomPos, (arg("traj"), arg("atom"), arg("stride")=1));
+  def("AnalyzeCenterOfMassPos",&AnalyzeCenterOfMassPos, (arg("traj"), arg("selection"), arg("stride")=1));
+  def("AnalyzeDistanceBetwAtoms",&AnalyzeDistanceBetwAtoms, (arg("traj"), arg("atom"), arg("atom"), arg("stride")=1));
+  def("AnalyzeAngle",&AnalyzeAngle, (arg("traj"), arg("atom"), arg("atom"), arg("atom"), arg("stride")=1));
+  def("AnalyzeDihedralAngle",&AnalyzeDihedralAngle, (arg("traj"), arg("atom"), arg("atom"), arg("atom"), arg("atom"), arg("stride")=1));
+  def("AnalyzeDistanceBetwCenterOfMass",&AnalyzeDistanceBetwCenterOfMass, (arg("traj"), arg("selection"), arg("selection"), arg("stride")=1));
+  def("AnalyzeRMSD",&AnalyzeRMSD, (arg("traj"), arg("reference_view"), arg("selection"), arg("stride")=1));
 }

modules/mol/alg/src/trajectory_analysis.cc

@@ -21,10 +21,8 @@
  * Author Niklaus Johner
  */
 #include <stdexcept>
-//#include <boost/bind.hpp>
 #include <ost/base.hh>
 #include <ost/geom/vec3.hh>
-//#include <ost/mol/alg/svd_superpose.hh>
 #include <ost/base.hh>
 #include <ost/geom/geom.hh>
 #include <ost/mol/entity_view.hh>
@@ -34,72 +32,71 @@
 
 namespace ost { namespace mol { namespace alg {
 
-geom::Vec3List ExtractAtomPosition(const CoordGroupHandle& traj, const AtomHandle& a1, unsigned int stride)
-//This function extracts the position of an atom returns it as a vector of geom::Vec3
-//Doesn't work in python, because it cannot create the vector of geom::Vec3
+geom::Vec3List AnalyzeAtomPos(const CoordGroupHandle& traj, const AtomHandle& a1, unsigned int stride)
+//This function extracts the position of an atom from a trajectory and returns it as a vector of geom::Vec3
 {
-  traj.CheckValidity();
+  CheckHandleValidity(traj);
   geom::Vec3List pos;
   pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
   int i1=a1.GetIndex();
   for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
     CoordFramePtr frame=traj.GetFrame(i);
-    pos.push_back(frame->GetAtomPosition(i1));
+    pos.push_back(frame->GetAtomPos(i1));
   }
   return pos;
 }
 
-std::vector<Real> ExtractDistance(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, 
+std::vector<Real> AnalyzeDistanceBetwAtoms(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, 
                                   unsigned int stride)
 //This function extracts the distance between two atoms from a trajectory and returns it as a vector
 {
-  traj.CheckValidity();
+  CheckHandleValidity(traj);
   std::vector<Real> dist;
   dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
   int i1=a1.GetIndex();
   int i2=a2.GetIndex();
   for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
     CoordFramePtr frame=traj.GetFrame(i);
-    dist.push_back((*frame).GetDistance(i1,i2));
+    dist.push_back(frame->GetDistanceBetwAtoms(i1,i2));
   }
   return dist;
 } 
 
-std::vector<Real> ExtractAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, 
+std::vector<Real> AnalyzeAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, 
                                const AtomHandle& a3, unsigned int stride)
 //This function extracts the angle between three atoms from a trajectory and returns it as a vector
 {
-  traj.CheckValidity();
+  CheckHandleValidity(traj);
   std::vector<Real> ang;
   ang.reserve(ceil(traj.GetFrameCount()/float(stride)));
   int i1=a1.GetIndex(),i2=a2.GetIndex(),i3=a3.GetIndex();
   for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
     CoordFramePtr frame=traj.GetFrame(i);
-    ang.push_back((*frame).GetAngle(i1,i2,i3));
+    ang.push_back(frame->GetAngle(i1,i2,i3));
   }
   return ang;
 }
 
-std::vector<Real> ExtractDihedral(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, 
+std::vector<Real> AnalyzeDihedralAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, 
                                   const AtomHandle& a3, const AtomHandle& a4, unsigned int stride)
 //This function extracts the dihedral angle between four atoms from a trajectory and returns it as a vector
 {
-  traj.CheckValidity();
+  CheckHandleValidity(traj);
   std::vector<Real> ang;
   ang.reserve(ceil(traj.GetFrameCount()/float(stride)));
   int i1=a1.GetIndex(),i2=a2.GetIndex(),i3=a3.GetIndex(),i4=a4.GetIndex();
   for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
     CoordFramePtr frame=traj.GetFrame(i);
-    ang.push_back((*frame).GetDihedralAngle(i1,i2,i3,i4));
+    ang.push_back(frame->GetDihedralAngle(i1,i2,i3,i4));
   }
   return ang;
 }
 
-geom::Vec3List ExtractCMPosition(const CoordGroupHandle& traj, const EntityView& Sele,unsigned int stride)
-//This function extracts the position of the CM of two selection (entity views) from a trajectory 
+geom::Vec3List AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const EntityView& Sele,unsigned int stride)
+//This function extracts the position of the CenterOfMass of a selection (entity view) from a trajectory 
 //and returns it as a vector. 
   {
-  traj.CheckValidity();
+  CheckHandleValidity(traj);
   geom::Vec3List pos;
   pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
   std::vector<unsigned long> indices;
@@ -107,17 +104,17 @@ geom::Vec3List ExtractCMPosition(const CoordGroupHandle& traj, const EntityView&
   GetIndicesAndMasses(Sele, indices, masses);
   for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
     CoordFramePtr frame=traj.GetFrame(i);
-    pos.push_back(frame->GetCMPosition(indices,masses));
+    pos.push_back(frame->GetCenterOfMassPos(indices,masses));
   }
   return pos;
 }
 
-std::vector<Real> ExtractCMDistance(const CoordGroupHandle& traj, const EntityView& Sele1,
+std::vector<Real> AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, const EntityView& Sele1,
                                     const EntityView& Sele2, unsigned int stride)
-//This function extracts the distance between the CM of two selection (entity views) from a trajectory 
+//This function extracts the distance between the CenterOfMass of two selection (entity views) from a trajectory 
 //and returns it as a vector. 
   {
-  traj.CheckValidity();
+  CheckHandleValidity(traj);
   std::vector<Real> dist;
   dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
   std::vector<unsigned long> indices1,indices2;
@@ -126,21 +123,21 @@ std::vector<Real> ExtractCMDistance(const CoordGroupHandle& traj, const EntityVi
   GetIndicesAndMasses(Sele2, indices2, masses2);
   for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
     CoordFramePtr frame=traj.GetFrame(i);
-    dist.push_back((*frame).GetCMDistance(indices1,masses1,indices2,masses2));
+    dist.push_back(frame->GetDistanceBetwCenterOfMass(indices1,masses1,indices2,masses2));
   }
   return dist;
 }
 
-std::vector<Real> ExtractRMSD(const CoordGroupHandle& traj, const EntityView& ReferenceView,
-                                    const EntityView& SeleView, unsigned int stride)
+std::vector<Real> AnalyzeRMSD(const CoordGroupHandle& traj, const EntityView& reference_view,
+                                    const EntityView& sele_view, unsigned int stride)
 // This function extracts the rmsd between two entity views and returns it as a vector
 // The views don't have to be from the same entity
 // If you want to compare to frame i of the trajectory t, first use t.CopyFrame(i) for example:
-// eh=io.LoadPDB(...),t=io.LoadCHARMMTraj(eh,...);Sele=eh.Select(...);t.CopyFrame(0);mol.alg.ExtractRMSD(t,Sele,Sele)
+// eh=io.LoadPDB(...),t=io.LoadCHARMMTraj(eh,...);Sele=eh.Select(...);t.CopyFrame(0);mol.alg.AnalyzeRMSD(t,Sele,Sele)
   {
-  traj.CheckValidity();
-  int count_ref=ReferenceView.GetAtomCount();
-  int count_sele=SeleView.GetAtomCount();
+  CheckHandleValidity(traj);
+  int count_ref=reference_view.GetAtomCount();
+  int count_sele=sele_view.GetAtomCount();
   if (count_ref!=count_sele){
     throw std::runtime_error("atom counts of the two views are not equal");
   }  
@@ -148,11 +145,11 @@ std::vector<Real> ExtractRMSD(const CoordGroupHandle& traj, const EntityView& Re
   rmsd.reserve(ceil(traj.GetFrameCount()/float(stride)));
   std::vector<unsigned long> sele_indices;
   std::vector<geom::Vec3> ref_pos;
-  GetIndices(ReferenceView, sele_indices);
-  GetPositions(SeleView, ref_pos);
+  GetIndices(reference_view, sele_indices);
+  GetPositions(sele_view, ref_pos);
   for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
     CoordFramePtr frame=traj.GetFrame(i);
-    rmsd.push_back((*frame).GetRMSD(ref_pos,sele_indices));
+    rmsd.push_back(frame->GetRMSD(ref_pos,sele_indices));
   }
   return rmsd;
 }

modules/mol/alg/src/trajectory_analysis.hh

@@ -33,13 +33,13 @@
 
 namespace ost { namespace mol { namespace alg {
   
-  geom::Vec3List DLLEXPORT_OST_MOL_ALG ExtractAtomPosition(const CoordGroupHandle& traj, const AtomHandle& a1,unsigned int stride=1);
-  geom::Vec3List DLLEXPORT_OST_MOL_ALG ExtractCMPosition(const CoordGroupHandle& traj, const EntityView& Sele,unsigned int stride=1);
-  std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractDistance(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,unsigned int stride=1);
-  std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3,unsigned int stride=1);
-  std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractCMDistance(const CoordGroupHandle& traj, const EntityView& Sele1, const EntityView& Sele2,unsigned int stride=1);
-  std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractDihedral(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4,unsigned int stride=1);
-  std::vector<Real> DLLEXPORT_OST_MOL_ALG ExtractRMSD(const CoordGroupHandle& traj, const EntityView& Reference_View, const EntityView& Sele,unsigned int stride=1);
+  geom::Vec3List DLLEXPORT_OST_MOL_ALG AnalyzeAtomPos(const CoordGroupHandle& traj, const AtomHandle& a1,unsigned int stride=1);
+  geom::Vec3List DLLEXPORT_OST_MOL_ALG AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const EntityView& sele,unsigned int stride=1);
+  std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDistanceBetwAtoms(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,unsigned int stride=1);
+  std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3,unsigned int stride=1);
+  std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, const EntityView& sele1, const EntityView& sele2,unsigned int stride=1);
+  std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeDihedralAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4,unsigned int stride=1);
+  std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeRMSD(const CoordGroupHandle& traj, const EntityView& reference_view, const EntityView& sele,unsigned int stride=1);
   
   
 }}}//ns

modules/mol/base/pymod/CMakeLists.txt

@@ -5,6 +5,7 @@ export_bond.cc
 export_chain.cc
 export_chain_view.cc
 export_coord_group.cc
+export_coord_frame.cc
 export_editors.cc
 export_entity.cc
 export_entity_view.cc

modules/mol/base/pymod/wrap_mol.cc

@@ -39,6 +39,7 @@ void export_AtomView();
 void export_ResidueView();
 void export_Editors();
 void export_CoordGroup();
+void export_CoordFrame();
 void export_PropertyID();
 void export_BoundingBox();
 void export_QueryViewWrapper();
@@ -65,6 +66,7 @@ BOOST_PYTHON_MODULE(_ost_mol)
   export_EntityView(); 
   export_Editors();
   export_CoordGroup();
+  export_CoordFrame();
   export_PropertyID();  
   export_BoundingBox();
   export_QueryViewWrapper();

modules/mol/base/src/coord_frame.cc

@@ -31,36 +31,53 @@
 
 namespace ost { namespace mol {
 
-  geom::Vec3 CoordFrame::GetAtomPosition(const AtomHandle& atom){
-    return this->GetAtomPosition(atom.GetIndex());
+  CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos)
+  {
+    return CoordFrame(atom_pos);
   }
-  geom::Vec3 CoordFrame::GetAtomPosition(const int& index){
-    return (*this)[index];
+  
+  geom::Vec3 CoordFrame::GetAtomPos(const AtomHandle& atom){
+  //Returns the position of the atom in this frame
+    return this->GetAtomPos(atom.GetIndex());
+  }
+  geom::Vec3 CoordFrame::GetAtomPos(int i1){
+  //Returns the position in this frame of the atom with index i1 
+    return (*this)[i1];
   }
   
-  Real CoordFrame::GetDistance(const AtomHandle& a1, const AtomHandle& a2){
-    return this->GetDistance(a1.GetIndex(),a2.GetIndex());
+  Real CoordFrame::GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2){
+  //Return the distance in this frame between two atoms
+    return this->GetDistanceBetwAtoms(a1.GetIndex(),a2.GetIndex());
   }
-  Real CoordFrame::GetDistance(const int& i1, const int& i2){
+  Real CoordFrame::GetDistanceBetwAtoms(int i1, int i2){
+  //Return the distance in this frame between the two atoms with indices i1 and i2 
     return geom::Distance((*this)[i1],(*this)[i2]);
   }
   
   Real CoordFrame::GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3){
+  //Returns the Angle between the three atoms, a2 being considered as the central atom
+  //i.e the angle between the vector (a1.pos-a2.pos) and (a3.pos-a2.pos)
     return this->GetAngle(a1.GetIndex(),a2.GetIndex(),a3.GetIndex());  
   }
-  Real CoordFrame::GetAngle(const int& i1, const int& i2, const int& i3){
+  Real CoordFrame::GetAngle(int i1, int i2, int i3){
+  //Returns the angl between the three atoms with indices i1,i2,i3
     return geom::Angle((*this)[i1]-(*this)[i2],(*this)[i3]-(*this)[i2]);
   }
   
-  Real CoordFrame::GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4){
+  Real CoordFrame::GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, 
+                                    const AtomHandle& a3, const AtomHandle& a4){
+  //Returns the Dihedral angle between the four atoms a1,a2,a3,a4
     return this->GetDihedralAngle(a1.GetIndex(),a2.GetIndex(),a3.GetIndex(),a4.GetIndex());
   }
-  Real CoordFrame::GetDihedralAngle(const int& i1, const int& i2, const int& i3, const int& i4){
+  Real CoordFrame::GetDihedralAngle(int i1, int i2, int i3, int i4){
+  //Returns the Dihedral angle between the four atoms with indices i1,i2,i3,i4
     return geom::DihedralAngle((*this)[i1],(*this)[i2],(*this)[i3],(*this)[i4]);
   }
 
-  void GetIndices(const EntityView& Sele, std::vector<unsigned long>& indices){
-    AtomViewList atoms=Sele.GetAtomList();
+  void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices){
+  // This function returns a vector containing the atom indices of the atoms in an EntityView
+  // It is used to accelerate the extraction of information from a trajectory
+    AtomViewList atoms=sele.GetAtomList();
     indices.reserve(atoms.size());
     for (AtomViewList::const_iterator i=atoms.begin(), 
          e=atoms.end(); i!=e; ++i) {
@@ -68,9 +85,11 @@ namespace ost { namespace mol {
     }
   }
   
-  void GetMasses(const EntityView& Sele, std::vector<Real>& masses){
+  void GetMasses(const EntityView& sele, std::vector<Real>& masses){
+  // This function returns a vector containing the atom masses of the atoms in an EntityView
+  // It is used together with GetIndices to accelerate the extraction of RMSD from a trajectory
     Real mass_tot=0.0;
-    AtomViewList atoms=Sele.GetAtomList();
+    AtomViewList atoms=sele.GetAtomList();
     masses.reserve(atoms.size());
     for (AtomViewList::const_iterator i=atoms.begin(), 
          e=atoms.end(); i!=e; ++i) {
@@ -84,26 +103,30 @@ namespace ost { namespace mol {
   }
   
   
-  void GetIndicesAndMasses(const EntityView& Sele, std::vector<unsigned long>& indices,std::vector<Real>& masses){
-    GetIndices(Sele, indices);
-    GetMasses(Sele, masses);
+  void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses){
+    GetIndices(sele, indices);
+    GetMasses(sele, masses);
   }
   
-  void GetPositions(const EntityView& Sele, std::vector<geom::Vec3>& ref_pos){
-    ref_pos.reserve(Sele.GetAtomCount());
-    for (mol::AtomViewIter a=Sele.AtomsBegin(),e=Sele.AtomsEnd(); a!=e; ++a) {
+  void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos){
+  //Returns the positions of all the atoms in the EntityView
+    ref_pos.reserve(sele.GetAtomCount());
+    for (mol::AtomViewIter a=sele.AtomsBegin(),e=sele.AtomsEnd(); a!=e; ++a) {
       ref_pos.push_back((*a).GetPos());
     }
   }
   
-  geom::Vec3 CoordFrame::GetCMPosition(const EntityView& Sele){
+  geom::Vec3 CoordFrame::GetCenterOfMassPos(const EntityView& sele){
+  //Returns the position of the centor of mass of the atoms in the EntityView
     std::vector<unsigned long> indices;
     std::vector<Real> masses;
-    GetIndicesAndMasses(Sele,indices,masses);
-    return this->GetCMPosition(indices,masses);
+    GetIndicesAndMasses(sele,indices,masses);
+    return this->GetCenterOfMassPos(indices,masses);
   }
   
-  geom::Vec3 CoordFrame::GetCMPosition(std::vector<unsigned long>& indices,std::vector<Real>& masses){
+  geom::Vec3 CoordFrame::GetCenterOfMassPos(std::vector<unsigned long>& indices,std::vector<Real>& masses){
+  //Returns the position of the centor of mass of the atoms from which the indices and masses are passed 
+  //as vectors in argument
     geom::Vec3 v;
     for (unsigned int i=0,e=indices.size();i!=e; i++) {
       v+=masses[i]*(*this)[indices[i]];
@@ -111,22 +134,27 @@ namespace ost { namespace mol {
     return v;
   }
   
-  Real CoordFrame::GetCMDistance(const EntityView& Sele1,const EntityView& Sele2){
+  Real CoordFrame::GetDistanceBetwCenterOfMass(const EntityView& sele1,const EntityView& sele2){
+  //Returns the distance between the centers of mass of the two EntityViews
     std::vector<unsigned long> indices1,indices2;
     std::vector<Real> masses1,masses2;
-    GetIndicesAndMasses(Sele1,indices1,masses1);
-    GetIndicesAndMasses(Sele2,indices2,masses2);
-    return this->GetCMDistance(indices1,masses1,indices2,masses2);
+    GetIndicesAndMasses(sele1,indices1,masses1);
+    GetIndicesAndMasses(sele2,indices2,masses2);
+    return this->GetDistanceBetwCenterOfMass(indices1,masses1,indices2,masses2);
   }
   
-  Real CoordFrame::GetCMDistance(std::vector<unsigned long>& indices1,std::vector<Real>& masses1,
+  Real CoordFrame::GetDistanceBetwCenterOfMass(std::vector<unsigned long>& indices1,std::vector<Real>& masses1,
                            std::vector<unsigned long>& indices2,std::vector<Real>& masses2){
-    geom::Vec3 v1=this->GetCMPosition(indices1, masses1);
-    geom::Vec3 v2=this->GetCMPosition(indices2, masses2);
+  //Returns the distance between the centers of mass of the two groups of atoms from which the
+  //indices and masses are given as vectors as argument
+    geom::Vec3 v1=this->GetCenterOfMassPos(indices1, masses1);
+    geom::Vec3 v2=this->GetCenterOfMassPos(indices2, masses2);
     return geom::Distance(v1,v2);
   }
   
   Real CoordFrame::GetRMSD(const std::vector<geom::Vec3>& ref_pos,const std::vector<unsigned long>& indices_sele){
+  //Returns the RMSD between the positions of the atoms whose indices are given in indices_sele and the positions
+  //given in ref_pos
     Real rmsd=0.0,val;
     for (unsigned int i1=0; i1!=indices_sele.size(); ++i1) {
       geom::Vec3 av_sele=(*this)[indices_sele[i1]];
@@ -134,19 +162,21 @@ namespace ost { namespace mol {
       val=geom::Length2(av_ref-av_sele);
       rmsd+=val;
     }
-    return rmsd/indices_sele.size();
+    return pow(rmsd/indices_sele.size(),0.5);
   }
   
-  Real CoordFrame::GetRMSD(const EntityView& Reference_View,const EntityView& Sele_View){
-    int count_ref=Reference_View.GetAtomCount();
-    int count_sele=Sele_View.GetAtomCount();
+  Real CoordFrame::GetRMSD(const EntityView& reference_view,const EntityView& sele_view){
+  //Return the rmsd between two EntityViews. The reference positions are taken directly from the reference_view
+  //whereas they are taken from this CoordFrame for the sele_view
+    int count_ref=reference_view.GetAtomCount();
+    int count_sele=sele_view.GetAtomCount();
     if (count_ref!=count_sele){
       throw std::runtime_error("atom counts of the two views are not equal");
     }
     std::vector<unsigned long> indices_sele;
     std::vector<geom::Vec3> ref_pos;
-    GetIndices(Sele_View,indices_sele);
-    GetPositions(Reference_View,ref_pos);
+    GetIndices(sele_view,indices_sele);
+    GetPositions(reference_view,ref_pos);
     return this->GetRMSD(ref_pos,indices_sele);
   }
   

modules/mol/base/src/coord_frame.hh

@@ -40,30 +40,34 @@ public:
   CoordFrame(const base_type& rhs) : base_type(rhs) { }
   CoordFrame(const std::vector<geom::Vec3>& rhs) : base_type(rhs) { }
   
-  geom::Vec3 GetAtomPosition(const AtomHandle& atom);
-  geom::Vec3 GetAtomPosition(const int& atom_index);
-  Real GetDistance(const AtomHandle& a1, const AtomHandle& a2);
-  Real GetDistance(const int& atom1_index, const int& atom2_index);
+  geom::Vec3 GetAtomPos(const AtomHandle& atom);
+  geom::Vec3 GetAtomPos(int atom_index);
+  Real GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2);
+  Real GetDistanceBetwAtoms(int atom1_index, int atom2_index);
   Real GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3);
-  Real GetAngle(const int& atom1_index, const int& atom2_index, const int& atom3_index);  
+  Real GetAngle(int atom1_index, int atom2_index, int atom3_index);  
   Real GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4);
-  Real GetDihedralAngle(const int& a1_index, const int& a2_index, const int& a3_index, const int& a4_index); 
-  geom::Vec3 GetCMPosition(const mol::EntityView& Sele);
-  geom::Vec3 GetCMPosition(std::vector<unsigned long>& indices,std::vector<Real>& masses);
-  Real GetCMDistance(const mol::EntityView& Sele1, const mol::EntityView& Sele2);
-  Real GetCMDistance(std::vector<unsigned long>& indices1,std::vector<Real>& masses1,std::vector<unsigned long>& indices2,std::vector<Real>& masses2);
-  Real GetRMSD(const std::vector<geom::Vec3>& ref_pos,const std::vector<unsigned long>& indices_sele);
-  Real GetRMSD(const mol::EntityView& Reference_View, const mol::EntityView& Sele_View);
+  Real GetDihedralAngle(int a1_index, int a2_index, int a3_index, int a4_index); 
+  geom::Vec3 GetCenterOfMassPos(const mol::EntityView& sele);
+  geom::Vec3 GetCenterOfMassPos(std::vector<unsigned long>& indices, std::vector<Real>& masses);
+  Real GetDistanceBetwCenterOfMass(const mol::EntityView& sele1, const mol::EntityView& sele2);
+  Real GetDistanceBetwCenterOfMass(std::vector<unsigned long>& indices1, std::vector<Real>& masses1, 
+                                   std::vector<unsigned long>& indices2, std::vector<Real>& masses2);
+  Real GetRMSD(const std::vector<geom::Vec3>& ref_pos, const std::vector<unsigned long>& indices_sele);
+  Real GetRMSD(const mol::EntityView& Reference_View, const mol::EntityView& sele_View);
 };
   
-  void GetIndices(const EntityView& Sele, std::vector<unsigned long>& indices);
-  void GetMasses(const EntityView& Sele,std::vector<Real>& masses);
-  void GetIndicesAndMasses(const EntityView& Sele, std::vector<unsigned long>& indices,std::vector<Real>& masses);
-  void GetPositions(const EntityView& Sele, std::vector<geom::Vec3>& ref_pos);
+  void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices);
+  void GetMasses(const EntityView& sele, std::vector<Real>& masses);
+  void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses);
+  void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos);
   
 typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
 typedef std::vector<CoordFramePtr> CoordFrameList;
 
+// factory method
+// create a frame froma Vec3List containing the positions of the atoms
+  DLLEXPORT_OST_MOL CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos);
 
 }}
 

modules/mol/base/src/coord_group.cc

@@ -84,8 +84,9 @@ CoordGroupHandle::operator bool() const
   return this->IsValid();
 }
 
-void CoordGroupHandle::AddFrame(const std::vector<geom::Vec3>& clist)
-{
+//void CoordGroupHandle::AddFrame(const std::vector<geom::Vec3>& clist)
+  void CoordGroupHandle::AddFrame(const geom::Vec3List& clist)
+  {
   this->CheckValidity();
   if (source_->IsMutable()) {
     source_->AddFrame(clist);    

modules/mol/base/src/coord_group.hh

@@ -25,7 +25,6 @@
 
 #include <vector>
 #include <boost/shared_array.hpp>
-//#include <ost/mol/alg/trajectory_analysis.hh>
 #include "atom_handle.hh"
 #include "coord_source.hh"
 
@@ -65,7 +64,8 @@ public:
   void Capture(uint frame);
   
   /// \brief add frame 
-  void AddFrame(const std::vector<geom::Vec3>& clist);
+  //void AddFrame(const std::vector<geom::Vec3>& clist);
+  void AddFrame(const geom::Vec3List& clist);
 
   void AddFrames(const CoordGroupHandle& cg);
   /// \brief set an indidivial atom position in the given frame
@@ -89,13 +89,9 @@ public:
   
   CoordGroupHandle(CoordSourcePtr source);
 
-  /*
-  //friend geom::Vec3List mol::alg::ExtractAtomPosition(const CoordGroupHandle& traj, const AtomHandle& a1,unsigned int stride=1);
-  */
 
-//private:
-  void CheckValidity() const;
 private:
+  void CheckValidity() const;
   CoordSourcePtr source_;
 };