2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 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.
37 * Implements PME force gathering tests.
39 * \author Aleksei Iupinov <a.yupinov@gmail.com>
40 * \ingroup module_ewald
47 #include <gmock/gmock.h>
49 #include "gromacs/mdtypes/inputrec.h"
50 #include "gromacs/utility/stringutil.h"
52 #include "testutils/refdata.h"
53 #include "testutils/testasserts.h"
55 #include "pmetestcommon.h"
64 /* Valid input instances */
66 //! A couple of valid inputs for boxes.
67 std::vector<Matrix3x3> const c_sampleBoxes{
69 Matrix3x3{ { 8.0F, 0.0F, 0.0F, 0.0F, 3.4F, 0.0F, 0.0F, 0.0F, 2.0F } },
71 Matrix3x3{ { 7.0F, 0.0F, 0.0F, 0.0F, 4.1F, 0.0F, 3.5F, 2.0F, 12.2F } },
74 //! A couple of valid inputs for grid sizes
75 std::vector<IVec> const c_sampleGridSizes{ IVec{ 16, 12, 14 }, IVec{ 13, 15, 11 } };
77 std::vector<real> const c_sampleChargesFull{ 4.95F, 3.11F, 3.97F, 1.08F, 2.09F, 1.1F,
78 4.13F, 3.31F, 2.8F, 5.83F, 5.09F, 6.1F,
79 2.86F, 0.24F, 5.76F, 5.19F, 0.72F };
81 //! All the input atom gridline indices
82 std::vector<IVec> const c_sampleGridLineIndicesFull{
83 IVec{ 4, 2, 6 }, IVec{ 1, 4, 10 }, IVec{ 0, 6, 6 }, IVec{ 0, 1, 4 }, IVec{ 6, 3, 0 },
84 IVec{ 7, 2, 2 }, IVec{ 8, 3, 1 }, IVec{ 4, 0, 3 }, IVec{ 0, 0, 0 }, IVec{ 8, 5, 8 },
85 IVec{ 4, 4, 2 }, IVec{ 7, 1, 7 }, IVec{ 8, 5, 5 }, IVec{ 2, 6, 5 }, IVec{ 1, 6, 2 },
86 IVec{ 7, 1, 8 }, IVec{ 3, 5, 1 },
89 // Spline values/derivatives below are also generated randomly, so they are bogus,
90 // but that should not affect the reproducibility, which we're after
92 //! A lot of bogus input spline values - should have at list (max PME order = 5) * (DIM = 3) * (total unique atom number in all test cases = 16) values
93 std::vector<real> const c_sampleSplineValuesFull{
94 0.12F, 0.81F, 0.29F, 0.22F, 0.13F, 0.19F, 0.12F, 0.8F, 0.44F, 0.38F, 0.32F, 0.36F, 0.27F,
95 0.11F, 0.17F, 0.94F, 0.07F, 0.9F, 0.98F, 0.96F, 0.07F, 0.94F, 0.77F, 0.24F, 0.84F, 0.16F,
96 0.77F, 0.57F, 0.52F, 0.27F, 0.39F, 0.45F, 0.6F, 0.59F, 0.44F, 0.91F, 0.97F, 0.43F, 0.24F,
97 0.52F, 0.73F, 0.55F, 0.99F, 0.39F, 0.97F, 0.35F, 0.1F, 0.68F, 0.19F, 0.1F, 0.77F, 0.2F,
98 0.43F, 0.69F, 0.76F, 0.32F, 0.31F, 0.94F, 0.53F, 0.6F, 0.93F, 0.57F, 0.94F, 0.88F, 0.75F,
99 0.77F, 0.91F, 0.72F, 0.07F, 0.78F, 0.09F, 0.02F, 0.48F, 0.97F, 0.89F, 0.39F, 0.48F, 0.19F,
100 0.02F, 0.92F, 0.8F, 0.41F, 0.53F, 0.32F, 0.38F, 0.58F, 0.36F, 0.46F, 0.92F, 0.91F, 0.01F,
101 0.86F, 0.54F, 0.86F, 0.94F, 0.37F, 0.35F, 0.81F, 0.89F, 0.48F, 0.34F, 0.18F, 0.11F, 0.02F,
102 0.87F, 0.95F, 0.66F, 0.67F, 0.38F, 0.45F, 0.04F, 0.94F, 0.54F, 0.76F, 0.58F, 0.83F, 0.31F,
103 0.73F, 0.71F, 0.06F, 0.35F, 0.32F, 0.35F, 0.61F, 0.27F, 0.98F, 0.83F, 0.11F, 0.3F, 0.42F,
104 0.95F, 0.69F, 0.58F, 0.29F, 0.1F, 0.68F, 0.94F, 0.62F, 0.51F, 0.47F, 0.04F, 0.47F, 0.34F,
105 0.71F, 0.52F, 0.19F, 0.69F, 0.5F, 0.59F, 0.05F, 0.74F, 0.11F, 0.4F, 0.81F, 0.24F, 0.53F,
106 0.71F, 0.07F, 0.17F, 0.41F, 0.23F, 0.78F, 0.27F, 0.1F, 0.71F, 0.36F, 0.67F, 0.6F, 0.94F,
107 0.69F, 0.19F, 0.58F, 0.68F, 0.5F, 0.62F, 0.38F, 0.29F, 0.44F, 0.04F, 0.89F, 0.0F, 0.76F,
108 0.22F, 0.16F, 0.08F, 0.62F, 0.51F, 0.62F, 0.83F, 0.72F, 0.96F, 0.99F, 0.4F, 0.79F, 0.83F,
109 0.21F, 0.43F, 0.32F, 0.44F, 0.72F, 0.21F, 0.4F, 0.93F, 0.07F, 0.11F, 0.41F, 0.24F, 0.04F,
110 0.36F, 0.15F, 0.92F, 0.08F, 0.99F, 0.35F, 0.42F, 0.7F, 0.17F, 0.39F, 0.69F, 0.0F, 0.86F,
111 0.89F, 0.59F, 0.81F, 0.77F, 0.15F, 0.89F, 0.17F, 0.76F, 0.67F, 0.58F, 0.78F, 0.26F, 0.19F,
112 0.69F, 0.18F, 0.46F, 0.6F, 0.69F, 0.23F, 0.34F, 0.3F, 0.64F, 0.34F, 0.6F, 0.99F, 0.69F,
113 0.57F, 0.75F, 0.07F, 0.36F, 0.75F, 0.81F, 0.8F, 0.42F, 0.09F, 0.94F, 0.66F, 0.35F, 0.67F,
114 0.34F, 0.66F, 0.02F, 0.47F, 0.78F, 0.21F, 0.02F, 0.18F, 0.42F, 0.2F, 0.46F, 0.34F, 0.4F,
115 0.46F, 0.96F, 0.86F, 0.25F, 0.25F, 0.22F, 0.37F, 0.59F, 0.19F, 0.45F, 0.61F, 0.04F, 0.71F,
116 0.77F, 0.51F, 0.77F, 0.15F, 0.78F, 0.36F, 0.62F, 0.24F, 0.86F, 0.2F, 0.77F, 0.08F, 0.09F,
117 0.3F, 0.0F, 0.6F, 0.99F, 0.69F,
120 //! A lot of bogus input spline derivatives - should have at list (max PME order = 5) * (DIM = 3) * (total unique atom number in all test cases = 16) values
121 std::vector<real> const c_sampleSplineDerivativesFull{
122 0.82F, 0.88F, 0.83F, 0.11F, 0.93F, 0.32F, 0.71F, 0.37F, 0.69F, 0.88F, 0.11F, 0.38F, 0.25F,
123 0.5F, 0.36F, 0.81F, 0.78F, 0.31F, 0.66F, 0.32F, 0.27F, 0.35F, 0.53F, 0.83F, 0.08F, 0.08F,
124 0.94F, 0.71F, 0.65F, 0.24F, 0.13F, 0.01F, 0.33F, 0.65F, 0.24F, 0.53F, 0.45F, 0.84F, 0.33F,
125 0.97F, 0.31F, 0.7F, 0.03F, 0.31F, 0.41F, 0.76F, 0.12F, 0.3F, 0.57F, 0.65F, 0.87F, 0.99F,
126 0.42F, 0.97F, 0.32F, 0.39F, 0.73F, 0.23F, 0.03F, 0.67F, 0.97F, 0.57F, 0.42F, 0.38F, 0.54F,
127 0.17F, 0.53F, 0.54F, 0.18F, 0.8F, 0.76F, 0.13F, 0.29F, 0.83F, 0.77F, 0.56F, 0.4F, 0.87F,
128 0.36F, 0.18F, 0.59F, 0.04F, 0.05F, 0.61F, 0.26F, 0.91F, 0.62F, 0.16F, 0.89F, 0.23F, 0.26F,
129 0.59F, 0.33F, 0.2F, 0.49F, 0.41F, 0.25F, 0.4F, 0.16F, 0.83F, 0.44F, 0.82F, 0.21F, 0.95F,
130 0.14F, 0.8F, 0.37F, 0.31F, 0.41F, 0.53F, 0.15F, 0.85F, 0.78F, 0.17F, 0.92F, 0.03F, 0.13F,
131 0.2F, 0.03F, 0.33F, 0.87F, 0.38F, 0, 0.08F, 0.79F, 0.36F, 0.53F, 0.05F, 0.07F, 0.94F,
132 0.23F, 0.85F, 0.13F, 0.27F, 0.23F, 0.22F, 0.26F, 0.38F, 0.15F, 0.48F, 0.18F, 0.33F, 0.23F,
133 0.62F, 0.1F, 0.36F, 0.99F, 0.07F, 0.02F, 0.04F, 0.09F, 0.29F, 0.52F, 0.29F, 0.83F, 0.97F,
134 0.61F, 0.81F, 0.49F, 0.56F, 0.08F, 0.09F, 0.03F, 0.65F, 0.46F, 0.1F, 0.06F, 0.06F, 0.39F,
135 0.29F, 0.04F, 0.03F, 0.1F, 0.83F, 0.94F, 0.59F, 0.97F, 0.82F, 0.2F, 0.66F, 0.23F, 0.11F,
136 0.03F, 0.16F, 0.27F, 0.53F, 0.94F, 0.46F, 0.43F, 0.29F, 0.97F, 0.64F, 0.46F, 0.37F, 0.43F,
137 0.48F, 0.37F, 0.93F, 0.5F, 0.2F, 0.92F, 0.09F, 0.74F, 0.55F, 0.44F, 0.05F, 0.13F, 0.17F,
138 0.79F, 0.44F, 0.11F, 0.6F, 0.64F, 0.05F, 0.96F, 0.3F, 0.45F, 0.47F, 0.42F, 0.74F, 0.91F,
139 0.06F, 0.89F, 0.24F, 0.26F, 0.68F, 0.4F, 0.88F, 0.5F, 0.65F, 0.48F, 0.15F, 0.0F, 0.41F,
140 0.67F, 0.4F, 0.31F, 0.73F, 0.77F, 0.36F, 0.26F, 0.74F, 0.46F, 0.56F, 0.78F, 0.92F, 0.32F,
141 0.9F, 0.06F, 0.55F, 0.6F, 0.13F, 0.38F, 0.93F, 0.5F, 0.92F, 0.96F, 0.82F, 0.0F, 0.04F,
142 0.9F, 0.55F, 0.97F, 1.0F, 0.23F, 0.46F, 0.52F, 0.49F, 0.0F, 0.32F, 0.16F, 0.4F, 0.62F,
143 0.36F, 0.03F, 0.63F, 0.16F, 0.58F, 0.97F, 0.03F, 0.44F, 0.07F, 0.22F, 0.75F, 0.32F, 0.61F,
144 0.94F, 0.33F, 0.7F, 0.57F, 0.5F, 0.84F, 0.7F, 0.47F, 0.18F, 0.09F, 0.25F, 0.77F, 0.94F,
145 0.85F, 0.09F, 0.83F, 0.02F, 0.91F,
148 //! 2 c_sample grids - only non-zero values have to be listed
149 std::vector<SparseRealGridValuesInput> const c_sampleGrids{ SparseRealGridValuesInput{
150 { IVec{ 0, 0, 0 }, 3.5F },
151 { IVec{ 7, 0, 0 }, -2.5F },
152 { IVec{ 3, 5, 7 }, -0.006F },
153 { IVec{ 1, 6, 7 }, -2.76F },
154 { IVec{ 3, 1, 2 }, 0.6F },
155 { IVec{ 6, 2, 4 }, 7.1F },
156 { IVec{ 4, 5, 6 }, 4.1F },
157 { IVec{ 4, 4, 6 }, -3.7F },
159 SparseRealGridValuesInput{
160 { IVec{ 0, 4, 0 }, 6.F },
161 { IVec{ 4, 2, 7 }, 13.76F },
162 { IVec{ 0, 6, 7 }, 3.6F },
163 { IVec{ 3, 2, 8 }, 0.61F },
164 { IVec{ 5, 4, 3 }, 2.1F },
165 { IVec{ 2, 5, 10 }, 3.6F },
166 { IVec{ 5, 3, 6 }, 2.1F },
167 { IVec{ 6, 6, 6 }, 6.1F },
170 //! Input forces for reduction
171 std::vector<RVec> const c_sampleForcesFull{
172 RVec{ 0.02F, 0.87F, 0.95F }, RVec{ 0.66F, 0.67F, 0.38F }, RVec{ 0.45F, 0.04F, 0.94F },
173 RVec{ 0.54F, 0.76F, 0.58F }, RVec{ 0.83F, 0.31F, 0.73F }, RVec{ 0.71F, 0.06F, 0.35F },
174 RVec{ 0.32F, 0.35F, 0.61F }, RVec{ 0.27F, 0.98F, 0.83F }, RVec{ 0.11F, 0.3F, 0.42F },
175 RVec{ 0.95F, 0.69F, 0.58F }, RVec{ 0.29F, 0.1F, 0.68F }, RVec{ 0.94F, 0.62F, 0.51F },
176 RVec{ 0.47F, 0.04F, 0.47F }, RVec{ 0.34F, 0.71F, 0.52F }
179 //! PME orders to test
180 std::vector<int> const pmeOrders{ 3, 4, 5 };
181 //! Atom counts to test
182 std::vector<size_t> const atomCounts{ 1, 2, 13 };
184 /* Helper structures for test input */
186 //! A structure for all the spline data which depends in size both on the PME order and atom count
187 struct AtomAndPmeOrderSizedData
190 SplineParamsVector splineValues;
191 //! Spline derivatives
192 SplineParamsVector splineDerivatives;
195 //! A structure for all the input atom data which depends in size on atom count - including a range of spline data for different PME orders
199 GridLineIndicesVector gridLineIndices;
201 ChargesVector charges;
203 CoordinatesVector coordinates;
204 //! Spline data for different orders
205 std::map<int, AtomAndPmeOrderSizedData> splineDataByPmeOrder;
208 //! A range of the test input data sets, uniquely identified by the atom count
209 typedef std::map<size_t, AtomSizedData> InputDataByAtomCount;
211 /*! \brief Convenience typedef of the test input parameters - unit cell box, PME interpolation
212 * order, grid dimensions, grid values, overwriting/reducing the input forces, atom count.
214 * The rest of the atom-related input data - gridline indices, spline theta values, spline dtheta
215 * values, atom charges - is looked up in the inputAtomDataSets_ test fixture variable.
217 typedef std::tuple<Matrix3x3, int, IVec, SparseRealGridValuesInput, size_t> GatherInputParameters;
220 class PmeGatherTest : public ::testing::TestWithParam<GatherInputParameters>
223 //! Storage of all the input atom datasets
224 static InputDataByAtomCount s_inputAtomDataSets_;
227 PmeGatherTest() = default;
228 //! Sets the input atom data references once
229 static void SetUpTestCase()
232 for (auto atomCount : atomCounts)
234 AtomSizedData atomData;
235 atomData.gridLineIndices =
236 GridLineIndicesVector(c_sampleGridLineIndicesFull).subArray(start, atomCount);
237 atomData.charges = ChargesVector(c_sampleChargesFull).subArray(start, atomCount);
239 atomData.coordinates.resize(atomCount, RVec{ 1e6, 1e7, -1e8 });
240 /* The coordinates are intentionally bogus in this test - only the size matters; the gridline indices are fed directly as inputs */
241 for (auto pmeOrder : pmeOrders)
243 AtomAndPmeOrderSizedData splineData;
244 const size_t dimSize = atomCount * pmeOrder;
245 for (int dimIndex = 0; dimIndex < DIM; dimIndex++)
247 splineData.splineValues[dimIndex] =
248 SplineParamsDimVector(c_sampleSplineValuesFull).subArray(dimIndex * dimSize, dimSize);
249 splineData.splineDerivatives[dimIndex] =
250 SplineParamsDimVector(c_sampleSplineDerivativesFull).subArray(dimIndex * dimSize, dimSize);
252 atomData.splineDataByPmeOrder[pmeOrder] = splineData;
254 s_inputAtomDataSets_[atomCount] = atomData;
260 /* Getting the input */
265 SparseRealGridValuesInput nonZeroGridValues;
266 std::tie(box, pmeOrder, gridSize, nonZeroGridValues, atomCount) = GetParam();
267 auto inputAtomData = s_inputAtomDataSets_[atomCount];
268 auto inputAtomSplineData = inputAtomData.splineDataByPmeOrder[pmeOrder];
270 /* Storing the input where it's needed, running the test */
272 inputRec.nkx = gridSize[XX];
273 inputRec.nky = gridSize[YY];
274 inputRec.nkz = gridSize[ZZ];
275 inputRec.pme_order = pmeOrder;
276 inputRec.coulombtype = eelPME;
277 inputRec.epsilon_r = 1.0;
279 TestReferenceData refData;
280 for (const auto& context : getPmeTestEnv()->getHardwareContexts())
282 CodePath codePath = context->getCodePath();
283 const bool supportedInput =
284 pmeSupportsInputForMode(*getPmeTestEnv()->hwinfo(), &inputRec, codePath);
287 /* Testing the failure for the unsupported input */
288 EXPECT_THROW_GMX(pmeInitWrapper(&inputRec, codePath, nullptr, nullptr, box),
289 NotImplementedError);
293 /* Describing the test uniquely */
295 formatString("Testing force gathering with %s %sfor PME grid size %d %d %d"
296 ", order %d, %zu atoms",
297 codePathToString(codePath), context->getDescription().c_str(),
298 gridSize[XX], gridSize[YY], gridSize[ZZ], pmeOrder, atomCount));
300 PmeSafePointer pmeSafe = pmeInitWrapper(&inputRec, codePath, context->getDeviceInfo(),
301 context->getPmeGpuProgram(), box);
302 std::unique_ptr<StatePropagatorDataGpu> stateGpu =
303 (codePath == CodePath::GPU)
304 ? makeStatePropagatorDataGpu(*pmeSafe.get(), context->deviceContext())
307 pmeInitAtoms(pmeSafe.get(), stateGpu.get(), codePath, inputAtomData.coordinates,
308 inputAtomData.charges);
310 /* Setting some more inputs */
311 pmeSetRealGrid(pmeSafe.get(), codePath, nonZeroGridValues);
312 pmeSetGridLineIndices(pmeSafe.get(), codePath, inputAtomData.gridLineIndices);
313 for (int dimIndex = 0; dimIndex < DIM; dimIndex++)
315 pmeSetSplineData(pmeSafe.get(), codePath, inputAtomSplineData.splineValues[dimIndex],
316 PmeSplineDataType::Values, dimIndex);
317 pmeSetSplineData(pmeSafe.get(), codePath, inputAtomSplineData.splineDerivatives[dimIndex],
318 PmeSplineDataType::Derivatives, dimIndex);
321 /* Explicitly copying the sample forces to be able to modify them */
322 auto inputForcesFull(c_sampleForcesFull);
323 GMX_RELEASE_ASSERT(inputForcesFull.size() >= atomCount, "Bad input forces size");
324 auto forces = ForcesVector(inputForcesFull).subArray(0, atomCount);
326 /* Running the force gathering itself */
327 pmePerformGather(pmeSafe.get(), codePath, forces);
328 pmeFinalizeTest(pmeSafe.get(), codePath);
330 /* Check the output forces correctness */
331 TestReferenceChecker forceChecker(refData.rootChecker());
332 const auto ulpTolerance = 3 * pmeOrder;
333 forceChecker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpTolerance));
334 forceChecker.checkSequence(forces.begin(), forces.end(), "Forces");
339 // An instance of static atom data
340 InputDataByAtomCount PmeGatherTest::s_inputAtomDataSets_;
342 //! Test for PME force gathering
343 TEST_P(PmeGatherTest, ReproducesOutputs)
345 EXPECT_NO_THROW_GMX(runTest());
348 //! Instantiation of the PME gathering test
349 INSTANTIATE_TEST_CASE_P(SaneInput,
351 ::testing::Combine(::testing::ValuesIn(c_sampleBoxes),
352 ::testing::ValuesIn(pmeOrders),
353 ::testing::ValuesIn(c_sampleGridSizes),
354 ::testing::ValuesIn(c_sampleGrids),
355 ::testing::ValuesIn(atomCounts)));