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