fixes

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index ce952939cd45ea33443dce0da33fa09ba6c15bae..6e5d624c5bff4f4d57b18457caa6747a945e3f16 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -77,14 +77,12 @@ using gmx::RVec;
 
 void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const char* file_name)
 {
-    FILE               *file;
+    FILE *file;
     file = std::fopen(file_name, "w+");
     for (int i = 1; i < trajectory.size(); i++) {
         std::fprintf(file, "MODEL%9d\n", i);
         for (int j = 0; j < trajectory[i].size(); j++) {
             std::fprintf(file, "ATOM  %5d   CA LYS     1    %8.3f%8.3f%8.3f  1.00  0.00\n", j, trajectory[i][j][0] * 10, trajectory[i][j][1] * 10, trajectory[i][j][2] * 10);
-            //                    123456  123456789012345678901    8    6   4567890123456
-            //std::fprintf(file, "1234567890123456789012345678901234567890123456789012345678901234567890\n");
         }
         std::fprintf(file, "ENDMDL\n");
     }
@@ -108,10 +106,9 @@ void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const ch
 79 - 80        LString(2)      Charge on the atom.
 */
 
-void make_correlation_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start)
+void make_correlation_file(std::vector< std::vector< std::vector< double > > > correlations, const char* file_name, int start)
 {
-
-    FILE               *file;
+    FILE *file;
     file = std::fopen(file_name, "w+");
     std::cout << "start printing\n";
     for (int i = start; i < correlations.size(); i++) {
@@ -132,8 +129,7 @@ void make_correlation_file(std::vector< std::vector< std::vector< float > > > co
 
 void make_rout_file(std::vector< int > dist, std::vector< std::vector< int > > rout, const char* file_name)
 {
-
-    FILE               *file;
+    FILE *file;
     file = std::fopen(file_name, "w+");
     std::cout << "start printing\n";
     for (int i = 0; i < rout.size(); i++) {
@@ -419,39 +415,65 @@ void
 Domains::finishAnalysis(int nframes)
 {
     //make_pdb_trajectory(frankenstein_trajectory, "/home/toluk/Develop/samples/reca_rd/frank_with_duo_test.pdb");
-    std::vector< std::vector< std::vector< float > > > crltns;
-    int k = 1000, m = 1;
-    crltns.resize(k);
+    std::vector< std::vector< std::vector< double > > > crltns;
+    int k = 950, m = 949;
+    crltns.resize(k + 1);
     for (int i = 0; i < crltns.size(); i++) {
         crltns[i].resize(index.size());
         for (int j = 0; j < crltns[i].size(); j++) {
             crltns[i][j].resize(index.size(), 0);
         }
     }
-    #pragma omp parallel
-    {
-        #pragma omp for schedule(dynamic)
-        for (int i = m; i < k; i += 1) {
+
+    //#pragma omp parallel
+    //{
+        //#pragma omp for schedule(dynamic)
+        for (int i = m; i <= k; i += 1) {
             std::cout << " \n " << " k = " << i << " \n";
-            std::vector< float > medx, medy;
-            medx.resize(index.size(), 0);
-            medy.resize(index.size(), 0);
+            //std::vector< float > medx, medy;
+            //medx.resize(index.size(), 0);
+            //medy.resize(index.size(), 0);
+            rvec *medx, *medy, temp_zero;
+            snew(medx, index.size());
+            snew(medy, index.size());
+            temp_zero[0] = 0;
+            temp_zero[1] = 0;
+            temp_zero[2] = 0;
+            for (int i = 0; i < index.size(); i++) {
+                copy_rvec(temp_zero, medx[i]);
+                copy_rvec(temp_zero, medy[i]);
+            }
             for (int j = 0; j < frankenstein_trajectory.size() - k; j++) {
                 for (int f = 0; f < index.size(); f++) {
                     rvec temp;
-                    rvec_sub(frankenstein_trajectory[j][f], frankenstein_trajectory[j + 1][f], temp);
-                    //rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + 1][f], temp);
-                    medx[f] += norm(temp);
-                    rvec_sub(frankenstein_trajectory[j + k - 1][f], frankenstein_trajectory[j + k][f], temp);
-                    //rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + k][f], temp);
-                    medy[f] += norm(temp);
+                    //rvec_sub(frankenstein_trajectory[j][f], frankenstein_trajectory[j + 1][f], temp);
+                    rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j][f], temp);
+                    rvec_inc(medx[f], temp);
+                    //medx[f] += norm(temp);
+                    //rvec_sub(frankenstein_trajectory[j + k - 1][f], frankenstein_trajectory[j + k][f], temp);
+                    rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + k][f], temp);
+                    rvec_inc(medy[f], temp);
+                    //medy[f] += norm(temp);
                 }
             }
+
             for (int j = 0; j < index.size(); j++) {
-                medx[j] /= frankenstein_trajectory.size();
-                medy[j] /= frankenstein_trajectory.size();
+                //medx[j] /= frankenstein_trajectory.size();
+                rvec temp;
+                real temp2 = frankenstein_trajectory.size() - k;
+                temp2 = 1 / temp2;
+                copy_rvec(medx[j], temp);
+                svmul(temp2, temp, medx[j]);
+                //medy[j] /= frankenstein_trajectory.size();
+                copy_rvec(medy[j], temp);
+                svmul(temp2, temp, medy[j]);
             }
-            std::vector< std::vector< float > > a, b, c;
+
+            /*for (int j = 0; j < index.size(); j++) {
+                std::cout << medx[j][0] << " " << medx[j][1] << " " << medx[j][2] << " " << medy[j][0] << " " << medy[j][1] << " " << medy[j][2] << "\n";
+            }*/
+
+            std::vector< std::vector< double > > a, b, c;
             a.resize(index.size());
             b.resize(index.size());
             c.resize(index.size());
@@ -460,32 +482,52 @@ Domains::finishAnalysis(int nframes)
                 b[j].resize(index.size(), 0);
                 c[j].resize(index.size(), 0);
             }
-            for (int j = 1; j < frankenstein_trajectory.size() - k; j++) {
+            for (int j = 0; j < frankenstein_trajectory.size() - k; j++) {
                 for (int f1 = 0; f1 < index.size(); f1++) {
-                    for (int f2 = f1; f2 < index.size(); f2++) {
+                    for (int f2 = 0; f2 < index.size(); f2++) {
                         rvec temp1, temp2;
-                        rvec_sub(frankenstein_trajectory[j][f1], frankenstein_trajectory[j + 1][f1], temp1);
-                        //rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j + 1][f1], temp1);
-                        rvec_sub(frankenstein_trajectory[j + k - 1][f2], frankenstein_trajectory[j + k][f2], temp2);
-                        //rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + k][f2], temp2);
-                        a[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp2) - medy[f2]);
-                        b[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp1) - medx[f1]);
-                        c[f1][f2] += (norm(temp2) - medy[f2]) * (norm(temp2) - medy[f2]);
+                        //rvec_sub(frankenstein_trajectory[j][f1], frankenstein_trajectory[j + 1][f1], temp1);
+                        rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j][f1], temp1);
+                        rvec_dec(temp1, medx[f1]);
+                        //rvec_sub(frankenstein_trajectory[j + k - 1][f2], frankenstein_trajectory[j + k][f2], temp2);
+                        rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + k][f2], temp2);
+                        rvec_dec(temp2, medy[f2]);
+                        //a[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp2) - medy[f2]);
+                        //b[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp1) - medx[f1]);
+                        //c[f1][f2] += (norm(temp2) - medy[f2]) * (norm(temp2) - medy[f2]);
+                        a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
+                        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]);
+                        //std::cout << a[f1][f2] << " " << b[f1][f2] << " " << c[f1][f2] << "\n";
+
                     }
                 }
             }
             for (int j = 0; j < index.size(); j++) {
-                for (int f = j; f < index.size(); f++) {
-                    crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f]) * std::sqrt(c[j][f]));
+                for (int f = 0; f < index.size(); f++) {
+                    crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
                 }
             }
