asdsad
authorMax <Infinity2573@gmail.com>
Sat, 16 Jul 2022 22:43:31 +0000 (01:43 +0300)
committerMax <Infinity2573@gmail.com>
Sat, 16 Jul 2022 22:43:31 +0000 (01:43 +0300)
src/dssp.cpp
src/dssptools.cpp
src/dssptools.h

index 51074df087404d3a7430f64c30921f05a725a5db..c56dd978f698ec9c84e4acbd67a59f6b39f2e913 100644 (file)
@@ -100,7 +100,7 @@ void Dssp::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* s
    options->addOption(
                EnumOption<modeType>("HBond Energy Calc Mode").store(&initParams.mode).required().defaultValue(modeType::Classique).enumValue(modeTypeNames).description("Mode for DSSP"));
    options->addOption(
-               EnumOption<NBSearchMethod>("NBSearchMethod Algo").store(&initParams.NBS).required().defaultValue(NBSearchMethod::Classique).enumValue(NBSearchMethodNames).description("Mode for DSSP"));
+               EnumOption<NBSearchMethod>("Algorithm for search residues' neighbours").store(&initParams.NBS).required().defaultValue(NBSearchMethod::DSSP).enumValue(NBSearchMethodNames).description("Mode for DSSP"));
    options->addOption(
                BooleanOption("v").store(&initParams.verbose).defaultValue(false).description(">:("));
    settings->setHelpText(desc);
index f5762465ddfa9dab695a3d480afe24ee25864007..8e21c3006a3ab8c2c48b656fe97814e0fdcfb8ae 100644 (file)
@@ -45,7 +45,8 @@
 
 
 /* Запилить версию поиска всех остатков против всех. Медленно, но как в дссп.
-   Расхождение может быть из-за тестбонда.
+   Расхождение может разного восприятия доноров и акцепторовю=. Поменять.
+   Пролин отдельное говно (где ассайн гидроген).
 */
 
 #include "dssptools.h"
@@ -79,18 +80,6 @@ std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const
    return _ResInfo[static_cast<std::size_t>(atomTypeName)];
 }
 
-t_resinfo* ResInfo::getInfo() const{
-    return resinfo;
-}
-
-void ResInfo::setHBondedPair(const t_resinfo &PairResinfo){
-    HBondedResidues.push_back(&PairResinfo);
-}
-
-std::vector<const t_resinfo*>   ResInfo::getHBondedPairs() const{
-    return HBondedResidues;
-}
-
 secondaryStructures::secondaryStructures(){
 }
 void secondaryStructures::initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez){
@@ -114,11 +103,11 @@ void secondaryStructures::secondaryStructuresData::setStatus(const HelixPosition
     TurnsStatusArray[static_cast<std::size_t>(turn)] = helixPosition;
 }
 
-bool secondaryStructures::secondaryStructuresData::getStatus(const secondaryStructureTypes secondaryStructureTypeName) const {
+bool secondaryStructures::secondaryStructuresData::getStatus(const secondaryStructureTypes secondaryStructureTypeName) const{
     return SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)];
 }
 
