src/gromacs/listed_forces/tests/bonded.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,2019,2020,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 /*! \internal \file
  36  * \brief
  37  * Implements test of bonded force routines
  38  *
  39  *
  40  * \todo We should re-work this. For example, a harmonic bond
  41  * has so few computations that force components should be
  42  * accurate to a small *and computed* relative error.
  43  *
  44  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  45  * \ingroup module_listed_forces
  46  */
  47 #include "gmxpre.h"
  48
  49 #include "gromacs/listed_forces/bonded.h"
  50
  51 #include <cmath>
  52
  53 #include <memory>
  54 #include <unordered_map>
  55
  56 #include <gtest/gtest.h>
  57
  58 #include "gromacs/listed_forces/listed_forces.h"
  59 #include "gromacs/math/paddedvector.h"
  60 #include "gromacs/math/units.h"
  61 #include "gromacs/math/vec.h"
  62 #include "gromacs/math/vectypes.h"
  63 #include "gromacs/mdtypes/mdatom.h"
  64 #include "gromacs/pbcutil/ishift.h"
  65 #include "gromacs/pbcutil/pbc.h"
  66 #include "gromacs/topology/idef.h"
  67 #include "gromacs/utility/enumerationhelpers.h"
  68 #include "gromacs/utility/strconvert.h"
  69 #include "gromacs/utility/stringstream.h"
  70 #include "gromacs/utility/textwriter.h"
  71
  72 #include "testutils/refdata.h"
  73 #include "testutils/testasserts.h"
  74
  75 namespace gmx
  76 {
  77 namespace test
  78 {
  79 namespace
  80 {
  81
  82 //! Number of atoms used in these tests.
  83 constexpr int c_numAtoms = 4;
  84
  85 /*! \brief Output from bonded kernels
  86  *
  87  * \todo Later this might turn into the actual output struct. */
  88 struct OutputQuantities
  89 {
  90     //! Energy of this interaction
  91     real energy = 0;
  92     //! Derivative with respect to lambda
  93     real dvdlambda = 0;
  94     //! Shift vectors
  95     rvec fshift[c_numShiftVectors] = { { 0 } };
  96     //! Forces
  97     alignas(GMX_REAL_MAX_SIMD_WIDTH * sizeof(real)) rvec4 f[c_numAtoms] = { { 0 } };
  98 };
  99
 100 /*! \brief Utility to check the output from bonded tests
 101  *
 102  * \param[in] checker Reference checker
 103  * \param[in] output  The output from the test to check
 104  * \param[in] bondedKernelFlavor  Flavor for determining what output to check
 105  */
 106 void checkOutput(TestReferenceChecker*    checker,
 107                  const OutputQuantities&  output,
 108                  const BondedKernelFlavor bondedKernelFlavor)
 109 {
 110     if (computeEnergy(bondedKernelFlavor))
 111     {
 112         checker->checkReal(output.energy, "Epot ");
 113         // Should still be zero when not doing FEP, so may as well test it.
 114         checker->checkReal(output.dvdlambda, "dVdlambda ");
 115     }
 116     checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
 117 }
 118
 119 /*! \brief Input structure for listed forces tests
 120  */
 121 struct iListInput
 122 {
 123 public:
 124     //! Function type
 125     int ftype = -1;
 126     //! Tolerance for float evaluation
 127     float ftoler = 1e-6;
 128     //! Tolerance for double evaluation
 129     double dtoler = 1e-8;
 130     //! Do free energy perturbation?
 131     bool fep = false;
 132     //! Interaction parameters
 133     t_iparams iparams = { { 0 } };
 134
 135     friend std::ostream& operator<<(std::ostream& out, const iListInput& input);
 136
 137     //! Constructor
 138     iListInput() {}
 139
 140     /*! \brief Constructor with tolerance
 141      *
 142      * \param[in] ftol Single precision tolerance
 143      * \param[in] dtol Double precision tolerance
 144      */
 145     iListInput(float ftol, double dtol)
 146     {
 147         ftoler = ftol;
 148         dtoler = dtol;
 149     }
 150     /*! \brief Set parameters for harmonic potential
 151      *
 152      * Free energy perturbation is turned on when A
 153      * and B parameters are different.
 154      * \param[in] ft  Function type
 155      * \param[in] rA  Equilibrium value A
 156      * \param[in] krA Force constant A
 157      * \param[in] rB  Equilibrium value B
 158      * \param[in] krB Force constant B
 159      * \return The structure itself.
 160      */
 161     iListInput setHarmonic(int ft, real rA, real krA, real rB, real krB)
 162     {
 163         iparams.harmonic.rA  = rA;
 164         iparams.harmonic.rB  = rB;
 165         iparams.harmonic.krA = krA;
 166         iparams.harmonic.krB = krB;
 167         ftype                = ft;
 168         fep                  = (rA != rB || krA != krB);
 169         return *this;
 170     }
 171     /*! \brief Set parameters for harmonic potential
 172      *
 173      * \param[in] ft  Function type
 174      * \param[in] rA  Equilibrium value
 175      * \param[in] krA Force constant
 176      * \return The structure itself.
 177      */
 178     iListInput setHarmonic(int ft, real rA, real krA) { return setHarmonic(ft, rA, krA, rA, krA); }
 179     /*! \brief Set parameters for cubic potential
 180      *
 181      * \param[in] b0   Equilibrium bond length
 182      * \param[in] kb   Harmonic force constant
 183      * \param[in] kcub Cubic force constant
 184      * \return The structure itself.
 185      */
 186     iListInput setCubic(real b0, real kb, real kcub)
 187     {
 188         ftype              = F_CUBICBONDS;
 189         iparams.cubic.b0   = b0;
 190         iparams.cubic.kb   = kb;
 191         iparams.cubic.kcub = kcub;
 192         return *this;
 193     }
 194     /*! \brief Set parameters for morse potential
 195      *
 196      * Free energy perturbation is turned on when A
 197      * and B parameters are different.
 198      * \param[in] b0A   Equilibrium value A
 199      * \param[in] cbA   Force constant A
 200      * \param[in] betaA Steepness parameter A
 201      * \param[in] b0B   Equilibrium value B
 202      * \param[in] cbB   Force constant B
 203      * \param[in] betaB Steepness parameter B
 204      * \return The structure itself.
 205      */
 206     iListInput setMorse(real b0A, real cbA, real betaA, real b0B, real cbB, real betaB)
 207     {
 208         ftype               = F_MORSE;
 209         iparams.morse.b0A   = b0A;
 210         iparams.morse.cbA   = cbA;
 211         iparams.morse.betaA = betaA;
 212         iparams.morse.b0B   = b0B;
 213         iparams.morse.cbB   = cbB;
 214         iparams.morse.betaB = betaB;
 215         fep                 = (b0A != b0B || cbA != cbB || betaA != betaB);
 216         return *this;
 217     }
 218     /*! \brief Set parameters for morse potential
 219      *
 220      * \param[in] b0A   Equilibrium value
 221      * \param[in] cbA   Force constant
 222      * \param[in] betaA Steepness parameter
 223      * \return The structure itself.
 224      */
 225     iListInput setMorse(real b0A, real cbA, real betaA)
 226     {
 227         return setMorse(b0A, cbA, betaA, b0A, cbA, betaA);
 228     }
 229     /*! \brief Set parameters for fene potential
 230      *
 231      * \param[in] bm Equilibrium bond length
 232      * \param[in] kb Force constant
 233      * \return The structure itself.
 234      */
 235     iListInput setFene(real bm, real kb)
 236     {
 237         ftype           = F_FENEBONDS;
 238         iparams.fene.bm = bm;
 239         iparams.fene.kb = kb;
 240         return *this;
 241     }
 242     /*! \brief Set parameters for linear angle potential
 243      *
 244      * Free energy perturbation is turned on when A
 245      * and B parameters are different.
 246      * \param[in] klinA Force constant A
 247      * \param[in] aA    The position of the central atom A
 248      * \param[in] klinB Force constant B
 249      * \param[in] aB    The position of the central atom B
 250      * \return The structure itself.
 251      */
 252     iListInput setLinearAngle(real klinA, real aA, real klinB, real aB)
 253     {
 254         ftype                  = F_LINEAR_ANGLES;
 255         iparams.linangle.klinA = klinA;
 256         iparams.linangle.aA    = aA;
 257         iparams.linangle.klinB = klinB;
 258         iparams.linangle.aB    = aB;
 259         fep                    = (klinA != klinB || aA != aB);
 260         return *this;
 261     }
 262     /*! \brief Set parameters for linear angle potential
 263      *
 264      * \param[in] klinA Force constant
 265      * \param[in] aA    The position of the central atom
 266      * \return The structure itself.
 267      */
 268     iListInput setLinearAngle(real klinA, real aA) { return setLinearAngle(klinA, aA, klinA, aA); }
 269     /*! \brief Set parameters for Urey Bradley potential
 270      *
 271      * Free energy perturbation is turned on when A
 272      * and B parameters are different.
 273      * \param[in] thetaA  Equilibrium angle A
 274      * \param[in] kthetaA Force constant A
 275      * \param[in] r13A    The distance between i and k atoms A
 276      * \param[in] kUBA    The force constant for 1-3 distance A
 277      * \param[in] thetaB  Equilibrium angle B
 278      * \param[in] kthetaB Force constant B
 279      * \param[in] r13B    The distance between i and k atoms B
 280      * \param[in] kUBB    The force constant for 1-3 distance B
 281      * \return The structure itself.
 282      */
 283     iListInput
 284     setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA, real thetaB, real kthetaB, real r13B, real kUBB)
 285     {
 286         ftype               = F_UREY_BRADLEY;
 287         iparams.u_b.thetaA  = thetaA;
 288         iparams.u_b.kthetaA = kthetaA;
 289         iparams.u_b.r13A    = r13A;
 290         iparams.u_b.kUBA    = kUBA;
 291         iparams.u_b.thetaB  = thetaB;
 292         iparams.u_b.kthetaB = kthetaB;
 293         iparams.u_b.r13B    = r13B;
 294         iparams.u_b.kUBB    = kUBB;
 295         fep = (thetaA != thetaB || kthetaA != kthetaB || r13A != r13B || kUBA != kUBB);
 296         return *this;
 297     }
 298     /*! \brief Set parameters for Urey Bradley potential
 299      *
 300      * \param[in] thetaA  Equilibrium angle
 301      * \param[in] kthetaA Force constant
 302      * \param[in] r13A    The distance between i and k atoms
 303      * \param[in] kUBA    The force constant for 1-3 distance
 304      * \return The structure itself.
 305      */
 306     iListInput setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA)
 307     {
 308         return setUreyBradley(thetaA, kthetaA, r13A, kUBA, thetaA, kthetaA, r13A, kUBA);
 309     }
 310     /*! \brief Set parameters for Cross Bond Bonds potential
 311      *
 312      * \param[in] r1e  First bond length i-j
 313      * \param[in] r2e  Second bond length i-k
 314      * \param[in] krr  The force constant
 315      * \return The structure itself.
 316      */
 317     iListInput setCrossBondBonds(real r1e, real r2e, real krr)
 318     {
 319         ftype                = F_CROSS_BOND_BONDS;
 320         iparams.cross_bb.r1e = r1e;
 321         iparams.cross_bb.r2e = r2e;
 322         iparams.cross_bb.krr = krr;
 323         return *this;
 324     }
 325     /*! \brief Set parameters for Cross Bond Angles potential
 326      *
 327      * \param[in] r1e  First bond length i-j
 328      * \param[in] r2e  Second bond length j-k
 329      * \param[in] r3e  Third bond length i-k
 330      * \param[in] krt  The force constant
 331      * \return The structure itself.
 332      */
 333     iListInput setCrossBondAngles(real r1e, real r2e, real r3e, real krt)
 334     {
 335         ftype                = F_CROSS_BOND_ANGLES;
 336         iparams.cross_ba.r1e = r1e;
 337         iparams.cross_ba.r2e = r2e;
 338         iparams.cross_ba.r3e = r3e;
 339         iparams.cross_ba.krt = krt;
 340         return *this;
 341     }
 342     /*! \brief Set parameters for Quartic Angles potential
 343      *
 344      * \param[in] theta Angle
 345      * \param[in] c     Array of parameters
 346      * \return The structure itself.
 347      */
 348     iListInput setQuarticAngles(real theta, const real c[5])
 349     {
 350         ftype                = F_QUARTIC_ANGLES;
 351         iparams.qangle.theta = theta;
 352         iparams.qangle.c[0]  = c[0];
 353         iparams.qangle.c[1]  = c[1];
 354         iparams.qangle.c[2]  = c[2];
 355         iparams.qangle.c[3]  = c[3];
 356         iparams.qangle.c[4]  = c[4];
 357         return *this;
 358     }
 359     /*! \brief Set parameters for proper dihedrals potential
 360      *
 361      * Free energy perturbation is turned on when A
 362      * and B parameters are different.
 363      * \param[in] ft   Function type
 364      * \param[in] phiA Dihedral angle A
 365      * \param[in] cpA  Force constant A
 366      * \param[in] mult Multiplicity of the angle
 367      * \param[in] phiB Dihedral angle B
 368      * \param[in] cpB  Force constant B
 369      * \return The structure itself.
 370      */
 371     iListInput setPDihedrals(int ft, real phiA, real cpA, int mult, real phiB, real cpB)
 372     {
 373         ftype              = ft;
 374         iparams.pdihs.phiA = phiA;
 375         iparams.pdihs.cpA  = cpA;
 376         iparams.pdihs.phiB = phiB;
 377         iparams.pdihs.cpB  = cpB;
 378         iparams.pdihs.mult = mult;
 379         fep                = (phiA != phiB || cpA != cpB);
 380         return *this;
 381     }
 382     /*! \brief Set parameters for proper dihedrals potential
 383      *
 384      * \param[in] ft   Function type
 385      * \param[in] phiA Dihedral angle
 386      * \param[in] cpA  Force constant
 387      * \param[in] mult Multiplicity of the angle
 388      * \return The structure itself.
 389      */
 390     iListInput setPDihedrals(int ft, real phiA, real cpA, int mult)
 391     {
 392         return setPDihedrals(ft, phiA, cpA, mult, phiA, cpA);
 393     }
 394     /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
 395      *
 396      * Free energy perturbation is turned on when A
 397      * and B parameters are different.
 398      * \param[in] rbcA Force constants A
 399      * \param[in] rbcB Force constants B
 400      * \return The structure itself.
 401      */
 402     iListInput setRbDihedrals(const real rbcA[NR_RBDIHS], const real rbcB[NR_RBDIHS])
 403     {
 404         ftype = F_RBDIHS;
 405         fep   = false;
 406         for (int i = 0; i < NR_RBDIHS; i++)
 407         {
 408             iparams.rbdihs.rbcA[i] = rbcA[i];
 409             iparams.rbdihs.rbcB[i] = rbcB[i];
 410             fep                    = fep || (rbcA[i] != rbcB[i]);
 411         }
 412         return *this;
 413     }
 414     /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
 415      *
 416      * \param[in] rbc Force constants
 417      * \return The structure itself.
 418      */
 419     iListInput setRbDihedrals(const real rbc[NR_RBDIHS]) { return setRbDihedrals(rbc, rbc); }
 420     /*! \brief Set parameters for Polarization
 421      *
 422      * \param[in] alpha Polarizability
 423      * \return The structure itself.
 424      */
 425     iListInput setPolarization(real alpha)
 426     {
 427         ftype                  = F_POLARIZATION;
 428         fep                    = false;
 429         iparams.polarize.alpha = alpha;
 430         return *this;
 431     }
 432     /*! \brief Set parameters for Anharmonic Polarization
 433      *
 434      * \param[in] alpha Polarizability (nm^3)
 435      * \param[in] drcut The cut-off distance (nm) after which the
 436      *                  fourth power kicks in
 437      * \param[in] khyp  The force constant for the fourth power
 438      * \return The structure itself.
 439      */
 440     iListInput setAnharmPolarization(real alpha, real drcut, real khyp)
 441     {
 442         ftype                         = F_ANHARM_POL;
 443         fep                           = false;
 444         iparams.anharm_polarize.alpha = alpha;
 445         iparams.anharm_polarize.drcut = drcut;
 446         iparams.anharm_polarize.khyp  = khyp;
 447         return *this;
 448     }
 449     /*! \brief Set parameters for Thole Polarization
 450      *
 451      * \param[in] a      Thole factor
 452      * \param[in] alpha1 Polarizability 1 (nm^3)
 453      * \param[in] alpha2 Polarizability 2 (nm^3)
 454      * \param[in] rfac   Distance factor
 455      * \return The structure itself.
 456      */
 457     iListInput setTholePolarization(real a, real alpha1, real alpha2, real rfac)
 458     {
 459         ftype                = F_THOLE_POL;
 460         fep                  = false;
 461         iparams.thole.a      = a;
 462         iparams.thole.alpha1 = alpha1;
 463         iparams.thole.alpha2 = alpha2;
 464         iparams.thole.rfac   = rfac;
 465         return *this;
 466     }
 467     /*! \brief Set parameters for Water Polarization
 468      *
 469      * \param[in] alpha_x Polarizability X (nm^3)
 470      * \param[in] alpha_y Polarizability Y (nm^3)
 471      * \param[in] alpha_z Polarizability Z (nm^3)
 472      * \param[in] rOH     Oxygen-Hydrogen distance
 473      * \param[in] rHH     Hydrogen-Hydrogen distance
 474      * \param[in] rOD     Oxygen-Dummy distance
 475      * \return The structure itself.
 476      */
 477     iListInput setWaterPolarization(real alpha_x, real alpha_y, real alpha_z, real rOH, real rHH, real rOD)
 478     {
 479         ftype             = F_WATER_POL;
 480         fep               = false;
 481         iparams.wpol.al_x = alpha_x;
 482         iparams.wpol.al_y = alpha_y;
 483         iparams.wpol.al_z = alpha_z;
 484         iparams.wpol.rOH  = rOH;
 485         iparams.wpol.rHH  = rHH;
 486         iparams.wpol.rOD  = rOD;
 487         return *this;
 488     }
 489 };
 490
 491 //! Prints the interaction and parameters to a stream
 492 std::ostream& operator<<(std::ostream& out, const iListInput& input)
 493 {
 494     using std::endl;
 495     out << "Function type " << input.ftype << " called " << interaction_function[input.ftype].name
 496         << " ie. labelled '" << interaction_function[input.ftype].longname << "' in an energy file"
 497         << endl;
 498
 499     // Organize to print the legacy C union t_iparams, whose
 500     // relevant contents vary with ftype.
 501     StringOutputStream stream;
 502     {
 503         TextWriter writer(&stream);
 504         printInteractionParameters(&writer, input.ftype, input.iparams);
 505     }
 506     out << "Function parameters " << stream.toString();
 507     out << "Parameters trigger FEP? " << (input.fep ? "true" : "false") << endl;
 508     return out;
 509 }
 510
 511 /*! \brief Utility to fill iatoms struct
 512  *
 513  * \param[in]  ftype  Function type
 514  * \param[out] iatoms Pointer to iatoms struct
 515  */
 516 void fillIatoms(int ftype, std::vector<t_iatom>* iatoms)
 517 {
 518     std::unordered_map<int, std::vector<int>> ia = { { 2, { 0, 0, 1, 0, 1, 2, 0, 2, 3 } },
 519                                                      { 3, { 0, 0, 1, 2, 0, 1, 2, 3 } },
 520                                                      { 4, { 0, 0, 1, 2, 3 } },
 521                                                      { 5, { 0, 0, 1, 2, 3, 0 } } };
 522     EXPECT_TRUE(ftype >= 0 && ftype < F_NRE);
 523     int nral = interaction_function[ftype].nratoms;
 524     for (auto& i : ia[nral])
 525     {
 526         iatoms->push_back(i);
 527     }
 528 }
 529
 530 class ListedForcesTest :
 531     public ::testing::TestWithParam<std::tuple<iListInput, PaddedVector<RVec>, PbcType>>
 532 {
 533 protected:
 534     matrix                 box_;
 535     t_pbc                  pbc_;
 536     PaddedVector<RVec>     x_;
 537     PbcType                pbcType_;
 538     iListInput             input_;
 539     TestReferenceData      refData_;
 540     TestReferenceChecker   checker_;
 541     FloatingPointTolerance shiftForcesTolerance_ = defaultRealTolerance();
 542     ListedForcesTest() : checker_(refData_.rootChecker())
 543     {
 544         input_   = std::get<0>(GetParam());
 545         x_       = std::get<1>(GetParam());
 546         pbcType_ = std::get<2>(GetParam());
 547         clear_mat(box_);
 548         box_[0][0] = box_[1][1] = box_[2][2] = 1.5;
 549         set_pbc(&pbc_, pbcType_, box_);
 550         // We need quite specific tolerances here since angle functions
 551         // etc. are not very precise and reproducible.
 552         test::FloatingPointTolerance tolerance(test::FloatingPointTolerance(
 553                 input_.ftoler, input_.dtoler, 1.0e-6, 1.0e-12, 10000, 100, false));
 554         checker_.setDefaultTolerance(tolerance);
 555         // The SIMD acos() is only accurate to 2-3 ULP, so the angles
 556         // computed by it and the non-SIMD code paths (that use
 557         // std::acos) differ by enough to require quite large
 558         // tolerances for the shift forces in mixed precision.
 559         float singleShiftForcesAbsoluteTolerance =
 560                 ((input_.ftype == F_POLARIZATION) || (input_.ftype == F_ANHARM_POL)
 561                                  || (IS_ANGLE(input_.ftype))
 562                          ? 5e-3
 563                          : 5e-5);
 564         // Note that std::numeric_limits isn't required by the standard to
 565         // have an implementation for uint64_t(!) but this is likely to
 566         // work because that type is likely to be a typedef for one of
 567         // the other numerical types that happens to be 64-bits wide.
 568         shiftForcesTolerance_ = FloatingPointTolerance(singleShiftForcesAbsoluteTolerance,
 569                                                        1e-8,
 570                                                        1e-6,
 571                                                        1e-12,
 572                                                        std::numeric_limits<uint64_t>::max(),
 573                                                        std::numeric_limits<uint64_t>::max(),
 574                                                        false);
 575     }
 576     void testOneIfunc(TestReferenceChecker* checker, const std::vector<t_iatom>& iatoms, const real lambda)
 577     {
 578         SCOPED_TRACE(std::string("Testing PBC type: ") + c_pbcTypeNames[pbcType_]);
 579         std::vector<int>  ddgatindex = { 0, 1, 2, 3 };
 580         std::vector<real> charge     = { 1.5, -2.0, 1.5, -1.0 };
 581         /* Here we run both the standard, plain-C force+shift-forces+energy+free-energy
 582          * kernel flavor and the potentially optimized, with SIMD and less output,
 583          * force only kernels. Note that we also run the optimized kernel for free-energy
 584          * input when lambda=0, as the force output should match the non free-energy case.
 585          */
 586         std::vector<BondedKernelFlavor> flavors = { BondedKernelFlavor::ForcesAndVirialAndEnergy };
 587         if (!input_.fep || lambda == 0)
 588         {
 589             flavors.push_back(BondedKernelFlavor::ForcesSimdWhenAvailable);
 590         }
 591         for (const auto flavor : flavors)
 592         {
 593             SCOPED_TRACE("Testing bonded kernel flavor: " + c_bondedKernelFlavorStrings[flavor]);
 594             OutputQuantities output;
 595             output.energy = calculateSimpleBond(input_.ftype,
 596                                                 iatoms.size(),
 597                                                 iatoms.data(),
 598                                                 &input_.iparams,
 599                                                 as_rvec_array(x_.data()),
 600                                                 output.f,
 601                                                 output.fshift,
 602                                                 &pbc_,
 603                                                 lambda,
 604                                                 &output.dvdlambda,
 605                                                 charge,
 606                                                 /* struct t_fcdata * */ nullptr,
 607                                                 nullptr,
 608                                                 nullptr,
 609                                                 ddgatindex.data(),
 610                                                 flavor);
 611             // Internal consistency test of both test input
 612             // and bonded functions.
 613             EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0))) << "dvdlambda was " << output.dvdlambda;
 614             checkOutput(checker, output, flavor);
 615             auto shiftForcesChecker = checker->checkCompound("Shift-Forces", "Shift-forces");
 616             if (computeVirial(flavor))
 617             {
 618                 shiftForcesChecker.setDefaultTolerance(shiftForcesTolerance_);
 619                 shiftForcesChecker.checkVector(output.fshift[c_centralShiftIndex], "Central");
 620             }
 621             else
 622             {
 623                 // Permit omitting to compare shift forces with
 624                 // reference data when that is useless.
 625                 shiftForcesChecker.disableUnusedEntriesCheck();
 626             }
 627         }
 628     }
 629     void testIfunc()
 630     {
 631         TestReferenceChecker thisChecker =
 632                 checker_.checkCompound("FunctionType", interaction_function[input_.ftype].name)
 633                         .checkCompound("FEP", (input_.fep ? "Yes" : "No"));
 634         std::vector<t_iatom> iatoms;
 635         fillIatoms(input_.ftype, &iatoms);
 636         if (input_.fep)
 637         {
 638             const int numLambdas = 3;
 639             for (int i = 0; i < numLambdas; ++i)
 640             {
 641                 const real lambda       = i / (numLambdas - 1.0);
 642                 auto       valueChecker = thisChecker.checkCompound("Lambda", toString(lambda));
 643                 testOneIfunc(&valueChecker, iatoms, lambda);
 644             }
 645         }
 646         else
 647         {
 648             testOneIfunc(&thisChecker, iatoms, 0.0);
 649         }
 650     }
 651 };
 652
 653 TEST_P(ListedForcesTest, Ifunc)
 654 {
 655     testIfunc();
 656 }
 657
 658 //! Function types for testing bonds. Add new terms at the end.
 659 std::vector<iListInput> c_InputBonds = {
 660     { iListInput(2e-6F, 1e-8).setHarmonic(F_BONDS, 0.15, 500.0) },
 661     { iListInput(2e-6F, 1e-8).setHarmonic(F_BONDS, 0.15, 500.0, 0.17, 400.0) },
 662     { iListInput(1e-4F, 1e-8).setHarmonic(F_G96BONDS, 0.15, 50.0) },
 663     { iListInput().setHarmonic(F_G96BONDS, 0.15, 50.0, 0.17, 40.0) },
 664     { iListInput().setCubic(0.16, 50.0, 2.0) },
 665     { iListInput(2e-6F, 1e-8).setMorse(0.15, 50.0, 2.0, 0.17, 40.0, 1.6) },
 666     { iListInput(2e-6F, 1e-8).setMorse(0.15, 30.0, 2.7) },
 667     { iListInput().setFene(0.4, 5.0) }
 668 };
 669
 670 //! Constants for Quartic Angles
 671 const real cQuarticAngles[5] = { 1.1, 2.3, 4.6, 7.8, 9.2 };
 672
 673 //! Function types for testing angles. Add new terms at the end.
 674 std::vector<iListInput> c_InputAngles = {
 675     { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 100.0, 50.0) },
 676     { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 100.15, 50.0, 95.0, 30.0) },
 677     { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES, 100.0, 50.0) },
 678     { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES, 100.0, 50.0, 95.0, 30.0) },
 679     { iListInput().setLinearAngle(50.0, 0.4) },
 680     { iListInput().setLinearAngle(50.0, 0.4, 40.0, 0.6) },
 681     { iListInput(2e-6, 1e-8).setCrossBondBonds(0.8, 0.7, 45.0) },
 682     { iListInput(3e-6, 1e-8).setCrossBondAngles(0.8, 0.7, 0.3, 45.0) },
 683     { iListInput(2e-2, 1e-8).setUreyBradley(950.0, 46.0, 0.3, 5.0) },
 684     { iListInput(2e-2, 1e-8).setUreyBradley(100.0, 45.0, 0.3, 5.0, 90.0, 47.0, 0.32, 7.0) },
 685     { iListInput(2e-3, 1e-8).setQuarticAngles(87.0, cQuarticAngles) }
 686 };
 687
 688 //! Constants for Ryckaert-Bellemans A
 689 const real rbcA[NR_RBDIHS] = { -5.35, 13.6, 8.4, -16.7, 0.3, 12.4 };
 690
 691 //! Constants for Ryckaert-Bellemans B
 692 const real rbcB[NR_RBDIHS] = { -6.35, 12.6, 8.1, -10.7, 0.9, 15.4 };
 693
 694 //! Constants for Ryckaert-Bellemans without FEP
 695 const real rbc[NR_RBDIHS] = { -7.35, 13.6, 8.4, -16.7, 1.3, 12.4 };
 696
 697 //! Function types for testing dihedrals. Add new terms at the end.
 698 std::vector<iListInput> c_InputDihs = {
 699     { iListInput(5e-4, 1e-8).setPDihedrals(F_PDIHS, -100.0, 10.0, 2, -80.0, 20.0) },
 700     { iListInput(1e-4, 1e-8).setPDihedrals(F_PDIHS, -105.0, 15.0, 2) },
 701     { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS, 100.0, 50.0) },
 702     { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS, 100.15, 50.0, 95.0, 30.0) },
 703     { iListInput(4e-4, 1e-8).setRbDihedrals(rbcA, rbcB) },
 704     { iListInput(4e-4, 1e-8).setRbDihedrals(rbc) }
 705 };
 706
 707 //! Function types for testing polarization. Add new terms at the end.
 708 std::vector<iListInput> c_InputPols = {
 709     { iListInput(2e-5, 1e-8).setPolarization(0.12) },
 710     { iListInput(2e-3, 1e-8).setAnharmPolarization(0.0013, 0.02, 1235.6) },
 711     { iListInput(1.4e-3, 1e-8).setTholePolarization(0.26, 0.07, 0.09, 1.6) },
 712     { iListInput(2e-3, 1e-8).setWaterPolarization(0.001, 0.0012, 0.0016, 0.095, 0.15, 0.02) },
 713 };
 714
 715 //! Function types for testing polarization. Add new terms at the end.
 716 std::vector<iListInput> c_InputRestraints = {
 717     { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRES, -100.0, 10.0, 2, -80.0, 20.0) },
 718     { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRES, -105.0, 15.0, 2) },
 719     { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRESZ, -100.0, 10.0, 2, -80.0, 20.0) },
 720     { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRESZ, -105.0, 15.0, 2) },
 721     { iListInput(2e-3, 1e-8).setHarmonic(F_RESTRANGLES, 100.0, 50.0) },
 722     { iListInput(2e-3, 1e-8).setHarmonic(F_RESTRANGLES, 100.0, 50.0, 110.0, 45.0) }
 723 };
 724
 725 //! Function types for testing bond with zero length, has zero reference length to make physical sense.
 726 std::vector<iListInput> c_InputBondsZeroLength = {
 727     { iListInput().setHarmonic(F_BONDS, 0.0, 500.0) },
 728 };
 729
 730 //! Function types for testing angles with zero angle, has zero reference angle to make physical sense.
 731 std::vector<iListInput> c_InputAnglesZeroAngle = {
 732     { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 0.0, 50.0) },
 733 };
 734
 735 } // namespace
 736 } // namespace test
 737
 738 //! Print an RVec to \c os
 739 static void PrintTo(const RVec& value, std::ostream* os)
 740 {
 741     *os << value[XX] << " " << value[YY] << " " << value[ZZ] << std::endl;
 742 }
 743
 744 //! Print a padded vector of RVec to \c os
 745 static void PrintTo(const PaddedVector<RVec>& vector, std::ostream* os)
 746 {
 747     if (vector.empty())
 748     {
 749         *os << "Empty vector" << std::endl;
 750     }
 751     else
 752     {
 753         *os << "Vector of RVec containing:" << std::endl;
 754         std::for_each(vector.begin(), vector.end(), [os](const RVec& v) { PrintTo(v, os); });
 755     }
 756 }
 757
 758 namespace test
 759 {
 760 namespace
 761 {
 762
 763 /*! \brief Coordinates for testing
 764  *
 765  * Taken from a butane molecule, so we have some
 766  * normal-sized bonds and angles to test.
 767  *
 768  * \todo Test also some weirder values */
 769 std::vector<PaddedVector<RVec>> c_coordinatesForTests = {
 770     { { 1.382, 1.573, 1.482 }, { 1.281, 1.559, 1.596 }, { 1.292, 1.422, 1.663 }, { 1.189, 1.407, 1.775 } }
 771 };
 772
 773 //! Coordinates for testing bonds with zero length
 774 std::vector<PaddedVector<RVec>> c_coordinatesForTestsZeroBondLength = {
 775     { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.005, 0.0, 0.1 }, { 0.005, 0.0, 0.1 } }
 776 };
 777
 778 //! Coordinates for testing bonds with zero length
 779 std::vector<PaddedVector<RVec>> c_coordinatesForTestsZeroAngle = {
 780     { { 0.005, 0.0, 0.1 }, { 0.0, 0.0, 0.0 }, { 0.005, 0.0, 0.1 }, { 0.5, 0.18, 0.22 } }
 781 };
 782
 783 //! PBC values for testing
 784 std::vector<PbcType> c_pbcForTests = { PbcType::No, PbcType::XY, PbcType::Xyz };
 785
 786 INSTANTIATE_TEST_SUITE_P(Bond,
 787                          ListedForcesTest,
 788                          ::testing::Combine(::testing::ValuesIn(c_InputBonds),
 789                                             ::testing::ValuesIn(c_coordinatesForTests),
 790                                             ::testing::ValuesIn(c_pbcForTests)));
 791
 792 INSTANTIATE_TEST_SUITE_P(Angle,
 793                          ListedForcesTest,
 794                          ::testing::Combine(::testing::ValuesIn(c_InputAngles),
 795                                             ::testing::ValuesIn(c_coordinatesForTests),
 796                                             ::testing::ValuesIn(c_pbcForTests)));
 797
 798 INSTANTIATE_TEST_SUITE_P(Dihedral,
 799                          ListedForcesTest,
 800                          ::testing::Combine(::testing::ValuesIn(c_InputDihs),
 801                                             ::testing::ValuesIn(c_coordinatesForTests),
 802                                             ::testing::ValuesIn(c_pbcForTests)));
 803
 804 INSTANTIATE_TEST_SUITE_P(Polarize,
 805                          ListedForcesTest,
 806                          ::testing::Combine(::testing::ValuesIn(c_InputPols),
 807                                             ::testing::ValuesIn(c_coordinatesForTests),
 808                                             ::testing::ValuesIn(c_pbcForTests)));
 809
 810 INSTANTIATE_TEST_SUITE_P(Restraints,
 811                          ListedForcesTest,
 812                          ::testing::Combine(::testing::ValuesIn(c_InputRestraints),
 813                                             ::testing::ValuesIn(c_coordinatesForTests),
 814                                             ::testing::ValuesIn(c_pbcForTests)));
 815
 816 INSTANTIATE_TEST_SUITE_P(BondZeroLength,
 817                          ListedForcesTest,
 818                          ::testing::Combine(::testing::ValuesIn(c_InputBondsZeroLength),
 819                                             ::testing::ValuesIn(c_coordinatesForTestsZeroBondLength),
 820                                             ::testing::ValuesIn(c_pbcForTests)));
 821
 822 INSTANTIATE_TEST_SUITE_P(AngleZero,
 823                          ListedForcesTest,
 824                          ::testing::Combine(::testing::ValuesIn(c_InputAnglesZeroAngle),
 825                                             ::testing::ValuesIn(c_coordinatesForTestsZeroAngle),
 826                                             ::testing::ValuesIn(c_pbcForTests)));
 827
 828 } // namespace
 829
 830 } // namespace test
 831
 832 } // namespace gmx