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