2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2021, 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::Trajectory.
39 * \author Sergey Gorelov <gorelov_sv@pnpi.nrcki.ru>
40 * \author Anatoly Titov <titov_ai@pnpi.nrcki.ru>
41 * \author Alexey Shvetsov <alexxyum@gmail.com>
42 * \ingroup module_trajectoryanalysis
46 There's something wrong with energy of the last residue
50 #include "dssptools.h"
53 #include "gromacs/math/units.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include <gromacs/trajectoryanalysis.h>
57 #include "gromacs/trajectoryanalysis/topologyinformation.h"
66 namespace analysismodules
69 //void ResInfo::setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex)
71 // _ResInfo.at(static_cast<std::size_t>(atomTypeName)) = atomIndex;
74 //std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const
76 // return _ResInfo[static_cast<std::size_t>(atomTypeName)];
79 std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const{
80 return _backboneIndices[static_cast<std::size_t>(atomTypeName)];
83 secondaryStructures::secondaryStructures(){
85 void secondaryStructures::initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez){
86 SecondaryStructuresStatusMap.resize(0);
87 SecondaryStructuresStringLine.resize(0);
88 std::vector<std::size_t> temp; temp.resize(0),
89 PiHelixPreference = PiHelicesPreferencez;
90 ResInfoMap = &ResInfoMatrix;
91 SecondaryStructuresStatusMap.resize(ResInfoMatrix.size());
92 SecondaryStructuresStringLine.resize(ResInfoMatrix.size(), '~');
95 void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName){
96 SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)] = true;
99 void secondaryStructures::secondaryStructuresData::setStatus(const HelixPositions helixPosition, const turnsTypes turn){
100 TurnsStatusArray[static_cast<std::size_t>(turn)] = helixPosition;
103 bool secondaryStructures::secondaryStructuresData::getStatus(const secondaryStructureTypes secondaryStructureTypeName) const{
104 return SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)];
107 bool secondaryStructures::secondaryStructuresData::isBreakPartnerWith(const secondaryStructuresData *partner) const{
108 return breakPartners[0] == partner || breakPartners[1] == partner;
111 HelixPositions secondaryStructures::secondaryStructuresData::getStatus(const turnsTypes turn) const{
112 return TurnsStatusArray[static_cast<std::size_t>(turn)];
115 void secondaryStructures::secondaryStructuresData::setBreak(secondaryStructuresData *breakPartner){
116 if (breakPartners[0] == nullptr){
117 breakPartners[0] = breakPartner;
120 breakPartners[1] = breakPartner;
122 setStatus(secondaryStructureTypes::Break);
125 bool secondaryStructures::hasHBondBetween(std::size_t Donor, std::size_t Acceptor) const{ // prob should add resi name comparison ?
126 if( (*ResInfoMap)[Donor].acceptor[0] == nullptr ||
127 (*ResInfoMap)[Donor].acceptor[1] == nullptr ||
128 (*ResInfoMap)[Acceptor].info == nullptr ){
133 // std::cout << "Comparing DONOR " << Donor << " And ACCEPTOR " << Acceptor << " :";
134 // std::cout << "DONOR's acceptor adresses are " << (*ResInfoMap)[Donor].acceptor[0] << ", " << (*ResInfoMap)[Donor].acceptor[1] << " and ACCEPTOR adress is " << (*ResInfoMap)[Acceptor].info << std::endl;
135 // std::cout << "DONOR's acceptors' nr are = " << (*ResInfoMap)[Donor].acceptor[0]->nr << ", " << (*ResInfoMap)[Donor].acceptor[1]->nr << " And ACCEPTOR's nr = " << (*ResInfoMap)[Acceptor].info->nr << std::endl;
136 // std::cout << "DONOR's acceptors' chainID are = " << (*ResInfoMap)[Donor].acceptor[0]->chainid << ", " << (*ResInfoMap)[Donor].acceptor[1]->chainid << " And ACCEPTOR's chainID = " << (*ResInfoMap)[Acceptor].info->chainid << std::endl;
138 // if( ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
139 // ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff ) ){
140 // std::cout << "HBond Exist" << std::endl;
142 return ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
143 ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff );
147 bool secondaryStructures::NoChainBreaksBetween(std::size_t Resi1, std::size_t Resi2) const{
148 std::size_t i{Resi1}, j{Resi2}; // From i to j → i <= j
154 if ( !(SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i])) ){
157 std::cout << "Patternsearch has detected a CHAINBREAK between " << Resi1 << " and " << Resi2 << std::endl;
163 bridgeTypes secondaryStructures::calculateBridge(std::size_t i, std::size_t j) const{
164 if( i < 1 || j < 1 || i + 1 >= ResInfoMap->size() || j + 1 >= ResInfoMap->size() ){
165 return bridgeTypes::None;
167 if(NoChainBreaksBetween(i - 1, i + 1) && NoChainBreaksBetween(j - 1, j + 1)){
168 if((hasHBondBetween(i + 1, j) && hasHBondBetween(j, i - 1)) || (hasHBondBetween(j + 1, i) && hasHBondBetween(i, j - 1)) ){ //possibly swap
169 return bridgeTypes::ParallelBridge;
171 else if((hasHBondBetween(i + 1, j - 1) && hasHBondBetween(j + 1, i - 1)) || (hasHBondBetween(j, i) && hasHBondBetween(i, j)) ){ //possibly swap
172 return bridgeTypes::AntiParallelBridge;
175 return bridgeTypes::None;
178 void secondaryStructures::analyzeBridgesAndLaddersPatterns(){
179 for(std::size_t i {1}; i + 4 < SecondaryStructuresStatusMap.size(); ++i){
180 for(std::size_t j {i + 3}; j + 1 < SecondaryStructuresStatusMap.size(); ++j ){
181 bridgeTypes type {calculateBridge(i, j)};
182 if (type == bridgeTypes::None){
200 // for (std::size_t i{ 1 }; i < HBondsMap.front().size() - 1; ++i){
201 // for (std::size_t j{ 1 }; j < HBondsMap.front().size() - 1; ++j){
202 // if (std::abs(static_cast<int>(i) - static_cast<int>(j)) > 2){
203 // if ((HBondsMap[i - 1][j] && HBondsMap[j][i + 1]) ||
204 // (HBondsMap[j - 1][i] && HBondsMap[i][j + 1])){
205 // Bridge[i].push_back(j);
207 // if ((HBondsMap[i][j] && HBondsMap[j][i]) ||
208 // (HBondsMap[i - 1][j + 1] && HBondsMap[j - 1][i + 1])){
209 // AntiBridge[i].push_back(j);
214 // for (std::size_t i{ 0 }; i < HBondsMap.front().size(); ++i){
215 // if ((!Bridge[i].empty() || !AntiBridge[i].empty())){
216 // setStatus(i, secondaryStructureTypes::Bulge);
219 // for (std::size_t i{ 2 }; i + 2 < HBondsMap.front().size(); ++i){
220 // for (std::size_t j { i - 2 }; j <= (i + 2); ++j){
225 // for (std::vector<bridgeTypes>::const_iterator bridge {Bridges.begin()}; bridge != Bridges.end(); ++bridge ){
226 // if (!getBridge(*bridge)[i].empty() || !getBridge(*bridge)[j].empty()){
227 // for (std::size_t i_resi{ 0 }; i_resi < getBridge(*bridge)[i].size(); ++i_resi){
228 // for (std::size_t j_resi{ 0 }; j_resi < getBridge(*bridge)[j].size(); ++j_resi){
229 // if (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
230 // - static_cast<int>(getBridge(*bridge)[j][j_resi]))
231 // && (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
232 // - static_cast<int>(getBridge(*bridge)[j][j_resi]))
235 // for (std::size_t k{ 0 }; k <= i - j; ++k){
236 // setStatus(i + k, secondaryStructureTypes::Ladder);
240 // for (std::size_t k{ 0 }; k <= j - i; ++k){
241 // setStatus(i + k, secondaryStructureTypes::Ladder);
254 void secondaryStructures::analyzeTurnsAndHelicesPatterns(){
255 for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
256 std::cout << "Testing Helix_" << static_cast<std::size_t>(i) + 3 << std::endl;
257 std::size_t stride {static_cast<std::size_t>(i) + 3};
258 for(std::size_t j {0}; j + static_cast<std::size_t>(i) < SecondaryStructuresStatusMap.size(); ++j){
259 std::cout << "Testing " << j << " and " << j + stride << std::endl;
260 if ( hasHBondBetween(j, j + (static_cast<std::size_t>(i))) && NoChainBreaksBetween(j, j + stride) ){
261 std::cout << j << " and " << j + stride << " has hbond!" << std::endl;
262 SecondaryStructuresStatusMap[j + static_cast<std::size_t>(i)].setStatus(HelixPositions::End, i);
264 for (std::size_t k {1}; k < (static_cast<std::size_t>(i)); ++k){
265 if( SecondaryStructuresStatusMap[j + k].getStatus(i) == HelixPositions::None ){
266 SecondaryStructuresStatusMap[j + k].setStatus(HelixPositions::Middle, i);
271 if( SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::End ){
272 SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start_AND_End, i);
275 SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start, i);
281 for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
282 std::size_t stride {static_cast<std::size_t>(i) + 3};
283 for(std::size_t j {1}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
284 if ( (SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start_AND_End ) &&
285 (SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start_AND_End ) ){
287 secondaryStructureTypes Helix;
289 case turnsTypes::Turn_3:
290 for (std::size_t k {0}; empty && k < stride; ++k){
291 empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_3);
293 Helix = secondaryStructureTypes::Helix_3;
295 case turnsTypes::Turn_5:
296 for (std::size_t k {0}; empty && k < stride; ++k){
297 empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_5) || (PiHelixPreference && SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_5));
299 Helix = secondaryStructureTypes::Helix_4;
302 Helix = secondaryStructureTypes::Helix_4;
305 std::cout << j << " is HELIX" << std::endl;
307 for(std::size_t k {0}; k < (static_cast<std::size_t>(i)); ++k ){
308 SecondaryStructuresStatusMap[j + k].setStatus(Helix);
315 for(std::size_t i {1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
316 if (SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Loop)){
318 for(const turnsTypes &j : {turnsTypes::Turn_3, turnsTypes::Turn_4, turnsTypes::Turn_5}){
319 std::size_t stride {static_cast<std::size_t>(i) + 3};
320 for(std::size_t k {1}; k < stride; ++k){
321 isTurn = (i >= k) && (SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start || SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start_AND_End) ;
326 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Turn);
328 else if (SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Bend)){
329 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bend);
336 void secondaryStructures::analyzePPHelicesPatterns(){}
338 std::string secondaryStructures::patternSearch(){
341 // analyzeBridgesAndLaddersPatterns();
342 analyzeTurnsAndHelicesPatterns();
343 // analyzePPHelicesPatterns();
344 NoChainBreaksBetween(80, 75);
346 // for(std::size_t i {0}; i < ResInfoMap->size(); ++i){
347 // std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) << std::endl;
350 // std::cout.precision(5);
351 // for(std::size_t i{0}; i < ResInfoMap->size(); ++i, std::cout << std::endl << std::endl){
352 // std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) ;
353 // if ( (*ResInfoMap)[i].donor[0] != nullptr ){
354 // std::cout << " has donor[0] = " << (*ResInfoMap)[i].donor[0]->nr << " " << *((*ResInfoMap)[i].donor[0]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[0] << " and" ;
357 // std::cout << " has no donor[0] and" ;
359 // if ( (*ResInfoMap)[i].acceptor[0] != nullptr ){
360 // std::cout << " has acceptor[0] = " << (*ResInfoMap)[i].acceptor[0]->nr << " " << *((*ResInfoMap)[i].acceptor[0]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[0] ;
363 // std::cout << " has no acceptor[0]" ;
365 // std::cout << std::endl << "Also, " << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name);
366 // if ( (*ResInfoMap)[i].donor[1] != nullptr ){
367 // std::cout << " has donor[1] = " << (*ResInfoMap)[i].donor[1]->nr << " " << *((*ResInfoMap)[i].donor[1]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[1] << " and" ;
370 // std::cout << " has no donor[1] and" ;
372 // if ( (*ResInfoMap)[i].acceptor[1] != nullptr ){
373 // std::cout << " has acceptor[1] = " << (*ResInfoMap)[i].acceptor[1]->nr << " " << *((*ResInfoMap)[i].acceptor[1]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[1] ;
376 // std::cout << " has no acceptor[1]" ;
382 for(std::size_t i {static_cast<std::size_t>(secondaryStructureTypes::Bend)}; i != static_cast<std::size_t>(secondaryStructureTypes::Count); ++i){
383 for(std::size_t j {0}; j < SecondaryStructuresStatusMap.size(); ++j){
384 if (SecondaryStructuresStatusMap[j].getStatus(static_cast<secondaryStructureTypes>(i))){
385 SecondaryStructuresStringLine[j] = secondaryStructureTypeNames[i] ;
392 if(SecondaryStructuresStatusMap.size() > 1){
393 for(std::size_t i {0}, linefactor{1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
394 if( SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) && SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break) ){
395 if(SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]) ){
396 SecondaryStructuresStringLine.insert(SecondaryStructuresStringLine.begin() + i + linefactor, secondaryStructureTypeNames[secondaryStructureTypes::Break]);
402 return SecondaryStructuresStringLine;
405 secondaryStructures::~secondaryStructures(){
406 SecondaryStructuresStatusMap.resize(0);
407 SecondaryStructuresStringLine.resize(0);
410 DsspTool::DsspStorage::DsspStorage(){
411 storaged_data.resize(0);
414 void DsspTool::DsspStorage::clearAll(){
415 storaged_data.resize(0);
418 std::mutex DsspTool::DsspStorage::mx;
420 void DsspTool::DsspStorage::storageData(int frnr, std::string data){
421 std::lock_guard<std::mutex> guardian(mx);
422 std::pair<int, std::string> datapair(frnr, data);
423 storaged_data.push_back(datapair);
426 std::vector<std::pair<int, std::string>> DsspTool::DsspStorage::returnData(){
427 std::sort(storaged_data.begin(), storaged_data.end());
428 return storaged_data;
431 void alternateNeighborhoodSearch::setCutoff(const real &cutoff_init){
432 cutoff = cutoff_init;
435 void alternateNeighborhoodSearch::FixAtomCoordinates(real &coordinate, const real vector_length){
436 while (coordinate < 0) {
437 coordinate += vector_length;
439 while (coordinate >= vector_length) {
440 coordinate -= vector_length;
444 void alternateNeighborhoodSearch::ReCalculatePBC(int &x, const int &x_max) {
453 void alternateNeighborhoodSearch::GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
454 rvec coordinates, box_vector_length;
455 num_of_miniboxes.resize(0);
456 num_of_miniboxes.resize(3);
457 for (std::size_t i{XX}; i <= ZZ; ++i) {
458 box_vector_length[i] = std::sqrt(
459 std::pow(fr.box[i][XX], 2) + std::pow(fr.box[i][YY], 2) + std::pow(fr.box[i][ZZ], 2));
460 num_of_miniboxes[i] = std::floor((box_vector_length[i] / cutoff)) + 1;
462 MiniBoxesMap.resize(0);
463 MiniBoxesReverseMap.resize(0);
464 MiniBoxesMap.resize(num_of_miniboxes[XX], std::vector<std::vector<std::vector<std::size_t> > >(
465 num_of_miniboxes[YY], std::vector<std::vector<std::size_t> >(
466 num_of_miniboxes[ZZ], std::vector<std::size_t>(
468 MiniBoxesReverseMap.resize(IndexMap.size(), std::vector<std::size_t>(3));
469 for (std::vector<ResInfo>::const_iterator i {IndexMap.begin()}; i != IndexMap.end(); ++i) {
470 for (std::size_t j{XX}; j <= ZZ; ++j) {
471 coordinates[j] = fr.x[i->getIndex(backboneAtomTypes::AtomCA)][j];
472 FixAtomCoordinates(coordinates[j], box_vector_length[j]);
474 MiniBoxesMap[std::floor(coordinates[XX] / cutoff)][std::floor(coordinates[YY] / cutoff)][std::floor(
475 coordinates[ZZ] / cutoff)].push_back(i - IndexMap.begin());
476 for (std::size_t j{XX}; j <= ZZ; ++j){
477 MiniBoxesReverseMap[i - IndexMap.begin()][j] = std::floor(coordinates[j] / cutoff);
482 void alternateNeighborhoodSearch::AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
483 GetMiniBoxesMap(fr, IndexMap);
484 MiniBoxSize[XX] = MiniBoxesMap.size();
485 MiniBoxSize[YY] = MiniBoxesMap.front().size();
486 MiniBoxSize[ZZ] = MiniBoxesMap.front().front().size();
488 PairMap.resize(IndexMap.size(), std::vector<bool>(IndexMap.size(), false));
489 ResiI = PairMap.begin();
490 ResiJ = ResiI->begin();
492 for (std::vector<ResInfo>::const_iterator i = IndexMap.begin(); i != IndexMap.end(); ++i){
493 for (offset[XX] = -1; offset[XX] <= 1; ++offset[XX]) {
494 for (offset[YY] = -1; offset[YY] <= 1; ++offset[YY]) {
495 for (offset[ZZ] = -1; offset[ZZ] <= 1; ++offset[ZZ]) {
496 for (std::size_t k{XX}; k <= ZZ; ++k) {
497 fixBox[k] = MiniBoxesReverseMap[i - IndexMap.begin()][k] + offset[k];
498 ReCalculatePBC(fixBox[k], MiniBoxSize[k]);
500 for (std::size_t j{0}; j < MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]].size(); ++j) {
501 if ( (i - IndexMap.begin()) != MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]){
502 PairMap[i - IndexMap.begin()][MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]] = true;
503 PairMap[MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]][i - IndexMap.begin()] = true;
512 bool alternateNeighborhoodSearch::findNextPair(){
518 for(; ResiI != PairMap.end(); ++ResiI, ResiJ = ResiI->begin() ){
519 for(; ResiJ != ResiI->end(); ++ResiJ){
521 resiIpos = ResiI - PairMap.begin();
522 resiJpos = ResiJ - ResiI->begin();
523 if ( ResiJ != ResiI->end() ){
526 else if (ResiI != PairMap.end()) {
528 ResiJ = ResiI->begin();
541 std::size_t alternateNeighborhoodSearch::getResiI() const {
545 std::size_t alternateNeighborhoodSearch::getResiJ() const {
550 DsspTool::DsspStorage DsspTool::Storage;
552 DsspTool::DsspTool(){
555 void DsspTool::calculateBends(const t_trxframe &fr, const t_pbc *pbc)
557 const float benddegree{ 70.0 }, maxdist{ 2.5 };
558 float degree{ 0 }, vdist{ 0 }, vprod{ 0 };
559 gmx::RVec a{ 0, 0, 0 }, b{ 0, 0, 0 };
560 for (std::size_t i{ 0 }; i + 1 < IndexMap.size(); ++i)
562 if (CalculateAtomicDistances(static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
563 static_cast<int>(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)),
568 PatternSearch.SecondaryStructuresStatusMap[i].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i + 1]);
569 PatternSearch.SecondaryStructuresStatusMap[i + 1].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i]);
571 // std::cout << "Break between " << i + 1 << " and " << i + 2 << std::endl;
574 for (std::size_t i{ 2 }; i + 2 < IndexMap.size() ; ++i)
576 if (PatternSearch.SecondaryStructuresStatusMap[i - 2].getStatus(secondaryStructureTypes::Break) ||
577 PatternSearch.SecondaryStructuresStatusMap[i - 1].getStatus(secondaryStructureTypes::Break) ||
578 PatternSearch.SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) ||
579 PatternSearch.SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break)
584 for (int j{ 0 }; j < 3; ++j)
586 a[j] = fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j]
587 - fr.x[IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA)][j];
588 b[j] = fr.x[IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA)][j]
589 - fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j];
591 vdist = (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
592 vprod = CalculateAtomicDistances(IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA),
593 IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
596 * gmx::c_angstrom / gmx::c_nano
597 * CalculateAtomicDistances(IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
598 IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA),
601 * gmx::c_angstrom / gmx::c_nano;
602 degree = std::acos(vdist / vprod) * gmx::c_rad2Deg;
603 if (degree > benddegree)
605 PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bend);
610 void DsspTool::calculateHBondEnergy(ResInfo& Donor,
612 const t_trxframe& fr,
616 * DSSP uses eq from dssp 2.x
617 * kCouplingConstant = 27.888, // = 332 * 0.42 * 0.2
618 * E = k * (1/rON + 1/rCH - 1/rOH - 1/rCN) where CO comes from one AA and NH from another
623 const float kCouplingConstant = 27.888;
624 const float minimalAtomDistance{ 0.5 },
626 float HbondEnergy{ 0 };
627 float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 };
629 if( !(Donor.is_proline) ){
630 if (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
631 && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || (initParams.addHydrogens) ) ) // Kinda ew
634 distanceNO = CalculateAtomicDistances(
635 Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
636 distanceNC = CalculateAtomicDistances(
637 Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
639 if (initParams.addHydrogens){
640 if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
641 // std::cout << "On donor " << Donor.info->nr << *(Donor.info->name) << std::endl;
642 // std::cout << "Prev donor is " << Donor.prevResi->info->nr << *(Donor.prevResi->info->name) << std::endl;
643 // std::cout << "Prev C index is " << Donor.prevResi->getIndex(backboneAtomTypes::AtomC) << std::endl;
644 // std::cout << "Prev O index is " << Donor.prevResi->getIndex(backboneAtomTypes::AtomO) << std::endl;
646 float prevCODist {CalculateAtomicDistances(Donor.prevResi->getIndex(backboneAtomTypes::AtomC), Donor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)};
647 for (int i{XX}; i <= ZZ; ++i){
648 float prevCO = fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
649 atomH[i] = prevCO / prevCODist;
651 distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
652 distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
655 distanceHO = distanceNO;
656 distanceHC = distanceNC;
660 distanceHO = CalculateAtomicDistances(
661 Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
662 distanceHC = CalculateAtomicDistances(
663 Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
666 if (CalculateAtomicDistances(
667 Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
670 if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
671 || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
673 HbondEnergy = minEnergy;
675 // std::cout << "HBOND exists cause of distance" << std::endl;
681 * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
682 HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
684 // std::cout.precision(5);
685 // std::cout << "Calculated ENERGY = " << HbondEnergy << std::endl;
687 // if ( HbondEnergy == 0){
688 // std::cout << "Calculated ENERGY = " << HbondEnergy << " For donor " << Donor.info->nr << " and acceptor " << Acceptor.info->nr << std::endl;
691 if ( HbondEnergy < minEnergy ){
692 HbondEnergy = minEnergy;
699 // std::cerr << "PRO DETECTED! THIS IS RESI № " << Donor.info->nr << std::endl; //IT WORKS JUST FINE
702 if (HbondEnergy < Donor.acceptorEnergy[0]){
703 Donor.acceptor[1] = Donor.acceptor[0];
704 Donor.acceptor[0] = Acceptor.info;
705 Donor.acceptorEnergy[0] = HbondEnergy;
707 else if (HbondEnergy < Donor.acceptorEnergy[1]){
708 Donor.acceptor[1] = Acceptor.info;
709 Donor.acceptorEnergy[1] = HbondEnergy;
712 if (HbondEnergy < Acceptor.donorEnergy[0]){
713 Acceptor.donor[1] = Acceptor.donor[0];
714 Acceptor.donor[0] = Donor.info;
715 Acceptor.donorEnergy[0] = HbondEnergy;
717 else if (HbondEnergy < Acceptor.donorEnergy[1]){
718 Acceptor.donor[1] = Donor.info;
719 Acceptor.donorEnergy[1] = HbondEnergy;
724 /* Calculate Distance From B to A */
725 float DsspTool::CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
727 gmx::RVec r{ 0, 0, 0 };
728 pbc_dx(pbc, fr.x[A], fr.x[B], r.as_vec());
729 // return r.norm(); // TODO * gmx::c_nm2A; if not PDB, i guess????
730 return r.norm() * gmx::c_nm2A; // TODO * gmx::c_nm2A; if not PDB, i guess????
733 /* Calculate Distance From B to A, where A is only fake coordinates */
734 float DsspTool::CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
736 gmx::RVec r{ 0, 0, 0 };
737 pbc_dx(pbc, A, fr.x[B], r.as_vec());
738 // return r.norm(); // TODO * gmx::c_nm2A; if not PDB, i guess????
739 return r.norm() * gmx::c_nm2A; // TODO * gmx::c_nm2A; if not PDB, i guess????
742 void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz)
745 std::cout << "Init started" << std::endl;
746 initParams = initParamz;
747 ResInfo _backboneAtoms;
750 int resicompare{ top.atoms()->atom[static_cast<std::size_t>(*(initParams.sel_.atomIndices().begin()))].resind };
752 IndexMap.push_back(_backboneAtoms);
753 IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
754 proLINE = *(IndexMap[i].info->name);
755 if( proLINE.compare("PRO") == 0 ){
756 IndexMap[i].is_proline = true;
759 for (gmx::ArrayRef<const int>::iterator ai{ initParams.sel_.atomIndices().begin() }; (ai != initParams.sel_.atomIndices().end()); ++ai){
760 if (resicompare != top.atoms()->atom[static_cast<std::size_t>(*ai)].resind)
763 resicompare = top.atoms()->atom[static_cast<std::size_t>(*ai)].resind;
764 IndexMap.emplace_back(_backboneAtoms);
765 IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
766 proLINE = *(IndexMap[i].info->name);
767 if( proLINE.compare("PRO") == 0 ){
768 IndexMap[i].is_proline = true;
772 std::string atomname(*(top.atoms()->atomname[static_cast<std::size_t>(*ai)]));
773 if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA])
775 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomCA)] = *ai;
777 else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC])
779 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomC)] = *ai;
781 else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO])
783 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomO)] = *ai;
785 else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN])
787 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomN)] = *ai;
789 else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH] && initParamz.addHydrogens == false) // Юзать водород в структуре
791 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
794 // if( atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO]
795 // || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH]){
796 // std::cout << "Atom " << atomname << " №" << *ai << " From Resi " << *(top.atoms()->resinfo[i].name) << " №" << resicompare << std::endl;
800 for (std::size_t j {1}; j < IndexMap.size(); ++j){
801 IndexMap[j].prevResi = &(IndexMap[j - 1]);
803 IndexMap[j - 1].nextResi = &(IndexMap[j]);
805 // std::cout << "Resi " << IndexMap[i].info->nr << *(IndexMap[i].info->name) << std::endl;
806 // std::cout << "Prev resi is " << IndexMap[i].prevResi->info->nr << *(IndexMap[i].prevResi->info->name) << std::endl;
807 // std::cout << "Prev resi's next resi is " << IndexMap[i - 1].nextResi->info->nr << *(IndexMap[i - 1].nextResi->info->name) << std::endl;
808 // std::cout << IndexMap[j].prevResi->info->nr;
809 // std::cout << *(IndexMap[j].prevResi->info->name) ;
810 // std::cout << " have CA = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomCA) ;
811 // std::cout << " C = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomC);
812 // std::cout << " O = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomO);
813 // std::cout << " N = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomN);
814 // std::cout << " H = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomH) << std::endl;
820 void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc)
823 switch(initParams.NBS){
824 case (NBSearchMethod::Classique): {
826 // store positions of CA atoms to use them for nbSearch
827 std::vector<gmx::RVec> positionsCA_;
828 for (std::size_t i{ 0 }; i < IndexMap.size(); ++i)
830 positionsCA_.emplace_back(fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)]);
833 AnalysisNeighborhood nb_;
834 nb_.setCutoff(initParams.cutoff_);
835 AnalysisNeighborhoodPositions nbPos_(positionsCA_);
836 gmx::AnalysisNeighborhoodSearch start = nb_.initSearch(pbc, nbPos_);
837 gmx::AnalysisNeighborhoodPairSearch pairSearch = start.startPairSearch(nbPos_);
838 gmx::AnalysisNeighborhoodPair pair;
839 while (pairSearch.findNextPair(&pair))
841 if(CalculateAtomicDistances(
842 IndexMap[pair.refIndex()].getIndex(backboneAtomTypes::AtomCA), IndexMap[pair.testIndex()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
843 < minimalCAdistance){
844 calculateHBondEnergy(IndexMap[pair.refIndex()], IndexMap[pair.testIndex()], fr, pbc);
845 if (IndexMap[pair.testIndex()].info != IndexMap[pair.refIndex() + 1].info){
846 calculateHBondEnergy(IndexMap[pair.testIndex()], IndexMap[pair.refIndex()], fr, pbc);
853 case (NBSearchMethod::Experimental): { // TODO FIX
855 alternateNeighborhoodSearch as_;
857 as_.setCutoff(initParams.cutoff_);
859 as_.AltPairSearch(fr, IndexMap);
861 while (as_.findNextPair()){
862 if(CalculateAtomicDistances(
863 IndexMap[as_.getResiI()].getIndex(backboneAtomTypes::AtomCA), IndexMap[as_.getResiJ()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
864 < minimalCAdistance){
865 calculateHBondEnergy(IndexMap[as_.getResiI()], IndexMap[as_.getResiJ()], fr, pbc);
866 if (IndexMap[as_.getResiJ()].info != IndexMap[as_.getResiI() + 1].info){
867 calculateHBondEnergy(IndexMap[as_.getResiJ()], IndexMap[as_.getResiI()], fr, pbc);
876 for(std::vector<ResInfo>::iterator Donor {IndexMap.begin()}; Donor != IndexMap.end() ; ++Donor){
877 for(std::vector<ResInfo>::iterator Acceptor {Donor + 1} ; Acceptor != IndexMap.end() ; ++Acceptor){
878 if(CalculateAtomicDistances(
879 Donor->getIndex(backboneAtomTypes::AtomCA), Acceptor->getIndex(backboneAtomTypes::AtomCA), fr, pbc)
880 < minimalCAdistance){
881 calculateHBondEnergy(*Donor, *Acceptor, fr, pbc);
882 if (Acceptor != Donor + 1){
883 calculateHBondEnergy(*Acceptor, *Donor, fr, pbc);
893 // for(std::size_t i {0}; i < IndexMap.size(); ++i){
894 // std::cout << IndexMap[i].info->nr << " " << *(IndexMap[i].info->name) << std::endl;
897 PatternSearch.initiateSearch(IndexMap, initParams.PPHelices);
898 calculateBends(fr, pbc);
899 Storage.storageData(frnr, PatternSearch.patternSearch());
903 std::vector<std::pair<int, std::string>> DsspTool::getData(){
904 return Storage.returnData();
907 } // namespace analysismodules