First blood

Signed-off-by: Anatoly Titov <toluk@omrb.pnpi.spb.ru>

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a948f068b156585573a11614a48c17efe1b01e56..62ec199fddd244bc8411da29f7045cd95d668f58 100644 (file)
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -31,9 +31,7 @@ endif()
 find_package(GROMACS 2016 REQUIRED)
 gromacs_check_double(GMX_DOUBLE)
 gromacs_check_compiler(CXX)
-include_directories(
-       ${GROMACS_INCLUDE_DIRS}
-       ${CMAKE_SOURCE_DIR}/include)
+
 add_definitions(${GROMACS_DEFINITIONS})
 
 # Use static linking on MSVC

diff --git a/include/new_fit.h b/include/new_fit.h
index 07148934acb8bf72d1fdeb5c4cb2ad2640fa7549..085220f8fd0b3503bf08dc2ea9bb5d817cc155a6 100644 (file)
--- a/include/new_fit.h
+++ b/include/new_fit.h
@@ -1,31 +1,12 @@
 #ifndef NEW_FIT_H
 #define NEW_FIT_H
 
-#include <iostream>
-#include <chrono>
-#include <omp.h>
-#include <thread>
-
-#include <gromacs/analysisdata/modules/plot.h>
-#include <gromacs/options/ioptionscontainer.h>
-#include <gromacs/pbcutil/pbc.h>
-#include <gromacs/math/functions.h>
-#include <gromacs/selection/selectionoption.h>
-#include <gromacs/utility/smalloc.h>
-#include <gromacs/trajectory/trajectoryframe.h>
-#include <gromacs/trajectoryanalysis/analysissettings.h>
 
 #include "gromacs/math/vectypes.h"
 
 using gmx::RVec;
 
-#ifdef __cplusplus
-extern "C" {
-#endif
 
-#ifdef __cplusplus
-}
-#endif
 
         void implement_fit(int if_size, std::vector< RVec > &if_frame, long double if_x, long double if_y, long double if_z, long double if_A, long double if_B, long double if_C);
 

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 9ff25d5f2a09c1a116e96be2af5f0e99be280efc..daa0132eacf705d58686c5e6196d08fee14163ba 100644 (file)
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,4 +1,7 @@
 add_executable(rcore  new_fit.cpp rcore.cpp)
+include_directories(
+        ${GROMACS_INCLUDE_DIRS}
+        ${CMAKE_SOURCE_DIR}/include)
 set_target_properties(rcore PROPERTIES
                       COMPILE_FLAGS "${GROMACS_CXX_FLAGS}")
 target_link_libraries(rcore ${GROMACS_LIBRARIES})

diff --git a/src/new_fit.cpp b/src/new_fit.cpp
index ce764eb9cff821acb978cd3007d47201bb9e63d3..cb3132942dda0be1210e5e501c54f5428a30e4ba 100644 (file)
--- a/src/new_fit.cpp
+++ b/src/new_fit.cpp
@@ -1,4 +1,16 @@
-#include "new_fit.h"
+#include <iostream>
+#include <chrono>
+#include <omp.h>
+#include <thread>
+
+#include <gromacs/analysisdata/modules/plot.h>
+#include <gromacs/options/ioptionscontainer.h>
+#include <gromacs/pbcutil/pbc.h>
+#include <gromacs/math/functions.h>
+#include <gromacs/selection/selectionoption.h>
+#include <gromacs/utility/smalloc.h>
+#include <gromacs/trajectory/trajectoryframe.h>
+#include <gromacs/trajectoryanalysis/analysissettings.h>
 
 using gmx::RVec;
 

diff --git a/src/rcore.cpp b/src/rcore.cpp
index 5b82cd5fdd1ce5b4a156de7da49f269d31ebd55b..d67df70e1a455f9bbcd8f5626fc220b162dd3a79 100644 (file)
--- a/src/rcore.cpp
+++ b/src/rcore.cpp
@@ -34,8 +34,13 @@
  */
 #include <string>
 #include <vector>
+#include <iostream>
+#include <omp.h>
+
+
 
 #include <gromacs/trajectoryanalysis.h>
+#include "new_fit.h"
 
 using namespace gmx;
 
@@ -64,12 +69,23 @@ class AnalysisTemplate : public TrajectoryAnalysisModule
         std::string                      fnDist_;
         double                           cutoff_;
         Selection                        refsel_;
-        SelectionList                    sel_;
 
         AnalysisNeighborhood             nb_;
 
         AnalysisData                     data_;
         AnalysisDataAverageModulePointer avem_;
+        
+        
+        Selection                                           selec;
+        std::vector< std::pair< int, int > >                fitting_pairs;
+        std::vector< std::vector < RVec > >                 trajectory;
+        std::vector< int >                                  noise;
+        std::vector< int >                                  index;
+        const long double                                   epsi = 0.15;
+        const long double                                   fitting_prec = 0.0000001;
+        int                                                 frames = 0;
+        
+        
 };
 
 
