src/dssptools.cpp

   1 /*
   2 * This file is part of the GROMACS molecular simulation package.
   3 *
   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.
   8 *
   9 * GROMACS is free software; you can redistribute it and/or
  10 * modify it under the terms of the GNU Lesser General Public License
  11 * as published by the Free Software Foundation; either version 2.1
  12 * of the License, or (at your option) any later version.
  13 *
  14 * GROMACS is distributed in the hope that it will be useful,
  15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17 * Lesser General Public License for more details.
  18 *
  19 * You should have received a copy of the GNU Lesser General Public
  20 * License along with GROMACS; if not, see
  21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23 *
  24 * If you want to redistribute modifications to GROMACS, please
  25 * consider that scientific software is very special. Version
  26 * control is crucial - bugs must be traceable. We will be happy to
  27 * consider code for inclusion in the official distribution, but
  28 * derived work must not be called official GROMACS. Details are found
  29 * in the README & COPYING files - if they are missing, get the
  30 * official version at http://www.gromacs.org.
  31 *
  32 * To help us fund GROMACS development, we humbly ask that you cite
  33 * the research papers on the package. Check out http://www.gromacs.org.
  34 */
  35 /*! \internal \file
  36 * \brief
  37 * Implements gmx::analysismodules::Trajectory.
  38 *
  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
  43 */
  44
  45 /*
  46     There's something wrong with energy calculations of redidues with E ≈ -0
  47 */
  48
  49
  50 #include "dssptools.h"
  51
  52 #include <algorithm>
  53 #include "gromacs/math/units.h"
  54
  55 #include "gromacs/pbcutil/pbc.h"
  56 #include <gromacs/trajectoryanalysis.h>
  57 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  58 #include <set>
  59 #include <fstream>
  60 #include <mutex>
  61 #include <iostream>
  62
  63 namespace gmx
  64 {
  65
  66 namespace analysismodules
  67 {
  68
  69 //void ResInfo::setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex)
  70 //{
  71 //   _ResInfo.at(static_cast<std::size_t>(atomTypeName)) = atomIndex;
  72 //}
  73
  74 //std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const
  75 //{
  76 //   return _ResInfo[static_cast<std::size_t>(atomTypeName)];
  77 //}
  78
  79 std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const{
  80     return _backboneIndices[static_cast<std::size_t>(atomTypeName)];
  81 }
  82
  83 secondaryStructures::secondaryStructures(){
  84 }
  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(), '~');
  93 }
  94
  95 void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName){
  96     SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)] = true;
  97     SecondaryStructuresStatus = secondaryStructureTypeName;
  98 }
  99
 100 void secondaryStructures::secondaryStructuresData::setStatus(const HelixPositions helixPosition, const turnsTypes turn){
 101     TurnsStatusArray[static_cast<std::size_t>(turn)] = helixPosition;
 102 }
 103
 104 bool secondaryStructures::secondaryStructuresData::getStatus(const secondaryStructureTypes secondaryStructureTypeName) const{
 105     return SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)];
 106 }
 107
 108 bool secondaryStructures::secondaryStructuresData::isBreakPartnerWith(const secondaryStructuresData *partner) const{
 109     return breakPartners[0] == partner || breakPartners[1] == partner;
 110 }
 111
 112 HelixPositions secondaryStructures::secondaryStructuresData::getStatus(const turnsTypes turn) const{
 113     return TurnsStatusArray[static_cast<std::size_t>(turn)];
 114 }
 115
 116 secondaryStructureTypes secondaryStructures::secondaryStructuresData::getStatus() const{
 117     return  SecondaryStructuresStatus;
 118 }
 119
 120 void secondaryStructures::secondaryStructuresData::setBreak(secondaryStructuresData *breakPartner){
 121     if (breakPartners[0] == nullptr){
 122         breakPartners[0] = breakPartner;
 123     }
 124     else{
 125         breakPartners[1] = breakPartner;
 126     }
 127     setStatus(secondaryStructureTypes::Break);
 128 }
 129
 130 void secondaryStructures::secondaryStructuresData::setBridge(secondaryStructuresData *bridgePartner, std::size_t bridgePartnerIndex, bridgeTypes bridgeType){
 131     if(bridgeType == bridgeTypes::ParallelBridge){
 132         bridgePartners[0] = bridgePartner;
 133         bridgePartnersIndexes[0] = bridgePartnerIndex;
 134     }
 135     else if (bridgeType == bridgeTypes::AntiParallelBridge){
 136         bridgePartners[1] = bridgePartner;
 137         bridgePartnersIndexes[1] = bridgePartnerIndex;
 138     }
 139 }
 140
 141 bool secondaryStructures::secondaryStructuresData::hasBridges() const{
 142     return bridgePartners[0] || bridgePartners[1];
 143
 144 }
 145
 146 bool secondaryStructures::secondaryStructuresData::hasBridges(bridgeTypes bridgeType) const{
 147     if(bridgeType == bridgeTypes::ParallelBridge){
 148         return bridgePartners[0] != nullptr;
 149     }
 150     else if (bridgeType == bridgeTypes::AntiParallelBridge){
 151         return bridgePartners[1] != nullptr;
 152     }
 153
 154     return false;
 155
 156 }
 157
 158 bool secondaryStructures::secondaryStructuresData::isBridgePartnerWith(secondaryStructuresData *bridgePartner, bridgeTypes bridgeType) const{
 159     if(bridgeType == bridgeTypes::ParallelBridge){
 160         return bridgePartners[0] == bridgePartner;
 161     }
 162     else if (bridgeType == bridgeTypes::AntiParallelBridge){
 163         return bridgePartners[1] == bridgePartner;
 164     }
 165     return false;
 166 }
 167
 168 std::size_t secondaryStructures::secondaryStructuresData::getBridgePartnerIndex(bridgeTypes bridgeType) const{
 169     if (bridgeType == bridgeTypes::ParallelBridge){
 170         return bridgePartnersIndexes[0];
 171     }
 172     return bridgePartnersIndexes[1];
 173 }
 174
 175 secondaryStructures::secondaryStructuresData secondaryStructures::secondaryStructuresData::getBridgePartner(bridgeTypes bridgeType) const{
 176     if(bridgeType == bridgeTypes::ParallelBridge){
 177         return *(bridgePartners[0]);
 178     }
 179     return *(bridgePartners[1]);
 180 }
 181
 182 bool secondaryStructures::hasHBondBetween(std::size_t Donor, std::size_t Acceptor) const{
 183     if( (*ResInfoMap)[Donor].acceptor[0] == nullptr ||
 184         (*ResInfoMap)[Donor].acceptor[1] == nullptr ||
 185         (*ResInfoMap)[Acceptor].info == nullptr ){
 186 //        std::cout << "Bad hbond check. Reason(s): " ;
 187 //        if ( (*ResInfoMap)[Donor].acceptor[0] == nullptr ){
 188 //            std::cout << "Donor has no acceptor[0]; ";
 189 //        }
 190 //        if ( (*ResInfoMap)[Donor].acceptor[1] == nullptr ){
 191 //            std::cout << "Donor has no acceptor[1]; ";
 192 //        }
 193 //        if ( (*ResInfoMap)[Acceptor].info == nullptr ){
 194 //            std::cout << "No info about acceptor; ";
 195 //        }
 196 //        std::cout << std::endl;
 197         return false;
 198     }
 199 //    else if (!( (*ResInfoMap)[Acceptor].donor[0] == nullptr ||
 200 //              (*ResInfoMap)[Acceptor].donor[1] == nullptr ||
 201 //              (*ResInfoMap)[Donor].info == nullptr )) {
 202 //        std::cout << "Comparing DONOR №" << (*ResInfoMap)[Donor].info->nr << " And ACCEPTOR №" << (*ResInfoMap)[Acceptor].info->nr << ": " << std::endl;
 203 //        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;
 204 //        std::cout << "Donor's acceptors' energy are = " << (*ResInfoMap)[Donor].acceptorEnergy[0] << ", " << (*ResInfoMap)[Donor].acceptorEnergy[1] << std::endl;
 205 //        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;
 206 //        std::cout << "Acceptors's donors' energy are = " << (*ResInfoMap)[Acceptor].donorEnergy[0] << ", " << (*ResInfoMap)[Acceptor].donorEnergy[1] << std::endl;
 207 //        if( ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
 208 //                ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff ) ){
 209 //            std::cout << "HBond Exist" << std::endl;
 210 //        }
 211 //    }
 212 //    else {
 213 //        std::cout << "Bad hbond check. Reason(s): " ;
 214 //        if ( (*ResInfoMap)[Acceptor].donor[0] == nullptr ){
 215 //            std::cout << "Acceptor has no donor[0]; ";
 216 //        }
 217 //        if ( (*ResInfoMap)[Acceptor].donor[1] == nullptr ){
 218 //            std::cout << "Acceptor has no donor[1]; ";
 219 //        }
 220 //        if ( (*ResInfoMap)[Donor].info == nullptr ){
 221 //            std::cout << "No info about donor; ";
 222 //        }
 223 //        std::cout << std::endl;
 224 //    }
 225
 226     return ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
 227            ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff );
 228
 229
 230 }
 231
 232 bool secondaryStructures::NoChainBreaksBetween(std::size_t Resi1, std::size_t Resi2) const{
 233     std::size_t i{Resi1}, j{Resi2}; // From i to j → i <= j
 234     if ( i > j ){
 235         std::swap(i, j);
 236     }
 237
 238     for (; i != j; ++i){
 239         if ( SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]) ){
 240             std::cout << "Patternsearch has detected a CHAINBREAK between " << Resi1 << " and " << Resi2 << std::endl;
 241             return false;
 242         }
 243     }
 244     return true;
 245 }
 246
 247 bridgeTypes secondaryStructures::calculateBridge(std::size_t i, std::size_t j) const{
 248     if( i < 1 || j < 1 || i + 1 >= ResInfoMap->size() || j + 1 >= ResInfoMap->size() ){ // Protection from idiotz, не обязательно нужен
 249         return bridgeTypes::None;
 250     }
 251
 252     ResInfo a{(*ResInfoMap)[i]}, b{(*ResInfoMap)[j]};
 253
 254     if(NoChainBreaksBetween(i - 1, i + 1) && NoChainBreaksBetween(j - 1, j + 1) && a.prevResi && a.nextResi && b.prevResi && b.nextResi){
 255         if((hasHBondBetween(i + 1, j) && hasHBondBetween(j, i - 1)) || (hasHBondBetween(j + 1, i) && hasHBondBetween(i, j - 1)) ){
 256             return bridgeTypes::ParallelBridge;
 257         }
 258         else if((hasHBondBetween(i + 1, j - 1) && hasHBondBetween(j + 1, i - 1)) || (hasHBondBetween(j, i) && hasHBondBetween(i, j)) ){
 259             return bridgeTypes::AntiParallelBridge;
 260         }
 261     }
 262     return bridgeTypes::None;
 263 }
 264
 265 void secondaryStructures::analyzeBridgesAndStrandsPatterns(){
 266 //  Calculate Bridgez
 267     for(std::size_t i {1}; i + 4 < SecondaryStructuresStatusMap.size(); ++i){
 268         for(std::size_t j {i + 3}; j + 1 < SecondaryStructuresStatusMap.size(); ++j ){
 269             switch(calculateBridge(i, j)){
 270             case bridgeTypes::ParallelBridge : {
 271                 SecondaryStructuresStatusMap[i].setBridge(&(SecondaryStructuresStatusMap[j]), j, bridgeTypes::ParallelBridge);
 272                 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bridge);
 273                 SecondaryStructuresStatusMap[j].setBridge(&(SecondaryStructuresStatusMap[i]), i, bridgeTypes::ParallelBridge);
 274                 SecondaryStructuresStatusMap[j].setStatus(secondaryStructureTypes::Bridge);
 275             }
 276             case bridgeTypes::AntiParallelBridge : {
 277                 SecondaryStructuresStatusMap[i].setBridge(&(SecondaryStructuresStatusMap[j]), j, bridgeTypes::AntiParallelBridge);
 278                 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bridge);
 279                 SecondaryStructuresStatusMap[j].setBridge(&(SecondaryStructuresStatusMap[i]), i, bridgeTypes::AntiParallelBridge);
 280                 SecondaryStructuresStatusMap[j].setStatus(secondaryStructureTypes::Bridge);
 281             }
 282             default :
 283                 continue;
 284             }
 285         }
 286     }
 287 //  Calculate Extended Strandz
 288     for(std::size_t i {1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
 289         bool is_Estrand {false};
 290         for(std::size_t j {2}; j != 0; --j){
 291             for(const bridgeTypes &bridgeType : {bridgeTypes::ParallelBridge, bridgeTypes::AntiParallelBridge}){
 292                 if (SecondaryStructuresStatusMap[i].hasBridges(bridgeType) && SecondaryStructuresStatusMap[i + j].hasBridges(bridgeType) ){
 293                     std::size_t i_partner{SecondaryStructuresStatusMap[i].getBridgePartnerIndex(bridgeType)}, j_partner{SecondaryStructuresStatusMap[i + j].getBridgePartnerIndex(bridgeType)}, second_strand{};
 294                     if ( abs(i_partner - j_partner) < 6){
 295                         if (i_partner < j_partner){
 296                             second_strand = i_partner;
 297                         }
 298                         else{
 299                             second_strand = j_partner;
 300                         }
 301                         for(std::size_t k{abs(i_partner - j_partner)}; k >= 0; --k){
 302                             SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Strand);
 303                             std::cout << second_strand + k << " is strand" << std::endl;
 304                         }
 305                         is_Estrand = true;
 306                     }
 307
 308                 }
 309             }
 310             if (is_Estrand){
 311                 for(;j >= 0; --j){
 312                     SecondaryStructuresStatusMap[i + j].setStatus(secondaryStructureTypes::Strand);
 313                     std::cout << i + j << " is strand" << std::endl;
 314                 }
 315                 break;
 316             }
 317         }
 318     }
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331 //    for (std::size_t i{ 1 }; i < HBondsMap.front().size() - 1; ++i){
 332 //        for (std::size_t j{ 1 }; j < HBondsMap.front().size() - 1; ++j){
 333 //            if (std::abs(static_cast<int>(i) - static_cast<int>(j)) > 2){
 334 //                if ((HBondsMap[i - 1][j] && HBondsMap[j][i + 1])    ||
 335 //                    (HBondsMap[j - 1][i] && HBondsMap[i][j + 1])){
 336 //                    Bridge[i].push_back(j);
 337 //                }
 338 //                if ((HBondsMap[i][j] && HBondsMap[j][i])    ||
 339 //                    (HBondsMap[i - 1][j + 1] && HBondsMap[j - 1][i + 1])){
 340 //                    AntiBridge[i].push_back(j);
 341 //                }
 342 //            }
 343 //        }
 344 //    }
 345 //    for (std::size_t i{ 0 }; i < HBondsMap.front().size(); ++i){
 346 //        if ((!Bridge[i].empty() || !AntiBridge[i].empty())){
 347 //            setStatus(i, secondaryStructureTypes::Bulge);
 348 //        }
 349 //    }
 350 //    for (std::size_t i{ 2 }; i + 2 < HBondsMap.front().size(); ++i){
 351 //        for (std::size_t j { i - 2 }; j <= (i + 2); ++j){
 352 //            if (j == i){
 353 //                continue;
 354 //            }
 355 //            else {
 356 //                for (std::vector<bridgeTypes>::const_iterator bridge {Bridges.begin()}; bridge != Bridges.end(); ++bridge ){
 357 //                    if (!getBridge(*bridge)[i].empty() || !getBridge(*bridge)[j].empty()){
 358 //                        for (std::size_t i_resi{ 0 }; i_resi < getBridge(*bridge)[i].size(); ++i_resi){
 359 //                            for (std::size_t j_resi{ 0 }; j_resi < getBridge(*bridge)[j].size(); ++j_resi){
 360 //                                if (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
 361 //                                        - static_cast<int>(getBridge(*bridge)[j][j_resi]))
 362 //                                        && (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
 363 //                                        - static_cast<int>(getBridge(*bridge)[j][j_resi]))
 364 //                                        < 5)){
 365 //                                    if (j < i){
 366 //                                        for (std::size_t k{ 0 }; k <= i - j; ++k){
 367 //                                            setStatus(i + k, secondaryStructureTypes::Ladder);
 368 //                                        }
 369 //                                    }
 370 //                                    else{
 371 //                                        for (std::size_t k{ 0 }; k <= j - i; ++k){
 372 //                                            setStatus(i + k, secondaryStructureTypes::Ladder);
 373 //                                        }
 374 //                                    }
 375 //                                }
 376 //                            }
 377 //                        }
 378 //                    }
 379 //                }
 380 //            }
 381 //        }
 382 //    }
 383 }
 384
 385 void secondaryStructures::analyzeTurnsAndHelicesPatterns(){
 386     for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
 387         std::size_t stride {static_cast<std::size_t>(i) + 3};
 388         std::cout << "Testing Helix_" << stride << std::endl;
 389         for(std::size_t j {0}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
 390             std::cout << "Testing " << j << " and " << j + stride << std::endl;
 391             if ( hasHBondBetween(j + stride, j) && NoChainBreaksBetween(j, j + stride) ){
 392                 std::cout << j << " and " << j + stride << " has hbond!" << std::endl;
 393                 SecondaryStructuresStatusMap[j + stride].setStatus(HelixPositions::End, i);
 394
 395                 for (std::size_t k {1}; k < stride; ++k){
 396                     if( SecondaryStructuresStatusMap[j + k].getStatus(i) == HelixPositions::None ){
 397                         SecondaryStructuresStatusMap[j + k].setStatus(HelixPositions::Middle, i);
 398                         SecondaryStructuresStatusMap[j + k].setStatus(secondaryStructureTypes::Turn);
 399                     }
 400                 }
 401
 402                 if( SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::End ){
 403                     SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start_AND_End, i);
 404                 }
 405                 else {
 406                     SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start, i);
 407                 }
 408             }
 409         }
 410     }
 411
 412     for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
 413         std::size_t stride {static_cast<std::size_t>(i) + 3};
 414         for(std::size_t j {1}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
 415             if ( (SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start_AND_End ) &&
 416                  (SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start_AND_End ) ){
 417                 bool empty = true;
 418                 secondaryStructureTypes Helix;
 419                 switch(i){
 420                 case turnsTypes::Turn_3:
 421                     for (std::size_t k {0}; empty && k < stride; ++k){
 422                         empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_3);
 423                     }
 424                     Helix = secondaryStructureTypes::Helix_3;
 425                     break;
 426                 case turnsTypes::Turn_5:
 427                     for (std::size_t k {0}; empty && k < stride; ++k){
 428                         empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_5) || (PiHelixPreference && SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_4)); //TODO
 429                     }
 430                     Helix = secondaryStructureTypes::Helix_5;
 431                     break;
 432                 default:
 433                     Helix = secondaryStructureTypes::Helix_4;
 434                     break;
 435                 }
 436                 if ( empty || Helix == secondaryStructureTypes::Helix_4 ){
 437                     for(std::size_t k {0}; k < stride; ++k ){
 438                         SecondaryStructuresStatusMap[j + k].setStatus(Helix);
 439                     }
 440                 }
 441             }
 442         }
 443     }
 444
 445     /* Не знаю зач они в дссп так сделали, этож полное говно */
 446
 447 //    for(std::size_t i {1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
 448 //        if (SecondaryStructuresStatusMap[i].getStatus() == secondaryStructureTypes::Loop || SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Bend)){
 449 //            bool isTurn = false;
 450 //            for(const turnsTypes &j : {turnsTypes::Turn_3, turnsTypes::Turn_4, turnsTypes::Turn_5}){
 451 //                std::size_t stride {static_cast<std::size_t>(j) + 3};
 452 //                for(std::size_t k {1}; k < stride and !isTurn; ++k){
 453 //                    isTurn = (i >= k) && (SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start || SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start_AND_End) ;
 454 //                }
 455 //            }
 456 //            if (isTurn){
 457 //                SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Turn);
 458 //            }
 459 //        }
 460 //    }
 461
 462 }
 463
 464 void secondaryStructures::analyzePPHelicesPatterns(){}
 465
 466 std::string secondaryStructures::patternSearch(){
 467
 468
 469     analyzeBridgesAndStrandsPatterns();
 470     analyzeTurnsAndHelicesPatterns();
 471 //    analyzePPHelicesPatterns();
 472
 473 //    for(std::size_t i {0}; i < ResInfoMap->size(); ++i){
 474 //        std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) << std::endl;
 475 //    }
 476
 477 //    std::cout.precision(5);
 478 //    for(std::size_t i{0}; i < ResInfoMap->size(); ++i, std::cout << std::endl << std::endl){
 479 //        std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) ;
 480 //        if ( (*ResInfoMap)[i].donor[0] != nullptr ){
 481 //            std::cout << " has donor[0] = " << (*ResInfoMap)[i].donor[0]->nr << " " << *((*ResInfoMap)[i].donor[0]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[0] << " and" ;
 482 //        }
 483 //        else {
 484 //            std::cout << " has no donor[0] and" ;
 485 //        }
 486 //        if ( (*ResInfoMap)[i].acceptor[0] != nullptr ){
 487 //            std::cout << " has acceptor[0] = " << (*ResInfoMap)[i].acceptor[0]->nr << " " << *((*ResInfoMap)[i].acceptor[0]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[0] ;
 488 //        }
 489 //        else {
 490 //            std::cout << " has no acceptor[0]" ;
 491 //        }
 492 //        std::cout << std::endl << "Also, " << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name);
 493 //        if ( (*ResInfoMap)[i].donor[1] != nullptr ){
 494 //            std::cout << " has donor[1] = " << (*ResInfoMap)[i].donor[1]->nr << " " << *((*ResInfoMap)[i].donor[1]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[1] << " and" ;
 495 //        }
 496 //        else {
 497 //            std::cout << " has no donor[1] and" ;
 498 //        }
 499 //        if ( (*ResInfoMap)[i].acceptor[1] != nullptr ){
 500 //            std::cout << " has acceptor[1] = " << (*ResInfoMap)[i].acceptor[1]->nr << " " << *((*ResInfoMap)[i].acceptor[1]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[1] ;
 501 //        }
 502 //        else {
 503 //            std::cout << " has no acceptor[1]" ;
 504 //        }
 505 //    }
 506
 507     /*Write Data*/
 508
 509     for(std::size_t i {static_cast<std::size_t>(secondaryStructureTypes::Bend)}; i != static_cast<std::size_t>(secondaryStructureTypes::Count); ++i){
 510         for(std::size_t j {0}; j < SecondaryStructuresStatusMap.size(); ++j){
 511             if (SecondaryStructuresStatusMap[j].getStatus(static_cast<secondaryStructureTypes>(i))){
 512                 SecondaryStructuresStringLine[j] = secondaryStructureTypeNames[i] ;
 513             }
 514         }
 515     }
 516
 517     /*Add breaks*/
 518
 519     if(SecondaryStructuresStatusMap.size() > 1){
 520         for(std::size_t i {0}, linefactor{1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
 521             if( SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) && SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break) ){
 522                 if(SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]) ){
 523                     SecondaryStructuresStringLine.insert(SecondaryStructuresStringLine.begin() + i + linefactor, secondaryStructureTypeNames[secondaryStructureTypes::Break]);
 524                     ++linefactor;
 525                 }
 526             }
 527         }
 528     }
 529     return SecondaryStructuresStringLine;
 530 }
 531
 532 secondaryStructures::~secondaryStructures(){
 533     SecondaryStructuresStatusMap.resize(0);
 534     SecondaryStructuresStringLine.resize(0);
 535 }
 536
 537 DsspTool::DsspStorage::DsspStorage(){
 538     storaged_data.resize(0);
 539 }
 540
 541 void DsspTool::DsspStorage::clearAll(){
 542     storaged_data.resize(0);
 543 }
 544
 545 std::mutex DsspTool::DsspStorage::mx;
 546
 547 void DsspTool::DsspStorage::storageData(int frnr, std::string data){
 548     std::lock_guard<std::mutex> guardian(mx);
 549     std::pair<int, std::string> datapair(frnr, data);
 550     storaged_data.push_back(datapair);
 551 }
 552
 553 std::vector<std::pair<int, std::string>> DsspTool::DsspStorage::returnData(){
 554     std::sort(storaged_data.begin(), storaged_data.end());
 555     return storaged_data;
 556 }
 557
 558 void alternateNeighborhoodSearch::setCutoff(const real &cutoff_init){
 559     cutoff = cutoff_init;
 560 }
 561
 562 void alternateNeighborhoodSearch::FixAtomCoordinates(real &coordinate, const real vector_length){
 563     while (coordinate < 0) {
 564         coordinate += vector_length;
 565     }
 566     while (coordinate >= vector_length) {
 567         coordinate -= vector_length;
 568     }
 569 }
 570
 571 void alternateNeighborhoodSearch::ReCalculatePBC(int &x, const int &x_max) {
 572     if (x < 0) {
 573         x += x_max;
 574     }
 575     if (x >= x_max) {
 576         x -= x_max;
 577     }
 578 }
 579
 580 void alternateNeighborhoodSearch::GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
 581     rvec coordinates, box_vector_length;
 582     num_of_miniboxes.resize(0);
 583     num_of_miniboxes.resize(3);
 584     for (std::size_t i{XX}; i <= ZZ; ++i) {
 585         box_vector_length[i] = std::sqrt(
 586                     std::pow(fr.box[i][XX], 2) + std::pow(fr.box[i][YY], 2) + std::pow(fr.box[i][ZZ], 2));
 587         num_of_miniboxes[i] = std::floor((box_vector_length[i] / cutoff)) + 1;
 588     }
 589     MiniBoxesMap.resize(0);
 590     MiniBoxesReverseMap.resize(0);
 591     MiniBoxesMap.resize(num_of_miniboxes[XX], std::vector<std::vector<std::vector<std::size_t> > >(
 592                              num_of_miniboxes[YY], std::vector<std::vector<std::size_t> >(
 593                              num_of_miniboxes[ZZ], std::vector<std::size_t>(
 594                              0))));
 595     MiniBoxesReverseMap.resize(IndexMap.size(), std::vector<std::size_t>(3));
 596     for (std::vector<ResInfo>::const_iterator i {IndexMap.begin()}; i != IndexMap.end(); ++i) {
 597         for (std::size_t j{XX}; j <= ZZ; ++j) {
 598             coordinates[j] = fr.x[i->getIndex(backboneAtomTypes::AtomCA)][j];
 599             FixAtomCoordinates(coordinates[j], box_vector_length[j]);
 600         }
 601         MiniBoxesMap[std::floor(coordinates[XX] / cutoff)][std::floor(coordinates[YY] / cutoff)][std::floor(
 602                           coordinates[ZZ] / cutoff)].push_back(i - IndexMap.begin());
 603         for (std::size_t j{XX}; j <= ZZ; ++j){
 604             MiniBoxesReverseMap[i - IndexMap.begin()][j] = std::floor(coordinates[j] / cutoff);
 605         }
 606     }
 607 }
 608
 609 void alternateNeighborhoodSearch::AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
 610     GetMiniBoxesMap(fr, IndexMap);
 611     MiniBoxSize[XX] = MiniBoxesMap.size();
 612     MiniBoxSize[YY] = MiniBoxesMap.front().size();
 613     MiniBoxSize[ZZ] = MiniBoxesMap.front().front().size();
 614     PairMap.resize(0);
 615     PairMap.resize(IndexMap.size(), std::vector<bool>(IndexMap.size(), false));
 616     ResiI = PairMap.begin();
 617     ResiJ = ResiI->begin();
 618
 619     for (std::vector<ResInfo>::const_iterator i = IndexMap.begin(); i != IndexMap.end(); ++i){
 620         for (offset[XX] = -1; offset[XX] <= 1; ++offset[XX]) {
 621             for (offset[YY] = -1; offset[YY] <= 1; ++offset[YY]) {
 622                 for (offset[ZZ] = -1; offset[ZZ] <= 1; ++offset[ZZ]) {
 623                     for (std::size_t k{XX}; k <= ZZ; ++k) {
 624                         fixBox[k] = MiniBoxesReverseMap[i - IndexMap.begin()][k] + offset[k];
 625                         ReCalculatePBC(fixBox[k], MiniBoxSize[k]);
 626                     }
 627                     for (std::size_t j{0}; j < MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]].size(); ++j) {
 628                         if ( (i - IndexMap.begin()) != MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]){
 629                             PairMap[i - IndexMap.begin()][MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]] = true;
 630                             PairMap[MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]][i - IndexMap.begin()] = true;
 631                         }
 632                     }
 633                 }
 634             }
 635         }
 636     }
 637 }
 638
 639 bool alternateNeighborhoodSearch::findNextPair(){
 640
 641     if(!PairMap.size()){
 642         return false;
 643     }
 644
 645     for(; ResiI != PairMap.end(); ++ResiI, ResiJ = ResiI->begin() ){
 646         for(; ResiJ != ResiI->end(); ++ResiJ){
 647             if(*ResiJ){
 648                 resiIpos = ResiI - PairMap.begin();
 649                 resiJpos = ResiJ - ResiI->begin();
 650                 if ( ResiJ != ResiI->end() ){
 651                     ++ResiJ;
 652                 }
 653                 else if (ResiI != PairMap.end()) {
 654                     ++ResiI;
 655                     ResiJ = ResiI->begin();
 656                 }
 657                 else {
 658                     return false; // ???
 659                 }
 660                 return true;
 661             }
 662         }
 663     }
 664
 665     return false;
 666 }
 667
 668 std::size_t alternateNeighborhoodSearch::getResiI() const {
 669     return resiIpos;
 670 }
 671
 672 std::size_t alternateNeighborhoodSearch::getResiJ() const {
 673     return resiJpos;
 674 }
 675
 676
 677 DsspTool::DsspStorage DsspTool::Storage;
 678
 679 DsspTool::DsspTool(){
 680 }
 681
 682 void DsspTool::calculateBends(const t_trxframe &fr, const t_pbc *pbc)
 683 {
 684    const float benddegree{ 70.0 }, maxdist{ 2.5 };
 685    float       degree{ 0 }, vdist{ 0 }, vprod{ 0 };
 686    gmx::RVec   a{ 0, 0, 0 }, b{ 0, 0, 0 };
 687    for (std::size_t i{ 0 }; i + 1 < IndexMap.size(); ++i)
 688    {
 689        if (CalculateAtomicDistances(static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
 690                                     static_cast<int>(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)),
 691                                     fr,
 692                                     pbc)
 693            > maxdist)
 694        {
 695            PatternSearch.SecondaryStructuresStatusMap[i].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i + 1]);
 696            PatternSearch.SecondaryStructuresStatusMap[i + 1].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i]);
 697
 698 //           std::cout << "Break between " << i + 1 << " and " << i + 2 << std::endl;
 699        }
 700    }
 701    for (std::size_t i{ 2 }; i + 2 < IndexMap.size() ; ++i)
 702    {
 703        if (PatternSearch.SecondaryStructuresStatusMap[i - 2].getStatus(secondaryStructureTypes::Break) ||
 704            PatternSearch.SecondaryStructuresStatusMap[i - 1].getStatus(secondaryStructureTypes::Break) ||
 705            PatternSearch.SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) ||
 706            PatternSearch.SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break)
 707           )
 708        {
 709            continue;
 710        }
 711        for (int j{ 0 }; j < 3; ++j)
 712        {
 713            a[j] = fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j]
 714                   - fr.x[IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA)][j];
 715            b[j] = fr.x[IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA)][j]
 716                   - fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j];
 717        }
 718        vdist = (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
 719        vprod = CalculateAtomicDistances(IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA),
 720                                         IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
 721                                         fr,
 722                                         pbc)
 723                * gmx::c_angstrom / gmx::c_nano
 724                * CalculateAtomicDistances(IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
 725                                           IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA),
 726                                           fr,
 727                                           pbc)
 728                * gmx::c_angstrom / gmx::c_nano;
 729        degree = std::acos(vdist / vprod) * gmx::c_rad2Deg;
 730        if (degree > benddegree)
 731        {
 732            PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bend);
 733        }
 734    }
 735 }
 736
 737 void DsspTool::calculateHBondEnergy(ResInfo& Donor,
 738                        ResInfo& Acceptor,
 739                        const t_trxframe&          fr,
 740                        const t_pbc*               pbc)
 741 {
 742    /*
 743     * DSSP uses eq from dssp 2.x
 744     * kCouplingConstant = 27.888,  //  = 332 * 0.42 * 0.2
 745     * E = k * (1/rON + 1/rCH - 1/rOH - 1/rCN) where CO comes from one AA and NH from another
 746     * if R is in A
 747     * Hbond if E < -0.5
 748     *
 749     * For the note, H-Bond Donor is N-H («Donor of H») and H-Bond Acceptor is C=O («Acceptor of H»)
 750     *
 751     */
 752
 753     if (CalculateAtomicDistances(
 754                 Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
 755         >= minimalCAdistance)
 756     {
 757         return void();
 758     }
 759
 760     const float kCouplingConstant = 27.888;
 761     const float minimalAtomDistance{ 0.5 },
 762             minEnergy{ -9.9 };
 763     float HbondEnergy{ 0 };
 764     float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 };
 765
 766 //    std::cout << "For Donor №" << Donor.info->nr - 1 << " and Accpetor №" << Acceptor.info->nr - 1 << std::endl;
 767
 768     if( !(Donor.is_proline) && (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
 769                                 && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || initParams.addHydrogens ) ) ){ // TODO
 770         distanceNO = CalculateAtomicDistances(
 771                Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 772         distanceNC = CalculateAtomicDistances(
 773                Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 774         if (initParams.addHydrogens){
 775             if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
 776                rvec atomH{};
 777                float prevCODist {CalculateAtomicDistances(Donor.prevResi->getIndex(backboneAtomTypes::AtomC), Donor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)};
 778                for (int i{XX}; i <= ZZ; ++i){
 779                    float prevCO = fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
 780                    atomH[i] = fr.x[Donor.getIndex(backboneAtomTypes::AtomH)][i]; // Но на самом деле берутся координаты N
 781                    atomH[i] += prevCO / prevCODist;
 782                }
 783                distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 784                distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 785             }
 786             else{
 787                 distanceHO = distanceNO;
 788                 distanceHC = distanceNC;
 789             }
 790        }
 791        else {
 792            distanceHO = CalculateAtomicDistances(
 793                    Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 794            distanceHC = CalculateAtomicDistances(
 795                    Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 796        }
 797        if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
 798         || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
 799        {
 800             HbondEnergy = minEnergy;
 801        }
 802        else{
 803            HbondEnergy =
 804                    kCouplingConstant
 805                    * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
 806        }
 807
 808 //       std::cout << "CA-CA distance: " << CalculateAtomicDistances(
 809 //                        Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc) << std::endl;
 810 //       std::cout << "N-O distance: " << distanceNO << std::endl;
 811 //       std::cout << "N-C distance: " << distanceNC << std::endl;
 812 //       std::cout << "H-O distance: " << distanceHO << std::endl;
 813 //       std::cout << "H-C distance: " << distanceHC << std::endl;
 814
 815        HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
 816
 817 //       if ( HbondEnergy < minEnergy ){ // I don't think that this is correct
 818 //            HbondEnergy = minEnergy;
 819 //       }
 820
 821 //       std::cout << "Calculated energy = " << HbondEnergy << std::endl;
 822     }
 823 //    else{
 824 //        std::cout << "Donor Is Proline" << std::endl;
 825 //    }
 826
 827     if (HbondEnergy < Donor.acceptorEnergy[0]){
 828            Donor.acceptor[1] = Donor.acceptor[0];
 829            Donor.acceptor[0] = Acceptor.info;
 830            Donor.acceptorEnergy[0] = HbondEnergy;
 831     }
 832     else if (HbondEnergy < Donor.acceptorEnergy[1]){
 833            Donor.acceptor[1] = Acceptor.info;
 834            Donor.acceptorEnergy[1] = HbondEnergy;
 835     }
 836
 837     if (HbondEnergy < Acceptor.donorEnergy[0]){
 838            Acceptor.donor[1] = Acceptor.donor[0];
 839            Acceptor.donor[0] = Donor.info;
 840            Acceptor.donorEnergy[0] = HbondEnergy;
 841     }
 842     else if (HbondEnergy < Acceptor.donorEnergy[1]){
 843            Acceptor.donor[1] = Donor.info;
 844            Acceptor.donorEnergy[1] = HbondEnergy;
 845     }
 846 }
 847
 848
 849 /* Calculate Distance From B to A */
 850 float DsspTool::CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
 851 {
 852    gmx::RVec r{ 0, 0, 0 };
 853    pbc_dx(pbc, fr.x[A], fr.x[B], r.as_vec());
 854    return r.norm() * gmx::c_nm2A; // НЕ ТРОГАТЬ
 855 }
 856
 857 /* Calculate Distance From B to A, where A is only fake coordinates */
 858 float DsspTool::CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
 859 {
 860    gmx::RVec r{ 0, 0, 0 };
 861    pbc_dx(pbc, A, fr.x[B], r.as_vec());
 862    return r.norm() * gmx::c_nm2A; // НЕ ТРОГАТЬ
 863 }
 864
 865 void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz)
 866 {
 867    initParams = initParamz;
 868    ResInfo _backboneAtoms;
 869    std::size_t                 i{ 0 };
 870    std::string proLINE;
 871    int resicompare{ top.atoms()->atom[static_cast<std::size_t>(*(initParams.sel_.atomIndices().begin()))].resind };
 872    IndexMap.resize(0);
 873    IndexMap.push_back(_backboneAtoms);
 874    IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
 875    proLINE = *(IndexMap[i].info->name);
 876    if( proLINE.compare("PRO") == 0 ){
 877        IndexMap[i].is_proline = true;
 878    }
 879
 880    for (gmx::ArrayRef<const int>::iterator ai{ initParams.sel_.atomIndices().begin() }; (ai != initParams.sel_.atomIndices().end()); ++ai){
 881        if (resicompare != top.atoms()->atom[static_cast<std::size_t>(*ai)].resind)
 882        {
 883            ++i;
 884            resicompare = top.atoms()->atom[static_cast<std::size_t>(*ai)].resind;
 885            IndexMap.emplace_back(_backboneAtoms);
 886            IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
 887            proLINE = *(IndexMap[i].info->name);
 888            if( proLINE.compare("PRO") == 0 ){
 889                IndexMap[i].is_proline = true;
 890            }
 891
 892        }
 893        std::string atomname(*(top.atoms()->atomname[static_cast<std::size_t>(*ai)]));
 894        if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA])
 895        {
 896            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomCA)] = *ai;
 897        }
 898        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC])
 899        {
 900            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomC)] = *ai;
 901        }
 902        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO])
 903        {
 904            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomO)] = *ai;
 905        }
 906        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN])
 907        {
 908            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomN)] = *ai;
 909            if (initParamz.addHydrogens == true){
 910                IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
 911            }
 912        }
 913        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH] && initParamz.addHydrogens == false) // Юзать водород в структуре
 914        {
 915            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
 916        }
 917
 918
 919
 920 //       if( atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO]
 921 //       || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH]){
 922 //           std::cout << "Atom " << atomname << " №" << *ai << " From Resi " << *(top.atoms()->resinfo[i].name) << " №" << resicompare << std::endl;
 923 //       }
 924    }
 925
 926    for (std::size_t j {1}; j < IndexMap.size(); ++j){
 927        IndexMap[j].prevResi = &(IndexMap[j - 1]);
 928
 929        IndexMap[j - 1].nextResi = &(IndexMap[j]);
 930
 931 //           std::cout << "Resi " << IndexMap[i].info->nr << *(IndexMap[i].info->name) << std::endl;
 932 //           std::cout << "Prev resi is " << IndexMap[i].prevResi->info->nr << *(IndexMap[i].prevResi->info->name) << std::endl;
 933 //           std::cout << "Prev resi's next resi is " << IndexMap[i - 1].nextResi->info->nr << *(IndexMap[i - 1].nextResi->info->name) << std::endl;
 934 //         std::cout << IndexMap[j].prevResi->info->nr;
 935 //         std::cout << *(IndexMap[j].prevResi->info->name) ;
 936 //         std::cout << " have CA = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomCA) ;
 937 //         std::cout << " C = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomC);
 938 //         std::cout << " O = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomO);
 939 //         std::cout << " N = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomN);
 940 //         std::cout << " H = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomH) << std::endl;
 941    }
 942
 943    nres = i + 1;
 944 }
 945
 946 void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc)
 947 {
 948
 949     switch(initParams.NBS){
 950     case (NBSearchMethod::Classique): {
 951
 952         // store positions of CA atoms to use them for nbSearch
 953         std::vector<gmx::RVec> positionsCA_;
 954         for (std::size_t i{ 0 }; i < IndexMap.size(); ++i)
 955         {
 956             positionsCA_.emplace_back(fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)]);
 957         }
 958
 959         AnalysisNeighborhood nb_;
 960         nb_.setCutoff(initParams.cutoff_);
 961         AnalysisNeighborhoodPositions       nbPos_(positionsCA_);
 962         gmx::AnalysisNeighborhoodSearch     start      = nb_.initSearch(pbc, nbPos_);
 963         gmx::AnalysisNeighborhoodPairSearch pairSearch = start.startPairSearch(nbPos_);
 964         gmx::AnalysisNeighborhoodPair       pair;
 965         while (pairSearch.findNextPair(&pair))
 966         {
 967             if(CalculateAtomicDistances(
 968                         IndexMap[pair.refIndex()].getIndex(backboneAtomTypes::AtomCA), IndexMap[pair.testIndex()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
 969                 < minimalCAdistance){
 970                 calculateHBondEnergy(IndexMap[pair.refIndex()], IndexMap[pair.testIndex()], fr, pbc);
 971                 if (IndexMap[pair.testIndex()].info != IndexMap[pair.refIndex() + 1].info){
 972                     calculateHBondEnergy(IndexMap[pair.testIndex()], IndexMap[pair.refIndex()], fr, pbc);
 973                 }
 974             }
 975         }
 976
 977         break;
 978     }
 979     case (NBSearchMethod::Experimental): { // TODO FIX
 980
 981         alternateNeighborhoodSearch as_;
 982
 983         as_.setCutoff(initParams.cutoff_);
 984
 985         as_.AltPairSearch(fr, IndexMap);
 986
 987         while (as_.findNextPair()){
 988             if(CalculateAtomicDistances(
 989                         IndexMap[as_.getResiI()].getIndex(backboneAtomTypes::AtomCA), IndexMap[as_.getResiJ()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
 990                 < minimalCAdistance){
 991                 calculateHBondEnergy(IndexMap[as_.getResiI()], IndexMap[as_.getResiJ()], fr, pbc);
 992                 if (IndexMap[as_.getResiJ()].info != IndexMap[as_.getResiI() + 1].info){
 993                     calculateHBondEnergy(IndexMap[as_.getResiJ()], IndexMap[as_.getResiI()], fr, pbc);
 994                 }
 995             }
 996         }
 997
 998         break;
 999     }
1000     default: {
1001
1002         for(std::vector<ResInfo>::iterator Donor {IndexMap.begin()}; Donor != IndexMap.end() ; ++Donor){
1003             for(std::vector<ResInfo>::iterator Acceptor {Donor + 1} ; Acceptor != IndexMap.end() ; ++Acceptor){
1004                 if(CalculateAtomicDistances(
1005                             Donor->getIndex(backboneAtomTypes::AtomCA), Acceptor->getIndex(backboneAtomTypes::AtomCA), fr, pbc)
1006                     < minimalCAdistance){
1007                     calculateHBondEnergy(*Donor, *Acceptor, fr, pbc);
1008                     if (Acceptor != Donor + 1){
1009                         calculateHBondEnergy(*Acceptor, *Donor, fr, pbc);
1010                     }
1011                 }
1012             }
1013         }
1014         break;
1015     }
1016     }
1017
1018
1019 //    for(std::size_t i {0}; i < IndexMap.size(); ++i){
1020 //        std::cout << IndexMap[i].info->nr << " " << *(IndexMap[i].info->name) << std::endl;
1021 //    }
1022
1023    PatternSearch.initiateSearch(IndexMap, initParams.PPHelices);
1024    calculateBends(fr, pbc);
1025    Storage.storageData(frnr, PatternSearch.patternSearch());
1026
1027 }
1028
1029 std::vector<std::pair<int, std::string>> DsspTool::getData(){
1030     return Storage.returnData();
1031 }
1032
1033 } // namespace analysismodules
1034
1035 } // namespace gmx