сохранение нормальной траектории

но нужно сделать force copy как-то

diff --git a/src/fitng.cpp b/src/fitng.cpp
index 6c30431ed31bacd72ebb225ed8a371ef6d15220b..9b3e6db936fede08ca790c0a120b70f670c5ad9f 100644 (file)
--- a/src/fitng.cpp
+++ b/src/fitng.cpp
@@ -51,6 +51,7 @@
 #include <gromacs/pbcutil/pbc.h>
 #include <gromacs/utility/smalloc.h>
 #include <gromacs/math/do_fit.h>
+#include <gromacs/fileio/trxio.h>
 
 using namespace gmx;
 
@@ -245,6 +246,17 @@ void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &mi
     midB[2] /= pairs.size();
 }
 
+float FrameMidLength(std::vector< RVec > a) {
+    RVec b;
+    b[0] = 0;
+    b[1] = 0;
+    b[2] = 0;
+    for (int i = 0; i < a.size(); i++) {
+        b += a[i];
+    }
+    return b.norm();
+}
+
 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;
@@ -399,7 +411,7 @@ void TrajFitting(std::vector< std::vector < RVec > > &trj, double FC, std::vecto
             }
         }
         clone.resize(0);
-        std::cout << "frame " << i << " fitted\n";
+        //std::cout << "frame " << i << " fitted\n";
     }
 }
 
@@ -438,9 +450,12 @@ class Fitng : public TrajectoryAnalysisModule
         SelectionList                                               sel_;
 
         std::vector< int >                                          index;
+        const int                                                   *ind;
         std::vector< std::vector< std::pair< int, int > > >         domains_index_number;
         std::vector< std::vector < RVec > >                         trajectory;
+        std::vector< t_trxframe >                                   saved_traj;
         //Selection                                                   selec;
+        std::vector< std::vector< std::pair< int, int > > >         pairs;
         int                                                         frames = 0;
         int                                                         method = 1;
         double                                                      FitConst = 0;
@@ -489,22 +504,14 @@ Fitng::initOptions(IOptionsContainer          *options,
     settings->setPBC(true);
 }
 
-// Fitng -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
-// Fitng -s '/home/toluk/Develop/FitSamples/cube_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/cube.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_termal_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_termal.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_1st_fr_D100.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_D100.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_1st_fr_D10.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_D10.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/garmoshka_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/garmoshka.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/garmoshka_2nd_fr.pdb' -f '/home/toluk/Develop/FitSamples/garmoshka.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/rubber_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/rubber.pdb' -select 'name NA'
-// rubber_1st_fr
-
 // /home/toluk/Develop/samples/reca_test_ground/
 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL0' -n '/home/toluk/Develop/samples/reca_test_ground/t0.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/FullCAFit.pdb' -e 1000
 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsFit.pdb' -e 1000 -method 2
 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsOnlyFit.pdb' -e 1000 -method 3
 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsWExtraOnlyFit.pdb' -e 1000 -method 4
 
+// -s '/home/titov_ai/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/titov_ai/Develop/samples/reca_rd/reca_rd.mono.xtc' -sf '/home/titov_ai/Develop/samples/reca_rd/SLCa' -n '/home/titov_ai/Develop/samples/reca_rd/CA.ndx' -pdb_out '/home/titov_ai/Develop/samples/reca_rd/test.xtc' -method 1
+
 void
 Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
 {
@@ -512,7 +519,8 @@ Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyIn
     for (ArrayRef< const int >::iterator ai = sel_[0].atomIndices().begin(); (ai < sel_[0].atomIndices().end()); ai++) {
         index.push_back(*ai);
     }
-    trajectory.resize(0);
+
+    saved_traj.resize(0);
 
     if (sel_.size() > 1 && method >= 2) {
         domain_reading(sel_, index, domains_index_number);
@@ -529,18 +537,35 @@ void
 Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                       TrajectoryAnalysisModuleData *pdata)
 {
-    trajectory.resize(trajectory.size() + 1);
-    trajectory.back().resize(index.size());
-    for (int i = 0; i < index.size(); i++) {
-        trajectory.back()[i] = fr.x[index[i]];
+    t_trxframe c = fr;
+    saved_traj.push_back(c); //тут нужно force copy сделать
+
+    std::vector< RVec > a, b;
+    a.resize(index.size());
+    b.resize(index.size());
+    for (int j = 0; j < index.size(); j++) {
+        a[j] = saved_traj.back().x[index[j]];
     }
+    for (int j = 0; j < index.size(); j++) {
+        b[j] = saved_traj.front().x[index[j]];
+    }
+    std::cout << FrameMidLength(a) << " " << FrameMidLength(b) << " ";
 }
 
 
 void
 Fitng::finishAnalysis(int nframes)
 {
-    std::vector< std::vector< std::pair< int, int > > > pairs;
+    //std::vector < RVec > a;
+    //a.resize(index.size());
+    trajectory.resize(saved_traj.size()/*, a*/);
+    for (int i = 0; i < saved_traj.size(); i++) {
+        for (int j = 0; j < index.size(); j++) {
+            //trajectory[i][j] = saved_traj[i].x[index[j]];
+            trajectory[i].push_back(saved_traj[i].x[index[j]]);
+        }
+        std::cout << FrameMidLength(trajectory[i]) << " " << i << " ";
+    }
 
     switch (method) {
         case 1:
@@ -548,7 +573,6 @@ Fitng::finishAnalysis(int nframes)
             for (int i = 0; i < index.size(); i++) {
                 pairs[0].push_back(std::make_pair(i, i));
             }
-            TrajFitting(trajectory, FitConst, pairs);
             break;
         case 2:
             pairs.resize(1);
@@ -557,27 +581,27 @@ Fitng::finishAnalysis(int nframes)
                     pairs[0].push_back(std::make_pair(domains_index_number[i][j].second, domains_index_number[i][j].second));
                 }
             }
-            TrajFitting(trajectory, FitConst, pairs);
             break;
         case 3:
             MakeDomainFitPairs(pairs, domains_index_number);
-            TrajFitting(trajectory, FitConst, pairs);
             break;
         case 4:
             domain_expansion(index, domains_index_number, trajectory[0]);
             MakeDomainFitPairs(pairs, domains_index_number);
-            TrajFitting(trajectory, FitConst, pairs);
             break;
     }
 
-
-    printPDBtraj(OutPutTrjName.c_str(), trajectory, pairs, index);
+    TrajFitting(trajectory, FitConst, pairs);
 }
 
 void
 Fitng::writeOutput()
 {
-
+    //printPDBtraj(OutPutTrjName.c_str(), trajectory, pairs, index);
+    std::cout << saved_traj.size();
+    for (int i = 0; i < saved_traj.size(); i++) {
+        std::cout << write_trxframe_indexed(open_trx(OutPutTrjName.c_str(), "a+"), static_cast<const t_trxframe *>(&saved_traj[i]), index.size(), static_cast<const int *>(&index.front())/*saved_traj[i].index*/, NULL) << " ";
+    }
 }