+            sfree(medx);
+            sfree(medy);
         }
-    }
+    //}
 //    std::cout << "\noutput file\n";
-//    make_correlation_file(crltns, "/home/toluk/Develop/samples/reca_rd/corrs2.txt", m);
+    make_correlation_file(crltns, "/home/toluk/Develop/samples/reca_rd/corrs2.txt", m);
 //    std::cout << "\n\n file done \n\n";
 
-    std::cout << "\n\n\n\n";
+
+
+
+
+
+
+
+
+
+
+
+    /*std::cout << "\n\n\n\n";
     for (int i2 = 0; i2 < 75; i2++) {
         long int count_here = 0;
         for (int i = 0; i < k; i++) {
@@ -499,7 +541,7 @@ Domains::finishAnalysis(int nframes)
         }
         std::cout << i2 << " " << count_here << "\n";
     }
-    std::cout << "\n\n\n\n";
+    std::cout << "\n\n\n\n";*/
 
 
 
@@ -535,11 +577,11 @@ Domains::finishAnalysis(int nframes)
     std::cout << "building routs\n" ;
 
 
-
+/*
     for (int i = 0; i < index.size(); i++) {
-        std::cout << dive_in(graph, i, 0, frames, 0);
+        std::cout << i << " " << dive_in(graph, i, 0, frames, 0) << " ";
     }
-
+*/
     /*while (flag) {
         flag = false;
         for (long int i = 0; i < rout_old.size(); i++) {