src/gromacs/applied_forces/awh/correlationhistory.cpp

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