src/colorvec.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements gmx::analysismodules::Freevolume.
  38  *
  39  * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42
  43 #include <gromacs/trajectoryanalysis.h>
  44 #include <gromacs/trajectoryanalysis/topologyinformation.h>
  45 #include <gromacs/selection/nbsearch.h>
  46 #include <gromacs/pbcutil/pbc.h>
  47
  48 #include <iostream>
  49 #include <fstream>
  50 //#include <chrono>
  51 #include <cfloat>
  52 #include <omp.h>
  53 #include <vector>
  54 #include <cstdio>
  55 #include <iomanip>
  56
  57 // структура углов для одной краски
  58 struct colorLocalAngles {
  59     gmx::RVec   a1{0., 0., 0.}, a2{0., 0., 0.};
  60     gmx::RVec   b1{0., 0., 0.}, b2{0., 0., 0.};
  61     gmx::RVec   n1{0., 0., 0.}, n2{0., 0., 0.};
  62     double      a12 {0}, b12 {0}, n12 {0};
  63     std::vector< double > betaAngles;
  64 };
  65
  66 // функция парсинга одной строки
  67 void parseBetaListDATLine(const std::string &currentLine, std::vector< std::vector< unsigned int > > &localInputBL) {
  68     size_t                      equalCount = 0;
  69     std::vector< unsigned int > a;
  70     a.resize(0);
  71     for (size_t i = 0; i < currentLine.size(); ++i) {
  72         if (currentLine[i] == '=') { // подсчитываем число "пустых" символов
  73             ++equalCount;
  74         } else {
  75             size_t temp {i};
  76             for (size_t j = i + 1; j < currentLine.size(); ++j) {
  77                 if (currentLine[j] == 'B' || currentLine[j] == 'E') {
  78                     ++temp;
  79                 } else {
  80                     break;
  81                 }
  82             }
  83             if (temp - i > 3) {
  84                 localInputBL.push_back(a);
  85                 for (size_t j = i; j <= temp; ++j) {
  86                     localInputBL.back().push_back(j - equalCount);
  87                 }
  88                 i = temp;
  89             }
  90         }
  91     }
  92 }
  93
  94 // функция нахождения бета-листов в структуре по файлу ДССП
  95 void betaListDigestion(const std::string &inputFile, std::vector< std::vector< std::vector< unsigned int > > > &inputBL) {
  96     inputBL.resize(0);
  97     std::ifstream   file(inputFile);
  98     std::string     line;
  99     getline(file, line); // считываем число в первой строке - кол-во осмысленных элементов в строках - нам не нужно
 100     getline(file, line);
 101     if (line.size() > 3) {
 102         do {
 103             inputBL.resize(inputBL.size() + 1);
 104             parseBetaListDATLine(line, inputBL.back());
 105             line = "";
 106             getline(file, line);
 107         } while (line.size() > 3);
 108     } else {
 109         throw "DSSP DAT FILE IS EMPTY";
 110     }
 111     file.close();
 112 }
 113
 114 // функция выделения индексов для бэталистов
 115 inline void aminoacidsIndexation(const std::vector< size_t > &inputIndex, const gmx::TopologyInformation &top, std::vector< std::vector< size_t > > &aminoacidsIndex) {
 116     aminoacidsIndex.resize(0);
 117     for (size_t i = 0; i < inputIndex.size(); ++i) {
 118         aminoacidsIndex.resize(std::max(aminoacidsIndex.size(), static_cast< size_t >(top.atoms()->atom[inputIndex[i]].resind + 1)));
 119         aminoacidsIndex[top.atoms()->atom[inputIndex[i]].resind].push_back(inputIndex[i]);
 120     }
 121 }
 122
 123 // функция поиска RVec в кадре по имени->индексу
 124 gmx::RVec returnRVec(const std::vector< gmx::RVec > &frame, const std::vector< std::pair< std::string, size_t > > &colorsIndex, const std::string &toFind) {
 125     for (auto &i : colorsIndex) {
 126         if (i.first == toFind) {
 127             return frame[i.second];
 128         }
 129     }
 130     throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS";
 131 }
 132
 133 // функция векторного произведения двух RVec'ов
 134 inline void RVecVecMultiply(gmx::RVec &n, const gmx::RVec &a, const gmx::RVec &b) {
 135     n[0] = a[1] * b[2] - a[2] * b[1];
 136     n[1] = a[2] * b[0] - a[0] * b[2];
 137     n[2] = a[0] * b[1] - a[2] * b[0];
 138 }
 139
 140 // поиск угла между двумя RVec'ами
 141 inline double RVecAngle(const gmx::RVec &a, const gmx::RVec &b) {
 142     return std::acos((a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) / (a.norm() * b.norm())) * 180.0 / 3.14159265;
 143 }
 144
 145 // вычисление внутренних углов в краске
 146 void colorsAnglesEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex,
 147                             std::vector< colorLocalAngles > &colorStruct) {
 148     colorStruct.resize(colorsIndex.size());
 149     #pragma omp parallel for ordered schedule(dynamic)
 150     for (size_t i = 0; i < colorsIndex.size(); ++i) {
 151         colorStruct[i].betaAngles.resize(0);
 152         // "левый" блок
 153         colorStruct[i].a1 = returnRVec(frame, colorsIndex[i], "CAF") - returnRVec(frame, colorsIndex[i], "CAJ") + returnRVec(frame, colorsIndex[i], "CAK") - returnRVec(frame, colorsIndex[i], "CAO");
 154         colorStruct[i].a1 /= colorStruct[i].a1.norm();
 155         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");
 156         colorStruct[i].b1 /= colorStruct[i].b1.norm();
 157         RVecVecMultiply(colorStruct[i].n1, colorStruct[i].a1, colorStruct[i].b1);
 158         colorStruct[i].n1 /= colorStruct[i].n1.norm();
 159         // "правый" блок
 160         colorStruct[i].a2 = returnRVec(frame, colorsIndex[i], "CAB") - returnRVec(frame, colorsIndex[i], "CBQ") + returnRVec(frame, colorsIndex[i], "CBP") - returnRVec(frame, colorsIndex[i], "CBL");
 161         colorStruct[i].a2 /= colorStruct[i].a2.norm();
 162         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");
 163         colorStruct[i].b2 /= colorStruct[i].b2.norm();
 164         RVecVecMultiply(colorStruct[i].n2, colorStruct[i].a2, colorStruct[i].b2);
 165         colorStruct[i].n2 *= -1; // для "сонаправленности" векторов нормали
 166         colorStruct[i].n2 /= colorStruct[i].n2.norm();
 167         colorStruct[i].a12 = RVecAngle(colorStruct[i].a1, colorStruct[i].a2);
 168         colorStruct[i].b12 = RVecAngle(colorStruct[i].b1, colorStruct[i].b2);
 169         colorStruct[i].n12 = RVecAngle(colorStruct[i].n1, colorStruct[i].n2);
 170     }
 171     #pragma omp barrier
 172 }
 173
 174 // вычисление направляющего вектора в бэта-листе
 175 inline void betaListsRVecsEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< std::vector< unsigned int > > &inputBetaLists,
 176                                      std::vector< gmx::RVec > &temp, const std::vector< size_t > &inputCA) {
 177     temp.resize(0);
 178     gmx::RVec tempA;
 179     for (const auto &i : inputBetaLists) {
 180         tempA = frame[inputCA[i[i.size() - 1]]] + frame[inputCA[i[i.size() - 2]]] - frame[inputCA[i[1]]] - frame[inputCA[i[0]]];
 181         tempA /= tempA.norm();
 182         temp.push_back(tempA);
 183     }
 184 }
 185
 186 //// определение близко ли к белку находится краска
 187 //bool isNearPeptide(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > &inputFrame,
 188 //                          /*gmx::AnalysisNeighborhood &nbhood, */const std::vector< size_t > &inputIndex,
 189 //                          const std::vector< std::pair< std::string, size_t > > &inputColor) {
 190 //    /*gmx::AnalysisNeighborhoodSearch     nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
 191 //    gmx::AnalysisNeighborhoodPair       pair;*/
 192
 193 //    for (size_t i = 0; i < inputColor.size(); ++i) {
 194 //        std::cout << inputColor[i].first << " " << i << " / "  << inputColor.size() << std::endl;
 195 //        gmx::AnalysisNeighborhood           nbhood;
 196 //        nbhood.setCutoff(0.8);
 197 //        gmx::AnalysisNeighborhoodSearch     nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
 198 //        gmx::AnalysisNeighborhoodPair       pair;
 199 //        gmx::AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(inputFrame[inputColor[i].second].as_vec());
 200 //        int count1 {0};
 201 //        while (pairSearch.findNextPair(&pair)) {
 202 //            std::cout << ++count1 << " ";
 203 //            for (size_t j = 0; j < inputIndex.size(); ++j) {
 204 //                if (pair.refIndex() == inputIndex[j]) {
 205 //                    return true;
 206 //                }
 207 //            }
 208 //        }
 209 //        std::cout << inputColor[i].first << " atom done" << std::endl;
 210 //    }
 211 //    return false;
 212 //}
 213
 214 // определение близко ли к белку находится краска
 215 bool isNearPeptide( const t_pbc *inputPBC,
 216                     const std::vector< gmx::RVec > &inputFrame, const std::vector< size_t > &inputIndex,
 217                     const std::vector< std::pair< std::string, size_t > > &inputColor, const double cutOff) {
 218     gmx::RVec tempRVec;
 219     for (size_t i {0}; i < inputColor.size(); ++i) {
 220         for (size_t j {0}; j < inputIndex.size(); ++j) {
 221             pbc_dx(inputPBC, inputFrame[inputIndex[j]], inputFrame[inputColor[i].second], tempRVec);
 222             if (tempRVec.norm() <= cutOff) {
 223                 return true;
 224             }
 225         }
 226     }
 227     return false;
 228 }
 229
 230 // поиск ближайших к краске бэта-листов
 231 //inline void searchNearBetaLists(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > &inputFrame,
 232 //                                gmx::AnalysisNeighborhood &nbhood,
 233 //                                const std::vector< std::vector< unsigned int > > &inputBLists,
 234 //                                const std::vector< std::pair< std::string, size_t > > &inputColor,
 235 //                                std::vector< bool > &outputList, const std::vector< std::vector< size_t > > &inputAminoacids) {
 236 //    outputList.resize(0);
 237 //    outputList.resize(inputBLists.size(), false);
 238 //    gmx::AnalysisNeighborhoodSearch      nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
 239 //    gmx::AnalysisNeighborhoodPair        pair;
 240 //    for (const auto &i : inputColor) {
 241 //        while (nbsearch.startPairSearch(inputFrame[i.second].as_vec()).findNextPair(&pair)) {
 242 //            for (size_t j = 0; j < inputBLists.size(); ++j) {
 243 //                for (size_t k = 0; k < inputBLists[j].size(); ++k) {
 244 //                    for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
 245 //                        if (pair.testIndex() == inputAminoacids[inputBLists[j][k]][m]) {
 246 //                            outputList[j] = true;
 247 //                        }
 248 //                    }
 249 //                }
 250 //            }
 251 //        }
 252 //    }
 253 //}
 254
 255 inline void searchNearBetaLists(const t_pbc *inputPBC,
 256                                 const std::vector< gmx::RVec > &inputFrame, const std::vector< std::vector< unsigned int > > &inputBLists,
 257                                 const std::vector< std::pair< std::string, size_t > > &inputColor, std::vector< bool > &outputList,
 258                                 const std::vector< std::vector< size_t > > &inputAminoacids, const double cutOff) {
 259     outputList.resize(0);
 260     outputList.resize(inputBLists.size(), false);
 261     gmx::RVec tempRVec;
 262     for (size_t i {0}; i < inputColor.size(); ++i) {
 263         for (size_t j = 0; j < inputBLists.size(); ++j) {
 264             for (size_t k = 0; k < inputBLists[j].size(); ++k) {
 265                 for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
 266                     pbc_dx(inputPBC, inputFrame[inputAminoacids[inputBLists[j][k]][m]], inputFrame[inputColor[i].second], tempRVec);
 267                     if (tempRVec.norm() <= cutOff) {
 268                         outputList[j] = true;
 269                         break;
 270                     }
 271                 }
 272             }
 273         }
 274     }
 275 }
 276
 277 // определение углов между краской и направляющим вектором бэта-листа
 278 inline void computeAnglesColorsVsBeta(const gmx::RVec &inputBetaRVec, colorLocalAngles &colorFormation) {
 279     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a1));
 280     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b1));
 281     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n1));
 282     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a2));
 283     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b2));
 284     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n2));
 285 }
 286
 287 // функция записи в файл значений углов для кадра
 288 void anglesFileDump(const int frameNum, const std::string &output, const std::vector< bool > &toPeptide, const std::vector< colorLocalAngles > &colorFormation) {
 289     std::ofstream   file(output, std::ofstream::app);
 290     file << "frame =" << std::setw(8) << frameNum << std::endl;
 291     for (size_t i = 0; i < colorFormation.size(); ++i) {
 292         file << "color #" << std::setw(3) << i;
 293         if (toPeptide[i]) {
 294             file << std::setw(4) << "yes";
 295         } else {
 296             file << std::setw(4) << "no";
 297         }
 298         file << std::fixed;
 299         file << std::setw(8) << std::setprecision(2) << colorFormation[i].a12;
 300         file << std::setw(8) << std::setprecision(2) << colorFormation[i].b12;
 301         file << std::setw(8) << std::setprecision(2) << colorFormation[i].n12;
 302         if (colorFormation[i].betaAngles.size() == 0) {
 303             file << std::setw(4) << 0;
 304         } else {
 305             int temp = colorFormation[i].betaAngles.size() / 6;
 306             std::vector< double > betaTemp;
 307             betaTemp.resize(0);
 308             betaTemp.resize(6, 0.); // magic number, meh
 309             for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
 310                 betaTemp[j % 6] += colorFormation[i].betaAngles[j];
 311             }
 312             for (size_t j = 0; j < betaTemp.size(); ++j) {
 313                 betaTemp[j] /= static_cast< double >(temp);
 314             }
 315             file << std::setw(4) << temp;
 316             for (size_t j = 0; j < betaTemp.size(); ++j) {
 317                 file << std::setw(8) << std::setprecision(2) << betaTemp[j];
 318             }
 319             for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
 320                 file << std::setw(8) << std::setprecision(2) << colorFormation[i].betaAngles[j];
 321             }
 322         }
 323         file << std::endl;
 324     }
 325     file.close();
 326 }
 327
 328 /*! \brief
 329  * \ingroup module_trajectoryanalysis
 330  */
 331 class colorVec : public gmx::TrajectoryAnalysisModule
 332 {
 333     public:
 334
 335         colorVec();
 336         virtual ~colorVec();
 337
 338         //! Set the options and setting
 339         virtual void initOptions(gmx::IOptionsContainer          *options,
 340                                  gmx::TrajectoryAnalysisSettings *settings);
 341
 342         //! First routine called by the analysis framework
 343         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 344         virtual void initAnalysis(const gmx::TrajectoryAnalysisSettings &settings,
 345                                   const gmx::TopologyInformation        &top);
 346
 347         //! Call for each frame of the trajectory
 348         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 349         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 350                                   gmx::TrajectoryAnalysisModuleData *pdata);
 351
 352         //! Last routine called by the analysis framework
 353         // virtual void finishAnalysis(t_pbc *pbc);
 354         virtual void finishAnalysis(int nframes);
 355
 356         //! Routine to write output, that is additional over the built-in
 357         virtual void writeOutput();
 358
 359     private:
 360
 361         gmx::SelectionList                                              sel_;
 362         std::string                                                     fnOut           {"name"};   // selectable
 363         std::string                                                     fnBetaListsDat;             // selectable
 364         float                                                           effRad          {0.8};      // selectable
 365         std::vector< size_t >                                           index;
 366         std::vector< size_t >                                           indexCA;
 367         std::vector< std::vector< size_t > >                            aminoacidsIndex;
 368         std::vector< std::vector< std::pair< std::string, size_t > > >  colorsNames;
 369         std::vector< std::vector< std::vector< unsigned int > > >       betaLists;
 370         gmx::AnalysisNeighborhood                                       nb_;
 371
 372         // Copy and assign disallowed by base.
 373 };
 374
 375 colorVec::colorVec(): TrajectoryAnalysisModule()
 376 {
 377 }
 378
 379 colorVec::~colorVec()
 380 {
 381 }
 382
 383 /*
 384  * ndx:
 385  * [peptide]
 386  * [CA]
 387  * [color_{i}]
 388  *
 389  */
 390
 391 void
 392 colorVec::initOptions(  gmx::IOptionsContainer          *options,
 393                         gmx::TrajectoryAnalysisSettings *settings)
 394 {
 395     static const char *const desc[] = {
 396         "[THISMODULE] to be done"
 397     };
 398     // Add the descriptive text (program help text) to the options
 399     settings->setHelpText(desc);
 400     // Add option for selection list
 401     options->addOption(gmx::SelectionOption("sel")
 402                             .storeVector(&sel_)
 403                             .required().dynamicMask().multiValue()
 404                             .description("select pepride and colors / -sf"));
 405     // Add option for input file names
 406     options->addOption(gmx::StringOption("dat")
 407                             .store(&fnBetaListsDat)
 408                             .description("a file to make dynamic beta lists"));
 409     // Add option for output file name
 410     options->addOption(gmx::StringOption("out")
 411                             .store(&fnOut)
 412                             .description("Index file for the algorithm output."));
 413     // Add option for effRad constant
 414     options->addOption(gmx::FloatOption("efRad")
 415                             .store(&effRad)
 416                             .description("max distance from colors to peptide in nm to consider to be \"near\""));
 417     // Control input settings
 418     settings->setFlags(gmx::TrajectoryAnalysisSettings::efNoUserPBC);
 419     settings->setFlag(gmx::TrajectoryAnalysisSettings::efUseTopX);
 420     //settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 421     settings->setPBC(true);
 422 }
 423
 424 void
 425 colorVec::initAnalysis( const gmx::TrajectoryAnalysisSettings   &settings,
 426                         const gmx::TopologyInformation          &top)
 427 {
 428     // считывание индекса
 429     index.resize(0);
 430     for (gmx::ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ai++) {
 431         index.push_back(static_cast< size_t >(*ai));
 432     }
 433     // считывание индекса
 434     indexCA.resize(0);
 435     for (gmx::ArrayRef< const int >::iterator ai {sel_[1].atomIndices().begin()}; (ai < sel_[1].atomIndices().end()); ai++) {
 436         indexCA.push_back(static_cast< size_t >(*ai));
 437     }
 438     // считывание красок
 439     colorsNames.resize(sel_.size() - 2);
 440     for (size_t i = 2; i < sel_.size(); ++i) {
 441         for (gmx::ArrayRef< const int >::iterator ai {sel_[i].atomIndices().begin()}; (ai < sel_[i].atomIndices().end()); ai++) {
 442             colorsNames[i - 2].push_back(std::make_pair(*(top.atoms()->atomname[*ai]), static_cast< size_t >(*ai)));
 443         }
 444     }
 445     // разбор dat файла для создания бэта листов
 446     betaListDigestion(fnBetaListsDat, betaLists);
 447     // разбор топологии для индексации бета листов
 448     aminoacidsIndexation(index, top, aminoacidsIndex);
 449     // задание радиуса рассматриваемых соседей
 450     nb_.setCutoff(effRad);
 451     // вывод базовых значений
 452     std::cout << std::endl;
 453     std::cout << "effective radius = " << std::setw(30) << effRad                   << std::endl;
 454     std::cout << "DAT file         = " << std::setw(30) << fnBetaListsDat           << std::endl;
 455     std::cout << "output file      = " << std::setw(30) << fnOut                    << std::endl;
 456     std::cout << "peptide atoms #  = " << std::setw(30) << index.size()             << std::endl;
 457     std::cout << "CA-atmos #       = " << std::setw(30) << indexCA.size()           << std::endl;
 458     std::cout << "Beta frames      = " << std::setw(30) << betaLists.size()         << std::endl;
 459     std::cout << "residue #        = " << std::setw(30) << aminoacidsIndex.size()   << std::endl;
 460     std::cout << "colors #         = " << std::setw(30) << colorsNames.size()       << std::endl;
 461
 462     /*
 463      * формально можно сделать:
 464      * найти соотношение между именами для углов и индексными номерами
 465      * заполнить std::vector< std::vector< size_t > > nameToIndex;
 466      * и в последствии считать:
 467      * a1 = frame[nameToIndex[0][1] + ... - ... - ...
 468      * b1 a2 b2 по аналогии только для [@<-[1, 2, 3]]
 469      *
 470      */
 471 }
 472
 473 void
 474 colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::TrajectoryAnalysisModuleData *pdata)
 475 {
 476     std::vector< gmx::RVec > trajectoryFrame;
 477     trajectoryFrame.resize(0);
 478     // считывания текущего фрейма траектории
 479     for (size_t i {0}; i < fr.natoms; ++i) {
 480         trajectoryFrame.push_back(fr.x[i]);
 481     }
 482     // подсчёт углов в красках
 483     std::vector< colorLocalAngles > colorFormation;
 484     colorsAnglesEvaluation(trajectoryFrame, colorsNames, colorFormation);
 485     // рассчёт положения относительно белка
 486     /*
 487      * формально можно совместить определение близости с белком и определение ближайших бэта-листов
 488      * для этого думаю нужно как-то соотнести или сделать соотношение между индексом и бэта-листами (т.е. хранить что бы за О(1) делать)
 489      *
 490      */
 491
 492
 493
 494     std::vector< bool >     colorsToPeptide;
 495     colorsToPeptide.resize(0);
 496     colorsToPeptide.resize(colorsNames.size(), false);
 497     for (size_t i = 0; i < colorsNames.size(); ++i) {
 498         colorsToPeptide[i] = isNearPeptide(pbc, trajectoryFrame, index, colorsNames[i], effRad);
 499     }
 500     // расчёт угла и среднего угла с ближайшими бета листами
 501     std::vector< std::vector< bool > > colorsToBeta;
 502     colorsToBeta.resize(0);
 503     colorsToBeta.resize(colorsNames.size());
 504     std::vector< std::vector< double > > colorsToBetaAngles;
 505     colorsToBetaAngles.resize(0);
 506     std::vector< gmx::RVec > betaListsRVecs;
 507     betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA);
 508     for (size_t i = 0; i < colorsToPeptide.size(); ++i) {
 509         if (colorsToPeptide[i]) {
 510             searchNearBetaLists(pbc, trajectoryFrame, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex, effRad);
 511         }
 512     }
 513     for (size_t i = 0; i < colorsToBeta.size(); ++i) {
 514         for (size_t j = 0; j < colorsToBeta[i].size(); ++j) {
 515             if (colorsToBeta[i][j]) {
 516                 computeAnglesColorsVsBeta(betaListsRVecs[j], colorFormation[i]);
 517             }
 518         }
 519     }
 520     // вывод в файл(ы) информацию для фрейма
 521     anglesFileDump(frnr, fnOut, colorsToPeptide, colorFormation);
 522     std::cout << "frame number = " << std::setw(8) << frnr << " trajectory time = " << std::setw(8) << fr.time << std::endl;
 523 }
 524
 525 void
 526 colorVec::finishAnalysis(int nframes)
 527 {
 528
 529 }
 530
 531 void
 532 colorVec::writeOutput()
 533 {
 534     /*
 535      *
 536      * frame #<current frame number>
 537      * 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>]
 538      *
 539      */
 540 }
 541
 542 /*! \brief
 543  * The main function for the analysis template.
 544  */
 545 int
 546 main(int argc, char *argv[])
 547 {
 548     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<colorVec>(argc, argv);
 549 }