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 <gorelov_sv@pnpi.nrcki.ru>
  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 };
  81
  82 enum class backboneAtomTypes : std::size_t
  83 {
  84    AtomCA,
  85    AtomC,
  86    AtomO,
  87    AtomN,
  88    AtomH,
  89    Count
  90 };
  91 //! String values corresponding to backbone atom types
  92 const gmx::EnumerationArray<backboneAtomTypes, const char*> backboneAtomTypeNames = {
  93    { "CA", "C", "O", "N", "H" }
  94 };
  95
  96 struct ResInfo {
  97     std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)>     _backboneIndices{ 0, 0, 0, 0, 0 }; // TODO something with zeroes
  98     std::size_t                                                                     getIndex(backboneAtomTypes atomTypeName) const;
  99     t_resinfo                                                                       *info{nullptr}, *donor[2]{nullptr, nullptr}, *acceptor[2]{nullptr, nullptr};
 100     ResInfo                                                                         *prevResi{nullptr}, *nextResi{nullptr};
 101     float                                                                           donorEnergy[2]{}, acceptorEnergy[2]{};
 102     bool                                                                            is_proline{false};
 103 };
 104
 105
 106 //class ResInfo
 107 //{
 108 //public:
 109 //   void   setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex);
 110 //   std::size_t getIndex(backboneAtomTypes atomTypeName) const;
 111
 112 //   ResInfo                                                                    *donor[2]{nullptr, nullptr}, *acceptor[2]{nullptr, nullptr};
 113 //   float                                                                        donorEnergy[2]{}, acceptorEnergy[2]{};
 114 //   std::string                                                                  resiName;
 115
 116 //private:
 117 //   std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)>  _ResInfo{ 0, 0, 0, 0, 0 };
 118 //};
 119
 120 enum class secondaryStructureTypes : std::size_t { // TODO
 121     Loop = 0, // ~
 122     Break, // =
 123     Bend, // S
 124     Turn, // T
 125     Helis_PPII, // TODO
 126     Helix_5, // I
 127     Helix_3, // G
 128     Bulge, // B
 129     Ladder, // E
 130     Helix_4, // H
 131     Count
 132
 133 };
 134
 135 enum class turnsTypes : std::size_t {
 136     Turn_3 = 0,
 137     Turn_4,
 138     Turn_5,
 139     Count
 140 };
 141
 142 enum class HelixPositions : std::size_t {
 143     None = 0,
 144     Start,
 145     Middle,
 146     End,
 147     Start_AND_End,
 148     Count
 149 };
 150
 151 enum class bridgeTypes : std::size_t {
 152     None = 0,
 153     AntiParallelBridge,
 154     ParallelBridge,
 155     Count
 156 };
 157
 158 class secondaryStructures{ // PatterSearch Wrapper
 159 public:
 160     secondaryStructures();
 161     void                    initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez = true);
 162     std::string             patternSearch();
 163     ~secondaryStructures();
 164
 165     class secondaryStructuresData{ // PatternSearch Tool
 166     public:
 167         void                                setStatus(const secondaryStructureTypes secondaryStructureTypeName);
 168         void                                setStatus(const HelixPositions helixPosition, const turnsTypes turn);
 169         bool                                getStatus(const secondaryStructureTypes secondaryStructureTypeName) const, isBreakPartnerWith(const secondaryStructuresData *partner) const;
 170         HelixPositions                      getStatus(const turnsTypes turn) const;
 171         secondaryStructureTypes             getStatus() const;
 172         void                                setBreak(secondaryStructuresData *breakPartner);
 173     private:
 174         std::array<bool, static_cast<std::size_t>(secondaryStructureTypes::Count)>  SecondaryStructuresStatusArray{ true, 0, 0, 0, 0, 0, 0 };
 175         secondaryStructuresData                                        *breakPartners[2]{nullptr, nullptr};
 176         secondaryStructureTypes                                         SecondaryStructuresStatus {secondaryStructureTypes::Loop};
 177         std::array<HelixPositions, static_cast<std::size_t>(turnsTypes::Count)>  TurnsStatusArray {HelixPositions::None, HelixPositions::None, HelixPositions::None};
 178         std::array<bool, static_cast<std::size_t>(bridgeTypes::Count)>  BridgesStatusArray{ 0, 0 };
 179     };
 180
 181     std::vector<secondaryStructuresData>     SecondaryStructuresStatusMap;
 182 private:
 183
 184     const std::vector<ResInfo>                                       *ResInfoMap;
 185
 186     const gmx::EnumerationArray<secondaryStructureTypes, const char> secondaryStructureTypeNames = {
 187        { '~', '=', 'S', 'T', '?', 'I', 'G', 'E', 'B', 'H'} // TODO
 188     };
 189     std::string     SecondaryStructuresStringLine;
 190
 191     bool            hasHBondBetween(std::size_t resi1, std::size_t resi2) const;
 192
 193     bool            NoChainBreaksBetween(std::size_t ResiStart, std::size_t ResiEnd) const, isLoop(const std::size_t resiIndex) const, PiHelixPreference;
 194     bridgeTypes     calculateBridge(std::size_t i, std::size_t j) const;
 195     void            analyzeBridgesAndLaddersPatterns(), analyzeTurnsAndHelicesPatterns(), analyzePPHelicesPatterns();
 196
 197     const float                           HBondEnergyCutOff{ -0.5 };
 198 };
 199
 200 class alternateNeighborhoodSearch{
 201 public:
 202     alternateNeighborhoodSearch(){}
 203     void setCutoff(const real &cutoff_init);
 204     void AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo>      &IndexMap);
 205     bool findNextPair();
 206     std::size_t getResiI() const, getResiJ() const;
 207 private:
 208     real                                                 cutoff = 4.0;
 209     std::size_t resiIpos{0}, resiJpos{0};
 210     std::vector<std::vector<bool>> PairMap;
 211     std::vector<std::size_t> MiniBoxMap, partners, num_of_miniboxes;
 212     std::vector<std::vector<std::size_t>> MiniBoxesReverseMap;
 213     std::vector<std::vector<std::vector<std::vector<std::size_t>>>> MiniBoxesMap;
 214     std::array<int, 3> MiniBoxSize{0, 0, 0}, offset{-1, -1, -1}, fixBox{0,0,0};
 215     void GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo>      &IndexMap);
 216     void FixAtomCoordinates(real &coordinate, const real vector_length);
 217     void ReCalculatePBC(int &x, const int &x_max);
 218     std::vector<std::vector<bool>>::iterator ResiI;
 219     std::vector<bool>::iterator ResiJ;
 220     bool init = true;
 221 };
 222
 223 class DsspTool
 224 {
 225 public:
 226    DsspTool();
 227    void initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz);
 228
 229    void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc);
 230
 231    int                                   nres{};
 232
 233    std::vector<std::pair<int, std::string>> getData();
 234 private:
 235    initParameters                        initParams;
 236    AtomsDataPtr                          atoms_;
 237    std::string                           filename_;
 238    void                                  calculateBends(const t_trxframe& fr, const t_pbc* pbc);
 239    std::vector<std::string>              ResiNames;
 240    std::vector<std::size_t>              AtomResi;
 241    std::vector<ResInfo>                  IndexMap;
 242    std::vector<std::vector<std::size_t>> nturnsmap, Bridges, AntiBridges;
 243    secondaryStructures                   PatternSearch;
 244    const float                           HBondEnergyCutOff{ -0.5 },
 245                                          minimalCAdistance{ 9.0 };
 246    void                                  calculateHBondEnergy(ResInfo& Donor,
 247                                                                    ResInfo& Acceptor,
 248                                                                    const t_trxframe&          fr,
 249                                                                    const t_pbc*               pbc);
 250    float CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
 251    float CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
 252
 253    class DsspStorage
 254    {
 255    public:
 256        DsspStorage();
 257
 258        /*Storages our dirty data, duh*/
 259        void storageData(int frnr, std::string data);
 260
 261        /*Clear data after usage for the next pdb file*/
 262        void clearAll();
 263
 264        /*Perform several checks and fucks and shits and returns sorted sexy data*/
 265        std::vector<std::pair<int, std::string>> returnData();
 266    private:
 267        std::vector<std::pair<int, std::string>> storaged_data;
 268        static std::mutex mx;
 269    };
 270
 271    static DsspStorage                           Storage;
 272 };
 273
 274 }
 275 //end of namespace analysismodules
 276
 277 }
 278 //end of namespace gmx
 279
 280
 281 #endif