- fixed setprec...(3 to 5)
authorAnatoly <Titov_AI@pnpi.nrcki.ru>
Fri, 23 Oct 2020 09:31:10 +0000 (12:31 +0300)
committerAnatoly <Titov_AI@pnpi.nrcki.ru>
Fri, 23 Oct 2020 09:31:10 +0000 (12:31 +0300)
- added pbc_dx()
- added some "safe-pillow" output

src/colorvec.cpp

index 9128e6726a0ecbf8b4eba3bcf80840d6c5b32e71..7eedff20607ca871b1ecac1f187ffd14e82f88fe 100644 (file)
@@ -43,6 +43,7 @@
 #include <gromacs/trajectoryanalysis.h>
 #include <gromacs/trajectoryanalysis/topologyinformation.h>
 #include <gromacs/selection/nbsearch.h>
+#include <gromacs/pbcutil/pbc.h>
 
 #include <iostream>
 #include <fstream>
@@ -211,11 +212,14 @@ inline void betaListsRVecsEvaluation(const std::vector< gmx::RVec > &frame, cons
 //}
 
 // определение близко ли к белку находится краска
-bool isNearPeptide( const std::vector< gmx::RVec > &inputFrame, const std::vector< size_t > &inputIndex,
+bool isNearPeptide( const t_pbc *inputPBC,
+                    const std::vector< gmx::RVec > &inputFrame, const std::vector< size_t > &inputIndex,
                     const std::vector< std::pair< std::string, size_t > > &inputColor, const double cutOff) {
+    gmx::RVec tempRVec;
     for (size_t i {0}; i < inputColor.size(); ++i) {
         for (size_t j {0}; j < inputIndex.size(); ++j) {
-            if ((inputFrame[inputIndex[j]] - inputFrame[inputColor[i].second]).norm() <= cutOff) {
+            pbc_dx(inputPBC, inputFrame[inputIndex[j]], inputFrame[inputColor[i].second], tempRVec);
+            if (tempRVec.norm() <= cutOff) {
                 return true;
             }
         }
@@ -248,16 +252,19 @@ bool isNearPeptide( const std::vector< gmx::RVec > &inputFrame, const std::vecto
 //    }
 //}
 
-inline void searchNearBetaLists(const std::vector< gmx::RVec > &inputFrame, const std::vector< std::vector< unsigned int > > &inputBLists,
+inline void searchNearBetaLists(const t_pbc *inputPBC,
+                                const std::vector< gmx::RVec > &inputFrame, const std::vector< std::vector< unsigned int > > &inputBLists,
                                 const std::vector< std::pair< std::string, size_t > > &inputColor, std::vector< bool > &outputList,
                                 const std::vector< std::vector< size_t > > &inputAminoacids, const double cutOff) {
     outputList.resize(0);
     outputList.resize(inputBLists.size(), false);
+    gmx::RVec tempRVec;
     for (size_t i {0}; i < inputColor.size(); ++i) {
         for (size_t j = 0; j < inputBLists.size(); ++j) {
             for (size_t k = 0; k < inputBLists[j].size(); ++k) {
                 for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
-                    if ((inputFrame[inputAminoacids[inputBLists[j][k]][m]] - inputFrame[inputColor[i].second]).norm() <= cutOff) {
+                    pbc_dx(inputPBC, inputFrame[inputAminoacids[inputBLists[j][k]][m]], inputFrame[inputColor[i].second], tempRVec);
+                    if (tempRVec.norm() <= cutOff) {
                         outputList[j] = true;
                         break;
                     }
@@ -288,9 +295,10 @@ void anglesFileDump(const int frameNum, const std::string &output, const std::ve
         } else {
             file << std::setw(4) << "no";
         }
-        file << std::setw(8) << std::setprecision(3) << colorFormation[i].a12;
-        file << std::setw(8) << std::setprecision(3) << colorFormation[i].b12;
-        file << std::setw(8) << std::setprecision(3) << colorFormation[i].n12;
+        file << std::fixed;
+        file << std::setw(8) << std::setprecision(5) << colorFormation[i].a12;
+        file << std::setw(8) << std::setprecision(5) << colorFormation[i].b12;
+        file << std::setw(8) << std::setprecision(5) << colorFormation[i].n12;
         if (colorFormation[i].betaAngles.size() == 0) {
             file << std::setw(4) << 0;
         } else {
@@ -306,10 +314,10 @@ void anglesFileDump(const int frameNum, const std::string &output, const std::ve
             }
             file << std::setw(4) << temp;
             for (size_t j = 0; j < betaTemp.size(); ++j) {
-                file << std::setw(8) << std::setprecision(3) << betaTemp[j];
+                file << std::setw(8) << std::setprecision(5) << betaTemp[j];
             }
             for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
-                file << std::setw(8) << std::setprecision(3) << colorFormation[i].betaAngles[j];
+                file << std::setw(8) << std::setprecision(5) << colorFormation[i].betaAngles[j];
             }
         }
         file << std::endl;
@@ -440,6 +448,17 @@ colorVec::initAnalysis( const gmx::TrajectoryAnalysisSettings   &settings,
     aminoacidsIndexation(index, top, aminoacidsIndex);
     // задание радиуса рассматриваемых соседей
     nb_.setCutoff(effRad);
+    // вывод базовых значений
+    std::cout << std::endl;
+    std::cout << "effective radius = " << std::setw(20) << effRad                   << std::endl;
+    std::cout << "DAT file         = " << std::setw(20) << fnBetaListsDat           << std::endl;
+    std::cout << "output file      = " << std::setw(20) << fnOut                    << std::endl;
+    std::cout << "peptide atoms #  = " << std::setw(20) << index.size()             << std::endl;
+    std::cout << "CA-atmos #       = " << std::setw(10) << indexCA.size()           << std::endl;
+    std::cout << "Beta frames      = " << std::setw(10) << betaLists.size()         << std::endl;
+    std::cout << "residue #        = " << std::setw(10) << aminoacidsIndex.size()   << std::endl;
+    std::cout << "colors #         = " << std::setw(10) << colorsNames.size()       << std::endl;
+
     /*
      * формально можно сделать:
      * найти соотношение между именами для углов и индексными номерами
@@ -470,11 +489,13 @@ colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::Trajecto
      *
      */
 
+
+
     std::vector< bool >     colorsToPeptide;
     colorsToPeptide.resize(0);
     colorsToPeptide.resize(colorsNames.size(), false);
     for (size_t i = 0; i < colorsNames.size(); ++i) {
-        colorsToPeptide[i] = isNearPeptide(trajectoryFrame, index, colorsNames[i], effRad);
+        colorsToPeptide[i] = isNearPeptide(pbc, trajectoryFrame, index, colorsNames[i], effRad);
     }
     // расчёт угла и среднего угла с ближайшими бета листами
     std::vector< std::vector< bool > > colorsToBeta;
@@ -486,7 +507,7 @@ colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::Trajecto
     betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA);
     for (size_t i = 0; i < colorsToPeptide.size(); ++i) {
         if (colorsToPeptide[i]) {
-            searchNearBetaLists(trajectoryFrame, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex, effRad);
+            searchNearBetaLists(pbc, trajectoryFrame, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex, effRad);
         }
     }
     for (size_t i = 0; i < colorsToBeta.size(); ++i) {