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(),
  81                                correlationGrid.tensorSize(), correlationGrid.blockDataListSize());
  82
  83     return correlationGridHistory;
  84 }
  85
  86 /* Update the correlation grid history for checkpointing. */
  87 void updateCorrelationGridHistory(CorrelationGridHistory* correlationGridHistory,
  88                                   const CorrelationGrid&  correlationGrid)
  89 {
  90     GMX_RELEASE_ASSERT(correlationGridHistory != nullptr, "We need a valid history object");
  91
  92     gmx::ArrayRef<CorrelationBlockDataHistory> blockDataBuffer = correlationGridHistory->blockDataBuffer;
  93
  94     /* Store the grid in a linear array */
  95     gmx::index bufferIndex = 0;
  96     for (const CorrelationTensor& tensor : correlationGrid.tensors())
  97     {
  98         const int numDims    = tensor.blockDataList()[0].coordData().size();
  99         const int tensorSize = tensor.blockDataList()[0].correlationIntegral().size();
 100
 101         /* Loop of the tensor elements, ignore the symmetric data */
 102         int d1 = 0;
 103         int d2 = 0;
 104         for (int k = 0; k < tensorSize; k++)
 105         {
 106             /* BlockData for each correlation element */
 107             for (const CorrelationBlockData& blockData : tensor.blockDataList())
 108             {
 109                 const CorrelationBlockData::CoordData& cx = blockData.coordData()[d1];
 110                 const CorrelationBlockData::CoordData& cy = blockData.coordData()[d2];
 111
 112                 CorrelationBlockDataHistory& bdh = blockDataBuffer[bufferIndex];
 113
 114                 bdh.blockSumWeight                 = blockData.blockSumWeight();
 115                 bdh.blockSumSquareWeight           = blockData.blockSumSquareWeight();
 116                 bdh.blockSumWeightX                = cx.blockSumWeightX;
 117                 bdh.blockSumWeightY                = cy.blockSumWeightX;
 118                 bdh.sumOverBlocksSquareBlockWeight = blockData.sumOverBlocksSquareBlockWeight();
 119                 bdh.sumOverBlocksBlockSquareWeight = blockData.sumOverBlocksBlockSquareWeight();
 120                 bdh.sumOverBlocksBlockWeightBlockWeightX = cx.sumOverBlocksBlockWeightBlockWeightX;
 121                 bdh.sumOverBlocksBlockWeightBlockWeightY = cy.sumOverBlocksBlockWeightBlockWeightX;
 122                 bdh.previousBlockIndex                   = blockData.previousBlockIndex();
 123                 bdh.blockLength                          = blockData.blockLength();
 124                 bdh.correlationIntegral                  = blockData.correlationIntegral()[k];
 125
 126                 bufferIndex++;
 127             }
 128
 129             d1++;
 130             if (d1 == numDims)
 131             {
 132                 d1 = 0;
 133                 d2++;
 134             }
 135         }
 136     }
 137
 138     GMX_RELEASE_ASSERT(bufferIndex == blockDataBuffer.ssize(),
 139                        "We should store exactly as many elements as the buffer size");
 140 }
 141
 142 void CorrelationBlockData::restoreFromHistory(const CorrelationBlockDataHistory& blockHistory,
 143                                               const std::vector<CorrelationBlockData::CoordData>& coordData,
 144                                               const std::vector<double>& correlationIntegral)
 145 {
 146     blockSumWeight_                 = blockHistory.blockSumWeight;
 147     blockSumSquareWeight_           = blockHistory.blockSumSquareWeight;
 148     sumOverBlocksSquareBlockWeight_ = blockHistory.sumOverBlocksSquareBlockWeight;
 149     sumOverBlocksBlockSquareWeight_ = blockHistory.sumOverBlocksBlockSquareWeight;
 150     previousBlockIndex_             = blockHistory.previousBlockIndex;
 151     blockLength_                    = blockHistory.blockLength;
 152     coordData_                      = coordData;
 153     correlationIntegral_            = correlationIntegral;
 154 }
 155
 156 /* Restore a correlation element from history. */
 157 void CorrelationTensor::restoreFromHistory(const std::vector<CorrelationBlockDataHistory>& blockDataBuffer,
 158                                            size_t* bufferIndex)
 159 {
 160     /* Blockdata for each correlation element */
 161     for (CorrelationBlockData& blockData : blockDataList_)
 162     {
 163         /* Correlation elements for each tensor */
 164         const int numDims    = blockDataList_[0].coordData().size();
 165         const int tensorSize = blockDataList_[0].correlationIntegral().size();
 166         int       d1         = 0;
 167         int       d2         = 0;
 168
 169         /* Temporary containers to collect data */
 170         std::vector<CorrelationBlockData::CoordData> coordData(numDims);
 171         std::vector<double>                          correlationIntegral(tensorSize);
 172         for (int k = 0; k < tensorSize; k++)
 173         {
 174             if (*bufferIndex >= blockDataBuffer.size())
 175             {
 176                 GMX_THROW(InvalidInputError(
 177                         "Mismatch of the correlation tensor size for the force correlation between "
 178                         "checkpoint and simulation. Likely you have provided a checkpoint from a "
 179                         "different simulation."));
 180             }
 181             const CorrelationBlockDataHistory& blockHistory = blockDataBuffer[*bufferIndex];
 182
 183             /* To simplify the checkpointing, CorrelationBlockDataHistory
 184              * duplicates some weight data for all tensor elements.
 185              * Here we collect the coordinate and tensor data
 186              * in temporary buffers.
 187              */
 188             coordData[d1].blockSumWeightX = blockHistory.blockSumWeightX;
 189             coordData[d2].blockSumWeightX = blockHistory.blockSumWeightY;
 190             coordData[d1].sumOverBlocksBlockWeightBlockWeightX =
 191                     blockHistory.sumOverBlocksBlockWeightBlockWeightX;
 192             coordData[d2].sumOverBlocksBlockWeightBlockWeightX =
 193                     blockHistory.sumOverBlocksBlockWeightBlockWeightY;
 194
 195             correlationIntegral[k] = blockHistory.correlationIntegral;
 196
 197             /* Check if we collected all data needed for blockData */
 198             if (k == tensorSize - 1)
 199             {
 200                 blockData.restoreFromHistory(blockHistory, coordData, correlationIntegral);
 201             }
 202
 203             (*bufferIndex)++;
 204
 205             d1++;
 206             if (d1 == numDims)
 207             {
 208                 d1 = 0;
 209                 d2++;
 210             }
 211         }
 212     }
 213 }
 214
 215 /* Restores the correlation grid state from the correlation grid history. */
 216 void CorrelationGrid::restoreStateFromHistory(const CorrelationGridHistory& correlationGridHistory)
 217 {
 218     if (tensors_.size() != static_cast<size_t>(correlationGridHistory.numCorrelationTensors))
 219     {
 220         GMX_THROW(InvalidInputError(
 221                 "Mismatch of the grid size for the force correlation between checkpoint and "
 222                 "simulation. Likely you have provided a checkpoint from a different simulation."));
 223     }
 224
 225     /* Extract the state from the linear history array */
 226     size_t bufferIndex = 0;
 227     for (CorrelationTensor& tensor : tensors_)
 228     {
 229         tensor.restoreFromHistory(correlationGridHistory.blockDataBuffer, &bufferIndex);
 230     }
 231
 232     if (bufferIndex != correlationGridHistory.blockDataBuffer.size())
 233     {
 234         GMX_THROW(
 235                 InvalidInputError("Mismatch of the correlation tensor size for the force "
 236                                   "correlation between checkpoint and simulation. Likely you have "
 237                                   "provided a checkpoint from a different simulation."));
 238     }
 239 }
 240
 241 } // namespace gmx