From 0e5232b6699a8bc3b74cb2b5d45ede509eefecfd Mon Sep 17 00:00:00 2001
From: Anatoly Titov <toluk@omrb.pnpi.spb.ru>
Date: Wed, 19 Dec 2018 13:59:13 +0300
Subject: [PATCH] minor fixes

---
 src/domains.cpp | 83 ++++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 71 insertions(+), 12 deletions(-)

diff --git a/src/domains.cpp b/src/domains.cpp
index 31709a7..6ecb70b 100644
--- a/src/domains.cpp
+++ b/src/domains.cpp
@@ -49,6 +49,8 @@
 #include <gromacs/utility/smalloc.h>
 #include <gromacs/math/do_fit.h>
 
+#include "fittingn.cpp"
+
 #define MAX_NTRICVEC 12
 
 using namespace gmx;
@@ -95,6 +97,26 @@ void update_graph(std::vector< std::vector< node > > &ugwi_graph, rvec *ugwi_x,
     }
 }
 
+void update_graph_v2(std::vector< std::vector< node > > &ugwi_graph, std::vector < RVec > ugwi_x, long double ugwi_epsi) {
+    rvec ugwi_temp;
+    int ugwi_for = ugwi_graph.size();
+    for (int i = 0; i < ugwi_for; i++) {
+        for (int j = i; j < ugwi_for; j++) {
+            rvec_sub(ugwi_x[i].as_vec(), ugwi_x[j].as_vec(), ugwi_temp);
+            rvec_dec(ugwi_temp, ugwi_graph[i][j].r.as_vec());
+            if (norm(ugwi_temp) <= ugwi_epsi) {
+                if (i == j) {
+                    ugwi_graph[i][j].n++;
+                }
+                else {
+                    ugwi_graph[i][j].n++;
+                    ugwi_graph[j][i].n++;
+                }
+            }
+        }
+    }
+}
+
 void check_domains(long double cd_delta, int cd_frames, std::vector< std::vector< std::vector< node > > > &cd_graph) {
     int cd_for1 = cd_graph.size(), cd_for2 = cd_graph[1].size();
     for (int k = 0; k < cd_for1; k++) {
@@ -212,7 +234,9 @@ void print_domains(std::vector< std::vector< int > > pd_domains, std::vector< in
                 std::fprintf(fpNdx_, "\n");
             }
             std::fprintf(fpNdx_, "%5d ", index[pd_domains[i][j]] + 1);
+            std::printf("%5d ", index[pd_domains[i][j]] + 1);
         }
+        std::printf("\n\n");
         std::fprintf(fpNdx_,"\n\n");
     }
     std::fprintf(fpNdx_,"\n");
@@ -270,6 +294,7 @@ class Domains : public TrajectoryAnalysisModule
         int                                                         domain_min_size     = 5; // selectable
 
         real                                                        **w_rls;
+        std::vector< std::vector< std::pair< int, int > > >         w_rls2;
         int                                                         bone;
         double                                                      delta               = 0.90; // selectable
         double                                                      epsi                = 0.15; // selectable
@@ -329,23 +354,14 @@ void
 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
                       const TopologyInformation        &top)
 {
-    domains_ePBC = top.ePBC();
 }
 
 void
 Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
                              const t_trxframe                       &fr)
 {
-    t_pbc  pbc;
-    t_pbc *ppbc = settings.hasPBC() ? &pbc : NULL;
-    matrix boxx;
-    copy_mat(fr.box, boxx);
-    if (ppbc != NULL) {
-        set_pbc(ppbc, domains_ePBC, boxx);
-    }
-    ConstArrayRef< int > atomind  = selec.atomIndices();
     index.resize(0);
-    for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
+    for (ArrayRef< const int >::iterator ai = selec.atomIndices().begin(); (ai < selec.atomIndices().end()); ai++) {
         index.push_back(*ai);
     }
     trajectory.resize(2);
@@ -355,9 +371,13 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
         trajectory[0][i] = fr.x[index[i]];
     }
 
+
+
+
     bone = index.size() - domain_min_size + 1;
     graph.resize(bone);
     snew(w_rls, bone);
+    w_rls2.resize(bone);
     for (int i = 0; i < bone; i++) {
         snew(w_rls[i], index.size());
         for (int j = 0; j < index.size(); j++) {
@@ -367,6 +387,10 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
                 w_rls[i][j] = 0;
             }
         }
+        w_rls2[i].resize(0);
+        for (int j = 0; j < domain_min_size; j++) {
+            w_rls2[i].push_back(std::make_pair(j + i, j + i));
+        }
         rvec *etalon;
         snew(etalon, index.size());
         for (int j = 0; j < index.size(); j++) {
@@ -386,13 +410,14 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
     for (int i = 0; i < index.size(); i++) {
         trajectory[1][i] = fr.x[index[i]];
     }
+    std::vector < RVec > temp;
     frames++;
 
     #pragma omp parallel
     {
         #pragma omp for schedule(dynamic)
         for (int j = 0; j < bone; j++) {
-            rvec *etalon, *traj;
+            /*rvec *etalon, *traj;
             snew(etalon, index.size());
             for (int k = 0; k < index.size(); k++) {
                 copy_rvec(trajectory[0][k].as_vec(), etalon[k]);
@@ -406,7 +431,11 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
             do_fit(index.size(), w_rls[j], etalon, traj);
             update_graph(graph[j], traj, epsi);
             sfree(etalon);
-            sfree(traj);
+            sfree(traj);*/
+
+            std::vector < RVec > temp = trajectory[1];
+            MyFitNew(trajectory[0], temp, w_rls2[j]);
+            update_graph_v2(graph[j], temp, epsi);
         }
     }
     std::cout << "frame: " << frames << " analyzed\n";
@@ -460,3 +489,33 @@ main(int argc, char *argv[])
 {
     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-- 
2.22.0