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, const int _pp_stretch){
  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     pp_stretch = _pp_stretch;
  92     SecondaryStructuresStatusMap.resize(ResInfoMatrix.size());
  93     SecondaryStructuresStringLine.resize(ResInfoMatrix.size(), '~');
  94 }
  95
  96 void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName, bool status){
  97     SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)] = status;
  98     if (status){
  99         SecondaryStructuresStatus = secondaryStructureTypeName;
 100     }
 101     else{
 102         SecondaryStructuresStatus = secondaryStructureTypes::Loop;
 103     }
 104 }
 105
 106 void secondaryStructures::secondaryStructuresData::setStatus(const HelixPositions helixPosition, const turnsTypes turn){
 107     TurnsStatusArray[static_cast<std::size_t>(turn)] = helixPosition;
 108 }
 109
 110 bool secondaryStructures::secondaryStructuresData::getStatus(const secondaryStructureTypes secondaryStructureTypeName) const{
 111     return SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)];
 112 }
 113
 114 bool secondaryStructures::secondaryStructuresData::isBreakPartnerWith(const secondaryStructuresData *partner) const{
 115     return breakPartners[0] == partner || breakPartners[1] == partner;
 116 }
 117
 118 HelixPositions secondaryStructures::secondaryStructuresData::getStatus(const turnsTypes turn) const{
 119     return TurnsStatusArray[static_cast<std::size_t>(turn)];
 120 }
 121
 122 secondaryStructureTypes secondaryStructures::secondaryStructuresData::getStatus() const{
 123     return  SecondaryStructuresStatus;
 124 }
 125
 126 void secondaryStructures::secondaryStructuresData::setBreak(secondaryStructuresData *breakPartner){
 127     if (breakPartners[0] == nullptr){
 128         breakPartners[0] = breakPartner;
 129     }
 130     else{
 131         breakPartners[1] = breakPartner;
 132     }
 133     setStatus(secondaryStructureTypes::Break);
 134 }
 135
 136 void secondaryStructures::secondaryStructuresData::setBridge(secondaryStructuresData *bridgePartner, std::size_t bridgePartnerIndex, bridgeTypes bridgeType){
 137     if(bridgeType == bridgeTypes::ParallelBridge){
 138         bridgePartners[0] = bridgePartner;
 139         bridgePartnersIndexes[0] = bridgePartnerIndex;
 140     }
 141     else if (bridgeType == bridgeTypes::AntiParallelBridge){
 142         bridgePartners[1] = bridgePartner;
 143         bridgePartnersIndexes[1] = bridgePartnerIndex;
 144     }
 145 }
 146
 147 bool secondaryStructures::secondaryStructuresData::hasBridges() const{
 148     return bridgePartners[0] || bridgePartners[1];
 149
 150 }
 151
 152 bool secondaryStructures::secondaryStructuresData::hasBridges(bridgeTypes bridgeType) const{
 153     if(bridgeType == bridgeTypes::ParallelBridge){
 154         return bridgePartners[0] != nullptr;
 155     }
 156     else if (bridgeType == bridgeTypes::AntiParallelBridge){
 157         return bridgePartners[1] != nullptr;
 158     }
 159
 160     return false;
 161
 162 }
 163
 164 bool secondaryStructures::secondaryStructuresData::isBridgePartnerWith(secondaryStructuresData *bridgePartner, bridgeTypes bridgeType) const{
 165     if(bridgeType == bridgeTypes::ParallelBridge){
 166         return bridgePartners[0] == bridgePartner;
 167     }
 168     else if (bridgeType == bridgeTypes::AntiParallelBridge){
 169         return bridgePartners[1] == bridgePartner;
 170     }
 171     return false;
 172 }
 173
 174 std::size_t secondaryStructures::secondaryStructuresData::getBridgePartnerIndex(bridgeTypes bridgeType) const{
 175     if (bridgeType == bridgeTypes::ParallelBridge){
 176         return bridgePartnersIndexes[0];
 177     }
 178     return bridgePartnersIndexes[1];
 179 }
 180
 181 secondaryStructures::secondaryStructuresData secondaryStructures::secondaryStructuresData::getBridgePartner(bridgeTypes bridgeType) const{
 182     if(bridgeType == bridgeTypes::ParallelBridge){
 183         return *(bridgePartners[0]);
 184     }
 185     return *(bridgePartners[1]);
 186 }
 187
 188 bool secondaryStructures::hasHBondBetween(std::size_t Donor, std::size_t Acceptor) const{
 189     if( (*ResInfoMap)[Donor].acceptor[0] == nullptr ||
 190         (*ResInfoMap)[Donor].acceptor[1] == nullptr ||
 191         (*ResInfoMap)[Acceptor].info == nullptr ){
 192 //        std::cout << "Bad hbond check. Reason(s): " ;
 193 //        if ( (*ResInfoMap)[Donor].acceptor[0] == nullptr ){
 194 //            std::cout << "Donor has no acceptor[0]; ";
 195 //        }
 196 //        if ( (*ResInfoMap)[Donor].acceptor[1] == nullptr ){
 197 //            std::cout << "Donor has no acceptor[1]; ";
 198 //        }
 199 //        if ( (*ResInfoMap)[Acceptor].info == nullptr ){
 200 //            std::cout << "No info about acceptor; ";
 201 //        }
 202 //        std::cout << std::endl;
 203         return false;
 204     }
 205 //    else if (!( (*ResInfoMap)[Acceptor].donor[0] == nullptr ||
 206 //              (*ResInfoMap)[Acceptor].donor[1] == nullptr ||
 207 //              (*ResInfoMap)[Donor].info == nullptr )) {
 208 //        std::cout << "Comparing DONOR №" << (*ResInfoMap)[Donor].info->nr << " And ACCEPTOR №" << (*ResInfoMap)[Acceptor].info->nr << ": " << std::endl;
 209 //        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;
 210 //        std::cout << "Donor's acceptors' energy are = " << (*ResInfoMap)[Donor].acceptorEnergy[0] << ", " << (*ResInfoMap)[Donor].acceptorEnergy[1] << std::endl;
 211 //        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;
 212 //        std::cout << "Acceptors's donors' energy are = " << (*ResInfoMap)[Acceptor].donorEnergy[0] << ", " << (*ResInfoMap)[Acceptor].donorEnergy[1] << std::endl;
 213 //        if( ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
 214 //                ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff ) ){
 215 //            std::cout << "HBond Exist" << std::endl;
 216 //        }
 217 //    }
 218 //    else {
 219 //        std::cout << "Bad hbond check. Reason(s): " ;
 220 //        if ( (*ResInfoMap)[Acceptor].donor[0] == nullptr ){
 221 //            std::cout << "Acceptor has no donor[0]; ";
 222 //        }
 223 //        if ( (*ResInfoMap)[Acceptor].donor[1] == nullptr ){
 224 //            std::cout << "Acceptor has no donor[1]; ";
 225 //        }
 226 //        if ( (*ResInfoMap)[Donor].info == nullptr ){
 227 //            std::cout << "No info about donor; ";
 228 //        }
 229 //        std::cout << std::endl;
 230 //    }
 231
 232     return ( (*ResInfoMap)[Donor].acceptor[0] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[0] < HBondEnergyCutOff ) ||
 233            ( (*ResInfoMap)[Donor].acceptor[1] == (*ResInfoMap)[Acceptor].info && (*ResInfoMap)[Donor].acceptorEnergy[1] < HBondEnergyCutOff );
 234
 235
 236 }
 237
 238 bool secondaryStructures::NoChainBreaksBetween(std::size_t Resi1, std::size_t Resi2) const{
 239     std::size_t i{Resi1}, j{Resi2}; // From i to j → i <= j
 240     if ( i > j ){
 241         std::swap(i, j);
 242     }
 243
 244     for (; i != j; ++i){
 245         if ( SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]) ){
 246 //            std::cout << "Patternsearch has detected a CHAINBREAK between " << Resi1 << " and " << Resi2 << std::endl;
 247             return false;
 248         }
 249     }
 250     return true;
 251 }
 252
 253 bridgeTypes secondaryStructures::calculateBridge(std::size_t i, std::size_t j) const{
 254     if( i < 1 || j < 1 || i + 1 >= ResInfoMap->size() || j + 1 >= ResInfoMap->size() ){ // Protection from idiotz, не обязательно нужен
 255         return bridgeTypes::None;
 256     }
 257
 258     ResInfo a{(*ResInfoMap)[i]}, b{(*ResInfoMap)[j]};
 259
 260     if(NoChainBreaksBetween(i - 1, i + 1) && NoChainBreaksBetween(j - 1, j + 1) && a.prevResi && a.nextResi && b.prevResi && b.nextResi){
 261         if((hasHBondBetween(i + 1, j) && hasHBondBetween(j, i - 1)) || (hasHBondBetween(j + 1, i) && hasHBondBetween(i, j - 1)) ){
 262             return bridgeTypes::ParallelBridge;
 263         }
 264         else if((hasHBondBetween(i + 1, j - 1) && hasHBondBetween(j + 1, i - 1)) || (hasHBondBetween(j, i) && hasHBondBetween(i, j)) ){
 265             return bridgeTypes::AntiParallelBridge;
 266         }
 267     }
 268     return bridgeTypes::None;
 269 }
 270
 271 void secondaryStructures::analyzeBridgesAndStrandsPatterns(){
 272 //  Calculate Bridgez
 273     for(std::size_t i {1}; i + 4 < SecondaryStructuresStatusMap.size(); ++i){
 274         for(std::size_t j {i + 3}; j + 1 < SecondaryStructuresStatusMap.size(); ++j ){
 275             switch(calculateBridge(i, j)){
 276             case bridgeTypes::ParallelBridge : {
 277                 SecondaryStructuresStatusMap[i].setBridge(&(SecondaryStructuresStatusMap[j]), j, bridgeTypes::ParallelBridge);
 278                 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bridge);
 279                 SecondaryStructuresStatusMap[j].setBridge(&(SecondaryStructuresStatusMap[i]), i, bridgeTypes::ParallelBridge);
 280                 SecondaryStructuresStatusMap[j].setStatus(secondaryStructureTypes::Bridge);
 281             }
 282             case bridgeTypes::AntiParallelBridge : {
 283                 SecondaryStructuresStatusMap[i].setBridge(&(SecondaryStructuresStatusMap[j]), j, bridgeTypes::AntiParallelBridge);
 284                 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bridge);
 285                 SecondaryStructuresStatusMap[j].setBridge(&(SecondaryStructuresStatusMap[i]), i, bridgeTypes::AntiParallelBridge);
 286                 SecondaryStructuresStatusMap[j].setStatus(secondaryStructureTypes::Bridge);
 287             }
 288             default :
 289                 continue;
 290             }
 291         }
 292     }
 293 //  Calculate Extended Strandz
 294     for(std::size_t i {1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
 295         bool is_Estrand {false};
 296         for(std::size_t j {2}; j != 0; --j){
 297             for(const bridgeTypes &bridgeType : {bridgeTypes::ParallelBridge, bridgeTypes::AntiParallelBridge}){
 298                 if (SecondaryStructuresStatusMap[i].hasBridges(bridgeType) && SecondaryStructuresStatusMap[i + j].hasBridges(bridgeType) ){
 299                     std::size_t i_partner{SecondaryStructuresStatusMap[i].getBridgePartnerIndex(bridgeType)}, j_partner{SecondaryStructuresStatusMap[i + j].getBridgePartnerIndex(bridgeType)}, second_strand{};
 300                     if ( abs(i_partner - j_partner) < 6){
 301                         if (i_partner < j_partner){
 302                             second_strand = i_partner;
 303                         }
 304                         else{
 305                             second_strand = j_partner;
 306                         }
 307                         for(int k{abs(i_partner - j_partner)}; k >= 0; --k){
 308                             if (SecondaryStructuresStatusMap[second_strand + k].getStatus(secondaryStructureTypes::Bridge)){
 309                                 SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Bridge, false);
 310                             }
 311
 312                             SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Strand);
 313                         }
 314                         is_Estrand = true;
 315                     }
 316                 }
 317             }
 318             if (is_Estrand){
 319                 for(std::size_t k{0}; k <= j; ++k){
 320                     if (SecondaryStructuresStatusMap[i + k].getStatus(secondaryStructureTypes::Bridge)){
 321                         SecondaryStructuresStatusMap[i + k].setStatus(secondaryStructureTypes::Bridge, false);
 322                     }
 323                     SecondaryStructuresStatusMap[i + k].setStatus(secondaryStructureTypes::Strand);
 324                 }
 325                 break;
 326             }
 327         }
 328     }
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341 //    for (std::size_t i{ 1 }; i < HBondsMap.front().size() - 1; ++i){
 342 //        for (std::size_t j{ 1 }; j < HBondsMap.front().size() - 1; ++j){
 343 //            if (std::abs(static_cast<int>(i) - static_cast<int>(j)) > 2){
 344 //                if ((HBondsMap[i - 1][j] && HBondsMap[j][i + 1])    ||
 345 //                    (HBondsMap[j - 1][i] && HBondsMap[i][j + 1])){
 346 //                    Bridge[i].push_back(j);
 347 //                }
 348 //                if ((HBondsMap[i][j] && HBondsMap[j][i])    ||
 349 //                    (HBondsMap[i - 1][j + 1] && HBondsMap[j - 1][i + 1])){
 350 //                    AntiBridge[i].push_back(j);
 351 //                }
 352 //            }
 353 //        }
 354 //    }
 355 //    for (std::size_t i{ 0 }; i < HBondsMap.front().size(); ++i){
 356 //        if ((!Bridge[i].empty() || !AntiBridge[i].empty())){
 357 //            setStatus(i, secondaryStructureTypes::Bulge);
 358 //        }
 359 //    }
 360 //    for (std::size_t i{ 2 }; i + 2 < HBondsMap.front().size(); ++i){
 361 //        for (std::size_t j { i - 2 }; j <= (i + 2); ++j){
 362 //            if (j == i){
 363 //                continue;
 364 //            }
 365 //            else {
 366 //                for (std::vector<bridgeTypes>::const_iterator bridge {Bridges.begin()}; bridge != Bridges.end(); ++bridge ){
 367 //                    if (!getBridge(*bridge)[i].empty() || !getBridge(*bridge)[j].empty()){
 368 //                        for (std::size_t i_resi{ 0 }; i_resi < getBridge(*bridge)[i].size(); ++i_resi){
 369 //                            for (std::size_t j_resi{ 0 }; j_resi < getBridge(*bridge)[j].size(); ++j_resi){
 370 //                                if (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
 371 //                                        - static_cast<int>(getBridge(*bridge)[j][j_resi]))
 372 //                                        && (abs(static_cast<int>(getBridge(*bridge)[i][i_resi])
 373 //                                        - static_cast<int>(getBridge(*bridge)[j][j_resi]))
 374 //                                        < 5)){
 375 //                                    if (j < i){
 376 //                                        for (std::size_t k{ 0 }; k <= i - j; ++k){
 377 //                                            setStatus(i + k, secondaryStructureTypes::Ladder);
 378 //                                        }
 379 //                                    }
 380 //                                    else{
 381 //                                        for (std::size_t k{ 0 }; k <= j - i; ++k){
 382 //                                            setStatus(i + k, secondaryStructureTypes::Ladder);
 383 //                                        }
 384 //                                    }
 385 //                                }
 386 //                            }
 387 //                        }
 388 //                    }
 389 //                }
 390 //            }
 391 //        }
 392 //    }
 393 }
 394
 395 void secondaryStructures::analyzeTurnsAndHelicesPatterns(){
 396     for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
 397         std::size_t stride {static_cast<std::size_t>(i) + 3};
 398         for(std::size_t j {0}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
 399             if ( hasHBondBetween(j + stride, j) && NoChainBreaksBetween(j, j + stride) ){
 400                 std::cout << "Resi " << j << " is Helix_" << stride << "start" << std::endl;
 401                 SecondaryStructuresStatusMap[j + stride].setStatus(HelixPositions::End, i);
 402
 403                 for (std::size_t k {1}; k < stride; ++k){
 404                     if( SecondaryStructuresStatusMap[j + k].getStatus(i) == HelixPositions::None ){
 405                         SecondaryStructuresStatusMap[j + k].setStatus(HelixPositions::Middle, i);
 406 //                        SecondaryStructuresStatusMap[j + k].setStatus(secondaryStructureTypes::Turn);
 407                     }
 408                 }
 409
 410                 if( SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::End ){
 411                     SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start_AND_End, i);
 412                 }
 413                 else {
 414                     SecondaryStructuresStatusMap[j].setStatus(HelixPositions::Start, i);
 415                 }
 416             }
 417         }
 418     }
 419
 420     for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
 421         std::size_t stride {static_cast<std::size_t>(i) + 3};
 422         for(std::size_t j {1}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
 423             if ( (SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j - 1].getStatus(i) == HelixPositions::Start_AND_End ) &&
 424                  (SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start || SecondaryStructuresStatusMap[j].getStatus(i) == HelixPositions::Start_AND_End ) ){
 425                 bool empty = true;
 426                 secondaryStructureTypes Helix;
 427                 switch(i){
 428                 case turnsTypes::Turn_3:
 429                     for (std::size_t k {0}; empty && k < stride; ++k){
 430                         empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_3);
 431                     }
 432                     Helix = secondaryStructureTypes::Helix_3;
 433                     break;
 434                 case turnsTypes::Turn_5:
 435                     for (std::size_t k {0}; empty && k < stride; ++k){
 436                         empty = SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Loop ) || SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_5) || (PiHelixPreference && SecondaryStructuresStatusMap[j + k].getStatus(secondaryStructureTypes::Helix_4)); //TODO
 437                     }
 438                     Helix = secondaryStructureTypes::Helix_5;
 439                     break;
 440                 default:
 441                     Helix = secondaryStructureTypes::Helix_4;
 442                     break;
 443                 }
 444                 if ( empty || Helix == secondaryStructureTypes::Helix_4 ){
 445                     for(std::size_t k {0}; k < stride; ++k ){
 446                         std::cout << "Resi " << j << " is Helix_" << static_cast<std::size_t>(Helix) << std::endl;
 447                         SecondaryStructuresStatusMap[j + k].setStatus(Helix);
 448                     }
 449                 }
 450             }
 451         }
 452     }
 453
 454     /* Не знаю зач они в дссп так сделали, этож полное говно */
 455
 456     for(std::size_t i {1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
 457         if (static_cast<int>(SecondaryStructuresStatusMap[i].getStatus()) <= static_cast<int>(secondaryStructureTypes::Turn)){
 458             bool isTurn = false;
 459             for(const turnsTypes &j : {turnsTypes::Turn_3, turnsTypes::Turn_4, turnsTypes::Turn_5}){
 460                 std::size_t stride {static_cast<std::size_t>(j) + 3};
 461                 for(std::size_t k {1}; k < stride and !isTurn; ++k){
 462                     isTurn = (i >= k) && (SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start || SecondaryStructuresStatusMap[i - k].getStatus(j) == HelixPositions::Start_AND_End) ;
 463                 }
 464             }
 465             if (isTurn){
 466                 std::cout << "Resi " << i << " is Turn" << std::endl;
 467                 SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Turn);
 468             }
 469         }
 470     }
 471
 472 }
 473
 474 std::string secondaryStructures::patternSearch(){
 475
 476
 477     analyzeBridgesAndStrandsPatterns();
 478     analyzeTurnsAndHelicesPatterns();
 479
 480 //    for(std::size_t i {0}; i < ResInfoMap->size(); ++i){
 481 //        std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) << std::endl;
 482 //    }
 483
 484 //    std::cout.precision(5);
 485 //    for(std::size_t i{0}; i < ResInfoMap->size(); ++i, std::cout << std::endl << std::endl){
 486 //        std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) ;
 487 //        if ( (*ResInfoMap)[i].donor[0] != nullptr ){
 488 //            std::cout << " has donor[0] = " << (*ResInfoMap)[i].donor[0]->nr << " " << *((*ResInfoMap)[i].donor[0]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[0] << " and" ;
 489 //        }
 490 //        else {
 491 //            std::cout << " has no donor[0] and" ;
 492 //        }
 493 //        if ( (*ResInfoMap)[i].acceptor[0] != nullptr ){
 494 //            std::cout << " has acceptor[0] = " << (*ResInfoMap)[i].acceptor[0]->nr << " " << *((*ResInfoMap)[i].acceptor[0]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[0] ;
 495 //        }
 496 //        else {
 497 //            std::cout << " has no acceptor[0]" ;
 498 //        }
 499 //        std::cout << std::endl << "Also, " << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name);
 500 //        if ( (*ResInfoMap)[i].donor[1] != nullptr ){
 501 //            std::cout << " has donor[1] = " << (*ResInfoMap)[i].donor[1]->nr << " " << *((*ResInfoMap)[i].donor[1]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[1] << " and" ;
 502 //        }
 503 //        else {
 504 //            std::cout << " has no donor[1] and" ;
 505 //        }
 506 //        if ( (*ResInfoMap)[i].acceptor[1] != nullptr ){
 507 //            std::cout << " has acceptor[1] = " << (*ResInfoMap)[i].acceptor[1]->nr << " " << *((*ResInfoMap)[i].acceptor[1]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[1] ;
 508 //        }
 509 //        else {
 510 //            std::cout << " has no acceptor[1]" ;
 511 //        }
 512 //    }
 513
 514     /*Write Data*/
 515
 516     for(std::size_t i {static_cast<std::size_t>(secondaryStructureTypes::Bend)}; i != static_cast<std::size_t>(secondaryStructureTypes::Count); ++i){
 517         for(std::size_t j {0}; j < SecondaryStructuresStatusMap.size(); ++j){
 518             if (SecondaryStructuresStatusMap[j].getStatus(static_cast<secondaryStructureTypes>(i))){
 519                 SecondaryStructuresStringLine[j] = secondaryStructureTypeNames[i] ;
 520             }
 521         }
 522     }
 523
 524     /*Add breaks*/
 525
 526     if(SecondaryStructuresStatusMap.size() > 1){
 527         for(std::size_t i {0}, linefactor{1}; i + 1 < SecondaryStructuresStatusMap.size(); ++i){
 528             if( SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) && SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break) ){
 529                 if(SecondaryStructuresStatusMap[i].isBreakPartnerWith(&SecondaryStructuresStatusMap[i + 1]) && SecondaryStructuresStatusMap[i + 1].isBreakPartnerWith(&SecondaryStructuresStatusMap[i]) ){
 530                     SecondaryStructuresStringLine.insert(SecondaryStructuresStringLine.begin() + i + linefactor, secondaryStructureTypeNames[secondaryStructureTypes::Break]);
 531                     ++linefactor;
 532                 }
 533             }
 534         }
 535     }
 536     return SecondaryStructuresStringLine;
 537 }
 538
 539 secondaryStructures::~secondaryStructures(){
 540     SecondaryStructuresStatusMap.resize(0);
 541     SecondaryStructuresStringLine.resize(0);
 542 }
 543
 544 DsspTool::DsspStorage::DsspStorage(){
 545     storaged_data.resize(0);
 546 }
 547
 548 void DsspTool::DsspStorage::clearAll(){
 549     storaged_data.resize(0);
 550 }
 551
 552 std::mutex DsspTool::DsspStorage::mx;
 553
 554 void DsspTool::DsspStorage::storageData(int frnr, std::string data){
 555     std::lock_guard<std::mutex> guardian(mx);
 556     std::pair<int, std::string> datapair(frnr, data);
 557     storaged_data.push_back(datapair);
 558 }
 559
 560 std::vector<std::pair<int, std::string>> DsspTool::DsspStorage::returnData(){
 561     std::sort(storaged_data.begin(), storaged_data.end());
 562     return storaged_data;
 563 }
 564
 565 void alternateNeighborhoodSearch::setCutoff(const real &cutoff_init){
 566     cutoff = cutoff_init;
 567 }
 568
 569 void alternateNeighborhoodSearch::FixAtomCoordinates(real &coordinate, const real vector_length){
 570     while (coordinate < 0) {
 571         coordinate += vector_length;
 572     }
 573     while (coordinate >= vector_length) {
 574         coordinate -= vector_length;
 575     }
 576 }
 577
 578 void alternateNeighborhoodSearch::ReCalculatePBC(int &x, const int &x_max) {
 579     if (x < 0) {
 580         x += x_max;
 581     }
 582     if (x >= x_max) {
 583         x -= x_max;
 584     }
 585 }
 586
 587 void alternateNeighborhoodSearch::GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
 588     rvec coordinates, box_vector_length;
 589     num_of_miniboxes.resize(0);
 590     num_of_miniboxes.resize(3);
 591     for (std::size_t i{XX}; i <= ZZ; ++i) {
 592         box_vector_length[i] = std::sqrt(
 593                     std::pow(fr.box[i][XX], 2) + std::pow(fr.box[i][YY], 2) + std::pow(fr.box[i][ZZ], 2));
 594         num_of_miniboxes[i] = std::floor((box_vector_length[i] / cutoff)) + 1;
 595     }
 596     MiniBoxesMap.resize(0);
 597     MiniBoxesReverseMap.resize(0);
 598     MiniBoxesMap.resize(num_of_miniboxes[XX], std::vector<std::vector<std::vector<std::size_t> > >(
 599                              num_of_miniboxes[YY], std::vector<std::vector<std::size_t> >(
 600                              num_of_miniboxes[ZZ], std::vector<std::size_t>(
 601                              0))));
 602     MiniBoxesReverseMap.resize(IndexMap.size(), std::vector<std::size_t>(3));
 603     for (std::vector<ResInfo>::const_iterator i {IndexMap.begin()}; i != IndexMap.end(); ++i) {
 604         for (std::size_t j{XX}; j <= ZZ; ++j) {
 605             coordinates[j] = fr.x[i->getIndex(backboneAtomTypes::AtomCA)][j];
 606             FixAtomCoordinates(coordinates[j], box_vector_length[j]);
 607         }
 608         MiniBoxesMap[std::floor(coordinates[XX] / cutoff)][std::floor(coordinates[YY] / cutoff)][std::floor(
 609                           coordinates[ZZ] / cutoff)].push_back(i - IndexMap.begin());
 610         for (std::size_t j{XX}; j <= ZZ; ++j){
 611             MiniBoxesReverseMap[i - IndexMap.begin()][j] = std::floor(coordinates[j] / cutoff);
 612         }
 613     }
 614 }
 615
 616 void alternateNeighborhoodSearch::AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap){
 617     GetMiniBoxesMap(fr, IndexMap);
 618     MiniBoxSize[XX] = MiniBoxesMap.size();
 619     MiniBoxSize[YY] = MiniBoxesMap.front().size();
 620     MiniBoxSize[ZZ] = MiniBoxesMap.front().front().size();
 621     PairMap.resize(0);
 622     PairMap.resize(IndexMap.size(), std::vector<bool>(IndexMap.size(), false));
 623     ResiI = PairMap.begin();
 624     ResiJ = ResiI->begin();
 625
 626     for (std::vector<ResInfo>::const_iterator i = IndexMap.begin(); i != IndexMap.end(); ++i){
 627         for (offset[XX] = -1; offset[XX] <= 1; ++offset[XX]) {
 628             for (offset[YY] = -1; offset[YY] <= 1; ++offset[YY]) {
 629                 for (offset[ZZ] = -1; offset[ZZ] <= 1; ++offset[ZZ]) {
 630                     for (std::size_t k{XX}; k <= ZZ; ++k) {
 631                         fixBox[k] = MiniBoxesReverseMap[i - IndexMap.begin()][k] + offset[k];
 632                         ReCalculatePBC(fixBox[k], MiniBoxSize[k]);
 633                     }
 634                     for (std::size_t j{0}; j < MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]].size(); ++j) {
 635                         if ( (i - IndexMap.begin()) != MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]){
 636                             PairMap[i - IndexMap.begin()][MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]] = true;
 637                             PairMap[MiniBoxesMap[fixBox[XX]][fixBox[YY]][fixBox[ZZ]][j]][i - IndexMap.begin()] = true;
 638                         }
 639                     }
 640                 }
 641             }
 642         }
 643     }
 644 }
 645
 646 bool alternateNeighborhoodSearch::findNextPair(){
 647
 648     if(!PairMap.size()){
 649         return false;
 650     }
 651
 652     for(; ResiI != PairMap.end(); ++ResiI, ResiJ = ResiI->begin() ){
 653         for(; ResiJ != ResiI->end(); ++ResiJ){
 654             if(*ResiJ){
 655                 resiIpos = ResiI - PairMap.begin();
 656                 resiJpos = ResiJ - ResiI->begin();
 657                 if ( ResiJ != ResiI->end() ){
 658                     ++ResiJ;
 659                 }
 660                 else if (ResiI != PairMap.end()) {
 661                     ++ResiI;
 662                     ResiJ = ResiI->begin();
 663                 }
 664                 else {
 665                     return false; // ???
 666                 }
 667                 return true;
 668             }
 669         }
 670     }
 671
 672     return false;
 673 }
 674
 675 std::size_t alternateNeighborhoodSearch::getResiI() const {
 676     return resiIpos;
 677 }
 678
 679 std::size_t alternateNeighborhoodSearch::getResiJ() const {
 680     return resiJpos;
 681 }
 682
 683
 684 DsspTool::DsspStorage DsspTool::Storage;
 685
 686 DsspTool::DsspTool(){
 687 }
 688
 689 void DsspTool::calculateBends(const t_trxframe &fr, const t_pbc *pbc)
 690 {
 691    const float benddegree{ 70.0 }, maxdist{ 2.5 };
 692    float       degree{ 0 }, vdist{ 0 }, vprod{ 0 };
 693    gmx::RVec   a{ 0, 0, 0 }, b{ 0, 0, 0 };
 694    for (std::size_t i{ 0 }; i + 1 < IndexMap.size(); ++i)
 695    {
 696        if (CalculateAtomicDistances(static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
 697                                     static_cast<int>(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)),
 698                                     fr,
 699                                     pbc)
 700            > maxdist)
 701        {
 702            PatternSearch.SecondaryStructuresStatusMap[i].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i + 1]);
 703            PatternSearch.SecondaryStructuresStatusMap[i + 1].setBreak(&PatternSearch.SecondaryStructuresStatusMap[i]);
 704
 705 //           std::cout << "Break between " << i + 1 << " and " << i + 2 << std::endl;
 706        }
 707    }
 708    for (std::size_t i{ 2 }; i + 2 < IndexMap.size() ; ++i)
 709    {
 710        if (PatternSearch.SecondaryStructuresStatusMap[i - 2].getStatus(secondaryStructureTypes::Break) ||
 711            PatternSearch.SecondaryStructuresStatusMap[i - 1].getStatus(secondaryStructureTypes::Break) ||
 712            PatternSearch.SecondaryStructuresStatusMap[i].getStatus(secondaryStructureTypes::Break) ||
 713            PatternSearch.SecondaryStructuresStatusMap[i + 1].getStatus(secondaryStructureTypes::Break)
 714           )
 715        {
 716            continue;
 717        }
 718        for (int j{ 0 }; j < 3; ++j)
 719        {
 720            a[j] = fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j]
 721                   - fr.x[IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA)][j];
 722            b[j] = fr.x[IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA)][j]
 723                   - fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)][j];
 724        }
 725        vdist = (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
 726        vprod = CalculateAtomicDistances(IndexMap[i - 2].getIndex(backboneAtomTypes::AtomCA),
 727                                         IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
 728                                         fr,
 729                                         pbc)
 730                * gmx::c_angstrom / gmx::c_nano
 731                * CalculateAtomicDistances(IndexMap[i].getIndex(backboneAtomTypes::AtomCA),
 732                                           IndexMap[i + 2].getIndex(backboneAtomTypes::AtomCA),
 733                                           fr,
 734                                           pbc)
 735                * gmx::c_angstrom / gmx::c_nano;
 736        degree = std::acos(vdist / vprod) * gmx::c_rad2Deg;
 737        if (degree > benddegree)
 738        {
 739            PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Bend);
 740        }
 741    }
 742 }
 743
 744 float DsspTool::CalculateDihedralAngle(const int &A, const int &B, const int &C, const int &D, const t_trxframe &fr, const t_pbc *pbc){
 745     float result{360}, u{}, v{};
 746     gmx::RVec v12{}, v43{}, z{}, p{}, x{}, y{};
 747     pbc_dx(pbc, fr.x[A], fr.x[B], v12.as_vec());
 748     pbc_dx(pbc, fr.x[D], fr.x[C], v43.as_vec());
 749     pbc_dx(pbc, fr.x[B], fr.x[C], z.as_vec());
 750
 751     for(std::size_t j {XX}; j <= ZZ; ++j){
 752         v12[j] *= gmx::c_nm2A;
 753         v43[j] *= gmx::c_nm2A;
 754         z[j] *= gmx::c_nm2A;
 755     }
 756     for(std::size_t i{XX}, j{j + 1}, k{i + 2}; i <= ZZ; ++i, ++j, ++k){
 757         if (j > 2){
 758             j -= 3;
 759         }
 760         if (k > 2){
 761             k -= 3;
 762         }
 763         p[i] = (z[j] * v12[k]) - (z[k] * v12[j]);
 764         x[i] = (z[j] * v43[k]) - (z[k] * v43[j]);
 765
 766     }
 767     for(std::size_t i{XX}, j{j + 1}, k{i + 2}; i <= ZZ; ++i, ++j, ++k){
 768         if (j > 2){
 769             j -= 3;
 770         }
 771         if (k > 2){
 772             k -= 3;
 773         }
 774         y[i] = (z[j] * x[k]) - (z[k] * x[j]);
 775     }
 776
 777 //    std::cout << "v12 = " << v12[0] << ", " << v12[1] << ", " << v12[2] << std::endl;
 778 //    std::cout << "v43 = " << v43[0] << ", " << v43[1] << ", " << v43[2] << std::endl;
 779 //    std::cout << "z = " << z[0] << ", " << z[1] << ", " << z[2] << std::endl;
 780 //    std::cout << "p = " << p[0] << ", " << p[1] << ", " << p[2] << std::endl;
 781 //    std::cout << "x = " << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
 782 //    std::cout << "y = " << y[0] << ", " << y[1] << ", " << y[2] << std::endl;
 783
 784     u = (x[XX] * x[XX]) + (x[YY] * x[YY]) + (x[ZZ] * x[ZZ]);
 785     v = (y[XX] * y[XX]) + (y[YY] * y[YY]) + (y[ZZ] * y[ZZ]);
 786
 787 //    std::cout << "u = " << u << std::endl;
 788 //    std::cout << "v = " << v << std::endl;
 789
 790     if (u > 0 and v > 0){
 791         u = ((p[XX] * x[XX]) + (p[YY] * x[YY]) + (p[ZZ] * x[ZZ])) / std::sqrt(u);
 792         v = ((p[XX] * y[XX]) + (p[YY] * y[YY]) + (p[ZZ] * y[ZZ])) / std::sqrt(v);
 793 //        std::cout << "new u = " << u << std::endl;
 794 //        std::cout << "new v = " << v << std::endl;
 795         if (u != 0 or v != 0){
 796             result = std::atan2(v, u) * gmx::c_rad2Deg;
 797 //            std::cout << "result = " << result << std::endl;
 798         }
 799     }
 800     return result;
 801 }
 802
 803 void DsspTool::calculateDihedrals(const t_trxframe &fr, const t_pbc *pbc){
 804     const float epsilon = 29;
 805     const float phi_min = -75 - epsilon; // -104
 806     const float phi_max = -75 + epsilon; // -46
 807     const float psi_min = 145 - epsilon; // 116
 808     const float psi_max = 145 + epsilon; // 176
 809     std::vector<float>          phi(IndexMap.size(), 360), psi(IndexMap.size(), 360);
 810 //    phi.resize(0);
 811 //    psi.resize(0);
 812 //    phi.resize(IndexMap.size(), 360);
 813 //    psi.resize(IndexMap.size(), 360);
 814
 815     for (std::size_t i = 1; i + 1 < IndexMap.size(); ++i){ // TODO add index verifictaion (check if those atom indexes exist)
 816 //        std::cout << "For resi " << i << " :" << std::endl;
 817         phi[i] = CalculateDihedralAngle(static_cast<int>(IndexMap[i - 1].getIndex(backboneAtomTypes::AtomC)),
 818                                         static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomN)),
 819                                         static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomCA)),
 820                                         static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
 821                                         fr,
 822                                         pbc);
 823         psi[i] = CalculateDihedralAngle(static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomN)),
 824                                         static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomCA)),
 825                                         static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
 826                                         static_cast<int>(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)),
 827                                         fr,
 828                                         pbc);
 829 //        std::cout << "For " << i << " phi = " << phi[i] << ", psi = " << psi[i] << std::endl;
 830 //        std::cout << "phi[" << i << "] = " << phi[i] << std::endl;
 831 //        std::cout << "psi[" << i << "] = " << psi[i] << std::endl;
 832     }
 833
 834     for (std::size_t i = 1; i + 3 < IndexMap.size(); ++i){
 835         switch (initParams.pp_stretch){
 836         case 2: {
 837             if (phi_min > phi[i] or phi[i] > phi_max or
 838                 phi_min > phi[i + 1] or phi[i + 1]> phi_max){
 839                 continue;
 840             }
 841
 842             if (psi_min > psi[i] or psi[i] > psi_max or
 843                 psi_min > psi[i + 1] or psi[i + 1] > psi_max){
 844                 continue;
 845             }
 846
 847             switch (PatternSearch.SecondaryStructuresStatusMap[i].getStatus(turnsTypes::Turn_PP)){
 848                 case HelixPositions::None:
 849                     PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start, turnsTypes::Turn_PP);
 850                     break;
 851
 852                 case HelixPositions::End:
 853                     PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start_AND_End, turnsTypes::Turn_PP);
 854                     break;
 855
 856                 default:
 857                     break;
 858             }
 859
 860             PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(HelixPositions::End, turnsTypes::Turn_PP);
 861             /* Пропустил проверку того, что заменяемая ак - петля */
 862             PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Helix_PP);
 863             PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(secondaryStructureTypes::Helix_PP);
 864             break;
 865         }
 866         case 3:{
 867             if (phi_min > phi[i] or phi[i] > phi_max or
 868                 phi_min > phi[i + 1] or phi[i + 1]> phi_max or
 869                 phi_min > phi[i + 2] or phi[i + 2]> phi_max){
 870                 continue;
 871             }
 872
 873             if (psi_min > psi[i] or psi[i] > psi_max or
 874                 psi_min > psi[i + 1] or psi[i + 1] > psi_max or
 875                 psi_min > psi[i + 2] or psi[i + 2] > psi_max){
 876                 continue;
 877             }
 878
 879             switch (PatternSearch.SecondaryStructuresStatusMap[i].getStatus(turnsTypes::Turn_PP)){
 880                 case HelixPositions::None:
 881                     PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start, turnsTypes::Turn_PP);
 882                     break;
 883
 884                 case HelixPositions::End:
 885                     PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start_AND_End, turnsTypes::Turn_PP);
 886                     break;
 887
 888                 default:
 889                     break;
 890             }
 891
 892             PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(HelixPositions::Middle, turnsTypes::Turn_PP);
 893             PatternSearch.SecondaryStructuresStatusMap[i + 2].setStatus(HelixPositions::End, turnsTypes::Turn_PP);
 894             /* Пропустил проверку того, что заменяемая ак - петля */
 895             PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Helix_PP);
 896             PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(secondaryStructureTypes::Helix_PP);
 897             PatternSearch.SecondaryStructuresStatusMap[i + 2].setStatus(secondaryStructureTypes::Helix_PP);
 898
 899             break;
 900         }
 901         default:
 902             throw std::runtime_error("Unsupported stretch length");
 903         }
 904     }
 905
 906 }
 907
 908 void DsspTool::calculateHBondEnergy(ResInfo& Donor,
 909                        ResInfo& Acceptor,
 910                        const t_trxframe&          fr,
 911                        const t_pbc*               pbc)
 912 {
 913    /*
 914     * DSSP uses eq from dssp 2.x
 915     * kCouplingConstant = 27.888,  //  = 332 * 0.42 * 0.2
 916     * E = k * (1/rON + 1/rCH - 1/rOH - 1/rCN) where CO comes from one AA and NH from another
 917     * if R is in A
 918     * Hbond if E < -0.5
 919     *
 920     * For the note, H-Bond Donor is N-H («Donor of H») and H-Bond Acceptor is C=O («Acceptor of H»)
 921     *
 922     */
 923
 924     if (CalculateAtomicDistances(
 925                 Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
 926         >= minimalCAdistance)
 927     {
 928         return void();
 929     }
 930
 931     const float kCouplingConstant = 27.888;
 932     const float minimalAtomDistance{ 0.5 },
 933             minEnergy{ -9.9 };
 934     float HbondEnergy{ 0 };
 935     float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 };
 936
 937 //    std::cout << "For Donor №" << Donor.info->nr - 1 << " and Accpetor №" << Acceptor.info->nr - 1 << std::endl;
 938
 939     if( !(Donor.is_proline) && (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
 940                                 && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || initParams.addHydrogens ) ) ){ // TODO
 941         distanceNO = CalculateAtomicDistances(
 942                Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 943         distanceNC = CalculateAtomicDistances(
 944                Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 945         if (initParams.addHydrogens){
 946             if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
 947                rvec atomH{};
 948                float prevCODist {CalculateAtomicDistances(Donor.prevResi->getIndex(backboneAtomTypes::AtomC), Donor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)};
 949                for (int i{XX}; i <= ZZ; ++i){
 950                    float prevCO = fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
 951                    atomH[i] = fr.x[Donor.getIndex(backboneAtomTypes::AtomH)][i]; // Но на самом деле берутся координаты N
 952                    atomH[i] += prevCO / prevCODist;
 953                }
 954                distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 955                distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 956             }
 957             else{
 958                 distanceHO = distanceNO;
 959                 distanceHC = distanceNC;
 960             }
 961        }
 962        else {
 963            distanceHO = CalculateAtomicDistances(
 964                    Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
 965            distanceHC = CalculateAtomicDistances(
 966                    Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
 967        }
 968        if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
 969         || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
 970        {
 971             HbondEnergy = minEnergy;
 972        }
 973        else{
 974            HbondEnergy =
 975                    kCouplingConstant
 976                    * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
 977        }
 978
 979 //       std::cout << "CA-CA distance: " << CalculateAtomicDistances(
 980 //                        Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc) << std::endl;
 981 //       std::cout << "N-O distance: " << distanceNO << std::endl;
 982 //       std::cout << "N-C distance: " << distanceNC << std::endl;
 983 //       std::cout << "H-O distance: " << distanceHO << std::endl;
 984 //       std::cout << "H-C distance: " << distanceHC << std::endl;
 985
 986        HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
 987
 988 //       if ( HbondEnergy < minEnergy ){ // I don't think that this is correct
 989 //            HbondEnergy = minEnergy;
 990 //       }
 991
 992 //       std::cout << "Calculated energy = " << HbondEnergy << std::endl;
 993     }
 994 //    else{
 995 //        std::cout << "Donor Is Proline" << std::endl;
 996 //    }
 997
 998     if (HbondEnergy < Donor.acceptorEnergy[0]){
 999            Donor.acceptor[1] = Donor.acceptor[0];
1000            Donor.acceptor[0] = Acceptor.info;
1001            Donor.acceptorEnergy[0] = HbondEnergy;
1002     }
1003     else if (HbondEnergy < Donor.acceptorEnergy[1]){
1004            Donor.acceptor[1] = Acceptor.info;
1005            Donor.acceptorEnergy[1] = HbondEnergy;
1006     }
1007
1008     if (HbondEnergy < Acceptor.donorEnergy[0]){
1009            Acceptor.donor[1] = Acceptor.donor[0];
1010            Acceptor.donor[0] = Donor.info;
1011            Acceptor.donorEnergy[0] = HbondEnergy;
1012     }
1013     else if (HbondEnergy < Acceptor.donorEnergy[1]){
1014            Acceptor.donor[1] = Donor.info;
1015            Acceptor.donorEnergy[1] = HbondEnergy;
1016     }
1017 }
1018
1019
1020 /* Calculate Distance From B to A */
1021 float DsspTool::CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
1022 {
1023    gmx::RVec r{ 0, 0, 0 };
1024    pbc_dx(pbc, fr.x[A], fr.x[B], r.as_vec());
1025    return r.norm() * gmx::c_nm2A; // НЕ ТРОГАТЬ
1026 }
1027
1028 /* Calculate Distance From B to A, where A is only fake coordinates */
1029 float DsspTool::CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
1030 {
1031    gmx::RVec r{ 0, 0, 0 };
1032    pbc_dx(pbc, A, fr.x[B], r.as_vec());
1033    return r.norm() * gmx::c_nm2A; // НЕ ТРОГАТЬ
1034 }
1035
1036 void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz)
1037 {
1038    initParams = initParamz;
1039    ResInfo _backboneAtoms;
1040    std::size_t                 i{ 0 };
1041    std::string proLINE;
1042    int resicompare{ top.atoms()->atom[static_cast<std::size_t>(*(initParams.sel_.atomIndices().begin()))].resind };
1043    IndexMap.resize(0);
1044    IndexMap.push_back(_backboneAtoms);
1045    IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
1046    proLINE = *(IndexMap[i].info->name);
1047    if( proLINE.compare("PRO") == 0 ){
1048        IndexMap[i].is_proline = true;
1049    }
1050
1051    for (gmx::ArrayRef<const int>::iterator ai{ initParams.sel_.atomIndices().begin() }; (ai != initParams.sel_.atomIndices().end()); ++ai){
1052        if (resicompare != top.atoms()->atom[static_cast<std::size_t>(*ai)].resind)
1053        {
1054            ++i;
1055            resicompare = top.atoms()->atom[static_cast<std::size_t>(*ai)].resind;
1056            IndexMap.emplace_back(_backboneAtoms);
1057            IndexMap[i].info = &(top.atoms()->resinfo[resicompare]);
1058            proLINE = *(IndexMap[i].info->name);
1059            if( proLINE.compare("PRO") == 0 ){
1060                IndexMap[i].is_proline = true;
1061            }
1062
1063        }
1064        std::string atomname(*(top.atoms()->atomname[static_cast<std::size_t>(*ai)]));
1065        if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA])
1066        {
1067            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomCA)] = *ai;
1068        }
1069        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC])
1070        {
1071            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomC)] = *ai;
1072        }
1073        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO])
1074        {
1075            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomO)] = *ai;
1076        }
1077        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN])
1078        {
1079            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomN)] = *ai;
1080            if (initParamz.addHydrogens == true){
1081                IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
1082            }
1083        }
1084        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH] && initParamz.addHydrogens == false) // Юзать водород в структуре
1085        {
1086            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
1087        }
1088
1089
1090
1091 //       if( atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO]
1092 //       || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH]){
1093 //           std::cout << "Atom " << atomname << " №" << *ai << " From Resi " << *(top.atoms()->resinfo[i].name) << " №" << resicompare << std::endl;
1094 //       }
1095    }
1096
1097    for (std::size_t j {1}; j < IndexMap.size(); ++j){
1098        IndexMap[j].prevResi = &(IndexMap[j - 1]);
1099
1100        IndexMap[j - 1].nextResi = &(IndexMap[j]);
1101
1102 //           std::cout << "Resi " << IndexMap[i].info->nr << *(IndexMap[i].info->name) << std::endl;
1103 //           std::cout << "Prev resi is " << IndexMap[i].prevResi->info->nr << *(IndexMap[i].prevResi->info->name) << std::endl;
1104 //           std::cout << "Prev resi's next resi is " << IndexMap[i - 1].nextResi->info->nr << *(IndexMap[i - 1].nextResi->info->name) << std::endl;
1105 //         std::cout << IndexMap[j].prevResi->info->nr;
1106 //         std::cout << *(IndexMap[j].prevResi->info->name) ;
1107 //         std::cout << " have CA = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomCA) ;
1108 //         std::cout << " C = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomC);
1109 //         std::cout << " O = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomO);
1110 //         std::cout << " N = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomN);
1111 //         std::cout << " H = " << IndexMap[j].prevResi->getIndex(backboneAtomTypes::AtomH) << std::endl;
1112    }
1113
1114    nres = i + 1;
1115 }
1116
1117 void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc)
1118 {
1119
1120     switch(initParams.NBS){
1121     case (NBSearchMethod::Classique): {
1122
1123         // store positions of CA atoms to use them for nbSearch
1124         std::vector<gmx::RVec> positionsCA_;
1125         for (std::size_t i{ 0 }; i < IndexMap.size(); ++i)
1126         {
1127             positionsCA_.emplace_back(fr.x[IndexMap[i].getIndex(backboneAtomTypes::AtomCA)]);
1128         }
1129
1130         AnalysisNeighborhood nb_;
1131         nb_.setCutoff(initParams.cutoff_);
1132         AnalysisNeighborhoodPositions       nbPos_(positionsCA_);
1133         gmx::AnalysisNeighborhoodSearch     start      = nb_.initSearch(pbc, nbPos_);
1134         gmx::AnalysisNeighborhoodPairSearch pairSearch = start.startPairSearch(nbPos_);
1135         gmx::AnalysisNeighborhoodPair       pair;
1136         while (pairSearch.findNextPair(&pair))
1137         {
1138             if(CalculateAtomicDistances(
1139                         IndexMap[pair.refIndex()].getIndex(backboneAtomTypes::AtomCA), IndexMap[pair.testIndex()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
1140                 < minimalCAdistance){
1141                 calculateHBondEnergy(IndexMap[pair.refIndex()], IndexMap[pair.testIndex()], fr, pbc);
1142                 if (IndexMap[pair.testIndex()].info != IndexMap[pair.refIndex() + 1].info){
1143                     calculateHBondEnergy(IndexMap[pair.testIndex()], IndexMap[pair.refIndex()], fr, pbc);
1144                 }
1145             }
1146         }
1147
1148         break;
1149     }
1150     case (NBSearchMethod::Experimental): { // TODO FIX
1151
1152         alternateNeighborhoodSearch as_;
1153
1154         as_.setCutoff(initParams.cutoff_);
1155
1156         as_.AltPairSearch(fr, IndexMap);
1157
1158         while (as_.findNextPair()){
1159             if(CalculateAtomicDistances(
1160                         IndexMap[as_.getResiI()].getIndex(backboneAtomTypes::AtomCA), IndexMap[as_.getResiJ()].getIndex(backboneAtomTypes::AtomCA), fr, pbc)
1161                 < minimalCAdistance){
1162                 calculateHBondEnergy(IndexMap[as_.getResiI()], IndexMap[as_.getResiJ()], fr, pbc);
1163                 if (IndexMap[as_.getResiJ()].info != IndexMap[as_.getResiI() + 1].info){
1164                     calculateHBondEnergy(IndexMap[as_.getResiJ()], IndexMap[as_.getResiI()], fr, pbc);
1165                 }
1166             }
1167         }
1168
1169         break;
1170     }
1171     default: {
1172
1173         for(std::vector<ResInfo>::iterator Donor {IndexMap.begin()}; Donor != IndexMap.end() ; ++Donor){
1174             for(std::vector<ResInfo>::iterator Acceptor {Donor + 1} ; Acceptor != IndexMap.end() ; ++Acceptor){
1175                 if(CalculateAtomicDistances(
1176                             Donor->getIndex(backboneAtomTypes::AtomCA), Acceptor->getIndex(backboneAtomTypes::AtomCA), fr, pbc)
1177                     < minimalCAdistance){
1178                     calculateHBondEnergy(*Donor, *Acceptor, fr, pbc);
1179                     if (Acceptor != Donor + 1){
1180                         calculateHBondEnergy(*Acceptor, *Donor, fr, pbc);
1181                     }
1182                 }
1183             }
1184         }
1185         break;
1186     }
1187     }
1188
1189
1190 //    for(std::size_t i {0}; i < IndexMap.size(); ++i){
1191 //        std::cout << IndexMap[i].info->nr << " " << *(IndexMap[i].info->name) << std::endl;
1192 //    }
1193
1194    PatternSearch.initiateSearch(IndexMap, initParams.PPHelices, initParams.pp_stretch);
1195    calculateBends(fr, pbc);
1196    calculateDihedrals(fr, pbc);
1197    Storage.storageData(frnr, PatternSearch.patternSearch());
1198
1199 }
1200
1201 std::vector<std::pair<int, std::string>> DsspTool::getData(){
1202     return Storage.returnData();
1203 }
1204
1205 } // namespace analysismodules
1206
1207 } // namespace gmx