added diffusion calculation

[alexxy/gromacs-testing.git] / src / spacetimecorr.cpp

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index 99b0cf3cb0beb39b443d308618518bfb0f73f927..98ef82d9ffd5bda8f853e00b7268172b8496c5a3 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -130,6 +130,20 @@ void make_best_corrs_graphics(std::vector< std::vector< std::vector< float > > >
     std::fclose(file);
 }
 
+void make_diffusion_file(const char* file_name, std::vector< std::vector< float > > D)
+{
+    FILE *file;
+    file = std::fopen(file_name, "w+");
+    for (int i = 0; i < D.size(); i++) {
+        std::fprintf(file, "diffusion for tau= %d\n", i);
+        for (int j = 0; j < D[i].size(); j++) {
+            std::fprintf(file, "%f ", D[i][j]);
+        }
+        std::fprintf(file, "\n");
+    }
+    std::fclose(file);
+}
+
 bool mysortfunc (std::vector< int > a, std::vector< int > b) {
     return (a.size() > b.size());
 }
@@ -351,6 +365,30 @@ void graph_back_bone_evaluation(std::vector< std::vector < std::pair< int, int >
     }
 }
 
+gmx::RVec evaluate_com(std::vector< RVec > frame) {
+    RVec temp;
+    temp[0] = 0;
+    temp[1] = 0;
+    temp[2] = 0;
+    for (int i = 0; i < frame.size(); i++) {
+        temp += frame[i];
+    }
+    temp[0] /= frame.size();
+    temp[1] /= frame.size();
+    temp[2] /= frame.size();
+    return temp;
+}
+
+void evaluate_diffusion(std::vector< std::vector< RVec > > trj, std::vector< std::vector< float > > &D/*, int max_frame_depth*/) {
+    D.resize(trj.size() * 0.9 /*max_frame_depth*/);
+    for (int i = 1; i < trj.size() * 0.9 /*max_frame_depth*/; i++) {
+        D[i].resize(trj.size() - 1 - i /*trj.size() - 1 - max_frame_depth*/);
+        for (int j = 0; j < trj.size() - 1 - i /*trj.size() - 1 - max_frame_depth*/; j++) {
+            D[i][j] = (evaluate_com(trj[j]) - evaluate_com(trj[j + i])).norm2() / (2 * i);
+        }
+    }
+}
+
 /*! \brief
  * \ingroup module_trajectoryanalysis
  */
@@ -508,6 +546,10 @@ SpaceTimeCorr::finishAnalysis(int nframes)
 
     std::cout << "Correlation's evaluation - end\n" << "graph evaluation: start\n";
 
+
+
+
+
     graph_calculation(graph, sub_graph, sub_graph_rbr, trajectory, basic_frame, crltns, crl_border, eff_rad, k);
 
     std::cout << "graph evaluation: end\n" << "routs evaluation: start\n";
@@ -521,6 +563,12 @@ SpaceTimeCorr::finishAnalysis(int nframes)
     make_best_corrs_graphics(crltns, rout_new, index, (OutPutName + "_routs_graphics.txt").c_str());
     std::cout << "corelation routs' pairs' graphics printed\n";
 
+
+    std::cout << "extra params\n";
+    std::vector< std::vector< float > > diffusion;
+    evaluate_diffusion(trajectory, diffusion);
+    make_diffusion_file((OutPutName + "_diffusion.txt").c_str(), diffusion);
+
     std::cout << "Finish Analysis - end\n\n";
 }