2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Declares trajectory analysis module for secondary structure asignment.
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
45 #ifndef GMX_TRAJECTORYANALYSIS_MODULES_DSSPTOOLS_H
46 #define GMX_TRAJECTORYANALYSIS_MODULES_DSSPTOOLS_H
49 #include "gromacs/math/units.h"
50 #include "gromacs/pbcutil/pbc.h"
51 #include <gromacs/trajectoryanalysis.h>
52 #include "gromacs/trajectoryanalysis/topologyinformation.h"
60 namespace analysismodules
63 enum class NBSearchMethod : std::size_t
71 const gmx::EnumerationArray<NBSearchMethod, const char*> NBSearchMethodNames = {
72 { "Classique", "Experimental", "DSSP" }
75 struct initParameters {
77 real cutoff_; // = 4.0; ???
79 bool verbose, PPHelices, addHydrogens;
82 enum class backboneAtomTypes : std::size_t
91 //! String values corresponding to backbone atom types
92 const gmx::EnumerationArray<backboneAtomTypes, const char*> backboneAtomTypeNames = {
93 { "CA", "C", "O", "N", "H" }
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};
109 // void setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex);
110 // std::size_t getIndex(backboneAtomTypes atomTypeName) const;
112 // ResInfo *donor[2]{nullptr, nullptr}, *acceptor[2]{nullptr, nullptr};
113 // float donorEnergy[2]{}, acceptorEnergy[2]{};
114 // std::string resiName;
117 // std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)> _ResInfo{ 0, 0, 0, 0, 0 };
120 enum class secondaryStructureTypes : std::size_t { // TODO
135 enum class turnsTypes : std::size_t {
142 enum class HelixPositions : std::size_t {
151 enum class bridgeTypes : std::size_t {
158 class secondaryStructures{ // PatterSearch Wrapper
160 secondaryStructures();
161 void initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez = true);
162 std::string patternSearch();
163 ~secondaryStructures();
165 class secondaryStructuresData{ // PatternSearch Tool
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, bridgeTypes bridgeType);
173 bool hasBridges() const, isBridgePartnerWith(secondaryStructuresData *bridgePartner, bridgeTypes bridgeType) const;
175 std::array<bool, static_cast<std::size_t>(secondaryStructureTypes::Count)> SecondaryStructuresStatusArray{ true, 0, 0, 0, 0, 0, 0 };
176 secondaryStructuresData *breakPartners[2]{nullptr, nullptr};
177 secondaryStructuresData *bridgePartners[2]{nullptr, nullptr};
178 secondaryStructureTypes SecondaryStructuresStatus {secondaryStructureTypes::Loop};
179 std::array<HelixPositions, static_cast<std::size_t>(turnsTypes::Count)> TurnsStatusArray {HelixPositions::None, HelixPositions::None, HelixPositions::None};
182 std::vector<secondaryStructuresData> SecondaryStructuresStatusMap;
185 const std::vector<ResInfo> *ResInfoMap;
187 const gmx::EnumerationArray<secondaryStructureTypes, const char> secondaryStructureTypeNames = {
188 { '~', '=', 'S', 'T', 'P', 'I', 'G', 'E', 'B', 'H'}
190 std::string SecondaryStructuresStringLine;
192 bool hasHBondBetween(std::size_t resi1, std::size_t resi2) const;
194 bool NoChainBreaksBetween(std::size_t ResiStart, std::size_t ResiEnd) const, isLoop(const std::size_t resiIndex) const, PiHelixPreference;
195 bridgeTypes calculateBridge(std::size_t i, std::size_t j) const;
196 void analyzeBridgesAndLaddersPatterns(), analyzeTurnsAndHelicesPatterns(), analyzePPHelicesPatterns();
198 const float HBondEnergyCutOff{ -0.5 };
201 class alternateNeighborhoodSearch{
203 alternateNeighborhoodSearch(){}
204 void setCutoff(const real &cutoff_init);
205 void AltPairSearch(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap);
207 std::size_t getResiI() const, getResiJ() const;
210 std::size_t resiIpos{0}, resiJpos{0};
211 std::vector<std::vector<bool>> PairMap;
212 std::vector<std::size_t> MiniBoxMap, partners, num_of_miniboxes;
213 std::vector<std::vector<std::size_t>> MiniBoxesReverseMap;
214 std::vector<std::vector<std::vector<std::vector<std::size_t>>>> MiniBoxesMap;
215 std::array<int, 3> MiniBoxSize{0, 0, 0}, offset{-1, -1, -1}, fixBox{0,0,0};
216 void GetMiniBoxesMap(const t_trxframe &fr, const std::vector<ResInfo> &IndexMap);
217 void FixAtomCoordinates(real &coordinate, const real vector_length);
218 void ReCalculatePBC(int &x, const int &x_max);
219 std::vector<std::vector<bool>>::iterator ResiI;
220 std::vector<bool>::iterator ResiJ;
228 void initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top, const initParameters &initParamz);
230 void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc);
234 std::vector<std::pair<int, std::string>> getData();
236 initParameters initParams;
238 std::string filename_;
239 void calculateBends(const t_trxframe& fr, const t_pbc* pbc);
240 std::vector<std::string> ResiNames;
241 std::vector<std::size_t> AtomResi;
242 std::vector<ResInfo> IndexMap;
243 std::vector<std::vector<std::size_t>> nturnsmap, Bridges, AntiBridges;
244 secondaryStructures PatternSearch;
245 const float HBondEnergyCutOff{ -0.5 },
246 minimalCAdistance{ 9.0 };
247 void calculateHBondEnergy(ResInfo& Donor,
249 const t_trxframe& fr,
251 float CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
252 float CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc);
259 /*Storages our dirty data, duh*/
260 void storageData(int frnr, std::string data);
262 /*Clear data after usage for the next pdb file*/
265 /*Perform several checks and fucks and shits and returns sorted sexy data*/
266 std::vector<std::pair<int, std::string>> returnData();
268 std::vector<std::pair<int, std::string>> storaged_data;
269 static std::mutex mx;
272 static DsspStorage Storage;
276 //end of namespace analysismodules
279 //end of namespace gmx