--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2013,2014,2015, 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.
+ */
+/*! \internal \file
+ * \brief
+ * Implements gmx::analysismodules::Freevolume.
+ *
+ * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
+ * \ingroup module_trajectoryanalysis
+ */
+
+#include "domaintype.h"
+#include "newfit.h"
+
+#include <gromacs/trajectoryanalysis.h>
+#include <gromacs/trajectoryanalysis/topologyinformation.h>
+#include <gromacs/selection/nbsearch.h>
+
+#include <iostream>
+#include <fstream>
+//#include <chrono>
+#include <cfloat>
+#include <omp.h>
+#include <vector>
+#include <cstdio>
+
+// структура углов для одной краски
+struct colorLocalAngles {
+ gmx::RVec a1, a2;
+ gmx::RVec b1, b2;
+ gmx::RVec n1, n2;
+ double a12, b12, n12;
+ std::vector< double > betaAngles;
+};
+
+// хрен пойми почему, но захотелось рекурсивно сделать | можно сделать и через вайл
+long long returnBetaEnd(const std::string ¤tLine, long long pos) {
+ switch (currentLine[pos]) {
+ case 'B' :
+ returnBetaEnd(currentLine, pos + 1);
+ break;
+ case 'E' :
+ returnBetaEnd(currentLine, pos + 1);
+ break;
+ default :
+ return (pos - 1);
+ break;
+ }
+ return (pos - 1); // формально логически лишняя строчка, но Qt ругается
+}
+
+// функция парсинга одной строки
+void parseBetaListDATLine(const std::string ¤tLine, std::vector< std::vector< unsigned int > > &localInputBL) {
+ size_t equalCount = 0;
+ std::vector< unsigned int > a;
+ a.resize(0);
+ for (size_t i = 0; i < currentLine.size(); ++i) {
+ if (currentLine[i] == '=') { // подсчитываем число "пустых" символов
+ ++equalCount;
+ } else {
+ long long temp = returnBetaEnd(currentLine, i);
+ if (temp - i > 3) {
+ localInputBL.push_back(a);
+ for (size_t j = i; j <= temp; ++j) {
+ localInputBL.back().push_back(j - equalCount);
+ }
+ i = temp;
+ }
+ }
+ }
+}
+
+// функция нахождения бета-листов в структуре по файлу ДССП
+void betaListDigestion(const std::string &inputFile, std::vector< std::vector< std::vector< unsigned int > > > &inputBL) {
+ inputBL.resize(0);
+ std::ifstream file(inputFile);
+ std::string line;
+ getline(file, line); // считываем число в первой строке - кол-во осмысленных элементов в строках - нам не нужно
+ getline(file, line);
+ if (line.size() > 0) {
+ do {
+ inputBL.resize(inputBL.size() + 1);
+ parseBetaListDATLine(line, inputBL.back());
+ getline(file, line);
+ } while (line.size() > 0);
+ } else {
+ throw "DSSP DAT FILE IS EMPTY";
+ }
+ file.close();
+}
+
+// функция поиска RVec в кадре по имени->индексу
+gmx::RVec returnRVec(const std::vector< RVec > &frame, const std::vector< std::pair< std::string, size_t > > &colorsIndex, std::string toFind) {
+ for (auto &i : colorsIndex) {
+ if (i.first == toFind) {
+ return frame[i.second];
+ }
+ }
+ throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS";
+}
+
+// функция векторного произведения двух RVec'ов
+inline void RVecVecMultiply(gmx::RVec &n, const gmx::RVec &a, const gmx::RVec &b) {
+ n[0] = a[1] * b[2] - a[2] * b[1];
+ n[1] = a[2] * b[0] - a[0] * b[2];
+ n[2] = a[0] * b[1] - a[2] * b[0];
+}
+
+// поиск угла между двумя RVec'ами
+inline double RVecAngle(const gmx::RVec &a, const gmx::RVec &b) {
+ return std::acos((a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) / (a.norm() * b.norm()));
+}
+
+// вычисление внутренних углов в краске
+void colorsAnglesEvaluation(const std::vector< RVec > &frame, const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex, std::vector< colorLocalAngles > &colorStruct) {
+ colorStruct.resize(colorsIndex.size());
+ #pragma omp parallel for ordered schedule(dynamic)
+ for (size_t i = 0; i < colorsIndex.size(); ++i) {
+ colorStruct[i].betaAngles.resize(0);
+ // "левый" блок
+ colorStruct[i].a1 = returnRVec(frame, colorsIndex[i], "CAF") - returnRVec(frame, colorsIndex[i], "CAJ") + returnRVec(frame, colorsIndex[i], "CAK") - returnRVec(frame, colorsIndex[i], "CAO");
+ colorStruct[i].a1 /= colorStruct[i].a1.norm();
+ colorStruct[i].b1 = returnRVec(frame, colorsIndex[i], "CAO") - returnRVec(frame, colorsIndex[i], "CAJ") + returnRVec(frame, colorsIndex[i], "CAM") - returnRVec(frame, colorsIndex[i], "CAH") + returnRVec(frame, colorsIndex[i], "CAK") - returnRVec(frame, colorsIndex[i], "CAF");
+ colorStruct[i].b1 /= colorStruct[i].b1.norm();
+ RVecVecMultiply(colorStruct[i].n1, colorStruct[i].a1, colorStruct[i].b1);
+ colorStruct[i].n1 /= colorStruct[i].n1.norm();
+ // "правый" блок
+ colorStruct[i].a2 = returnRVec(frame, colorsIndex[i], "CAB") - returnRVec(frame, colorsIndex[i], "CBQ") + returnRVec(frame, colorsIndex[i], "CBP") - returnRVec(frame, colorsIndex[i], "CBL");
+ colorStruct[i].a2 /= colorStruct[i].a2.norm();
+ colorStruct[i].b2 = returnRVec(frame, colorsIndex[i], "CBL") - returnRVec(frame, colorsIndex[i], "CBQ") + returnRVec(frame, colorsIndex[i], "CBN") - returnRVec(frame, colorsIndex[i], "CBS") + returnRVec(frame, colorsIndex[i], "CBP") - returnRVec(frame, colorsIndex[i], "CAB");
+ colorStruct[i].b2 /= colorStruct[i].b2.norm();
+ RVecVecMultiply(colorStruct[i].n2, colorStruct[i].a2, colorStruct[i].b2);
+ colorStruct[i].n2 *= -1; // для "сонаправленности" векторов нормали
+ colorStruct[i].n2 /= colorStruct[i].n2.norm();
+ colorStruct[i].a12 = RVecAngle(colorStruct[i].a1, colorStruct[i].a2);
+ colorStruct[i].b12 = RVecAngle(colorStruct[i].b1, colorStruct[i].b2);
+ colorStruct[i].n12 = RVecAngle(colorStruct[i].n1, colorStruct[i].n2);
+ }
+ #pragma omp barrier
+}
+
+// вычисление направляющего вектора в бэта-листе
+inline void betaListsRVecsEvaluation(const std::vector< RVec > &frame, const std::vector< std::vector< unsigned int > > &inputBetaLists, std::vector< gmx::RVec > &temp) {
+ temp.resize(0);
+ gmx::RVec tempA;
+ for (const auto &i : inputBetaLists) {
+ tempA = frame[i[i.size() - 1]] + frame[i[i.size() - 2]] - frame[i[1]] - frame[i[0]];
+ tempA /= tempA.norm();
+ temp.push_back(tempA);
+ }
+}
+
+// определение близко ли к белку находится краска
+bool isNearPeptide(const t_trxframe &fr, const t_pbc *pbc, const std::vector< RVec > inputFrame,
+ gmx::AnalysisNeighborhood &nbhood, const std::vector< size_t > &inputIndex,
+ const std::vector< std::pair< std::string, size_t > > &inputColor) {
+ gmx::AnalysisNeighborhoodSearch nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
+ gmx::AnalysisNeighborhoodPair pair;
+ for (const auto &i : inputColor) {
+ while (nbsearch.startPairSearch(inputFrame[i.second].as_vec()).findNextPair(&pair)) {
+ for (const auto &j : inputIndex) {
+ if (pair.testIndex() == j) {
+ return true;
+ }
+ }
+ }
+ }
+ return false;
+}
+
+// поиск ближайших к краске бэта-листов
+inline void searchNearBetaLists(const t_trxframe &fr, const t_pbc *pbc, const std::vector< RVec > inputFrame,
+ gmx::AnalysisNeighborhood &nbhood,
+ const std::vector< std::vector< unsigned int > > inputBLists,
+ const std::vector< std::pair< std::string, size_t > > &inputColor,
+ std::vector< bool > &outputList) {
+ outputList.resize(0);
+ outputList.resize(inputBLists.size(), false);
+ gmx::AnalysisNeighborhoodSearch nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
+ gmx::AnalysisNeighborhoodPair pair;
+ for (const auto &i : inputColor) {
+ while (nbsearch.startPairSearch(inputFrame[i.second].as_vec()).findNextPair(&pair)) {
+ for (size_t j = 0; j < inputBLists.size(); ++j) {
+ for (size_t k = 0; k < inputBLists[j].size(); ++k) {
+ if (pair.testIndex() == inputBLists[j][k]) {
+ outputList[j] = true;
+ }
+ }
+ }
+ }
+ }
+
+}
+
+// определение углов между краской и направляющим вектором бэта-листа
+inline void computeAnglesColorsVsBeta(const gmx::RVec &inputBetaRVec, colorLocalAngles &colorFormation) {
+ colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a1));
+ colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b1));
+ colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n1));
+ colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a2));
+ colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b2));
+ colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n2));
+}
+
+// функция записи в файл значений углов для кадра
+void anglesFileDump(const int frameNum, const std::string &output, const std::vector< bool > &toPeptide, const std::vector< colorLocalAngles > &colorFormation) {
+ std::ofstream file(output);
+ int temp = 0;
+ std::vector< double > betaTemp;
+ betaTemp.resize(0);
+ file << "frame =" << std::setw(8) << frameNum << std::endl;
+ for (size_t i = 0; i < colorFormation.size(); ++i) {
+ file << "color #" << std::setw(3) << i;
+ if (toPeptide[i]) {
+ file << std::setw(4) << "yes";
+ } else {
+ file << std::setw(4) << "no";
+ }
+ file << std::setw(7) << std::setprecision(2) << colorFormation[i].a12 << colorFormation[i].b12 << colorFormation[i].n12;
+ temp = colorFormation[i].betaAngles.size() / 6;
+ if (temp == 0) {
+ file << std::setw(4) << 0;
+ } else {
+ betaTemp.resize(6, 0.); // magic number, meh
+ for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
+ betaTemp[j % 6] += colorFormation[i].betaAngles[j];
+ }
+ for (size_t j = 0; j < betaTemp.size(); ++j) {
+ betaTemp[j] /= temp;
+ }
+ file << std::setw(4) << temp;
+ for (size_t j = 0; j < betaTemp.size(); ++j) {
+ file << std::setw(7) << std::setprecision(2) << betaTemp[j];
+ }
+ for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
+ file << std::setw(7) << std::setprecision(2) << colorFormation[i].betaAngles[j];
+ }
+ }
+ }
+ file.close();
+}
+
+/*! \brief
+ * \ingroup module_trajectoryanalysis
+ */
+class colorVec : public TrajectoryAnalysisModule
+{
+ public:
+
+ colorVec();
+ virtual ~colorVec();
+
+ //! Set the options and setting
+ virtual void initOptions(IOptionsContainer *options,
+ TrajectoryAnalysisSettings *settings);
+
+ //! First routine called by the analysis framework
+ // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
+ virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
+ const TopologyInformation &top);
+
+ //! Call for each frame of the trajectory
+ // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
+ virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+ TrajectoryAnalysisModuleData *pdata);
+
+ //! Last routine called by the analysis framework
+ // virtual void finishAnalysis(t_pbc *pbc);
+ virtual void finishAnalysis(int nframes);
+
+ //! Routine to write output, that is additional over the built-in
+ virtual void writeOutput();
+
+ private:
+
+ SelectionList sel_;
+ std::string fnOut; // selectable
+ std::string fnBetaListsDat; // selectable
+ float effRad {0.8}; // selectable
+ std::vector< size_t > index;
+ std::vector< std::vector< std::pair< std::string, size_t > > > colorsNames;
+ std::vector< std::vector< std::vector< unsigned int > > > betaLists;
+ AnalysisNeighborhood nb_;
+
+ // Copy and assign disallowed by base.
+};
+
+colorVec::colorVec(): TrajectoryAnalysisModule()
+{
+}
+
+colorVec::~colorVec()
+{
+}
+
+/*
+ * ndx:
+ * [peptide]
+ * [color_{i}]
+ *
+ */
+
+void
+colorVec::initOptions(IOptionsContainer *options,
+ TrajectoryAnalysisSettings *settings)
+{
+ static const char *const desc[] = {
+ "[THISMODULE] to be done"
+ };
+ // Add the descriptive text (program help text) to the options
+ settings->setHelpText(desc);
+ // Add option for selection list
+ options->addOption(SelectionOption("sel")
+ .storeVector(&sel_)
+ .required().dynamicMask().multiValue()
+ .description("select residue and domains to form a rigid skeleton / -sf"));
+ // Add option for output file name
+ options->addOption(FileNameOption("dat").filetype(eftOptionFileType_NR).inputFile()
+ .store(&fnBetaListsDat).description("a file to make dynamic beta lists"));
+ // Add option for output file name
+ options->addOption(FileNameOption("out").filetype(eftOptionFileType_NR).outputFile()
+ .store(&fnOut).defaultBasename("colorAngles")
+ .description("Index file for the algorithm output."));
+ // Add option for effRad constant
+ options->addOption(FloatOption("efRad")
+ .store(&effRad)
+ .description("max distance from colors to peptide in nm to consider to be \"near\""));
+ // Control input settings
+ settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
+ settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
+ //settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
+ settings->setPBC(true);
+}
+
+void
+colorVec::initAnalysis(const TrajectoryAnalysisSettings &settings,
+ const TopologyInformation &top)
+{
+ // считывание индекса
+ index.resize(0);
+ for (ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ai++) {
+ index.push_back(static_cast< size_t >(*ai));
+ }
+ // считывание красок
+ colorsNames.resize(sel_.size() - 1);
+ for (size_t i = 1; i < sel_.size(); ++i) {
+ for (ArrayRef< const int >::iterator ai {sel_[i].atomIndices().begin()}; (ai < sel_[i].atomIndices().end()); ai++) {
+ colorsNames[i - 1].push_back(std::make_pair(*(top.atoms()->resinfo[top.atoms()->atom[*ai].resind].name), static_cast< size_t >(*ai)));
+ }
+ }
+ // разбор dat файла для создания бэта листов
+ betaListDigestion(fnBetaListsDat, betaLists);
+ // задание радиуса рассматриваемых соседей
+ nb_.setCutoff(effRad);
+ /*
+ * формально можно сделать:
+ * найти соотношение между именами для углов и индексными номерами
+ * заполнить std::vector< std::vector< size_t > > nameToIndex;
+ * и в последствии считать:
+ * a1 = frame[nameToIndex[0][1] + ... - ... - ...
+ * b1 a2 b2 по аналогии только для [@<-[1, 2, 3]]
+ *
+ */
+}
+
+void
+colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
+{
+ std::vector< RVec > trajectoryFrame;
+ trajectoryFrame.resize(0);
+ // считывания текущего фрейма траектории
+ for (size_t i {0}; i < index.size(); ++i) {
+ trajectoryFrame.push_back(fr.x[index[i]]);
+ }
+ // подсчёт углов в красках
+ std::vector< colorLocalAngles > colorFormation;
+ colorsAnglesEvaluation(trajectoryFrame, colorsNames, colorFormation);
+ // рассчёт положения относительно белка
+ /*
+ * формально можно совместить определение близости с белком и определение ближайших бэта-листов
+ * для этого думаю нужно как-то соотнести или сделать соотношение между индексом и бэта-листами (т.е. хранить что бы за О(1) делать)
+ *
+ */
+
+ std::vector< bool > colorsToPeptide;
+ colorsToPeptide.resize(0);
+ colorsToPeptide.resize(colorsNames.size(), false);
+ for (size_t i = 0; i < colorsNames.size(); ++i) {
+ colorsToPeptide[i] = isNearPeptide(fr, pbc, trajectoryFrame, nb_, index, colorsNames[i]);
+ }
+ // расчёт угла и среднего угла с ближайшими бета листами
+ std::vector< std::vector< bool > > colorsToBeta;
+ colorsToBeta.resize(0);
+ colorsToBeta.resize(colorsNames.size());
+ std::vector< std::vector< double > > colorsToBetaAngles;
+ colorsToBetaAngles.resize(0);
+ std::vector< gmx::RVec > betaListsRVecs;
+ betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs);
+ for (size_t i = 0; i < colorsToPeptide.size(); ++i) {
+ if (colorsToPeptide[i]) {
+ searchNearBetaLists(fr, pbc, trajectoryFrame, nb_, betaLists[frnr], colorsNames[i], colorsToBeta[i]);
+ }
+ }
+ for (size_t i = 0; i < colorsToBeta.size(); ++i) {
+ for (size_t j = 0; j < colorsToBeta[i].size(); ++j) {
+ if (colorsToBeta[i][j]) {
+ computeAnglesColorsVsBeta(betaListsRVecs[j], colorFormation[i]);
+ }
+ }
+ }
+ // вывод в файл(ы) информацию для фрейма
+ anglesFileDump(frnr, fnOut, colorsToPeptide, colorFormation);
+}
+
+void
+colorVec::finishAnalysis(int nframes)
+{
+
+}
+
+void
+colorVec::writeOutput()
+{
+ /*
+ *
+ * frame #<current frame number>
+ * color #<color number> <yes/no for "close" to peptide> a12 b12 n12 <number of beta lists> <6 avg beta angles> [<6 angles for each beta list>]
+ *
+ */
+ std::cout << "\n\t colorAngles finished successfully\n";
+}
+
+/*! \brief
+ * The main function for the analysis template.
+ */
+int
+main(int argc, char *argv[])
+{
+ return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<colorVec>(argc, argv);
+}
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2013,2014,2015, 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.
- */
-/*! \internal \file
- * \brief
- * Implements gmx::analysismodules::Freevolume.
- *
- * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
- * \ingroup module_trajectoryanalysis
- */
-
-#include "domaintype.h"
-#include "newfit.h"
-
-#include <gromacs/trajectoryanalysis.h>
-#include <gromacs/trajectoryanalysis/topologyinformation.h>
-//#include "/home/titov_ai/Develop/gromacs-original/src/gromacs/trajectoryanalysis/topologyinformation.h"
-
-#include <iostream>
-//#include <chrono>
-#include <cfloat>
-#include <omp.h>
-#include <vector>
-#include <cstdio>
-
-using namespace gmx;
-
-using gmx::RVec;
-
-/*! \brief
- * \ingroup module_trajectoryanalysis
- */
-class Domains : public TrajectoryAnalysisModule
-{
- public:
-
- Domains();
- virtual ~Domains();
-
- //! Set the options and setting
- virtual void initOptions(IOptionsContainer *options,
- TrajectoryAnalysisSettings *settings);
-
- //! First routine called by the analysis framework
- // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
- virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
- const TopologyInformation &top);
-
- //! Call for each frame of the trajectory
- // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
- virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
- TrajectoryAnalysisModuleData *pdata);
-
- //! Last routine called by the analysis framework
- // virtual void finishAnalysis(t_pbc *pbc);
- virtual void finishAnalysis(int nframes);
-
- //! Routine to write output, that is additional over the built-in
- virtual void writeOutput();
-
- private:
-
- std::string fnNdx_;
- domainType testSubject;
- std::vector< size_t > index;
- std::vector< RVec > trajectory;
- std::vector< RVec > reference;
- Selection selec;
- std::vector< std::vector< std::pair< size_t, size_t > > > fitPairs;
- unsigned int bone;
- int window {-1}; // selectable
- int domain_min_size {-1}; // selectable
- double delta {-1}; // selectable
- double epsi {-1}; // selectable
- int twStep {-1}; // selectable
- int DomainSearchingAlgorythm {-1}; // selectable
- // Copy and assign disallowed by base.
-};
-
-Domains::Domains(): TrajectoryAnalysisModule()
-{
-}
-
-Domains::~Domains()
-{
-}
-
-void
-Domains::initOptions(IOptionsContainer *options,
- TrajectoryAnalysisSettings *settings)
-{
- static const char *const desc[] = {
- "[THISMODULE] to be done"
- "-f trajectory.xtc -s structure.tpr -n index.ndx -sf selectionListFile -on outPutName -dms 4+ -DSA 0/1/2 -eps 0.2 -dlt 0.95 -wSize 1000 -twStep 100"
- };
- // Add the descriptive text (program help text) to the options
- 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."));
- // Add option for output file name
- options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
- .store(&fnNdx_).defaultBasename("domains")
- .description("Index file for the algorithm output."));
- // Add option for domain min size constant
- options->addOption(gmx::IntegerOption("dms")
- .store(&domain_min_size)
- .description("minimum domain size, should be >= 4."));
- // Add option for Domain's Searching Algorythm
- options->addOption(gmx::IntegerOption("DSA")
- .store(&DomainSearchingAlgorythm)
- .description("Domain's Searching Algorythm: 0 == (default) from bigger to smaller | 1 == just the first biggest domain | 2 == from smaller to bigger."));
- // Add option for epsi constant
- options->addOption(DoubleOption("eps")
- .store(&epsi)
- .description("thermal vibrations' constant."));
- // Add option for delta constant
- options->addOption(DoubleOption("dlt")
- .store(&delta)
- .description("domain membership probability."));
- // Add option for domain min size constant
- options->addOption(gmx::IntegerOption("wSize")
- .store(&window)
- .description("flowing window to evaluate domains from."));
- // Add option for time step between windows' starts
- options->addOption(gmx::IntegerOption("twStep")
- .store(&twStep)
- .description("time step between windows' starting positions."));
- // Control input settings
- settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
- settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
- //settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
- settings->setPBC(true);
-}
-
-// -f trajectory.xtc -s structure.tpr -n index.ndx -sf selectionListFile -on outPutName -dms 4+ -DSA 0/1/2 -eps 0.2 -dlt 0.95 -wSize 1000 -twStep 100
-
-void
-Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
- const TopologyInformation &top)
-{
- // считывание индекса
- index.resize(0);
- for (ArrayRef< const int >::iterator ai {selec.atomIndices().begin()}; (ai < selec.atomIndices().end()); ai++) {
- index.push_back(static_cast< size_t >(*ai));
- }
-
- // считывание референсной структуры
- reference.resize(0);
-
- if (top.hasFullTopology()) {
- for (const auto &i : index) {
- reference.push_back(top.x().at(i));
- }
- }
-
- // вычисление числа основных (от основание) последовательностей
- // в идеале нужно брать все сочетания, но ограничиваемся последовательными наборами
- bone = static_cast< size_t >(index.size() - static_cast< size_t >(domain_min_size) + 1);
-
- // создание пар для фитирования структур
- fitPairs.resize(0);
- fitPairs.resize(bone);
- for (size_t i {0}; i < bone; i++) {
- fitPairs[i].resize(0);
- for (size_t j {0}; j < domain_min_size; j++) {
- fitPairs[i].push_back(std::make_pair(j + i, j + i));
- }
- }
-
- trajectory.resize(0);
- trajectory.resize(index.size());
-
- // задание необходимых параметров для вычисления доменов
- testSubject.setDefaults(index, reference, window, domain_min_size, DomainSearchingAlgorythm, twStep, bone, epsi, delta, fnNdx_);
-}
-
-void
-Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
-{
- // считывания текущего фрейма траектории
- for (size_t i {0}; i < index.size(); i++) {
- trajectory[i] = fr.x[index[i]];
- }
-
- // создание копий фрейма для каждой основной последовательности
- std::vector< std::vector< RVec > > tsTemp;
- tsTemp.resize(0);
- tsTemp.resize(bone, trajectory);
-
- // фитирование каждой копии
- #pragma omp parallel for ordered schedule(dynamic) firstprivate(reference)
- for (size_t i = 0; i < bone; i++) {
- MyFitNew(reference, tsTemp[i], fitPairs[i], 0.000'001);
- }
- #pragma omp barrier
-
- // (обновление/потенциальное выделение) доменов
- testSubject.update(tsTemp, frnr);
-}
-
-void
-Domains::finishAnalysis(int nframes)
-{
-
-}
-
-void
-Domains::writeOutput()
-{
- std::cout << "\n END \n";
-}
-
-/*! \brief
- * The main function for the analysis template.
- */
-int
-main(int argc, char *argv[])
-{
- return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
-}
biod.pnpi.spb.ru / alexxy / gromacs-colorvec.git / commitdiff