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>
46 #include <gromacs/pbcutil/pbc.h>
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;
66 struct colorAnglesPairs
68 std::vector< std::pair<size_t, size_t> > a1, a2;
69 std::vector< std::pair<size_t, size_t> > b1, b2;
72 // функция парсинга одной строки
73 void parseBetaListDATLine(const std::string ¤tLine, std::vector< std::vector< unsigned int > > &localInputBL) {
74 size_t equalCount {0};
75 std::vector< unsigned int > a;
77 for (size_t i {0}; i < currentLine.size(); ++i) {
78 if (currentLine[i] == '=') { // подсчитываем число "пустых" символов
82 for (size_t j {i + 1}; j < currentLine.size(); ++j) {
83 if (currentLine[j] == 'B' || currentLine[j] == 'E') {
90 localInputBL.push_back(a);
91 for (size_t j {i}; j <= temp; ++j) {
92 localInputBL.back().push_back(j - equalCount);
100 // функция нахождения бета-листов в структуре по файлу ДССП
101 void betaListDigestion(const std::string &inputFile, std::vector< std::vector< std::vector< unsigned int > > > &inputBL) {
103 std::ifstream file(inputFile);
105 getline(file, line); // считываем число в первой строке - кол-во осмысленных элементов в строках - нам не нужно
107 if (line.size() > 3) {
109 inputBL.resize(inputBL.size() + 1);
110 parseBetaListDATLine(line, inputBL.back());
113 } while (line.size() > 3);
115 throw "DSSP DAT FILE IS EMPTY";
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]);
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];
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;
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);
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];
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];
155 RVecVecMultiply(colorStruct[i].n1, colorStruct[i].a1, colorStruct[i].b1);
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];
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];
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);
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) {
176 if (j.first == toFind) {
181 throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS";
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")));
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) {
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);
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;*/
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());
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]) {
240 // std::cout << inputColor[i].first << " atom done" << std::endl;
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) {
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) {
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;
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);
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;
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));
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);
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) {
346 output.resize(colorsIndex.size());
347 for (size_t i {0}; i < colorsIndex.size(); ++i) {
350 for (size_t i {0}; i < outputTemp.size(); ++i) {
351 for (size_t j {i + 1}; 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);
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;
368 file << std::setw(4) << "yes";
370 file << std::setw(4) << "no";
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;
379 int temp = colorFormation[i].betaAngles.size() / 6;
380 std::vector< double > betaTemp;
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];
386 for (size_t j {0}; j < betaTemp.size(); ++j) {
387 betaTemp[j] /= static_cast< double >(temp);
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];
393 for (size_t j {0}; j < colorFormation[i].betaAngles.size(); ++j) {
394 file << std::setw(8) << std::setprecision(2) << colorFormation[i].betaAngles[j];
397 if (inputStack[i].size() == 0) {
398 file << std::setw(4) << "no";
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];
412 * \ingroup module_trajectoryanalysis
414 class colorVec : public gmx::TrajectoryAnalysisModule
421 //! Set the options and setting
422 virtual void initOptions(gmx::IOptionsContainer *options,
423 gmx::TrajectoryAnalysisSettings *settings);
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);
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);
435 //! Last routine called by the analysis framework
436 // virtual void finishAnalysis(t_pbc *pbc);
437 virtual void finishAnalysis(int nframes);
439 //! Routine to write output, that is additional over the built-in
440 virtual void writeOutput();
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_;
458 // Copy and assign disallowed by base.
461 colorVec::colorVec(): TrajectoryAnalysisModule()
465 colorVec::~colorVec()
478 colorVec::initOptions( gmx::IOptionsContainer *options,
479 gmx::TrajectoryAnalysisSettings *settings)
481 static const char *const desc[] = {
482 "[THISMODULE] to be done"
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")
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")
498 .description("Index file for the algorithm output."));
499 // Add option for effRad constant
500 options->addOption(gmx::FloatOption("efRad")
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")
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")
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);
519 colorVec::initAnalysis( const gmx::TrajectoryAnalysisSettings &settings,
520 const gmx::TopologyInformation &top)
522 // считывание индекса
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));
527 // считывание индекса
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));
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)));
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;
559 * формально можно сделать:
560 * найти соотношение между именами для углов и индексными номерами
561 * заполнить std::vector< std::vector< size_t > > nameToIndex;
562 * и в последствии считать:
563 * a1 = frame[nameToIndex[0][1] + ... - ... - ...
564 * b1 a2 b2 по аналогии только для [@<-[1, 2, 3]]
570 colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::TrajectoryAnalysisModuleData *pdata)
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]);
578 // подсчёт углов в красках
579 std::vector< colorLocalAngles > colorFormation;
580 colorsAnglesEvaluation(trajectoryFrame, colorsVectorsIndex, colorFormation);
581 // рассчёт положения относительно белка
583 * формально можно совместить определение близости с белком и определение ближайших бэта-листов
584 * для этого думаю нужно как-то соотнести или сделать соотношение между индексом и бэта-листами (т.е. хранить что бы за О(1) делать)
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);
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);
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]);
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;
624 colorVec::finishAnalysis(int nframes)
630 colorVec::writeOutput()
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>]
641 * The main function for the analysis template.
644 main(int argc, char *argv[])
646 return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<colorVec>(argc, argv);