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