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