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 of the last residue
  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 }
  98
  99 void secondaryStructures::secondaryStructuresData::setStatus(const HelixPositions helixPosition, const turnsTypes turn){
 100     TurnsStatusArray[static_cast<std::size_t>(turn)] = helixPosition;
 101 }
 102
 103 bool secondaryStructures::secondaryStructuresData::getStatus(const secondaryStructureTypes secondaryStructureTypeName) const{
 104     return SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)];
 105 }
 106
 107 bool secondaryStructures::secondaryStructuresData::isBreakPartnerWith(const secondaryStructuresData *partner) const{
 108     return breakPartners[0] == partner || breakPartners[1] == partner;
 109 }
 110
 111 HelixPositions secondaryStructures::secondaryStructuresData::getStatus(const turnsTypes turn) const{
 112     return TurnsStatusArray[static_cast<std::size_t>(turn)];
 113 }
 114
 115 void secondaryStructures::secondaryStructuresData::setBreak(secondaryStructuresData *breakPartner){
 116     if (breakPartners[0] == nullptr){
 117         breakPartners[0] = breakPartner;
 118     }
 119     else{
 120         breakPartners[1] = breakPartner;
 121     }
 122     setStatus(secondaryStructureTypes::Break);
 123 }
 124
 125 bool secondaryStructures::hasHBondBetween(std::size_t Donor, std::size_t Acceptor) const{ // prob should add resi name comparison ?
 126     if( (*ResInfoMap)[Donor].acceptor[0] == nullptr ||
 127         (*ResInfoMap)[Donor].acceptor[1] == nullptr ||
 128         (*ResInfoMap)[Acceptor].info == nullptr ){
 129         return false;
 130     }
 131     else {
 132
 133 //        std::cout << "Comparing DONOR " << Donor << " And ACCEPTOR " << Acceptor << " :";
 134 //        std::cout << "DONOR's acceptor adresses are " << (*ResInfoMap)[Donor].acceptor[0] << ", " << (*ResInfoMap)[Donor].acceptor[1] << " and ACCEPTOR adress is " << (*ResInfoMap)[Acceptor].info << std::endl;
 135 //        std::cout << "DONOR's acceptors' nr are = " << (*ResInfoMap)[Donor].acceptor[0]->nr << ", " << (*ResInfoMap)[Donor].acceptor[1]->nr << " And ACCEPTOR's nr = " << (*ResInfoMap)[Acceptor].info->nr << std::endl;
 136 //        std::cout << "DONOR's acceptors' chainID are = " << (*ResInfoMap)[Donor].acceptor[0]->chainid << ", " << (*ResInfoMap)[Donor].acceptor[1]->chainid << " And ACCEPTOR's chainID = " << (*ResInfoMap)[Acceptor].info->chainid << std::endl;
 137
 138 //        if( ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
 139 //                ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff ) ){
 140 //            std::cout << "HBond Exist" << std::endl;
 141 //        }
 142         return ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
 143                 ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff );
 144     }
 145 }
 146
 147 bool secondaryStructures::NoChainBreaksBetween(std::size_t Resi1, std::size_t Resi2) const{
 148     bool flag{true};
 149     std::size_t i{Resi1}, j[Resi2]; // From i to j
 150     if ( Resi2 < Resi1){
 151         std::swap(i, j);
 152     }
 153
 154     for (; i != j; ++i){
 155         flag = !(SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]));
 156     }
 157
 158     if (!flag){
 159         std::cout << "Patternsearch has detected a CHAINBREAK between " << Resi1 << " and " << Resi2 << std::endl;
 160     }
 161     return flag;
 162 }
 163
 164 bridgeTypes secondaryStructures::calculateBridge(std::size_t i, std::size_t j) const{
 165     if( i < 1 || j < 1 || i + 1 >= ResInfoMap->size() || j + 1 >= ResInfoMap->size() ){
 166         return bridgeTypes::None;
 167     }
 168     if(NoChainBreaksBetween(i - 1, i + 1) && NoChainBreaksBetween(j - 1, j + 1)){
 169         if((hasHBondBetween(i + 1, j) && hasHBondBetween(j, i - 1)) || (hasHBondBetween(j + 1, i) && hasHBondBetween(i, j - 1)) ){ //possibly swap
 170             return bridgeTypes::ParallelBridge;
 171         }
 172         else if((hasHBondBetween(i + 1, j - 1) && hasHBondBetween(j + 1, i - 1)) || (hasHBondBetween(j, i) && hasHBondBetween(i, j)) ){ //possibly swap
 173             return bridgeTypes::AntiParallelBridge;
 174         }
 175     }
 176     return bridgeTypes::None;
 177 }
 178
 179 void secondaryStructures::analyzeBridgesAndLaddersPatterns(){
 180     for(std::size_t i {1}; i + 4 < SecondaryStructuresStatusMap.size(); ++i){
 181         for(std::size_t j {i + 3}; j + 1 < SecondaryStructuresStatusMap.size(); ++j ){
 182             bridgeTypes type {calculateBridge(i, j)};
 183             if (type == bridgeTypes::None){
 184                 continue;
 185             }
 186             bool found {false};
 187         }
 188     }
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201 //    for (std::size_t i{ 1 }; i < HBondsMap.front().size() - 1; ++i){
 202 //        for (std::size_t j{ 1 }; j < HBondsMap.front().size() - 1; ++j){
 203 //            if (std::abs(static_cast<int>(i) - static_cast<int>(j)) > 2){
 204 //                if ((HBondsMap[i - 1][j] && HBondsMap[j][i + 1])    ||
 205 //                    (HBondsMap[j - 1][i] && HBondsMap[i][j + 1])){
 206 //                    Bridge[i].push_back(j);
 207 //                }
 208 //                if ((HBondsMap[i][j] && HBondsMap[j][i])    ||
 209 //                    (HBondsMap[i - 1][j + 1] && HBondsMap[j - 1][i + 1])){
 210 //                    AntiBridge[i].push_back(j);
 211 //                }
 212 //            }
 213 //        }
 214 //    }
 215 //    for (std::size_t i{ 0 }; i < HBondsMap.front().size(); ++i){
 216 //        if ((!Bridge[i].empty() || !AntiBridge[i].empty())){
 217 //            setStatus(i, secondaryStructureTypes::Bulge);
 218 //        }
 219 //    }
 220 //    for (std::size_t i{ 2 }; i + 2 < HBondsMap.front().size(); ++i){
 221 //        for (std::size_t j { i - 2 }; j <= (i + 2); ++j){
 222 //            if (j == i){
 223 //                continue;
 224 //            }
 225 //            else {
 226 //                for (std::vector<bridgeTypes>::const_iterator bridge {Bridges.begin()}; bridge != Bridges.end(); ++bridge ){
 227 //                    if (!getBridge(*bridge)[i].empty() || !getBridge(*bridge)[j].empty()){
 228 //                        for (std::size_t i_resi{ 0 }; i_resi < getBridge(*bridge)[i].size(); ++i_resi){
 229 //                            for (std::size_t j_resi{ 0 }; j_resi < getBridge(*bridge)[j].size(); ++j_resi){
 230 //                                if (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
 231 //                                        - static_cast<int>(getBridge(*bridge)[j][j_resi]))
 232 //                                        && (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
 233 //                                        - static_cast<int>(getBridge(*bridge)[j][j_resi]))
 234 //                                        < 5)){
 235 //                                    if (j < i){
 236 //                                        for (std::size_t k{ 0 }; k <= i - j; ++k){
 237 //                                            setStatus(i + k, secondaryStructureTypes::Ladder);
 238 //                                        }
 239 //                                    }
 240 //                                    else{
 241 //                                        for (std::size_t k{ 0 }; k <= j - i; ++k){
 242 //                                            setStatus(i + k, secondaryStructureTypes::Ladder);
 243 //                                        }
 244 //                                    }
 245 //                                }
 246 //                            }
 247 //                        }
 248 //                    }
 249 //                }
 250 //            }
 251 //        }
 252 //    }
 253 }
 254
 255 void secondaryStructures::analyzeTurnsAndHelicesPatterns(){
 256     for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
 257         std::cout << "Testing Helix_" << static_cast<std::size_t>(i) + 3 << std::endl;
 258         std::size_t stride {static_cast<std::size_t>(i) + 3};
 259         for(std::size_t j {0}; j + static_cast<std::size_t>(i) < SecondaryStructuresStatusMap.size(); ++j){
 260             std::cout << "Testing " << j << " and " << j + stride << std::endl;
 261             if ( hasHBondBetween(j, j + (static_cast<std::size_t>(i))) && NoChainBreaksBetween(j, j + stride) ){
 262                 std::cout << j << " and " << j + stride << " has hbond!" << std::endl;
 263                 SecondaryStructuresStatusMap[j + static_cast<std::size_t>(i)].setStatus(HelixPositions::End, i);
 264
 265                 for (std::size_t k {1}; k < (static_cast<std::size_t>(i)); ++k){
 266                     if( SecondaryStructuresStatusMap[j + k].getStatus(i) == HelixPositions::None ){
 267                         SecondaryStructuresStatusMap[j + k].setStatus(HelixPositions::Middle, i);
 268                     }
 269
 270                 }
 271
 272                 if( SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::End ){
 273                     SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start_AND_End, i);
 274                 }
 275                 else {
 276                     SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start, i);
 277                 }
 278             }
 279         }
 280     }
 281
 282     for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
 283         std::size_t stride {static_cast<std::size_t>(i) + 3};
 284         for(std::size_t j {1}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
 285             if ( (SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start_AND_End ) &&
 286                  (SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start_AND_End ) ){
 287                 bool empty = true;
 288                 secondaryStructureTypes Helix;
 289                 switch(i){
 290                 case turnsTypes::Turn_3:
 291                     for (std::size_t k {0}; empty && k < stride; ++k){
 292                         empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_3);
 293                     }
 294                     Helix = secondaryStructureTypes::Helix_3;
 295                     break;
 296                 case turnsTypes::Turn_5:
 297                     for (std::size_t k {0}; empty && k < stride; ++k){
 298                         empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_5) || (PiHelixPreference && SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_5));
 299                     }
 300                     Helix = secondaryStructureTypes::Helix_4;
 301                     break;
 302                 default:
 303                     Helix = secondaryStructureTypes::Helix_4;
 304                     break;
 305                 }
 306                 std::cout << j << " is HELIX" << std::endl;
 307                 if ( empty ){
 308                     for(std::size_t k {0}; k < (static_cast<std::size_t>(i)); ++k ){
 309                         SecondaryStructuresStatusMap[j + k].setStatus(Helix);
 310                     }
 311                 }
 312             }
 313         }
 314     }
 315
 316     for(std::size_t i {1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
 317         if (SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Loop)){
 318             bool isTurn = false;
 319             for(const turnsTypes &j : {turnsTypes::Turn_3, turnsTypes::Turn_4, turnsTypes::Turn_5}){
 320                 std::size_t stride {static_cast<std::size_t>(i) + 3};
 321                 for(std::size_t k {1}; k < stride; ++k){
 322                     isTurn = (i >= k) && (SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start || SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start_AND_End) ;
 323                 }
 324             }
 325
 326             if (isTurn){
 327                 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Turn);
 328             }
 329             else if (SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Bend)){
 330                 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bend);
 331             }
 332         }
 333     }
 334
 335 }
 336
 337 void secondaryStructures::analyzePPHelicesPatterns(){}
 338
 339 std::string secondaryStructures::patternSearch(){
 340
 341 //    analyzeBridgesAndLaddersPatterns();
 342     analyzeTurnsAndHelicesPatterns();
 343 //    analyzePPHelicesPatterns();
 344
 345 //    for(std::size_t i {0}; i < ResInfoMap->size(); ++i){
 346 //        std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) << std::endl;
 347 //    }
 348
 349 //    std::cout.precision(5);
 350 //    for(std::size_t i{0}; i < ResInfoMap->size(); ++i, std::cout << std::endl << std::endl){
 351 //        std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) ;
 352 //        if ( (*ResInfoMap)[i].donor[0] != nullptr ){
 353 //            std::cout << " has donor[0] = " << (*ResInfoMap)[i].donor[0]->nr << " " << *((*ResInfoMap)[i].donor[0]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[0] << " and" ;
 354 //        }
 355 //        else {
 356 //            std::cout << " has no donor[0] and" ;
 357 //        }
 358 //        if ( (*ResInfoMap)[i].acceptor[0] != nullptr ){
 359 //            std::cout << " has acceptor[0] = " << (*ResInfoMap)[i].acceptor[0]->nr << " " << *((*ResInfoMap)[i].acceptor[0]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[0] ;
 360 //        }
 361 //        else {
 362 //            std::cout << " has no acceptor[0]" ;
 363 //        }
 364 //        std::cout << std::endl << "Also, " << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name);
 365 //        if ( (*ResInfoMap)[i].donor[1] != nullptr ){
 366 //            std::cout << " has donor[1] = " << (*ResInfoMap)[i].donor[1]->nr << " " << *((*ResInfoMap)[i].donor[1]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[1] << " and" ;
 367 //        }
 368 //        else {
 369 //            std::cout << " has no donor[1] and" ;
 370 //        }
 371 //        if ( (*ResInfoMap)[i].acceptor[1] != nullptr ){
 372 //            std::cout << " has acceptor[1] = " << (*ResInfoMap)[i].acceptor[1]->nr << " " << *((*ResInfoMap)[i].acceptor[1]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[1] ;
 373 //        }
 374 //        else {
 375 //            std::cout << " has no acceptor[1]" ;
 376 //        }
 377 //    }
 378
 379     /*Write Data*/
 380
 381     for(std::size_t i {static_cast<std::size_t>(secondaryStructureTypes::Bend)}; i != static_cast<std::size_t>(secondaryStructureTypes::Count); ++i){
 382         for(std::size_t j {0}; j < SecondaryStructuresStatusMap.size(); ++j){
 383             if (SecondaryStructuresStatusMap[j].getStatus(static_cast<secondaryStructureTypes>(i))){
 384                 SecondaryStructuresStringLine[j] = secondaryStructureTypeNames[i] ;
 385             }
 386         }
 387     }
 388
 389     /*Add breaks*/
 390
 391     if(SecondaryStructuresStatusMap.size() > 1){
 392         for(std::size_t i {0}, linefactor{1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
 393             if( SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) && SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break) ){
 394                 if(SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]) ){
 395                     SecondaryStructuresStringLine.insert(SecondaryStructuresStringLine.begin() + i + linefactor, secondaryStructureTypeNames[secondaryStructureTypes::Break]);
 396                     ++linefactor;
 397                 }
 398             }
 399         }
 400     }
 401     return SecondaryStructuresStringLine;
 402 }
 403
 404 secondaryStructures::~secondaryStructures(){
 405     SecondaryStructuresStatusMap.resize(0);
 406     SecondaryStructuresStringLine.resize(0);
 407 }
 408
 409 DsspTool::DsspStorage::DsspStorage(){
 410     storaged_data.resize(0);
 411 }
 412
 413 void DsspTool::DsspStorage::clearAll(){
 414     storaged_data.resize(0);
 415 }
 416
 417 std::mutex DsspTool::DsspStorage::mx;
 418
 419 void DsspTool::DsspStorage::storageData(int frnr, std::string data){
 420     std::lock_guard<std::mutex> guardian(mx);
 421     std::pair<int, std::string> datapair(frnr, data);
 422     storaged_data.push_back(datapair);
 423 }
 424
 425 std::vector<std::pair<int, std::string>> DsspTool::DsspStorage::returnData(){
 426     std::sort(storaged_data.begin(), storaged_data.end());
 427     return storaged_data;
 428 }
 429
 430 void alternateNeighborhoodSearch::setCutoff(const real &cutoff_init){
 431     cutoff = cutoff_init;
 432 }
 433
 434 void alternateNeighborhoodSearch::FixAtomCoordinates(real &coordinate, const real vector_length){
 435     while (coordinate < 0) {
 436         coordinate += vector_length;
 437     }
 438     while (coordinate >= vector_length) {
 439         coordinate -= vector_length;
 440     }
 441 }
 442
 443 void alternateNeighborhoodSearch::ReCalculatePBC(int &x, const int &x_max) {
 444     if (x < 0) {
 445         x += x_max;
 446     }
 447     if (x >= x_max) {
 448         x -= x_max;
 449     }
 450 }
 451
 452 void alternateNeighborhoodSearch::GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
 453     rvec coordinates, box_vector_length;
 454     num_of_miniboxes.resize(0);
 455     num_of_miniboxes.resize(3);
 456     for (std::size_t i{XX}; i <= ZZ; ++i) {
 457         box_vector_length[i] = std::sqrt(
 458                     std::pow(fr.box[i][XX], 2) + std::pow(fr.box[i][YY], 2) + std::pow(fr.box[i][ZZ], 2));
 459         num_of_miniboxes[i] = std::floor((box_vector_length[i] / cutoff)) + 1;
 460     }
 461     MiniBoxesMap.resize(0);
 462     MiniBoxesReverseMap.resize(0);
 463     MiniBoxesMap.resize(num_of_miniboxes[XX], std::vector<std::vector<std::vector<std::size_t> > >(
 464                              num_of_miniboxes[YY], std::vector<std::vector<std::size_t> >(
 465                              num_of_miniboxes[ZZ], std::vector<std::size_t>(
 466                              0))));
 467     MiniBoxesReverseMap.resize(IndexMap.size(), std::vector<std::size_t>(3));
 468     for (std::vector<ResInfo>::const_iterator i {IndexMap.begin()}; i != IndexMap.end(); ++i) {
 469         for (std::size_t j{XX}; j <= ZZ; ++j) {
 470             coordinates[j] = fr.x[i->getIndex(backboneAtomTypes::AtomCA)][j];
 471             FixAtomCoordinates(coordinates[j], box_vector_length[j]);
 472         }
 473         MiniBoxesMap[std::floor(coordinates[XX] / cutoff)][std::floor(coordinates[YY] / cutoff)][std::floor(
 474                           coordinates[ZZ] / cutoff)].push_back(i - IndexMap.begin());
 475         for (std::size_t j{XX}; j <= ZZ; ++j){
 476             MiniBoxesReverseMap[i - IndexMap.begin()][j] = std::floor(coordinates[j] / cutoff);
 477         }
 478     }
 479 }
 480
 481 void alternateNeighborhoodSearch::AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
 482     GetMiniBoxesMap(fr, IndexMap);
 483     MiniBoxSize[XX] = MiniBoxesMap.size();
 484     MiniBoxSize[YY] = MiniBoxesMap.front().size();
 485     MiniBoxSize[ZZ] = MiniBoxesMap.front().front().size();
 486     PairMap.resize(0);
 487     PairMap.resize(IndexMap.size(), std::vector<bool>(IndexMap.size(), false));
 488     ResiI = PairMap.begin();
 489     ResiJ = ResiI->begin();
 490
 491     for (std::vector<ResInfo>::const_iterator i = IndexMap.begin(); i != IndexMap.end(); ++i){
 492         for (offset[XX] = -1; offset[XX] <= 1; ++offset[XX]) {
 493             for (offset[YY] = -1; offset[YY] <= 1; ++offset[YY]) {
 494                 for (offset[ZZ] = -1; offset[ZZ] <= 1; ++offset[ZZ]) {
 495                     for (std::size_t k{XX}; k <= ZZ; ++k) {
 496                         fixBox[k] = MiniBoxesReverseMap[i - IndexMap.begin()][k] + offset[k];
 497                         ReCalculatePBC(fixBox[k], MiniBoxSize[k]);
 498                     }
 499                     for (std::size_t j{0}; j < MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]].size(); ++j) {
 500                         if ( (i - IndexMap.begin()) != MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]){
 501                             PairMap[i - IndexMap.begin()][MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]] = true;
 502                             PairMap[MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]][i - IndexMap.begin()] = true;
 503                         }
 504                     }
 505                 }
 506             }
 507         }
 508     }
 509 }
 510
 511 bool alternateNeighborhoodSearch::findNextPair(){
 512
 513     if(!PairMap.size()){
 514         return false;
 515     }
 516
 517     for(; ResiI != PairMap.end(); ++ResiI, ResiJ = ResiI->begin() ){
 518         for(; ResiJ != ResiI->end(); ++ResiJ){
 519             if(*ResiJ){
 520                 resiIpos = ResiI - PairMap.begin();
 521                 resiJpos = ResiJ - ResiI->begin();
 522                 if ( ResiJ != ResiI->end() ){
 523                     ++ResiJ;
 524                 }
 525                 else if (ResiI != PairMap.end()) {
 526                     ++ResiI;
 527                     ResiJ = ResiI->begin();
 528                 }
 529                 else {
 530                     return false; // ???
 531                 }
 532                 return true;
 533             }
 534         }
 535     }
 536
 537     return false;
 538 }
 539
 540 std::size_t alternateNeighborhoodSearch::getResiI() const {
 541     return resiIpos;
 542 }
 543
 544 std::size_t alternateNeighborhoodSearch::getResiJ() const {
 545     return resiJpos;
 546 }
 547
 548
 549 DsspTool::DsspStorage DsspTool::Storage;
 550
 551 DsspTool::DsspTool(){
 552 }
 553
 554 void DsspTool::calculateBends(const t_trxframe &fr, const t_pbc *pbc)
 555 {
 556    const float benddegree{ 70.0 }, maxdist{ 2.5 };
 557    float       degree{ 0 }, vdist{ 0 }, vprod{ 0 };
 558    gmx::RVec   a{ 0, 0, 0 }, b{ 0, 0, 0 };
 559    for (std::size_t i{ 0 }; i + 1 < IndexMap.size(); ++i)
 560    {
 561        if (CalculateAtomicDistances(static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
 562                                     static_cast<int>(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)),
 563                                     fr,
 564                                     pbc)
 565            > maxdist)
 566        {
 567            PatternSearch.SecondaryStructuresStatusMap[i].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i + 1]);
 568            PatternSearch.SecondaryStructuresStatusMap[i + 1].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i]);
 569
 570 //           std::cout << "Break between " << i + 1 << " and " << i + 2 << std::endl;
 571        }
 572    }
 573    for (std::size_t i{ 2 }; i + 2 < IndexMap.size() ; ++i)
 574    {
 575        if (PatternSearch.SecondaryStructuresStatusMap[i - 2].getStatus(secondaryStructureTypes::Break) ||
 576            PatternSearch.SecondaryStructuresStatusMap[i - 1].getStatus(secondaryStructureTypes::Break) ||
 577            PatternSearch.SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) ||
 578            PatternSearch.SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break)
 579           )
 580        {
 581            continue;
 582        }
 583        for (int j{ 0 }; j < 3; ++j)
 584        {
 585            a[j] = fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j]
 586                   - fr.x[IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA)][j];
 587            b[j] = fr.x[IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA)][j]
 588                   - fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j];
 589        }
 590        vdist = (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
 591        vprod = CalculateAtomicDistances(IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA),
 592                                         IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
 593                                         fr,
 594                                         pbc)
 595                * gmx::c_angstrom / gmx::c_nano
 596                * CalculateAtomicDistances(IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
 597                                           IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA),
 598                                           fr,
 599                                           pbc)
 600                * gmx::c_angstrom / gmx::c_nano;
 601        degree = std::acos(vdist / vprod) * gmx::c_rad2Deg;
 602        if (degree > benddegree)
 603        {
 604            PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bend);
 605        }
 606    }
 607 }
 608
 609 void DsspTool::calculateHBondEnergy(ResInfo& Donor,
 610                        ResInfo& Acceptor,
 611                        const t_trxframe&          fr,
 612                        const t_pbc*               pbc)
 613 {
 614    /*
 615     * DSSP uses eq from dssp 2.x
 616     * kCouplingConstant = 27.888,  //  = 332 * 0.42 * 0.2
 617     * E = k * (1/rON + 1/rCH - 1/rOH - 1/rCN) where CO comes from one AA and NH from another
 618     * if R is in A
 619     * Hbond if E < -0.5
 620     */
 621
 622     const float kCouplingConstant = 27.888;
 623     const float minimalAtomDistance{ 0.5 },
 624             minEnergy{ -9.9 };
 625     float HbondEnergy{ 0 };
 626     float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 };
 627
 628    if( !(Donor.is_proline) ){
 629        if (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
 630            && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || (initParams.addHydrogens) ) )  // Kinda ew
 631        {
 632
 633            distanceNO = CalculateAtomicDistances(
 634                    Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 635            distanceNC = CalculateAtomicDistances(
 636                    Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 637
 638            if (initParams.addHydrogens){
 639                if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
 640 //                   std::cout << "On donor " << Donor.info->nr << *(Donor.info->name) << std::endl;
 641 //                   std::cout << "Prev donor is " << Donor.prevResi->info->nr << *(Donor.prevResi->info->name) << std::endl;
 642 //                   std::cout << "Prev C index is " << Donor.prevResi->getIndex(backboneAtomTypes::AtomC) << std::endl;
 643 //                   std::cout << "Prev O index is " << Donor.prevResi->getIndex(backboneAtomTypes::AtomO) << std::endl;
 644                    rvec atomH{};
 645                    float prevCODist {CalculateAtomicDistances(Donor.prevResi->getIndex(backboneAtomTypes::AtomC), Donor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)};
 646                    for (int i{XX}; i <= ZZ; ++i){
 647                        float prevCO = fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
 648                        atomH[i] = prevCO / prevCODist;
 649                    }
 650                    distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 651                    distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 652                }
 653                else{
 654                    distanceHO = distanceNO;
 655                    distanceHC = distanceNC;
 656                }
 657            }
 658            else {
 659                distanceHO = CalculateAtomicDistances(
 660                        Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 661                distanceHC = CalculateAtomicDistances(
 662                        Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 663            }
 664
 665            if (CalculateAtomicDistances(
 666                        Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
 667                < minimalCAdistance)
 668            {
 669                if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
 670                    || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
 671                {
 672                    HbondEnergy = minEnergy;
 673
 674 //                   std::cout << "HBOND exists cause of distance" << std::endl;
 675                }
 676                else
 677                {
 678                    HbondEnergy =
 679                            kCouplingConstant
 680                            * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
 681                    HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
 682
 683 //                   std::cout.precision(5);
 684 //                   std::cout << "Calculated ENERGY = " << HbondEnergy << std::endl;
 685
 686 //                   if ( HbondEnergy == 0){
 687 //                       std::cout << "Calculated ENERGY = " << HbondEnergy << " For donor " << Donor.info->nr << " and acceptor " << Acceptor.info->nr << std::endl;
 688 //                   }
 689
 690                    if ( HbondEnergy < minEnergy ){
 691                        HbondEnergy = minEnergy;
 692                    }
 693                }
 694            }
 695        }
 696    }
 697 //   else {
 698 //       std::cerr << "PRO DETECTED! THIS IS RESI № " << Donor.info->nr << std::endl; //IT WORKS JUST FINE
 699 //   }
 700
 701    if (HbondEnergy < Donor.acceptorEnergy[0]){
 702        Donor.acceptor[1] = Donor.acceptor[0];
 703        Donor.acceptor[0] = Acceptor.info;
 704        Donor.acceptorEnergy[0] = HbondEnergy;
 705    }
 706    else if (HbondEnergy < Donor.acceptorEnergy[1]){
 707        Donor.acceptor[1] = Acceptor.info;
 708        Donor.acceptorEnergy[1] = HbondEnergy;
 709    }
 710
 711    if (HbondEnergy < Acceptor.donorEnergy[0]){
 712        Acceptor.donor[1] = Acceptor.donor[0];
 713        Acceptor.donor[0] = Donor.info;
 714        Acceptor.donorEnergy[0] = HbondEnergy;
 715    }
 716    else if (HbondEnergy < Acceptor.donorEnergy[1]){
 717        Acceptor.donor[1] = Donor.info;
 718        Acceptor.donorEnergy[1] = HbondEnergy;
 719    }
 720 }
 721
 722
 723 /* Calculate Distance From B to A */
 724 float DsspTool::CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
 725 {
 726    gmx::RVec r{ 0, 0, 0 };
 727    pbc_dx(pbc, fr.x[A], fr.x[B], r.as_vec());
 728 //   return r.norm(); // TODO * gmx::c_nm2A; if not PDB, i guess????
 729    return r.norm() * gmx::c_nm2A; // TODO * gmx::c_nm2A; if not PDB, i guess????
 730 }
 731
 732 /* Calculate Distance From B to A, where A is only fake coordinates */
 733 float DsspTool::CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
 734 {
 735    gmx::RVec r{ 0, 0, 0 };
 736    pbc_dx(pbc, A, fr.x[B], r.as_vec());
 737 //   return r.norm(); // TODO * gmx::c_nm2A; if not PDB, i guess????
 738    return r.norm() * gmx::c_nm2A; // TODO * gmx::c_nm2A; if not PDB, i guess????
 739 }
 740
 741 void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz)
 742 {
 743
 744    std::cout << "Init started" << std::endl;
 745    initParams = initParamz;
 746    ResInfo _backboneAtoms;
 747    std::size_t                 i{ 0 };
 748    std::string proLINE;
 749    int resicompare{ top.atoms()->atom[static_cast<std::size_t>(*(initParams.sel_.atomIndices().begin()))].resind };
 750    IndexMap.resize(0);
 751    IndexMap.push_back(_backboneAtoms);
 752    IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
 753    proLINE = *(IndexMap[i].info->name);
 754    if( proLINE.compare("PRO") == 0 ){
 755        IndexMap[i].is_proline = true;
 756    }
 757
 758    for (gmx::ArrayRef<const int>::iterator ai{ initParams.sel_.atomIndices().begin() }; (ai != initParams.sel_.atomIndices().end()); ++ai){
 759        if (resicompare != top.atoms()->atom[static_cast<std::size_t>(*ai)].resind)
 760        {
 761            ++i;
 762            resicompare = top.atoms()->atom[static_cast<std::size_t>(*ai)].resind;
 763            IndexMap.emplace_back(_backboneAtoms);
 764            IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
 765            proLINE = *(IndexMap[i].info->name);
 766            if( proLINE.compare("PRO") == 0 ){
 767                IndexMap[i].is_proline = true;
 768            }
 769
 770        }
 771        std::string atomname(*(top.atoms()->atomname[static_cast<std::size_t>(*ai)]));
 772        if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA])
 773        {
 774            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomCA)] = *ai;
 775        }
 776        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC])
 777        {
 778            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomC)] = *ai;
 779        }
 780        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO])
 781        {
 782            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomO)] = *ai;
 783        }
 784        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN])
 785        {
 786            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomN)] = *ai;
 787        }
 788        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH] && initParamz.addHydrogens == false) // Юзать водород в структуре
 789        {
 790            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
 791        }
 792
 793 //       if( atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO]
 794 //       || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH]){
 795 //           std::cout << "Atom " << atomname << " №" << *ai << " From Resi " << *(top.atoms()->resinfo[i].name) << " №" << resicompare << std::endl;
 796 //       }
 797    }
 798
 799    for (std::size_t j {1}; j < IndexMap.size(); ++j){
 800        IndexMap[j].prevResi = &(IndexMap[j - 1]);
 801
 802        IndexMap[j - 1].nextResi = &(IndexMap[j]);
 803
 804 //           std::cout << "Resi " << IndexMap[i].info->nr << *(IndexMap[i].info->name) << std::endl;
 805 //           std::cout << "Prev resi is " << IndexMap[i].prevResi->info->nr << *(IndexMap[i].prevResi->info->name) << std::endl;
 806 //           std::cout << "Prev resi's next resi is " << IndexMap[i - 1].nextResi->info->nr << *(IndexMap[i - 1].nextResi->info->name) << std::endl;
 807 //         std::cout << IndexMap[j].prevResi->info->nr;
 808 //         std::cout << *(IndexMap[j].prevResi->info->name) ;
 809 //         std::cout << " have CA = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomCA) ;
 810 //         std::cout << " C = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomC);
 811 //         std::cout << " O = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomO);
 812 //         std::cout << " N = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomN);
 813 //         std::cout << " H = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomH) << std::endl;
 814    }
 815
 816    nres = i + 1;
 817 }
 818
 819 void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc)
 820 {
 821
 822     switch(initParams.NBS){
 823     case (NBSearchMethod::Classique): {
 824
 825         // store positions of CA atoms to use them for nbSearch
 826         std::vector<gmx::RVec> positionsCA_;
 827         for (std::size_t i{ 0 }; i < IndexMap.size(); ++i)
 828         {
 829             positionsCA_.emplace_back(fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)]);
 830         }
 831
 832         AnalysisNeighborhood nb_;
 833         nb_.setCutoff(initParams.cutoff_);
 834         AnalysisNeighborhoodPositions       nbPos_(positionsCA_);
 835         gmx::AnalysisNeighborhoodSearch     start      = nb_.initSearch(pbc, nbPos_);
 836         gmx::AnalysisNeighborhoodPairSearch pairSearch = start.startPairSearch(nbPos_);
 837         gmx::AnalysisNeighborhoodPair       pair;
 838         while (pairSearch.findNextPair(&pair))
 839         {
 840             if(CalculateAtomicDistances(
 841                         IndexMap[pair.refIndex()].getIndex(backboneAtomTypes::AtomCA), IndexMap[pair.testIndex()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
 842                 < minimalCAdistance){
 843                 calculateHBondEnergy(IndexMap[pair.refIndex()], IndexMap[pair.testIndex()], fr, pbc);
 844                 if (IndexMap[pair.testIndex()].info != IndexMap[pair.refIndex() + 1].info){
 845                     calculateHBondEnergy(IndexMap[pair.testIndex()], IndexMap[pair.refIndex()], fr, pbc);
 846                 }
 847             }
 848         }
 849
 850         break;
 851     }
 852     case (NBSearchMethod::Experimental): { // TODO FIX
 853
 854         alternateNeighborhoodSearch as_;
 855
 856         as_.setCutoff(initParams.cutoff_);
 857
 858         as_.AltPairSearch(fr, IndexMap);
 859
 860         while (as_.findNextPair()){
 861             if(CalculateAtomicDistances(
 862                         IndexMap[as_.getResiI()].getIndex(backboneAtomTypes::AtomCA), IndexMap[as_.getResiJ()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
 863                 < minimalCAdistance){
 864                 calculateHBondEnergy(IndexMap[as_.getResiI()], IndexMap[as_.getResiJ()], fr, pbc);
 865                 if (IndexMap[as_.getResiJ()].info != IndexMap[as_.getResiI() + 1].info){
 866                     calculateHBondEnergy(IndexMap[as_.getResiJ()], IndexMap[as_.getResiI()], fr, pbc);
 867                 }
 868             }
 869         }
 870
 871         break;
 872     }
 873     default: {
 874
 875         for(std::vector<ResInfo>::iterator Donor {IndexMap.begin()}; Donor != IndexMap.end() ; ++Donor){
 876             for(std::vector<ResInfo>::iterator Acceptor {Donor + 1} ; Acceptor != IndexMap.end() ; ++Acceptor){
 877                 if(CalculateAtomicDistances(
 878                             Donor->getIndex(backboneAtomTypes::AtomCA), Acceptor->getIndex(backboneAtomTypes::AtomCA), fr, pbc)
 879                     < minimalCAdistance){
 880                     calculateHBondEnergy(*Donor, *Acceptor, fr, pbc);
 881                     if (Acceptor != Donor + 1){
 882                         calculateHBondEnergy(*Acceptor, *Donor, fr, pbc);
 883                     }
 884                 }
 885             };
 886         };
 887         break;
 888     }
 889     }
 890
 891
 892 //    for(std::size_t i {0}; i < IndexMap.size(); ++i){
 893 //        std::cout << IndexMap[i].info->nr << " " << *(IndexMap[i].info->name) << std::endl;
 894 //    }
 895
 896    PatternSearch.initiateSearch(IndexMap, initParams.PPHelices);
 897    calculateBends(fr, pbc);
 898    Storage.storageData(frnr, PatternSearch.patternSearch());
 899
 900 }
 901
 902 std::vector<std::pair<int, std::string>> DsspTool::getData(){
 903     return Storage.returnData();
 904 }
 905
 906 } // namespace analysismodules
 907
 908 } // namespace gmx