- fixed an issue with peptide being only of CA

diff --git a/src/colorvec.cpp b/src/colorvec.cpp
index 31fc750d3190fe4625be37d0e3f3f291391741c6..b562ba91a900178192202cab8b6a345aa2d4c944 100644 (file)
--- a/src/colorvec.cpp
+++ b/src/colorvec.cpp
@@ -120,6 +120,15 @@ void betaListDigestion(const std::string &inputFile, std::vector< std::vector< s
     file.close();
 }
 
+// функция выделения индексов для бэталистов
+inline void aminoacidsIndexation(const std::vector< size_t > inputIndex, const TopologyInformation &top, std::vector< std::vector< size_t > > &aminoacidsIndex) {
+    aminoacidsIndex.resize(0);
+    for (size_t i = 0; i < inputIndex.size(); ++i) {
+        aminoacidsIndex.resize(std::max(aminoacidsIndex.size(), static_cast< size_t >(top.atoms()->atom[inputIndex[i]].resind)));
+        aminoacidsIndex[top.atoms()->atom[inputIndex[i]].resind].push_back(inputIndex[i]);
+    }
+}
+
 // функция поиска 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) {
@@ -171,11 +180,11 @@ void colorsAnglesEvaluation(const std::vector< RVec > &frame, const std::vector<
 }
 
 // вычисление направляющего вектора в бэта-листе
-inline void betaListsRVecsEvaluation(const std::vector< RVec > &frame, const std::vector< std::vector< unsigned int > > &inputBetaLists, std::vector< gmx::RVec > &temp) {
+inline void betaListsRVecsEvaluation(const std::vector< RVec > &frame, const std::vector< std::vector< unsigned int > > &inputBetaLists, std::vector< gmx::RVec > &temp, const std::vector< size_t > &inputCA) {
     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 = frame[inputCA[i[i.size() - 1]]] + frame[inputCA[i[i.size() - 2]]] - frame[inputCA[i[1]]] - frame[inputCA[i[0]]];
         tempA /= tempA.norm();
         temp.push_back(tempA);
     }
@@ -204,7 +213,7 @@ inline void searchNearBetaLists(const t_trxframe &fr, const t_pbc *pbc, const st
                                 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) {
+                                std::vector< bool > &outputList, const std::vector< std::vector< size_t > > &inputAminoacids) {
     outputList.resize(0);
     outputList.resize(inputBLists.size(), false);
     gmx::AnalysisNeighborhoodSearch      nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
@@ -213,8 +222,10 @@ inline void searchNearBetaLists(const t_trxframe &fr, const t_pbc *pbc, const st
         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;
+                    for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
+                        if (pair.testIndex() == inputAminoacids[inputBLists[j][k]][m]) {
+                            outputList[j] = true;
+                        }
                     }
                 }
             }
@@ -309,6 +320,8 @@ class colorVec : public TrajectoryAnalysisModule
         std::string                                                     fnBetaListsDat;         // selectable
         float                                                           effRad          {0.8};  // selectable
         std::vector< size_t >                                           index;
+        std::vector< size_t >                                           indexCA;
+        std::vector< std::vector< size_t > >                            aminoacidsIndex;
         std::vector< std::vector< std::pair< std::string, size_t > > >  colorsNames;
         std::vector< std::vector< std::vector< unsigned int > > >       betaLists;
         AnalysisNeighborhood                                            nb_;
@@ -327,6 +340,7 @@ colorVec::~colorVec()
 /*
  * ndx:
  * [peptide]
+ * [CA]
  * [color_{i}]
  *
  */
@@ -372,15 +386,22 @@ colorVec::initAnalysis(const TrajectoryAnalysisSettings &settings,
     for (ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ai++) {
         index.push_back(static_cast< size_t >(*ai));
     }
+    // считывание индекса
+    indexCA.resize(0);
+    for (ArrayRef< const int >::iterator ai {sel_[1].atomIndices().begin()}; (ai < sel_[2].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) {
+    colorsNames.resize(sel_.size() - 2);
+    for (size_t i = 2; 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)));
+            colorsNames[i - 2].push_back(std::make_pair(*(top.atoms()->resinfo[top.atoms()->atom[*ai].resind].name), static_cast< size_t >(*ai)));
         }
     }
     // разбор dat файла для создания бэта листов
     betaListDigestion(fnBetaListsDat, betaLists);
+    // разбор топологии для индексации бета листов
+    aminoacidsIndexation(index, top, aminoacidsIndex);
     // задание радиуса рассматриваемых соседей
     nb_.setCutoff(effRad);
     /*
@@ -426,10 +447,10 @@ colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAna
     std::vector< std::vector< double > > colorsToBetaAngles;
     colorsToBetaAngles.resize(0);
     std::vector< gmx::RVec > betaListsRVecs;
-    betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs);
+    betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA);
     for (size_t i = 0; i < colorsToPeptide.size(); ++i) {
         if (colorsToPeptide[i]) {
-            searchNearBetaLists(fr, pbc, trajectoryFrame, nb_, betaLists[frnr], colorsNames[i], colorsToBeta[i]);
+            searchNearBetaLists(fr, pbc, trajectoryFrame, nb_, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex);
         }
     }
     for (size_t i = 0; i < colorsToBeta.size(); ++i) {