src/gromacs/trajectoryanalysis/modules/distance.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010,2011,2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * 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::Distance.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42 #include "distance.h"
  43
  44 #include "gromacs/legacyheaders/pbc.h"
  45 #include "gromacs/legacyheaders/vec.h"
  46
  47 #include "gromacs/analysisdata/analysisdata.h"
  48 #include "gromacs/analysisdata/modules/plot.h"
  49 #include "gromacs/options/basicoptions.h"
  50 #include "gromacs/options/filenameoption.h"
  51 #include "gromacs/options/options.h"
  52 #include "gromacs/selection/selection.h"
  53 #include "gromacs/selection/selectionoption.h"
  54 #include "gromacs/trajectoryanalysis/analysissettings.h"
  55 #include "gromacs/utility/exceptions.h"
  56 #include "gromacs/utility/stringutil.h"
  57
  58 namespace gmx
  59 {
  60
  61 namespace analysismodules
  62 {
  63
  64 const char Distance::name[]             = "distance";
  65 const char Distance::shortDescription[] =
  66     "Calculate distances";
  67
  68 Distance::Distance()
  69     : TrajectoryAnalysisModule(name, shortDescription),
  70       avem_(new AnalysisDataAverageModule())
  71 {
  72     data_.setColumnCount(4);
  73     registerAnalysisDataset(&data_, "distance");
  74 }
  75
  76
  77 Distance::~Distance()
  78 {
  79 }
  80
  81
  82 void
  83 Distance::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
  84 {
  85     static const char *const desc[] = {
  86         "g_dist can calculate the distance between two positions as",
  87         "a function of time. The total distance and its",
  88         "x, y and z components are plotted."
  89     };
  90
  91     options->setDescription(concatenateStrings(desc));
  92
  93     options->addOption(FileNameOption("o").filetype(eftPlot).outputFile()
  94                            .store(&fnDist_).defaultBasename("dist")
  95                            .description("Computed distances"));
  96     options->addOption(SelectionOption("select").required().valueCount(2)
  97                            .store(sel_));
  98 }
  99
 100
 101 void
 102 Distance::initAnalysis(const TrajectoryAnalysisSettings &settings,
 103                        const TopologyInformation         & /*top*/)
 104 {
 105     if (sel_[0].posCount() != 1)
 106     {
 107         GMX_THROW(InvalidInputError("The first selection does not define a single position"));
 108     }
 109     if (sel_[1].posCount() != 1)
 110     {
 111         GMX_THROW(InvalidInputError("The second selection does not define a single position"));
 112     }
 113
 114     data_.addModule(avem_);
 115     AnalysisDataPlotModulePointer plotm_(new AnalysisDataPlotModule());
 116     plotm_->setSettings(settings.plotSettings());
 117     plotm_->setFileName(fnDist_);
 118     plotm_->setTitle("Distance");
 119     plotm_->setXAxisIsTime();
 120     plotm_->setYLabel("Distance (nm)");
 121     data_.addModule(plotm_);
 122 }
 123
 124
 125 void
 126 Distance::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 127                        TrajectoryAnalysisModuleData *pdata)
 128 {
 129     AnalysisDataHandle       dh   = pdata->dataHandle(data_);
 130     const Selection         &sel1 = pdata->parallelSelection(sel_[0]);
 131     const Selection         &sel2 = pdata->parallelSelection(sel_[1]);
 132     rvec                     dx;
 133     real                     r;
 134     const SelectionPosition &p1 = sel1.position(0);
 135     const SelectionPosition &p2 = sel2.position(0);
 136
 137     if (pbc != NULL)
 138     {
 139         pbc_dx(pbc, p1.x(), p2.x(), dx);
 140     }
 141     else
 142     {
 143         rvec_sub(p1.x(), p2.x(), dx);
 144     }
 145     r = norm(dx);
 146     dh.startFrame(frnr, fr.time);
 147     dh.setPoint(0, r);
 148     dh.setPoints(1, 3, dx);
 149     dh.finishFrame();
 150 }
 151
 152
 153 void
 154 Distance::finishAnalysis(int /*nframes*/)
 155 {
 156 }
 157
 158
 159 void
 160 Distance::writeOutput()
 161 {
 162     fprintf(stderr, "Average distance: %f\n", avem_->average(0));
 163     fprintf(stderr, "Std. deviation:   %f\n", avem_->stddev(0));
 164 }
 165
 166 } // namespace analysismodules
 167
 168 } // namespace gmx