saddasd

diff --git a/src/dssp.cpp b/src/dssp.cpp
index 65734156b8a1c8a9ce27020eaa77d5d0b1a35f74..2757404ec25a17e99c72bc9b39522cdce343acb8 100644 (file)
--- a/src/dssp.cpp
+++ b/src/dssp.cpp
@@ -102,7 +102,9 @@ void Dssp::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* s
    options->addOption(
                EnumOption<NBSearchMethod>("algo").store(&initParams.NBS).required().defaultValue(NBSearchMethod::DSSP).enumValue(NBSearchMethodNames).description("Algorithm for search residues' neighbours"));
    options->addOption(
-               BooleanOption("pp").store(&initParams.PPHelices).defaultValue(true).description("Prefer Pi Helices"));
+               BooleanOption("h").store(&initParams.addHydrogens).defaultValue(false).description("Add hydrogens or not"));
+   options->addOption(
+               BooleanOption("p").store(&initParams.PPHelices).defaultValue(true).description("Prefer Pi Helices"));
    options->addOption(
                BooleanOption("v").store(&initParams.verbose).defaultValue(false).description(">:("));
    settings->setHelpText(desc);

diff --git a/src/dssptools.cpp b/src/dssptools.cpp
index ff761a8856a0e3ef2081000b02a891628dea23bb..5902024925b373af06368d1acfe0c6abe4d2e088 100644 (file)
--- a/src/dssptools.cpp
+++ b/src/dssptools.cpp
@@ -601,39 +601,47 @@ void DsspTool::calculateHBondEnergy(ResInfo& Donor,
     const float minimalAtomDistance{ 0.5 },
             minEnergy{ -9.9 };
     float HbondEnergy{ 0 };
-    float distanceON{ 0 }, distanceCH{ 0 }, distanceOH{ 0 }, distanceCN{ 0 };
+    float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 };
 
    if( !(Donor.is_proline) ){
+       if (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
+           && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || (initParams.addHydrogens) ) )  // Kinda ew
+       {
 
-       distanceON = CalculateAtomicDistances(
-               Donor.getIndex(backboneAtomTypes::AtomO), Acceptor.getIndex(backboneAtomTypes::AtomN), fr, pbc);
-       distanceCN = CalculateAtomicDistances(
-               Donor.getIndex(backboneAtomTypes::AtomC), Acceptor.getIndex(backboneAtomTypes::AtomN), fr, pbc);
-
-       switch(initParams.mode){
-       case modeType::Classique: {
-           distanceOH = distanceON;
-           distanceCH = distanceCN;
-           break;
-       }
-       default: {
-           distanceOH = CalculateAtomicDistances(
-                   Donor.getIndex(backboneAtomTypes::AtomO), Acceptor.getIndex(backboneAtomTypes::AtomH), fr, pbc);
-           distanceCH = CalculateAtomicDistances(
-                   Donor.getIndex(backboneAtomTypes::AtomC), Acceptor.getIndex(backboneAtomTypes::AtomH), fr, pbc);
-           break;
-       }
-       }
+           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 (Acceptor.prevResi != nullptr){
+                   rvec atomH{};
+                   float prevCODist {CalculateAtomicDistances(Acceptor.prevResi->getIndex(backboneAtomTypes::AtomC), Acceptor.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];
+                       atomH[i] = prevCO / prevCODist;
+                   }
+                   distanceHO = CalculateAtomicDistances(Donor.getIndex(backboneAtomTypes::AtomO), atomH, fr, pbc);
+                   distanceHC = CalculateAtomicDistances(Donor.getIndex(backboneAtomTypes::AtomC), atomH, 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 (Donor.getIndex(backboneAtomTypes::AtomC) && Donor.getIndex(backboneAtomTypes::AtomO)
-           && Acceptor.getIndex(backboneAtomTypes::AtomN) && ( Acceptor.getIndex(backboneAtomTypes::AtomH) || (initParams.mode == modeType::Classique) ) )  // Kinda ew
-       {
            if (CalculateAtomicDistances(
                        Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
                < minimalCAdistance)
            {
-               if ((distanceON < minimalAtomDistance) || (distanceCH < minimalAtomDistance)
-                   || (distanceOH < minimalAtomDistance) || (distanceCN < minimalAtomDistance))
+               if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
+                   || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
                {
                    HbondEnergy = minEnergy;
 
@@ -643,15 +651,17 @@ void DsspTool::calculateHBondEnergy(ResInfo& Donor,
                {
                    HbondEnergy =
                            kCouplingConstant
-                           * ((1 / distanceON) + (1 / distanceCH) - (1 / distanceOH) - (1 / distanceCN));
+                           * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
                    HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
 
+                   std::cout << "Calculated ENERGY - " << HbondEnergy << std::endl;
+
                    if ( HbondEnergy < minEnergy ){
                        HbondEnergy = minEnergy;
                    }
 
 //                   if ( HbondEnergy < -0.5 ){
-//                       std::cout << "HBOND exists cause of energy" << std::endl; // NOT working in original DSSP
+//                       std::cout << "HBOND exists cause of energy" << std::endl;
 //                   }
 
 
@@ -693,6 +703,14 @@ float DsspTool::CalculateAtomicDistances(const int &A, const int &B, const t_trx
    return r.norm(); // TODO * gmx::c_nm2A; if not PDB, i guess????
 }
 
+/* Calculate Distance From B to A, where A is only fake coordinates */
+float DsspTool::CalculateAtomicDistances(const int &A, const rvec &B, const t_trxframe &fr, const t_pbc *pbc)
+{
+   gmx::RVec r{ 0, 0, 0 };
+   pbc_dx(pbc, fr.x[A], B, r.as_vec());
+   return r.norm(); // TODO * gmx::c_nm2A; if not PDB, i guess????
+}
+
 void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz)
 {
 
@@ -721,6 +739,9 @@ void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const
            if( proLINE.compare("PRO") == 0 ){
                IndexMap[i].is_proline = true;
            }
+           IndexMap[i - 1].nextResi = &IndexMap[i];
+
+           IndexMap[i].prevResi = &IndexMap[i - 1];
        }
        std::string atomname(*(top.atoms()->atomname[static_cast<std::size_t>(*ai)]));
        if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomCA])
@@ -739,7 +760,7 @@ void DsspTool::initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const
        {
            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomN)] = *ai;
        }
-       else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH])
+       else if (atomname == backboneAtomTypeNames[backboneAtomTypes::AtomH] && initParamz.addHydrogens == false) // Юзать водород в структуре
        {
            IndexMap[i]._backboneIndices[static_cast<std::size_t>(backboneAtomTypes::AtomH)] = *ai;
        }

diff --git a/src/dssptools.h b/src/dssptools.h
index 1f1f3a522d912337a60e8a9843f42ba0d453d3b2..08adfafc78dea0dce107d557abae439a5ced75ba 100644 (file)
--- a/src/dssptools.h
+++ b/src/dssptools.h
@@ -89,7 +89,7 @@ struct initParameters {
     real                                                 cutoff_; // = 4.0; ???
     modeType                                             mode;
     NBSearchMethod                                       NBS;
-    bool                                                 verbose, PPHelices;
+    bool                                                 verbose, PPHelices, addHydrogens;
 };
 
 enum class backboneAtomTypes : std::size_t
@@ -107,11 +107,12 @@ const gmx::EnumerationArray<backboneAtomTypes, const char*> backboneAtomTypeName
 };
 
 struct ResInfo {
-    std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)>     _backboneIndices{ 0, 0, 0, 0, 0 }; // todo something with zeroes
+    std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)>     _backboneIndices{ 0, 0, 0, 0, 0 }; // TODO something with zeroes
     std::size_t                                                                     getIndex(backboneAtomTypes atomTypeName) const;
     t_resinfo                                                                       *info{nullptr}, *donor[2]{nullptr, nullptr}, *acceptor[2]{nullptr, nullptr};
+    ResInfo                                                                         *prevResi{nullptr}, *nextResi{nullptr};
     float                                                                           donorEnergy[2]{}, acceptorEnergy[2]{};
-    bool                                                                            is_proline{false}, is_deformed_backbone{false};
+    bool                                                                            is_proline{false};
 };
 
 
@@ -256,7 +257,8 @@ private:
                                                                    ResInfo& Acceptor,
                                                                    const t_trxframe&          fr,
                                                                    const t_pbc*               pbc);
-   float CalculateAtomicDistances(const int& A, const int& B, const t_trxframe& fr, const t_pbc* pbc);
+   float CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
+   float CalculateAtomicDistances(const int &A, const rvec &B, const t_trxframe &fr, const t_pbc *pbc);
 
    class DsspStorage
    {