src/programs/mdrun/tests/virtualsites.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35
  36 /*! \internal \file
  37  * \brief Sanity checks for virtual sites
  38  *
  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.
  54  *
  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.
  61  *
  62  * \author Pascal Merz <pascal.merz@me.com>
  63  * \ingroup module_mdrun_integration_tests
  64  */
  65 #include "gmxpre.h"
  66
  67 #include "config.h"
  68
  69 #include "gromacs/topology/idef.h"
  70 #include "gromacs/topology/ifunc.h"
  71 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  72 #include "gromacs/utility/stringutil.h"
  73
  74 #include "testutils/mpitest.h"
  75 #include "testutils/simulationdatabase.h"
  76 #include "testutils/testmatchers.h"
  77
  78 #include "gromacs/utility/strconvert.h"
  79
  80 #include "moduletest.h"
  81 #include "simulatorcomparison.h"
  82 #include "trajectoryreader.h"
  83
  84 namespace gmx::test
  85 {
  86 namespace
  87 {
  88 using VirtualSiteTestParams = std::tuple<std::string, std::string, std::string>;
  89 class VirtualSiteTest : public MdrunTestFixture, public ::testing::WithParamInterface<VirtualSiteTestParams>
  90 {
  91 public:
  92     struct VirtualSite;
  93
  94     /*! \brief Check against reference implementation
  95      *
  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.
 100      *
 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.
 108      *
 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.
 112      */
 113     static void checkVirtualSitesAgainstReferenceImplementation(const std::string& trajectoryName,
 114                                                                 ArrayRef<const VirtualSite> virtualSites,
 115                                                                 FloatingPointTolerance tolerance)
 116     {
 117         SCOPED_TRACE(
 118                 "Checking virtual site positions and velocities against reference implementation.");
 119         TrajectoryFrameReader trajectoryFrameReader(trajectoryName);
 120         while (trajectoryFrameReader.readNextFrame())
 121         {
 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)
 129             {
 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));
 135             }
 136             EXPECT_THAT(refPositions, Pointwise(RVecEq(tolerance), positions));
 137             EXPECT_THAT(refVelocities, Pointwise(RVecEq(tolerance), velocities));
 138         }
 139     }
 140
 141     /*! \brief Check the reference implementation
 142      *
 143      * This tests that the reference implementation position and velocities correspond.
 144      * This is done by
 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.
 156      *
 157      * Maybe unused because this test runs only in double precision.
 158      */
 159     [[maybe_unused]] static void checkReferenceImplementation()
 160     {
 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 } }
 181         };
 182
 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++)
 191         {
 192             endPositions.emplace_back(startPositions[idx] + timeStep * halfStepVelocities[idx]);
 193             halfStepPositions.emplace_back(startPositions[idx]
 194                                            + real(0.5) * timeStep * halfStepVelocities[idx]);
 195         }
 196
 197         // Check that displacement equals the calculated velocities
 198         for (const auto& vSite : virtualSites)
 199         {
 200             SCOPED_TRACE(formatString("Checking %s", interaction_function[vSite.type].longname));
 201
 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
 204             // clang-format off
 205             GCC_DIAGNOSTIC_IGNORE(-Wunused-variable)
 206             // clang-format on
 207             /* Calculate start and end virtual position
 208              *
 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.
 212              */
 213             const auto [startVPosition, vVelocityUnused1] =
 214                     vSite.calculate(startPositions, halfStepVelocities);
 215             const auto [endVPosition1, vVelocityUnused2] =
 216                     vSite.calculate(endPositions, halfStepVelocities);
 217
 218             /* Calculate virtual velocity at half step using reference implementation
 219              *
 220              * The virtual positions are exact, but we don't need them.
 221              */
 222             const auto [halfStepVPositionUnused, halfStepVVelocity] =
 223                     vSite.calculate(halfStepPositions, halfStepVelocities);
 224             GCC_DIAGNOSTIC_RESET
 225
 226             // We can now integrate the virtual positions using the half step velocity
 227             const auto endVPosition2 = startVPosition + timeStep * halfStepVVelocity;
 228
 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!
 237              */
 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.";
 246
 247             EXPECT_REAL_EQ_TOL(displacement1[XX],
 248                                displacement2[XX],
 249                                relativeToleranceAsFloatingPoint(displacement1[XX], 1e-9));
 250             EXPECT_REAL_EQ_TOL(displacement1[YY],
 251                                displacement2[YY],
 252                                relativeToleranceAsFloatingPoint(displacement1[YY], 1e-9));
 253             EXPECT_REAL_EQ_TOL(displacement1[ZZ],
 254                                displacement2[ZZ],
 255                                relativeToleranceAsFloatingPoint(displacement1[ZZ], 1e-9));
 256         }
 257     }
 258
 259     //! Holds parameters of virtual site and allows calculation
 260     struct VirtualSite
 261     {
 262         //! Type of virtual site
 263         int type;
 264         //! Index of virtual site
 265         int atomIdx;
 266         //! Indices of constructing atoms
 267         std::vector<int> constructingAtomIdx;
 268         //! Construction parameters
 269         std::vector<real> parameters;
 270
 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;
 274
 275     private:
 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;
 280     };
 281
 282     /*! \brief Helper function returning a list of virtual sites from the topology
 283      *
 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.
 286      */
 287     static std::vector<VirtualSite> vSiteList(const TopologyInformation& topologyInformation)
 288     {
 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++)
 293         {
 294             const auto& interactionList = localTopology.idef.il.at(vsiteType);
 295             if (vsiteType == F_VSITE4FD || interactionList.empty())
 296             {
 297                 // 4FD is deprecated. Interaction list empty means system doesn't contain this type.
 298                 continue;
 299             }
 300             const int   numConstructingAtoms = interaction_function[vsiteType].nratoms - 1;
 301             const int   defaultIncrement     = numConstructingAtoms + 2;
 302             std::string indexString;
 303
 304             for (int i = 0; i < interactionList.size();)
 305             {
 306                 const int parameterIdx   = interactionList.iatoms[i];
 307                 const int virtualSiteIdx = interactionList.iatoms[i + 1];
 308                 if (!indexString.empty())
 309                 {
 310                     indexString += ", ";
 311                 }
 312                 indexString += toString(virtualSiteIdx);
 313
 314                 if (vsiteType == F_VSITEN)
 315                 {
 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)
 320                     {
 321                         vSite.constructingAtomIdx.push_back(interactionList.iatoms[j + 2]);
 322                         vSite.parameters.push_back(
 323                                 localTopology.idef.iparams[interactionList.iatoms[j]].vsiten.a);
 324                     }
 325                     virtualSites.push_back(std::move(vSite));
 326                     i += 3 * vSiteNConstructingAtoms;
 327                 }
 328                 else
 329                 {
 330                     virtualSites.emplace_back(VirtualSite{
 331                             vsiteType,
 332                             virtualSiteIdx,
 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;
 338                 }
 339             }
 340         }
 341         return virtualSites;
 342     }
 343 };
 344
 345 // check-source gets confused by these
 346 //! \cond
 347 template<>
 348 [[nodiscard]] std::tuple<RVec, RVec>
 349 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE1>(ArrayRef<const RVec> positions,
 350                                                        ArrayRef<const RVec> velocities) const
 351 {
 352     return { positions[constructingAtomIdx.at(0)], velocities[constructingAtomIdx.at(0)] };
 353 }
 354 template<>
 355 [[nodiscard]] std::tuple<RVec, RVec>
 356 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE2>(ArrayRef<const RVec> positions,
 357                                                        ArrayRef<const RVec> velocities) const
 358 {
 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)] };
 362 }
 363 template<>
 364 [[nodiscard]] std::tuple<RVec, RVec>
 365 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE2FD>(ArrayRef<const RVec> positions,
 366                                                          ArrayRef<const RVec> velocities) const
 367 {
 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;
 375
 376     return { ri + (a / rij.norm()) * rij,
 377              vi + (a / rij.norm()) * (vij - rij * (vij.dot(rij) / rij.norm2())) };
 378 }
 379 template<>
 380 [[nodiscard]] std::tuple<RVec, RVec>
 381 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE3>(ArrayRef<const RVec> positions,
 382                                                        ArrayRef<const RVec> velocities) const
 383 {
 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)];
 392
 393     return { (1 - a - b) * ri + a * rj + b * rk, (1 - a - b) * vi + a * vj + b * vk };
 394 }
 395 template<>
 396 [[nodiscard]] std::tuple<RVec, RVec>
 397 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE3FD>(ArrayRef<const RVec> positions,
 398                                                          ArrayRef<const RVec> velocities) const
 399 {
 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;
 412
 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;
 416
 417     return { ri + (b / rijk.norm()) * rijk,
 418              vi + (b / rijk.norm()) * (vijk - rijk * (vijk.dot(rijk) / rijk.norm2())) };
 419 }
 420 template<>
 421 [[nodiscard]] std::tuple<RVec, RVec>
 422 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE3FAD>(ArrayRef<const RVec> positions,
 423                                                           ArrayRef<const RVec> velocities) const
 424 {
 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;
 439
 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
 443                        - rij
 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()));
 449
 450     return { ri + (a / rij.norm()) * rij + (b / rPerp.norm()) * rPerp,
 451              vi + a * dtRijNormRij + b * dtRPerpNormRPerp };
 452 }
 453 template<>
 454 [[nodiscard]] std::tuple<RVec, RVec>
 455 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE3OUT>(ArrayRef<const RVec> positions,
 456                                                           ArrayRef<const RVec> velocities) const
 457 {
 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;
 471
 472     return { ri + a * rij + b * rik + c * rij.cross(rik),
 473              vi + a * vij + b * vik + c * (vij.cross(rik) + rij.cross(vik)) };
 474 }
 475 template<>
 476 [[nodiscard]] std::tuple<RVec, RVec>
 477 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITE4FDN>(ArrayRef<const RVec> positions,
 478                                                           ArrayRef<const RVec> velocities) const
 479 {
 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;
 497
 498     const auto rja = a * rik - rij;
 499     const auto rjb = b * ril - rij;
 500     const auto rm  = rja.cross(rjb);
 501
 502     const auto vja = a * vik - vij;
 503     const auto vjb = b * vil - vij;
 504     const auto vm  = vja.cross(rjb) + rja.cross(vjb);
 505
 506     return { ri + (c / rm.norm()) * rm, vi + (c / rm.norm()) * (vm - rm * (vm.dot(rm) / rm.norm2())) };
 507 }
 508
 509 template<>
 510 [[nodiscard]] std::tuple<RVec, RVec>
 511 VirtualSiteTest::VirtualSite::calculateVSite<F_VSITEN>(ArrayRef<const RVec> positions,
 512                                                        ArrayRef<const RVec> velocities) const
 513 {
 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.");
 518
 519     RVec rSum(0, 0, 0);
 520     RVec vSum(0, 0, 0);
 521
 522     const auto parameterSize = parameters.size();
 523     for (auto idx = decltype(parameterSize){ 0 }; idx < parameterSize; idx++)
 524     {
 525         rSum += parameters[idx] * (positions[constructingAtomIdx[idx + 1]] - ri);
 526         vSum += parameters[idx] * (velocities[constructingAtomIdx[idx + 1]] - vi);
 527     }
 528
 529     return { ri + rSum, vi + vSum };
 530 }
 531
 532 [[nodiscard]] std::tuple<RVec, RVec>
 533 VirtualSiteTest::VirtualSite::calculate(ArrayRef<const RVec> constructingPositions,
 534                                         ArrayRef<const RVec> constructingVelocities) const
 535 {
 536     switch (type)
 537     {
 538         case F_VSITE1:
 539             return calculateVSite<F_VSITE1>(constructingPositions, constructingVelocities);
 540         case F_VSITE2:
 541             return calculateVSite<F_VSITE2>(constructingPositions, constructingVelocities);
 542         case F_VSITE2FD:
 543             return calculateVSite<F_VSITE2FD>(constructingPositions, constructingVelocities);
 544         case F_VSITE3:
 545             return calculateVSite<F_VSITE3>(constructingPositions, constructingVelocities);
 546         case F_VSITE3FD:
 547             return calculateVSite<F_VSITE3FD>(constructingPositions, constructingVelocities);
 548         case F_VSITE3FAD:
 549             return calculateVSite<F_VSITE3FAD>(constructingPositions, constructingVelocities);
 550         case F_VSITE3OUT:
 551             return calculateVSite<F_VSITE3OUT>(constructingPositions, constructingVelocities);
 552         case F_VSITE4FDN:
 553             return calculateVSite<F_VSITE4FDN>(constructingPositions, constructingVelocities);
 554         case F_VSITEN:
 555             return calculateVSite<F_VSITEN>(constructingPositions, constructingVelocities);
 556         default: throw NotImplementedError("Unknown virtual site type");
 557     }
 558 }
 559 //! \endcond
 560
 561 TEST(VirtualSiteVelocityTest, ReferenceIsCorrect)
 562 {
 563     // Test is too sensitive to run in single precision
 564     if constexpr (GMX_DOUBLE)
 565     {
 566         VirtualSiteTest::checkReferenceImplementation();
 567     }
 568 }
 569
 570 TEST_P(VirtualSiteTest, WithinToleranceOfReference)
 571 {
 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";
 578
 579     if (integrator == "md-vv" && pcoupling == "parrinello-rahman")
 580     {
 581         // Parrinello-Rahman is not implemented in md-vv
 582         return;
 583     }
 584
 585     if ((integrator == "sd" || integrator == "bd") && tcoupling != "no")
 586     {
 587         // bd and sd handle temperature coupling implicitly and would set tcoupling to "no" anyway
 588         return;
 589     }
 590
 591     // Prepare mdp input
 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")
 599     {
 600         mdpFieldValues["tau-p"] = "2";
 601     }
 602
 603     // Run grompp
 604     runner_.useTopGroAndNdxFromDatabase(simulationName);
 605     runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
 606     runGrompp(&runner_);
 607     // Run mdrun
 608     runMdrun(&runner_);
 609
 610     TopologyInformation topologyInformation;
 611     topologyInformation.fillFromInputFile(runner_.tprFileName_);
 612     const auto virtualSites = vSiteList(topologyInformation);
 613
 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();
 618
 619     checkVirtualSitesAgainstReferenceImplementation(
 620             runner_.fullPrecisionTrajectoryFileName_, virtualSites, tolerance);
 621 }
 622
 623 INSTANTIATE_TEST_CASE_P(
 624         VelocitiesConformToExpectations,
 625         VirtualSiteTest,
 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")));
 629
 630 } // namespace
 631 } // namespace gmx::test