src/gromacs/applied_forces/awh/correlationtensor.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source 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
  36 /*! \internal \file
  37  *
  38  * \brief
  39  * Declares the CorrelationTensor class for correlation tensor statistics.
  40  *
  41  * \author Viveca Lindahl
  42  * \author Berk Hess <hess@kth.se>
  43  * \ingroup module_awh
  44  */
  45
  46 #ifndef GMX_AWH_CORRELATIONTENSOR_H
  47 #define GMX_AWH_CORRELATIONTENSOR_H
  48
  49 #include <cstddef>
  50
  51 #include <vector>
  52
  53 #include "gromacs/utility/arrayref.h"
  54 #include "gromacs/utility/basedefinitions.h"
  55 #include "gromacs/utility/gmxassert.h"
  56
  57 namespace gmx
  58 {
  59
  60 struct CorrelationBlockDataHistory;
  61
  62 /*! \internal \brief Correlation block averaging weight-only data.
  63  */
  64 class CorrelationBlockData
  65 {
  66 public:
  67     /*! \internal \brief Correlation block averaging data.
  68      */
  69     struct CoordData
  70     {
  71         /*! \brief Constructor.
  72          */
  73         CoordData() : blockSumWeightX(0), sumOverBlocksBlockWeightBlockWeightX(0) {}
  74
  75         double blockSumWeightX; /**< Weighted sum of x for current block. */
  76         double sumOverBlocksBlockWeightBlockWeightX; /**< Sum over all blocks in the simulation of block weight times sum_wx. */
  77     };
  78
  79     /*! \brief Constructor.
  80      *
  81      * \param[in] numDim           The dimensionality.
  82      * \param[in] blockLengthInit  The initial block length.
  83      */
  84     CorrelationBlockData(int numDim, double blockLengthInit) :
  85         blockSumWeight_(0),
  86         blockSumSquareWeight_(0),
  87         sumOverBlocksSquareBlockWeight_(0),
  88         sumOverBlocksBlockSquareWeight_(0),
  89         blockLength_(blockLengthInit),
  90         previousBlockIndex_(-1),
  91         coordData_(numDim),
  92         correlationIntegral_(numDim * (numDim + 1) / 2)
  93     {
  94     }
  95
  96     /*! \brief Restore the state from history.
  97      *
  98      * \param[in] blockHistory         The block data history containing the weight sums.
  99      * \param[in] coordData            The coordinate data.
 100      * \param[in] correlationIntegral  The correlation integral for all tensor elements.
 101      */
 102     void restoreFromHistory(const CorrelationBlockDataHistory& blockHistory,
 103                             const std::vector<CoordData>&      coordData,
 104                             const std::vector<double>&         correlationIntegral);
 105
 106     /*! \brief Adds a weighted data vector to one point in the correlation grid.
 107      *
 108      * \param[in] weight  The weight of the data.
 109      * \param[in] data    One data point for each grid dimension.
 110      */
 111     void addData(double weight, gmx::ArrayRef<const double> data)
 112     {
 113         GMX_ASSERT(data.size() == coordData_.size(),
 114                    "Size of data should match the size of coordData");
 115
 116         blockSumWeight_ += weight;
 117         blockSumSquareWeight_ += weight * weight;
 118
 119         for (size_t d = 0; d < coordData_.size(); d++)
 120         {
 121             coordData_[d].blockSumWeightX += weight * data[d];
 122         }
 123     }
 124
 125     /*! \brief Adds a filled data block to correlation time integral.
 126      */
 127     void addBlockToCorrelationIntegral();
 128
 129     /*! \brief Returns the sum weights for current block. */
 130     double blockSumWeight() const { return blockSumWeight_; }
 131
 132     /*! \brief Returns the sum weights^2 for current block. */
 133     double blockSumSquareWeight() const { return blockSumSquareWeight_; }
 134
 135     /*! \brief Returns the sum over blocks of block weight^2. */
 136     double sumOverBlocksSquareBlockWeight() const { return sumOverBlocksSquareBlockWeight_; }
 137
 138     /*! \brief Returns the sum over blocks of weight^2. */
 139     double sumOverBlocksBlockSquareWeight() const { return sumOverBlocksBlockSquareWeight_; }
 140
 141     /*! \brief Returns the length of each block used for block averaging. */
 142     double blockLength() const { return blockLength_; }
 143
 144     /*! \brief Double the length of each block used for block averaging. */
 145     void doubleBlockLength() { blockLength_ *= 2; }
 146
 147     /*! \brief Return the last block index data was added to (needed only for block length in terms of time). */
 148     int previousBlockIndex() const { return previousBlockIndex_; }
 149
 150     /*! \brief Set the last block index data was added to.
 151      *
 152      * \param[in] blockIndex  The block index.
 153      */
 154     void setPreviousBlockIndex(int blockIndex) { previousBlockIndex_ = blockIndex; }
 155
 156     /*! \brief Return sums for each coordinate dimension. */
 157     const std::vector<CoordData>& coordData() const { return coordData_; }
 158
 159     /*! \brief Return the correlation integral tensor. */
 160     const std::vector<double>& correlationIntegral() const { return correlationIntegral_; }
 161
 162 private:
 163     /* Weight sum data, indentical for all dimensions */
 164     double blockSumWeight_;                 /**< Sum weights for current block. */
 165     double blockSumSquareWeight_;           /**< Sum weights^2 for current block. */
 166     double sumOverBlocksSquareBlockWeight_; /**< Sum over all blocks in the simulation of block weight^2. */
 167     double sumOverBlocksBlockSquareWeight_; /**< Sum over all blocks in the simulation of weight^2. */
 168     double blockLength_; /**< The length of each block used for block averaging */
 169     int    previousBlockIndex_; /**< The last block index data was added to (needed only for block length in terms of time). */
 170
 171     /* Sums for each coordinate dimension. */
 172     std::vector<CoordData> coordData_; /**< Array with sums for each coordinate dimension. */
 173
 174     /* Correlation tensor. */
 175     std::vector<double> correlationIntegral_; /**< Array with the correlation elements corr(x, y) in the tensor, where x, y are vector components. */
 176 };
 177
 178 /*! \internal
 179  * \brief Correlation data for computing the correlation tensor of one grid point.
 180  *
 181  * The time integrated autocorrelation of the desired quantity is computed using
 182  * block averages, which is a computationally efficient and low memory method.
 183  * Most of the work here goes into computing the block averages for weights
 184  * and the coordinate quantity. This is done for a number of blocks in
 185  * the range of \p c_numCorrelationBlocks/2 + 1 to \p c_numCorrelationBlocks,
 186  * depending on the current simulation length. As the simulation time
 187  * progresses, the blocks get longer. This is implemented in an efficient
 188  * way by keeping track of log2(\p c_numCorrelationBlocks) \p BlockData
 189  * data blocks with block length increasing progressively by a factor of 2.
 190  * Once \p c_numCorrelationBlocks are reached, all block lengths are doubled.
 191  */
 192 class CorrelationTensor
 193 {
 194 public:
 195     /*! \brief 64 blocks is a good trade-off between signal and noise */
 196     static constexpr int c_numCorrelationBlocks = 64;
 197
 198     /*! \brief Constructor.
 199      *
 200      * \param[in] numDim          The dimensionality.
 201      * \param[in] numBlockData     The number of data block data structs.
 202      * \param[in] blockLengthInit  The initial block length.
 203      */
 204     CorrelationTensor(int numDim, int numBlockData, double blockLengthInit);
 205
 206     /*! \brief Get a const reference to the list of block data.
 207      */
 208     const std::vector<CorrelationBlockData>& blockDataList() const { return blockDataList_; }
 209
 210     /*! \brief
 211      * Get the total weight of the data in the correlation matrix.
 212      *
 213      * \returns the weight of the added data.
 214      */
 215     double getWeight() const
 216     {
 217         /* The last blockdata has only 1 block containing all data */
 218         return blockDataList().back().blockSumWeight();
 219     }
 220
 221     /*! \brief Restore a correlation element from history.
 222      *
 223      * \param[in]     blockDataBuffer  The linear correlation grid data history buffer.
 224      * \param[in,out] bufferIndex      The index in \p blockDataBuffer to start reading, is increased with the number of blocks read.
 225      */
 226     void restoreFromHistory(const std::vector<CorrelationBlockDataHistory>& blockDataBuffer,
 227                             size_t*                                         bufferIndex);
 228
 229 private:
 230     /*! \brief Updates the block length by doubling.
 231      *
 232      * The length of all blocks is doubled. This is achieved by removing
 233      * the shortest block, moving all other blocks and duplicating
 234      * the data of longest block to a nw block of double length (but
 235      * currenly only half filled with data).
 236      */
 237     void doubleBlockLengths();
 238
 239     /*! \brief Updates the block length such that data fits.
 240      *
 241      * \param[in] samplingLength  Sampling length of all data, in time or weight.
 242      */
 243     void updateBlockLengths(double samplingLength);
 244
 245 public:
 246     /*! \brief Adds a weighted data vector to one point in the correlation grid.
 247      *
 248      * \note To avoid rounding noise, data with weight smaller than 1e-6
 249      *       is ignored.
 250      *
 251      * \param[in] weight               The weight of the data.
 252      * \param[in] data                 One data point for each grid dimension.
 253      * \param[in] blockLengthInWeight  If true, a block is measured in probability weight, otherwise
 254      * in time. \param[in] t                    The simulation time.
 255      */
 256     void addData(double weight, gmx::ArrayRef<const double> data, bool blockLengthInWeight, double t);
 257
 258     /*! \brief Returns an element of the time integrated correlation tensor at a given point in the grid.
 259      *
 260      * The units of the integral are time*(units of data)^2. This will be friction units time/length^2
 261      * if the data unit is 1/length.
 262      *
 263      * The correlation index lists the elements of the upper-triangular
 264      * correlation matrix row-wise, so e.g. in 3D:
 265      * 0 (0,0), 1 (1,0), 2 (1,1), 3 (2,0), 4 (2,1), 5 (2,2).
 266      * (TODO: this should ideally not have to be known by the caller.)
 267      *
 268      * \param[in] tensorIndex  Index in the tensor.
 269      * \param[in] dtSample     The sampling interval length.
 270      * \returns the integral.
 271      */
 272     double getTimeIntegral(int tensorIndex, double dtSample) const;
 273
 274     /*! \brief
 275      * Returns the volume element of the correlation metric.
 276      *
 277      * The matrix of the metric equals the time-integrated correlation matrix. The volume element of
 278      * the metric therefore equals the square-root of the absolute value of its determinant
 279      * according to the standard formula for a volume element in a metric space.
 280      *
 281      * Since the units of the metric matrix elements are time*(units of data)^2, the volume element
 282      * has units of (sqrt(time)*(units of data))^(ndim of data).
 283      *
 284      * \param[in] dtSample  The sampling interval length.
 285      * \returns the volume element.
 286      */
 287     double getVolumeElement(double dtSample) const;
 288
 289 private:
 290     std::vector<CorrelationBlockData> blockDataList_; /**< The block data for different, consecutively doubling block lengths. */
 291 };
 292
 293 } // namespace gmx
 294
 295 #endif /* GMX_AWH_CORRELATIONTENSOR_H */