src/gromacs/awh/correlationtensor.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2018,2019, 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  * \brief
  38  * Implements the CorrelationTensor class for correlation tensor statistics.
  39  *
  40  * \author Viveca Lindahl
  41  * \author Berk Hess <hess@kth.se>
  42  * \ingroup module_awh
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "correlationtensor.h"
  48
  49 #include <cassert>
  50
  51 #include "gromacs/math/functions.h"
  52 #include "gromacs/math/utilities.h"
  53 #include "gromacs/utility/gmxassert.h"
  54
  55 namespace gmx
  56 {
  57
  58 namespace
  59 {
  60
  61 /*! \brief
  62  * Get the block index for the current block and simulation length.
  63  *
  64  * This is simply how many blocks of length \p blockLength fit completely
  65  * into \p totalAccumulatedLength, which is either the current time minus
  66  * the start time, which time weighting, or the total weight.
  67  *
  68  * \param[in] blockLength               Block length.
  69  * \param[in] currentAccumulatedLength  Sampling length of all data collected up to the current
  70  * step, in time or weight. \returns the block index.
  71  */
  72 int getBlockIndex(double blockLength, double currentAccumulatedLength)
  73 {
  74     return static_cast<int>(currentAccumulatedLength / blockLength);
  75 }
  76
  77 } // namespace
  78
  79 double CorrelationTensor::getTimeIntegral(int tensorIndex, double dtSample) const
  80 {
  81     const CorrelationBlockData& blockData           = blockDataList_[0];
  82     double                      weight              = blockData.sumOverBlocksBlockSquareWeight();
  83     double                      correlationIntegral = 0;
  84     if (weight > 0)
  85     {
  86         correlationIntegral = blockData.correlationIntegral()[tensorIndex] / weight;
  87     }
  88
  89     return 0.5 * correlationIntegral * dtSample;
  90 }
  91
  92 double CorrelationTensor::getVolumeElement(double dtSample) const
  93 {
  94     double det;
  95
  96     switch (blockDataList_[0].correlationIntegral().size())
  97     {
  98         case 1:
  99             /* 1-dimensional tensor: [a] */
 100             det = getTimeIntegral(0, dtSample);
 101             break;
 102         case 3:
 103         {
 104             /* 2-dimensional tensor: [a b; b c] */
 105             double a = getTimeIntegral(0, dtSample);
 106             double b = getTimeIntegral(1, dtSample);
 107             double c = getTimeIntegral(2, dtSample);
 108
 109             det = a * c - b * b;
 110         }
 111         break;
 112         case 6:
 113         {
 114             /* 3-dimensional tensor: [a b d; b c e; d e f] */
 115             double a = getTimeIntegral(0, dtSample);
 116             double b = getTimeIntegral(1, dtSample);
 117             double c = getTimeIntegral(2, dtSample);
 118             double d = getTimeIntegral(3, dtSample);
 119             double e = getTimeIntegral(4, dtSample);
 120             double f = getTimeIntegral(5, dtSample);
 121
 122             det = a * c * f + 2 * b * d * e - d * c * d - b * b * f - a * e * e;
 123         }
 124         break;
 125         default:
 126             det = 0;
 127             /* meh */
 128             break;
 129     }
 130
 131     /* Returns 0 if no data, not supported number of dims
 132        or not enough data to give a positive determinant (as it should be) */
 133     return det > 0 ? std::sqrt(det) : 0;
 134 }
 135
 136 void CorrelationTensor::doubleBlockLengths()
 137 {
 138     /* We need to shift the data so that a given blockdata gets the data for double the block
 139        length. The data for the shortest block length is not needed anymore. */
 140
 141     for (size_t i = 0; i < blockDataList_.size() - 1; i++)
 142     {
 143         blockDataList_[i] = blockDataList_[i + 1];
 144     }
 145
 146     /* The blockdata which has 1 block is the same as the old one but with double the block length */
 147     blockDataList_.back().doubleBlockLength();
 148 }
 149
 150 void CorrelationTensor::updateBlockLengths(double samplingLength)
 151 {
 152     /* How many times do we need to double the longest block length to fit the data? */
 153     double longestLength = blockDataList_.back().blockLength();
 154     int    numDoublings  = 0;
 155     while (samplingLength > longestLength)
 156     {
 157         numDoublings++;
 158         longestLength *= 2;
 159     }
 160
 161     while (numDoublings > 0)
 162     {
 163         doubleBlockLengths();
 164         numDoublings--;
 165     }
 166 }
 167
 168 /* Adds a filled data block to correlation time integral. */
 169 void CorrelationBlockData::addBlockToCorrelationIntegral()
 170 {
 171     const bool firstBlock = (sumOverBlocksSquareBlockWeight_ == 0);
 172
 173     if (!firstBlock)
 174     {
 175         const int numDim      = coordData_.size();
 176         int       tensorIndex = 0;
 177         for (int d1 = 0; d1 < numDim; d1++)
 178         {
 179             const CoordData& c1 = coordData_[d1];
 180
 181             for (int d2 = 0; d2 <= d1; d2++)
 182             {
 183                 const CoordData& c2 = coordData_[d2];
 184
 185                 /* Compute change in correlaion integral due to adding
 186                  * the block through computing the difference of the block
 187                  * average with the old average for one component (we use x)
 188                  * and with the new component (we use y).
 189                  */
 190                 /* Need the old average, before the data of this block was added */
 191                 GMX_ASSERT(sumOverBlocksSquareBlockWeight_,
 192                            "Denominator should be > 0 (should be guaranteed by the conditional "
 193                            "above)");
 194                 double oldAverageX =
 195                         c1.sumOverBlocksBlockWeightBlockWeightX / sumOverBlocksSquareBlockWeight_;
 196
 197                 double newSumWeightBlockWeightY =
 198                         c2.sumOverBlocksBlockWeightBlockWeightX + blockSumWeight_ * c2.blockSumWeightX;
 199                 double newSumSquareBlockWeight =
 200                         sumOverBlocksSquareBlockWeight_ + gmx::square(blockSumWeight_);
 201
 202                 GMX_ASSERT(newSumSquareBlockWeight > 0, "Denominator should be > 0");
 203                 double newAverageY = newSumWeightBlockWeightY / newSumSquareBlockWeight;
 204
 205                 double diffBlockWithOldAverageX = c1.blockSumWeightX - oldAverageX * blockSumWeight_;
 206                 double diffBlockWithNewAverageY = c2.blockSumWeightX - newAverageY * blockSumWeight_;
 207
 208                 /* Update the correlation integral using the changes in averages. */
 209                 correlationIntegral_[tensorIndex] += diffBlockWithOldAverageX * diffBlockWithNewAverageY;
 210                 tensorIndex++;
 211             }
 212         }
 213     }
 214
 215     /* Add the weights of the block to the block sums and clear the weights */
 216     sumOverBlocksSquareBlockWeight_ += gmx::square(blockSumWeight_);
 217     sumOverBlocksBlockSquareWeight_ += blockSumSquareWeight_;
 218     for (auto& c : coordData_)
 219     {
 220         c.sumOverBlocksBlockWeightBlockWeightX += blockSumWeight_ * c.blockSumWeightX;
 221         /* Reset */
 222         c.blockSumWeightX = 0;
 223     }
 224     /* Reset */
 225     blockSumWeight_       = 0;
 226     blockSumSquareWeight_ = 0;
 227 }
 228
 229 void CorrelationTensor::addData(double weight, gmx::ArrayRef<const double> data, bool blockLengthInWeight, double t)
 230 {
 231     /* We should avoid adding data with very small weight to avoid
 232      * divergence close to 0/0. The total spread weight for each sample is 1,
 233      * so 1e-6 is a completely negligible amount.
 234      */
 235     if (weight < 1e-6)
 236     {
 237         /* Nothing to add */
 238         return;
 239     }
 240
 241     /*  The sampling length is measured either in the total (local) weight or the current time */
 242     double samplingLength = blockLengthInWeight ? getWeight() + weight : t;
 243
 244     /* Make sure the blocks are long enough to fit all data */
 245     updateBlockLengths(samplingLength);
 246
 247     /* Store the data for each block length considered. First update the longest block which has all
 248        data since it's used for updating the correlation function for the other block lengths. */
 249     for (size_t i = 0; i < blockDataList_.size() - 1; i++)
 250     {
 251         CorrelationBlockData& bd = blockDataList_[i];
 252
 253         /* Find current block index given the block length. */
 254         int blockIndex = getBlockIndex(bd.blockLength(), samplingLength);
 255
 256         if (bd.previousBlockIndex() >= 0 && bd.previousBlockIndex() != blockIndex)
 257         {
 258             /* Changed block. Update correlation with old data before adding to new block. */
 259             bd.addBlockToCorrelationIntegral();
 260         }
 261
 262         /* Keep track of which block index data was last added to */
 263         bd.setPreviousBlockIndex(blockIndex);
 264
 265         /* Store the data */
 266         bd.addData(weight, data);
 267     }
 268
 269     /* The last blockdata has only 1 block which contains all data so far.
 270      * Add the data for the largest block length.
 271      */
 272     blockDataList_.back().addData(weight, data);
 273 }
 274
 275 CorrelationTensor::CorrelationTensor(int numDim, int numBlockData, double blockLengthInit)
 276 {
 277     unsigned int scaling = 1;
 278
 279     GMX_RELEASE_ASSERT(numBlockData < static_cast<int>(sizeof(scaling) * 8),
 280                        "numBlockData should we smaller than the number of bits in scaling");
 281
 282     for (int n = 0; n < numBlockData; n++)
 283     {
 284         blockDataList_.emplace_back(numDim, scaling * blockLengthInit);
 285         scaling <<= 1; /* Double block length */
 286     }
 287 }
 288
 289 } // namespace gmx