From d2ab957721e853eec809419bfec04608b7d427c3 Mon Sep 17 00:00:00 2001
From: Anatoly Titov <toluk@omrb.pnpi.spb.ru>
Date: Thu, 7 Sep 2017 20:35:50 +0300
Subject: [PATCH] - added effective radius for corr search for atom's
 interaction - added a selection option for the constant, mentioned above -
 added output for entire corr function for the atom pairs which reach highest
 numbers along the time interval

---
 src/spacetimecorr.cpp | 48 ++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 47 insertions(+), 1 deletion(-)

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index ff9dd7d..9699791 100644
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -43,6 +43,7 @@
 #include <iostream>
 #include <algorithm>
 #include <omp.h>
+#include <set>
 
 #include <gromacs/trajectoryanalysis.h>
 #include <gromacs/utility/smalloc.h>
@@ -100,6 +101,23 @@ void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std
     std::fclose(file);
 }
 
+void make_best_corrs_graphics(std::vector< std::vector< std::vector< double > > > correlations,
+                              std::vector< std::pair< int, int > > corr_pairs,
+                              std::vector< int > indx,
+                              const char* file_name)
+{
+    FILE *file;
+    file = std::fopen(file_name, "w+");
+    for (int i = 0; i < corr_pairs.size(); i++) {
+        std::fprintf(file, "%3d %3d", indx[corr_pairs[i].first] + 1, indx[corr_pairs[i].second] + 1);
+        for (int j = 0; j < correlations.size(); j++) {
+            std::fprintf(file, "%3.2f ", correlations[j][corr_pairs[i].first][corr_pairs[i].second]);
+        }
+        std::fprintf(file, "\n");
+    }
+    std::fclose(file);
+}
+
 bool mysortfunc (std::vector< int > a, std::vector< int > b) { return (a.size() > b.size()); }
 
 /*! \brief
@@ -153,6 +171,7 @@ class Domains : public TrajectoryAnalysisModule
         int                                                         basic_frame         = 0;
         int                                                         tau                 = 0;
         double                                                      crl_border          = 0;
+        double                                                      eff_rad             = 1.5;
         real                                                        **w_rls;
 
         int                                                         domains_ePBC;
@@ -188,6 +207,10 @@ Domains::initOptions(IOptionsContainer          *options,
     options->addOption(DoubleOption("crl")
                             .store(&crl_border)
                             .description("make graph based on corrs > constant"));
+    // Add option for eff_rad constant
+    options->addOption(DoubleOption("ef_rad")
+                            .store(&eff_rad)
+                            .description("effective radius for atoms to evaluate corrs"));
     // Add option for selection list
     options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
                            .required().dynamicMask().multiValue()
@@ -452,6 +475,7 @@ Domains::finishAnalysis(int nframes)
     make_correlation_file(crltns, "corrs_vectors.txt", m);
     std::cout << "done\n" ;
 
+    int imax = 0;
     for (int i1 = 0; i1 < 100; i1++) {
         int i5 = 0;
         for (int i2 = 1; i2 < crltns.size(); i2++) {
@@ -463,6 +487,9 @@ Domains::finishAnalysis(int nframes)
                 }
             }
         }
+        if (i5 > 0) {
+            imax = i1;
+        }
         std::cout << i1 << " - " << i5 << " | ";
         if (i1 % 10 == 0) {
             std::cout << "\n" ;
@@ -470,13 +497,32 @@ Domains::finishAnalysis(int nframes)
     }
     std::cout << "\n";
 
+    std::set< std::pair< int, int > > crr_prs;
+    std::vector< std::pair< int, int > > corr_prs;
+    for (int i2 = 0; i2 < crltns.size(); i2++) {
+        for (int i3 = 0; i3 < crltns[i2].size(); i3++) {
+            for (int i4 = i3 + 1; i4 < crltns[i2][i3].size(); i4++) {
+                if (crltns[i2][i3][i4] >= (double)imax / 100) {
+                    crr_prs.insert(std::make_pair(i3, i4));
+                }
+            }
+        }
+    }
+    for (std::set< std::pair< int, int > >::iterator it = crr_prs.begin(); it != crr_prs.end(); it++) {
+        corr_prs.push_back(*it);
+    }
+    make_best_corrs_graphics(crltns, corr_prs, index, "best_graphics.txt");
+
     std::vector< std::vector< int > > graph;
     graph.resize(index.size());
+    rvec temp;
 
     for (int i = 1; i <= k; i++) {
         for (int j = 0; j < index.size(); j++) {
             for (int f = 0; f < index.size(); f++) {
-                if (std::abs(crltns[i][j][f]) >= crl_border) {
+                copy_rvec(frankenstein_trajectory[basic_frame][j], temp);
+                rvec_dec(temp, frankenstein_trajectory[basic_frame][f]);
+                if (std::abs(crltns[i][j][f]) >= crl_border && norm(temp) <= eff_rad) {
                     bool local_flag = true;
                     for (int f2 = 0; f2 < graph[j].size(); f2++) {
                         if (graph[j][f2] == f) {
-- 
2.22.0