src/gromacs/applied_forces/awh/correlationgrid.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2019,2020,2021, 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 CorrelationGrid class to collect correlation statistics on a grid, using several block lengths.
  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 "correlationgrid.h"
  48
  49 #include "gromacs/utility/gmxassert.h"
  50 #include "gromacs/utility/real.h"
  51
  52 namespace gmx
  53 {
  54
  55 namespace
  56 {
  57
  58 /*! \brief
  59  * Return the number of block data structs needed for keeping a certain number of blocks.
  60  *
  61  * The list start with 1 block and doubles, so we need 1 + 2log(numBlocks).
  62  *
  63  * \param[in] numBlocks  Number of blocks.
  64  * \returns the number of block data structs.
  65  */
  66 int getBlockDataListSize(int numBlocks)
  67 {
  68     int blockDataListSize = 1;
  69
  70     while (numBlocks > (1 << (blockDataListSize - 1)))
  71     {
  72         blockDataListSize++;
  73     }
  74
  75     GMX_RELEASE_ASSERT((1 << (blockDataListSize - 1)) == numBlocks,
  76                        "numBlocks should be a power of 2");
  77
  78     return blockDataListSize;
  79 }
  80
  81 } // namespace
  82
  83 CorrelationGrid::CorrelationGrid(int                numPoints,
  84                                  int                numDim,
  85                                  double             blockLengthInit,
  86                                  BlockLengthMeasure blockLengthMeasure,
  87                                  double             dtSample) :
  88     dtSample(dtSample), blockLengthMeasure(blockLengthMeasure)
  89 {
  90     /* Set the initial block length for the block averaging. The length doesn't really matter
  91        after the block length has been doubled a few times, as long as it's set small enough */
  92     if (blockLengthMeasure == BlockLengthMeasure::Weight)
  93     {
  94         blockLengthInit = blockLengthInit > 0 ? blockLengthInit : 1;
  95     }
  96     else
  97     {
  98         blockLengthInit = blockLengthInit > 0 ? blockLengthInit : dtSample;
  99     }
 100
 101     /* Set the number of blocks. The number of blocks determines the current span of the data
 102        and how many different block lengths (nblockdata) we need to keep track of to be able to
 103        increase the block length later */
 104     int numBlocks         = CorrelationTensor::c_numCorrelationBlocks;
 105     int BlockDataListSize = getBlockDataListSize(numBlocks);
 106
 107     tensors_.resize(numPoints, CorrelationTensor(numDim, BlockDataListSize, blockLengthInit));
 108 }
 109
 110 int CorrelationGrid::getNumBlocks() const
 111 {
 112     const auto& blockDataList  = tensors()[0].blockDataList();
 113     double      maxBlockLength = blockDataList.back().blockLength();
 114     double      minBlockLength = blockDataList[0].blockLength();
 115
 116     /* If we have a finite block span we have a constant number of blocks, otherwise we are always adding more blocks (and we don't keep track of the number) */
 117     if (maxBlockLength < GMX_DOUBLE_MAX)
 118     {
 119         return static_cast<int>(maxBlockLength / minBlockLength);
 120     }
 121     else
 122     {
 123         return -1;
 124     }
 125 }
 126
 127 double CorrelationGrid::getBlockLength() const
 128 {
 129     /* Return the  minimum blocklength */
 130     return tensors()[0].blockDataList()[0].blockLength();
 131 }
 132
 133 } // namespace gmx