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