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
  47 #include <iostream>
  48 #include <fstream>
  49 //#include <chrono>
  50 #include <cfloat>
  51 #include <omp.h>
  52 #include <vector>
  53 #include <cstdio>
  54 #include <iomanip>
  55
  56 // структура углов для одной краски
  57 struct colorLocalAngles {
  58     gmx::RVec   a1, a2;
  59     gmx::RVec   b1, b2;
  60     gmx::RVec   n1, n2;
  61     double      a12, b12, n12;
  62     std::vector< double > betaAngles;
  63 };
  64
  65 // хрен пойми почему, но захотелось рекурсивно сделать | можно сделать и через вайл
  66 long long returnBetaEnd(const std::string &currentLine, long long pos) {
  67     switch (currentLine[pos]) {
  68         case 'B' :
  69             returnBetaEnd(currentLine, pos + 1);
  70             //std::cout << "B" << std::endl;
  71             break;
  72         case 'E' :
  73             returnBetaEnd(currentLine, pos + 1);
  74             //std::cout << "E" << std::endl;
  75             break;
  76         default :
  77             return (pos - 1);
  78             break;
  79     }
  80     return (pos - 1); // формально логически лишняя строчка, но Qt ругается
  81 }
  82
  83 // функция парсинга одной строки
  84 void parseBetaListDATLine(const std::string &currentLine, std::vector< std::vector< unsigned int > > &localInputBL) {
  85     size_t                      equalCount = 0;
  86     std::vector< unsigned int > a;
  87     a.resize(0);
  88     //std::cout << "parsing " << currentLine.size() << std::endl;
  89     for (size_t i = 0; i < currentLine.size(); ++i) {
  90         //std::cout << i << std::endl;
  91         if (currentLine[i] == '=') { // подсчитываем число "пустых" символов
  92             ++equalCount;
  93         } else {
  94             long long temp = returnBetaEnd(currentLine, i);
  95             //std::cout << "|" << std::endl;
  96             if (temp - static_cast< long long >(i) > 3) {
  97                 localInputBL.push_back(a);
  98                 for (size_t j = i; j <= temp; ++j) {
  99                     localInputBL.back().push_back(j - equalCount);
 100                 }
 101                 i = temp;
 102             }
 103         }
 104     }
 105     //std::cout << std::endl;
 106 }
 107
 108 // функция нахождения бета-листов в структуре по файлу ДССП
 109 void betaListDigestion(const std::string &inputFile, std::vector< std::vector< std::vector< unsigned int > > > &inputBL) {
 110     inputBL.resize(0);
 111     //std::cout << "\t\tbetaListDigestion" << std::endl;
 112     std::ifstream   file(inputFile);
 113     std::string     line;
 114     getline(file, line); // считываем число в первой строке - кол-во осмысленных элементов в строках - нам не нужно
 115     getline(file, line);
 116     if (line.size() > 3) {
 117         do {
 118             //std::cout << "\t\tNew line = " << line << " | ";
 119             inputBL.resize(inputBL.size() + 1);
 120             parseBetaListDATLine(line, inputBL.back());
 121             line = "";
 122             getline(file, line);
 123         } while (line.size() > 3);
 124     } else {
 125         throw "DSSP DAT FILE IS EMPTY";
 126     }
 127     file.close();
 128 }
 129
 130 // функция выделения индексов для бэталистов
 131 inline void aminoacidsIndexation(const std::vector< size_t > &inputIndex, const gmx::TopologyInformation &top, std::vector< std::vector< size_t > > &aminoacidsIndex) {
 132     aminoacidsIndex.resize(0);
 133     for (size_t i = 0; i < inputIndex.size(); ++i) {
 134         std::cout << top.atoms()->atom[inputIndex[i]].resind << std::endl;
 135         aminoacidsIndex.resize(std::max(aminoacidsIndex.size(), static_cast< size_t >(top.atoms()->atom[inputIndex[i]].resind + 1)));
 136         aminoacidsIndex[top.atoms()->atom[inputIndex[i]].resind].push_back(inputIndex[i]);
 137     }
 138 }
 139
 140 // функция поиска RVec в кадре по имени->индексу
 141 gmx::RVec returnRVec(const std::vector< gmx::RVec > &frame, const std::vector< std::pair< std::string, size_t > > &colorsIndex, std::string toFind) {
 142     for (auto &i : colorsIndex) {
 143         if (i.first == toFind) {
 144             return frame[i.second];
 145         }
 146     }
 147     throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS";
 148 }
 149
 150 // функция векторного произведения двух RVec'ов
 151 inline void RVecVecMultiply(gmx::RVec &n, const gmx::RVec &a, const gmx::RVec &b) {
 152     n[0] = a[1] * b[2] - a[2] * b[1];
 153     n[1] = a[2] * b[0] - a[0] * b[2];
 154     n[2] = a[0] * b[1] - a[2] * b[0];
 155 }
 156
 157 // поиск угла между двумя RVec'ами
 158 inline double RVecAngle(const gmx::RVec &a, const gmx::RVec &b) {
 159     return std::acos((a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) / (a.norm() * b.norm()));
 160 }
 161
 162 // вычисление внутренних углов в краске
 163 void colorsAnglesEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex, std::vector< colorLocalAngles > &colorStruct) {
 164     colorStruct.resize(colorsIndex.size());
 165     #pragma omp parallel for ordered schedule(dynamic)
 166     for (size_t i = 0; i < colorsIndex.size(); ++i) {
 167         colorStruct[i].betaAngles.resize(0);
 168         // "левый" блок
 169         colorStruct[i].a1 = returnRVec(frame, colorsIndex[i], "CAF") - returnRVec(frame, colorsIndex[i], "CAJ") + returnRVec(frame, colorsIndex[i], "CAK") - returnRVec(frame, colorsIndex[i], "CAO");
 170         colorStruct[i].a1 /= colorStruct[i].a1.norm();
 171         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");
 172         colorStruct[i].b1 /= colorStruct[i].b1.norm();
 173         RVecVecMultiply(colorStruct[i].n1, colorStruct[i].a1, colorStruct[i].b1);
 174         colorStruct[i].n1 /= colorStruct[i].n1.norm();
 175         // "правый" блок
 176         colorStruct[i].a2 = returnRVec(frame, colorsIndex[i], "CAB") - returnRVec(frame, colorsIndex[i], "CBQ") + returnRVec(frame, colorsIndex[i], "CBP") - returnRVec(frame, colorsIndex[i], "CBL");
 177         colorStruct[i].a2 /= colorStruct[i].a2.norm();
 178         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");
 179         colorStruct[i].b2 /= colorStruct[i].b2.norm();
 180         RVecVecMultiply(colorStruct[i].n2, colorStruct[i].a2, colorStruct[i].b2);
 181         colorStruct[i].n2 *= -1; // для "сонаправленности" векторов нормали
 182         colorStruct[i].n2 /= colorStruct[i].n2.norm();
 183         colorStruct[i].a12 = RVecAngle(colorStruct[i].a1, colorStruct[i].a2);
 184         colorStruct[i].b12 = RVecAngle(colorStruct[i].b1, colorStruct[i].b2);
 185         colorStruct[i].n12 = RVecAngle(colorStruct[i].n1, colorStruct[i].n2);
 186     }
 187     #pragma omp barrier
 188 }
 189
 190 // вычисление направляющего вектора в бэта-листе
 191 inline void betaListsRVecsEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< std::vector< unsigned int > > &inputBetaLists, std::vector< gmx::RVec > &temp, const std::vector< size_t > &inputCA) {
 192     temp.resize(0);
 193     gmx::RVec tempA;
 194     for (const auto &i : inputBetaLists) {
 195         tempA = frame[inputCA[i[i.size() - 1]]] + frame[inputCA[i[i.size() - 2]]] - frame[inputCA[i[1]]] - frame[inputCA[i[0]]];
 196         tempA /= tempA.norm();
 197         temp.push_back(tempA);
 198     }
 199 }
 200
 201 // определение близко ли к белку находится краска
 202 bool isNearPeptide(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > inputFrame,
 203                           gmx::AnalysisNeighborhood &nbhood, const std::vector< size_t > &inputIndex,
 204                           const std::vector< std::pair< std::string, size_t > > &inputColor) {
 205     gmx::AnalysisNeighborhoodSearch      nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
 206     gmx::AnalysisNeighborhoodPair        pair;
 207     for (const auto &i : inputColor) {
 208         while (nbsearch.startPairSearch(inputFrame[i.second].as_vec()).findNextPair(&pair)) {
 209             for (const auto &j : inputIndex) {
 210                 if (pair.testIndex() == j) {
 211                     return true;
 212                 }
 213             }
 214         }
 215     }
 216     return false;
 217 }
 218
 219 // поиск ближайших к краске бэта-листов
 220 inline void searchNearBetaLists(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > inputFrame,
 221                                 gmx::AnalysisNeighborhood &nbhood,
 222                                 const std::vector< std::vector< unsigned int > > inputBLists,
 223                                 const std::vector< std::pair< std::string, size_t > > &inputColor,
 224                                 std::vector< bool > &outputList, const std::vector< std::vector< size_t > > &inputAminoacids) {
 225     outputList.resize(0);
 226     outputList.resize(inputBLists.size(), false);
 227     gmx::AnalysisNeighborhoodSearch      nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
 228     gmx::AnalysisNeighborhoodPair        pair;
 229     for (const auto &i : inputColor) {
 230         while (nbsearch.startPairSearch(inputFrame[i.second].as_vec()).findNextPair(&pair)) {
 231             for (size_t j = 0; j < inputBLists.size(); ++j) {
 232                 for (size_t k = 0; k < inputBLists[j].size(); ++k) {
 233                     for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
 234                         if (pair.testIndex() == inputAminoacids[inputBLists[j][k]][m]) {
 235                             outputList[j] = true;
 236                         }
 237                     }
 238                 }
 239             }
 240         }
 241     }
 242
 243 }
 244
 245 // определение углов между краской и направляющим вектором бэта-листа
 246 inline void computeAnglesColorsVsBeta(const gmx::RVec &inputBetaRVec, colorLocalAngles &colorFormation) {
 247     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a1));
 248     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b1));
 249     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n1));
 250     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a2));
 251     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b2));
 252     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n2));
 253 }
 254
 255 // функция записи в файл значений углов для кадра
 256 void anglesFileDump(const int frameNum, const std::string &output, const std::vector< bool > &toPeptide, const std::vector< colorLocalAngles > &colorFormation) {
 257     std::ofstream   file(output);
 258     int temp = 0;
 259     std::vector< double > betaTemp;
 260     betaTemp.resize(0);
 261     file << "frame =" << std::setw(8) << frameNum << std::endl;
 262     for (size_t i = 0; i < colorFormation.size(); ++i) {
 263         file << "color #" << std::setw(3) << i;
 264         if (toPeptide[i]) {
 265             file << std::setw(4) << "yes";
 266         } else {
 267             file << std::setw(4) << "no";
 268         }
 269         file << std::setw(7) << std::setprecision(2) << colorFormation[i].a12 << colorFormation[i].b12 << colorFormation[i].n12;
 270         temp = colorFormation[i].betaAngles.size() / 6;
 271         if (temp == 0) {
 272             file << std::setw(4) << 0;
 273         } else {
 274             betaTemp.resize(6, 0.); // magic number, meh
 275             for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
 276                 betaTemp[j % 6] += colorFormation[i].betaAngles[j];
 277             }
 278             for (size_t j = 0; j < betaTemp.size(); ++j) {
 279                 betaTemp[j] /= temp;
 280             }
 281             file << std::setw(4) << temp;
 282             for (size_t j = 0; j < betaTemp.size(); ++j) {
 283                 file << std::setw(7) << std::setprecision(2) << betaTemp[j];
 284             }
 285             for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
 286                 file << std::setw(7) << std::setprecision(2) << colorFormation[i].betaAngles[j];
 287             }
 288         }
 289     }
 290     file.close();
 291 }
 292
 293 /*! \brief
 294  * \ingroup module_trajectoryanalysis
 295  */
 296 class colorVec : public gmx::TrajectoryAnalysisModule
 297 {
 298     public:
 299
 300         colorVec();
 301         virtual ~colorVec();
 302
 303         //! Set the options and setting
 304         virtual void initOptions(gmx::IOptionsContainer          *options,
 305                                  gmx::TrajectoryAnalysisSettings *settings);
 306
 307         //! First routine called by the analysis framework
 308         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 309         virtual void initAnalysis(const gmx::TrajectoryAnalysisSettings &settings,
 310                                   const gmx::TopologyInformation        &top);
 311
 312         //! Call for each frame of the trajectory
 313         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 314         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 315                                   gmx::TrajectoryAnalysisModuleData *pdata);
 316
 317         //! Last routine called by the analysis framework
 318         // virtual void finishAnalysis(t_pbc *pbc);
 319         virtual void finishAnalysis(int nframes);
 320
 321         //! Routine to write output, that is additional over the built-in
 322         virtual void writeOutput();
 323
 324     private:
 325
 326         gmx::SelectionList                                              sel_;
 327         std::string                                                     fnOut           {"name"};   // selectable
 328         std::string                                                     fnBetaListsDat;             // selectable
 329         float                                                           effRad          {0.8};      // selectable
 330         std::vector< size_t >                                           index;
 331         std::vector< size_t >                                           indexCA;
 332         std::vector< std::vector< size_t > >                            aminoacidsIndex;
 333         std::vector< std::vector< std::pair< std::string, size_t > > >  colorsNames;
 334         std::vector< std::vector< std::vector< unsigned int > > >       betaLists;
 335         gmx::AnalysisNeighborhood                                       nb_;
 336
 337         // Copy and assign disallowed by base.
 338 };
 339
 340 colorVec::colorVec(): TrajectoryAnalysisModule()
 341 {
 342 }
 343
 344 colorVec::~colorVec()
 345 {
 346 }
 347
 348 /*
 349  * ndx:
 350  * [peptide]
 351  * [CA]
 352  * [color_{i}]
 353  *
 354  */
 355
 356 void
 357 colorVec::initOptions(  gmx::IOptionsContainer          *options,
 358                         gmx::TrajectoryAnalysisSettings *settings)
 359 {
 360     static const char *const desc[] = {
 361         "[THISMODULE] to be done"
 362     };
 363     // Add the descriptive text (program help text) to the options
 364     settings->setHelpText(desc);
 365     // Add option for selection list
 366     options->addOption(gmx::SelectionOption("sel")
 367                             .storeVector(&sel_)
 368                             .required().dynamicMask().multiValue()
 369                             .description("select pepride and colors / -sf"));
 370     // Add option for input file names
 371     options->addOption(gmx::StringOption("dat")
 372                             .store(&fnBetaListsDat)
 373                             .description("a file to make dynamic beta lists"));
 374     // Add option for output file name
 375     options->addOption(gmx::StringOption("out")
 376                             .store(&fnOut)
 377                             .description("Index file for the algorithm output."));
 378     // Add option for effRad constant
 379     options->addOption(gmx::FloatOption("efRad")
 380                             .store(&effRad)
 381                             .description("max distance from colors to peptide in nm to consider to be \"near\""));
 382     // Control input settings
 383     settings->setFlags(gmx::TrajectoryAnalysisSettings::efNoUserPBC);
 384     settings->setFlag(gmx::TrajectoryAnalysisSettings::efUseTopX);
 385     //settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 386     settings->setPBC(true);
 387 }
 388
 389 void
 390 colorVec::initAnalysis( const gmx::TrajectoryAnalysisSettings   &settings,
 391                         const gmx::TopologyInformation          &top)
 392 {
 393     // считывание индекса
 394     std::cout << "\tReading indexes." << std::endl;
 395     index.resize(0);
 396     for (gmx::ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ai++) {
 397         index.push_back(static_cast< size_t >(*ai));
 398     }
 399     // считывание индекса
 400     indexCA.resize(0);
 401     for (gmx::ArrayRef< const int >::iterator ai {sel_[1].atomIndices().begin()}; (ai < sel_[2].atomIndices().end()); ai++) {
 402         index.push_back(static_cast< size_t >(*ai));
 403     }
 404     // считывание красок
 405     colorsNames.resize(sel_.size() - 2);
 406     for (size_t i = 2; i < sel_.size(); ++i) {
 407         for (gmx::ArrayRef< const int >::iterator ai {sel_[i].atomIndices().begin()}; (ai < sel_[i].atomIndices().end()); ai++) {
 408             colorsNames[i - 2].push_back(std::make_pair(*(top.atoms()->resinfo[top.atoms()->atom[*ai].resind].name), static_cast< size_t >(*ai)));
 409         }
 410     }
 411     std::cout << "\tBeta-lists digestion." << std::endl;
 412     // разбор dat файла для создания бэта листов
 413     betaListDigestion(fnBetaListsDat, betaLists);
 414     // разбор топологии для индексации бета листов
 415     std::cout << "\tAminoacids indexation." << std::endl;
 416     aminoacidsIndexation(index, top, aminoacidsIndex);
 417     // задание радиуса рассматриваемых соседей
 418     nb_.setCutoff(effRad);
 419     /*
 420      * формально можно сделать:
 421      * найти соотношение между именами для углов и индексными номерами
 422      * заполнить std::vector< std::vector< size_t > > nameToIndex;
 423      * и в последствии считать:
 424      * a1 = frame[nameToIndex[0][1] + ... - ... - ...
 425      * b1 a2 b2 по аналогии только для [@<-[1, 2, 3]]
 426      *
 427      */
 428 }
 429
 430 void
 431 colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::TrajectoryAnalysisModuleData *pdata)
 432 {
 433     std::cout << "\tFrame analisis start ..." << std::endl;
 434     std::vector< gmx::RVec > trajectoryFrame;
 435     trajectoryFrame.resize(0);
 436     // считывания текущего фрейма траектории
 437     for (size_t i {0}; i < index.size(); ++i) {
 438         trajectoryFrame.push_back(fr.x[index[i]]);
 439     }
 440     // подсчёт углов в красках
 441     std::vector< colorLocalAngles > colorFormation;
 442     std::cout << "\t\tColors' angles evaluation." << std::endl;
 443     colorsAnglesEvaluation(trajectoryFrame, colorsNames, colorFormation);
 444     // рассчёт положения относительно белка
 445     /*
 446      * формально можно совместить определение близости с белком и определение ближайших бэта-листов
 447      * для этого думаю нужно как-то соотнести или сделать соотношение между индексом и бэта-листами (т.е. хранить что бы за О(1) делать)
 448      *
 449      */
 450
 451     std::vector< bool >     colorsToPeptide;
 452     colorsToPeptide.resize(0);
 453     colorsToPeptide.resize(colorsNames.size(), false);
 454     std::cout << "\t\tWhich colors are \"near\" the pepride mass." << std::endl;
 455     for (size_t i = 0; i < colorsNames.size(); ++i) {
 456         colorsToPeptide[i] = isNearPeptide(fr, pbc, trajectoryFrame, nb_, index, colorsNames[i]);
 457     }
 458     // расчёт угла и среднего угла с ближайшими бета листами
 459     std::vector< std::vector< bool > > colorsToBeta;
 460     colorsToBeta.resize(0);
 461     colorsToBeta.resize(colorsNames.size());
 462     std::vector< std::vector< double > > colorsToBetaAngles;
 463     colorsToBetaAngles.resize(0);
 464     std::vector< gmx::RVec > betaListsRVecs;
 465     std::cout << "\t\tBeta-lists' RVecs search." << std::endl;
 466     betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA);
 467     std::cout << "\t\tSearching nearby beta-lists." << std::endl;
 468     for (size_t i = 0; i < colorsToPeptide.size(); ++i) {
 469         if (colorsToPeptide[i]) {
 470             searchNearBetaLists(fr, pbc, trajectoryFrame, nb_, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex);
 471         }
 472     }
 473     std::cout << "\t\tComputing angles colors vs betas." << std::endl;
 474     for (size_t i = 0; i < colorsToBeta.size(); ++i) {
 475         for (size_t j = 0; j < colorsToBeta[i].size(); ++j) {
 476             if (colorsToBeta[i][j]) {
 477                 computeAnglesColorsVsBeta(betaListsRVecs[j], colorFormation[i]);
 478             }
 479         }
 480     }
 481     // вывод в файл(ы) информацию для фрейма
 482     std::cout << "\t\tDumping data." << std::endl;
 483     anglesFileDump(frnr, fnOut, colorsToPeptide, colorFormation);
 484     std::cout << "\tFrame analysed." << std::endl;
 485 }
 486
 487 void
 488 colorVec::finishAnalysis(int nframes)
 489 {
 490
 491 }
 492
 493 void
 494 colorVec::writeOutput()
 495 {
 496     /*
 497      *
 498      * frame #<current frame number>
 499      * 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>]
 500      *
 501      */
 502     std::cout << "\n\t colorAngles finished successfully" << std::endl;
 503 }
 504
 505 /*! \brief
 506  * The main function for the analysis template.
 507  */
 508 int
 509 main(int argc, char *argv[])
 510 {
 511     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<colorVec>(argc, argv);
 512 }