2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2021, 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 * \brief Sanity checks for virtual sites
39 * The tests in this file test the virtual site implementation in two ways.
40 * 1) An artificial test system containing all virtual site types is run
41 * end-to-end. The virtual sites are recalculated using a reference
42 * implementation, and compared to the trajectory values. This ensures
43 * that no relevant (real or virtual) coordinates were changed between
44 * virtual site computation, and that the mdrun implementation agrees
45 * with the reference implementation. The latter has the advantage to
46 * be written closer to the analytical expressions, hence easier to
47 * check by eye. Unlike the mdrun implementation, it can also be unit-
48 * tested (see 2)). Since this is an end-to-end test, it also ensures
49 * that virtual sites can be processed by grompp and run by mdrun.
50 * 2) The reference implementation is tested to have corresponding positions
51 * and velocities. This is achieved by comparing virtual site positions
52 * calculated from propagated real positions to virtual site positions
53 * calculated by propagating virtual sites using the virtual velocities.
55 * Note, this only ensures that the position and velocity implementation
56 * match, not that they are actually correct. Some regression test systems
57 * include virtual sites, so there is some testing that no bugs are introduced.
58 * It would be good to have unit tests, though. This can either be achieved by
59 * refactoring the mdrun implementation, or adding them to the reference
60 * implementation here. See also #3911.
62 * \author Pascal Merz <pascal.merz@me.com>
63 * \ingroup module_mdrun_integration_tests
69 #include "gromacs/topology/idef.h"
70 #include "gromacs/topology/ifunc.h"
71 #include "gromacs/trajectoryanalysis/topologyinformation.h"
72 #include "gromacs/utility/stringutil.h"
74 #include "testutils/mpitest.h"
75 #include "testutils/simulationdatabase.h"
76 #include "testutils/testmatchers.h"
78 #include "gromacs/utility/strconvert.h"
80 #include "moduletest.h"
81 #include "simulatorcomparison.h"
82 #include "trajectoryreader.h"
88 using VirtualSiteTestParams = std::tuple<std::string, std::string, std::string>;
89 class VirtualSiteTest : public MdrunTestFixture, public ::testing::WithParamInterface<VirtualSiteTestParams>
94 /*! \brief Check against reference implementation
96 * Check that the positions and velocities of virtual sites are equal to the reference
97 * implementation, within a given tolerance. This loops over the trajectory, calculates
98 * the virtual site positions and velocities using a reference implementation, and
99 * compares it with the reported virtual site positions and velocities.
101 * The expectation is that the trajectory and the reference calculations are identical.
102 * This ensures that neither the real atoms nor the virtual sites are changed between
103 * virtual site calculation and trajectory printing. The reference implementation is
104 * also closer to the analytically derived equations than the mdrun implementation,
105 * making it easier to verify. Finally, the reference implementation is tested (within
106 * this file), which is a reasonable work-around for the fact that the actual implementation
107 * doesn't have unit tests.
109 * Note that the reference implementation does not take into account PBC. It's intended
110 * to be used with simple test cases in which molecules are ensured not to be broken
111 * across periodic boundaries.
113 static void checkVirtualSitesAgainstReferenceImplementation(const std::string& trajectoryName,
114 ArrayRef<const VirtualSite> virtualSites,
115 FloatingPointTolerance tolerance)
118 "Checking virtual site positions and velocities against reference implementation.");
119 TrajectoryFrameReader trajectoryFrameReader(trajectoryName);
120 while (trajectoryFrameReader.readNextFrame())
122 const auto frame = trajectoryFrameReader.frame();
123 SCOPED_TRACE(formatString("Checking frame %s", frame.frameName().c_str()));
124 std::vector<RVec> positions;
125 std::vector<RVec> velocities;
126 std::vector<RVec> refPositions;
127 std::vector<RVec> refVelocities;
128 for (const auto& vSite : virtualSites)
130 auto [refPosition, refVelocity] = vSite.calculate(frame.x(), frame.v());
131 refPositions.emplace_back(refPosition);
132 refVelocities.emplace_back(refVelocity);
133 positions.emplace_back(frame.x().at(vSite.atomIdx));
134 velocities.emplace_back(frame.v().at(vSite.atomIdx));
136 EXPECT_THAT(refPositions, Pointwise(RVecEq(tolerance), positions));
137 EXPECT_THAT(refVelocities, Pointwise(RVecEq(tolerance), velocities));
141 /*! \brief Check the reference implementation
143 * This tests that the reference implementation position and velocities correspond.
145 * a) generating real atom starting positions (x(t)) and half-step velocities (v(t+dt/2))
146 * b) propagating the real atoms positions by one time step x(t+dt) = x(t) + dt*v(t+dt/2)
147 * c) calculating the half-step positions x(t+dt/2) = 0.5 * (x(t) + x(t+dt))
148 * d) calculating the virtual position xv(t) := xv(x(t)) and xv1(t+dt) := xv(x(t+dt))
149 * using the reference implementation
150 * e) calculating the virtual velocity vv(t+dt/2) := vv(x(t+dt/2), v(t+dt/2))
151 * using the reference implementation
152 * f) calculating the virtual position xv2(t+dt) = xv(t) + dt*xv(t+dt/2)
153 * g) comparing xv1(t+dt) and xv2(t+dt)
154 * If the calculation of the virtual positions and velocities correspond, xv1 and xv2 will
155 * be identical up to some integration error.
157 * Maybe unused because this test runs only in double precision.
159 [[maybe_unused]] static void checkReferenceImplementation()
161 SCOPED_TRACE("Checking virtual site reference implementation.");
162 // Randomly generated real atom positions and velocities
163 std::vector<RVec> startPositions = { { 2.641321, 2.076298, 2.138602 },
164 { 3.776765, 3.154901, 1.556379 },
165 { 2.376669, 1.166706, 2.457044 },
166 { 3.242320, 2.142465, 2.023578 } };
167 std::vector<RVec> halfStepVelocities = { { 0.154667, 0.319010, 0.458749 },
168 { -0.010590, -0.191858, -0.096820 },
169 { -0.008609, 0.004656, 0.448852 },
170 { 0.411874, -0.038205, -0.151459 } };
171 // Virtual site definitions with randomly generated parameters
172 std::vector<VirtualSite> virtualSites = {
173 { F_VSITE1, 6, { 0 }, {} },
174 { F_VSITE2, 6, { 0, 1 }, { 0.710573 } },
175 { F_VSITE2FD, 6, { 0, 1 }, { 0.292430 } },
176 { F_VSITE3, 6, { 0, 1, 2 }, { 0.060990, 0.543636 } },
177 { F_VSITE3FD, 6, { 0, 1, 2 }, { 0.125024, 0.444587 } },
178 { F_VSITE3FAD, 6, { 0, 1, 2 }, { 0.414850, 0.349767 } },
179 { F_VSITE3OUT, 6, { 0, 1, 2 }, { 0.779323, 0.093773, 0.743164 } },
180 { F_VSITE4FDN, 6, { 0, 1, 2, 3 }, { 0.975111, 0.952180, 0.757594 } }
183 // Make integration step
184 const real timeStep = 1e-5;
185 std::vector<RVec> endPositions;
186 std::vector<RVec> halfStepPositions;
187 GMX_RELEASE_ASSERT(startPositions.size() == halfStepVelocities.size(),
188 "Need positions and velocities for every real atom.");
189 const auto numRealAtoms = startPositions.size();
190 for (auto idx = decltype(numRealAtoms){ 0 }; idx < numRealAtoms; idx++)
192 endPositions.emplace_back(startPositions[idx] + timeStep * halfStepVelocities[idx]);
193 halfStepPositions.emplace_back(startPositions[idx]
194 + real(0.5) * timeStep * halfStepVelocities[idx]);
197 // Check that displacement equals the calculated velocities
198 for (const auto& vSite : virtualSites)
200 SCOPED_TRACE(formatString("Checking %s", interaction_function[vSite.type].longname));
202 // GCC 7 falsely flags unused "variables" in structured bindings, GCC 8+ fixed this
203 // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=81767
205 GCC_DIAGNOSTIC_IGNORE(-Wunused-variable)
207 /* Calculate start and end virtual position
209 * The reference implementation always calculates the virtual velocity, but
210 * since we don't have real velocities at full steps, these virtual velocities
211 * are unprecise. We don't need them anyway, so we'll just ignore them.
213 const auto [startVPosition, vVelocityUnused1] =
214 vSite.calculate(startPositions, halfStepVelocities);
215 const auto [endVPosition1, vVelocityUnused2] =
216 vSite.calculate(endPositions, halfStepVelocities);
218 /* Calculate virtual velocity at half step using reference implementation
220 * The virtual positions are exact, but we don't need them.
222 const auto [halfStepVPositionUnused, halfStepVVelocity] =
223 vSite.calculate(halfStepPositions, halfStepVelocities);
226 // We can now integrate the virtual positions using the half step velocity
227 const auto endVPosition2 = startVPosition + timeStep * halfStepVVelocity;
229 /* We can now calculate the displacement of the virtual site in two ways:
230 * (1) Using endVPosition1, calculated by advancing the real positions
231 * and calculating the virtual position from them.
232 * (2) Using endVPosition2, calculated by advancing the virtual position
233 * by the virtual velocity.
234 * Comparing the difference of the displacement with relative tolerance makes
235 * this at least somewhat independent of the choice of time step and coordinates -
236 * assuming that the displacement doesn't go to zero, of course!
238 const auto displacement1 = endVPosition1 - startVPosition;
239 const auto displacement2 = endVPosition2 - startVPosition;
240 ASSERT_GT(std::abs(displacement1[XX]), 0)
241 << "adjust choice of coordinates or time step.";
242 ASSERT_GT(std::abs(displacement1[YY]), 0)
243 << "adjust choice of coordinates or time step.";
244 ASSERT_GT(std::abs(displacement1[ZZ]), 0)
245 << "adjust choice of coordinates or time step.";
247 EXPECT_REAL_EQ_TOL(displacement1[XX],
249 relativeToleranceAsFloatingPoint(displacement1[XX], 1e-9));
250 EXPECT_REAL_EQ_TOL(displacement1[YY],
252 relativeToleranceAsFloatingPoint(displacement1[YY], 1e-9));
253 EXPECT_REAL_EQ_TOL(displacement1[ZZ],
255 relativeToleranceAsFloatingPoint(displacement1[ZZ], 1e-9));
259 //! Holds parameters of virtual site and allows calculation
262 //! Type of virtual site
264 //! Index of virtual site
266 //! Indices of constructing atoms
267 std::vector<int> constructingAtomIdx;
268 //! Construction parameters
269 std::vector<real> parameters;
271 //! Dispatch function to compute position and velocity of virtual site from reference implementation based on \p type
272 [[nodiscard]] std::tuple<RVec, RVec> calculate(ArrayRef<const RVec> constructingPositions,
273 ArrayRef<const RVec> constructingVelocities) const;
276 //! Templated reference implementation of virtual site position and velocity calculation
277 template<int vsiteType>
278 [[nodiscard]] std::tuple<RVec, RVec> calculateVSite(ArrayRef<const RVec> positions,
279 ArrayRef<const RVec> velocities) const;
282 /*! \brief Helper function returning a list of virtual sites from the topology
284 * This also prints the indices of the virtual sites. If any tests fail, this
285 * can be used to understand which type is failing.
287 static std::vector<VirtualSite> vSiteList(const TopologyInformation& topologyInformation)
289 std::vector<VirtualSite> virtualSites;
290 const auto& localTopology = *topologyInformation.expandedTopology();
291 printf("Reading virtual site types...\n");
292 for (int vsiteType = F_VSITE1; vsiteType <= F_VSITEN; vsiteType++)
294 const auto& interactionList = localTopology.idef.il.at(vsiteType);
295 if (vsiteType == F_VSITE4FD || interactionList.empty())
297 // 4FD is deprecated. Interaction list empty means system doesn't contain this type.
300 const int numConstructingAtoms = interaction_function[vsiteType].nratoms - 1;
301 const int defaultIncrement = numConstructingAtoms + 2;
302 std::string indexString;
304 for (int i = 0; i < interactionList.size();)
306 const int parameterIdx = interactionList.iatoms[i];
307 const int virtualSiteIdx = interactionList.iatoms[i + 1];
308 if (!indexString.empty())
312 indexString += toString(virtualSiteIdx);
314 if (vsiteType == F_VSITEN)
316 const int vSiteNConstructingAtoms =
317 localTopology.idef.iparams[parameterIdx].vsiten.n;
318 VirtualSite vSite{ vsiteType, virtualSiteIdx, { interactionList.iatoms[i + 2] }, {} };
319 for (int j = 3; j < 3 * vSiteNConstructingAtoms; j += 3)
321 vSite.constructingAtomIdx.push_back(interactionList.iatoms[j + 2]);
322 vSite.parameters.push_back(
323 localTopology.idef.iparams[interactionList.iatoms[j]].vsiten.a);
325 virtualSites.push_back(std::move(vSite));
326 i += 3 * vSiteNConstructingAtoms;
330 virtualSites.emplace_back(VirtualSite{
333 { &interactionList.iatoms[i + 2],
334 &interactionList.iatoms[i + 2 + numConstructingAtoms] },
335 { localTopology.idef.iparams[parameterIdx].generic.buf,
336 localTopology.idef.iparams[parameterIdx].generic.buf + MAXFORCEPARAM } });
337 i += defaultIncrement;
345 // check-source gets confused by these
348 [[nodiscard]] std::tuple<RVec, RVec>
349 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE1>(ArrayRef<const RVec> positions,
350 ArrayRef<const RVec> velocities) const
352 return { positions[constructingAtomIdx.at(0)], velocities[constructingAtomIdx.at(0)] };
355 [[nodiscard]] std::tuple<RVec, RVec>
356 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE2>(ArrayRef<const RVec> positions,
357 ArrayRef<const RVec> velocities) const
359 const auto& a = parameters[0];
360 return { (1 - a) * positions[constructingAtomIdx.at(0)] + a * positions[constructingAtomIdx.at(1)],
361 (1 - a) * velocities[constructingAtomIdx.at(0)] + a * velocities[constructingAtomIdx.at(1)] };
364 [[nodiscard]] std::tuple<RVec, RVec>
365 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE2FD>(ArrayRef<const RVec> positions,
366 ArrayRef<const RVec> velocities) const
368 const auto& a = parameters[0];
369 const auto& ri = positions[constructingAtomIdx.at(0)];
370 const auto& rj = positions[constructingAtomIdx.at(1)];
371 const auto& vi = velocities[constructingAtomIdx.at(0)];
372 const auto& vj = velocities[constructingAtomIdx.at(1)];
373 const auto rij = rj - ri;
374 const auto vij = vj - vi;
376 return { ri + (a / rij.norm()) * rij,
377 vi + (a / rij.norm()) * (vij - rij * (vij.dot(rij) / rij.norm2())) };
380 [[nodiscard]] std::tuple<RVec, RVec>
381 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE3>(ArrayRef<const RVec> positions,
382 ArrayRef<const RVec> velocities) const
384 const auto& a = parameters[0];
385 const auto& b = parameters[1];
386 const auto& ri = positions[constructingAtomIdx.at(0)];
387 const auto& rj = positions[constructingAtomIdx.at(1)];
388 const auto& rk = positions[constructingAtomIdx.at(2)];
389 const auto& vi = velocities[constructingAtomIdx.at(0)];
390 const auto& vj = velocities[constructingAtomIdx.at(1)];
391 const auto& vk = velocities[constructingAtomIdx.at(2)];
393 return { (1 - a - b) * ri + a * rj + b * rk, (1 - a - b) * vi + a * vj + b * vk };
396 [[nodiscard]] std::tuple<RVec, RVec>
397 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE3FD>(ArrayRef<const RVec> positions,
398 ArrayRef<const RVec> velocities) const
400 const auto& a = parameters[0];
401 const auto& b = parameters[1];
402 const auto& ri = positions[constructingAtomIdx.at(0)];
403 const auto& rj = positions[constructingAtomIdx.at(1)];
404 const auto& rk = positions[constructingAtomIdx.at(2)];
405 const auto& vi = velocities[constructingAtomIdx.at(0)];
406 const auto& vj = velocities[constructingAtomIdx.at(1)];
407 const auto& vk = velocities[constructingAtomIdx.at(2)];
408 const auto rij = rj - ri;
409 const auto rjk = rk - rj;
410 const auto vij = vj - vi;
411 const auto vjk = vk - vj;
413 // TODO: Should be uncommented after resolution of #3909
414 const auto rijk = /*(1 - a) **/ rij + a * rjk;
415 const auto vijk = /*(1 - a) **/ vij + a * vjk;
417 return { ri + (b / rijk.norm()) * rijk,
418 vi + (b / rijk.norm()) * (vijk - rijk * (vijk.dot(rijk) / rijk.norm2())) };
421 [[nodiscard]] std::tuple<RVec, RVec>
422 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE3FAD>(ArrayRef<const RVec> positions,
423 ArrayRef<const RVec> velocities) const
425 // Note: a = d * cos(theta)
426 // b = d * sin(theta)
427 const auto& a = parameters[0];
428 const auto& b = parameters[1];
429 const auto& ri = positions[constructingAtomIdx.at(0)];
430 const auto& rj = positions[constructingAtomIdx.at(1)];
431 const auto& rk = positions[constructingAtomIdx.at(2)];
432 const auto& vi = velocities[constructingAtomIdx.at(0)];
433 const auto& vj = velocities[constructingAtomIdx.at(1)];
434 const auto& vk = velocities[constructingAtomIdx.at(2)];
435 const auto rij = rj - ri;
436 const auto rjk = rk - rj;
437 const auto vij = vj - vi;
438 const auto vjk = vk - vj;
440 const auto rPerp = rjk - rij * (rij.dot(rjk) / rij.norm2());
441 const auto dtRijNormRij = (1 / rij.norm()) * (vij - rij * (vij.dot(rij) / rij.norm2()));
442 const auto vPerp = vjk
444 * ((vij.dot(rjk) + rij.dot(vjk)) / rij.norm2()
445 - 2 * rij.dot(rjk) * rij.dot(vij) / rij.norm2() / rij.norm2())
446 - vij * (rij.dot(rjk) / rij.norm2());
447 const auto dtRPerpNormRPerp =
448 (1 / rPerp.norm()) * (vPerp - rPerp * (vPerp.dot(rPerp) / rPerp.norm2()));
450 return { ri + (a / rij.norm()) * rij + (b / rPerp.norm()) * rPerp,
451 vi + a * dtRijNormRij + b * dtRPerpNormRPerp };
454 [[nodiscard]] std::tuple<RVec, RVec>
455 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE3OUT>(ArrayRef<const RVec> positions,
456 ArrayRef<const RVec> velocities) const
458 const auto& a = parameters[0];
459 const auto& b = parameters[1];
460 const auto& c = parameters[2];
461 const auto& ri = positions[constructingAtomIdx.at(0)];
462 const auto& rj = positions[constructingAtomIdx.at(1)];
463 const auto& rk = positions[constructingAtomIdx.at(2)];
464 const auto& vi = velocities[constructingAtomIdx.at(0)];
465 const auto& vj = velocities[constructingAtomIdx.at(1)];
466 const auto& vk = velocities[constructingAtomIdx.at(2)];
467 const auto rij = rj - ri;
468 const auto rik = rk - ri;
469 const auto vij = vj - vi;
470 const auto vik = vk - vi;
472 return { ri + a * rij + b * rik + c * rij.cross(rik),
473 vi + a * vij + b * vik + c * (vij.cross(rik) + rij.cross(vik)) };
476 [[nodiscard]] std::tuple<RVec, RVec>
477 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE4FDN>(ArrayRef<const RVec> positions,
478 ArrayRef<const RVec> velocities) const
480 const auto& a = parameters[0];
481 const auto& b = parameters[1];
482 const auto& c = parameters[2];
483 const auto& ri = positions[constructingAtomIdx.at(0)];
484 const auto& rj = positions[constructingAtomIdx.at(1)];
485 const auto& rk = positions[constructingAtomIdx.at(2)];
486 const auto& rl = positions[constructingAtomIdx.at(3)];
487 const auto& vi = velocities[constructingAtomIdx.at(0)];
488 const auto& vj = velocities[constructingAtomIdx.at(1)];
489 const auto& vk = velocities[constructingAtomIdx.at(2)];
490 const auto& vl = velocities[constructingAtomIdx.at(3)];
491 const auto rij = rj - ri;
492 const auto rik = rk - ri;
493 const auto ril = rl - ri;
494 const auto vij = vj - vi;
495 const auto vik = vk - vi;
496 const auto vil = vl - vi;
498 const auto rja = a * rik - rij;
499 const auto rjb = b * ril - rij;
500 const auto rm = rja.cross(rjb);
502 const auto vja = a * vik - vij;
503 const auto vjb = b * vil - vij;
504 const auto vm = vja.cross(rjb) + rja.cross(vjb);
506 return { ri + (c / rm.norm()) * rm, vi + (c / rm.norm()) * (vm - rm * (vm.dot(rm) / rm.norm2())) };
510 [[nodiscard]] std::tuple<RVec, RVec>
511 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITEN>(ArrayRef<const RVec> positions,
512 ArrayRef<const RVec> velocities) const
514 const auto& ri = positions[constructingAtomIdx.at(0)];
515 const auto& vi = velocities[constructingAtomIdx.at(0)];
516 GMX_RELEASE_ASSERT(constructingAtomIdx.size() == parameters.size() - 1,
517 "VSITEN atom / parameters mismatch.");
522 const auto parameterSize = parameters.size();
523 for (auto idx = decltype(parameterSize){ 0 }; idx < parameterSize; idx++)
525 rSum += parameters[idx] * (positions[constructingAtomIdx[idx + 1]] - ri);
526 vSum += parameters[idx] * (velocities[constructingAtomIdx[idx + 1]] - vi);
529 return { ri + rSum, vi + vSum };
532 [[nodiscard]] std::tuple<RVec, RVec>
533 VirtualSiteTest::VirtualSite::calculate(ArrayRef<const RVec> constructingPositions,
534 ArrayRef<const RVec> constructingVelocities) const
539 return calculateVSite<F_VSITE1>(constructingPositions, constructingVelocities);
541 return calculateVSite<F_VSITE2>(constructingPositions, constructingVelocities);
543 return calculateVSite<F_VSITE2FD>(constructingPositions, constructingVelocities);
545 return calculateVSite<F_VSITE3>(constructingPositions, constructingVelocities);
547 return calculateVSite<F_VSITE3FD>(constructingPositions, constructingVelocities);
549 return calculateVSite<F_VSITE3FAD>(constructingPositions, constructingVelocities);
551 return calculateVSite<F_VSITE3OUT>(constructingPositions, constructingVelocities);
553 return calculateVSite<F_VSITE4FDN>(constructingPositions, constructingVelocities);
555 return calculateVSite<F_VSITEN>(constructingPositions, constructingVelocities);
556 default: throw NotImplementedError("Unknown virtual site type");
561 TEST(VirtualSiteVelocityTest, ReferenceIsCorrect)
563 // Test is too sensitive to run in single precision
564 if constexpr (GMX_DOUBLE)
566 VirtualSiteTest::checkReferenceImplementation();
570 TEST_P(VirtualSiteTest, WithinToleranceOfReference)
572 const auto& params = GetParam();
573 const auto& integrator = std::get<0>(params);
574 const auto& tcoupling = std::get<1>(params);
575 const auto& pcoupling = std::get<2>(params);
576 const real timeStep = 0.001;
577 const auto& simulationName = "vsite_test";
579 if (integrator == "md-vv" && pcoupling == "parrinello-rahman")
581 // Parrinello-Rahman is not implemented in md-vv
585 if ((integrator == "sd" || integrator == "bd") && tcoupling != "no")
587 // bd and sd handle temperature coupling implicitly and would set tcoupling to "no" anyway
592 auto mdpFieldValues = prepareMdpFieldValues(simulationName, integrator, tcoupling, pcoupling);
593 mdpFieldValues["nsteps"] = "8";
594 mdpFieldValues["nstxout"] = "4";
595 mdpFieldValues["nstvout"] = "4";
596 mdpFieldValues["dt"] = toString(timeStep);
597 mdpFieldValues["constraints"] = "none";
598 if (pcoupling == "parrinello-rahman")
600 mdpFieldValues["tau-p"] = "2";
604 runner_.useTopGroAndNdxFromDatabase(simulationName);
605 runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
610 TopologyInformation topologyInformation;
611 topologyInformation.fillFromInputFile(runner_.tprFileName_);
612 const auto virtualSites = vSiteList(topologyInformation);
614 // This is in line with other tests (e.g. exact continuation, rerun), which
615 // never reach the same reproducibility for BD as for the other integrators.
616 const auto tolerance =
617 (integrator == "bd") ? relativeToleranceAsUlp(1.0, 100) : defaultRealTolerance();
619 checkVirtualSitesAgainstReferenceImplementation(
620 runner_.fullPrecisionTrajectoryFileName_, virtualSites, tolerance);
623 INSTANTIATE_TEST_CASE_P(
624 VelocitiesConformToExpectations,
626 ::testing::Combine(::testing::Values("md", "md-vv", "sd", "bd"),
627 ::testing::Values("no", "v-rescale", "nose-hoover"),
628 ::testing::Values("no", "c-rescale", "parrinello-rahman")));
631 } // namespace gmx::test