- added reference choice
authorAnatoly <Titov_AI@pnpi.nrcki.ru>
Wed, 27 Mar 2019 13:44:16 +0000 (16:44 +0300)
committerAnatoly <Titov_AI@pnpi.nrcki.ru>
Wed, 27 Mar 2019 13:44:16 +0000 (16:44 +0300)
- structures are reseted to speed up fitting

src/fitng.cpp

index 835e1a1b98e527ca7b30bad9627f9d57035332bf..f4e41de6f222e52b1ae6fd33084e3997adca3971 100644 (file)
@@ -53,6 +53,7 @@
 #include <gromacs/utility/smalloc.h>
 #include <gromacs/math/do_fit.h>
 #include <gromacs/fileio/trxio.h>
+#include <gromacs/trajectoryanalysis/topologyinformation.h>
 
 using namespace gmx;
 
@@ -257,11 +258,22 @@ void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &mi
     return b.norm();
 }
 
+void reset_coord(std::vector< RVec > &frame, RVec move) {
+    for (int i = 0; i < frame.size(); i++) {
+        frame[i] -= move;
+    }
+}
+
 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs, double FtCnst) {
     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
     RVec ma, mb;
     CalcMid(a, b, ma, mb, pairs);
+
+    reset_coord(a, ma);
+    reset_coord(b, mb);
+
     rvec_dec(ma, mb);
+
     double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
     for (int i = 0; i < pairs.size(); i++) {
         f1 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
@@ -457,6 +469,8 @@ class Fitng : public TrajectoryAnalysisModule
         const int                                                   *ind;
         std::vector< std::vector< std::pair< int, int > > >         domains_index_number;
         std::vector< std::vector < RVec > >                         trajectory;
+        std::vector< RVec >                                         reference;
+
         std::vector< t_trxframe >                                   saved_traj;
         //Selection                                                   selec;
         std::vector< std::vector< std::pair< int, int > > >         pairs;
@@ -497,7 +511,7 @@ Fitng::initOptions(IOptionsContainer          *options,
                             .store(&FitConst)
                             .description("Fitting untill diff <= FitConst, by default == 0.00001"));
     // Add option for output pdb traj (meh edition)
-    options->addOption(StringOption("pdb_out")
+    options->addOption(StringOption("trj_out")
                             .store(&OutPutTrjName)
                             .description("transformed trajectory"));
     // Add option for trajectory Fit transformation
@@ -506,6 +520,7 @@ Fitng::initOptions(IOptionsContainer          *options,
                             .description("1 - fitting all -> all | 2 - keep only domains, all -> all | 3 - keep only domains, domain -> domain | 4 - extend domains to full coverage, domain -> domain"));
     // Control input settings
     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
+    settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
     settings->setPBC(true);
 }
 
@@ -519,6 +534,15 @@ Fitng::initOptions(IOptionsContainer          *options,
 
 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/test.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/testfit.xtc' -method 1
 
+// -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/local_test_trjconvfit_time.xtc' -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.pdb' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/local_test_prefit_time.xtc' -method 1 -e 500 -FitConst 0.000001
+// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test_Done.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1
+
+// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.md_npt.tpr' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test_Done.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -trj_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1 -h
+
+
+// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.md_npt.tpr' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test1k.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -trj_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1 -h
+
+
 void
 Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
 {
@@ -527,12 +551,19 @@ Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyIn
         index.push_back(*ai);
     }
 
+    reference.resize(0);
+    if (top.hasFullTopology()) {
+        for (int i = 0; i < index.size(); i++) {
+            reference.push_back(top.x().at(index[i]));
+        }
+    }
+    //std::cout << reference.size() << "\n important \n";
+
     saved_traj.resize(0);
 
     if (sel_.size() > 1 && method >= 2) {
         domain_reading(sel_, index, domains_index_number);
     }
-
 }
 
 void
@@ -541,10 +572,10 @@ Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_t
     trajectory.resize(0);
     std::vector < RVec > a;
     a.resize(0);
-    trajectory.push_back(a);
-    for (int i = 0; i < index.size(); i++) {
+    trajectory.push_back(reference);
+    /*for (int i = 0; i < index.size(); i++) {
         trajectory[0].push_back(fr.x[index[i]]);
-    }
+    }*/
     trajectory.push_back(a);
     switch (method) {
         case 1: