2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2018, 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 #include "gromacs/mdlib/settle.h"
42 #include <gtest/gtest.h>
44 #include "gromacs/math/vec.h"
45 #include "gromacs/math/vectypes.h"
46 #include "gromacs/mdtypes/mdatom.h"
47 #include "gromacs/pbcutil/pbc.h"
48 #include "gromacs/topology/idef.h"
49 #include "gromacs/topology/ifunc.h"
50 #include "gromacs/topology/topology.h"
51 #include "gromacs/utility/smalloc.h"
52 #include "gromacs/utility/stringutil.h"
53 #include "gromacs/utility/unique_cptr.h"
55 #include "testutils/testasserts.h"
63 //! Database of 51 water atom input positions (DIM reals per atom, taken from spc216.gro) for use as test inputs.
64 const double g_positions[] = {
118 //! Simple cubic simulation box to use in tests
119 matrix g_box = {{real(1.86206), 0, 0}, {0, real(1.86206), 0}, {0, 0, real(1.86206)}};
121 //! Convenience typedef
122 typedef std::tuple<int, bool, bool, bool> SettleTestParameters;
124 /*! \brief Test fixture for testing SETTLE position updates
126 * \todo This also tests that if the calling code requires velocities
127 * and virial updates, that those outputs do change, but does not test
128 * that those changes are correct.
130 * \todo Only no-PBC and cubic-PBC are tested here, but the correct
131 * function of the SIMD version of set_pbx_auic in all cases should be
134 class SettleTest : public ::testing::TestWithParam<SettleTestParameters>
137 //! Updated water atom positions to constrain (DIM reals per atom)
138 std::vector<real> updatedPositions_;
139 //! Water atom velocities to constrain (DIM reals per atom)
140 std::vector<real> velocities_;
141 //! PBC option to test
143 //! PBC option to test
148 updatedPositions_(std::begin(g_positions), std::end(g_positions)),
149 velocities_(updatedPositions_.size(), 0)
151 set_pbc(&pbcNone_, epbcNONE, g_box);
152 set_pbc(&pbcXYZ_, epbcXYZ, g_box);
154 // Perturb the atom positions, to appear like an
155 // "update," and where there is definitely constraining
157 for (size_t i = 0; i != updatedPositions_.size(); ++i)
161 updatedPositions_[i] += 0.01;
165 updatedPositions_[i] -= 0.01;
169 updatedPositions_[i] += 0.02;
173 updatedPositions_[i] -= 0.02;
179 TEST_P(SettleTest, SatisfiesConstraints)
182 bool usePbc, useVelocities, calcVirial;
183 // Make some symbolic names for the parameter combination under
185 std::tie(numSettles, usePbc, useVelocities, calcVirial) = GetParam();
187 // Make a string that describes which parameter combination is
188 // being tested, to help make failing tests comprehensible.
189 std::string testDescription = formatString("while testing %d SETTLEs, %sPBC, %svelocities and %scalculating the virial",
191 usePbc ? "with " : "without ",
192 useVelocities ? "with " : "without ",
193 calcVirial ? "" : "not ");
195 const int settleType = 0;
196 const int atomsPerSettle = NRAL(F_SETTLE);
197 ASSERT_LE(numSettles, updatedPositions_.size() / (atomsPerSettle * DIM)) << "cannot test that many SETTLEs " << testDescription;
199 // Set up the topology.
201 mtop.moltype.resize(1);
202 mtop.molblock.resize(1);
203 mtop.molblock[0].type = 0;
204 const int settleStride = 1 + atomsPerSettle;
206 snew(iatoms, numSettles*settleStride);
207 for (int i = 0; i < numSettles; ++i)
209 iatoms[i*settleStride + 0] = settleType;
210 iatoms[i*settleStride + 1] = i*atomsPerSettle + 0;
211 iatoms[i*settleStride + 2] = i*atomsPerSettle + 1;
212 iatoms[i*settleStride + 3] = i*atomsPerSettle + 2;
214 mtop.moltype[0].ilist[F_SETTLE].iatoms = iatoms;
215 mtop.moltype[0].ilist[F_SETTLE].nr = numSettles*settleStride;
217 // Set up the SETTLE parameters.
218 mtop.ffparams.ntypes = 1;
219 snew(mtop.ffparams.iparams, mtop.ffparams.ntypes);
220 const real dOH = 0.09572;
221 const real dHH = 0.15139;
222 mtop.ffparams.iparams[settleType].settle.doh = dOH;
223 mtop.ffparams.iparams[settleType].settle.dhh = dHH;
225 // Set up the masses.
227 std::vector<real> mass, massReciprocal;
228 const real oxygenMass = 15.9994, hydrogenMass = 1.008;
229 for (int i = 0; i < numSettles; ++i)
231 mass.push_back(oxygenMass);
232 mass.push_back(hydrogenMass);
233 mass.push_back(hydrogenMass);
234 massReciprocal.push_back(1./oxygenMass);
235 massReciprocal.push_back(1./hydrogenMass);
236 massReciprocal.push_back(1./hydrogenMass);
238 mdatoms.massT = mass.data();
239 mdatoms.invmass = massReciprocal.data();
240 mdatoms.homenr = numSettles * atomsPerSettle;
242 // Finally make the settle data structures
243 settledata *settled = settle_init(mtop);
244 settle_set_constraints(settled, &mtop.moltype[0].ilist[F_SETTLE], mdatoms);
246 // Copy the original positions from the array of doubles to a vector of reals
247 std::vector<real> startingPositions(std::begin(g_positions), std::end(g_positions));
251 tensor virial = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
252 int numThreads = 1, threadIndex = 0;
253 const real reciprocalTimeStep = 1.0/0.002;
254 csettle(settled, numThreads, threadIndex,
255 usePbc ? &pbcXYZ_ : &pbcNone_,
256 startingPositions.data(), updatedPositions_.data(), reciprocalTimeStep,
257 useVelocities ? velocities_.data() : nullptr,
258 calcVirial, virial, &errorOccured);
259 settle_free(settled);
260 EXPECT_FALSE(errorOccured) << testDescription;
262 // The necessary tolerances for the test to pass were determined
263 // empirically. This isn't nice, but the required behaviour that
264 // SETTLE produces constrained coordinates consistent with
265 // sensible sampling needs to be tested at a much higher level.
266 FloatingPointTolerance tolerance =
267 relativeToleranceAsPrecisionDependentUlp(dOH*dOH, 80, 380);
269 // Verify the updated coordinates match the requirements
270 int positionIndex = 0, velocityIndex = 0;
271 for (int i = 0; i < numSettles; ++i)
274 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
277 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
280 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
283 EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH1), tolerance) << formatString("for water %d ", i) << testDescription;
284 EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
285 EXPECT_REAL_EQ_TOL(dHH*dHH, distance2(positionH1, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
287 // This merely tests whether the velocities were
288 // updated from the starting values of zero (or not),
289 // but not whether the update was correct.
290 for (int j = 0; j < atomsPerSettle * DIM; ++j, ++velocityIndex)
292 EXPECT_TRUE(useVelocities == (0. != velocities_[velocityIndex])) << formatString("for water %d velocity coordinate %d ", i, j) << testDescription;
296 // This merely tests whether the viral was updated from
297 // the starting values of zero (or not), but not whether
298 // any update was correct.
299 for (int d = 0; d < DIM; ++d)
301 for (int dd = 0; dd < DIM; ++dd)
303 EXPECT_TRUE(calcVirial == (0. != virial[d][dd])) << formatString("for virial component[%d][%d] ", d, dd) << testDescription;
308 // Scan the full Cartesian product of numbers of SETTLE interactions
309 // (4 and 17 are chosen to test cases that do and do not match
310 // hardware SIMD widths), and whether or not we use PBC, velocities or
311 // calculate the virial contribution.
312 INSTANTIATE_TEST_CASE_P(WithParameters, SettleTest,
313 ::testing::Combine(::testing::Values(1, 4, 7),