src/gromacs/analysisdata/modules/displacement.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   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.
  13
  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.
  18  *
  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.
  25  *
  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.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements gmx::AnalysisDataDisplacementModule.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_analysisdata
  37  */
  38 #include "gromacs/analysisdata/modules/displacement.h"
  39
  40 #ifdef HAVE_CONFIG_H
  41 #include <config.h>
  42 #endif
  43
  44 // Legacy include.
  45 #include "smalloc.h"
  46
  47 #include "gromacs/basicmath.h"
  48 #include "gromacs/analysisdata/dataframe.h"
  49 #include "gromacs/analysisdata/modules/histogram.h"
  50 #include "gromacs/fatalerror/exceptions.h"
  51 #include "gromacs/fatalerror/gmxassert.h"
  52
  53 #include "displacement-impl.h"
  54
  55 namespace gmx
  56 {
  57
  58 /********************************************************************
  59  * AnalysisDataDisplacementModule::Impl
  60  */
  61
  62 AnalysisDataDisplacementModule::Impl::Impl()
  63     : nmax(0), tmax(0.0), ndim(3),
  64       bFirst(true), t0(0.0), dt(0.0), t(0.0),
  65       max_store(-1), nstored(0), oldval(NULL),
  66       histm(NULL)
  67 {
  68 }
  69
  70 AnalysisDataDisplacementModule::Impl::~Impl()
  71 {
  72     sfree(oldval);
  73 }
  74
  75 /********************************************************************
  76  * AnalysisDataDisplacementModule
  77  */
  78
  79 AnalysisDataDisplacementModule::AnalysisDataDisplacementModule()
  80     : _impl(new Impl())
  81 {
  82     setMultipoint(true);
  83 }
  84
  85
  86 AnalysisDataDisplacementModule::~AnalysisDataDisplacementModule()
  87 {
  88 }
  89
  90
  91 void
  92 AnalysisDataDisplacementModule::setMaxTime(real tmax)
  93 {
  94     _impl->tmax = tmax;
  95 }
  96
  97
  98 void
  99 AnalysisDataDisplacementModule::setMSDHistogram(AnalysisDataBinAverageModule *histm)
 100 {
 101     GMX_RELEASE_ASSERT(!_impl->histm, "Can only set MSD histogram once");
 102     _impl->histm = histm;
 103     addModule(histm);
 104 }
 105
 106
 107 AnalysisDataFrameRef
 108 AnalysisDataDisplacementModule::tryGetDataFrameInternal(int /*index*/) const
 109 {
 110     return AnalysisDataFrameRef();
 111 }
 112
 113
 114 bool
 115 AnalysisDataDisplacementModule::requestStorageInternal(int /*nframes*/)
 116 {
 117     return false;
 118 }
 119
 120
 121 int
 122 AnalysisDataDisplacementModule::flags() const
 123 {
 124     return efAllowMulticolumn;
 125 }
 126
 127
 128 void
 129 AnalysisDataDisplacementModule::dataStarted(AbstractAnalysisData *data)
 130 {
 131     if (data->columnCount() % _impl->ndim != 0)
 132     {
 133         GMX_THROW(APIError("Data has incorrect number of columns"));
 134     }
 135     _impl->nmax = data->columnCount();
 136     snew(_impl->oldval, _impl->nmax);
 137     _impl->ci = -_impl->nmax;
 138
 139     int ncol = _impl->nmax / _impl->ndim + 1;
 140     _impl->currValues_.reserve(ncol);
 141     setColumnCount(ncol);
 142 }
 143
 144
 145 void
 146 AnalysisDataDisplacementModule::frameStarted(const AnalysisDataFrameHeader &header)
 147 {
 148     // Initialize times.
 149     if (_impl->bFirst)
 150     {
 151         _impl->t0 = header.x();
 152     }
 153     else if (_impl->dt <= 0)
 154     {
 155         _impl->dt = header.x() - _impl->t0;
 156         if (_impl->dt < 0 || gmx_within_tol(_impl->dt, 0.0, GMX_REAL_EPS))
 157         {
 158             GMX_THROW(APIError("Identical or decreasing frame times"));
 159         }
 160     }
 161     else
 162     {
 163         if (!gmx_within_tol(header.x() - _impl->t, _impl->dt, GMX_REAL_EPS))
 164         {
 165             GMX_THROW(APIError("Frames not evenly spaced"));
 166         }
 167     }
 168     _impl->t = header.x();
 169
 170     // Allocate memory for all the positions once it is possible.
 171     if (_impl->max_store == -1 && !_impl->bFirst)
 172     {
 173         _impl->max_store = _impl->nmax * (int)(_impl->tmax/_impl->dt + 1);
 174         srenew(_impl->oldval, _impl->max_store);
 175     }
 176
 177     // Increment the index where current positions are stored.
 178     _impl->ci += _impl->nmax;
 179     if (_impl->ci >= _impl->max_store)
 180     {
 181         _impl->ci = 0;
 182     }
 183
 184 /*
 185     for (int i = 0; i < _impl->nmax; ++i)
 186     {
 187         _impl->p[_impl->ci + i].bPres = false;
 188     }
 189 */
 190     _impl->nstored++;
 191     _impl->bFirst = false;
 192 }
 193
 194
 195 void
 196 AnalysisDataDisplacementModule::pointsAdded(const AnalysisDataPointSetRef &points)
 197 {
 198     if (points.firstColumn() % _impl->ndim != 0
 199         || points.columnCount() % _impl->ndim != 0)
 200     {
 201         GMX_THROW(APIError("Partial data points"));
 202     }
 203     for (int i = 0; i < points.columnCount(); ++i)
 204     {
 205         _impl->oldval[_impl->ci + points.firstColumn() + i] = points.y(i);
 206     }
 207 }
 208
 209
 210 void
 211 AnalysisDataDisplacementModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 212 {
 213     if (_impl->nstored <= 1)
 214     {
 215         return;
 216     }
 217
 218     int step, i;
 219     int rc;
 220
 221     if (_impl->nstored == 2)
 222     {
 223         if (_impl->histm)
 224         {
 225             _impl->histm->init(histogramFromBins(0, _impl->max_store / _impl->nmax,
 226                                                  _impl->dt).integerBins());
 227         }
 228         notifyDataStart();
 229     }
 230     AnalysisDataFrameHeader header(_impl->nstored - 2, _impl->t, 0);
 231     notifyFrameStart(header);
 232
 233     for (i = _impl->ci - _impl->nmax, step = 1;
 234          step < _impl->nstored && i != _impl->ci;
 235          i -= _impl->nmax, ++step)
 236     {
 237         if (i < 0)
 238         {
 239             i += _impl->max_store;
 240         }
 241         _impl->currValues_.clear();
 242         _impl->currValues_.push_back(AnalysisDataValue(step * _impl->dt));
 243         int k = 1;
 244         for (int j = 0; j < _impl->nmax; j += _impl->ndim, ++k)
 245         {
 246             real dist2 = 0.0;
 247
 248             for (int d = 0; d < _impl->ndim; ++d)
 249             {
 250                 dist2 += sqr(_impl->oldval[_impl->ci + j + d]
 251                              - _impl->oldval[i + j + d]);
 252             }
 253             _impl->currValues_.push_back(AnalysisDataValue(dist2));
 254         }
 255         notifyPointsAdd(AnalysisDataPointSetRef(header, _impl->currValues_));
 256     }
 257
 258     notifyFrameFinish(header);
 259 }
 260
 261
 262 void
 263 AnalysisDataDisplacementModule::dataFinished()
 264 {
 265     if (_impl->nstored >= 2)
 266     {
 267         notifyDataFinish();
 268     }
 269 }
 270
 271 } // namespace gmx