added 2 new correlations

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

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index 5a787cf18637fbfa365a2a92f09b5cdd1c9b1d9d..6727257d88690b829c158e7b3e1ff782f0763718 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -442,16 +442,22 @@ Domains::finishAnalysis(int nframes)
 
     //frankenstein_trajectory = trajectory;
 
-    std::vector< std::vector< std::vector< double > > > crltns;
+    std::vector< std::vector< std::vector< double > > > crltns, crltns_scalars, crltns_mini;
     int k = 1000, m = 0;
     if (tau > 0) {
         k = tau;
     }
     crltns.resize(k + 1);
+    crltns_scalars.resize(k + 1);
+    crltns_mini.resize(k + 1);
     for (int i = 0; i < crltns.size(); i++) {
         crltns[i].resize(index.size());
+        crltns_scalars[i].resize(index.size());
+        crltns_mini[i].resize(index.size());
         for (int j = 0; j < crltns[i].size(); j++) {
             crltns[i][j].resize(index.size(), 0);
+            crltns_scalars[i][j].resize(index.size(), 0);
+            crltns_mini[i][j].resize(index.size(), 0);
         }
     }
 
@@ -462,49 +468,89 @@ Domains::finishAnalysis(int nframes)
         for (int i = m; i <= k; i += 1) {
             std::cout << " k = " << i << " ";
 
-            rvec *medx, *medy, temp_zero;
+            rvec *medx, *medy, temp_zero, *medx_mini, *medy_mini;
+            std::vector< double > medx_scalar, medy_scalar;
             snew(medx, index.size());
             snew(medy, index.size());
+            medx_scalar.resize(index.size(), 0);
+            medy_scalar.resize(index.size(), 0);
+            snew(medx_mini, index.size());
+            snew(medy_mini, 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]);
+                copy_rvec(temp_zero, medx_mini[i]);
+                copy_rvec(temp_zero, medy_mini[i]);
             }
 
-            for (int j = 0; j < frankenstein_trajectory.size() - i; j++) {
+            for (int j = 0; j < frankenstein_trajectory.size() - i - 1; j++) {
                 for (int f = 0; f < index.size(); f++) {
                     rvec temp;
 
                     rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j][f], temp);
                     rvec_inc(medx[f], temp);
+
+                    medx_scalar[f] += norm(temp);
+
                     rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + i][f], temp);
                     rvec_inc(medy[f], temp);
+
+                    medy_scalar[f] += norm(temp);
+
+                    rvec_sub(frankenstein_trajectory[j][f], frankenstein_trajectory[j + 1][f], temp);
+                    rvec_inc(medx_mini[f], temp);
+                    rvec_sub(frankenstein_trajectory[j + i][f], frankenstein_trajectory[j + i + 1][f], temp);
+                    rvec_inc(medy_mini[f], temp);
+
                 }
             }
 
             for (int j = 0; j < index.size(); j++) {
                 rvec temp;
-                real temp2 = frankenstein_trajectory.size();
+                real temp2 = frankenstein_trajectory.size() - 1;
                 temp2 -= i;
                 temp2 = 1 / temp2;
+
                 copy_rvec(medx[j], temp);
                 svmul(temp2, temp, medx[j]);
                 copy_rvec(medy[j], temp);
                 svmul(temp2, temp, medy[j]);
+
+                medx_scalar[j] *= temp2;
+                medy_scalar[j] *= temp2;
+
+                copy_rvec(medx_mini[j], temp);
+                svmul(temp2, temp, medx_mini[j]);
+                copy_rvec(medy_mini[j], temp);
+                svmul(temp2, temp, medy_mini[j]);
+
             }
 
-            std::vector< std::vector< double > > a, b, c;
+            std::vector< std::vector< double > > a, b, c, as, bs, cs, am, bm, cm;
             a.resize(index.size());
             b.resize(index.size());
             c.resize(index.size());
+            as.resize(index.size());
+            bs.resize(index.size());
+            cs.resize(index.size());
+            am.resize(index.size());
+            bm.resize(index.size());
+            cm.resize(index.size());
             for (int j = 0; j < index.size(); j++) {
                 a[j].resize(index.size(), 0);
                 b[j].resize(index.size(), 0);
                 c[j].resize(index.size(), 0);
+                as[j].resize(index.size(), 0);
+                bs[j].resize(index.size(), 0);
+                cs[j].resize(index.size(), 0);
+                am[j].resize(index.size(), 0);
+                bm[j].resize(index.size(), 0);
+                cm[j].resize(index.size(), 0);
             }
-            for (int j = 0; j < frankenstein_trajectory.size() - i; j++) {
+            for (int j = 0; j < frankenstein_trajectory.size() - i - 1; j++) {
                 for (int f1 = 0; f1 < index.size(); f1++) {
                     for (int f2 = 0; f2 < index.size(); f2++) {
                         rvec temp1, temp2;
@@ -515,22 +561,42 @@ Domains::finishAnalysis(int nframes)
                         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]);
+
+                        rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j][f1], temp1);
+                        rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + i][f2], temp2);
+                        as[f1][f2] += (norm(temp1) - medx_scalar[f1]) * (norm(temp2) - medy_scalar[f2]);
+                        bs[f1][f2] += (norm(temp1) - medx_scalar[f1]) * (norm(temp1) - medx_scalar[f1]);
+                        cs[f1][f2] += (norm(temp2) - medy_scalar[f2]) * (norm(temp2) - medy_scalar[f2]);
+
+                        rvec_sub(frankenstein_trajectory[j][f1], frankenstein_trajectory[j + 1][f1], temp1);
+                        rvec_dec(temp1, medx_mini[f1]);
+                        rvec_sub(frankenstein_trajectory[j + i][f2], frankenstein_trajectory[j + i + 1][f2], temp2);
+                        rvec_dec(temp2, medy_mini[f2]);
+                        am[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
+                        bm[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
+                        cm[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
                     }
                 }
             }
             for (int j = 0; j < index.size(); j++) {
                 for (int f = 0; f < index.size(); f++) {
                     crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
+                    crltns_scalars[i][j][f] = as[j][f] / (std::sqrt(bs[j][f] * cs[j][f]));
+                    crltns_mini[i][j][f] = am[j][f] / (std::sqrt(bm[j][f] * cm[j][f]));
                 }
             }
             sfree(medx);
             sfree(medy);
+            sfree(medx_mini);
+            sfree(medy_mini);
         }
     }
     std::cout << "\n" ;
 
     std::cout << "printing correlation's file\n" ;
-    make_correlation_file(crltns, "corrs.txt", m);
+    make_correlation_file(crltns, "corrs_vectors.txt", m);
+    make_correlation_file(crltns_scalars, "corrs_vectors_scalars.txt", m);
+    make_correlation_file(crltns_mini, "corrs_mini_vectors.txt", m);
     std::cout << "done\n" ;
 
     for (int i1 = 0; i1 < 100; i1++) {