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 struct colorAnglesPairs
  67 {
  68     std::vector< std::pair<size_t, size_t> > a1, a2;
  69     std::vector< std::pair<size_t, size_t> > b1, b2;
  70 };
  71
  72 struct massChargePair {
  73     double m{0.}, ch{0.};
  74 };
  75
  76 // функция парсинга одной строки
  77 void parseBetaListDATLine(const std::string                             &currentLine,
  78                           std::vector< std::vector< unsigned int > >    &localInputBL) {
  79     size_t                      equalCount {0};
  80     std::vector< unsigned int > a;
  81     a.resize(0);
  82     for (size_t i {0}; i < currentLine.size(); ++i) {
  83         if (currentLine[i] == '=') { // подсчитываем число "пустых" символов
  84             ++equalCount;
  85         } else {
  86             size_t temp {i};
  87             for (size_t j {i + 1}; j < currentLine.size(); ++j) {
  88                 if (currentLine[j] == 'B' || currentLine[j] == 'E') {
  89                     ++temp;
  90                 } else {
  91                     break;
  92                 }
  93             }
  94             if (temp - i > 3) {
  95                 localInputBL.push_back(a);
  96                 for (size_t j {i}; j <= temp; ++j) {
  97                     localInputBL.back().push_back(j - equalCount);
  98                 }
  99                 i = temp;
 100             }
 101         }
 102     }
 103 }
 104
 105 // функция нахождения бета-листов в структуре по файлу ДССП
 106 void betaListDigestion(const std::string                                            &inputFile,
 107                        std::vector< std::vector< std::vector< unsigned int > > >    &inputBL) {
 108     inputBL.resize(0);
 109     std::ifstream   file(inputFile);
 110     std::string     line;
 111     getline(file, line); // считываем число в первой строке - кол-во осмысленных элементов в строках - нам не нужно
 112     getline(file, line);
 113     if (line.size() > 3) {
 114         do {
 115             inputBL.resize(inputBL.size() + 1);
 116             parseBetaListDATLine(line, inputBL.back());
 117             line = "";
 118             getline(file, line);
 119         } while (line.size() > 3);
 120     } else {
 121         throw "DSSP DAT FILE IS EMPTY";
 122     }
 123     file.close();
 124 }
 125
 126 // функция выделения индексов для бэталистов
 127 inline void aminoacidsIndexation(const std::vector< size_t >            &inputIndex,
 128                                  const gmx::TopologyInformation         &top,
 129                                  std::vector< std::vector< size_t > >   &aminoacidsIndex) {
 130     aminoacidsIndex.resize(0);
 131     for (size_t i {0}; i < inputIndex.size(); ++i) {
 132         aminoacidsIndex.resize(std::max(aminoacidsIndex.size(), static_cast< size_t >(top.atoms()->atom[inputIndex[i]].resind + 1)));
 133         aminoacidsIndex[top.atoms()->atom[inputIndex[i]].resind].push_back(inputIndex[i]);
 134     }
 135 }
 136
 137 // функция векторного произведения двух RVec'ов
 138 inline void RVecVecMultiply(gmx::RVec &n, const gmx::RVec &a, const gmx::RVec &b) {
 139     n[0] = a[1] * b[2] - a[2] * b[1];
 140     n[1] = a[2] * b[0] - a[0] * b[2];
 141     n[2] = a[0] * b[1] - a[2] * b[0];
 142 }
 143
 144 // поиск угла между двумя RVec'ами
 145 inline double RVecAngle(const gmx::RVec &a, const gmx::RVec &b) {
 146     return std::acos((a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) / (a.norm() * b.norm())) * 180.0 / 3.14159265;
 147 }
 148
 149 // вычисление внутренних углов в краске
 150 void colorsAnglesEvaluation(const std::vector< gmx::RVec >          &frame,
 151                             const std::vector< colorAnglesPairs >   &inputPairs,
 152                             std::vector< colorLocalAngles >         &colorStruct) {
 153     colorStruct.resize(0);
 154     colorStruct.resize(inputPairs.size());
 155     #pragma omp parallel for ordered schedule(dynamic)
 156     for (size_t i = 0; i < inputPairs.size(); ++i) {
 157         colorStruct[i].betaAngles.resize(0);
 158         // "левый" блок
 159         for (size_t j {0}; j < inputPairs[i].a1.size(); ++j) {
 160             colorStruct[i].a1 += frame[inputPairs[i].a1[j].first] - frame[inputPairs[i].a1[j].second];
 161         }
 162         for (size_t j {0}; j < inputPairs[i].b1.size(); ++j) {
 163             colorStruct[i].b1 += frame[inputPairs[i].b1[j].first] - frame[inputPairs[i].b1[j].second];
 164         }
 165         RVecVecMultiply(colorStruct[i].n1, colorStruct[i].a1, colorStruct[i].b1);
 166         // "правый" блок
 167         for (size_t j {0}; j < inputPairs[i].a2.size(); ++j) {
 168             colorStruct[i].a2 += frame[inputPairs[i].a2[j].first] - frame[inputPairs[i].a2[j].second];
 169         }
 170         for (size_t j {0}; j < inputPairs[i].b2.size(); ++j) {
 171             colorStruct[i].b2 += frame[inputPairs[i].b2[j].first] - frame[inputPairs[i].b2[j].second];
 172         }
 173         RVecVecMultiply(colorStruct[i].n2, colorStruct[i].a2, colorStruct[i].b2);
 174         colorStruct[i].n2 *= -1; // для "сонаправленности" векторов нормали
 175         colorStruct[i].a12 = RVecAngle(colorStruct[i].a1, colorStruct[i].a2);
 176         colorStruct[i].b12 = RVecAngle(colorStruct[i].b1, colorStruct[i].b2);
 177         colorStruct[i].n12 = RVecAngle(colorStruct[i].n1, colorStruct[i].n2);
 178     }
 179     #pragma omp barrier
 180 }
 181
 182 // функция поиска индекса в краске по имени->индексу
 183 size_t returnRVecIndex(const std::vector< std::pair< std::string, size_t > > &colorsIndex, const std::string &toFind) {
 184     for (auto &i : colorsIndex) {
 185         if (i.first == toFind) {
 186             return i.second;
 187         }
 188     }
 189     throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS";
 190 }
 191
 192 // заполнение пресетов для поиска средних векторов в краске
 193 void colorAnglesPairsEvaluation(const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex,
 194                                 std::vector< colorAnglesPairs > &output) {
 195     for (size_t i {0}; i < output.size(); ++i) {
 196         output[i].a1.resize(0);
 197         output[i].a1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAF"), returnRVecIndex(colorsIndex[i], "CAJ")));
 198         output[i].a1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAK"), returnRVecIndex(colorsIndex[i], "CAO")));
 199         output[i].b1.resize(0);
 200         output[i].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAO"), returnRVecIndex(colorsIndex[i], "CAJ")));
 201         output[i].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAM"), returnRVecIndex(colorsIndex[i], "CAH")));
 202         output[i].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAK"), returnRVecIndex(colorsIndex[i], "CAF")));
 203         output[i].a2.resize(0);
 204         output[i].a2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAB"), returnRVecIndex(colorsIndex[i], "CBQ")));
 205         output[i].a2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CBP"), returnRVecIndex(colorsIndex[i], "CBL")));
 206         output[i].b2.resize(0);
 207         output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CBL"), returnRVecIndex(colorsIndex[i], "CBQ")));
 208         output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CBN"), returnRVecIndex(colorsIndex[i], "CBS")));
 209         output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CBP"), returnRVecIndex(colorsIndex[i], "CAB")));
 210     }
 211 }
 212
 213 // вычисление направляющего вектора в бэта-листе
 214 inline void betaListsRVecsEvaluation(const std::vector< gmx::RVec > &frame,
 215                                      const std::vector< std::vector< unsigned int > > &inputBetaLists,
 216                                      std::vector< gmx::RVec >       &temp,
 217                                      const std::vector< size_t >    &inputCA) {
 218     temp.resize(0);
 219     gmx::RVec tempA;
 220     for (const auto &i : inputBetaLists) {
 221         tempA = frame[inputCA[i[i.size() - 1]]] + frame[inputCA[i[i.size() - 2]]] - frame[inputCA[i[1]]] - frame[inputCA[i[0]]];
 222         tempA /= tempA.norm();
 223         temp.push_back(tempA);
 224     }
 225 }
 226
 227 // определение близко ли к белку находится краска
 228 bool isNearPeptide( const t_pbc *inputPBC,
 229                     const std::vector< gmx::RVec > &inputFrame,
 230                     const std::vector< size_t > &inputIndex,
 231                     const std::vector< std::pair< std::string, size_t > > &inputColor, const double cutOff) {
 232     gmx::RVec tempRVec;
 233     for (size_t i {0}; i < inputColor.size(); ++i) {
 234         for (size_t j {0}; j < inputIndex.size(); ++j) {
 235             pbc_dx(inputPBC, inputFrame[inputIndex[j]], inputFrame[inputColor[i].second], tempRVec);
 236             if (tempRVec.norm() <= cutOff) {
 237                 return true;
 238             }
 239         }
 240     }
 241     return false;
 242 }
 243
 244 // поиск ближайших к краске бэта-листов
 245 inline void searchNearBetaLists(const t_pbc *inputPBC,
 246                                 const std::vector< gmx::RVec > &inputFrame, const std::vector< std::vector< unsigned int > > &inputBLists,
 247                                 const std::vector< std::pair< std::string, size_t > > &inputColor, std::vector< bool > &outputList,
 248                                 const std::vector< std::vector< size_t > > &inputAminoacids, const double cutOff) {
 249     outputList.resize(0);
 250     outputList.resize(inputBLists.size(), false);
 251     gmx::RVec tempRVec;
 252     for (size_t i {0}; i < inputColor.size(); ++i) {
 253         for (size_t j {0}; j < inputBLists.size(); ++j) {
 254             for (size_t k {0}; k < inputBLists[j].size(); ++k) {
 255                 for (size_t m {0}; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
 256                     pbc_dx(inputPBC, inputFrame[inputAminoacids[inputBLists[j][k]][m]], inputFrame[inputColor[i].second], tempRVec);
 257                     if (tempRVec.norm() <= cutOff) {
 258                         outputList[j] = true;
 259                         break;
 260                     }
 261                 }
 262             }
 263         }
 264     }
 265 }
 266
 267 // определение углов между краской и направляющим вектором бэта-листа
 268 inline void computeAnglesColorsVsBeta(const gmx::RVec &inputBetaRVec, colorLocalAngles &colorFormation) {
 269     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a1));
 270     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b1));
 271     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n1));
 272     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a2));
 273     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b2));
 274     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n2));
 275 }
 276
 277 // поиск стаков красок
 278 void colorsStackingSearch(const t_pbc *inputPBC, const std::vector< gmx::RVec > &inputFrame,
 279                           const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex,
 280                           const float radius, const float epsi,
 281                           std::vector< std::vector< int > > &output) {
 282     std::vector< std::vector< int > > outputTemp;
 283     outputTemp.resize(0);
 284     outputTemp.resize(colorsIndex.size());
 285     for (size_t i {0}; i < colorsIndex.size(); ++i) {
 286         outputTemp[i].resize(0);
 287         outputTemp[i].resize(colorsIndex.size(), 0);
 288     }
 289     gmx::RVec temp;
 290     #pragma omp parallel for ordered schedule(dynamic)
 291     for (size_t i = 0; i < colorsIndex.size(); ++i) {
 292         for (size_t j {i + 1}; j < colorsIndex.size(); ++j) {
 293             for (size_t k1 {0}; k1 < colorsIndex[i].size(); ++k1) {
 294                 for (size_t k2 {0}; k2 < colorsIndex[j].size(); ++k2) {
 295                     pbc_dx(inputPBC, inputFrame[colorsIndex[i][k1].second], inputFrame[colorsIndex[j][k2].second], temp);
 296                     if (temp.norm() <= radius) {
 297                         ++outputTemp[i][j];
 298                         ++outputTemp[j][i];
 299                     }
 300                 }
 301             }
 302         }
 303     }
 304     #pragma omp barrier
 305     output.resize(0);
 306     output.resize(colorsIndex.size());
 307     for (size_t i {0}; i < colorsIndex.size(); ++i) {
 308         output[i].resize(0);
 309     }
 310     for (size_t i {0}; i < outputTemp.size(); ++i) {
 311         for (size_t j {i + 1}; j < outputTemp[i].size(); ++j) {
 312             if (static_cast< float >(outputTemp[i][j]) >= static_cast< float >(colorsIndex.size()) * epsi) {
 313                 output[i].push_back(j);
 314                 output[j].push_back(i);
 315             }
 316         }
 317     }
 318 }
 319
 320 // запоминаем массу/заряды атомов красок
 321 void colorMassChargesEval(  const std::vector< std::pair< std::string, size_t > >   &color,
 322                             const gmx::TopologyInformation                          &top,
 323                             std::vector< std::pair< size_t, massChargePair > >      &outputIndexMassChargePair) {
 324     outputIndexMassChargePair.resize(0);
 325     massChargePair temp;
 326     for (size_t i {0}; i < color.size(); ++i) {
 327         temp.m = top.atoms()->atom[color[i].second].m;
 328         temp.ch = top.atoms()->atom[color[i].second].q;
 329         outputIndexMassChargePair.push_back(std::make_pair(color[i].second, temp));
 330     }
 331 }
 332
 333 // вычисляем дипольный вектор краски
 334 gmx::RVec colorDipoleEval(const std::vector< std::pair< size_t, massChargePair > >  &inputIndexMassChargePair,
 335                           const std::vector< gmx::RVec >                            &inputFrame) {
 336     gmx::RVec temp {0., 0., 0.};
 337     double temp2 {0.};
 338     for (auto &i : inputIndexMassChargePair) {
 339         temp += inputFrame[i.first] * static_cast< float >(i.second.m);
 340         temp2 += i.second.m;
 341     }
 342     temp /= static_cast< float >(temp2);
 343     return temp;
 344 }
 345
 346 // функция записи в файл значений углов для кадра
 347 void anglesFileDump(const int frameNum, const std::string       &output,
 348                     const std::vector< bool >                   &toPeptide,
 349                     const std::vector< colorLocalAngles >       &colorFormation,
 350                     const std::vector< std::vector< int > >     &inputStack,
 351                     const std::vector< std::vector< double > >  &dipoleBetaAngles) {
 352     std::ofstream   file(output, std::ofstream::app);
 353     file << "frame =" << std::setw(8) << frameNum << std::endl;
 354     for (size_t i {0}; i < colorFormation.size(); ++i) {
 355         file << "color #" << std::setw(3) << i;
 356         if (toPeptide[i]) {
 357             file << std::setw(4) << "yes";
 358         } else {
 359             file << std::setw(4) << "no";
 360         }
 361         file << std::fixed;
 362         file << std::setw(8) << std::setprecision(2) << colorFormation[i].a12;
 363         file << std::setw(8) << std::setprecision(2) << colorFormation[i].b12;
 364         file << std::setw(8) << std::setprecision(2) << colorFormation[i].n12;
 365         if (colorFormation[i].betaAngles.size() == 0) {
 366             file << std::setw(4) << 0 << " [] ";
 367         } else {
 368             int temp = colorFormation[i].betaAngles.size() / 6;
 369             std::vector< double > betaTemp;
 370             betaTemp.resize(0);
 371             betaTemp.resize(6, 0.); // magic number, meh
 372             for (size_t j {0}; j < colorFormation[i].betaAngles.size(); ++j) {
 373                 betaTemp[j % 6] += colorFormation[i].betaAngles[j];
 374             }
 375             for (size_t j {0}; j < betaTemp.size(); ++j) {
 376                 betaTemp[j] /= static_cast< double >(temp);
 377             }
 378             file << std::setw(4) << temp << " [ ";
 379             for (size_t j {0}; j < betaTemp.size(); ++j) {
 380                 file << std::setw(8) << std::setprecision(2) << betaTemp[j];
 381             }
 382             for (size_t j {0}; j < colorFormation[i].betaAngles.size(); ++j) {
 383                 file << std::setw(8) << std::setprecision(2) << colorFormation[i].betaAngles[j];
 384             }
 385             file << " ] ";
 386         }
 387         if (inputStack[i].size() == 0) {
 388             file << std::setw(4) << "no []";
 389         } else {
 390             file << std::setw(4) << "yes [ ";
 391             file << std::setw(4) << inputStack[i].size();
 392             for (size_t j {0}; j < inputStack[i].size(); ++j) {
 393                 file << std::setw(4) << inputStack[i][j];
 394             }
 395             file << " ] ";
 396         }
 397         if (dipoleBetaAngles[i].size() == 0) {
 398             file << "no []";
 399         } else {
 400             double tempDBangle {0.};
 401             for (size_t j {0}; j < dipoleBetaAngles[i].size(); ++j) {
 402                 tempDBangle += dipoleBetaAngles[i][j];
 403             }
 404             file << std::setw(8) << std::setprecision(2) << tempDBangle / dipoleBetaAngles[i].size();
 405             for (size_t j {0}; j < dipoleBetaAngles[i].size(); ++j) {
 406                 file << std::setw(8) << std::setprecision(2) << dipoleBetaAngles[i][j];
 407             }
 408         }
 409         file << std::endl;
 410     }
 411     file.close();
 412 }
 413
 414 /*! \brief
 415  * \ingroup module_trajectoryanalysis
 416  */
 417 class colorVec : public gmx::TrajectoryAnalysisModule
 418 {
 419     public:
 420
 421         colorVec();
 422         virtual ~colorVec();
 423
 424         //! Set the options and setting
 425         virtual void initOptions(gmx::IOptionsContainer          *options,
 426                                  gmx::TrajectoryAnalysisSettings *settings);
 427
 428         //! First routine called by the analysis framework
 429         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 430         virtual void initAnalysis(const gmx::TrajectoryAnalysisSettings &settings,
 431                                   const gmx::TopologyInformation        &top);
 432
 433         //! Call for each frame of the trajectory
 434         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 435         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 436                                   gmx::TrajectoryAnalysisModuleData *pdata);
 437
 438         //! Last routine called by the analysis framework
 439         // virtual void finishAnalysis(t_pbc *pbc);
 440         virtual void finishAnalysis(int nframes);
 441
 442         //! Routine to write output, that is additional over the built-in
 443         virtual void writeOutput();
 444
 445     private:
 446
 447         gmx::SelectionList                                                  sel_;
 448         std::string                                                         fnOut           {"name"};   // selectable
 449         std::string                                                         fnBetaListsDat;             // selectable
 450         float                                                               effRad          {0.4};      // selectable
 451         float                                                               stackRad        {0.5};      // selectable
 452         float                                                               stackCoef       {0.66};     // selectable
 453         std::vector< size_t >                                               index;
 454         std::vector< size_t >                                               indexCA;
 455         std::vector< std::vector< size_t > >                                aminoacidsIndex;
 456         std::vector< std::vector< std::pair< std::string, size_t > > >      colorsNames;
 457         std::vector< colorAnglesPairs >                                     colorsVectorsIndex;
 458         std::vector< colorLocalAngles >                                     colorsAngles;
 459         std::vector< std::vector< std::pair< size_t, massChargePair > > >   colorsMassCharges;
 460         std::vector< std::vector< std::vector< unsigned int > > >           betaLists;
 461
 462         std::vector< gmx::RVec >                                            trajectoryFrame;
 463
 464         // Copy and assign disallowed by base.
 465 };
 466
 467 colorVec::colorVec(): TrajectoryAnalysisModule()
 468 {
 469 }
 470
 471 colorVec::~colorVec()
 472 {
 473 }
 474
 475 /*
 476  * ndx:
 477  * [peptide]
 478  * [CA]
 479  * [color_{i}]
 480  *
 481  */
 482
 483 void
 484 colorVec::initOptions(  gmx::IOptionsContainer          *options,
 485                         gmx::TrajectoryAnalysisSettings *settings)
 486 {
 487     static const char *const desc[] = {
 488         "[THISMODULE] to be done"
 489     };
 490     // Add the descriptive text (program help text) to the options
 491     settings->setHelpText(desc);
 492     // Add option for selection list
 493     options->addOption(gmx::SelectionOption("sel")
 494                             .storeVector(&sel_)
 495                             .required().dynamicMask().multiValue()
 496                             .description("select pepride and colors / -sf"));
 497     // Add option for input file names
 498     options->addOption(gmx::StringOption("dat")
 499                             .store(&fnBetaListsDat)
 500                             .description("a file to make dynamic beta lists"));
 501     // Add option for output file name
 502     options->addOption(gmx::StringOption("out")
 503                             .store(&fnOut)
 504                             .description("Index file for the algorithm output."));
 505     // Add option for effRad constant
 506     options->addOption(gmx::FloatOption("efRad")
 507                             .store(&effRad)
 508                             .description("max distance from colors to peptide in nm to consider to be \"near\" - default == 0.4"));
 509     // Add option for effRad constant
 510     options->addOption(gmx::FloatOption("stackRad")
 511                             .store(&stackRad)
 512                             .description("max distance from color to color in nm to consider to be \"stacked\" - default == 0.5"));
 513     // Add option for effRad constant
 514     options->addOption(gmx::FloatOption("stackCoef")
 515                             .store(&stackCoef)
 516                             .description("min % (1 == 100%) stacked mass of a colorto be considered \"stacked\" - default == 0.66"));
 517     // Control input settings
 518     settings->setFlags(gmx::TrajectoryAnalysisSettings::efNoUserPBC);
 519     settings->setFlag(gmx::TrajectoryAnalysisSettings::efUseTopX);
 520     //settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 521     settings->setPBC(true);
 522 }
 523
 524 void
 525 colorVec::initAnalysis( const gmx::TrajectoryAnalysisSettings   &settings,
 526                         const gmx::TopologyInformation          &top)
 527 {
 528     // считывание индекса
 529     index.resize(0);
 530     for (gmx::ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ++ai) {
 531         index.push_back(static_cast< size_t >(*ai));
 532     }
 533     // считывание индекса
 534     indexCA.resize(0);
 535     for (gmx::ArrayRef< const int >::iterator ai {sel_[1].atomIndices().begin()}; (ai < sel_[1].atomIndices().end()); ++ai) {
 536         indexCA.push_back(static_cast< size_t >(*ai));
 537     }
 538     // считывание красок
 539     colorsNames.resize(0);
 540     colorsNames.resize(sel_.size() - 2);
 541     for (size_t i {2}; i < sel_.size(); ++i) {
 542         for (gmx::ArrayRef< const int >::iterator ai {sel_[i].atomIndices().begin()}; (ai < sel_[i].atomIndices().end()); ++ai) {
 543             colorsNames[i - 2].push_back(std::make_pair(*(top.atoms()->atomname[*ai]), static_cast< size_t >(*ai)));
 544         }
 545     }
 546     colorsVectorsIndex.resize(0);
 547     colorsVectorsIndex.resize(colorsNames.size());
 548     colorAnglesPairsEvaluation(colorsNames, colorsVectorsIndex);
 549     // разбор dat файла для создания бэта листов
 550     betaListDigestion(fnBetaListsDat, betaLists);
 551     // разбор топологии для индексации бета листов
 552     aminoacidsIndexation(index, top, aminoacidsIndex);
 553
 554     // запоминание распределения зарядов в красках
 555     colorsMassCharges.resize(0);
 556     colorsMassCharges.resize(colorsNames.size());
 557     for (size_t i {0}; i < colorsNames.size(); ++i) {
 558         colorMassChargesEval(colorsNames[i], top, colorsMassCharges[i]);
 559     }
 560
 561     // вывод базовых значений
 562     std::cout << std::endl;
 563     std::cout << "effective radius = " << std::setw(30) << effRad                   << std::endl;
 564     std::cout << "DAT file         = " << std::setw(30) << fnBetaListsDat           << std::endl;
 565     std::cout << "output file      = " << std::setw(30) << fnOut                    << std::endl;
 566     std::cout << "peptide atoms #  = " << std::setw(30) << index.size()             << std::endl;
 567     std::cout << "CA-atmos #       = " << std::setw(30) << indexCA.size()           << std::endl;
 568     std::cout << "Beta frames      = " << std::setw(30) << betaLists.size()         << std::endl;
 569     std::cout << "residue #        = " << std::setw(30) << aminoacidsIndex.size()   << std::endl;
 570     std::cout << "colors #         = " << std::setw(30) << colorsNames.size()       << std::endl;
 571 }
 572
 573 void
 574 colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::TrajectoryAnalysisModuleData *pdata)
 575 {
 576     trajectoryFrame.resize(0);
 577     // считывания текущего фрейма траектории
 578     for (size_t i {0}; i < fr.natoms; ++i) {
 579         trajectoryFrame.push_back(fr.x[i]);
 580     }
 581     // подсчёт углов в красках
 582     colorsAngles.resize(0);
 583     colorsAnglesEvaluation(trajectoryFrame, colorsVectorsIndex, colorsAngles);
 584     // рассчёт положения относительно белка
 585     std::vector< bool > colorsToPeptide;
 586     colorsToPeptide.resize(0);
 587     colorsToPeptide.resize(colorsNames.size(), false);
 588     for (size_t i {0}; i < colorsNames.size(); ++i) {
 589         colorsToPeptide[i] = isNearPeptide(pbc, trajectoryFrame, index, colorsNames[i], effRad);
 590     }
 591     // расчёт угла и среднего угла с ближайшими бета листами
 592     std::vector< std::vector< bool > > colorsToBeta;
 593     colorsToBeta.resize(0);
 594     colorsToBeta.resize(colorsNames.size());
 595     std::vector< std::vector< double > > colorsToBetaAngles;
 596     colorsToBetaAngles.resize(0);
 597     std::vector< gmx::RVec > betaListsRVecs;
 598     // подсчёт направляющей beta листа
 599     betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA);
 600     // поиск ближайших beta листов
 601     for (size_t i {0}; i < colorsToPeptide.size(); ++i) {
 602         if (colorsToPeptide[i]) {
 603             searchNearBetaLists(pbc, trajectoryFrame, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex, effRad);
 604         }
 605     }
 606     // подсчёт углов между векторами краски и beta листами
 607     for (size_t i {0}; i < colorsToBeta.size(); ++i) {
 608         for (size_t j {0}; j < colorsToBeta[i].size(); ++j) {
 609             if (colorsToBeta[i][j]) {
 610                 computeAnglesColorsVsBeta(betaListsRVecs[j], colorsAngles[i]);
 611             }
 612         }
 613     }
 614     // поиск стакированных красок
 615     std::vector< std::vector< int > > colorStack;
 616     colorsStackingSearch(pbc, trajectoryFrame, colorsNames, stackRad, stackCoef, colorStack);
 617     // рассчёт дипольных моментов для красок
 618     std::vector< gmx::RVec > colorsDipole;
 619     std::vector< std::vector< double > > dipoleBetaAngles;
 620     colorsDipole.resize(0);
 621     dipoleBetaAngles.resize(0);
 622     dipoleBetaAngles.resize(colorsNames.size());
 623     for (auto &i : colorsMassCharges) {
 624         colorsDipole.push_back(colorDipoleEval(i, trajectoryFrame));
 625     }
 626     for (size_t i {0}; i < colorsToBeta.size(); ++i) {
 627         dipoleBetaAngles[i].resize(0);
 628         for (size_t j {0}; j < colorsToBeta[i].size(); ++j) {
 629             if (colorsToBeta[i][j]) {
 630                 dipoleBetaAngles[i].push_back(RVecAngle(colorsDipole[i], betaListsRVecs[j]));
 631             }
 632         }
 633     }
 634     // вывод в файл(ы) информацию для фрейма
 635     anglesFileDump(frnr, fnOut, colorsToPeptide, colorsAngles, colorStack, dipoleBetaAngles);
 636     std::cout << "frame number = " << std::setw(8) << frnr << " trajectory time = " << std::setw(8) << fr.time << std::endl;
 637 }
 638
 639 void
 640 colorVec::finishAnalysis(int nframes)
 641 {
 642
 643 }
 644
 645 void
 646 colorVec::writeOutput()
 647 {
 648     /*
 649      *
 650      * frame #<current frame number>
 651      * 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>]
 652      *
 653      */
 654 }
 655
 656 /*! \brief
 657  * The main function for the analysis template.
 658  */
 659 int
 660 main(int argc, char *argv[])
 661 {
 662     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<colorVec>(argc, argv);
 663 }