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(
 100         AnalysisDataBinAverageModulePointer histm)
 101 {
 102     GMX_RELEASE_ASSERT(_impl->histm == NULL, "Can only set MSD histogram once");
 103     _impl->histm = histm.get();
 104     addModule(histm);
 105 }
 106
 107
 108 AnalysisDataFrameRef
 109 AnalysisDataDisplacementModule::tryGetDataFrameInternal(int /*index*/) const
 110 {
 111     return AnalysisDataFrameRef();
 112 }
 113
 114
 115 bool
 116 AnalysisDataDisplacementModule::requestStorageInternal(int /*nframes*/)
 117 {
 118     return false;
 119 }
 120
 121
 122 int
 123 AnalysisDataDisplacementModule::flags() const
 124 {
 125     return efAllowMulticolumn;
 126 }
 127
 128
 129 void
 130 AnalysisDataDisplacementModule::dataStarted(AbstractAnalysisData *data)
 131 {
 132     if (data->columnCount() % _impl->ndim != 0)
 133     {
 134         GMX_THROW(APIError("Data has incorrect number of columns"));
 135     }
 136     _impl->nmax = data->columnCount();
 137     snew(_impl->oldval, _impl->nmax);
 138     _impl->ci = -_impl->nmax;
 139
 140     int ncol = _impl->nmax / _impl->ndim + 1;
 141     _impl->currValues_.reserve(ncol);
 142     setColumnCount(ncol);
 143 }
 144
 145
 146 void
 147 AnalysisDataDisplacementModule::frameStarted(const AnalysisDataFrameHeader &header)
 148 {
 149     // Initialize times.
 150     if (_impl->bFirst)
 151     {
 152         _impl->t0 = header.x();
 153     }
 154     else if (_impl->dt <= 0)
 155     {
 156         _impl->dt = header.x() - _impl->t0;
 157         if (_impl->dt < 0 || gmx_within_tol(_impl->dt, 0.0, GMX_REAL_EPS))
 158         {
 159             GMX_THROW(APIError("Identical or decreasing frame times"));
 160         }
 161     }
 162     else
 163     {
 164         if (!gmx_within_tol(header.x() - _impl->t, _impl->dt, GMX_REAL_EPS))
 165         {
 166             GMX_THROW(APIError("Frames not evenly spaced"));
 167         }
 168     }
 169     _impl->t = header.x();
 170
 171     // Allocate memory for all the positions once it is possible.
 172     if (_impl->max_store == -1 && !_impl->bFirst)
 173     {
 174         _impl->max_store = _impl->nmax * (int)(_impl->tmax/_impl->dt + 1);
 175         srenew(_impl->oldval, _impl->max_store);
 176     }
 177
 178     // Increment the index where current positions are stored.
 179     _impl->ci += _impl->nmax;
 180     if (_impl->ci >= _impl->max_store)
 181     {
 182         _impl->ci = 0;
 183     }
 184
 185 /*
 186     for (int i = 0; i < _impl->nmax; ++i)
 187     {
 188         _impl->p[_impl->ci + i].bPres = false;
 189     }
 190 */
 191     _impl->nstored++;
 192     _impl->bFirst = false;
 193 }
 194
 195
 196 void
 197 AnalysisDataDisplacementModule::pointsAdded(const AnalysisDataPointSetRef &points)
 198 {
 199     if (points.firstColumn() % _impl->ndim != 0
 200         || points.columnCount() % _impl->ndim != 0)
 201     {
 202         GMX_THROW(APIError("Partial data points"));
 203     }
 204     for (int i = 0; i < points.columnCount(); ++i)
 205     {
 206         _impl->oldval[_impl->ci + points.firstColumn() + i] = points.y(i);
 207     }
 208 }
 209
 210
 211 void
 212 AnalysisDataDisplacementModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 213 {
 214     if (_impl->nstored <= 1)
 215     {
 216         return;
 217     }
 218
 219     int step, i;
 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