project(gromacs-rcore CXX)
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
if (NOT CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel." FORCE)
endif()
set(GROMACS_SUFFIX ${GMX_SUFFIX})
endif()
-find_package(GROMACS 2016 REQUIRED)
+if (GMX_OPENMP)
+ find_package(OpenMP REQUIRED)
+else()
+ find_package(OpenMP)
+endif()
+
+find_package(GROMACS 2019 REQUIRED)
gromacs_check_double(GMX_DOUBLE)
gromacs_check_compiler(CXX)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \file
+ * \brief
+ * Declares gmx::TopologyInformation.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \inlibraryapi
+ * \ingroup module_trajectoryanalysis
+ */
+#ifndef GMX_TRAJECTORYANALYSIS_TOPOLOGYINFORMATION_H
+#define GMX_TRAJECTORYANALYSIS_TOPOLOGYINFORMATION_H
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "gromacs/math/vectypes.h"
+#include "gromacs/topology/atoms.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/classhelpers.h"
+
+//! Forward declaration
+typedef struct gmx_rmpbc* gmx_rmpbc_t;
+
+namespace gmx
+{
+
+template<typename T>
+class ArrayRef;
+
+class TopologyInformation;
+class TrajectoryAnalysisRunnerCommon;
+
+/*! \libinternal
+ * \brief Topology information available to a trajectory analysis module.
+ *
+ * This class is used to provide topology information to trajectory
+ * analysis modules and to manage memory for them. Having a single
+ * wrapper object instead of passing each item separately makes
+ * TrajectoryAnalysisModule interface simpler, and also reduces the
+ * need to change existing code if additional information is added.
+ *
+ * It is intended that eventually most clients of this class will be
+ * analysis tools ported to the new analysis framework, but we will
+ * use this infrastructure also from the legacy analysis tools during
+ * the transition period. That will make it easier to put those tools
+ * under tests, and eventually port them.
+ *
+ * Methods in this class do not throw if not explicitly stated.
+ *
+ * The main data content is constant once loaded, but some content is
+ * constructed only when required (e.g. atoms_ and
+ * expandedTopology_). Their data members are mutable, so that the
+ * lazy construction idiom works properly. Some clients wish to modify
+ * the t_atoms, so there is an efficient mechanism for them to get a
+ * copy they can modify without disturbing this class. (The
+ * implementation releases the cached lazily constructed atoms_, but
+ * from the point of view of the caller, that is a copy.) The getters
+ * and copy functions are const for callers, which correctly expresses
+ * that the topology information is not being changed, merely copied
+ * and presented in a different form.
+ *
+ * \ingroup module_trajectoryanalysis
+ */
+class TopologyInformation
+{
+public:
+ //! Returns true if a topology file was loaded.
+ bool hasTopology() const { return hasLoadedMtop_; }
+ //! Returns true if a full topology file was loaded.
+ bool hasFullTopology() const { return bTop_; }
+ /*! \brief Builder function to fill the contents of
+ * TopologyInformation in \c topInfo from \c filename.
+ *
+ * Different tools require, might need, would benefit from, or
+ * do not need topology information. This functions implements
+ * the two-phase construction that is currently needed to
+ * support that.
+ *
+ * Any coordinate or run input file format will work, but the
+ * kind of data available from the getter methods afterwards
+ * will vary. For example, the mtop() available after reading
+ * a plain structure file will have a single molecule block and
+ * molecule type, regardless of contents.
+ *
+ * After reading, this object can return many kinds of primary
+ * and derived data structures to its caller.
+ *
+ * \todo This should throw upon error but currently does
+ * not. */
+ void fillFromInputFile(const std::string& filename);
+ /*! \brief Returns the loaded topology, or nullptr if not loaded. */
+ gmx_mtop_t* mtop() const { return mtop_.get(); }
+ //! Returns the loaded topology fully expanded, or nullptr if no topology is available.
+ const gmx_localtop_t* expandedTopology() const;
+ /*! \brief Returns a read-only handle to the fully expanded
+ * atom data arrays, which might be valid but empty if no
+ * topology is available. */
+ const t_atoms* atoms() const;
+ /*! \brief Copies the fully expanded atom data arrays, which
+ * might be valid but empty if no topology is available. */
+ AtomsDataPtr copyAtoms() const;
+ //! Returns the ePBC field from the topology.
+ int ePBC() const { return ePBC_; }
+ /*! \brief
+ * Gets the configuration positions from the topology file.
+ *
+ * If TrajectoryAnalysisSettings::efUseTopX has not been specified,
+ * this method should not be called.
+ *
+ * \throws APIError if topology position coordinates are not available
+ */
+ ArrayRef<const RVec> x() const;
+ /*! \brief
+ * Gets the configuration velocities from the topology file.
+ *
+ * If TrajectoryAnalysisSettings::efUseTopV has not been specified,
+ * this method should not be called.
+ *
+ * \throws APIError if topology velocity coordinates are not available
+ */
+ ArrayRef<const RVec> v() const;
+ /*! \brief
+ * Gets the configuration box from the topology file.
+ *
+ * \param[out] box Box size from the topology file, must not be nullptr.
+ */
+ void getBox(matrix box) const;
+ /*! \brief Returns a name for the topology.
+ *
+ * If a full topology was read from a a file, returns the name
+ * it contained, otherwise the empty string. */
+ const char* name() const;
+
+ TopologyInformation();
+ ~TopologyInformation();
+
+private:
+ //! The topology structure, or nullptr if no topology loaded.
+ std::unique_ptr<gmx_mtop_t> mtop_;
+ //! Whether a topology has been loaded.
+ bool hasLoadedMtop_;
+ //! The fully expanded topology structure, nullptr if not yet constructed.
+ mutable ExpandedTopologyPtr expandedTopology_;
+ //! The fully expanded atoms data structure, nullptr if not yet constructed.
+ mutable AtomsDataPtr atoms_;
+ //! true if full tpx file was loaded, false otherwise.
+ bool bTop_;
+ //! Position coordinates from the topology (can be nullptr).
+ std::vector<RVec> xtop_;
+ //! Velocity coordinates from the topology (can be nullptr).
+ std::vector<RVec> vtop_;
+ //! The box loaded from the topology file.
+ matrix boxtop_{};
+ //! The ePBC field loaded from the topology file.
+ int ePBC_;
+
+ // TODO This type is probably movable if we need that.
+ GMX_DISALLOW_COPY_AND_ASSIGN(TopologyInformation);
+
+ /*! \brief
+ * Needed to initialize the data.
+ */
+ friend class TrajectoryAnalysisRunnerCommon;
+};
+
+//! Convenience overload useful for implementing legacy tools.
+gmx_rmpbc_t gmx_rmpbc_init(const gmx::TopologyInformation& topInfo);
+
+} // namespace gmx
+
+#endif
-add_executable(rcore rcore.cpp)
+add_executable(rcore rcore.cpp newfit.cpp)
include_directories(
${GROMACS_INCLUDE_DIRS}
- ${CMAKE_SOURCE_DIR}/include)
+ ${CMAKE_SOURCE_DIR}/include
+ ${CMAKE_SOURCE_DIR}/src)
set_target_properties(rcore PROPERTIES
COMPILE_FLAGS "${GROMACS_CXX_FLAGS}")
target_link_libraries(rcore ${GROMACS_LIBRARIES})
--- /dev/null
+#include "newfit.h"
+
+// вычисление модуля расстояния между двумя точками, с учётом геометрического преобразования
+double F (const DVec &ai, const DVec &biR, const DVec &angl) {
+ const double cosA {cos(angl[0])}, cosB {cos(angl[1])}, cosC {cos(angl[2])}, sinA {sin(angl[0])}, sinB {sin(angl[1])}, sinC {sin(angl[2])};
+ return sqrt( pow(ai[0] + biR[1] * (cosA * sinC - cosC * sinA * sinB) - biR[2] * (sinA * sinC + cosA * cosC * sinB) - cosB * cosC * biR[0], 2) +
+ pow(ai[1] - biR[1] * (cosA * cosC + sinA * sinB * sinC) + biR[2] * (cosC * sinA - cosA * sinB * sinC) - cosB * sinC * biR[0], 2) +
+ pow(ai[2] + sinB * biR[0] - cosA * cosB * biR[2] - cosB * sinA * biR[1], 2) );
+}
+
+template< typename T1, typename T2 >
+void checkIfNan(T1 &a, const T2 &b) {
+ if (!std::isnan(b)) {
+ a += b;
+ } else {
+ a += 0.;
+ }
+}
+
+// вычисление функции F и её частных производных
+void searchF0xyzabc (long double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, const RVec &ai, const RVec &biR, const DVec &angl) {
+ const double cosA {cos(angl[0])}, cosB {cos(angl[1])}, cosC {cos(angl[2])}, sinA {sin(angl[0])}, sinB {sin(angl[1])}, sinC {sin(angl[2])};
+ double t01 {pow(ai[0] + biR[1] * (cosA * sinC - cosC * sinA * sinB) - biR[2] * (sinA * sinC + cosA * cosC * sinB) - cosB * cosC * biR[0], 2)};
+ double t02 {pow(ai[1] - biR[1] * (cosA * cosC + sinA * sinB * sinC) + biR[2] * (cosC * sinA - cosA * sinB * sinC) - cosB * sinC * biR[0], 2)};
+ double t03 {pow(ai[2] + sinB * biR[0] - cosA * cosB * biR[2] - cosB * sinA * biR[1], 2)};
+ double t04 {sqrt(t01 + t02 + t03)};
+ double t05 {(ai[0] + biR[1] * (cosA * sinC - cosC * sinA * sinB) - biR[2] * (sinA * sinC + cosA * cosC * sinB) - cosB * cosC * biR[0])};
+ double t06 {(ai[2] + sinB * biR[0] - cosA * cosB * biR[2] - cosB * sinA * biR[1])};
+ double t07 {(ai[1] - biR[1] * (cosA * cosC + sinA * sinB * sinC) + biR[2] * (cosC * sinA - cosA * sinB * sinC) - cosB * sinC * biR[0])};
+
+ checkIfNan(F, t04);
+ checkIfNan(Fx, -(2 * cosB * cosC * t05 - 2 * sinB * t06 + 2 * cosB * sinC * t07) / (2 * t04));
+ checkIfNan(Fy, -(2 * (cosA * cosC + sinA * sinB * sinC) * t07 - 2 * (cosA * sinC - cosC * sinA * sinB) * t05 + 2 * cosB * sinA * t06) / (2 * t04));
+ checkIfNan(Fz, -(2 * (sinA * sinC + cosA * cosC * sinB) * t05 - 2 * (cosC * sinA - cosA * sinB * sinC) * t07 + 2 * cosA * cosB * t06) / (2 * t04));
+ checkIfNan(Fa, -(2 * (cosA * cosB * biR[1] - cosB * sinA * biR[2]) * t06 -
+ 2 * (biR[1] * (cosC * sinA - cosA * sinB * sinC) + biR[2] * (cosA * cosC + sinA * sinB * sinC)) * t07 +
+ 2 * (biR[1] * (sinA * sinC + cosA * cosC * sinB) + biR[2] * (cosA * sinC - cosC * sinA * sinB)) * t05) / (2 * t04));
+ checkIfNan(Fb, -(2 * (cosA * cosB * sinC * biR[2] - sinB * sinC * biR[0] + cosB * sinA * sinC * biR[1]) * t07 +
+ 2 * (cosA * cosB * cosC * biR[2] - cosC * sinB * biR[0] + cosB * cosC * sinA * biR[1]) * t05 -
+ 2 * (cosB * biR[0] + sinA * sinB * biR[1] + cosA * sinB * biR[2]) * t06) / (2 * t04));
+ checkIfNan(Fc, (2 * (biR[1] * (cosA * cosC + sinA * sinB * sinC) - biR[2] * (cosC * sinA - cosA * sinB * sinC) + cosB * sinC * biR[0]) * t05 -
+ 2 * (biR[2] * (sinA * sinC + cosA * cosC * sinB) - biR[1] * (cosA * sinC - cosC * sinA * sinB) + cosB * cosC * biR[0]) * t07) / (2 * t04));
+}
+
+// применения геометрического преобразования: смещение на вектор (x, y, z) и повороты на градусы A, B, C относительно осей (x, y, z)
+//
+// А порядок поворотов какой?
+void ApplyFit (std::vector< RVec > &b, const DVec &R, const DVec &angl) {
+ DVec t(0., 0., 0.);
+ const double cosA {cos(angl[0])}, cosB {cos(angl[1])}, cosC {cos(angl[2])}, sinA {sin(angl[0])}, sinB {sin(angl[1])}, sinC {sin(angl[2])};
+ for (auto &i : b) {
+ t = i.toDVec();
+ i[0] = static_cast< float >((t[2] + R[2]) * (sinA * sinC + cosA * cosC * sinB) - (t[1] + R[1]) * (cosA * sinC - cosC * sinA * sinB) + cosB * cosC * (t[0] + R[0]));
+ i[1] = static_cast< float >((t[1] + R[1]) * (cosA * cosC + sinA * sinB * sinC) - (t[2] + R[2]) * (cosC * sinA - cosA * sinB * sinC) + cosB * sinC * (t[0] + R[0]));
+ i[2] = static_cast< float >(cosA * cosB * (t[2] + R[2]) - sinB * (t[0] + R[0]) + cosB * sinA * (t[1] + R[1]));
+ }
+}
+
+// подсчёт геометрических центров множеств a и b на основе пар pairs
+void CalcMid (const std::vector< RVec > &a, const std::vector< RVec > &b, DVec &midA, DVec &midB, const std::vector< std::pair< size_t, size_t > > &pairs) {
+ for (const auto &i : pairs) {
+ midA += a[i.first].toDVec();
+ midB += b[i.second].toDVec();
+ }
+ midA[0] /= pairs.size();
+ midA[1] /= pairs.size();
+ midA[2] /= pairs.size();
+
+ midB[0] /= pairs.size();
+ midB[1] /= pairs.size();
+ midB[2] /= pairs.size();
+}
+
+// функция фитирования фрейма b на фрейм a на основе пар pairs с "точностью" FtCnst
+void MyFitNew (const std::vector< RVec > &a, std::vector< RVec > &b, const std::vector< std::pair< size_t, size_t > > &pairs, long double FtCnst) {
+ double fx {0.}, fy {0.}, fz {0.}, fa {0.}, fb {0.}, fc {0.};
+ long double f1 {0.}, f2 {0.}, l {1};
+ DVec tempA(0., 0., 0.), tempB(0., 0., 0.), tempC(0., 0., 0.), tempD(0., 0., 0.);
+ CalcMid(a, b, tempA, tempB, pairs);
+ tempA -= tempB;
+ // поиск стартового отклонения множеств
+ for (const auto &i : pairs) {
+ f1 += F(a[i.first].toDVec(), b[i.second].toDVec() + tempA, tempC);
+ }
+ if (FtCnst == 0) {
+ FtCnst = 0.000'001; // experimental
+ }
+ if (f1 == 0) {
+ ApplyFit(b, tempA, tempC);
+ return;
+ } else {
+ // поиск оптимального геометрического преобразования методом градиентного спуска
+ while (f1 < f2 || f1 - f2 > FtCnst) {
+ f1 = 0.; f2 = 0.; fx = 0.; fy = 0.; fz = 0.; fa = 0.; fb = 0.; fc = 0.; l = 1.;
+ // поиск производных и отклонения при заданных параметрах
+ for (const auto &i : pairs) {
+ searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc, a[i.first], b[i.second] + tempA.toRVec(), tempC);
+ }
+ while (true) {
+ f2 = 0;
+ // поиск отклонения при новых параметрах
+ tempB[0] = tempA[0] - l * fx;
+ tempB[1] = tempA[1] - l * fy;
+ tempB[2] = tempA[2] - l * fz;
+ tempD[0] = tempC[0] - l * fa;
+ tempD[1] = tempC[1] - l * fb;
+ tempD[2] = tempC[2] - l * fc;
+ for (const auto &i : pairs) {
+ f2 += F(a[i.first].toDVec(), b[i.second].toDVec() + tempB, tempD);
+ }
+ if (f2 < f1) {
+ tempA = tempB;
+ tempC = tempD;
+ break;
+ } else {
+ // по факту существует минимальное значение переменной типа double
+ // в случае его достижения дважды останавливаем цикл т.к. упёрлись в предел допустимой точности
+ // таким образом гарантируем выход из цикла
+ if (l == DBL_MIN || l == 0.) { //DBL_TRUE_MIN
+ f2 = f1;
+ break;
+ }
+ l *= 0.1;
+ }
+ }
+ }
+ ApplyFit(b, tempA, tempC);
+ }
+}
--- /dev/null
+#ifndef NEWFIT_H
+#define NEWFIT_H
+
+#include <gromacs/trajectoryanalysis.h>
+#include <gromacs/trajectoryanalysis/topologyinformation.h>
+//#include "/home/titov_ai/Develop/gromacs-original/src/gromacs/trajectoryanalysis/topologyinformation.h"
+#include <gromacs/math/vectypes.h>
+#include <gromacs/math/vec.h>
+
+#include <iostream>
+#include <cfloat>
+#include <omp.h>
+#include <vector>
+#include <cstdio>
+
+using gmx::RVec;
+using gmx::DVec;
+using namespace gmx;
+
+// вычисление модуля расстояния между двумя точками, с учётом геометрического преобразования
+double F (const DVec &ai, const DVec &bi_plusR, const DVec &angl);
+
+template< typename T1, typename T2 > void checkIfNan(T1 &a, const T2 &b);
+
+// вычисление функции F и её частных производных
+void searchF0xyzabc (long double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, const RVec &ai, const RVec &biR, const DVec &angl);
+
+// применения геометрического преобразования: смещение на вектор (x, y, z) и повороты на градусы A, B, C относительно осей (x, y, z)
+void ApplyFit (std::vector< RVec > &b, const DVec &R, const DVec &angl);
+
+// подсчёт геометрических центров множеств a и b на основе пар pairs
+void CalcMid (const std::vector< RVec > &a, const std::vector< RVec > &b, DVec &midA, DVec &midB, const std::vector< std::pair< size_t, size_t > > &pairs);
+
+// функция фитирования фрейма b на фрейм a на основе пар pairs с "точностью" FtCnst
+void MyFitNew (const std::vector< RVec > &a, std::vector< RVec > &b, const std::vector< std::pair< size_t, size_t > > &pairs, long double FtCnst);
+
+#endif // NEWFIT_H
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
+
+#include <fstream>
#include <omp.h>
#include <iostream>
+#include <iomanip>
-#include <gromacs/trajectoryanalysis.h>
-#include <gromacs/math/do_fit.h>
-#include <gromacs/utility/smalloc.h>
+#include <newfit.h>
-using namespace gmx;
+#include <gromacs/trajectoryanalysis/topologyinformation.h>
/*! \brief
* Template class to serve as a basis for user analysis tools.
class AnalysisTemplate : public TrajectoryAnalysisModule
{
public:
- AnalysisTemplate();
virtual void initOptions(IOptionsContainer *options,
TrajectoryAnalysisSettings *settings);
virtual void writeOutput();
private:
- class ModuleData;
-
- std::string fnNdx_;
- double cutoff_;
- Selection refsel_;
-
- AnalysisNeighborhood nb_;
-
- AnalysisData data_;
- AnalysisDataAverageModulePointer avem_;
-
-
- Selection selec;
- std::vector< std::vector < RVec > > trajectory;
- std::vector< double > noise;
- std::vector< int > index;
- double epsi = 0.1; // selectable
- int etalon_frame = 0; // selectable
- int frames = 0;
+ //class ModuleData;
+
+ std::string outputFName; // selectable
+ Selection sel_;
+ std::vector< RVec > reference;
+ std::vector< std::vector < RVec > > trajectory;
+ std::vector< int > index;
+ double epsi {0.25}; // selectable
+ double prt {0.66}; // selectable
};
-AnalysisTemplate::AnalysisTemplate()
- : cutoff_(0.0)
-{
- registerAnalysisDataset(&data_, "avedist");
-}
-
void
AnalysisTemplate::initOptions(IOptionsContainer *options,
TrajectoryAnalysisSettings *settings)
settings->setHelpText(desc);
// Add option for selecting a subset of atoms
options->addOption(SelectionOption("select")
- .store(&selec).required()
- .description("Atoms that are considered as part of the excluded volume"));
+ .store(&sel_).required()
+ .description("Atoms that are considered as part of the excluded volume"));
// Add option for output file name
options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
- .store(&fnNdx_).defaultBasename("rcore")
- .description("Index file from the rcore"));
+ .store(&outputFName).defaultBasename("rcore")
+ .description("Index file from the rcore"));
// Add option for epsi constant
- options->addOption(DoubleOption("epsilon")
- .store(&epsi)
- .description("thermal vibrations' constant"));
- // Add option for etalon_frame constant
- options->addOption(gmx::IntegerOption("etalon_frame")
- .store(&etalon_frame)
- .description("basic frame to base evaluation on"));
+ options->addOption(DoubleOption("eps")
+ .store(&epsi)
+ .description("thermal vibrations' constant"));
+ // Add option for particion constant
+ options->addOption(DoubleOption("eps")
+ .store(&epsi)
+ .description("thermal vibrations' constant"));
}
void
AnalysisTemplate::initAnalysis(const TrajectoryAnalysisSettings &settings,
- const TopologyInformation & /*top*/)
+ const TopologyInformation &top)
{
index.resize(0);
- ConstArrayRef< int > atomind = selec.atomIndices();
- for (ConstArrayRef< int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++)
- {
+ for (ArrayRef< const int >::iterator ai {sel_.atomIndices().begin()}; ai < sel_.atomIndices().end(); ++ai) {
index.push_back(*ai);
}
+
+ reference.resize(0);
+ if (top.hasFullTopology()) {
+ for (size_t i {0}; i < index.size(); ++i) {
+ reference.push_back(top.x().at(index[i]));
+ }
+ }
+
+ trajectory.resize(0);
}
void
-AnalysisTemplate::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+AnalysisTemplate::analyzeFrame(int frnr,
+ const t_trxframe &fr,
+ t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
{
- trajectory.resize(frames + 1);
- trajectory[frames].resize(index.size());
- for (int i = 0; i < index.size(); i++)
+ trajectory.resize(trajectory.size() + 1);
+ trajectory.back().resize(0);
+ for (size_t i {0}; i < index.size(); ++i)
{
- trajectory[frames][i] = fr.x[index[i]];
+ trajectory.back().push_back(fr.x[index[i]]);
}
- frames++;
}
void
-AnalysisTemplate::finishAnalysis(int /*nframes*/)
+AnalysisTemplate::finishAnalysis(int nframes)
{
- real *w_rls;
- rvec *x1, *x2, temp;
- std::vector< bool > flag;
- double noise_mid = 9999;
- int count = 0;
- FILE *fpNdx_;
-
- snew(x1, index.size());
- snew(x2, index.size());
- snew(w_rls, index.size());
- flag.resize(index.size(), true);
- noise.resize(index.size(), 0);
-
- while (noise_mid > epsi)
- {
- noise_mid = 0;
- count = 0;
- noise.assign(index.size(), 0);
+ std::vector< std::pair< size_t, size_t > > fitPairs;
+ fitPairs.resize(0);
+ for (size_t i {0}; i < index.size(); ++i) {
+ fitPairs.push_back(std::make_pair(i, i));
+ }
- for (int i = 0; i < index.size(); i++)
- {
- if (flag[i])
- {
- w_rls[i] = 1;
- }
- else
- {
- w_rls[i] = 0;
- }
- }
+ std::vector< std::pair< size_t, double > > noise;
- for (int i = 0; i < frames; i++)
- {
- for (int j = 0; j < trajectory[i].size(); j++)
- {
- copy_rvec(trajectory[i][j].as_vec(), x2[j]);
- copy_rvec(trajectory[etalon_frame][j].as_vec(), x1[j]);
- }
- reset_x(index.size(), NULL, index.size(), NULL, x1, w_rls);
- reset_x(index.size(), NULL, index.size(), NULL, x2, w_rls);
- do_fit(index.size(), w_rls, x1, x2);
- if (i % 100 == 0)
- {
- std::cout << i << " / " << frames << "\n";
- }
- for (int j = 0; j < index.size(); j++)
- {
- rvec_sub(x1[j], x2[j], temp);
- noise[j] += norm(temp);
- }
+ size_t old = index.size(), nw = 0;
+ while (old != nw) {
+ noise.resize(0);
+ for (size_t i {0}; i < index.size(); ++i) {
+ noise.push_back(std::make_pair(i, 0.));
}
-
- for (int i = 0; i < index.size(); i++)
- {
- noise[i] /= frames;
+ for (size_t i {0}; i < trajectory.size(); ++i) {
+ MyFitNew(reference, trajectory[i], fitPairs, 0.000'001);
}
-
- for (int i = 0; i < index.size(); i++)
- {
- if (flag[i])
- {
- count++;
- noise_mid += noise[i];
+ for (size_t i {0}; i < trajectory.size(); ++i) {
+ for (size_t j {0}; j < trajectory[i].size(); ++j) {
+ noise[j].second += (reference[j] - trajectory[i][j]).norm();
}
}
- noise_mid /= count;
-
- for (int i = 0; i < index.size(); i++)
- {
- if (noise[i] > noise_mid)
- {
- flag[i] = false;
- }
+ for (size_t i {0}; i < noise.size(); ++i) {
+ noise[i].second /= trajectory.size();
}
- }
-
- int write_count = 0;
- fpNdx_ = std::fopen(fnNdx_.c_str(), "w+");
- std::fprintf(fpNdx_, "[ rcore ]\n");
- for (int i = 0; i < noise.size(); i++)
- {
- if (noise[i] <= epsi)
- {
- write_count++;
- if (write_count > 20) {
- write_count -= 20;
- std::fprintf(fpNdx_, "\n");
+ std::sort(noise.begin(), noise.end(), [](const std::pair< size_t, double > a, const std::pair< size_t, double > b){ return a.second < b.second;});
+ if (size_t temp {fitPairs.size() * prt}; noise[temp].second > epsi) {
+ fitPairs.resize(0);
+ for (size_t i {0}; i <= temp; ++i) {
+ fitPairs.push_back(std::make_pair(noise[i].first, noise[i].first));
+ }
+ old = nw;
+ nw = temp;
+ } else {
+ old = nw;
+ std::ofstream file(outputFName, std::ofstream::app);
+ file << "[ rcore ]\n";
+ for (size_t i {0}; i < temp; ++i) {
+ file << std::setw(7) << noise[i].first;
+ if (temp % 7 == 0) {
+ file << "\n";
+ }
}
- std::fprintf(fpNdx_, "%5d ", index[i] + 1);
}
}
- std::fprintf(fpNdx_,"\n");
- std::fclose(fpNdx_);
- sfree(x1);
- sfree(x2);
- sfree(w_rls);
}
void
biod.pnpi.spb.ru / alexxy / gromacs-rcore.git / commitdiff