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 < unsigned int > > > sels;
 100         std::vector< int >                                          index;
 101
 102
 103         int                                                         frames              = 0;
 104         int                                                         tau_jump            = 100;
 105
 106         int                                                         window              = 1000; // selectable
 107         int                                                         tau                 = 500;  // selectable
 108         float                                                       crl_border          = 0.5;  // selectable
 109         float                                                       eff_rad             = 1.5;  // selectable
 110         std::string                                                 OutPutName;                 // selectable
 111         int                                                         mode                = 0;    // selectable
 112
 113         correlationType                                             corrs;
 114         // Copy and assign disallowed by base.
 115 };
 116
 117 SpaceTimeCorr::SpaceTimeCorr(): TrajectoryAnalysisModule()
 118 {
 119 }
 120
 121 SpaceTimeCorr::~SpaceTimeCorr()
 122 {
 123 }
 124
 125 void
 126 SpaceTimeCorr::initOptions(IOptionsContainer          *options,
 127                      TrajectoryAnalysisSettings *settings)
 128 {
 129     static const char *const desc[] = {
 130         "[THISMODULE] to be done"
 131         "-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"
 132     };
 133     // Add the descriptive text (program help text) to the options
 134     settings->setHelpText(desc);
 135     // Add option for output file name
 136     //options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 137     //                        .store(&fnNdx_).defaultBasename("domains")
 138     //                        .description("Index file from the domains"));
 139     // Add option for working mode
 140     options->addOption(gmx::IntegerOption("mode")
 141                             .store(&mode)
 142                             .description("default 0 - creation of sigle matrixes' file | 1 - reading prepared correlation file to make graphs and etc."));
 143     // Add option for tau constant
 144     options->addOption(gmx::IntegerOption("tau")
 145                             .store(&tau)
 146                             .description("max number of frames to shift to evaluate correlations | all starts from array [0 ; tau]"));
 147     // Add option for window constant
 148     options->addOption(gmx::IntegerOption("window")
 149                             .store(&window)
 150                             .description("size of window in frames that correlation evaluation based on"));
 151     // Add option for crl_border constant
 152     options->addOption(FloatOption("crl")
 153                             .store(&crl_border)
 154                             .description("make graphs based on corrs > crl const"));
 155     // Add option for eff_rad constant
 156     options->addOption(FloatOption("efRad")
 157                             .store(&eff_rad)
 158                             .description("effective radius for atoms to see each other to evaluate corrs"));
 159     // Add option for selection list
 160     options->addOption(SelectionOption("selRaDs").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(&tau_jump)
 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
 183     std::vector< std::string >                                  resNms;
 184     std::vector< RVec >                                         reference;
 185
 186     ArrayRef< const int > atomind = sel_.front().atomIndices();
 187     index.resize(0);
 188     for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 189         index.push_back(*ai);
 190         resNms.push_back(*(top.atoms()->resinfo[top.atoms()->atom[*ai].resind].name) + std::to_string((top.atoms()->atom[*ai].resind + 1)));
 191     }
 192
 193     std::string str(sel_.back().name());
 194     unsigned long tempSize = std::stoul(str.substr(13, 6)) / static_cast< unsigned long >(tau_jump);
 195
 196     sels.resize(tempSize);
 197
 198     for (unsigned int i = 1; i < sel_.size(); i++) {
 199         std::string str(sel_[i].name());
 200         atomind = sel_[i].atomIndices();
 201         std::vector< unsigned int > a;
 202         a.resize(0);
 203         sels[std::stoul(str.substr(13, 6)) / static_cast< unsigned long >(tau_jump)].push_back(a);
 204         for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 205             sels[std::stoul(str.substr(13, 6)) / static_cast< unsigned long >(tau_jump)].back().push_back(static_cast< unsigned int >(*ai));
 206         }
 207     }
 208
 209     reference.resize(0);
 210     if (top.hasFullTopology()) {
 211         for (unsigned int i = 0; i < index.size(); i++) {
 212             reference.push_back(top.x().at(static_cast< unsigned int >(index[i])));
 213         }
 214     }
 215
 216     corrs.setDefaults(reference, window, tau, tau_jump, crl_border, eff_rad, 0, OutPutName, sels, resNms);
 217 }
 218
 219 void
 220 SpaceTimeCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
 221 {
 222 }
 223
 224 // -s *.tpr -f *.xtc -n *.ndx -sf 'name' -mode <> -tau <> -window <> -crl <> -ef_rad <> -out_put <>
 225
 226 void
 227 SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 228 {
 229     trajectoryFrame.resize(0);
 230     trajectoryFrame.resize(index.size());
 231     for (unsigned int i = 0; i < index.size(); i++) {
 232         trajectoryFrame[i] = fr.x[index[i]];
 233     }
 234     frames++;
 235
 236     corrs.update(frnr, trajectoryFrame);
 237 }
 238
 239 void
 240 SpaceTimeCorr::finishAnalysis(int nframes)
 241 {
 242     std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > rout_new;
 243     unsigned int k = 500;
 244     if (tau > 0) {
 245         k = static_cast< unsigned int>(tau);
 246     }
 247     if (window == 0) {
 248         window = static_cast< int >(k * 2);
 249     }
 250
 251     corrs.printData();
 252     std::cout << "Finish Analysis - end\n\n";
 253 }
 254
 255 void
 256 SpaceTimeCorr::writeOutput()
 257 {
 258
 259 }
 260
 261 /*! \brief
 262  * The main function for the analysis template.
 263  */
 264 int
 265 main(int argc, char *argv[])
 266 {
 267     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SpaceTimeCorr>(argc, argv);
 268 }