- 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     gmx::AnalysisNeighborhoodSearch     nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
199     gmx::AnalysisNeighborhoodPair       pair;
200     for (const auto &i : inputColor) {
201         gmx::AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(inputFrame[i.second].as_vec());
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         std::cout << " " << i.first;
213     }
214     return false;
215 }
216
217 // поиск ближайших к краске бэта-листов
218 inline void searchNearBetaLists(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > &inputFrame,
219                                 gmx::AnalysisNeighborhood &nbhood,
220                                 const std::vector< std::vector< unsigned int > > &inputBLists,
221                                 const std::vector< std::pair< std::string, size_t > > &inputColor,
222                                 std::vector< bool > &outputList, const std::vector< std::vector< size_t > > &inputAminoacids) {
223     outputList.resize(0);
224     outputList.resize(inputBLists.size(), false);
225     gmx::AnalysisNeighborhoodSearch      nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
226     gmx::AnalysisNeighborhoodPair        pair;
227     for (const auto &i : inputColor) {
228         while (nbsearch.startPairSearch(inputFrame[i.second].as_vec()).findNextPair(&pair)) {
229             for (size_t j = 0; j < inputBLists.size(); ++j) {
230                 for (size_t k = 0; k < inputBLists[j].size(); ++k) {
231                     for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
232                         if (pair.testIndex() == inputAminoacids[inputBLists[j][k]][m]) {
233                             outputList[j] = true;
234                         }
235                     }
236                 }
237             }
238         }
239     }
240
241 }
242
243 // определение углов между краской и направляющим вектором бэта-листа
244 inline void computeAnglesColorsVsBeta(const gmx::RVec &inputBetaRVec, colorLocalAngles &colorFormation) {
245     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a1));
246     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b1));
247     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n1));
248     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a2));
249     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b2));
250     colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n2));
251 }
252
253 // функция записи в файл значений углов для кадра
254 void anglesFileDump(const int frameNum, const std::string &output, const std::vector< bool > &toPeptide, const std::vector< colorLocalAngles > &colorFormation) {
255     std::ofstream   file(output);
256     int temp = 0;
257     std::vector< double > betaTemp;
258     betaTemp.resize(0);
259     file << "frame =" << std::setw(8) << frameNum << std::endl;
260     for (size_t i = 0; i < colorFormation.size(); ++i) {
261         file << "color #" << std::setw(3) << i;
262         if (toPeptide[i]) {
263             file << std::setw(4) << "yes";
264         } else {
265             file << std::setw(4) << "no";
266         }
267         file << std::setw(7) << std::setprecision(2) << colorFormation[i].a12 << colorFormation[i].b12 << colorFormation[i].n12;
268         temp = colorFormation[i].betaAngles.size() / 6;
269         if (temp == 0) {
270             file << std::setw(4) << 0;
271         } else {
272             betaTemp.resize(6, 0.); // magic number, meh
273             for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
274                 betaTemp[j % 6] += colorFormation[i].betaAngles[j];
275             }
276             for (size_t j = 0; j < betaTemp.size(); ++j) {
277                 betaTemp[j] /= temp;
278             }
279             file << std::setw(4) << temp;
280             for (size_t j = 0; j < betaTemp.size(); ++j) {
281                 file << std::setw(7) << std::setprecision(2) << betaTemp[j];
282             }
283             for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
284                 file << std::setw(7) << std::setprecision(2) << colorFormation[i].betaAngles[j];
285             }
286         }
287     }
288     file.close();
289 }
290
291 /*! \brief
292  * \ingroup module_trajectoryanalysis
293  */
294 class colorVec : public gmx::TrajectoryAnalysisModule
295 {
296     public:
297
298         colorVec();
299         virtual ~colorVec();
300
301         //! Set the options and setting
302         virtual void initOptions(gmx::IOptionsContainer          *options,
303                                  gmx::TrajectoryAnalysisSettings *settings);
304
305         //! First routine called by the analysis framework
306         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
307         virtual void initAnalysis(const gmx::TrajectoryAnalysisSettings &settings,
308                                   const gmx::TopologyInformation        &top);
309
310         //! Call for each frame of the trajectory
311         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
312         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
313                                   gmx::TrajectoryAnalysisModuleData *pdata);
314
315         //! Last routine called by the analysis framework
316         // virtual void finishAnalysis(t_pbc *pbc);
317         virtual void finishAnalysis(int nframes);
318
319         //! Routine to write output, that is additional over the built-in
320         virtual void writeOutput();
321
322     private:
323
324         gmx::SelectionList                                              sel_;
325         std::string                                                     fnOut           {"name"};   // selectable
326         std::string                                                     fnBetaListsDat;             // selectable
327         float                                                           effRad          {0.8};      // selectable
328         std::vector< size_t >                                           index;
329         std::vector< size_t >                                           indexCA;
330         std::vector< std::vector< size_t > >                            aminoacidsIndex;
331         std::vector< std::vector< std::pair< std::string, size_t > > >  colorsNames;
332         std::vector< std::vector< std::vector< unsigned int > > >       betaLists;
333         gmx::AnalysisNeighborhood                                       nb_;
334
335         // Copy and assign disallowed by base.
336 };
337
338 colorVec::colorVec(): TrajectoryAnalysisModule()
339 {
340 }
341
342 colorVec::~colorVec()
343 {
344 }
345
346 /*
347  * ndx:
348  * [peptide]
349  * [CA]
350  * [color_{i}]
351  *
352  */
353
354 void
355 colorVec::initOptions(  gmx::IOptionsContainer          *options,
356                         gmx::TrajectoryAnalysisSettings *settings)
357 {
358     static const char *const desc[] = {
359         "[THISMODULE] to be done"
360     };
361     // Add the descriptive text (program help text) to the options
362     settings->setHelpText(desc);
363     // Add option for selection list
364     options->addOption(gmx::SelectionOption("sel")
365                             .storeVector(&sel_)
366                             .required().dynamicMask().multiValue()
367                             .description("select pepride and colors / -sf"));
368     // Add option for input file names
369     options->addOption(gmx::StringOption("dat")
370                             .store(&fnBetaListsDat)
371                             .description("a file to make dynamic beta lists"));
372     // Add option for output file name
373     options->addOption(gmx::StringOption("out")
374                             .store(&fnOut)
375                             .description("Index file for the algorithm output."));
376     // Add option for effRad constant
377     options->addOption(gmx::FloatOption("efRad")
378                             .store(&effRad)
379                             .description("max distance from colors to peptide in nm to consider to be \"near\""));
380     // Control input settings
381     settings->setFlags(gmx::TrajectoryAnalysisSettings::efNoUserPBC);
382     settings->setFlag(gmx::TrajectoryAnalysisSettings::efUseTopX);
383     //settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
384     settings->setPBC(true);
385 }
386
387 void
388 colorVec::initAnalysis( const gmx::TrajectoryAnalysisSettings   &settings,
389                         const gmx::TopologyInformation          &top)
390 {
391     // считывание индекса
392     index.resize(0);
393     for (gmx::ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ai++) {
394         index.push_back(static_cast< size_t >(*ai));
395     }
396     // считывание индекса
397     indexCA.resize(0);
398     for (gmx::ArrayRef< const int >::iterator ai {sel_[1].atomIndices().begin()}; (ai < sel_[2].atomIndices().end()); ai++) {
399         index.push_back(static_cast< size_t >(*ai));
400     }
401     // считывание красок
402     colorsNames.resize(sel_.size() - 2);
403     for (size_t i = 2; i < sel_.size(); ++i) {
404         for (gmx::ArrayRef< const int >::iterator ai {sel_[i].atomIndices().begin()}; (ai < sel_[i].atomIndices().end()); ai++) {
405             colorsNames[i - 2].push_back(std::make_pair(*(top.atoms()->atomname[*ai]), static_cast< size_t >(*ai)));
406         }
407     }
408     // разбор dat файла для создания бэта листов
409     betaListDigestion(fnBetaListsDat, betaLists);
410     // разбор топологии для индексации бета листов
411     aminoacidsIndexation(index, top, aminoacidsIndex);
412     // задание радиуса рассматриваемых соседей
413     nb_.setCutoff(effRad);
414     /*
415      * формально можно сделать:
416      * найти соотношение между именами для углов и индексными номерами
417      * заполнить std::vector< std::vector< size_t > > nameToIndex;
418      * и в последствии считать:
419      * a1 = frame[nameToIndex[0][1] + ... - ... - ...
420      * b1 a2 b2 по аналогии только для [@<-[1, 2, 3]]
421      *
422      */
423 }
424
425 void
426 colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::TrajectoryAnalysisModuleData *pdata)
427 {
428     std::cout << "\tFrame analisis start ..." << std::endl;
429     std::vector< gmx::RVec > trajectoryFrame;
430     trajectoryFrame.resize(0);
431     // считывания текущего фрейма траектории
432     for (size_t i {0}; i < index.size(); ++i) {
433         trajectoryFrame.push_back(fr.x[index[i]]);
434     }
435     // подсчёт углов в красках
436     std::vector< colorLocalAngles > colorFormation;
437     std::cout << "\t\tColors' angles evaluation." << std::endl;
438     colorsAnglesEvaluation(trajectoryFrame, colorsNames, colorFormation);
439     // рассчёт положения относительно белка
440     /*
441      * формально можно совместить определение близости с белком и определение ближайших бэта-листов
442      * для этого думаю нужно как-то соотнести или сделать соотношение между индексом и бэта-листами (т.е. хранить что бы за О(1) делать)
443      *
444      */
445
446     std::vector< bool >     colorsToPeptide;
447     colorsToPeptide.resize(0);
448     colorsToPeptide.resize(colorsNames.size(), false);
449     std::cout << "\t\tWhich colors are \"near\" the pepride mass." << std::endl;
450     for (size_t i = 0; i < colorsNames.size(); ++i) {
451         colorsToPeptide[i] = isNearPeptide(fr, pbc, trajectoryFrame, nb_, index, colorsNames[i]);
452     }
453     // расчёт угла и среднего угла с ближайшими бета листами
454     std::vector< std::vector< bool > > colorsToBeta;
455     colorsToBeta.resize(0);
456     colorsToBeta.resize(colorsNames.size());
457     std::vector< std::vector< double > > colorsToBetaAngles;
458     colorsToBetaAngles.resize(0);
459     std::vector< gmx::RVec > betaListsRVecs;
460     std::cout << "\t\tBeta-lists' RVecs search." << std::endl;
461     betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA);
462     std::cout << "\t\tSearching nearby beta-lists." << std::endl;
463     for (size_t i = 0; i < colorsToPeptide.size(); ++i) {
464         if (colorsToPeptide[i]) {
465             searchNearBetaLists(fr, pbc, trajectoryFrame, nb_, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex);
466         }
467     }
468     std::cout << "\t\tComputing angles colors vs betas." << std::endl;
469     for (size_t i = 0; i < colorsToBeta.size(); ++i) {
470         for (size_t j = 0; j < colorsToBeta[i].size(); ++j) {
471             if (colorsToBeta[i][j]) {
472                 computeAnglesColorsVsBeta(betaListsRVecs[j], colorFormation[i]);
473             }
474         }
475     }
476     // вывод в файл(ы) информацию для фрейма
477     std::cout << "\t\tDumping data." << std::endl;
478     anglesFileDump(frnr, fnOut, colorsToPeptide, colorFormation);
479     std::cout << "\tFrame analysed." << std::endl;
480 }
481
482 void
483 colorVec::finishAnalysis(int nframes)
484 {
485
486 }
487
488 void
489 colorVec::writeOutput()
490 {
491     /*
492      *
493      * frame #<current frame number>
494      * 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>]
495      *
496      */
497     std::cout << "\n\t colorAngles finished successfully" << std::endl;
498 }
499
500 /*! \brief
501  * The main function for the analysis template.
502  */
503 int
504 main(int argc, char *argv[])
505 {
506     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<colorVec>(argc, argv);
507 }