src/dssptools.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Declares trajectory analysis module for secondary structure asignment.
  38  *
  39  * \author Sergey Gorelov <Infinity2573@gmail.com>
  40  * \author Anatoly Titov <titov_ai@pnpi.nrcki.ru>
  41  * \author Alexey Shvetsov <alexxyum@gmail.com>
  42  * \ingroup module_trajectoryanalysis
  43  */
  44
  45 #ifndef GMX_TRAJECTORYANALYSIS_MODULES_DSSPTOOLS_H
  46 #define GMX_TRAJECTORYANALYSIS_MODULES_DSSPTOOLS_H
  47
  48 #include <algorithm>
  49 #include "gromacs/math/units.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include <gromacs/trajectoryanalysis.h>
  52 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  53 #include <set>
  54 #include <fstream>
  55 #include <mutex>
  56
  57 namespace gmx
  58 {
  59
  60 namespace analysismodules
  61 {
  62
  63 enum class NBSearchMethod : std::size_t
  64 {
  65     Classique = 0,
  66     Experimental,
  67     DSSP,
  68     Count
  69 };
  70
  71 const gmx::EnumerationArray<NBSearchMethod, const char*> NBSearchMethodNames = {
  72     { "Classique", "Experimental", "DSSP" }
  73 };
  74
  75 struct initParameters {
  76     Selection                                            sel_;
  77     real                                                 cutoff_; // = 4.0; ???
  78     NBSearchMethod                                       NBS;
  79     bool                                                 verbose, PPHelices, addHydrogens;
  80     int                                                  pp_stretch;
  81 };
  82
  83 enum class backboneAtomTypes : std::size_t
  84 {
  85    AtomCA,
  86    AtomC,
  87    AtomO,
  88    AtomN,
  89    AtomH,
  90    Count
  91 };
  92 //! String values corresponding to backbone atom types
  93 const gmx::EnumerationArray<backboneAtomTypes, const char*> backboneAtomTypeNames = {
  94    { "CA", "C", "O", "N", "H" }
  95 };
  96
  97 struct ResInfo {
  98     std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)>     _backboneIndices{ 0, 0, 0, 0, 0 }; // TODO something with zeroes
  99     std::size_t                                                                     getIndex(backboneAtomTypes atomTypeName) const;
 100     t_resinfo                                                                       *info{nullptr}, *donor[2]{nullptr, nullptr}, *acceptor[2]{nullptr, nullptr};
 101     ResInfo                                                                         *prevResi{nullptr}, *nextResi{nullptr};
 102     double                                                                           donorEnergy[2]{}, acceptorEnergy[2]{};
 103     bool                                                                            is_proline{false};
 104 };
 105
 106
 107 //class ResInfo
 108 //{
 109 //public:
 110 //   void   setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex);
 111 //   std::size_t getIndex(backboneAtomTypes atomTypeName) const;
 112
 113 //   ResInfo                                                                    *donor[2]{nullptr, nullptr}, *acceptor[2]{nullptr, nullptr};
 114 //   float                                                                        donorEnergy[2]{}, acceptorEnergy[2]{};
 115 //   std::string                                                                  resiName;
 116
 117 //private:
 118 //   std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)>  _ResInfo{ 0, 0, 0, 0, 0 };
 119 //};
 120
 121 enum class secondaryStructureTypes : std::size_t { // TODO
 122     Loop = 0, // ~
 123     Break, // =
 124     Bend, // S
 125     Turn, // T
 126     Helix_PP, // P
 127     Helix_5, // I
 128     Helix_3, // G
 129     Strand, // E
 130     Bridge, // B
 131     Helix_4, // H
 132     Count
 133
 134 };
 135
 136 enum class turnsTypes : std::size_t {
 137     Turn_3 = 0,
 138     Turn_4,
 139     Turn_5,
 140     Turn_PP,
 141     Count
 142 };
 143
 144 enum class HelixPositions : std::size_t {
 145     None = 0,
 146     Start,
 147     Middle,
 148     End,
 149     Start_AND_End,
 150     Count
 151 };
 152
 153 enum class bridgeTypes : std::size_t {
 154     None = 0,
 155     AntiParallelBridge,
 156     ParallelBridge,
 157     Count
 158 };
 159
 160 class secondaryStructures{ // PatterSearch Wrapper
 161 public:
 162     secondaryStructures();
 163     void                    initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez, const int _pp_stretch);
 164     std::string             patternSearch();
 165     ~secondaryStructures();
 166
 167     class secondaryStructuresData{ // PatternSearch Tool
 168     public:
 169         void                                setStatus(const secondaryStructureTypes secondaryStructureTypeName, bool status = true);
 170         void                                setStatus(const HelixPositions helixPosition, const turnsTypes turn);
 171         bool                                getStatus(const secondaryStructureTypes secondaryStructureTypeName) const, isBreakPartnerWith(const secondaryStructuresData *partner) const;
 172         HelixPositions                      getStatus(const turnsTypes turn) const;
 173         secondaryStructureTypes             getStatus() const;
 174         void                                setBreak(secondaryStructuresData *breakPartner), setBridge(secondaryStructuresData *bridgePartner, std::size_t bridgePartnerIndex, bridgeTypes bridgeType);
 175         bool                                hasBridges() const, hasBridges(bridgeTypes bridgeType) const, isBridgePartnerWith(secondaryStructuresData *bridgePartner, bridgeTypes bridgeType) const;
 176         std::size_t                         getBridgePartnerIndex(bridgeTypes bridgeType) const;
 177         secondaryStructuresData             getBridgePartner(bridgeTypes bridgeType) const;
 178         void                                set_phi(const ResInfo resi), set_psi(const ResInfo resi);
 179     private:
 180         std::array<bool, static_cast<std::size_t>(secondaryStructureTypes::Count)>  SecondaryStructuresStatusArray{ true, 0, 0, 0, 0, 0, 0 };
 181         secondaryStructuresData                                        *breakPartners[2]{nullptr, nullptr};
 182         secondaryStructuresData                                        *bridgePartners[2]{nullptr, nullptr};
 183         std::size_t                                                     bridgePartnersIndexes[2]{};
 184         secondaryStructureTypes                                         SecondaryStructuresStatus {secondaryStructureTypes::Loop};
 185         std::array<HelixPositions, static_cast<std::size_t>(turnsTypes::Count)>  TurnsStatusArray {HelixPositions::None, HelixPositions::None, HelixPositions::None, HelixPositions::None};
 186     };
 187
 188     std::vector<secondaryStructuresData>     SecondaryStructuresStatusMap;
 189 private:
 190
 191     const std::vector<ResInfo>                                       *ResInfoMap;
 192
 193     const gmx::EnumerationArray<secondaryStructureTypes, const char> secondaryStructureTypeNames = {
 194        { '~', '=', 'S', 'T', 'P', 'I', 'G', 'E', 'B', 'H'}
 195     };
 196     std::string     SecondaryStructuresStringLine;
 197
 198     bool            hasHBondBetween(std::size_t resi1, std::size_t resi2) const;
 199
 200     bool            NoChainBreaksBetween(std::size_t ResiStart, std::size_t ResiEnd) const, isLoop(const std::size_t resiIndex) const, PiHelixPreference;
 201     bridgeTypes     calculateBridge(std::size_t i, std::size_t j) const;
 202     void            analyzeBridgesAndStrandsPatterns(), analyzeTurnsAndHelicesPatterns(), analyzePPHelicesPatterns();
 203
 204     const float                           HBondEnergyCutOff{ -0.5 };
 205 };
 206
 207 class alternateNeighborhoodSearch{
 208 public:
 209     alternateNeighborhoodSearch(){}
 210     void setCutoff(const real &cutoff_init);
 211     void AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo>      &IndexMap);
 212     bool findNextPair();
 213     std::size_t getResiI() const, getResiJ() const;
 214 private:
 215     real                                                 cutoff = 4.0;
 216     std::size_t resiIpos{0}, resiJpos{0};
 217     std::vector<std::vector<bool>> PairMap;
 218     std::vector<std::size_t> MiniBoxMap, partners, num_of_miniboxes;
 219     std::vector<std::vector<std::size_t>> MiniBoxesReverseMap;
 220     std::vector<std::vector<std::vector<std::vector<std::size_t>>>> MiniBoxesMap;
 221     std::array<int, 3> MiniBoxSize{0, 0, 0}, offset{-1, -1, -1}, fixBox{0,0,0};
 222     void GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo>      &IndexMap);
 223     void FixAtomCoordinates(real &coordinate, const real vector_length);
 224     void ReCalculatePBC(int &x, const int &x_max);
 225     std::vector<std::vector<bool>>::iterator ResiI;
 226     std::vector<bool>::iterator ResiJ;
 227     bool init = true;
 228 };
 229
 230 class DsspTool
 231 {
 232 public:
 233    DsspTool();
 234    void initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz);
 235
 236    void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc);
 237
 238    int                                   nres{};
 239
 240    std::vector<std::pair<int, std::string>> getData();
 241 private:
 242    initParameters                        initParams;
 243    AtomsDataPtr                          atoms_;
 244    std::string                           filename_;
 245    void                                  calculateBends(const t_trxframe& fr, const t_pbc* pbc), calculateDihedrals(const t_trxframe& fr, const t_pbc* pbc);
 246    std::vector<std::string>              ResiNames;
 247    std::vector<std::size_t>              AtomResi;
 248    std::vector<ResInfo>                  IndexMap;
 249    std::vector<std::vector<std::size_t>> nturnsmap, Bridges, AntiBridges;
 250    secondaryStructures                   PatternSearch;
 251    const float                           HBondEnergyCutOff{ -0.5 },
 252                                          minimalCAdistance{ 9.0 };
 253    void                                  calculateHBondEnergy(ResInfo& Donor,
 254                                                                    ResInfo& Acceptor,
 255                                                                    const t_trxframe&          fr,
 256                                                                    const t_pbc*               pbc);
 257    float CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
 258    float CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
 259    float CalculateDihedralAngle(const int &A, const int &B, const int &C, const int &D, const t_trxframe &fr, const t_pbc *pbc);
 260
 261    class DsspStorage
 262    {
 263    public:
 264        DsspStorage();
 265
 266        /*Storages our dirty data, duh*/
 267        void storageData(int frnr, std::string data);
 268
 269        /*Clear data after usage for the next pdb file*/
 270        void clearAll();
 271
 272        /*Perform several checks and fucks and shits and returns sorted sexy data*/
 273        std::vector<std::pair<int, std::string>> returnData();
 274    private:
 275        std::vector<std::pair<int, std::string>> storaged_data;
 276        static std::mutex mx;
 277    };
 278
 279    static DsspStorage                           Storage;
 280 };
 281
 282 }
 283 //end of namespace analysismodules
 284
 285 }
 286 //end of namespace gmx
 287
 288
 289 #endif