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, 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/math/utilities.h"
  50 #include "gromacs/utility/gmxassert.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),
  89     blockLengthMeasure(blockLengthMeasure)
  90 {
  91     /* Set the initial block length for the block averaging. The length doesn't really matter
  92        after the block length has been doubled a few times, as long as it's set small enough */
  93     if (blockLengthMeasure == BlockLengthMeasure::Weight)
  94     {
  95         blockLengthInit = blockLengthInit > 0 ? blockLengthInit : 1;
  96     }
  97     else
  98     {
  99         blockLengthInit = blockLengthInit > 0 ? blockLengthInit : dtSample;
 100     }
 101
 102     /* Set the number of blocks. The number of blocks determines the current span of the data
 103        and how many different block lengths (nblockdata) we need to keep track of to be able to
 104        increase the block length later */
 105     int numBlocks         = CorrelationTensor::c_numCorrelationBlocks;
 106     int BlockDataListSize = getBlockDataListSize(numBlocks);
 107
 108     tensors_.resize(numPoints, CorrelationTensor(numDim, BlockDataListSize, blockLengthInit));
 109 }
 110
 111 int CorrelationGrid::getNumBlocks() const
 112 {
 113     const auto& blockDataList  = tensors()[0].blockDataList();
 114     double      maxBlockLength = blockDataList.back().blockLength();
 115     double      minBlockLength = blockDataList[0].blockLength();
 116
 117     /* 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) */
 118     if (maxBlockLength < GMX_DOUBLE_MAX)
 119     {
 120         return static_cast<int>(maxBlockLength / minBlockLength);
 121     }
 122     else
 123     {
 124         return -1;
 125     }
 126 }
 127
 128 double CorrelationGrid::getBlockLength() const
 129 {
 130     /* Return the  minimum blocklength */
 131     return tensors()[0].blockDataList()[0].blockLength();
 132 }
 133
 134 } // namespace gmx