@@ -102,46 +118,37 @@ AnalysisTemplate::initOptions(IOptionsContainer          *options,
 
     settings->setHelpText(desc);
 
-    options->addOption(FileNameOption("o")
+    /*options->addOption(FileNameOption("o")
                            .filetype(eftPlot).outputFile()
                            .store(&fnDist_).defaultBasename("avedist")
                            .description("Average distances from reference group"));
 
     options->addOption(SelectionOption("reference")
                            .store(&refsel_).required()
-                           .description("Reference group to calculate distances from"));
+                           .description("Reference group to calculate distances from"));*/
     options->addOption(SelectionOption("select")
-                           .storeVector(&sel_).required().multiValue()
-                           .description("Groups to calculate distances to"));
+                           .store(&selec).required()
+                           .description("Atoms that are considered as part of the excluded volume"));
 
-    options->addOption(DoubleOption("cutoff").store(&cutoff_)
+   /* options->addOption(DoubleOption("cutoff").store(&cutoff_)
                            .description("Cutoff for distance calculation (0 = no cutoff)"));
 
-    settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
+    settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);*/
 }
 
+//domains -s '/home/toluk/Рабочий стол/reca_rd_2008/reca_rd.md.non-sol.tpr' -f '/home/toluk/Рабочий стол/reca_rd_2008/reca_rd.md.non-sol.xtc' -select 'name CA' -dt 100
 
 void
 AnalysisTemplate::initAnalysis(const TrajectoryAnalysisSettings &settings,
                                const TopologyInformation         & /*top*/)
 {
-    nb_.setCutoff(cutoff_);
-
-    data_.setColumnCount(0, sel_.size());
-
-    avem_.reset(new AnalysisDataAverageModule());
-    data_.addModule(avem_);
-
-    if (!fnDist_.empty())
-    {
-        AnalysisDataPlotModulePointer plotm(
-                new AnalysisDataPlotModule(settings.plotSettings()));
-        plotm->setFileName(fnDist_);
-        plotm->setTitle("Average distance");
-        plotm->setXAxisIsTime();
-        plotm->setYLabel("Distance (nm)");
-        data_.addModule(plotm);
+    std::cout << "select start\n";
+    index.resize(0);
+    ConstArrayRef< int > atomind = selec.atomIndices();
+    for (ConstArrayRef< int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
+        index.push_back(*ai);
     }
+    std::cout << "select finish\n";
 }
 
 
@@ -149,43 +156,57 @@ void
 AnalysisTemplate::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                                TrajectoryAnalysisModuleData *pdata)
 {
-    AnalysisDataHandle         dh     = pdata->dataHandle(data_);
-    const Selection           &refsel = pdata->parallelSelection(refsel_);
-
-    AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, refsel);
-    dh.startFrame(frnr, fr.time);
-    for (size_t g = 0; g < sel_.size(); ++g)
-    {
-        const Selection &sel   = pdata->parallelSelection(sel_[g]);
-        int              nr    = sel.posCount();
-        real             frave = 0.0;
-        for (int i = 0; i < nr; ++i)
-        {
-            SelectionPosition p = sel.position(i);
-            frave += nbsearch.minimumDistance(p.x());
-        }
-        frave /= nr;
-        dh.setPoint(g, frave);
+    std::cout << "trajectory start\n";
+    trajectory.resize(frames + 1);
+    trajectory[frames].resize(index.size());
+    for (int i = 0; i < index.size(); i++) {
+        trajectory[frames][i] = fr.x[index[i]];
     }
-    dh.finishFrame();
+    frames++;
+    std::cout << "trajectory finish\n";
 }
 
 
 void
 AnalysisTemplate::finishAnalysis(int /*nframes*/)
 {
+    std::cout << "analys start\n";
+    for (int i = 0; i < index.size(); i++) {
+        fitting_pairs.push_back(std::make_pair (i, i));
+    }
+        long double noise_mid = 9000;
+    while (noise_mid > epsi) {
+        for (int i = 1; i < frames; i++) {
+            new_fit(trajectory[0], trajectory[i], fitting_pairs, index.size(), fitting_prec);
+            for (int j = 0; j < fitting_pairs.size(); j++) {
+                noise[fitting_pairs[j].first] += norm(trajectory[0][fitting_pairs[j].first] - trajectory[i][fitting_pairs[j].first]);
+            }
+        }
+        noise_mid = 0;
+        for (int i = 0; i < fitting_pairs.size(); i++) {
+            noise[fitting_pairs[i].first] /= (frames - 1);
+            noise_mid += noise[fitting_pairs[i].first];
+        }
+        noise_mid /= fitting_pairs.size();
+        for (int i = 0; i < fitting_pairs.size(); i++) {
+            if (noise[fitting_pairs[i].first] > noise_mid) {
+                fitting_pairs.erase(fitting_pairs.begin() + i);
+                i--;
+            }
+        }
+    }
+    std::cout << "analys finish\n";
 }
 
 
 void
 AnalysisTemplate::writeOutput()
 {
-    // We print out the average of the mean distances for each group.
-    for (size_t g = 0; g < sel_.size(); ++g)
-    {
-        fprintf(stderr, "Average mean distance for '%s': %.3f nm\n",
-                sel_[g].name(), avem_->average(0, g));
+    std::cout << "output start\n";
+    for (int i = 0; i < fitting_pairs.size(); i++) {
+        std::cout << fitting_pairs[i].first << " ";
     }
+    std::cout << "output finish\n";
 }
 
 /*! \brief