AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAa
authorMax <Infinity2573@gmail.com>
Sun, 25 Sep 2022 11:28:49 +0000 (14:28 +0300)
committerMax <Infinity2573@gmail.com>
Sun, 25 Sep 2022 11:28:49 +0000 (14:28 +0300)
src/dssptools.cpp

index 26b6e00174991510535a52413fc7a7a7f065ea27..f13755a12d78d1e91e1e5a526febbcc488a3fcd3 100644 (file)
@@ -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<std::size_t>(backboneAtomTypes::AtomN)] = *ai;
+           if (initParamz.addHydrogens == true){
+               IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
+           }
        }
        else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH] && initParamz.addHydrogens == false) // Юзать водород в структуре
        {
            IndexMap[i]._backboneIndices[static_cast<std::size_t>(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;