From: Max Date: Sun, 25 Sep 2022 11:28:49 +0000 (+0300) Subject: AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAa X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=36f473c9132414e6bcfe4edbde3da16b29c5d51f;p=alexxy%2Fgromacs-dssp.git AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAa --- diff --git a/src/dssptools.cpp b/src/dssptools.cpp index 26b6e00..f13755a 100644 --- a/src/dssptools.cpp +++ b/src/dssptools.cpp @@ -345,35 +345,35 @@ std::string secondaryStructures::patternSearch(){ // std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) << std::endl; // } -// std::cout.precision(5); -// for(std::size_t i{0}; i < ResInfoMap->size(); ++i, std::cout << std::endl << std::endl){ -// std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) ; -// if ( (*ResInfoMap)[i].donor[0] != nullptr ){ -// std::cout << " has donor[0] = " << (*ResInfoMap)[i].donor[0]->nr << " " << *((*ResInfoMap)[i].donor[0]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[0] << " and" ; -// } -// else { -// std::cout << " has no donor[0] and" ; -// } -// if ( (*ResInfoMap)[i].acceptor[0] != nullptr ){ -// std::cout << " has acceptor[0] = " << (*ResInfoMap)[i].acceptor[0]->nr << " " << *((*ResInfoMap)[i].acceptor[0]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[0] ; -// } -// else { -// std::cout << " has no acceptor[0]" ; -// } -// std::cout << std::endl << "Also, " << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name); -// if ( (*ResInfoMap)[i].donor[1] != nullptr ){ -// std::cout << " has donor[1] = " << (*ResInfoMap)[i].donor[1]->nr << " " << *((*ResInfoMap)[i].donor[1]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[1] << " and" ; -// } -// else { -// std::cout << " has no donor[1] and" ; -// } -// if ( (*ResInfoMap)[i].acceptor[1] != nullptr ){ -// std::cout << " has acceptor[1] = " << (*ResInfoMap)[i].acceptor[1]->nr << " " << *((*ResInfoMap)[i].acceptor[1]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[1] ; -// } -// else { -// std::cout << " has no acceptor[1]" ; -// } -// } + std::cout.precision(5); + for(std::size_t i{0}; i < ResInfoMap->size(); ++i, std::cout << std::endl << std::endl){ + std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) ; + if ( (*ResInfoMap)[i].donor[0] != nullptr ){ + std::cout << " has donor[0] = " << (*ResInfoMap)[i].donor[0]->nr << " " << *((*ResInfoMap)[i].donor[0]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[0] << " and" ; + } + else { + std::cout << " has no donor[0] and" ; + } + if ( (*ResInfoMap)[i].acceptor[0] != nullptr ){ + std::cout << " has acceptor[0] = " << (*ResInfoMap)[i].acceptor[0]->nr << " " << *((*ResInfoMap)[i].acceptor[0]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[0] ; + } + else { + std::cout << " has no acceptor[0]" ; + } + std::cout << std::endl << "Also, " << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name); + if ( (*ResInfoMap)[i].donor[1] != nullptr ){ + std::cout << " has donor[1] = " << (*ResInfoMap)[i].donor[1]->nr << " " << *((*ResInfoMap)[i].donor[1]->name) << " with E = " << (*ResInfoMap)[i].donorEnergy[1] << " and" ; + } + else { + std::cout << " has no donor[1] and" ; + } + if ( (*ResInfoMap)[i].acceptor[1] != nullptr ){ + std::cout << " has acceptor[1] = " << (*ResInfoMap)[i].acceptor[1]->nr << " " << *((*ResInfoMap)[i].acceptor[1]->name) << " with E = " << (*ResInfoMap)[i].acceptorEnergy[1] ; + } + else { + std::cout << " has no acceptor[1]" ; + } + } /*Write Data*/ @@ -616,6 +616,9 @@ void DsspTool::calculateHBondEnergy(ResInfo& Donor, * E = k * (1/rON + 1/rCH - 1/rOH - 1/rCN) where CO comes from one AA and NH from another * if R is in A * Hbond if E < -0.5 + * + * For the note, H-Bond Donor is N-H («Donor of H») and H-Bond Acceptor is C=O («Acceptor of H») + * */ const float kCouplingConstant = 27.888; @@ -624,26 +627,26 @@ void DsspTool::calculateHBondEnergy(ResInfo& Donor, float HbondEnergy{ 0 }; float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 }; - if( !(Acceptor.is_proline) ){ - if (Donor.getIndex(backboneAtomTypes::AtomC) && Donor.getIndex(backboneAtomTypes::AtomO) - && Acceptor.getIndex(backboneAtomTypes::AtomN) && ( Acceptor.getIndex(backboneAtomTypes::AtomH) || (initParams.addHydrogens) ) ) // Kinda ew + if( !(Donor.is_proline) ){ + if (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO) + && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || (initParams.addHydrogens) ) ) { distanceNO = CalculateAtomicDistances( - Acceptor.getIndex(backboneAtomTypes::AtomN), Donor.getIndex(backboneAtomTypes::AtomO), fr, pbc); + Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc); distanceNC = CalculateAtomicDistances( - Acceptor.getIndex(backboneAtomTypes::AtomN), Donor.getIndex(backboneAtomTypes::AtomC), fr, pbc); + Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc); if (initParams.addHydrogens){ - if (Acceptor.prevResi != nullptr && Acceptor.prevResi->getIndex(backboneAtomTypes::AtomC) && Acceptor.prevResi->getIndex(backboneAtomTypes::AtomO)){ + if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){ rvec atomH{}; - float prevCODist {CalculateAtomicDistances(Acceptor.prevResi->getIndex(backboneAtomTypes::AtomC), Acceptor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)}; + float prevCODist {CalculateAtomicDistances(Donor.prevResi->getIndex(backboneAtomTypes::AtomC), Donor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)}; for (int i{XX}; i <= ZZ; ++i){ - float prevCO = fr.x[Acceptor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Acceptor.prevResi->getIndex(backboneAtomTypes::AtomO)][i]; + float prevCO = fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i]; atomH[i] = prevCO / prevCODist; } - distanceHO = CalculateAtomicDistances(atomH, Donor.getIndex(backboneAtomTypes::AtomO), fr, pbc); - distanceHC = CalculateAtomicDistances(atomH, Donor.getIndex(backboneAtomTypes::AtomC), fr, pbc); + distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc); + distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc); } else{ distanceHO = distanceNO; @@ -652,13 +655,13 @@ void DsspTool::calculateHBondEnergy(ResInfo& Donor, } else { distanceHO = CalculateAtomicDistances( - Acceptor.getIndex(backboneAtomTypes::AtomH), Donor.getIndex(backboneAtomTypes::AtomO), fr, pbc); + Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc); distanceHC = CalculateAtomicDistances( - Acceptor.getIndex(backboneAtomTypes::AtomH), Donor.getIndex(backboneAtomTypes::AtomC), fr, pbc); + Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc); } if (CalculateAtomicDistances( - Acceptor.getIndex(backboneAtomTypes::AtomCA), Donor.getIndex(backboneAtomTypes::AtomCA), fr, pbc) + Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc) < minimalCAdistance) { if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance) @@ -689,9 +692,6 @@ void DsspTool::calculateHBondEnergy(ResInfo& Donor, } } } -// else { -// std::cerr << "PRO DETECTED! THIS IS RESI № " << Donor.info->nr << std::endl; //IT WORKS JUST FINE -// } if (HbondEnergy < Donor.acceptorEnergy[0]){ Donor.acceptor[1] = Donor.acceptor[0]; @@ -777,12 +777,17 @@ void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN]) { IndexMap[i]._backboneIndices[static_cast(backboneAtomTypes::AtomN)] = *ai; + if (initParamz.addHydrogens == true){ + IndexMap[i]._backboneIndices[static_cast(backboneAtomTypes::AtomH)] = *ai; + } } else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH] && initParamz.addHydrogens == false) // Юзать водород в структуре { IndexMap[i]._backboneIndices[static_cast(backboneAtomTypes::AtomH)] = *ai; } + + // if( atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomC] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomO] // || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomN] || atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH]){ // std::cout << "Atom " << atomname << " №" << *ai << " From Resi " << *(top.atoms()->resinfo[i].name) << " №" << resicompare << std::endl;