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, 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 step, in time or weight.
  70  * \returns the block index.
  71  */
  72 int getBlockIndex(double blockLength,
  73                   double currentAccumulatedLength)
  74 {
  75     return static_cast<int>(currentAccumulatedLength/blockLength);
  76 }
  77
  78 }   // namespace
  79
  80 double CorrelationTensor::getTimeIntegral(int    tensorIndex,
  81                                           double dtSample) const
  82 {
  83     const CorrelationBlockData &blockData           = blockDataList_[0];
  84     double                      weight              = blockData.sumOverBlocksBlockSquareWeight();
  85     double                      correlationIntegral = 0;
  86     if (weight > 0)
  87     {
  88         correlationIntegral = blockData.correlationIntegral()[tensorIndex]/weight;
  89     }
  90
  91     return 0.5*correlationIntegral*dtSample;
  92 }
  93
  94 double CorrelationTensor::getVolumeElement(double dtSample) const
  95 {
  96     double det;
  97
  98     switch (blockDataList_[0].correlationIntegral().size())
  99     {
 100         case 1:
 101             /* 1-dimensional tensor: [a] */
 102             det = getTimeIntegral(0, dtSample);
 103             break;
 104         case 3:
 105         {
 106             /* 2-dimensional tensor: [a b; b c] */
 107             double a = getTimeIntegral(0, dtSample);
 108             double b = getTimeIntegral(1, dtSample);
 109             double c = getTimeIntegral(2, dtSample);
 110
 111             det      = a*c - b*b;
 112         }
 113         break;
 114         case 6:
 115         {
 116             /* 3-dimensional tensor: [a b d; b c e; d e f] */
 117             double a = getTimeIntegral(0, dtSample);
 118             double b = getTimeIntegral(1, dtSample);
 119             double c = getTimeIntegral(2, dtSample);
 120             double d = getTimeIntegral(3, dtSample);
 121             double e = getTimeIntegral(4, dtSample);
 122             double f = getTimeIntegral(5, dtSample);
 123
 124             det      = a*c*f + 2*b*d*e - d*c*d - b*b*f - a*e*e;
 125         }
 126         break;
 127         default:
 128             det = 0;
 129             /* meh */
 130             break;
 131     }
 132
 133     /* Returns 0 if no data, not supported number of dims
 134        or not enough data to give a positive determinant (as it should be) */
 135     return det > 0 ? std::sqrt(det) : 0;
 136 }
 137
 138 void CorrelationTensor::doubleBlockLengths()
 139 {
 140     /* We need to shift the data so that a given blockdata gets the data for double the block length.
 141        The data for the shortest block length is not needed anymore. */
 142
 143     for (size_t i = 0; i < blockDataList_.size() - 1; i++)
 144     {
 145         blockDataList_[i] = blockDataList_[i + 1];
 146     }
 147
 148     /* The blockdata which has 1 block is the same as the old one but with double the block length */
 149     blockDataList_.back().doubleBlockLength();
 150 }
 151
 152 void CorrelationTensor::updateBlockLengths(double samplingLength)
 153 {
 154     /* How many times do we need to double the longest block length to fit the data? */
 155     double longestLength = blockDataList_.back().blockLength();
 156     int    numDoublings  = 0;
 157     while (samplingLength > longestLength)
 158     {
 159         numDoublings++;
 160         longestLength *= 2;
 161     }
 162
 163     while (numDoublings > 0)
 164     {
 165         doubleBlockLengths();
 166         numDoublings--;
 167     }
 168 }
 169
 170 /* Adds a filled data block to correlation time integral. */
 171 void CorrelationBlockData::addBlockToCorrelationIntegral()
 172 {
 173     const bool firstBlock = (sumOverBlocksSquareBlockWeight_ == 0);
 174
 175     if (!firstBlock)
 176     {
 177         const int numDim      = coordData_.size();
 178         int       tensorIndex = 0;
 179         for (int d1 = 0; d1 < numDim; d1++)
 180         {
 181             const CoordData &c1 = coordData_[d1];
 182
 183             for (int d2 = 0; d2 <= d1; d2++)
 184             {
 185                 const CoordData &c2 = coordData_[d2];
 186
 187                 /* Compute change in correlaion integral due to adding
 188                  * the block through computing the difference of the block
 189                  * average with the old average for one component (we use x)
 190                  * and with the new component (we use y).
 191                  */
 192                 /* Need the old average, before the data of this block was added */
 193                 GMX_ASSERT(sumOverBlocksSquareBlockWeight_, "Denominator should be > 0 (should be guaranteed by the conditional above)");
 194                 double oldAverageX              = c1.sumOverBlocksBlockWeightBlockWeightX/sumOverBlocksSquareBlockWeight_;
 195
 196                 double newSumWeightBlockWeightY = c2.sumOverBlocksBlockWeightBlockWeightX + blockSumWeight_*c2.blockSumWeightX;
 197                 double newSumSquareBlockWeight  = sumOverBlocksSquareBlockWeight_ + gmx::square(blockSumWeight_);
 198
 199                 GMX_ASSERT(newSumSquareBlockWeight > 0, "Denominator should be > 0");
 200                 double newAverageY              = newSumWeightBlockWeightY/newSumSquareBlockWeight;
 201
 202                 double diffBlockWithOldAverageX = c1.blockSumWeightX - oldAverageX*blockSumWeight_;
 203                 double diffBlockWithNewAverageY = c2.blockSumWeightX - newAverageY*blockSumWeight_;
 204
 205                 /* Update the correlation integral using the changes in averages. */
 206                 correlationIntegral_[tensorIndex] += diffBlockWithOldAverageX*diffBlockWithNewAverageY;
 207                 tensorIndex++;
 208             }
 209         }
 210     }
 211
 212     /* Add the weights of the block to the block sums and clear the weights */
 213     sumOverBlocksSquareBlockWeight_ += gmx::square(blockSumWeight_);
 214     sumOverBlocksBlockSquareWeight_ += blockSumSquareWeight_;
 215     for (auto &c : coordData_)
 216     {
 217         c.sumOverBlocksBlockWeightBlockWeightX += blockSumWeight_*c.blockSumWeightX;
 218         /* Reset */
 219         c.blockSumWeightX                           = 0;
 220     }
 221     /* Reset */
 222     blockSumWeight_       = 0;
 223     blockSumSquareWeight_ = 0;
 224 }
 225
 226 void CorrelationTensor::addData(double                       weight,
 227                                 gmx::ArrayRef<const double>  data,
 228                                 bool                         blockLengthInWeight,
 229                                 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 data since it's
 248        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 &&
 257             bd.previousBlockIndex() != blockIndex)
 258         {
 259             /* Changed block. Update correlation with old data before adding to new block. */
 260             bd.addBlockToCorrelationIntegral();
 261         }
 262
 263         /* Keep track of which block index data was last added to */
 264         bd.setPreviousBlockIndex(blockIndex);
 265
 266         /* Store the data */
 267         bd.addData(weight, data);
 268     }
 269
 270     /* The last blockdata has only 1 block which contains all data so far.
 271      * Add the data for the largest block length.
 272      */
 273     blockDataList_.back().addData(weight, data);
 274 }
 275
 276 CorrelationTensor::CorrelationTensor(int    numDim,
 277                                      int    numBlockData,
 278                                      double blockLengthInit)
 279 {
 280     unsigned int scaling = 1;
 281
 282     GMX_RELEASE_ASSERT(numBlockData < static_cast<int>(sizeof(scaling)*8), "numBlockData should we smaller than the number of bits in scaling");
 283
 284     for (int n = 0; n < numBlockData; n++)
 285     {
 286         blockDataList_.emplace_back(numDim, scaling*blockLengthInit);
 287         scaling <<= 1; /* Double block length */
 288     }
 289 }
 290
 291 } // namespace gmx