2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Implements gmx::analysismodules::Freevolume.
39 * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
40 * \ingroup module_trajectoryanalysis
43 #include <gromacs/trajectoryanalysis.h>
44 #include <gromacs/trajectoryanalysis/topologyinformation.h>
45 #include <gromacs/selection/nbsearch.h>
56 // структура углов для одной краски
57 struct colorLocalAngles {
58 gmx::RVec a1{0., 0., 0.}, a2{0., 0., 0.};
59 gmx::RVec b1{0., 0., 0.}, b2{0., 0., 0.};
60 gmx::RVec n1{0., 0., 0.}, n2{0., 0., 0.};
61 double a12 {0}, b12 {0}, n12 {0};
62 std::vector< double > betaAngles;
65 // хрен пойми почему, но захотелось рекурсивно сделать | можно сделать и через вайл
66 void returnBetaEnd(const std::string ¤tLine, long long &pos) {
67 switch (currentLine[pos]) {
69 returnBetaEnd(currentLine, ++pos);
72 returnBetaEnd(currentLine, ++pos);
81 // функция парсинга одной строки
82 void parseBetaListDATLine(const std::string ¤tLine, std::vector< std::vector< unsigned int > > &localInputBL) {
83 size_t equalCount = 0;
84 std::vector< unsigned int > a;
86 for (size_t i = 0; i < currentLine.size(); ++i) {
87 if (currentLine[i] == '=') { // подсчитываем число "пустых" символов
91 returnBetaEnd(currentLine, temp);
92 std::cout << temp << " ";
93 if (temp - static_cast< long long >(i) > 3) {
94 localInputBL.push_back(a);
95 for (size_t j = i; j <= temp; ++j) {
96 localInputBL.back().push_back(j - equalCount);
104 // функция нахождения бета-листов в структуре по файлу ДССП
105 void betaListDigestion(const std::string &inputFile, std::vector< std::vector< std::vector< unsigned int > > > &inputBL) {
107 std::ifstream file(inputFile);
109 getline(file, line); // считываем число в первой строке - кол-во осмысленных элементов в строках - нам не нужно
111 if (line.size() > 3) {
113 inputBL.resize(inputBL.size() + 1);
114 parseBetaListDATLine(line, inputBL.back());
117 } while (line.size() > 3);
119 throw "DSSP DAT FILE IS EMPTY";
124 // функция выделения индексов для бэталистов
125 inline void aminoacidsIndexation(const std::vector< size_t > &inputIndex, const gmx::TopologyInformation &top, std::vector< std::vector< size_t > > &aminoacidsIndex) {
126 aminoacidsIndex.resize(0);
127 for (size_t i = 0; i < inputIndex.size(); ++i) {
128 aminoacidsIndex.resize(std::max(aminoacidsIndex.size(), static_cast< size_t >(top.atoms()->atom[inputIndex[i]].resind + 1)));
129 aminoacidsIndex[top.atoms()->atom[inputIndex[i]].resind].push_back(inputIndex[i]);
133 // функция поиска RVec в кадре по имени->индексу
134 gmx::RVec returnRVec(const std::vector< gmx::RVec > &frame, const std::vector< std::pair< std::string, size_t > > &colorsIndex, const std::string &toFind) {
135 for (auto &i : colorsIndex) {
136 if (i.first == toFind) {
137 return frame[i.second];
140 throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS";
143 // функция векторного произведения двух RVec'ов
144 inline void RVecVecMultiply(gmx::RVec &n, const gmx::RVec &a, const gmx::RVec &b) {
145 n[0] = a[1] * b[2] - a[2] * b[1];
146 n[1] = a[2] * b[0] - a[0] * b[2];
147 n[2] = a[0] * b[1] - a[2] * b[0];
150 // поиск угла между двумя RVec'ами
151 inline double RVecAngle(const gmx::RVec &a, const gmx::RVec &b) {
152 return std::acos((a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) / (a.norm() * b.norm())) * 180.0 / 3.14159265;
155 // вычисление внутренних углов в краске
156 void colorsAnglesEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex,
157 std::vector< colorLocalAngles > &colorStruct) {
158 colorStruct.resize(colorsIndex.size());
159 #pragma omp parallel for ordered schedule(dynamic)
160 for (size_t i = 0; i < colorsIndex.size(); ++i) {
161 colorStruct[i].betaAngles.resize(0);
163 colorStruct[i].a1 = returnRVec(frame, colorsIndex[i], "CAF") - returnRVec(frame, colorsIndex[i], "CAJ") + returnRVec(frame, colorsIndex[i], "CAK") - returnRVec(frame, colorsIndex[i], "CAO");
164 colorStruct[i].a1 /= colorStruct[i].a1.norm();
165 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");
166 colorStruct[i].b1 /= colorStruct[i].b1.norm();
167 RVecVecMultiply(colorStruct[i].n1, colorStruct[i].a1, colorStruct[i].b1);
168 colorStruct[i].n1 /= colorStruct[i].n1.norm();
170 colorStruct[i].a2 = returnRVec(frame, colorsIndex[i], "CAB") - returnRVec(frame, colorsIndex[i], "CBQ") + returnRVec(frame, colorsIndex[i], "CBP") - returnRVec(frame, colorsIndex[i], "CBL");
171 colorStruct[i].a2 /= colorStruct[i].a2.norm();
172 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");
173 colorStruct[i].b2 /= colorStruct[i].b2.norm();
174 RVecVecMultiply(colorStruct[i].n2, colorStruct[i].a2, colorStruct[i].b2);
175 colorStruct[i].n2 *= -1; // для "сонаправленности" векторов нормали
176 colorStruct[i].n2 /= colorStruct[i].n2.norm();
177 colorStruct[i].a12 = RVecAngle(colorStruct[i].a1, colorStruct[i].a2);
178 colorStruct[i].b12 = RVecAngle(colorStruct[i].b1, colorStruct[i].b2);
179 colorStruct[i].n12 = RVecAngle(colorStruct[i].n1, colorStruct[i].n2);
184 // вычисление направляющего вектора в бэта-листе
185 inline void betaListsRVecsEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< std::vector< unsigned int > > &inputBetaLists,
186 std::vector< gmx::RVec > &temp, const std::vector< size_t > &inputCA) {
189 for (const auto &i : inputBetaLists) {
190 tempA = frame[inputCA[i[i.size() - 1]]] + frame[inputCA[i[i.size() - 2]]] - frame[inputCA[i[1]]] - frame[inputCA[i[0]]];
191 tempA /= tempA.norm();
192 temp.push_back(tempA);
196 //// определение близко ли к белку находится краска
197 //bool isNearPeptide(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > &inputFrame,
198 // /*gmx::AnalysisNeighborhood &nbhood, */const std::vector< size_t > &inputIndex,
199 // const std::vector< std::pair< std::string, size_t > > &inputColor) {
200 // /*gmx::AnalysisNeighborhoodSearch nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
201 // gmx::AnalysisNeighborhoodPair pair;*/
203 // for (size_t i = 0; i < inputColor.size(); ++i) {
204 // std::cout << inputColor[i].first << " " << i << " / " << inputColor.size() << std::endl;
205 // gmx::AnalysisNeighborhood nbhood;
206 // nbhood.setCutoff(0.8);
207 // gmx::AnalysisNeighborhoodSearch nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
208 // gmx::AnalysisNeighborhoodPair pair;
209 // gmx::AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(inputFrame[inputColor[i].second].as_vec());
211 // while (pairSearch.findNextPair(&pair)) {
212 // std::cout << ++count1 << " ";
213 // for (size_t j = 0; j < inputIndex.size(); ++j) {
214 // if (pair.refIndex() == inputIndex[j]) {
219 // std::cout << inputColor[i].first << " atom done" << std::endl;
224 // определение близко ли к белку находится краска
225 bool isNearPeptide( const std::vector< gmx::RVec > &inputFrame, const std::vector< size_t > &inputIndex,
226 const std::vector< std::pair< std::string, size_t > > &inputColor, const double cutOff) {
227 for (size_t i {0}; i < inputColor.size(); ++i) {
228 for (size_t j {0}; j < inputIndex.size(); ++j) {
229 if ((inputFrame[inputIndex[j]] - inputFrame[inputColor[i].second]).norm() <= cutOff) {
237 // поиск ближайших к краске бэта-листов
238 //inline void searchNearBetaLists(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > &inputFrame,
239 // gmx::AnalysisNeighborhood &nbhood,
240 // const std::vector< std::vector< unsigned int > > &inputBLists,
241 // const std::vector< std::pair< std::string, size_t > > &inputColor,
242 // std::vector< bool > &outputList, const std::vector< std::vector< size_t > > &inputAminoacids) {
243 // outputList.resize(0);
244 // outputList.resize(inputBLists.size(), false);
245 // gmx::AnalysisNeighborhoodSearch nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
246 // gmx::AnalysisNeighborhoodPair pair;
247 // for (const auto &i : inputColor) {
248 // while (nbsearch.startPairSearch(inputFrame[i.second].as_vec()).findNextPair(&pair)) {
249 // for (size_t j = 0; j < inputBLists.size(); ++j) {
250 // for (size_t k = 0; k < inputBLists[j].size(); ++k) {
251 // for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
252 // if (pair.testIndex() == inputAminoacids[inputBLists[j][k]][m]) {
253 // outputList[j] = true;
262 inline void searchNearBetaLists(const std::vector< gmx::RVec > &inputFrame, const std::vector< std::vector< unsigned int > > &inputBLists,
263 const std::vector< std::pair< std::string, size_t > > &inputColor, std::vector< bool > &outputList,
264 const std::vector< std::vector< size_t > > &inputAminoacids, const double cutOff) {
265 outputList.resize(0);
266 outputList.resize(inputBLists.size(), false);
267 //std::cout << inputColor.size() << " " << inputBLists.size() << " " << inputAminoacids.size() << "\n";
268 for (size_t i {0}; i < inputColor.size(); ++i) {
269 for (size_t j = 0; j < inputBLists.size(); ++j) {
270 for (size_t k = 0; k < inputBLists[j].size(); ++k) {
271 for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) {
272 std::cout << inputFrame[inputAminoacids[inputBLists[j][k]][m]][0] << " " << inputFrame[inputAminoacids[inputBLists[j][k]][m]][1] << " " <<
273 inputFrame[inputAminoacids[inputBLists[j][k]][m]][2] << "\n";
274 if ((inputFrame[inputAminoacids[inputBLists[j][k]][m]] - inputFrame[inputColor[i].second]).norm() <= cutOff) {
275 outputList[j] = true;
277 m = inputAminoacids[inputBLists[j][k]].size() + 1;
278 k = inputBLists[j].size() + 1;
286 // определение углов между краской и направляющим вектором бэта-листа
287 inline void computeAnglesColorsVsBeta(const gmx::RVec &inputBetaRVec, colorLocalAngles &colorFormation) {
288 colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a1));
289 colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b1));
290 colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n1));
291 colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.a2));
292 colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.b2));
293 colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n2));
296 // функция записи в файл значений углов для кадра
297 void anglesFileDump(const int frameNum, const std::string &output, const std::vector< bool > &toPeptide, const std::vector< colorLocalAngles > &colorFormation) {
298 std::ofstream file(output, std::ofstream::app);
300 std::vector< double > betaTemp;
302 file << "frame =" << std::setw(8) << frameNum << std::endl;
303 for (size_t i = 0; i < colorFormation.size(); ++i) {
304 file << "color #" << std::setw(3) << i;
306 file << std::setw(4) << "yes";
308 file << std::setw(4) << "no";
310 file << std::setw(7) << std::setprecision(2) << colorFormation[i].a12 << std::setw(7) << colorFormation[i].b12 << std::setw(7) << colorFormation[i].n12;
311 temp = colorFormation[i].betaAngles.size() / 6;
313 file << std::setw(4) << 0;
315 betaTemp.resize(6, 0.); // magic number, meh
316 for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
317 betaTemp[j % 6] += colorFormation[i].betaAngles[j];
319 for (size_t j = 0; j < betaTemp.size(); ++j) {
322 file << std::setw(4) << temp;
323 for (size_t j = 0; j < betaTemp.size(); ++j) {
324 file << std::setw(7) << std::setprecision(2) << betaTemp[j];
326 for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) {
327 file << std::setw(7) << std::setprecision(2) << colorFormation[i].betaAngles[j];
336 * \ingroup module_trajectoryanalysis
338 class colorVec : public gmx::TrajectoryAnalysisModule
345 //! Set the options and setting
346 virtual void initOptions(gmx::IOptionsContainer *options,
347 gmx::TrajectoryAnalysisSettings *settings);
349 //! First routine called by the analysis framework
350 // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
351 virtual void initAnalysis(const gmx::TrajectoryAnalysisSettings &settings,
352 const gmx::TopologyInformation &top);
354 //! Call for each frame of the trajectory
355 // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
356 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
357 gmx::TrajectoryAnalysisModuleData *pdata);
359 //! Last routine called by the analysis framework
360 // virtual void finishAnalysis(t_pbc *pbc);
361 virtual void finishAnalysis(int nframes);
363 //! Routine to write output, that is additional over the built-in
364 virtual void writeOutput();
368 gmx::SelectionList sel_;
369 std::string fnOut {"name"}; // selectable
370 std::string fnBetaListsDat; // selectable
371 float effRad {0.8}; // selectable
372 std::vector< size_t > index;
373 std::vector< size_t > indexCA;
374 std::vector< std::vector< size_t > > aminoacidsIndex;
375 std::vector< std::vector< std::pair< std::string, size_t > > > colorsNames;
376 std::vector< std::vector< std::vector< unsigned int > > > betaLists;
377 gmx::AnalysisNeighborhood nb_;
379 // Copy and assign disallowed by base.
382 colorVec::colorVec(): TrajectoryAnalysisModule()
386 colorVec::~colorVec()
399 colorVec::initOptions( gmx::IOptionsContainer *options,
400 gmx::TrajectoryAnalysisSettings *settings)
402 static const char *const desc[] = {
403 "[THISMODULE] to be done"
405 // Add the descriptive text (program help text) to the options
406 settings->setHelpText(desc);
407 // Add option for selection list
408 options->addOption(gmx::SelectionOption("sel")
410 .required().dynamicMask().multiValue()
411 .description("select pepride and colors / -sf"));
412 // Add option for input file names
413 options->addOption(gmx::StringOption("dat")
414 .store(&fnBetaListsDat)
415 .description("a file to make dynamic beta lists"));
416 // Add option for output file name
417 options->addOption(gmx::StringOption("out")
419 .description("Index file for the algorithm output."));
420 // Add option for effRad constant
421 options->addOption(gmx::FloatOption("efRad")
423 .description("max distance from colors to peptide in nm to consider to be \"near\""));
424 // Control input settings
425 settings->setFlags(gmx::TrajectoryAnalysisSettings::efNoUserPBC);
426 settings->setFlag(gmx::TrajectoryAnalysisSettings::efUseTopX);
427 //settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
428 settings->setPBC(true);
432 colorVec::initAnalysis( const gmx::TrajectoryAnalysisSettings &settings,
433 const gmx::TopologyInformation &top)
435 // считывание индекса
437 for (gmx::ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ai++) {
438 index.push_back(static_cast< size_t >(*ai));
440 // считывание индекса
442 for (gmx::ArrayRef< const int >::iterator ai {sel_[1].atomIndices().begin()}; (ai < sel_[2].atomIndices().end()); ai++) {
443 index.push_back(static_cast< size_t >(*ai));
446 colorsNames.resize(sel_.size() - 2);
447 for (size_t i = 2; i < sel_.size(); ++i) {
448 for (gmx::ArrayRef< const int >::iterator ai {sel_[i].atomIndices().begin()}; (ai < sel_[i].atomIndices().end()); ai++) {
449 colorsNames[i - 2].push_back(std::make_pair(*(top.atoms()->atomname[*ai]), static_cast< size_t >(*ai)));
452 // разбор dat файла для создания бэта листов
453 betaListDigestion(fnBetaListsDat, betaLists);
454 for (size_t i = 0; i < betaLists.size(); ++i) {
455 std::cout << betaLists[i].size() << " ";
458 // разбор топологии для индексации бета листов
459 aminoacidsIndexation(index, top, aminoacidsIndex);
460 // задание радиуса рассматриваемых соседей
461 nb_.setCutoff(effRad);
463 * формально можно сделать:
464 * найти соотношение между именами для углов и индексными номерами
465 * заполнить std::vector< std::vector< size_t > > nameToIndex;
466 * и в последствии считать:
467 * a1 = frame[nameToIndex[0][1] + ... - ... - ...
468 * b1 a2 b2 по аналогии только для [@<-[1, 2, 3]]
474 colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::TrajectoryAnalysisModuleData *pdata)
476 std::vector< gmx::RVec > trajectoryFrame;
477 trajectoryFrame.resize(0);
478 // считывания текущего фрейма траектории
479 for (size_t i {0}; i < fr.natoms; ++i) {
480 trajectoryFrame.push_back(fr.x[i]);
482 // подсчёт углов в красках
483 std::vector< colorLocalAngles > colorFormation;
484 colorsAnglesEvaluation(trajectoryFrame, colorsNames, colorFormation);
485 // рассчёт положения относительно белка
487 * формально можно совместить определение близости с белком и определение ближайших бэта-листов
488 * для этого думаю нужно как-то соотнести или сделать соотношение между индексом и бэта-листами (т.е. хранить что бы за О(1) делать)
492 std::vector< bool > colorsToPeptide;
493 colorsToPeptide.resize(0);
494 colorsToPeptide.resize(colorsNames.size(), false);
495 for (size_t i = 0; i < colorsNames.size(); ++i) {
496 colorsToPeptide[i] = isNearPeptide(trajectoryFrame, index, colorsNames[i], effRad * 9000);
498 // расчёт угла и среднего угла с ближайшими бета листами
499 std::vector< std::vector< bool > > colorsToBeta;
500 colorsToBeta.resize(0);
501 colorsToBeta.resize(colorsNames.size());
502 std::vector< std::vector< double > > colorsToBetaAngles;
503 colorsToBetaAngles.resize(0);
504 std::vector< gmx::RVec > betaListsRVecs;
505 betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA);
506 for (size_t i = 0; i < colorsToPeptide.size(); ++i) {
507 if (colorsToPeptide[i]) {
508 searchNearBetaLists(trajectoryFrame, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex, effRad * 9000);
511 for (size_t i = 0; i < colorsToBeta.size(); ++i) {
512 for (size_t j = 0; j < colorsToBeta[i].size(); ++j) {
513 if (colorsToBeta[i][j]) {
514 computeAnglesColorsVsBeta(betaListsRVecs[j], colorFormation[i]);
518 // вывод в файл(ы) информацию для фрейма
519 anglesFileDump(frnr, fnOut, colorsToPeptide, colorFormation);
523 colorVec::finishAnalysis(int nframes)
529 colorVec::writeOutput()
533 * frame #<current frame number>
534 * 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>]
540 * The main function for the analysis template.
543 main(int argc, char *argv[])
545 return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<colorVec>(argc, argv);