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, // P
 126     Helix_5, // I
 127     Helix_3, // G
 128     Bridge, // B
 129     Strand, // 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), setBridge(secondaryStructuresData *bridgePartner, std::size_t bridgePartnerIndex, bridgeTypes bridgeType);
 173         bool                                hasBridges() const, hasBridges(bridgeTypes bridgeType) const, isBridgePartnerWith(secondaryStructuresData *bridgePartner, bridgeTypes bridgeType) const;
 174         std::size_t                         getBridgePartnerIndex(bridgeTypes bridgeType) const;
 175         secondaryStructuresData             getBridgePartner(bridgeTypes bridgeType) const;
 176     private:
 177         std::array<bool, static_cast<std::size_t>(secondaryStructureTypes::Count)>  SecondaryStructuresStatusArray{ true, 0, 0, 0, 0, 0, 0 };
 178         secondaryStructuresData                                        *breakPartners[2]{nullptr, nullptr};
 179         secondaryStructuresData                                        *bridgePartners[2]{nullptr, nullptr};
 180         std::size_t                                                     bridgePartnersIndexes[2]{};
 181         secondaryStructureTypes                                         SecondaryStructuresStatus {secondaryStructureTypes::Loop};
 182         std::array<HelixPositions, static_cast<std::size_t>(turnsTypes::Count)>  TurnsStatusArray {HelixPositions::None, HelixPositions::None, HelixPositions::None};
 183     };
 184
 185     std::vector<secondaryStructuresData>     SecondaryStructuresStatusMap;
 186 private:
 187
 188     const std::vector<ResInfo>                                       *ResInfoMap;
 189
 190     const gmx::EnumerationArray<secondaryStructureTypes, const char> secondaryStructureTypeNames = {
 191        { '~', '=', 'S', 'T', 'P', 'I', 'G', 'B', 'E', 'H'}
 192     };
 193     std::string     SecondaryStructuresStringLine;
 194
 195     bool            hasHBondBetween(std::size_t resi1, std::size_t resi2) const;
 196
 197     bool            NoChainBreaksBetween(std::size_t ResiStart, std::size_t ResiEnd) const, isLoop(const std::size_t resiIndex) const, PiHelixPreference;
 198     bridgeTypes     calculateBridge(std::size_t i, std::size_t j) const;
 199     void            analyzeBridgesAndStrandsPatterns(), analyzeTurnsAndHelicesPatterns(), analyzePPHelicesPatterns();
 200
 201     const float                           HBondEnergyCutOff{ -0.5 };
 202 };
 203
 204 class alternateNeighborhoodSearch{
 205 public:
 206     alternateNeighborhoodSearch(){}
 207     void setCutoff(const real &cutoff_init);
 208     void AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo>      &IndexMap);
 209     bool findNextPair();
 210     std::size_t getResiI() const, getResiJ() const;
 211 private:
 212     real                                                 cutoff = 4.0;
 213     std::size_t resiIpos{0}, resiJpos{0};
 214     std::vector<std::vector<bool>> PairMap;
 215     std::vector<std::size_t> MiniBoxMap, partners, num_of_miniboxes;
 216     std::vector<std::vector<std::size_t>> MiniBoxesReverseMap;
 217     std::vector<std::vector<std::vector<std::vector<std::size_t>>>> MiniBoxesMap;
 218     std::array<int, 3> MiniBoxSize{0, 0, 0}, offset{-1, -1, -1}, fixBox{0,0,0};
 219     void GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo>      &IndexMap);
 220     void FixAtomCoordinates(real &coordinate, const real vector_length);
 221     void ReCalculatePBC(int &x, const int &x_max);
 222     std::vector<std::vector<bool>>::iterator ResiI;
 223     std::vector<bool>::iterator ResiJ;
 224     bool init = true;
 225 };
 226
 227 class DsspTool
 228 {
 229 public:
 230    DsspTool();
 231    void initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz);
 232
 233    void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc);
 234
 235    int                                   nres{};
 236
 237    std::vector<std::pair<int, std::string>> getData();
 238 private:
 239    initParameters                        initParams;
 240    AtomsDataPtr                          atoms_;
 241    std::string                           filename_;
 242    void                                  calculateBends(const t_trxframe& fr, const t_pbc* pbc);
 243    std::vector<std::string>              ResiNames;
 244    std::vector<std::size_t>              AtomResi;
 245    std::vector<ResInfo>                  IndexMap;
 246    std::vector<std::vector<std::size_t>> nturnsmap, Bridges, AntiBridges;
 247    secondaryStructures                   PatternSearch;
 248    const float                           HBondEnergyCutOff{ -0.5 },
 249                                          minimalCAdistance{ 9.0 };
 250    void                                  calculateHBondEnergy(ResInfo& Donor,
 251                                                                    ResInfo& Acceptor,
 252                                                                    const t_trxframe&          fr,
 253                                                                    const t_pbc*               pbc);
 254    float CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
 255    float CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
 256
 257    class DsspStorage
 258    {
 259    public:
 260        DsspStorage();
 261
 262        /*Storages our dirty data, duh*/
 263        void storageData(int frnr, std::string data);
 264
 265        /*Clear data after usage for the next pdb file*/
 266        void clearAll();
 267
 268        /*Perform several checks and fucks and shits and returns sorted sexy data*/
 269        std::vector<std::pair<int, std::string>> returnData();
 270    private:
 271        std::vector<std::pair<int, std::string>> storaged_data;
 272        static std::mutex mx;
 273    };
 274
 275    static DsspStorage                           Storage;
 276 };
 277
 278 }
 279 //end of namespace analysismodules
 280
 281 }
 282 //end of namespace gmx
 283
 284
 285 #endif