2 * This file is part of the GROMACS molecular simulation package.
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.
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, correlationGrid.tensors().size(),
81 correlationGrid.tensorSize(), correlationGrid.blockDataListSize());
83 return correlationGridHistory;
86 /* Update the correlation grid history for checkpointing. */
87 void updateCorrelationGridHistory(CorrelationGridHistory* correlationGridHistory,
88 const CorrelationGrid& correlationGrid)
90 GMX_RELEASE_ASSERT(correlationGridHistory != nullptr, "We need a valid history object");
92 gmx::ArrayRef<CorrelationBlockDataHistory> blockDataBuffer = correlationGridHistory->blockDataBuffer;
94 /* Store the grid in a linear array */
95 gmx::index bufferIndex = 0;
96 for (const CorrelationTensor& tensor : correlationGrid.tensors())
98 const int numDims = tensor.blockDataList()[0].coordData().size();
99 const int tensorSize = tensor.blockDataList()[0].correlationIntegral().size();
101 /* Loop of the tensor elements, ignore the symmetric data */
104 for (int k = 0; k < tensorSize; k++)
106 /* BlockData for each correlation element */
107 for (const CorrelationBlockData& blockData : tensor.blockDataList())
109 const CorrelationBlockData::CoordData& cx = blockData.coordData()[d1];
110 const CorrelationBlockData::CoordData& cy = blockData.coordData()[d2];
112 CorrelationBlockDataHistory& bdh = blockDataBuffer[bufferIndex];
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];
138 GMX_RELEASE_ASSERT(bufferIndex == blockDataBuffer.ssize(),
139 "We should store exactly as many elements as the buffer size");
142 void CorrelationBlockData::restoreFromHistory(const CorrelationBlockDataHistory& blockHistory,
143 const std::vector<CorrelationBlockData::CoordData>& coordData,
144 const std::vector<double>& correlationIntegral)
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;
156 /* Restore a correlation element from history. */
157 void CorrelationTensor::restoreFromHistory(const std::vector<CorrelationBlockDataHistory>& blockDataBuffer,
160 /* Blockdata for each correlation element */
161 for (CorrelationBlockData& blockData : blockDataList_)
163 /* Correlation elements for each tensor */
164 const int numDims = blockDataList_[0].coordData().size();
165 const int tensorSize = blockDataList_[0].correlationIntegral().size();
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++)
174 if (*bufferIndex >= blockDataBuffer.size())
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."));
181 const CorrelationBlockDataHistory& blockHistory = blockDataBuffer[*bufferIndex];
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.
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;
195 correlationIntegral[k] = blockHistory.correlationIntegral;
197 /* Check if we collected all data needed for blockData */
198 if (k == tensorSize - 1)
200 blockData.restoreFromHistory(blockHistory, coordData, correlationIntegral);
215 /* Restores the correlation grid state from the correlation grid history. */
216 void CorrelationGrid::restoreStateFromHistory(const CorrelationGridHistory& correlationGridHistory)
218 if (tensors_.size() != static_cast<size_t>(correlationGridHistory.numCorrelationTensors))
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."));
225 /* Extract the state from the linear history array */
226 size_t bufferIndex = 0;
227 for (CorrelationTensor& tensor : tensors_)
229 tensor.restoreFromHistory(correlationGridHistory.blockDataBuffer, &bufferIndex);
232 if (bufferIndex != correlationGridHistory.blockDataBuffer.size())
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."));