src/gromacs/awh/correlationhistory.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  *
  38  * \brief
  39  * Implements helper functions for checkpointing the AWH state and observables history.
  40  *
  41  * \author Viveca Lindahl
  42  * \author Berk Hess <hess@kth.se>
  43  * \ingroup module_awh
  44  */
  45
  46 #include "gmxpre.h"
  47
  48 #include "correlationhistory.h"
  49
  50 #include <cassert>
  51
  52 #include "gromacs/gmxlib/network.h"
  53 #include "gromacs/mdtypes/awh_correlation_history.h"
  54 #include "gromacs/mdtypes/commrec.h"
  55 #include "gromacs/utility/arrayref.h"
  56 #include "gromacs/utility/exceptions.h"
  57 #include "gromacs/utility/gmxassert.h"
  58
  59 #include "correlationgrid.h"
  60
  61 namespace gmx
  62 {
  63
  64 void initCorrelationGridHistory(CorrelationGridHistory *correlationGridHistory,
  65                                 int                     numCorrelationTensors,
  66                                 int                     tensorSize,
  67                                 int                     blockDataListSize)
  68 {
  69     correlationGridHistory->numCorrelationTensors = numCorrelationTensors;
  70     correlationGridHistory->tensorSize            = tensorSize;
  71     correlationGridHistory->blockDataListSize     = blockDataListSize;
  72
  73     correlationGridHistory->blockDataBuffer.resize(numCorrelationTensors*tensorSize*blockDataListSize);
  74 }
  75
  76 CorrelationGridHistory initCorrelationGridHistoryFromState(const CorrelationGrid &correlationGrid)
  77 {
  78     CorrelationGridHistory correlationGridHistory;
  79
  80     initCorrelationGridHistory(&correlationGridHistory, correlationGrid.tensors().size(), correlationGrid.tensorSize(), correlationGrid.blockDataListSize());
  81
  82     return correlationGridHistory;
  83 }
  84
  85 /* Update the correlation grid history for checkpointing. */
  86 void updateCorrelationGridHistory(CorrelationGridHistory *correlationGridHistory,
  87                                   const CorrelationGrid  &correlationGrid)
  88 {
  89     GMX_RELEASE_ASSERT(correlationGridHistory != nullptr, "We need a valid history object");
  90
  91     gmx::ArrayRef<CorrelationBlockDataHistory> blockDataBuffer = correlationGridHistory->blockDataBuffer;
  92
  93     /* Store the grid in a linear array */
  94     gmx::index bufferIndex = 0;
  95     for (const CorrelationTensor &tensor : correlationGrid.tensors())
  96     {
  97         const int                numDims    = tensor.blockDataList()[0].coordData().size();
  98         const int                tensorSize = tensor.blockDataList()[0].correlationIntegral().size();
  99
 100         /* Loop of the tensor elements, ignore the symmetric data */
 101         int       d1         = 0;
 102         int       d2         = 0;
 103         for (int k = 0; k < tensorSize; k++)
 104         {
 105             /* BlockData for each correlation element */
 106             for (const CorrelationBlockData &blockData : tensor.blockDataList())
 107             {
 108                 const CorrelationBlockData::CoordData &cx  = blockData.coordData()[d1];
 109                 const CorrelationBlockData::CoordData &cy  = blockData.coordData()[d2];
 110
 111                 CorrelationBlockDataHistory           &bdh = blockDataBuffer[bufferIndex];
 112
 113                 bdh.blockSumWeight                       = blockData.blockSumWeight();
 114                 bdh.blockSumSquareWeight                 = blockData.blockSumSquareWeight();
 115                 bdh.blockSumWeightX                      = cx.blockSumWeightX;
 116                 bdh.blockSumWeightY                      = cy.blockSumWeightX;
 117                 bdh.sumOverBlocksSquareBlockWeight       = blockData.sumOverBlocksSquareBlockWeight();
 118                 bdh.sumOverBlocksBlockSquareWeight       = blockData.sumOverBlocksBlockSquareWeight();
 119                 bdh.sumOverBlocksBlockWeightBlockWeightX = cx.sumOverBlocksBlockWeightBlockWeightX;
 120                 bdh.sumOverBlocksBlockWeightBlockWeightY = cy.sumOverBlocksBlockWeightBlockWeightX;
 121                 bdh.previousBlockIndex                   = blockData.previousBlockIndex();
 122                 bdh.blockLength                          = blockData.blockLength();
 123                 bdh.correlationIntegral                  = blockData.correlationIntegral()[k];
 124
 125                 bufferIndex++;
 126             }
 127
 128             d1++;
 129             if (d1 == numDims)
 130             {
 131                 d1 = 0;
 132                 d2++;
 133             }
 134         }
 135     }
 136
 137     GMX_RELEASE_ASSERT(bufferIndex == blockDataBuffer.ssize(), "We should store exactly as many elements as the buffer size");
 138 }
 139
 140 void CorrelationBlockData::restoreFromHistory(const CorrelationBlockDataHistory                  &blockHistory,
 141                                               const std::vector<CorrelationBlockData::CoordData> &coordData,
 142                                               const std::vector<double>                          &correlationIntegral)
 143 {
 144     blockSumWeight_                 = blockHistory.blockSumWeight;
 145     blockSumSquareWeight_           = blockHistory.blockSumSquareWeight;
 146     sumOverBlocksSquareBlockWeight_ = blockHistory.sumOverBlocksSquareBlockWeight;
 147     sumOverBlocksBlockSquareWeight_ = blockHistory.sumOverBlocksBlockSquareWeight;
 148     previousBlockIndex_             = blockHistory.previousBlockIndex;
 149     blockLength_                    = blockHistory.blockLength;
 150     coordData_                      = coordData;
 151     correlationIntegral_            = correlationIntegral;
 152 }
 153
 154 /* Restore a correlation element from history. */
 155 void CorrelationTensor::restoreFromHistory(const std::vector<CorrelationBlockDataHistory> &blockDataBuffer,
 156                                            size_t                                         *bufferIndex)
 157 {
 158     /* Blockdata for each correlation element */
 159     for (CorrelationBlockData &blockData : blockDataList_)
 160     {
 161         /* Correlation elements for each tensor */
 162         const int numDims    = blockDataList_[0].coordData().size();
 163         const int tensorSize = blockDataList_[0].correlationIntegral().size();
 164         int       d1         = 0;
 165         int       d2         = 0;
 166
 167         /* Temporary containers to collect data */
 168         std::vector<CorrelationBlockData::CoordData> coordData(numDims);
 169         std::vector<double>                          correlationIntegral(tensorSize);
 170         for (int k = 0; k < tensorSize; k++)
 171         {
 172             if (*bufferIndex >= blockDataBuffer.size())
 173             {
 174                 GMX_THROW(InvalidInputError("Mismatch of the correlation tensor size for the force correlation between checkpoint and simulation. Likely you have provided a checkpoint from a different simulation."));
 175             }
 176             const CorrelationBlockDataHistory &blockHistory = blockDataBuffer[*bufferIndex];
 177
 178             /* To simplify the checkpointing, CorrelationBlockDataHistory
 179              * duplicates some weight data for all tensor elements.
 180              * Here we collect the coordinate and tensor data
 181              * in temporary buffers.
 182              */
 183             coordData[d1].blockSumWeightX                      = blockHistory.blockSumWeightX;
 184             coordData[d2].blockSumWeightX                      = blockHistory.blockSumWeightY;
 185             coordData[d1].sumOverBlocksBlockWeightBlockWeightX = blockHistory.sumOverBlocksBlockWeightBlockWeightX;
 186             coordData[d2].sumOverBlocksBlockWeightBlockWeightX = blockHistory.sumOverBlocksBlockWeightBlockWeightY;
 187
 188             correlationIntegral[k] = blockHistory.correlationIntegral;
 189
 190             /* Check if we collected all data needed for blockData */
 191             if (k == tensorSize - 1)
 192             {
 193                 blockData.restoreFromHistory(blockHistory, coordData, correlationIntegral);
 194             }
 195
 196             (*bufferIndex)++;
 197
 198             d1++;
 199             if (d1 == numDims)
 200             {
 201                 d1 = 0;
 202                 d2++;
 203             }
 204         }
 205     }
 206 }
 207
 208 /* Restores the correlation grid state from the correlation grid history. */
 209 void CorrelationGrid::restoreStateFromHistory(const CorrelationGridHistory &correlationGridHistory)
 210 {
 211     if (tensors_.size() != static_cast<size_t>(correlationGridHistory.numCorrelationTensors))
 212     {
 213         GMX_THROW(InvalidInputError("Mismatch of the grid size for the force correlation between checkpoint and simulation. Likely you have provided a checkpoint from a different simulation."));
 214     }
 215
 216     /* Extract the state from the linear history array */
 217     size_t bufferIndex = 0;
 218     for (CorrelationTensor &tensor : tensors_)
 219     {
 220         tensor.restoreFromHistory(correlationGridHistory.blockDataBuffer,
 221                                   &bufferIndex);
 222     }
 223
 224     if (bufferIndex != correlationGridHistory.blockDataBuffer.size())
 225     {
 226         GMX_THROW(InvalidInputError("Mismatch of the correlation tensor size for the force correlation between checkpoint and simulation. Likely you have provided a checkpoint from a different simulation."));
 227     }
 228 }
 229
 230 } // namespace gmx