RMSF analysis on trajectories.

Added a function to extract the root mean square fluctuation
of an entity on a trajectory.

modules/mol/alg/pymod/export_trajectory_analysis.cc

@@ -43,6 +43,7 @@ void export_TrajectoryAnalysis()
   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));
+  def("AnalyzeRMSF",&AnalyzeRMSF, (arg("traj"), arg("selection"), arg("first")=0, arg("last")=-1, arg("stride")=1));
   def("AnalyzeMinDistance", &AnalyzeMinDistance, (arg("traj"), arg("view1"), arg("view2"), arg("stride")=1));
   def("AnalyzeMinDistanceBetwCenterOfMassAndView", &AnalyzeMinDistanceBetwCenterOfMassAndView, (arg("traj"), arg("view_cm"), arg("view_atoms"), arg("stride")=1));
   def("AnalyzeAromaticRingInteraction", &AnalyzeAromaticRingInteraction, (arg("traj"), arg("view_ring1"), arg("view_ring2"), arg("stride")=1));

modules/mol/alg/src/trajectory_analysis.cc

@@ -331,11 +331,11 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
     if (to<from) {
       throw std::runtime_error("to smaller than from");
     }
-    unsigned int n_frames=to-from;
+    unsigned int n_frames=ceil((to-from)/stride);
     if (n_atoms==0){
       throw std::runtime_error("EntityView is empty");
     }
-    if (n_frames<=stride) {
+    if (n_frames<=1) {
       throw std::runtime_error("number of frames is too small");
     }
     std::vector<unsigned long> indices;
@@ -358,4 +358,45 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
     return eh;
   }
 
+  Real AnalyzeRMSF(const CoordGroupHandle& traj, const EntityView& selection, int from, int to, unsigned int stride)
+  // This function extracts the rmsf between two entity views and assigns it 
+  // 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.AnalyzeRMSD(t,Sele,Sele)
+  {
+    CheckHandleValidity(traj);
+    if (to==-1)to=traj.GetFrameCount();
+    unsigned int n_atoms=selection.GetAtomCount();
+    if (to<from) {
+      throw std::runtime_error("to smaller than from");
+    }
+    unsigned int n_frames=ceil((to-from)/stride);
+    if (n_atoms==0){
+      throw std::runtime_error("EntityView is empty");
+    }
+    if (n_frames<=1) {
+      throw std::runtime_error("number of frames is too small");
+    }
+    Real rmsf=0.0;
+    geom::Vec3 v;
+    std::vector<unsigned long> sele_indices;
+    std::vector<geom::Vec3> ref_pos(n_atoms,geom::Vec3(0.,0.,0.));
+    GetIndices(selection, sele_indices);
+    for (unsigned int j=0; j<n_atoms; ++j) {
+      for (int i=from; i<to; i+=stride) {
+        CoordFramePtr frame=traj.GetFrame(i);
+        ref_pos[j]+=frame->GetAtomPos(sele_indices[j]);
+      }
+      ref_pos[j]/=n_frames;
+    }
+    for (int i=from; i<to; i+=stride) {
+      CoordFramePtr frame=traj.GetFrame(i);
+      for (unsigned int j=0; j<n_atoms; ++j) {
+        v=frame->GetAtomPos(sele_indices[j])-ref_pos[j];
+        rmsf+=geom::Dot(v,v);
+      }
+    }
+    return pow(rmsf/float(n_atoms*n_frames),0.5);
+  }
+  
 }}} //ns

modules/mol/alg/src/trajectory_analysis.hh

@@ -40,6 +40,7 @@ namespace ost { namespace mol { namespace alg {
   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);
+  Real DLLEXPORT_OST_MOL_ALG AnalyzeRMSF(const CoordGroupHandle& traj, const EntityView& selection,int from=0, int to=-1, unsigned int stride=1);
   std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeMinDistance(const CoordGroupHandle& traj, const EntityView& view1, const EntityView& view2,unsigned int stride=1);
   std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeMinDistanceBetwCenterOfMassAndView(const CoordGroupHandle& traj, const EntityView& view_cm, const EntityView& view_atoms,unsigned int stride=1);
   std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, const EntityView& view_ring1, const EntityView& view_ring2,unsigned int stride=1);