src/gromacs/analysisdata/modules/displacement.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010-2018, The GROMACS development team.
   5  * Copyright (c) 2019,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements gmx::AnalysisDataDisplacementModule.
  39  *
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  * \ingroup module_analysisdata
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "displacement.h"
  46
  47 #include "gromacs/analysisdata/dataframe.h"
  48 #include "gromacs/analysisdata/datamodulemanager.h"
  49 #include "gromacs/analysisdata/modules/histogram.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/utility/exceptions.h"
  52 #include "gromacs/utility/gmxassert.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 namespace gmx
  56 {
  57
  58 /********************************************************************
  59  * AnalysisDataDisplacementModule::Impl
  60  */
  61
  62 /*! \internal \brief
  63  * Private implementation class for AnalysisDataDisplacementModule.
  64  *
  65  * \ingroup module_analysisdata
  66  */
  67 class AnalysisDataDisplacementModule::Impl
  68 {
  69 public:
  70     Impl();
  71     ~Impl();
  72
  73     //! Maximum number of particles for which the displacements are calculated.
  74     int nmax;
  75     //! Maximum time for which the displacements are needed.
  76     real tmax;
  77     //! Number of dimensions per data point.
  78     int ndim;
  79
  80     //! true if no frames have been read.
  81     bool bFirst;
  82     //! Stores the time of the first frame.
  83     real t0;
  84     //! Stores the time interval between frames.
  85     real dt;
  86     //! Stores the time of the current frame.
  87     real t;
  88     //! Stores the index in the store for the current positions.
  89     int ci;
  90
  91     //! Maximum number of positions to store for a particle.
  92     int max_store;
  93     //! The total number of positions ever stored (can be larger than \p max_store).
  94     int nstored;
  95     //! Old values.
  96     real* oldval;
  97     //! The most recently calculated displacements.
  98     std::vector<AnalysisDataValue> currValues_;
  99
 100     //! Histogram module for calculating MSD histograms, or NULL if not set.
 101     AnalysisDataBinAverageModule* histm;
 102 };
 103
 104 AnalysisDataDisplacementModule::Impl::Impl() :
 105     nmax(0),
 106     tmax(0.0),
 107     ndim(3),
 108     bFirst(true),
 109     t0(0.0),
 110     dt(0.0),
 111     t(0.0),
 112     ci(0),
 113     max_store(-1),
 114     nstored(0),
 115     oldval(nullptr),
 116     histm(nullptr)
 117 {
 118 }
 119
 120 AnalysisDataDisplacementModule::Impl::~Impl()
 121 {
 122     sfree(oldval);
 123 }
 124
 125 /********************************************************************
 126  * AnalysisDataDisplacementModule
 127  */
 128
 129 AnalysisDataDisplacementModule::AnalysisDataDisplacementModule() : _impl(new Impl())
 130 {
 131     setMultipoint(true);
 132 }
 133
 134
 135 AnalysisDataDisplacementModule::~AnalysisDataDisplacementModule() {}
 136
 137
 138 void AnalysisDataDisplacementModule::setMaxTime(real tmax)
 139 {
 140     _impl->tmax = tmax;
 141 }
 142
 143
 144 void AnalysisDataDisplacementModule::setMSDHistogram(const AnalysisDataBinAverageModulePointer& histm)
 145 {
 146     GMX_RELEASE_ASSERT(_impl->histm == nullptr, "Can only set MSD histogram once");
 147     _impl->histm = histm.get();
 148     addModule(histm);
 149 }
 150
 151
 152 AnalysisDataFrameRef AnalysisDataDisplacementModule::tryGetDataFrameInternal(int /*index*/) const
 153 {
 154     return AnalysisDataFrameRef();
 155 }
 156
 157
 158 bool AnalysisDataDisplacementModule::requestStorageInternal(int /*nframes*/)
 159 {
 160     return false;
 161 }
 162
 163
 164 int AnalysisDataDisplacementModule::flags() const
 165 {
 166     return efAllowMulticolumn;
 167 }
 168
 169
 170 void AnalysisDataDisplacementModule::dataStarted(AbstractAnalysisData* data)
 171 {
 172     if (data->columnCount() % _impl->ndim != 0)
 173     {
 174         GMX_THROW(APIError("Data has incorrect number of columns"));
 175     }
 176     _impl->nmax = data->columnCount();
 177     snew(_impl->oldval, _impl->nmax);
 178     _impl->ci = -_impl->nmax;
 179
 180     int ncol = _impl->nmax / _impl->ndim + 1;
 181     _impl->currValues_.reserve(ncol);
 182     setColumnCount(0, ncol);
 183 }
 184
 185
 186 void AnalysisDataDisplacementModule::frameStarted(const AnalysisDataFrameHeader& header)
 187 {
 188     // Initialize times.
 189     if (_impl->bFirst)
 190     {
 191         _impl->t0 = header.x();
 192     }
 193     else if (_impl->dt <= 0)
 194     {
 195         _impl->dt = header.x() - _impl->t0;
 196         if (_impl->dt < 0 || gmx_within_tol(_impl->dt, 0.0, GMX_REAL_EPS))
 197         {
 198             GMX_THROW(APIError("Identical or decreasing frame times"));
 199         }
 200     }
 201     else
 202     {
 203         if (!gmx_within_tol(header.x() - _impl->t, _impl->dt, GMX_REAL_EPS))
 204         {
 205             GMX_THROW(APIError("Frames not evenly spaced"));
 206         }
 207     }
 208     _impl->t = header.x();
 209
 210     // Allocate memory for all the positions once it is possible.
 211     if (_impl->max_store == -1 && !_impl->bFirst)
 212     {
 213         _impl->max_store = _impl->nmax * static_cast<int>(_impl->tmax / _impl->dt + 1);
 214         srenew(_impl->oldval, _impl->max_store);
 215     }
 216
 217     // Increment the index where current positions are stored.
 218     _impl->ci += _impl->nmax;
 219     if (_impl->ci >= _impl->max_store)
 220     {
 221         _impl->ci = 0;
 222     }
 223
 224     /*
 225         for (int i = 0; i < _impl->nmax; ++i)
 226         {
 227             _impl->p[_impl->ci + i].bPres = false;
 228         }
 229      */
 230     _impl->nstored++;
 231     _impl->bFirst = false;
 232 }
 233
 234
 235 void AnalysisDataDisplacementModule::pointsAdded(const AnalysisDataPointSetRef& points)
 236 {
 237     if (points.firstColumn() % _impl->ndim != 0 || points.columnCount() % _impl->ndim != 0)
 238     {
 239         GMX_THROW(APIError("Partial data points"));
 240     }
 241     for (int i = 0; i < points.columnCount(); ++i)
 242     {
 243         _impl->oldval[_impl->ci + points.firstColumn() + i] = points.y(i);
 244     }
 245 }
 246
 247
 248 void AnalysisDataDisplacementModule::frameFinished(const AnalysisDataFrameHeader& /*header*/)
 249 {
 250     if (_impl->nstored <= 1)
 251     {
 252         return;
 253     }
 254
 255     if (_impl->nstored == 2)
 256     {
 257         if (_impl->histm)
 258         {
 259             _impl->histm->init(
 260                     histogramFromBins(0, _impl->max_store / _impl->nmax, _impl->dt).integerBins());
 261         }
 262         moduleManager().notifyDataStart(this);
 263     }
 264     AnalysisDataFrameHeader header(_impl->nstored - 2, _impl->t, 0);
 265     moduleManager().notifyFrameStart(header);
 266
 267     for (int i = _impl->ci - _impl->nmax, step = 1; step < _impl->nstored && i != _impl->ci;
 268          i -= _impl->nmax, ++step)
 269     {
 270         if (i < 0)
 271         {
 272             i += _impl->max_store;
 273         }
 274         _impl->currValues_.clear();
 275         _impl->currValues_.emplace_back(step * _impl->dt);
 276         int k = 1;
 277         for (int j = 0; j < _impl->nmax; j += _impl->ndim, ++k)
 278         {
 279             real dist2 = 0.0;
 280
 281             for (int d = 0; d < _impl->ndim; ++d)
 282             {
 283                 real displ = _impl->oldval[_impl->ci + j + d] - _impl->oldval[i + j + d];
 284                 dist2 += displ * displ;
 285             }
 286             _impl->currValues_.emplace_back(dist2);
 287         }
 288         moduleManager().notifyPointsAdd(AnalysisDataPointSetRef(header, _impl->currValues_));
 289     }
 290
 291     moduleManager().notifyFrameFinish(header);
 292 }
 293
 294
 295 void AnalysisDataDisplacementModule::dataFinished()
 296 {
 297     if (_impl->nstored >= 2)
 298     {
 299         moduleManager().notifyDataFinish();
 300     }
 301 }
 302
 303 } // namespace gmx