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
147 updatedPositions_(std::begin(g_positions), std::end(g_positions)),
148 velocities_(updatedPositions_.size(), 0)
150 set_pbc(&pbcNone_, epbcNONE, g_box);
151 set_pbc(&pbcXYZ_, epbcXYZ, g_box);
153 // Perturb the atom positions, to appear like an
154 // "update," and where there is definitely constraining
156 for (size_t i = 0; i != updatedPositions_.size(); ++i)
160 updatedPositions_[i] += 0.01;
164 updatedPositions_[i] -= 0.01;
168 updatedPositions_[i] += 0.02;
172 updatedPositions_[i] -= 0.02;
178 TEST_P(SettleTest, SatisfiesConstraints)
181 bool usePbc, useVelocities, calcVirial;
182 // Make some symbolic names for the parameter combination under
184 std::tie(numSettles, usePbc, useVelocities, calcVirial) = GetParam();
186 // Make a string that describes which parameter combination is
187 // being tested, to help make failing tests comprehensible.
188 std::string testDescription = formatString("while testing %d SETTLEs, %sPBC, %svelocities and %scalculating the virial",
190 usePbc ? "with " : "without ",
191 useVelocities ? "with " : "without ",
192 calcVirial ? "" : "not ");
194 const int settleType = 0;
195 const int atomsPerSettle = NRAL(F_SETTLE);
196 ASSERT_LE(numSettles, updatedPositions_.size() / (atomsPerSettle * DIM)) << "cannot test that many SETTLEs " << testDescription;
198 // Set up the topology.
200 mtop.moltype.resize(1);
201 mtop.molblock.resize(1);
202 mtop.molblock[0].type = 0;
203 std::vector<int> &iatoms = mtop.moltype[0].ilist[F_SETTLE].iatoms;
204 for (int i = 0; i < numSettles; ++i)
206 iatoms.push_back(settleType);
207 iatoms.push_back(i*atomsPerSettle + 0);
208 iatoms.push_back(i*atomsPerSettle + 1);
209 iatoms.push_back(i*atomsPerSettle + 2);
212 // Set up the SETTLE parameters.
213 mtop.ffparams.ntypes = 1;
214 snew(mtop.ffparams.iparams, mtop.ffparams.ntypes);
215 const real dOH = 0.09572;
216 const real dHH = 0.15139;
217 mtop.ffparams.iparams[settleType].settle.doh = dOH;
218 mtop.ffparams.iparams[settleType].settle.dhh = dHH;
220 // Set up the masses.
222 std::vector<real> mass, massReciprocal;
223 const real oxygenMass = 15.9994, hydrogenMass = 1.008;
224 for (int i = 0; i < numSettles; ++i)
226 mass.push_back(oxygenMass);
227 mass.push_back(hydrogenMass);
228 mass.push_back(hydrogenMass);
229 massReciprocal.push_back(1./oxygenMass);
230 massReciprocal.push_back(1./hydrogenMass);
231 massReciprocal.push_back(1./hydrogenMass);
233 mdatoms.massT = mass.data();
234 mdatoms.invmass = massReciprocal.data();
235 mdatoms.homenr = numSettles * atomsPerSettle;
237 // Finally make the settle data structures
238 settledata *settled = settle_init(mtop);
239 const t_ilist ilist = { mtop.moltype[0].ilist[F_SETTLE].size(), 0, mtop.moltype[0].ilist[F_SETTLE].iatoms.data(), 0 };
240 settle_set_constraints(settled, &ilist, mdatoms);
242 // Copy the original positions from the array of doubles to a vector of reals
243 std::vector<real> startingPositions(std::begin(g_positions), std::end(g_positions));
247 tensor virial = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
248 int numThreads = 1, threadIndex = 0;
249 const real reciprocalTimeStep = 1.0/0.002;
250 csettle(settled, numThreads, threadIndex,
251 usePbc ? &pbcXYZ_ : &pbcNone_,
252 startingPositions.data(), updatedPositions_.data(), reciprocalTimeStep,
253 useVelocities ? velocities_.data() : nullptr,
254 calcVirial, virial, &errorOccured);
255 settle_free(settled);
256 EXPECT_FALSE(errorOccured) << testDescription;
258 // The necessary tolerances for the test to pass were determined
259 // empirically. This isn't nice, but the required behaviour that
260 // SETTLE produces constrained coordinates consistent with
261 // sensible sampling needs to be tested at a much higher level.
262 FloatingPointTolerance tolerance =
263 relativeToleranceAsPrecisionDependentUlp(dOH*dOH, 80, 380);
265 // Verify the updated coordinates match the requirements
266 int positionIndex = 0, velocityIndex = 0;
267 for (int i = 0; i < numSettles; ++i)
270 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
273 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
276 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
279 EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH1), tolerance) << formatString("for water %d ", i) << testDescription;
280 EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
281 EXPECT_REAL_EQ_TOL(dHH*dHH, distance2(positionH1, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
283 // This merely tests whether the velocities were
284 // updated from the starting values of zero (or not),
285 // but not whether the update was correct.
286 for (int j = 0; j < atomsPerSettle * DIM; ++j, ++velocityIndex)
288 EXPECT_TRUE(useVelocities == (0. != velocities_[velocityIndex])) << formatString("for water %d velocity coordinate %d ", i, j) << testDescription;
292 // This merely tests whether the viral was updated from
293 // the starting values of zero (or not), but not whether
294 // any update was correct.
295 for (int d = 0; d < DIM; ++d)
297 for (int dd = 0; dd < DIM; ++dd)
299 EXPECT_TRUE(calcVirial == (0. != virial[d][dd])) << formatString("for virial component[%d][%d] ", d, dd) << testDescription;
304 // Scan the full Cartesian product of numbers of SETTLE interactions
305 // (4 and 17 are chosen to test cases that do and do not match
306 // hardware SIMD widths), and whether or not we use PBC, velocities or
307 // calculate the virial contribution.
308 INSTANTIATE_TEST_CASE_P(WithParameters, SettleTest,
309 ::testing::Combine(::testing::Values(1, 4, 7),