simple corr

+ reference frame for corrs

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 3f8c5343c4a571346552247d313ea5dd21ced5b0..a825c44c69f1da4521328420864d3095fb505772 100644 (file)
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,9 +1,9 @@
 set(GROMACS_CXX_FLAGS "${GROMACS_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
 set( CMAKE_CXX_FLAGS "-fopenmp") #-fsanitize=thread -O2 -g
-add_executable(spacetimecorr spacetimecorr.cpp)
+add_executable(simplecorr spacetimecorr.cpp)
 include_directories(
         ${GROMACS_INCLUDE_DIRS}
         ${CMAKE_SOURCE_DIR}/include)
-set_target_properties(spacetimecorr PROPERTIES
+set_target_properties(simplecorr PROPERTIES
        COMPILE_FLAGS "${GROMACS_CXX_FLAGS}")
-target_link_libraries(spacetimecorr ${GROMACS_LIBRARIES})
+target_link_libraries(simplecorr ${GROMACS_LIBRARIES})

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index c9fe008bfe5ad0699615ba3712ff9d1f89f5f902..1f0dafb701f543632b5eedfc37ca2c19ecf04711 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -51,6 +51,7 @@
 #include <gromacs/trajectoryanalysis.h>
 #include <gromacs/utility/smalloc.h>
 #include <gromacs/math/do_fit.h>
+#include <gromacs/trajectoryanalysis/topologyinformation.h>
 
 using namespace gmx;
 
@@ -62,10 +63,10 @@ void make_correlation_matrix_file(std::vector< std::vector< std::vector< float >
     file = std::fopen(file_name, "w+");
     for (int i = start; i < correlations.size(); i++) {
         if (correlations.size() - start > 1) {
-            std::fprintf(file, "correlation between 'n' and 'n + %d' frames\n", i);
+            std::fprintf(file, "correlation frame '%d'\n", i);
         }
-        if (i % 50 == 0) {
-            std::cout << "correlation between 'n' and 'n + " << i << "' frames\n";
+        if (i % 100 == 0) {
+            std::cout << "correlation " << i << " done\n";
         }
         for (int j = 0; j < correlations[i].size(); j++) {
             for (int f = 0; f < correlations[i][j].size(); f++) {
@@ -90,15 +91,13 @@ void make_correlation_pairs_file(std::vector< std::vector< std::vector< float >
             }
             std::fprintf(file, "\n");
         }
-        if (i % 10 == 0) {
-            std::cout << "correlations in row " << i << " out of " << correlations.front().size() << " completed\n";
-        }
+        std::cout << "correlations in row " << i << " out of " << correlations.front().size() << " completed\n";
     }
     std::fclose(file);
 }
 
-void correlation_evaluation(std::vector< std::vector< RVec > > traj, int window_size, std::vector< std::vector< std::vector< float > > > &crl) {
-    crl.resize(traj.size() - window_size + 1);
+void correlation_evaluation(std::vector< RVec > ref, std::vector< std::vector< RVec > > traj, int window_size, std::vector< std::vector< std::vector< float > > > &crl) {
+    crl.resize(traj.size() - window_size);
     for (int i = 0; i < crl.size(); i++) {
         crl[i].resize(traj.front().size());
         for (int j = 0; j < crl[i].size(); j++) {
@@ -121,39 +120,36 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int window_
     //{
     #pragma omp parallel for schedule(dynamic)
     for (int i = 0; i < crl.size(); i++) {
-        std::vector< RVec > medx, medy;
+        //std::vector< RVec > medx, medy;
         RVec temp1, temp2;
-        float tmp;
-        medx.resize(traj.front().size(), temp_zero);
-        medy.resize(traj.front().size(), temp_zero);
+        //medx.resize(traj.front().size(), temp_zero);
+        //medy.resize(traj.front().size(), temp_zero);
             //#pragma omp for schedule(dynamic)
-        for (int j = i; j < i + window_size - 1; j++) {
+        /*for (int j = i; j < i + window_size; j++) {
             for (int f = 0; f < traj.front().size(); f++) {
-                medx[f] += (traj[i][f] - traj[j][f]);
-                medy[f] += (traj[i][f] - traj[j][f]);
+                /*medx[f] += (traj[i][f] - traj[j][f]);
+                medy[f] += (traj[i][f] - traj[j][f]);*/
+                /*medx[f] += traj[j][f];
+                medy[f] += traj[j][f];
             }
-        }
+        }*/
             //#pragma omp barrier
             //#pragma omp for schedule(dynamic)
-        for (int j = 0; j < traj.front().size(); j++) {
-            tmp = traj.size() - 1;
-            tmp -= i;
-            tmp = 1 / tmp;
-                //tmp = 1 / (traj.size() - 1 - i);
-            medx[j] *= tmp;
-            medy[j] *= tmp;
-        }
+        /*for (int j = 0; j < traj.front().size(); j++) {
+            medx[j] /= window_size;
+            medy[j] /= window_size;
+        }*/
             //#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 = i; j < i + window_size - 1; j++) {
+        for (int j = i; j < i + window_size; 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[i][f1] - traj[j][f1] - medx[f1];
-                    temp2 = traj[i][f2] - traj[j][f2] - medy[f2];
+                    temp1 = /*traj[i][f1] - */traj[j][f1] - /*medx[f1]*/ ref[f1];
+                    temp2 = /*traj[i][f2] - */traj[j][f2] - /*medy[f2]*/ ref[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]);
@@ -169,8 +165,8 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int window_
             }
         }
             //#pragma omp barrier
-        medx.resize(0);
-        medy.resize(0);
+        //medx.resize(0);
+        //medy.resize(0);
         std::cout << i << " corr done\n";
     }
     #pragma omp barrier
@@ -217,6 +213,7 @@ class SimpleCorr : public TrajectoryAnalysisModule
         SelectionList                                               sel_;
 
         std::vector< std::vector< RVec > >                          trajectory;
+        std::vector< RVec >                                         reference;
 
         std::vector< int >                                          index;
         int                                                         frames              = 0;
@@ -259,6 +256,7 @@ SimpleCorr::initOptions(IOptionsContainer          *options,
                             .description("<your name here> + <local file tag>.txt"));
     // Control input settings
     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
+    settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
     settings->setPBC(true);
 }
 
@@ -272,6 +270,13 @@ SimpleCorr::initAnalysis(const TrajectoryAnalysisSettings &settings,
         index.push_back(*ai);
     }
     trajectory.resize(0);
+
+    reference.resize(0);
+    if (top.hasFullTopology()) {
+        for (int i = 0; i < index.size(); i++) {
+            reference.push_back(top.x().at(index[i]));
+        }
+    }
 }
 
 void
@@ -289,6 +294,13 @@ SimpleCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings
 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.md_npt.no-sol.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/SLca' -Window 100000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/window100'
 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/trp-cage.md_npt.no-sol.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/SLca' -Window 100000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/window100'
 
+// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/800-900k-stFit.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/CA1-20.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/SLca' -Window 10000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/window10new'
+// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/900-1000k-stFit.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/CA1-20.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/SLca' -Window 10000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/w-test9001000'
+// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/800-1000k-stFit.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/CA1-20.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/SLca' -Window 100000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/w100-test8001000'
+
+// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.tpr' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/900-1000k-stFit.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/SLca' -Window 10000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/w10-test9001000'
+
+
 void
 SimpleCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                       TrajectoryAnalysisModuleData *pdata)
@@ -307,7 +319,7 @@ SimpleCorr::finishAnalysis(int nframes)
     std::vector< std::vector< std::vector< float > > > crltns;
 
     std::cout << "\nCorrelation's evaluation - start\n";
-    correlation_evaluation(trajectory, W_size, crltns);
+    correlation_evaluation(reference, trajectory, W_size, crltns);
     std::cout << "\nCorrelation's evaluation - end\n";
 
     std::cout << "corelation matrix - start printing\n";