fixed diffusion

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index 98ef82d9ffd5bda8f853e00b7268172b8496c5a3..19a426d49da57e900c4c893d4b6c469eeedc94cb 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -130,16 +130,12 @@ 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)
+void make_diffusion_file(const char* file_name, 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::fprintf(file, "%f ", D[i]);
     }
     std::fclose(file);
 }
@@ -182,22 +178,26 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame
 
     std::vector< float > d;
     d.resize(traj.front().size(), 0);
+    //int nth = omp_get_thread_num();
 
     #pragma omp parallel shared(crl, temp_zero, d)
     {
-        #pragma omp for schedule(dynamic)
+
         for (int i = tauS; i <= tauE; i += 1) {
             std::vector< RVec > medx, medy;
             RVec temp1, temp2;
             float tmp;
             medx.resize(traj.front().size(), temp_zero);
             medy.resize(traj.front().size(), temp_zero);
+            #pragma omp for schedule(dynamic)
             for (int j = 0; j < traj.size() - i - 1; j++) {
                 for (int f = 0; f < traj.front().size(); f++) {
                     medx[f] += (traj[b_frame][f] - traj[j][f]);
                     medy[f] += (traj[b_frame][f] - traj[j + i][f]);
                 }
             }
+            #pragma omp barrier
+            #pragma omp for schedule(dynamic)
             for (int j = 0; j < traj.front().size(); j++) {
                 tmp = traj.size() - 1;
                 tmp -= i;
@@ -206,12 +206,14 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame
                 medx[j] *= tmp;
                 medy[j] *= tmp;
             }
+            #pragma omp barrier
             std::vector< std::vector< float > > a, b, c;
             a.resize(traj.front().size(), d);
             b.resize(traj.front().size(), d);
             c.resize(traj.front().size(), d);
             for (int j = 0; j < traj.size() - i - 1; j++) {
                 for (int f1 = 0; f1 < traj.front().size(); f1++) {
+                    #pragma omp for schedule(dynamic)
                     for (int f2 = 0; f2 < traj.front().size(); f2++) {
                         temp1 = traj[b_frame][f1] - traj[j][f1] - medx[f1];
                         temp2 = traj[b_frame][f2] - traj[j + i][f2] - medy[f2];
@@ -219,19 +221,23 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame
                         b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
                         c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
                     }
+                    #pragma omp barrier
                 }
             }
+
+            #pragma omp for schedule(dynamic)
             for (int j = 0; j < traj.front().size(); j++) {
                 for (int f = 0; f < traj.front().size(); f++) {
                     crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
                 }
             }
+            #pragma omp barrier
             medx.resize(0);
             medy.resize(0);
             std::cout << i << " corr done\n";
         }
     }
-    #pragma omp barrier
+
 }
 
 void graph_calculation(std::vector< std::vector< std::pair< float, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
@@ -379,13 +385,16 @@ gmx::RVec evaluate_com(std::vector< RVec > frame) {
     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*/);
+void evaluate_diffusion(std::vector< std::vector< RVec > > trj, std::vector< float > &D/*, int max_frame_depth*/) {
+    D.resize(trj.size() * 0.9 - 1/*max_frame_depth*/);
+    float temp;
     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*/);
+        temp = 0;
         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);
+            temp += (evaluate_com(trj[j]) - evaluate_com(trj[j + i])).norm2();
+            //D[i][j] = (evaluate_com(trj[j]) - evaluate_com(trj[j + i])).norm2() / (2 * i);
         }
+        D[i - 1] = temp / (trj.size() - 1 - i) / (2 * i * 0.000001);
     }
 }
 
@@ -565,7 +574,7 @@ SpaceTimeCorr::finishAnalysis(int nframes)
 
 
     std::cout << "extra params\n";
-    std::vector< std::vector< float > > diffusion;
+    std::vector< float > diffusion;
     evaluate_diffusion(trajectory, diffusion);
     make_diffusion_file((OutPutName + "_diffusion.txt").c_str(), diffusion);