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 <chrono>
  45 #include <omp.h>
  46 #include <thread>
  47 #include <string>
  48
  49 #include <gromacs/trajectoryanalysis.h>
  50 #include <gromacs/pbcutil/pbc.h>
  51 #include <gromacs/pbcutil/rmpbc.h>
  52 #include <gromacs/utility/smalloc.h>
  53 #include <gromacs/math/vectypes.h>
  54 #include <gromacs/math/vec.h>
  55 #include <gromacs/math/do_fit.h>
  56
  57
  58
  59
  60
  61
  62
  63
  64 #include <string>
  65
  66 #include "gromacs/analysisdata/analysisdata.h"
  67 #include "gromacs/analysisdata/modules/average.h"
  68 #include "gromacs/analysisdata/modules/histogram.h"
  69 #include "gromacs/analysisdata/modules/plot.h"
  70 #include "gromacs/math/vec.h"
  71 #include "gromacs/options/basicoptions.h"
  72 #include "gromacs/options/filenameoption.h"
  73 #include "gromacs/options/ioptionscontainer.h"
  74 #include "gromacs/pbcutil/pbc.h"
  75 #include "gromacs/selection/selection.h"
  76 #include "gromacs/selection/selectionoption.h"
  77 #include "gromacs/trajectory/trajectoryframe.h"
  78 #include "gromacs/trajectoryanalysis/analysissettings.h"
  79 #include "gromacs/utility/arrayref.h"
  80 #include "gromacs/utility/exceptions.h"
  81 #include "gromacs/utility/stringutil.h"
  82
  83
  84
  85
  86
  87 using namespace gmx;
  88
  89 using gmx::RVec;
  90
  91 /*! \brief
  92  * Class used to compute free volume in a simulations box.
  93  *
  94  * Inherits TrajectoryAnalysisModule and all functions from there.
  95  * Does not implement any new functionality.
  96  *
  97  * \ingroup module_trajectoryanalysis
  98  */
  99 class Domains : public TrajectoryAnalysisModule
 100 {
 101     public:
 102
 103         Domains();
 104         virtual ~Domains();
 105
 106         //! Set the options and setting
 107         virtual void initOptions(IOptionsContainer          *options,
 108                                  TrajectoryAnalysisSettings *settings);
 109
 110         //! First routine called by the analysis framework
 111         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 112         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 113                                   const TopologyInformation        &top);
 114
 115         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 116                                          const t_trxframe                   &fr);
 117
 118         //! Call for each frame of the trajectory
 119         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 120         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 121                                   TrajectoryAnalysisModuleData *pdata);
 122
 123         //! Last routine called by the analysis framework
 124         // virtual void finishAnalysis(t_pbc *pbc);
 125         virtual void finishAnalysis(int nframes);
 126
 127         //! Routine to write output, that is additional over the built-in
 128         virtual void writeOutput();
 129
 130     private:
 131
 132         std::string                                                 fnNdx_;
 133         SelectionList                                               sel_;
 134         Selection                                                   selec;
 135
 136         std::vector< std::vector< int > >                           domains;
 137         std::vector< std::vector< RVec > >                          trajectory;
 138         std::vector< std::vector< RVec > >                          frankenstein_trajectory;
 139
 140
 141         std::vector< int >                                          index;
 142         std::vector< int >                                          numbers;
 143         int                                                         frames              = 0;
 144         int                                                         basic_frame         = 0;
 145         real                                                        **w_rls;
 146
 147         int                                                         domains_ePBC;
 148         // Copy and assign disallowed by base.
 149 };
 150
 151 Domains::Domains(): TrajectoryAnalysisModule()
 152 {
 153 }
 154
 155 Domains::~Domains()
 156 {
 157 }
 158
 159 void
 160 Domains::initOptions(IOptionsContainer          *options,
 161                      TrajectoryAnalysisSettings *settings)
 162 {
 163     static const char *const desc[] = {
 164         "[THISMODULE] to be done"
 165     };
 166     // Add the descriptive text (program help text) to the options
 167     settings->setHelpText(desc);
 168     // Add option for selecting a subset of atoms
 169     options->addOption(SelectionOption("select_type")
 170                            .store(&selec).required()
 171                            .description("Atoms that are considered as part of the excluded volume"));
 172     // Add option for output file name
 173     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 174                             .store(&fnNdx_).defaultBasename("domains")
 175                             .description("Index file from the domains"));
 176     options->addOption(SelectionOption("select_domains").storeVector(&sel_)
 177                            .required().dynamicMask().multiValue()
 178                            .description("Domains to form rigid skeleton"));
 179     // Control input settings
 180     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 181     settings->setPBC(true);
 182 }
 183
 184 void
 185 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 186                       const TopologyInformation        &top)
 187 {
 188     domains_ePBC = top.ePBC();
 189 }
 190
 191 void
 192 Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 193                              const t_trxframe                       &fr)
 194 {
 195     ConstArrayRef< int > atomind  = selec.atomIndices();
 196     index.resize(0);
 197     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 198         index.push_back(*ai);
 199     }
 200
 201     //checkSelections(sel_);
 202     domains.resize(sel_.size());
 203     for (int i = 0; i < sel_.size(); i++) {
 204         for (int j = 0; j < sel_[i].size(); j++) { //distance / modules / trajectoryanalyziz
 205             domains[i].push_back(sel_[i][j]);
 206         }
 207     }
 208
 209     trajectory.resize(1);
 210     trajectory[0].resize(selec.atomCount());
 211
 212     for (int i = 0; i < selec.atomCount(); i++) {
 213         trajectory[0][i] = fr.x[index[i]];
 214     }
 215
 216     frames++;
 217 }
 218
 219 void
 220 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 221                       TrajectoryAnalysisModuleData *pdata)
 222 {
 223     std::vector< std::vector< int > >   domains_local;
 224     std::vector< bool >                 flags;
 225
 226     flags.resize(index.size(), true);
 227     domains_local = domains; // технически оно за О(1) должно делаться, но не точно
 228
 229     trajectory.resize(trajectory.size() + 1);
 230     trajectory.back().resize(index.size());
 231
 232     for (int i = 0; i < index.size(); i++) {
 233         trajectory.back()[i] = fr.x[index[i]];
 234     }
 235
 236     for (int i = 0; i < domains.size(); i++) {
 237         for (int j = 0; j < domains[i].size(); j++) {
 238             flags[domains[i][j]] = false;
 239         }
 240     }
 241
 242     for (int i = 0; i < index.size(); i++) {
 243         if (flags[i]) {
 244             int a, b, dist = 9001;
 245             rvec temp;
 246             for (int j = 0; j < domains.size(); j++) {
 247                 for (int k = 0; k < domains[j].size(); k++) {
 248                     rvec_sub(trajectory.back()[i], trajectory.back()[domains[j][k]], temp);
 249                     if (norm(temp) <= dist) {
 250                         dist = norm(temp);
 251                         a = j;
 252                         b = k;
 253                     }
 254                 }
 255             }
 256             domains_local[a].push_back(i);
 257             flags[i] = false;
 258         }
 259     }
 260
 261     snew(w_rls, domains_local.size());
 262     for (int i = 0; i < domains_local.size(); i++) {
 263         snew(w_rls[i], index.size());
 264         for (int j = 0; j < index.size(); j++) {
 265             w_rls[i][j] = 0;
 266         }
 267         for (int j = 0; j < domains_local[i].size(); j++) {
 268             w_rls[i][domains_local[i][j]] = 1;
 269         }
 270     }
 271
 272     frankenstein_trajectory.resize(trajectory.size());
 273     frankenstein_trajectory.back().resize(index.size());
 274
 275     for (int i = 0; i < domains_local.size(); i++) {
 276         rvec *basic, *traj;
 277         snew(basic, index.size());
 278         for (int k = 0; k < index.size(); k++) {
 279             copy_rvec(trajectory[basic_frame][k].as_vec(), basic[k]);
 280         }
 281         snew(traj, index.size());
 282         for (int k = 0; k < index.size(); k++) {
 283             copy_rvec(trajectory.back()[k].as_vec(), traj[k]);
 284         }
 285         reset_x(index.size(), NULL, index.size(), NULL, basic, w_rls[i]);
 286         reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[i]);
 287         do_fit(index.size(), w_rls[i], basic, traj);
 288
 289         for (int j = 0; j < domains_local[i].size(); j++) {
 290             frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
 291         }
 292
 293         sfree(basic);
 294         sfree(traj);
 295     }
 296
 297     frames++;
 298 }
 299
 300 //domains -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
 301
 302 void
 303 Domains::finishAnalysis(int nframes)
 304 {
 305
 306 }
 307
 308 void
 309 Domains::writeOutput()
 310 {
 311
 312 }
 313
 314
 315 /*! \brief
 316  * The main function for the analysis template.
 317  */
 318 int
 319 main(int argc, char *argv[])
 320 {
 321     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 322 }