2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
39 * Implements helper functions for checkpointing the AWH state and observables history.
41 * \author Viveca Lindahl
42 * \author Berk Hess <hess@kth.se>
48 #include "correlationhistory.h"
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"
59 #include "correlationgrid.h"
64 void initCorrelationGridHistory(CorrelationGridHistory* correlationGridHistory,
65 int numCorrelationTensors,
67 int blockDataListSize)
69 correlationGridHistory->numCorrelationTensors = numCorrelationTensors;
70 correlationGridHistory->tensorSize = tensorSize;
71 correlationGridHistory->blockDataListSize = blockDataListSize;
73 correlationGridHistory->blockDataBuffer.resize(numCorrelationTensors * tensorSize * blockDataListSize);
76 CorrelationGridHistory initCorrelationGridHistoryFromState(const CorrelationGrid& correlationGrid)
78 CorrelationGridHistory correlationGridHistory;
80 initCorrelationGridHistory(&correlationGridHistory,
81 correlationGrid.tensors().size(),
82 correlationGrid.tensorSize(),
83 correlationGrid.blockDataListSize());
85 return correlationGridHistory;
88 /* Update the correlation grid history for checkpointing. */
89 void updateCorrelationGridHistory(CorrelationGridHistory* correlationGridHistory,
90 const CorrelationGrid& correlationGrid)
92 GMX_RELEASE_ASSERT(correlationGridHistory != nullptr, "We need a valid history object");
94 gmx::ArrayRef<CorrelationBlockDataHistory> blockDataBuffer = correlationGridHistory->blockDataBuffer;
96 /* Store the grid in a linear array */
97 gmx::index bufferIndex = 0;
98 for (const CorrelationTensor& tensor : correlationGrid.tensors())
100 const int numDims = tensor.blockDataList()[0].coordData().size();
101 const int tensorSize = tensor.blockDataList()[0].correlationIntegral().size();
103 /* Loop of the tensor elements, ignore the symmetric data */
106 for (int k = 0; k < tensorSize; k++)
108 /* BlockData for each correlation element */
109 for (const CorrelationBlockData& blockData : tensor.blockDataList())
111 const CorrelationBlockData::CoordData& cx = blockData.coordData()[d1];
112 const CorrelationBlockData::CoordData& cy = blockData.coordData()[d2];
114 CorrelationBlockDataHistory& bdh = blockDataBuffer[bufferIndex];
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];
140 GMX_RELEASE_ASSERT(bufferIndex == blockDataBuffer.ssize(),
141 "We should store exactly as many elements as the buffer size");
144 void CorrelationBlockData::restoreFromHistory(const CorrelationBlockDataHistory& blockHistory,
145 const std::vector<CorrelationBlockData::CoordData>& coordData,
146 const std::vector<double>& correlationIntegral)
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;
158 /* Restore a correlation element from history. */
159 void CorrelationTensor::restoreFromHistory(const std::vector<CorrelationBlockDataHistory>& blockDataBuffer,
162 /* Blockdata for each correlation element */
163 for (CorrelationBlockData& blockData : blockDataList_)
165 /* Correlation elements for each tensor */
166 const int numDims = blockDataList_[0].coordData().size();
167 const int tensorSize = blockDataList_[0].correlationIntegral().size();
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++)
176 if (*bufferIndex >= blockDataBuffer.size())
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."));
183 const CorrelationBlockDataHistory& blockHistory = blockDataBuffer[*bufferIndex];
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.
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;
197 correlationIntegral[k] = blockHistory.correlationIntegral;
199 /* Check if we collected all data needed for blockData */
200 if (k == tensorSize - 1)
202 blockData.restoreFromHistory(blockHistory, coordData, correlationIntegral);
217 /* Restores the correlation grid state from the correlation grid history. */
218 void CorrelationGrid::restoreStateFromHistory(const CorrelationGridHistory& correlationGridHistory)
220 if (tensors_.size() != static_cast<size_t>(correlationGridHistory.numCorrelationTensors))
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."));
227 /* Extract the state from the linear history array */
228 size_t bufferIndex = 0;
229 for (CorrelationTensor& tensor : tensors_)
231 tensor.restoreFromHistory(correlationGridHistory.blockDataBuffer, &bufferIndex);
234 if (bufferIndex != correlationGridHistory.blockDataBuffer.size())
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."));