-HelixPositions secondaryStructures::secondaryStructuresData::getStatus(const turnsTypes turn) const {
+HelixPositions secondaryStructures::secondaryStructuresData::getStatus(const turnsTypes turn) const{
     return TurnsStatusArray[static_cast<std::size_t>(turn)];
 }
 
@@ -127,6 +116,7 @@ secondaryStructureTypes secondaryStructures::secondaryStructuresData::getStatus(
 }
 
 bool secondaryStructures::hasHBondBetween(std::size_t resi1, std::size_t resi2) const{
+    if(  ){}
     for(const t_resinfo* &i : (*ResInfoMap)[resi1].getHBondedPairs()){
         if (i == (*ResInfoMap)[resi2].getInfo()){
             return true;
@@ -575,10 +565,10 @@ bool DsspTool::isHbondExist(const ResInfo& resA,
     * Hbond if E < -0.5
     */
 
-   std::string proLINE {*(resA.getInfo()->name)};
+   std::string proLINE {*(resA.resinfo->name)};
 
    if( proLINE.compare("PRO")){
-       std::cerr << "PRO detected! This is Resi№ " << *(resA.getInfo()->name) << std::endl;
+       std::cerr << "PRO-donor detected! This is Resi№ " << *(resA.resinfo->name) << std::endl;
        return false;
    }
    const float kCouplingConstant = 27.888;
@@ -594,21 +584,19 @@ bool DsspTool::isHbondExist(const ResInfo& resA,
            resA.getIndex(backboneAtomTypes::AtomC), resB.getIndex(backboneAtomTypes::AtomN), fr, pbc);
 
    switch(initParams.mode){
-   case modeType::Classique:
+   case modeType::Classique: {
        distanceOH = distanceON;
        distanceCH = distanceCN;
        break;
-   default:
+   }
+   default: {
        distanceOH = CalculateAtomicDistances(
                resA.getIndex(backboneAtomTypes::AtomO), resB.getIndex(backboneAtomTypes::AtomH), fr, pbc);
        distanceCH = CalculateAtomicDistances(
                resA.getIndex(backboneAtomTypes::AtomC), resB.getIndex(backboneAtomTypes::AtomH), fr, pbc);
        break;
    }
-
-   std::cout << "distance ON = " << distanceON << " distance CH = " <<  distanceCH << " distance OH = " << distanceOH << " distance CN = " << distanceCN << std::endl;
-   std::cout << "CACA distance = " << CalculateAtomicDistances(
-                    resA.getIndex(backboneAtomTypes::AtomCA), resB.getIndex(backboneAtomTypes::AtomCA), fr, pbc) << std::endl;
+   }
 
    if (resA.getIndex(backboneAtomTypes::AtomC) && resA.getIndex(backboneAtomTypes::AtomO)
        && resB.getIndex(backboneAtomTypes::AtomN) && ( resB.getIndex(backboneAtomTypes::AtomH) || (initParams.mode == modeType::Classique) ) )
@@ -635,6 +623,8 @@ bool DsspTool::isHbondExist(const ResInfo& resA,
    return false;
 }
 
+
+/* Calculate Distance From B to A */
 float DsspTool::CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc)
 {
    gmx::RVec r{ 0, 0, 0 };
@@ -746,8 +736,10 @@ void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc)
         gmx::AnalysisNeighborhoodPair       pair;
         while (pairSearch.findNextPair(&pair))
         {
+            /* refIndex - Donor, testIndex - Acceptor */
             if(isHbondExist(IndexMap[pair.refIndex()], IndexMap[pair.testIndex()], fr, pbc)){
-                IndexMap[pair.refIndex()].setHBondedPair(*(IndexMap[pair.testIndex()].getInfo()));
+                IndexMap[pair.refIndex()].setAcceptor(*(IndexMap[pair.testIndex()].getInfo()));
+                IndexMap[pair.testIndex()].setDonor(*(IndexMap[pair.refIndex()].getInfo()));
             }
         }
 
@@ -762,18 +754,27 @@ void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc)
         as_.AltPairSearch(fr, IndexMap);
 
         while (as_.findNextPair()){
+            /* getResiI - Donor, getResiJ - Acceptor */
             if(isHbondExist(IndexMap[as_.getResiI()], IndexMap[as_.getResiJ()], fr, pbc)){
-                IndexMap[as_.getResiI()].setHBondedPair(*(IndexMap[as_.getResiJ()].getInfo()));
+                IndexMap[as_.getResiI()].setAcceptor(*(IndexMap[as_.getResiJ()].getInfo()));
+                IndexMap[as_.getResiJ()].setDonor(*(IndexMap[as_.getResiI()].getInfo()));
             }
         }
 
         break;
     }
     default: {
-        for(std::size_t i {0}; i < IndexMap.size(); ++i){
-            for(std::size_t j {0}; j < IndexMap.size() && j != i; ++j){
+        for(std::size_t i {0}; i + 1 < IndexMap.size(); ++i){
+            for(std::size_t j {i + 1}; j < IndexMap.size(); ++j){
+                /* I - Donor, J - Acceptor */
                 if(isHbondExist(IndexMap[i], IndexMap[j], fr, pbc)){
-                    IndexMap[i].setHBondedPair(*(IndexMap[j].getInfo()));
+                    IndexMap[i].setAcceptor(*(IndexMap[j].getInfo()));
+                    IndexMap[j].setDonor(*(IndexMap[i].getInfo()));
+                }
+
+                if (j != i + 1 && isHbondExist(IndexMap[j], IndexMap[i], fr, pbc)){
+                    IndexMap[j].setAcceptor(*(IndexMap[i].getInfo()));
+                    IndexMap[i].setDonor(*(IndexMap[j].getInfo()));
                 }
             }
         }
index d3e984040dde1371b247e722b846dbf8df174681..c8782da1d71902cad5267d2499eb5c28c8bdcd36 100644 (file)
@@ -76,12 +76,12 @@ enum class NBSearchMethod : std::size_t
 {
     Classique = 0,
     Experimental,
-    Stoopid,
+    DSSP,
     Count
 };
 
 const gmx::EnumerationArray<NBSearchMethod, const char*> NBSearchMethodNames = {
-    { "Classique", "Experimental", "Stoopid" }
+    { "Classique", "Experimental", "DSSP" }
 };
 
 struct initParameters {
@@ -113,14 +113,12 @@ public:
    void   setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex);
    void   setInfo(const t_resinfo &info);
    std::size_t getIndex(backboneAtomTypes atomTypeName) const;
-   t_resinfo*   getInfo() const;
-   void   setHBondedPair(const t_resinfo &PairResinfo);
-   std::vector<const t_resinfo*>   getHBondedPairs() const;
+
+   t_resinfo                                                                    *resinfo, *donor[2]{nullptr, nullptr}, *acceptor[2]{nullptr, nullptr};
+   float                                                                        donorEnergy[2]{}, acceptorEnergy[2]{};
 
 private:
    std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)>  _ResInfo{ 0, 0, 0, 0, 0 };
-   t_resinfo                                                              *resinfo;
-   std::vector<const t_resinfo*>                                                      HBondedResidues;
 };
 
 enum class secondaryStructureTypes : std::size_t {
@@ -190,7 +188,6 @@ private:
     const gmx::EnumerationArray<secondaryStructureTypes, const char> secondaryStructureTypeNames = {
        { '=', 'S', 'T', 'I', 'G', 'E', 'B', 'H'} // TODO
     };
-    const std::vector<bridgeTypes> Bridges{ bridgeTypes::Bridge, bridgeTypes::AntiBridge};
     std::string     SecondaryStructuresStringLine;
     const std::size_t     strideFactor{3};