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