src/spacetimecorr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015, 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  * Implements gmx::analysismodules::Freevolume.
  38  *
  39  * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42
  43 #include <iostream>
  44 #include <vector>
  45 #include <cfloat>
  46
  47 #include <gromacs/trajectoryanalysis.h>
  48 #include <gromacs/trajectoryanalysis/topologyinformation.h>
  49
  50 //#include "/home/titov_ai/Develop/gromacs-original/src/gromacs/trajectoryanalysis/topologyinformation.h"
  51
  52 #include "newfit.h"
  53 #include "corrtype.h"
  54
  55 using namespace gmx;
  56
  57 using gmx::RVec;
  58
  59 /*! \brief
  60  * \ingroup module_trajectoryanalysis
  61  */
  62
  63 class SpaceTimeCorr : public TrajectoryAnalysisModule
  64 {
  65     public:
  66
  67         SpaceTimeCorr();
  68         virtual ~SpaceTimeCorr();
  69
  70         //! Set the options and setting
  71         virtual void initOptions(IOptionsContainer          *options,
  72                                  TrajectoryAnalysisSettings *settings);
  73
  74         //! First routine called by the analysis framework
  75         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
  76         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
  77                                   const TopologyInformation        &top);
  78
  79         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
  80                                          const t_trxframe                   &fr);
  81
  82         //! Call for each frame of the trajectory
  83         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
  84         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
  85                                   TrajectoryAnalysisModuleData *pdata);
  86
  87         //! Last routine called by the analysis framework
  88         // virtual void finishAnalysis(t_pbc *pbc);
  89         virtual void finishAnalysis(int nframes);
  90
  91         //! Routine to write output, that is additional over the built-in
  92         virtual void writeOutput();
  93
  94     private:
  95
  96         SelectionList                                           sel_;
  97         std::vector< RVec >                                     trajectoryFrame;
  98
  99         std::vector< std::vector < std::vector < size_t > > >   sels;
 100         std::vector< int >                                      index;
 101
 102         int                                                     tauJump            {100};  // selectable
 103
 104         int                                                     window             {1000}; // selectable
 105         int                                                     tau                {500};  // selectable
 106         float                                                   crlBorder          {0.5};  // selectable
 107         float                                                   effRad             {1.5};  // selectable
 108         std::string                                             outPutName;                 // selectable
 109         // может лучше в виде enum?
 110         int                                                     mode               {1};    // selectable
 111
 112         correlationType                                         corrs;
 113         // Copy and assign disallowed by base.
 114 };
 115
 116 SpaceTimeCorr::SpaceTimeCorr(): TrajectoryAnalysisModule()
 117 {
 118 }
 119
 120 SpaceTimeCorr::~SpaceTimeCorr()
 121 {
 122 }
 123
 124 void
 125 SpaceTimeCorr::initOptions(IOptionsContainer          *options,
 126                      TrajectoryAnalysisSettings *settings)
 127 {
 128     static const char *const desc[] = {
 129         "[THISMODULE] to be done"
 130         "-f trajectory.xtc -s structure.tpr -selRaDs residuePlusDomains.ndx -sf residuePlusDomains -mode 0/1 -tau 500 -window 1000 -crl 0.5 -efRad 1 -twStep 100 -outPut outputFilename"
 131     };
 132     // Add the descriptive text (program help text) to the options
 133     settings->setHelpText(desc);
 134     // Add option for output file name
 135     //options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 136     //                        .store(&fnNdx_).defaultBasename("domains")
 137     //                        .description("Index file from the domains"));
 138     // Add option for working mode
 139     options->addOption(gmx::IntegerOption("mode")
 140                             .store(&mode)
 141                             .description("default 0 - creation of sigle matrixes' file | 1 - reading prepared correlation file to make graphs and etc."));
 142     // Add option for tau constant
 143     options->addOption(gmx::IntegerOption("tau")
 144                             .store(&tau)
 145                             .description("max number of frames to shift to evaluate correlations | all starts from array [0 ; tau]"));
 146     // Add option for window constant
 147     options->addOption(gmx::IntegerOption("window")
 148                             .store(&window)
 149                             .description("size of window in frames that correlation evaluation based on"));
 150     // Add option for crl_border constant
 151     options->addOption(FloatOption("crl")
 152                             .store(&crlBorder)
 153                             .description("make graphs based on corrs > crl const"));
 154     // Add option for effRad constant
 155     options->addOption(FloatOption("efRad")
 156                             .store(&effRad)
 157                             .description("effective radius for atoms to see each other to evaluate corrs"));
 158     // Add option for selection list
 159     options->addOption(SelectionOption("sel")
 160                             .storeVector(&sel_)
 161                             .required().dynamicMask().multiValue()
 162                             .description("select residue and domains to form a rigid skeleton"));
 163     // Add option for output file names
 164     options->addOption(StringOption("outPut")
 165                             .store(&outPutName)
 166                             .description("<your name here> + <local file tag>.txt"));
 167     // Add option for time step between windows' starts
 168     options->addOption(gmx::IntegerOption("twStep")
 169                             .store(&tauJump)
 170                             .description("time step between windows' starting positions."));
 171     // Control input settings
 172     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 173     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 174     settings->setPBC(true);
 175 }
 176
 177 // -f trajectory.xtc -s structure.tpr -selRaDs residuePlusDomains.ndx -sf residuePlusDomains -mode 0/1 -tau 500 -window 1000 -crl 0.5 -efRad 1 -twStep 100 -outPut outputFilename
 178
 179 void
 180 SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
 181 {
 182     std::cout << "initAnalysis - start\n";
 183     std::vector< std::string >                                  resNms;
 184     std::vector< RVec >                                         reference;
 185     ArrayRef< const int > atomind {sel_.front().atomIndices()};
 186     index.resize(0);
 187     for (ArrayRef< const int >::iterator ai {atomind.begin()}; (ai < atomind.end()); ai++) {
 188         index.push_back(*ai);
 189         resNms.push_back(*(top.atoms()->resinfo[top.atoms()->atom[*ai].resind].name) + std::to_string((top.atoms()->atom[*ai].resind + 1)));
 190     }
 191     std::string str(sel_.back().name());
 192     size_t tempSize {std::stoul(str.substr(13, 6)) / static_cast< size_t >(tauJump) + 1};
 193     sels.resize(tempSize);
 194     std::vector< size_t > a;
 195     a.resize(0);
 196     for (size_t i {1}; i < sel_.size(); i++) {
 197         std::string str(sel_[i].name());
 198         atomind = sel_[i].atomIndices();
 199         sels[std::stoul(str.substr(13, 6)) / static_cast< size_t >(tauJump)].push_back(a);
 200         for (ArrayRef< const int >::iterator ai {atomind.begin()}; (ai < atomind.end()); ai++) {
 201             sels[std::stoul(str.substr(13, 6)) / static_cast< size_t >(tauJump)].back().push_back(static_cast< size_t >(*ai));
 202         }
 203     }
 204     reference.resize(0);
 205     if (top.hasFullTopology()) {
 206         for (const auto &i : index) {
 207             reference.push_back(top.x().at(static_cast< size_t >(i)));
 208         }
 209     }
 210     std::cout << "maybe ???\n";
 211     corrs.setDefaults(reference, window, tau, tauJump, crlBorder, effRad, mode, outPutName, index, sels, resNms);
 212     std::cout << "initAnalysis - finish\n";
 213 }
 214
 215 void
 216 SpaceTimeCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
 217 {
 218 }
 219
 220 // -s *.tpr -f *.xtc -n *.ndx -sf 'name' -mode <> -tau <> -window <> -crl <> -ef_rad <> -out_put <>
 221
 222 void
 223 SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 224 {
 225     trajectoryFrame.resize(0);
 226     trajectoryFrame.resize(index.size());
 227     for (size_t i {0}; i < index.size(); i++) {
 228         trajectoryFrame[i] = fr.x[index[i]];
 229     }
 230     corrs.update(frnr, trajectoryFrame);
 231     std::cout << "\t" << frnr << std::endl;
 232 }
 233
 234 void
 235 SpaceTimeCorr::finishAnalysis(int nframes)
 236 {
 237     std::vector< std::vector< std::vector< std::pair< size_t, size_t > > > > rout_new;
 238     size_t k {500};
 239     if (tau > 0) {
 240         k = static_cast< size_t >(tau);
 241     }
 242     if (window == 0) {
 243         window = static_cast< int >(k * 2);
 244     }
 245     corrs.readEval();
 246     corrs.printData();
 247     std::cout << "Finish Analysis - end\n\n";
 248 }
 249
 250 void
 251 SpaceTimeCorr::writeOutput()
 252 {
 253
 254 }
 255
 256 /*! \brief
 257  * The main function for the analysis template.
 258  */
 259 int
 260 main(int argc, char *argv[])
 261 {
 262     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SpaceTimeCorr>(argc, argv);
 263 }