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