From: Max <Infinity2573@gmail.com>
Date: Mon, 26 Sep 2022 11:10:46 +0000 (+0300)
Subject: Fixes fixes fixes
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=34d6a0c40712a6d356e629b4f880a193b9360eb0;p=alexxy%2Fgromacs-dssp.git

Fixes fixes fixes
---

diff --git a/src/dssptools.cpp b/src/dssptools.cpp
index 544ef07..939d29c 100644
--- a/src/dssptools.cpp
+++ b/src/dssptools.cpp
@@ -43,7 +43,7 @@
 */
 
 /*
-    There's something wrong with energy of the last residue
+    There's something wrong with energy calculations of redidues with E ≈ -0
 */
 
 
@@ -639,79 +639,60 @@ void DsspTool::calculateHBondEnergy(ResInfo& Donor,
     *
     */
 
+    if (CalculateAtomicDistances(
+                Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
+        >= minimalCAdistance)
+    {
+        return void();
+    }
+
     const float kCouplingConstant = 27.888;
     const float minimalAtomDistance{ 0.5 },
             minEnergy{ -9.9 };
     float HbondEnergy{ 0 };
     float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 };
 
-    std::cout << "For Donor №" << Donor.info->nr - 1 << " and accpetor №" << Acceptor.info->nr - 1 << std::endl;
-
-   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(
-                   Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
-           distanceNC = CalculateAtomicDistances(
-                   Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
-
-           if (initParams.addHydrogens){
-               if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
-                   rvec atomH{};
-                   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[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
-                       atomH[i] = fr.x[Donor.getIndex(backboneAtomTypes::AtomH)][i]; // Но на самом деле берутся координаты N
-                       atomH[i] += prevCO / prevCODist;
-                   }
-                   distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
-                   distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
+    std::cout << "For Donor №" << Donor.info->nr - 1 << " and Accpetor №" << Acceptor.info->nr - 1 << std::endl;
+
+    if( !(Donor.is_proline) && (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
+                                && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || initParams.addHydrogens ) ) ){ // TODO
+        distanceNO = CalculateAtomicDistances(
+               Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
+        distanceNC = CalculateAtomicDistances(
+               Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
+        if (initParams.addHydrogens){
+            if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
+               rvec atomH{};
+               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[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
+                   atomH[i] = fr.x[Donor.getIndex(backboneAtomTypes::AtomH)][i]; // Но на самом деле берутся координаты N
+                   atomH[i] += prevCO / prevCODist;
                }
-               else{
-                   distanceHO = distanceNO;
-                   distanceHC = distanceNC;
-               }
-           }
-           else {
-               distanceHO = CalculateAtomicDistances(
-                       Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
-               distanceHC = CalculateAtomicDistances(
-                       Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
-           }
-
-           if (CalculateAtomicDistances(
-                       Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
-               < minimalCAdistance)
-           {
-               if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
-                   || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
-               {
-                   HbondEnergy = minEnergy;
-               }
-               else
-               {
-                   HbondEnergy =
-                           kCouplingConstant
-                           * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
-//                   HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
-
-//                   std::cout.precision(5);
-//                   std::cout << "Calculated ENERGY = " << HbondEnergy << std::endl;
-
-//                   if ( HbondEnergy == 0){
-//                       std::cout << "Calculated ENERGY = " << HbondEnergy << " For donor " << Donor.info->nr << " and acceptor " << Acceptor.info->nr << std::endl;
-//                   }
-
-                   if ( HbondEnergy < minEnergy ){
-                       HbondEnergy = minEnergy;
-                   }
-               }
-           }
+               distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
+               distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
+            }
+            else{
+                distanceHO = distanceNO;
+                distanceHC = distanceNC;
+            }
+       }
+       else {
+           distanceHO = CalculateAtomicDistances(
+                   Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
+           distanceHC = CalculateAtomicDistances(
+                   Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
+       }
+       if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
+        || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
+       {
+            HbondEnergy = minEnergy;
+       }
+       else{
+           HbondEnergy =
+                   kCouplingConstant
+                   * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
        }
-
-
 
        std::cout << "CA-CA distance: " << CalculateAtomicDistances(
                         Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc) << std::endl;
@@ -720,31 +701,37 @@ void DsspTool::calculateHBondEnergy(ResInfo& Donor,
        std::cout << "H-O distance: " << distanceHO << std::endl;
        std::cout << "H-C distance: " << distanceHC << std::endl;
 
-   }
-
+       HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
 
-   std::cout << "Calculated energy = " << HbondEnergy << std::endl;
+       if ( HbondEnergy < minEnergy ){
+            HbondEnergy = minEnergy;
+       }
 
+       std::cout << "Calculated energy = " << HbondEnergy << std::endl;
+    }
+    else{
+        std::cout << "Donor Is Proline" << std::endl;
+    }
 
-   if (HbondEnergy < Donor.acceptorEnergy[0]){
-       Donor.acceptor[1] = Donor.acceptor[0];
-       Donor.acceptor[0] = Acceptor.info;
-       Donor.acceptorEnergy[0] = HbondEnergy;
-   }
-   else if (HbondEnergy < Donor.acceptorEnergy[1]){
-       Donor.acceptor[1] = Acceptor.info;
-       Donor.acceptorEnergy[1] = HbondEnergy;
-   }
+    if (HbondEnergy < Donor.acceptorEnergy[0]){
+           Donor.acceptor[1] = Donor.acceptor[0];
+           Donor.acceptor[0] = Acceptor.info;
+           Donor.acceptorEnergy[0] = HbondEnergy;
+    }
+    else if (HbondEnergy < Donor.acceptorEnergy[1]){
+           Donor.acceptor[1] = Acceptor.info;
+           Donor.acceptorEnergy[1] = HbondEnergy;
+    }
 
-   if (HbondEnergy < Acceptor.donorEnergy[0]){
-       Acceptor.donor[1] = Acceptor.donor[0];
-       Acceptor.donor[0] = Donor.info;
-       Acceptor.donorEnergy[0] = HbondEnergy;
-   }
-   else if (HbondEnergy < Acceptor.donorEnergy[1]){
-       Acceptor.donor[1] = Donor.info;
-       Acceptor.donorEnergy[1] = HbondEnergy;
-   }
+    if (HbondEnergy < Acceptor.donorEnergy[0]){
+           Acceptor.donor[1] = Acceptor.donor[0];
+           Acceptor.donor[0] = Donor.info;
+           Acceptor.donorEnergy[0] = HbondEnergy;
+    }
+    else if (HbondEnergy < Acceptor.donorEnergy[1]){
+           Acceptor.donor[1] = Donor.info;
+           Acceptor.donorEnergy[1] = HbondEnergy;
+    }
 }