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 <Infinity2573@gmail.com>
40 * \author Anatoly Titov <titov_ai@pnpi.nrcki.ru>
41 * \author Alexey Shvetsov <alexxyum@gmail.com>
42 * \ingroup module_trajectoryanalysis
45 #include "dssptools.h"
48 #include "gromacs/math/units.h"
50 #include "gromacs/pbcutil/pbc.h"
51 #include <gromacs/trajectoryanalysis.h>
52 #include "gromacs/trajectoryanalysis/topologyinformation.h"
61 namespace analysismodules
64 //void ResInfo::setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex)
66 // _ResInfo.at(static_cast<std::size_t>(atomTypeName)) = atomIndex;
69 //std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const
71 // return _ResInfo[static_cast<std::size_t>(atomTypeName)];
74 std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const{
75 return _backboneIndices[static_cast<std::size_t>(atomTypeName)];
78 secondaryStructures::secondaryStructures(){
80 void secondaryStructures::initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez, const int _pp_stretch){
81 SecondaryStructuresStatusMap.resize(0);
82 SecondaryStructuresStringLine.resize(0);
83 std::vector<std::size_t> temp; temp.resize(0),
84 PiHelixPreference = PiHelicesPreferencez;
85 ResInfoMap = &ResInfoMatrix;
86 pp_stretch = _pp_stretch;
87 SecondaryStructuresStatusMap.resize(ResInfoMatrix.size());
88 SecondaryStructuresStringLine.resize(ResInfoMatrix.size(), '~');
91 void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName, bool status){
92 SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)] = status;
94 SecondaryStructuresStatus = secondaryStructureTypeName;
97 SecondaryStructuresStatus = secondaryStructureTypes::Loop;
101 void secondaryStructures::secondaryStructuresData::setStatus(const HelixPositions helixPosition, const turnsTypes turn){
102 TurnsStatusArray[static_cast<std::size_t>(turn)] = helixPosition;
105 bool secondaryStructures::secondaryStructuresData::getStatus(const secondaryStructureTypes secondaryStructureTypeName) const{
106 return SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)];
109 bool secondaryStructures::secondaryStructuresData::isBreakPartnerWith(const secondaryStructuresData *partner) const{
110 return breakPartners[0] == partner || breakPartners[1] == partner;
113 HelixPositions secondaryStructures::secondaryStructuresData::getStatus(const turnsTypes turn) const{
114 return TurnsStatusArray[static_cast<std::size_t>(turn)];
117 secondaryStructureTypes secondaryStructures::secondaryStructuresData::getStatus() const{
118 return SecondaryStructuresStatus;
121 void secondaryStructures::secondaryStructuresData::setBreak(secondaryStructuresData *breakPartner){
122 if (breakPartners[0] == nullptr){
123 breakPartners[0] = breakPartner;
126 breakPartners[1] = breakPartner;
128 setStatus(secondaryStructureTypes::Break);
131 void secondaryStructures::secondaryStructuresData::setBridge(secondaryStructuresData *bridgePartner, std::size_t bridgePartnerIndex, bridgeTypes bridgeType){
132 if(bridgeType == bridgeTypes::ParallelBridge){
133 bridgePartners[0] = bridgePartner;
134 bridgePartnersIndexes[0] = bridgePartnerIndex;
136 else if (bridgeType == bridgeTypes::AntiParallelBridge){
137 bridgePartners[1] = bridgePartner;
138 bridgePartnersIndexes[1] = bridgePartnerIndex;
142 bool secondaryStructures::secondaryStructuresData::hasBridges() const{
143 return bridgePartners[0] || bridgePartners[1];
147 bool secondaryStructures::secondaryStructuresData::hasBridges(bridgeTypes bridgeType) const{
148 if(bridgeType == bridgeTypes::ParallelBridge){
149 return bridgePartners[0] != nullptr;
151 else if (bridgeType == bridgeTypes::AntiParallelBridge){
152 return bridgePartners[1] != nullptr;
159 bool secondaryStructures::secondaryStructuresData::isBridgePartnerWith(secondaryStructuresData *bridgePartner, bridgeTypes bridgeType) const{
160 if(bridgeType == bridgeTypes::ParallelBridge){
161 return bridgePartners[0] == bridgePartner;
163 else if (bridgeType == bridgeTypes::AntiParallelBridge){
164 return bridgePartners[1] == bridgePartner;
169 std::size_t secondaryStructures::secondaryStructuresData::getBridgePartnerIndex(bridgeTypes bridgeType) const{
170 if (bridgeType == bridgeTypes::ParallelBridge){
171 return bridgePartnersIndexes[0];
173 return bridgePartnersIndexes[1];
176 secondaryStructures::secondaryStructuresData secondaryStructures::secondaryStructuresData::getBridgePartner(bridgeTypes bridgeType) const{
177 if(bridgeType == bridgeTypes::ParallelBridge){
178 return *(bridgePartners[0]);
180 return *(bridgePartners[1]);
183 bool secondaryStructures::hasHBondBetween(std::size_t Donor, std::size_t Acceptor) const{
184 // if( (*ResInfoMap)[Donor].acceptor[0] == nullptr ||
185 // (*ResInfoMap)[Donor].acceptor[1] == nullptr ||
186 // (*ResInfoMap)[Acceptor].info == nullptr ){
187 // std::cout << "Bad hbond check. Reason(s): " ;
188 // if ( (*ResInfoMap)[Donor].acceptor[0] == nullptr ){
189 // std::cout << "Donor has no acceptor[0]; ";
191 // if ( (*ResInfoMap)[Donor].acceptor[1] == nullptr ){
192 // std::cout << "Donor has no acceptor[1]; ";
194 // if ( (*ResInfoMap)[Acceptor].info == nullptr ){
195 // std::cout << "No info about acceptor; ";
197 // std::cout << std::endl;
200 // else if (!( (*ResInfoMap)[Acceptor].donor[0] == nullptr ||
201 // (*ResInfoMap)[Acceptor].donor[1] == nullptr ||
202 // (*ResInfoMap)[Donor].info == nullptr )) {
203 // std::cout << "Comparing DONOR №" << (*ResInfoMap)[Donor].info->nr << " And ACCEPTOR №" << (*ResInfoMap)[Acceptor].info->nr << ": " << std::endl;
204 // std::cout << "Donor's acceptors' nr are = " << (*ResInfoMap)[Donor].acceptor[0]->nr << " (chain " << (*ResInfoMap)[Donor].acceptor[0]->chainid << ") , " << (*ResInfoMap)[Donor].acceptor[1]->nr << " (chain " << (*ResInfoMap)[Donor].acceptor[1]->chainid << ")" << std::endl;
205 // std::cout << "Donor's acceptors' energy are = " << (*ResInfoMap)[Donor].acceptorEnergy[0] << ", " << (*ResInfoMap)[Donor].acceptorEnergy[1] << std::endl;
206 // std::cout << "Acceptors's donors' nr are = " << (*ResInfoMap)[Acceptor].donor[0]->nr << " (chain " << (*ResInfoMap)[Acceptor].donor[0]->chainid << ") , " << (*ResInfoMap)[Acceptor].donor[1]->nr << " (chain " << (*ResInfoMap)[Acceptor].donor[1]->chainid << ")" << std::endl;
207 // std::cout << "Acceptors's donors' energy are = " << (*ResInfoMap)[Acceptor].donorEnergy[0] << ", " << (*ResInfoMap)[Acceptor].donorEnergy[1] << std::endl;
208 // if( ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
209 // ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff ) ){
210 // std::cout << "HBond Exist" << std::endl;
214 // std::cout << "Bad hbond check. Reason(s): " ;
215 // if ( (*ResInfoMap)[Acceptor].donor[0] == nullptr ){
216 // std::cout << "Acceptor has no donor[0]; ";
218 // if ( (*ResInfoMap)[Acceptor].donor[1] == nullptr ){
219 // std::cout << "Acceptor has no donor[1]; ";
221 // if ( (*ResInfoMap)[Donor].info == nullptr ){
222 // std::cout << "No info about donor; ";
224 // std::cout << std::endl;
227 return ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
228 ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff );
233 bool secondaryStructures::NoChainBreaksBetween(std::size_t Resi1, std::size_t Resi2) const{
234 std::size_t i{Resi1}, j{Resi2}; // From i to j → i <= j
240 if ( SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]) ){
241 // std::cout << "Patternsearch has detected a CHAINBREAK between " << Resi1 << " and " << Resi2 << std::endl;
248 bridgeTypes secondaryStructures::calculateBridge(std::size_t i, std::size_t j) const{
249 if( i < 1 || j < 1 || i + 1 >= ResInfoMap->size() || j + 1 >= ResInfoMap->size() ){ // Protection from idiotz, не обязательно нужен
250 return bridgeTypes::None;
253 ResInfo a{(*ResInfoMap)[i]}, b{(*ResInfoMap)[j]};
255 if(NoChainBreaksBetween(i - 1, i + 1) && NoChainBreaksBetween(j - 1, j + 1) && a.prevResi && a.nextResi && b.prevResi && b.nextResi){
256 if((hasHBondBetween(i + 1, j) && hasHBondBetween(j, i - 1)) || (hasHBondBetween(j + 1, i) && hasHBondBetween(i, j - 1)) ){
257 return bridgeTypes::ParallelBridge;
259 else if((hasHBondBetween(i + 1, j - 1) && hasHBondBetween(j + 1, i - 1)) || (hasHBondBetween(j, i) && hasHBondBetween(i, j)) ){
260 return bridgeTypes::AntiParallelBridge;
263 return bridgeTypes::None;
266 void secondaryStructures::analyzeBridgesAndStrandsPatterns(){
268 for(std::size_t i {1}; i + 4 < SecondaryStructuresStatusMap.size(); ++i){
269 for(std::size_t j {i + 3}; j + 1 < SecondaryStructuresStatusMap.size(); ++j ){
270 switch(calculateBridge(i, j)){
271 case bridgeTypes::ParallelBridge : {
272 SecondaryStructuresStatusMap[i].setBridge(&(SecondaryStructuresStatusMap[j]), j, bridgeTypes::ParallelBridge);
273 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bridge);
274 SecondaryStructuresStatusMap[j].setBridge(&(SecondaryStructuresStatusMap[i]), i, bridgeTypes::ParallelBridge);
275 SecondaryStructuresStatusMap[j].setStatus(secondaryStructureTypes::Bridge);
277 case bridgeTypes::AntiParallelBridge : {
278 SecondaryStructuresStatusMap[i].setBridge(&(SecondaryStructuresStatusMap[j]), j, bridgeTypes::AntiParallelBridge);
279 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bridge);
280 SecondaryStructuresStatusMap[j].setBridge(&(SecondaryStructuresStatusMap[i]), i, bridgeTypes::AntiParallelBridge);
281 SecondaryStructuresStatusMap[j].setStatus(secondaryStructureTypes::Bridge);
288 // Calculate Extended Strandz
289 for(std::size_t i {1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
290 bool is_Estrand {false};
291 for(std::size_t j {2}; j != 0; --j){
292 for(const bridgeTypes &bridgeType : {bridgeTypes::ParallelBridge, bridgeTypes::AntiParallelBridge}){
293 if (SecondaryStructuresStatusMap[i].hasBridges(bridgeType) && SecondaryStructuresStatusMap[i + j].hasBridges(bridgeType) ){
294 std::size_t i_partner{SecondaryStructuresStatusMap[i].getBridgePartnerIndex(bridgeType)}, j_partner{SecondaryStructuresStatusMap[i + j].getBridgePartnerIndex(bridgeType)}, second_strand{};
295 if ( abs(i_partner - j_partner) < 6){
296 if (i_partner < j_partner){
297 second_strand = i_partner;
300 second_strand = j_partner;
302 for(int k{abs(i_partner - j_partner)}; k >= 0; --k){
303 if (SecondaryStructuresStatusMap[second_strand + k].getStatus(secondaryStructureTypes::Bridge)){
304 SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Bridge, false);
307 SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Strand);
314 for(std::size_t k{0}; k <= j; ++k){
315 if (SecondaryStructuresStatusMap[i + k].getStatus(secondaryStructureTypes::Bridge)){
316 SecondaryStructuresStatusMap[i + k].setStatus(secondaryStructureTypes::Bridge, false);
318 SecondaryStructuresStatusMap[i + k].setStatus(secondaryStructureTypes::Strand);
336 // for (std::size_t i{ 1 }; i < HBondsMap.front().size() - 1; ++i){
337 // for (std::size_t j{ 1 }; j < HBondsMap.front().size() - 1; ++j){
338 // if (std::abs(static_cast<int>(i) - static_cast<int>(j)) > 2){
339 // if ((HBondsMap[i - 1][j] && HBondsMap[j][i + 1]) ||
340 // (HBondsMap[j - 1][i] && HBondsMap[i][j + 1])){
341 // Bridge[i].push_back(j);
343 // if ((HBondsMap[i][j] && HBondsMap[j][i]) ||
344 // (HBondsMap[i - 1][j + 1] && HBondsMap[j - 1][i + 1])){
345 // AntiBridge[i].push_back(j);
350 // for (std::size_t i{ 0 }; i < HBondsMap.front().size(); ++i){
351 // if ((!Bridge[i].empty() || !AntiBridge[i].empty())){
352 // setStatus(i, secondaryStructureTypes::Bulge);
355 // for (std::size_t i{ 2 }; i + 2 < HBondsMap.front().size(); ++i){
356 // for (std::size_t j { i - 2 }; j <= (i + 2); ++j){
361 // for (std::vector<bridgeTypes>::const_iterator bridge {Bridges.begin()}; bridge != Bridges.end(); ++bridge ){
362 // if (!getBridge(*bridge)[i].empty() || !getBridge(*bridge)[j].empty()){
363 // for (std::size_t i_resi{ 0 }; i_resi < getBridge(*bridge)[i].size(); ++i_resi){
364 // for (std::size_t j_resi{ 0 }; j_resi < getBridge(*bridge)[j].size(); ++j_resi){
365 // if (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
366 // - static_cast<int>(getBridge(*bridge)[j][j_resi]))
367 // && (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
368 // - static_cast<int>(getBridge(*bridge)[j][j_resi]))
371 // for (std::size_t k{ 0 }; k <= i - j; ++k){
372 // setStatus(i + k, secondaryStructureTypes::Ladder);
376 // for (std::size_t k{ 0 }; k <= j - i; ++k){
377 // setStatus(i + k, secondaryStructureTypes::Ladder);
390 void secondaryStructures::analyzeTurnsAndHelicesPatterns(){
391 for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
392 std::size_t stride {static_cast<std::size_t>(i) + 3};
393 for(std::size_t j {0}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
394 if (hasHBondBetween(j + stride, j)){
395 std::cout << "Bond between " << j << " and " << j + stride << " exists" << std::endl;
397 if (!NoChainBreaksBetween(j, j + stride)){
398 std::cout << "ChainBreak between " << j << " and " << j + stride << std::endl;
400 if ( hasHBondBetween(j + stride, j) && NoChainBreaksBetween(j, j + stride) ){
401 // std::cout << "Resi " << j << " is Helix_" << stride << " start" << std::endl;
402 SecondaryStructuresStatusMap[j + stride].setStatus(HelixPositions::End, i);
404 for (std::size_t k {1}; k < stride; ++k){
405 if( SecondaryStructuresStatusMap[j + k].getStatus(i) == HelixPositions::None ){
406 SecondaryStructuresStatusMap[j + k].setStatus(HelixPositions::Middle, i);
407 // SecondaryStructuresStatusMap[j + k].setStatus(secondaryStructureTypes::Turn);
411 if( SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::End ){
412 SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start_AND_End, i);
415 SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start, i);
421 for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
422 std::size_t stride {static_cast<std::size_t>(i) + 3};
423 for(std::size_t j {1}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
424 if ( (SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start_AND_End ) &&
425 (SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start_AND_End ) ){
427 secondaryStructureTypes Helix;
429 case turnsTypes::Turn_3:
430 for (std::size_t k {0}; empty && k < stride; ++k){
431 empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_3);
433 Helix = secondaryStructureTypes::Helix_3;
435 case turnsTypes::Turn_5:
436 for (std::size_t k {0}; empty && k < stride; ++k){
437 empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_5) || (PiHelixPreference && SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_4)); //TODO
439 Helix = secondaryStructureTypes::Helix_5;
442 Helix = secondaryStructureTypes::Helix_4;
445 if ( empty || Helix == secondaryStructureTypes::Helix_4 ){
446 for(std::size_t k {0}; k < stride; ++k ){
447 // std::cout << "Resi " << j << " is Helix_" << static_cast<std::size_t>(Helix) - 5 << std::endl;
448 SecondaryStructuresStatusMap[j + k].setStatus(Helix);
455 /* Не знаю зач они в дссп так сделали, этож полное говно */
457 for(std::size_t i {1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
458 if (static_cast<int>(SecondaryStructuresStatusMap[i].getStatus()) <= static_cast<int>(secondaryStructureTypes::Turn)){
460 for(const turnsTypes &j : {turnsTypes::Turn_3, turnsTypes::Turn_4, turnsTypes::Turn_5}){
461 std::size_t stride {static_cast<std::size_t>(j) + 3};
462 for(std::size_t k {1}; k < stride and !isTurn; ++k){
463 isTurn = (i >= k) && (SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start || SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start_AND_End) ;
467 // std::cout << "Resi " << i << " is Turn" << std::endl;
468 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Turn);
475 std::string secondaryStructures::patternSearch(){
478 analyzeBridgesAndStrandsPatterns();
479 analyzeTurnsAndHelicesPatterns();
481 // for(std::size_t i {0}; i < ResInfoMap->size(); ++i){
482 // std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) << std::endl;
485 // std::cout.precision(5);
486 // for(std::size_t i{0}; i < ResInfoMap->size(); ++i, std::cout << std::endl << std::endl){
487 // std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) ;
488 // if ( (*ResInfoMap)[i].donor[0] != nullptr ){
489 // std::cout << " has donor[0] = " << (*ResInfoMap)[i].donor[0]->nr << " " << *((*ResInfoMap)[i].donor[0]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[0] << " and" ;
492 // std::cout << " has no donor[0] and" ;
494 // if ( (*ResInfoMap)[i].acceptor[0] != nullptr ){
495 // std::cout << " has acceptor[0] = " << (*ResInfoMap)[i].acceptor[0]->nr << " " << *((*ResInfoMap)[i].acceptor[0]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[0] ;
498 // std::cout << " has no acceptor[0]" ;
500 // std::cout << std::endl << "Also, " << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name);
501 // if ( (*ResInfoMap)[i].donor[1] != nullptr ){
502 // std::cout << " has donor[1] = " << (*ResInfoMap)[i].donor[1]->nr << " " << *((*ResInfoMap)[i].donor[1]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[1] << " and" ;
505 // std::cout << " has no donor[1] and" ;
507 // if ( (*ResInfoMap)[i].acceptor[1] != nullptr ){
508 // std::cout << " has acceptor[1] = " << (*ResInfoMap)[i].acceptor[1]->nr << " " << *((*ResInfoMap)[i].acceptor[1]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[1] ;
511 // std::cout << " has no acceptor[1]" ;
517 for(std::size_t i {static_cast<std::size_t>(secondaryStructureTypes::Bend)}; i != static_cast<std::size_t>(secondaryStructureTypes::Count); ++i){
518 for(std::size_t j {0}; j < SecondaryStructuresStatusMap.size(); ++j){
519 if (SecondaryStructuresStatusMap[j].getStatus(static_cast<secondaryStructureTypes>(i))){
520 SecondaryStructuresStringLine[j] = secondaryStructureTypeNames[i] ;
527 if(SecondaryStructuresStatusMap.size() > 1){
528 for(std::size_t i {0}, linefactor{1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
529 if( SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) && SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break) ){
530 if(SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]) ){
531 SecondaryStructuresStringLine.insert(SecondaryStructuresStringLine.begin() + i + linefactor, secondaryStructureTypeNames[secondaryStructureTypes::Break]);
537 return SecondaryStructuresStringLine;
540 secondaryStructures::~secondaryStructures(){
541 SecondaryStructuresStatusMap.resize(0);
542 SecondaryStructuresStringLine.resize(0);
545 DsspTool::DsspStorage::DsspStorage(){
546 storaged_data.resize(0);
549 void DsspTool::DsspStorage::clearAll(){
550 storaged_data.resize(0);
553 std::mutex DsspTool::DsspStorage::mx;
555 void DsspTool::DsspStorage::storageData(int frnr, std::string data){
556 std::lock_guard<std::mutex> guardian(mx);
557 std::pair<int, std::string> datapair(frnr, data);
558 storaged_data.push_back(datapair);
561 std::vector<std::pair<int, std::string>> DsspTool::DsspStorage::returnData(){
562 std::sort(storaged_data.begin(), storaged_data.end());
563 return storaged_data;
566 void alternateNeighborhoodSearch::setCutoff(const real &cutoff_init){
567 cutoff = cutoff_init;
570 void alternateNeighborhoodSearch::FixAtomCoordinates(real &coordinate, const real vector_length){
571 while (coordinate < 0) {
572 coordinate += vector_length;
574 while (coordinate >= vector_length) {
575 coordinate -= vector_length;
579 void alternateNeighborhoodSearch::ReCalculatePBC(int &x, const int &x_max) {
588 void alternateNeighborhoodSearch::GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
589 rvec coordinates, box_vector_length;
590 num_of_miniboxes.resize(0);
591 num_of_miniboxes.resize(3);
592 for (std::size_t i{XX}; i <= ZZ; ++i) {
593 box_vector_length[i] = std::sqrt(
594 std::pow(fr.box[i][XX], 2) + std::pow(fr.box[i][YY], 2) + std::pow(fr.box[i][ZZ], 2));
595 num_of_miniboxes[i] = std::floor((box_vector_length[i] / cutoff)) + 1;
597 MiniBoxesMap.resize(0);
598 MiniBoxesReverseMap.resize(0);
599 MiniBoxesMap.resize(num_of_miniboxes[XX], std::vector<std::vector<std::vector<std::size_t> > >(
600 num_of_miniboxes[YY], std::vector<std::vector<std::size_t> >(
601 num_of_miniboxes[ZZ], std::vector<std::size_t>(
603 MiniBoxesReverseMap.resize(IndexMap.size(), std::vector<std::size_t>(3));
604 for (std::vector<ResInfo>::const_iterator i {IndexMap.begin()}; i != IndexMap.end(); ++i) {
605 for (std::size_t j{XX}; j <= ZZ; ++j) {
606 coordinates[j] = fr.x[i->getIndex(backboneAtomTypes::AtomCA)][j];
607 FixAtomCoordinates(coordinates[j], box_vector_length[j]);
609 MiniBoxesMap[std::floor(coordinates[XX] / cutoff)][std::floor(coordinates[YY] / cutoff)][std::floor(
610 coordinates[ZZ] / cutoff)].push_back(i - IndexMap.begin());
611 for (std::size_t j{XX}; j <= ZZ; ++j){
612 MiniBoxesReverseMap[i - IndexMap.begin()][j] = std::floor(coordinates[j] / cutoff);
617 void alternateNeighborhoodSearch::AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
618 GetMiniBoxesMap(fr, IndexMap);
619 MiniBoxSize[XX] = MiniBoxesMap.size();
620 MiniBoxSize[YY] = MiniBoxesMap.front().size();
621 MiniBoxSize[ZZ] = MiniBoxesMap.front().front().size();
623 PairMap.resize(IndexMap.size(), std::vector<bool>(IndexMap.size(), false));
624 ResiI = PairMap.begin();
625 ResiJ = ResiI->begin();
627 for (std::vector<ResInfo>::const_iterator i = IndexMap.begin(); i != IndexMap.end(); ++i){
628 for (offset[XX] = -1; offset[XX] <= 1; ++offset[XX]) {
629 for (offset[YY] = -1; offset[YY] <= 1; ++offset[YY]) {
630 for (offset[ZZ] = -1; offset[ZZ] <= 1; ++offset[ZZ]) {
631 for (std::size_t k{XX}; k <= ZZ; ++k) {
632 fixBox[k] = MiniBoxesReverseMap[i - IndexMap.begin()][k] + offset[k];
633 ReCalculatePBC(fixBox[k], MiniBoxSize[k]);
635 for (std::size_t j{0}; j < MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]].size(); ++j) {
636 if ( (i - IndexMap.begin()) != MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]){
637 PairMap[i - IndexMap.begin()][MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]] = true;
638 PairMap[MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]][i - IndexMap.begin()] = true;
647 bool alternateNeighborhoodSearch::findNextPair(){
653 for(; ResiI != PairMap.end(); ++ResiI, ResiJ = ResiI->begin() ){
654 for(; ResiJ != ResiI->end(); ++ResiJ){
656 resiIpos = ResiI - PairMap.begin();
657 resiJpos = ResiJ - ResiI->begin();
658 if ( ResiJ != ResiI->end() ){
661 else if (ResiI != PairMap.end()) {
663 ResiJ = ResiI->begin();
676 std::size_t alternateNeighborhoodSearch::getResiI() const {
680 std::size_t alternateNeighborhoodSearch::getResiJ() const {
685 DsspTool::DsspStorage DsspTool::Storage;
687 DsspTool::DsspTool(){
690 void DsspTool::calculateBends(const t_trxframe &fr, const t_pbc *pbc)
692 const float benddegree{ 70.0 }, maxdist{ 2.5 };
693 float degree{ 0 }, vdist{ 0 }, vprod{ 0 };
694 gmx::RVec a{ 0, 0, 0 }, b{ 0, 0, 0 };
695 for (std::size_t i{ 0 }; i + 1 < IndexMap.size(); ++i)
697 if (CalculateAtomicDistances(static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
698 static_cast<int>(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)),
703 PatternSearch.SecondaryStructuresStatusMap[i].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i + 1]);
704 PatternSearch.SecondaryStructuresStatusMap[i + 1].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i]);
706 // std::cout << "Break between " << i + 1 << " and " << i + 2 << std::endl;
709 for (std::size_t i{ 2 }; i + 2 < IndexMap.size() ; ++i)
711 if (PatternSearch.SecondaryStructuresStatusMap[i - 2].getStatus(secondaryStructureTypes::Break) ||
712 PatternSearch.SecondaryStructuresStatusMap[i - 1].getStatus(secondaryStructureTypes::Break) ||
713 PatternSearch.SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) ||
714 PatternSearch.SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break)
719 for (int j{ 0 }; j < 3; ++j)
721 a[j] = fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j]
722 - fr.x[IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA)][j];
723 b[j] = fr.x[IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA)][j]
724 - fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j];
726 vdist = (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
727 vprod = CalculateAtomicDistances(IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA),
728 IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
731 * gmx::c_angstrom / gmx::c_nano
732 * CalculateAtomicDistances(IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
733 IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA),
736 * gmx::c_angstrom / gmx::c_nano;
737 degree = std::acos(vdist / vprod) * gmx::c_rad2Deg;
738 if (degree > benddegree)
740 PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bend);
745 float DsspTool::CalculateDihedralAngle(const int &A, const int &B, const int &C, const int &D, const t_trxframe &fr, const t_pbc *pbc){
746 float result{360}, u{}, v{};
747 gmx::RVec v12{}, v43{}, z{}, p{}, x{}, y{};
748 pbc_dx(pbc, fr.x[A], fr.x[B], v12.as_vec());
749 pbc_dx(pbc, fr.x[D], fr.x[C], v43.as_vec());
750 pbc_dx(pbc, fr.x[B], fr.x[C], z.as_vec());
752 for(std::size_t j {XX}; j <= ZZ; ++j){
753 v12[j] *= gmx::c_nm2A;
754 v43[j] *= gmx::c_nm2A;
757 for(std::size_t i{XX}, j{j + 1}, k{i + 2}; i <= ZZ; ++i, ++j, ++k){
764 p[i] = (z[j] * v12[k]) - (z[k] * v12[j]);
765 x[i] = (z[j] * v43[k]) - (z[k] * v43[j]);
768 for(std::size_t i{XX}, j{j + 1}, k{i + 2}; i <= ZZ; ++i, ++j, ++k){
775 y[i] = (z[j] * x[k]) - (z[k] * x[j]);
778 // std::cout << "v12 = " << v12[0] << ", " << v12[1] << ", " << v12[2] << std::endl;
779 // std::cout << "v43 = " << v43[0] << ", " << v43[1] << ", " << v43[2] << std::endl;
780 // std::cout << "z = " << z[0] << ", " << z[1] << ", " << z[2] << std::endl;
781 // std::cout << "p = " << p[0] << ", " << p[1] << ", " << p[2] << std::endl;
782 // std::cout << "x = " << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
783 // std::cout << "y = " << y[0] << ", " << y[1] << ", " << y[2] << std::endl;
785 u = (x[XX] * x[XX]) + (x[YY] * x[YY]) + (x[ZZ] * x[ZZ]);
786 v = (y[XX] * y[XX]) + (y[YY] * y[YY]) + (y[ZZ] * y[ZZ]);
788 // std::cout << "u = " << u << std::endl;
789 // std::cout << "v = " << v << std::endl;
791 if (u > 0 and v > 0){
792 u = ((p[XX] * x[XX]) + (p[YY] * x[YY]) + (p[ZZ] * x[ZZ])) / std::sqrt(u);
793 v = ((p[XX] * y[XX]) + (p[YY] * y[YY]) + (p[ZZ] * y[ZZ])) / std::sqrt(v);
794 // std::cout << "new u = " << u << std::endl;
795 // std::cout << "new v = " << v << std::endl;
796 if (u != 0 or v != 0){
797 result = std::atan2(v, u) * gmx::c_rad2Deg;
798 // std::cout << "result = " << result << std::endl;
804 void DsspTool::calculateDihedrals(const t_trxframe &fr, const t_pbc *pbc){
805 const float epsilon = 29;
806 const float phi_min = -75 - epsilon; // -104
807 const float phi_max = -75 + epsilon; // -46
808 const float psi_min = 145 - epsilon; // 116
809 const float psi_max = 145 + epsilon; // 176
810 std::vector<float> phi(IndexMap.size(), 360), psi(IndexMap.size(), 360);
813 // phi.resize(IndexMap.size(), 360);
814 // psi.resize(IndexMap.size(), 360);
816 for (std::size_t i = 1; i + 1 < IndexMap.size(); ++i){ // TODO add index verifictaion (check if those atom indexes exist)
817 // std::cout << "For resi " << i << " :" << std::endl;
818 phi[i] = CalculateDihedralAngle(static_cast<int>(IndexMap[i - 1].getIndex(backboneAtomTypes::AtomC)),
819 static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomN)),
820 static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomCA)),
821 static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
824 psi[i] = CalculateDihedralAngle(static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomN)),
825 static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomCA)),
826 static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
827 static_cast<int>(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)),
830 // std::cout << "For " << i << " phi = " << phi[i] << ", psi = " << psi[i] << std::endl;
831 // std::cout << "phi[" << i << "] = " << phi[i] << std::endl;
832 // std::cout << "psi[" << i << "] = " << psi[i] << std::endl;
835 for (std::size_t i = 1; i + 3 < IndexMap.size(); ++i){
836 switch (initParams.pp_stretch){
838 if (phi_min > phi[i] or phi[i] > phi_max or
839 phi_min > phi[i + 1] or phi[i + 1]> phi_max){
843 if (psi_min > psi[i] or psi[i] > psi_max or
844 psi_min > psi[i + 1] or psi[i + 1] > psi_max){
848 switch (PatternSearch.SecondaryStructuresStatusMap[i].getStatus(turnsTypes::Turn_PP)){
849 case HelixPositions::None:
850 PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start, turnsTypes::Turn_PP);
853 case HelixPositions::End:
854 PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start_AND_End, turnsTypes::Turn_PP);
861 PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(HelixPositions::End, turnsTypes::Turn_PP);
862 /* Пропустил проверку того, что заменяемая ак - петля */
863 PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Helix_PP);
864 PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(secondaryStructureTypes::Helix_PP);
868 if (phi_min > phi[i] or phi[i] > phi_max or
869 phi_min > phi[i + 1] or phi[i + 1]> phi_max or
870 phi_min > phi[i + 2] or phi[i + 2]> phi_max){
874 if (psi_min > psi[i] or psi[i] > psi_max or
875 psi_min > psi[i + 1] or psi[i + 1] > psi_max or
876 psi_min > psi[i + 2] or psi[i + 2] > psi_max){
880 switch (PatternSearch.SecondaryStructuresStatusMap[i].getStatus(turnsTypes::Turn_PP)){
881 case HelixPositions::None:
882 PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start, turnsTypes::Turn_PP);
885 case HelixPositions::End:
886 PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start_AND_End, turnsTypes::Turn_PP);
893 PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(HelixPositions::Middle, turnsTypes::Turn_PP);
894 PatternSearch.SecondaryStructuresStatusMap[i + 2].setStatus(HelixPositions::End, turnsTypes::Turn_PP);
895 /* Пропустил проверку того, что заменяемая ак - петля */
896 PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Helix_PP);
897 PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(secondaryStructureTypes::Helix_PP);
898 PatternSearch.SecondaryStructuresStatusMap[i + 2].setStatus(secondaryStructureTypes::Helix_PP);
903 throw std::runtime_error("Unsupported stretch length");
909 void DsspTool::calculateHBondEnergy(ResInfo& Donor,
911 const t_trxframe& fr,
915 * DSSP uses eq from dssp 2.x
916 * kCouplingConstant = 27.888, // = 332 * 0.42 * 0.2
917 * E = k * (1/rON + 1/rCH - 1/rOH - 1/rCN) where CO comes from one AA and NH from another
921 * For the note, H-Bond Donor is N-H («Donor of H») and H-Bond Acceptor is C=O («Acceptor of H»)
925 if (CalculateAtomicDistances(
926 Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
927 >= minimalCAdistance)
932 const float kCouplingConstant = 27.888;
933 const float minimalAtomDistance{ 0.5 },
935 float HbondEnergy{ 0 };
936 float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 };
938 // std::cout << "For Donor №" << Donor.info->nr - 1 << " and Accpetor №" << Acceptor.info->nr - 1 << std::endl;
940 if( !(Donor.is_proline) && (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
941 && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || initParams.addHydrogens ) ) ){ // TODO
942 distanceNO = CalculateAtomicDistances(
943 Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
944 distanceNC = CalculateAtomicDistances(
945 Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
946 if (initParams.addHydrogens){
947 if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
949 float prevCODist {CalculateAtomicDistances(Donor.prevResi->getIndex(backboneAtomTypes::AtomC), Donor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)};
950 for (int i{XX}; i <= ZZ; ++i){
951 float prevCO = fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
952 atomH[i] = fr.x[Donor.getIndex(backboneAtomTypes::AtomH)][i]; // Но на самом деле берутся координаты N
953 atomH[i] += prevCO / prevCODist;
955 distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
956 distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
959 distanceHO = distanceNO;
960 distanceHC = distanceNC;
964 distanceHO = CalculateAtomicDistances(
965 Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
966 distanceHC = CalculateAtomicDistances(
967 Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
969 if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
970 || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
972 HbondEnergy = minEnergy;
977 * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
980 // std::cout << "CA-CA distance: " << CalculateAtomicDistances(
981 // Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc) << std::endl;
982 // std::cout << "N-O distance: " << distanceNO << std::endl;
983 // std::cout << "N-C distance: " << distanceNC << std::endl;
984 // std::cout << "H-O distance: " << distanceHO << std::endl;
985 // std::cout << "H-C distance: " << distanceHC << std::endl;
987 HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
989 if ( HbondEnergy < minEnergy ){
990 HbondEnergy = minEnergy;
993 // std::cout << "Calculated energy = " << HbondEnergy << std::endl;
996 // std::cout << "Donor Is Proline" << std::endl;
999 if (HbondEnergy < Donor.acceptorEnergy[0]){
1000 Donor.acceptor[1] = Donor.acceptor[0];
1001 Donor.acceptor[0] = Acceptor.info;
1002 Donor.acceptorEnergy[0] = HbondEnergy;
1004 else if (HbondEnergy < Donor.acceptorEnergy[1]){
1005 Donor.acceptor[1] = Acceptor.info;
1006 Donor.acceptorEnergy[1] = HbondEnergy;
1009 if (HbondEnergy < Acceptor.donorEnergy[0]){
1010 Acceptor.donor[1] = Acceptor.donor[0];
1011 Acceptor.donor[0] = Donor.info;
1012 Acceptor.donorEnergy[0] = HbondEnergy;
1014 else if (HbondEnergy < Acceptor.donorEnergy[1]){
1015 Acceptor.donor[1] = Donor.info;
1016 Acceptor.donorEnergy[1] = HbondEnergy;
1021 /* Calculate Distance From B to A */
1022 float DsspTool::CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
1024 gmx::RVec r{ 0, 0, 0 };
1025 pbc_dx(pbc, fr.x[A], fr.x[B], r.as_vec());
1026 return r.norm() * gmx::c_nm2A; // НЕ ТРОГАТЬ
1029 /* Calculate Distance From B to A, where A is only fake coordinates */
1030 float DsspTool::CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
1032 gmx::RVec r{ 0, 0, 0 };
1033 pbc_dx(pbc, A, fr.x[B], r.as_vec());
1034 return r.norm() * gmx::c_nm2A; // НЕ ТРОГАТЬ
1037 void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz)
1039 initParams = initParamz;
1040 ResInfo _backboneAtoms;
1042 std::string proLINE;
1043 int resicompare{ top.atoms()->atom[static_cast<std::size_t>(*(initParams.sel_.atomIndices().begin()))].resind };
1045 IndexMap.push_back(_backboneAtoms);
1046 IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
1047 proLINE = *(IndexMap[i].info->name);
1048 if( proLINE.compare("PRO") == 0 ){
1049 IndexMap[i].is_proline = true;
1052 for (gmx::ArrayRef<const int>::iterator ai{ initParams.sel_.atomIndices().begin() }; (ai != initParams.sel_.atomIndices().end()); ++ai){
1053 if (resicompare != top.atoms()->atom[static_cast<std::size_t>(*ai)].resind)
1056 resicompare = top.atoms()->atom[static_cast<std::size_t>(*ai)].resind;
1057 IndexMap.emplace_back(_backboneAtoms);
1058 IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
1059 proLINE = *(IndexMap[i].info->name);
1060 if( proLINE.compare("PRO") == 0 ){
1061 IndexMap[i].is_proline = true;
1065 std::string atomname(*(top.atoms()->atomname[static_cast<std::size_t>(*ai)]));
1066 if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA])
1068 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomCA)] = *ai;
1070 else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC])
1072 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomC)] = *ai;
1074 else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO])
1076 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomO)] = *ai;
1078 else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN])
1080 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomN)] = *ai;
1081 if (initParamz.addHydrogens == true){
1082 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
1085 else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH] && initParamz.addHydrogens == false) // Юзать водород в структуре
1087 IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
1092 // if( atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO]
1093 // || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH]){
1094 // std::cout << "Atom " << atomname << " №" << *ai << " From Resi " << *(top.atoms()->resinfo[i].name) << " №" << resicompare << std::endl;
1098 for (std::size_t j {1}; j < IndexMap.size(); ++j){
1099 IndexMap[j].prevResi = &(IndexMap[j - 1]);
1101 IndexMap[j - 1].nextResi = &(IndexMap[j]);
1103 // std::cout << "Resi " << IndexMap[i].info->nr << *(IndexMap[i].info->name) << std::endl;
1104 // std::cout << "Prev resi is " << IndexMap[i].prevResi->info->nr << *(IndexMap[i].prevResi->info->name) << std::endl;
1105 // std::cout << "Prev resi's next resi is " << IndexMap[i - 1].nextResi->info->nr << *(IndexMap[i - 1].nextResi->info->name) << std::endl;
1106 // std::cout << IndexMap[j].prevResi->info->nr;
1107 // std::cout << *(IndexMap[j].prevResi->info->name) ;
1108 // std::cout << " have CA = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomCA) ;
1109 // std::cout << " C = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomC);
1110 // std::cout << " O = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomO);
1111 // std::cout << " N = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomN);
1112 // std::cout << " H = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomH) << std::endl;
1118 void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc)
1121 switch(initParams.NBS){
1122 case (NBSearchMethod::Classique): {
1124 // store positions of CA atoms to use them for nbSearch
1125 std::vector<gmx::RVec> positionsCA_;
1126 for (std::size_t i{ 0 }; i < IndexMap.size(); ++i)
1128 positionsCA_.emplace_back(fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)]);
1131 AnalysisNeighborhood nb_;
1132 nb_.setCutoff(initParams.cutoff_);
1133 AnalysisNeighborhoodPositions nbPos_(positionsCA_);
1134 gmx::AnalysisNeighborhoodSearch start = nb_.initSearch(pbc, nbPos_);
1135 gmx::AnalysisNeighborhoodPairSearch pairSearch = start.startPairSearch(nbPos_);
1136 gmx::AnalysisNeighborhoodPair pair;
1137 while (pairSearch.findNextPair(&pair))
1139 if(CalculateAtomicDistances(
1140 IndexMap[pair.refIndex()].getIndex(backboneAtomTypes::AtomCA), IndexMap[pair.testIndex()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
1141 < minimalCAdistance){
1142 calculateHBondEnergy(IndexMap[pair.refIndex()], IndexMap[pair.testIndex()], fr, pbc);
1143 if (IndexMap[pair.testIndex()].info != IndexMap[pair.refIndex() + 1].info){
1144 calculateHBondEnergy(IndexMap[pair.testIndex()], IndexMap[pair.refIndex()], fr, pbc);
1151 case (NBSearchMethod::Experimental): { // TODO FIX
1153 alternateNeighborhoodSearch as_;
1155 as_.setCutoff(initParams.cutoff_);
1157 as_.AltPairSearch(fr, IndexMap);
1159 while (as_.findNextPair()){
1160 if(CalculateAtomicDistances(
1161 IndexMap[as_.getResiI()].getIndex(backboneAtomTypes::AtomCA), IndexMap[as_.getResiJ()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
1162 < minimalCAdistance){
1163 calculateHBondEnergy(IndexMap[as_.getResiI()], IndexMap[as_.getResiJ()], fr, pbc);
1164 if (IndexMap[as_.getResiJ()].info != IndexMap[as_.getResiI() + 1].info){
1165 calculateHBondEnergy(IndexMap[as_.getResiJ()], IndexMap[as_.getResiI()], fr, pbc);
1174 for(std::vector<ResInfo>::iterator Donor {IndexMap.begin()}; Donor != IndexMap.end() ; ++Donor){
1175 for(std::vector<ResInfo>::iterator Acceptor {Donor + 1} ; Acceptor != IndexMap.end() ; ++Acceptor){
1176 if(CalculateAtomicDistances(
1177 Donor->getIndex(backboneAtomTypes::AtomCA), Acceptor->getIndex(backboneAtomTypes::AtomCA), fr, pbc)
1178 < minimalCAdistance){
1179 calculateHBondEnergy(*Donor, *Acceptor, fr, pbc);
1180 if (Acceptor != Donor + 1){
1181 calculateHBondEnergy(*Acceptor, *Donor, fr, pbc);
1191 // for(std::size_t i {0}; i < IndexMap.size(); ++i){
1192 // std::cout << IndexMap[i].info->nr << " " << *(IndexMap[i].info->name) << std::endl;
1195 PatternSearch.initiateSearch(IndexMap, initParams.PPHelices, initParams.pp_stretch);
1196 calculateBends(fr, pbc);
1197 calculateDihedrals(fr, pbc);
1198 Storage.storageData(frnr, PatternSearch.patternSearch());
1202 std::vector<std::pair<int, std::string>> DsspTool::getData(){
1203 return Storage.returnData();
1206 } // namespace analysismodules