3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2009, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
33 * Implements gmx::AnalysisDataDisplacementModule.
35 * \author Teemu Murtola <teemu.murtola@cbr.su.se>
36 * \ingroup module_analysisdata
38 #include "gromacs/analysisdata/modules/displacement.h"
40 #include "gromacs/legacyheaders/smalloc.h"
41 #include "gromacs/legacyheaders/maths.h"
43 #include "gromacs/analysisdata/dataframe.h"
44 #include "gromacs/analysisdata/modules/histogram.h"
45 #include "gromacs/utility/exceptions.h"
46 #include "gromacs/utility/gmxassert.h"
51 /********************************************************************
52 * AnalysisDataDisplacementModule::Impl
56 * Private implementation class for AnalysisDataDisplacementModule.
58 * \ingroup module_analysisdata
60 class AnalysisDataDisplacementModule::Impl
66 //! Maximum number of particles for which the displacements are calculated.
68 //! Maximum time for which the displacements are needed.
70 //! Number of dimensions per data point.
73 //! true if no frames have been read.
75 //! Stores the time of the first frame.
77 //! Stores the time interval between frames.
79 //! Stores the time of the current frame.
81 //! Stores the index in the store for the current positions.
84 //! Maximum number of positions to store for a particle.
86 //! The total number of positions ever stored (can be larger than \p max_store).
90 //! The most recently calculated displacements.
91 std::vector<AnalysisDataValue> currValues_;
93 //! Histogram module for calculating MSD histograms, or NULL if not set.
94 AnalysisDataBinAverageModule *histm;
97 AnalysisDataDisplacementModule::Impl::Impl()
98 : nmax(0), tmax(0.0), ndim(3),
99 bFirst(true), t0(0.0), dt(0.0), t(0.0),
100 max_store(-1), nstored(0), oldval(NULL),
105 AnalysisDataDisplacementModule::Impl::~Impl()
110 /********************************************************************
111 * AnalysisDataDisplacementModule
114 AnalysisDataDisplacementModule::AnalysisDataDisplacementModule()
121 AnalysisDataDisplacementModule::~AnalysisDataDisplacementModule()
127 AnalysisDataDisplacementModule::setMaxTime(real tmax)
134 AnalysisDataDisplacementModule::setMSDHistogram(
135 AnalysisDataBinAverageModulePointer histm)
137 GMX_RELEASE_ASSERT(_impl->histm == NULL, "Can only set MSD histogram once");
138 _impl->histm = histm.get();
144 AnalysisDataDisplacementModule::tryGetDataFrameInternal(int /*index*/) const
146 return AnalysisDataFrameRef();
151 AnalysisDataDisplacementModule::requestStorageInternal(int /*nframes*/)
158 AnalysisDataDisplacementModule::flags() const
160 return efAllowMulticolumn;
165 AnalysisDataDisplacementModule::dataStarted(AbstractAnalysisData *data)
167 if (data->columnCount() % _impl->ndim != 0)
169 GMX_THROW(APIError("Data has incorrect number of columns"));
171 _impl->nmax = data->columnCount();
172 snew(_impl->oldval, _impl->nmax);
173 _impl->ci = -_impl->nmax;
175 int ncol = _impl->nmax / _impl->ndim + 1;
176 _impl->currValues_.reserve(ncol);
177 setColumnCount(ncol);
182 AnalysisDataDisplacementModule::frameStarted(const AnalysisDataFrameHeader &header)
187 _impl->t0 = header.x();
189 else if (_impl->dt <= 0)
191 _impl->dt = header.x() - _impl->t0;
192 if (_impl->dt < 0 || gmx_within_tol(_impl->dt, 0.0, GMX_REAL_EPS))
194 GMX_THROW(APIError("Identical or decreasing frame times"));
199 if (!gmx_within_tol(header.x() - _impl->t, _impl->dt, GMX_REAL_EPS))
201 GMX_THROW(APIError("Frames not evenly spaced"));
204 _impl->t = header.x();
206 // Allocate memory for all the positions once it is possible.
207 if (_impl->max_store == -1 && !_impl->bFirst)
209 _impl->max_store = _impl->nmax * (int)(_impl->tmax/_impl->dt + 1);
210 srenew(_impl->oldval, _impl->max_store);
213 // Increment the index where current positions are stored.
214 _impl->ci += _impl->nmax;
215 if (_impl->ci >= _impl->max_store)
221 for (int i = 0; i < _impl->nmax; ++i)
223 _impl->p[_impl->ci + i].bPres = false;
227 _impl->bFirst = false;
232 AnalysisDataDisplacementModule::pointsAdded(const AnalysisDataPointSetRef &points)
234 if (points.firstColumn() % _impl->ndim != 0
235 || points.columnCount() % _impl->ndim != 0)
237 GMX_THROW(APIError("Partial data points"));
239 for (int i = 0; i < points.columnCount(); ++i)
241 _impl->oldval[_impl->ci + points.firstColumn() + i] = points.y(i);
247 AnalysisDataDisplacementModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
249 if (_impl->nstored <= 1)
256 if (_impl->nstored == 2)
260 _impl->histm->init(histogramFromBins(0, _impl->max_store / _impl->nmax,
261 _impl->dt).integerBins());
265 AnalysisDataFrameHeader header(_impl->nstored - 2, _impl->t, 0);
266 notifyFrameStart(header);
268 for (i = _impl->ci - _impl->nmax, step = 1;
269 step < _impl->nstored && i != _impl->ci;
270 i -= _impl->nmax, ++step)
274 i += _impl->max_store;
276 _impl->currValues_.clear();
277 _impl->currValues_.push_back(AnalysisDataValue(step * _impl->dt));
279 for (int j = 0; j < _impl->nmax; j += _impl->ndim, ++k)
283 for (int d = 0; d < _impl->ndim; ++d)
285 real displ = _impl->oldval[_impl->ci + j + d]
286 - _impl->oldval[i + j + d];
287 dist2 += displ * displ;
289 _impl->currValues_.push_back(AnalysisDataValue(dist2));
291 notifyPointsAdd(AnalysisDataPointSetRef(header, _impl->currValues_));
294 notifyFrameFinish(header);
299 AnalysisDataDisplacementModule::dataFinished()
301 if (_impl->nstored >= 2)