Remove unused thole polarization rfac parameter
[alexxy/gromacs.git] / src / gromacs / listed_forces / tests / bonded.cpp
index bb3a503d37a9e3983e2b0a36e43682d37bb6887d..04bf2e1eb7a9adf4951edb18783b09fd838871c4 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
  * \brief
  * Implements test of bonded force routines
  *
+ *
+ * \todo We should re-work this. For example, a harmonic bond
+ * has so few computations that force components should be
+ * accurate to a small *and computed* relative error.
+ *
  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  * \ingroup module_listed_forces
  */
 #include <gtest/gtest.h>
 
 #include "gromacs/listed_forces/listed_forces.h"
+#include "gromacs/math/paddedvector.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/math/vectypes.h"
+#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/pbcutil/ishift.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/idef.h"
+#include "gromacs/utility/enumerationhelpers.h"
 #include "gromacs/utility/strconvert.h"
+#include "gromacs/utility/stringstream.h"
+#include "gromacs/utility/textwriter.h"
 
 #include "testutils/refdata.h"
 #include "testutils/testasserts.h"
 
 namespace gmx
 {
+namespace test
+{
 namespace
 {
 
@@ -76,27 +88,31 @@ constexpr int c_numAtoms = 4;
 struct OutputQuantities
 {
     //! Energy of this interaction
-    real  energy         = 0;
+    real energy = 0;
     //! Derivative with respect to lambda
-    real  dvdlambda      = 0;
+    real dvdlambda = 0;
     //! Shift vectors
-    rvec  fshift[N_IVEC] = {{0}};
+    rvec fshift[c_numShiftVectors] = { { 0 } };
     //! Forces
-    rvec4 f[c_numAtoms]  = {{0}};
+    alignas(GMX_REAL_MAX_SIMD_WIDTH * sizeof(real)) rvec4 f[c_numAtoms] = { { 0 } };
 };
 
 /*! \brief Utility to check the output from bonded tests
  *
  * \param[in] checker Reference checker
  * \param[in] output  The output from the test to check
+ * \param[in] bondedKernelFlavor  Flavor for determining what output to check
  */
-void checkOutput(test::TestReferenceChecker *checker,
-                 const OutputQuantities     &output)
+void checkOutput(TestReferenceChecker*    checker,
+                 const OutputQuantities&  output,
+                 const BondedKernelFlavor bondedKernelFlavor)
 {
-    checker->checkReal(output.energy, "Epot ");
-    // Should still be zero if not doing FEP, so may as well test it.
-    checker->checkReal(output.dvdlambda, "dVdlambda ");
-    checker->checkVector(output.fshift[CENTRAL], "Central shift forces");
+    if (computeEnergy(bondedKernelFlavor))
+    {
+        checker->checkReal(output.energy, "Epot ");
+        // Should still be zero when not doing FEP, so may as well test it.
+        checker->checkReal(output.dvdlambda, "dVdlambda ");
+    }
     checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
 }
 
@@ -104,419 +120,548 @@ void checkOutput(test::TestReferenceChecker *checker,
  */
 struct iListInput
 {
-    public:
-        //! Function type
-        int       ftype = -1;
-        //! Tolerance for float evaluation
-        float     ftoler = 1e-6;
-        //! Tolerance for double evaluation
-        double    dtoler = 1e-8;
-        //! Do free energy perturbation?
-        bool      fep = false;
-        //! Interaction parameters
-        t_iparams iparams = {{ 0 }};
-
-        //! Constructor
-        iListInput() {}
-
-        /*! \brief Constructor with tolerance
-         *
-         * \param[in] ftol Single precision tolerance
-         * \param[in] dtol Double precision tolerance
-         */
-        iListInput(float ftol, double dtol)
-        {
-            ftoler = ftol;
-            dtoler = dtol;
-        }
-        /*! \brief Set parameters for harmonic potential
-         *
-         * Free energy perturbation is turned on when A
-         * and B parameters are different.
-         * \param[in] ft  Function type
-         * \param[in] rA  Equilibrium value A
-         * \param[in] krA Force constant A
-         * \param[in] rB  Equilibrium value B
-         * \param[in] krB Force constant B
-         * \return The structure itself.
-         */
-        iListInput setHarmonic(int ft, real rA, real krA, real rB, real krB)
-        {
-            iparams.harmonic.rA  = rA;
-            iparams.harmonic.rB  = rB;
-            iparams.harmonic.krA = krA;
-            iparams.harmonic.krB = krB;
-            ftype                = ft;
-            fep                  = (rA != rB || krA != krB);
-            return *this;
-        }
-        /*! \brief Set parameters for harmonic potential
-         *
-         * \param[in] ft  Function type
-         * \param[in] rA  Equilibrium value
-         * \param[in] krA Force constant
-         * \return The structure itself.
-         */
-        iListInput setHarmonic(int ft, real rA, real krA)
-        {
-            return setHarmonic(ft, rA, krA, rA, krA);
-        }
-        /*! \brief Set parameters for cubic potential
-         *
-         * \param[in] b0   Equilibrium bond length
-         * \param[in] kb   Harmonic force constant
-         * \param[in] kcub Cubic force constant
-         * \return The structure itself.
-         */
-        iListInput setCubic(real b0, real kb, real kcub)
-        {
-            ftype              = F_CUBICBONDS;
-            iparams.cubic.b0   = b0;
-            iparams.cubic.kb   = kb;
-            iparams.cubic.kcub = kcub;
-            return *this;
-        }
-        /*! \brief Set parameters for morse potential
-         *
-         * Free energy perturbation is turned on when A
-         * and B parameters are different.
-         * \param[in] b0A   Equilibrium value A
-         * \param[in] cbA   Force constant A
-         * \param[in] betaA Steepness parameter A
-         * \param[in] b0B   Equilibrium value B
-         * \param[in] cbB   Force constant B
-         * \param[in] betaB Steepness parameter B
-         * \return The structure itself.
-         */
-        iListInput setMorse(real b0A, real cbA, real betaA, real b0B, real cbB, real betaB)
-        {
-            ftype               = F_MORSE;
-            iparams.morse.b0A   = b0A;
-            iparams.morse.cbA   = cbA;
-            iparams.morse.betaA = betaA;
-            iparams.morse.b0B   = b0B;
-            iparams.morse.cbB   = cbB;
-            iparams.morse.betaB = betaB;
-            fep                 = (b0A != b0B || cbA != cbB || betaA != betaB);
-            return *this;
-        }
-        /*! \brief Set parameters for morse potential
-         *
-         * \param[in] b0A   Equilibrium value
-         * \param[in] cbA   Force constant
-         * \param[in] betaA Steepness parameter
-         * \return The structure itself.
-         */
-        iListInput setMorse(real b0A, real cbA, real betaA)
-        {
-            return setMorse(b0A, cbA, betaA, b0A, cbA, betaA);
-        }
-        /*! \brief Set parameters for fene potential
-         *
-         * \param[in] bm Equilibrium bond length
-         * \param[in] kb Force constant
-         * \return The structure itself.
-         */
-        iListInput setFene(real bm, real kb)
-        {
-            ftype           = F_FENEBONDS;
-            iparams.fene.bm = bm;
-            iparams.fene.kb = kb;
-            return *this;
-        }
-        /*! \brief Set parameters for linear angle potential
-         *
-         * Free energy perturbation is turned on when A
-         * and B parameters are different.
-         * \param[in] klinA Force constant A
-         * \param[in] aA    The position of the central atom A
-         * \param[in] klinB Force constant B
-         * \param[in] aB    The position of the central atom B
-         * \return The structure itself.
-         */
-        iListInput setLinearAngle(real klinA, real aA, real klinB, real aB)
-        {
-            ftype                  = F_LINEAR_ANGLES;
-            iparams.linangle.klinA = klinA;
-            iparams.linangle.aA    = aA;
-            iparams.linangle.klinB = klinB;
-            iparams.linangle.aB    = aB;
-            fep                    = (klinA != klinB || aA != aB);
-            return *this;
-        }
-        /*! \brief Set parameters for linear angle potential
-         *
-         * \param[in] klinA Force constant
-         * \param[in] aA    The position of the central atom
-         * \return The structure itself.
-         */
-        iListInput setLinearAngle(real klinA, real aA)
-        {
-            return setLinearAngle(klinA, aA, klinA, aA);
-        }
-        /*! \brief Set parameters for Urey Bradley potential
-         *
-         * Free energy perturbation is turned on when A
-         * and B parameters are different.
-         * \param[in] thetaA  Equilibrium angle A
-         * \param[in] kthetaA Force constant A
-         * \param[in] r13A    The distance between i and k atoms A
-         * \param[in] kUBA    The force constant for 1-3 distance A
-         * \param[in] thetaB  Equilibrium angle B
-         * \param[in] kthetaB Force constant B
-         * \param[in] r13B    The distance between i and k atoms B
-         * \param[in] kUBB    The force constant for 1-3 distance B
-         * \return The structure itself.
-         */
-        iListInput setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA,
-                                  real thetaB, real kthetaB, real r13B, real kUBB)
-        {
-            ftype               = F_UREY_BRADLEY;
-            iparams.u_b.thetaA  = thetaA;
-            iparams.u_b.kthetaA = kthetaA;
-            iparams.u_b.r13A    = r13A;
-            iparams.u_b.kUBA    = kUBA;
-            iparams.u_b.thetaB  =  thetaB;
-            iparams.u_b.kthetaB = kthetaB;
-            iparams.u_b.r13B    = r13B;
-            iparams.u_b.kUBB    = kUBB;
-            fep                 = (thetaA != thetaB || kthetaA != kthetaB || r13A != r13B || kUBA != kUBB);
-            return *this;
-        }
-        /*! \brief Set parameters for Urey Bradley potential
-         *
-         * \param[in] thetaA  Equilibrium angle
-         * \param[in] kthetaA Force constant
-         * \param[in] r13A    The distance between i and k atoms
-         * \param[in] kUBA    The force constant for 1-3 distance
-         * \return The structure itself.
-         */
-        iListInput setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA)
-        {
-            return setUreyBradley(thetaA, kthetaA, r13A,  kUBA,
-                                  thetaA, kthetaA, r13A,  kUBA);
-        }
-        /*! \brief Set parameters for Cross Bond Bonds potential
-         *
-         * \param[in] r1e  First bond length i-j
-         * \param[in] r2e  Second bond length i-k
-         * \param[in] krr  The force constant
-         * \return The structure itself.
-         */
-        iListInput setCrossBondBonds(real r1e, real r2e, real krr)
-        {
-            ftype                = F_CROSS_BOND_BONDS;
-            iparams.cross_bb.r1e = r1e;
-            iparams.cross_bb.r2e = r2e;
-            iparams.cross_bb.krr = krr;
-            return *this;
-        }
-        /*! \brief Set parameters for Cross Bond Angles potential
-         *
-         * \param[in] r1e  First bond length i-j
-         * \param[in] r2e  Second bond length j-k
-         * \param[in] r3e  Third bond length i-k
-         * \param[in] krt  The force constant
-         * \return The structure itself.
-         */
-        iListInput setCrossBondAngles(real r1e, real r2e, real r3e, real krt)
-        {
-            ftype                = F_CROSS_BOND_ANGLES;
-            iparams.cross_ba.r1e = r1e;
-            iparams.cross_ba.r2e = r2e;
-            iparams.cross_ba.r3e = r3e;
-            iparams.cross_ba.krt = krt;
-            return *this;
-        }
-        /*! \brief Set parameters for Quartic Angles potential
-         *
-         * \param[in] theta Angle
-         * \param[in] c     Array of parameters
-         * \return The structure itself.
-         */
-        iListInput setQuarticAngles(real theta, const real c[5])
-        {
-            ftype                = F_QUARTIC_ANGLES;
-            iparams.qangle.theta = theta;
-            iparams.qangle.c[0]  = c[0];
-            iparams.qangle.c[1]  = c[1];
-            iparams.qangle.c[2]  = c[2];
-            iparams.qangle.c[3]  = c[3];
-            iparams.qangle.c[4]  = c[4];
-            return *this;
-        }
-        /*! \brief Set parameters for proper dihedrals potential
-         *
-         * Free energy perturbation is turned on when A
-         * and B parameters are different.
-         * \param[in] phiA Dihedral angle A
-         * \param[in] cpA  Force constant A
-         * \param[in] mult Multiplicity of the angle
-         * \param[in] phiB Dihedral angle B
-         * \param[in] cpB  Force constant B
-         * \return The structure itself.
-         */
-        iListInput setProperDihedrals(real phiA, real cpA, int mult, real phiB, real cpB)
-        {
-            ftype              = F_PDIHS;
-            iparams.pdihs.phiA = phiA;
-            iparams.pdihs.cpA  = cpA;
-            iparams.pdihs.phiB = phiB;
-            iparams.pdihs.cpB  = cpB;
-            iparams.pdihs.mult = mult;
-            fep                = (phiA != phiB || cpA != cpB);
-            return *this;
-        }
-        /*! \brief Set parameters for proper dihedrals potential
-         *
-         * \param[in] phiA Dihedral angle
-         * \param[in] cpA  Force constant
-         * \param[in] mult Multiplicity of the angle
-         * \return The structure itself.
-         */
-        iListInput setProperDihedrals(real phiA, real cpA, int mult)
-        {
-            return setProperDihedrals(phiA, cpA, mult, phiA, cpA);
-        }
-        /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
-         *
-         * Free energy perturbation is turned on when A
-         * and B parameters are different.
-         * \param[in] rbcA Force constants A
-         * \param[in] rbcB Force constants B
-         * \return The structure itself.
-         */
-        iListInput setRbDihedrals(const real rbcA[NR_RBDIHS], const real rbcB[NR_RBDIHS])
-        {
-            ftype = F_RBDIHS;
-            fep   = false;
-            for (int i = 0; i < NR_RBDIHS; i++)
-            {
-                iparams.rbdihs.rbcA[i] = rbcA[i];
-                iparams.rbdihs.rbcB[i] = rbcB[i];
-                fep                    = fep || (rbcA[i] != rbcB[i]);
-            }
-            return *this;
-        }
-        /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
-         *
-         * \param[in] rbc Force constants
-         * \return The structure itself.
-         */
-        iListInput setRbDihedrals(const real rbc[NR_RBDIHS])
+public:
+    //! Function type
+    int ftype = -1;
+    //! Tolerance for float evaluation
+    float ftoler = 1e-6;
+    //! Tolerance for double evaluation
+    double dtoler = 1e-8;
+    //! Do free energy perturbation?
+    bool fep = false;
+    //! Interaction parameters
+    t_iparams iparams = { { 0 } };
+
+    friend std::ostream& operator<<(std::ostream& out, const iListInput& input);
+
+    //! Constructor
+    iListInput() {}
+
+    /*! \brief Constructor with tolerance
+     *
+     * \param[in] ftol Single precision tolerance
+     * \param[in] dtol Double precision tolerance
+     */
+    iListInput(float ftol, double dtol)
+    {
+        ftoler = ftol;
+        dtoler = dtol;
+    }
+    /*! \brief Set parameters for harmonic potential
+     *
+     * Free energy perturbation is turned on when A
+     * and B parameters are different.
+     * \param[in] ft  Function type
+     * \param[in] rA  Equilibrium value A
+     * \param[in] krA Force constant A
+     * \param[in] rB  Equilibrium value B
+     * \param[in] krB Force constant B
+     * \return The structure itself.
+     */
+    iListInput setHarmonic(int ft, real rA, real krA, real rB, real krB)
+    {
+        iparams.harmonic.rA  = rA;
+        iparams.harmonic.rB  = rB;
+        iparams.harmonic.krA = krA;
+        iparams.harmonic.krB = krB;
+        ftype                = ft;
+        fep                  = (rA != rB || krA != krB);
+        return *this;
+    }
+    /*! \brief Set parameters for harmonic potential
+     *
+     * \param[in] ft  Function type
+     * \param[in] rA  Equilibrium value
+     * \param[in] krA Force constant
+     * \return The structure itself.
+     */
+    iListInput setHarmonic(int ft, real rA, real krA) { return setHarmonic(ft, rA, krA, rA, krA); }
+    /*! \brief Set parameters for cubic potential
+     *
+     * \param[in] b0   Equilibrium bond length
+     * \param[in] kb   Harmonic force constant
+     * \param[in] kcub Cubic force constant
+     * \return The structure itself.
+     */
+    iListInput setCubic(real b0, real kb, real kcub)
+    {
+        ftype              = F_CUBICBONDS;
+        iparams.cubic.b0   = b0;
+        iparams.cubic.kb   = kb;
+        iparams.cubic.kcub = kcub;
+        return *this;
+    }
+    /*! \brief Set parameters for morse potential
+     *
+     * Free energy perturbation is turned on when A
+     * and B parameters are different.
+     * \param[in] b0A   Equilibrium value A
+     * \param[in] cbA   Force constant A
+     * \param[in] betaA Steepness parameter A
+     * \param[in] b0B   Equilibrium value B
+     * \param[in] cbB   Force constant B
+     * \param[in] betaB Steepness parameter B
+     * \return The structure itself.
+     */
+    iListInput setMorse(real b0A, real cbA, real betaA, real b0B, real cbB, real betaB)
+    {
+        ftype               = F_MORSE;
+        iparams.morse.b0A   = b0A;
+        iparams.morse.cbA   = cbA;
+        iparams.morse.betaA = betaA;
+        iparams.morse.b0B   = b0B;
+        iparams.morse.cbB   = cbB;
+        iparams.morse.betaB = betaB;
+        fep                 = (b0A != b0B || cbA != cbB || betaA != betaB);
+        return *this;
+    }
+    /*! \brief Set parameters for morse potential
+     *
+     * \param[in] b0A   Equilibrium value
+     * \param[in] cbA   Force constant
+     * \param[in] betaA Steepness parameter
+     * \return The structure itself.
+     */
+    iListInput setMorse(real b0A, real cbA, real betaA)
+    {
+        return setMorse(b0A, cbA, betaA, b0A, cbA, betaA);
+    }
+    /*! \brief Set parameters for fene potential
+     *
+     * \param[in] bm Equilibrium bond length
+     * \param[in] kb Force constant
+     * \return The structure itself.
+     */
+    iListInput setFene(real bm, real kb)
+    {
+        ftype           = F_FENEBONDS;
+        iparams.fene.bm = bm;
+        iparams.fene.kb = kb;
+        return *this;
+    }
+    /*! \brief Set parameters for linear angle potential
+     *
+     * Free energy perturbation is turned on when A
+     * and B parameters are different.
+     * \param[in] klinA Force constant A
+     * \param[in] aA    The position of the central atom A
+     * \param[in] klinB Force constant B
+     * \param[in] aB    The position of the central atom B
+     * \return The structure itself.
+     */
+    iListInput setLinearAngle(real klinA, real aA, real klinB, real aB)
+    {
+        ftype                  = F_LINEAR_ANGLES;
+        iparams.linangle.klinA = klinA;
+        iparams.linangle.aA    = aA;
+        iparams.linangle.klinB = klinB;
+        iparams.linangle.aB    = aB;
+        fep                    = (klinA != klinB || aA != aB);
+        return *this;
+    }
+    /*! \brief Set parameters for linear angle potential
+     *
+     * \param[in] klinA Force constant
+     * \param[in] aA    The position of the central atom
+     * \return The structure itself.
+     */
+    iListInput setLinearAngle(real klinA, real aA) { return setLinearAngle(klinA, aA, klinA, aA); }
+    /*! \brief Set parameters for Urey Bradley potential
+     *
+     * Free energy perturbation is turned on when A
+     * and B parameters are different.
+     * \param[in] thetaA  Equilibrium angle A
+     * \param[in] kthetaA Force constant A
+     * \param[in] r13A    The distance between i and k atoms A
+     * \param[in] kUBA    The force constant for 1-3 distance A
+     * \param[in] thetaB  Equilibrium angle B
+     * \param[in] kthetaB Force constant B
+     * \param[in] r13B    The distance between i and k atoms B
+     * \param[in] kUBB    The force constant for 1-3 distance B
+     * \return The structure itself.
+     */
+    iListInput
+    setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA, real thetaB, real kthetaB, real r13B, real kUBB)
+    {
+        ftype               = F_UREY_BRADLEY;
+        iparams.u_b.thetaA  = thetaA;
+        iparams.u_b.kthetaA = kthetaA;
+        iparams.u_b.r13A    = r13A;
+        iparams.u_b.kUBA    = kUBA;
+        iparams.u_b.thetaB  = thetaB;
+        iparams.u_b.kthetaB = kthetaB;
+        iparams.u_b.r13B    = r13B;
+        iparams.u_b.kUBB    = kUBB;
+        fep = (thetaA != thetaB || kthetaA != kthetaB || r13A != r13B || kUBA != kUBB);
+        return *this;
+    }
+    /*! \brief Set parameters for Urey Bradley potential
+     *
+     * \param[in] thetaA  Equilibrium angle
+     * \param[in] kthetaA Force constant
+     * \param[in] r13A    The distance between i and k atoms
+     * \param[in] kUBA    The force constant for 1-3 distance
+     * \return The structure itself.
+     */
+    iListInput setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA)
+    {
+        return setUreyBradley(thetaA, kthetaA, r13A, kUBA, thetaA, kthetaA, r13A, kUBA);
+    }
+    /*! \brief Set parameters for Cross Bond Bonds potential
+     *
+     * \param[in] r1e  First bond length i-j
+     * \param[in] r2e  Second bond length i-k
+     * \param[in] krr  The force constant
+     * \return The structure itself.
+     */
+    iListInput setCrossBondBonds(real r1e, real r2e, real krr)
+    {
+        ftype                = F_CROSS_BOND_BONDS;
+        iparams.cross_bb.r1e = r1e;
+        iparams.cross_bb.r2e = r2e;
+        iparams.cross_bb.krr = krr;
+        return *this;
+    }
+    /*! \brief Set parameters for Cross Bond Angles potential
+     *
+     * \param[in] r1e  First bond length i-j
+     * \param[in] r2e  Second bond length j-k
+     * \param[in] r3e  Third bond length i-k
+     * \param[in] krt  The force constant
+     * \return The structure itself.
+     */
+    iListInput setCrossBondAngles(real r1e, real r2e, real r3e, real krt)
+    {
+        ftype                = F_CROSS_BOND_ANGLES;
+        iparams.cross_ba.r1e = r1e;
+        iparams.cross_ba.r2e = r2e;
+        iparams.cross_ba.r3e = r3e;
+        iparams.cross_ba.krt = krt;
+        return *this;
+    }
+    /*! \brief Set parameters for Quartic Angles potential
+     *
+     * \param[in] theta Angle
+     * \param[in] c     Array of parameters
+     * \return The structure itself.
+     */
+    iListInput setQuarticAngles(real theta, const real c[5])
+    {
+        ftype                = F_QUARTIC_ANGLES;
+        iparams.qangle.theta = theta;
+        iparams.qangle.c[0]  = c[0];
+        iparams.qangle.c[1]  = c[1];
+        iparams.qangle.c[2]  = c[2];
+        iparams.qangle.c[3]  = c[3];
+        iparams.qangle.c[4]  = c[4];
+        return *this;
+    }
+    /*! \brief Set parameters for proper dihedrals potential
+     *
+     * Free energy perturbation is turned on when A
+     * and B parameters are different.
+     * \param[in] ft   Function type
+     * \param[in] phiA Dihedral angle A
+     * \param[in] cpA  Force constant A
+     * \param[in] mult Multiplicity of the angle
+     * \param[in] phiB Dihedral angle B
+     * \param[in] cpB  Force constant B
+     * \return The structure itself.
+     */
+    iListInput setPDihedrals(int ft, real phiA, real cpA, int mult, real phiB, real cpB)
+    {
+        ftype              = ft;
+        iparams.pdihs.phiA = phiA;
+        iparams.pdihs.cpA  = cpA;
+        iparams.pdihs.phiB = phiB;
+        iparams.pdihs.cpB  = cpB;
+        iparams.pdihs.mult = mult;
+        fep                = (phiA != phiB || cpA != cpB);
+        return *this;
+    }
+    /*! \brief Set parameters for proper dihedrals potential
+     *
+     * \param[in] ft   Function type
+     * \param[in] phiA Dihedral angle
+     * \param[in] cpA  Force constant
+     * \param[in] mult Multiplicity of the angle
+     * \return The structure itself.
+     */
+    iListInput setPDihedrals(int ft, real phiA, real cpA, int mult)
+    {
+        return setPDihedrals(ft, phiA, cpA, mult, phiA, cpA);
+    }
+    /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
+     *
+     * Free energy perturbation is turned on when A
+     * and B parameters are different.
+     * \param[in] rbcA Force constants A
+     * \param[in] rbcB Force constants B
+     * \return The structure itself.
+     */
+    iListInput setRbDihedrals(const real rbcA[NR_RBDIHS], const real rbcB[NR_RBDIHS])
+    {
+        ftype = F_RBDIHS;
+        fep   = false;
+        for (int i = 0; i < NR_RBDIHS; i++)
         {
-            return setRbDihedrals(rbc, rbc);
+            iparams.rbdihs.rbcA[i] = rbcA[i];
+            iparams.rbdihs.rbcB[i] = rbcB[i];
+            fep                    = fep || (rbcA[i] != rbcB[i]);
         }
+        return *this;
+    }
+    /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
+     *
+     * \param[in] rbc Force constants
+     * \return The structure itself.
+     */
+    iListInput setRbDihedrals(const real rbc[NR_RBDIHS]) { return setRbDihedrals(rbc, rbc); }
+    /*! \brief Set parameters for Polarization
+     *
+     * \param[in] alpha Polarizability
+     * \return The structure itself.
+     */
+    iListInput setPolarization(real alpha)
+    {
+        ftype                  = F_POLARIZATION;
+        fep                    = false;
+        iparams.polarize.alpha = alpha;
+        return *this;
+    }
+    /*! \brief Set parameters for Anharmonic Polarization
+     *
+     * \param[in] alpha Polarizability (nm^3)
+     * \param[in] drcut The cut-off distance (nm) after which the
+     *                  fourth power kicks in
+     * \param[in] khyp  The force constant for the fourth power
+     * \return The structure itself.
+     */
+    iListInput setAnharmPolarization(real alpha, real drcut, real khyp)
+    {
+        ftype                         = F_ANHARM_POL;
+        fep                           = false;
+        iparams.anharm_polarize.alpha = alpha;
+        iparams.anharm_polarize.drcut = drcut;
+        iparams.anharm_polarize.khyp  = khyp;
+        return *this;
+    }
+    /*! \brief Set parameters for Thole Polarization
+     *
+     * \param[in] a      Thole factor
+     * \param[in] alpha1 Polarizability 1 (nm^3)
+     * \param[in] alpha2 Polarizability 2 (nm^3)
+     * \return The structure itself.
+     */
+    iListInput setTholePolarization(real a, real alpha1, real alpha2)
+    {
+        ftype                = F_THOLE_POL;
+        fep                  = false;
+        iparams.thole.a      = a;
+        iparams.thole.alpha1 = alpha1;
+        iparams.thole.alpha2 = alpha2;
+        return *this;
+    }
+    /*! \brief Set parameters for Water Polarization
+     *
+     * \param[in] alpha_x Polarizability X (nm^3)
+     * \param[in] alpha_y Polarizability Y (nm^3)
+     * \param[in] alpha_z Polarizability Z (nm^3)
+     * \param[in] rOH     Oxygen-Hydrogen distance
+     * \param[in] rHH     Hydrogen-Hydrogen distance
+     * \param[in] rOD     Oxygen-Dummy distance
+     * \return The structure itself.
+     */
+    iListInput setWaterPolarization(real alpha_x, real alpha_y, real alpha_z, real rOH, real rHH, real rOD)
+    {
+        ftype             = F_WATER_POL;
+        fep               = false;
+        iparams.wpol.al_x = alpha_x;
+        iparams.wpol.al_y = alpha_y;
+        iparams.wpol.al_z = alpha_z;
+        iparams.wpol.rOH  = rOH;
+        iparams.wpol.rHH  = rHH;
+        iparams.wpol.rOD  = rOD;
+        return *this;
+    }
 };
 
+//! Prints the interaction and parameters to a stream
+std::ostream& operator<<(std::ostream& out, const iListInput& input)
+{
+    using std::endl;
+    out << "Function type " << input.ftype << " called " << interaction_function[input.ftype].name
+        << " ie. labelled '" << interaction_function[input.ftype].longname << "' in an energy file"
+        << endl;
+
+    // Organize to print the legacy C union t_iparams, whose
+    // relevant contents vary with ftype.
+    StringOutputStream stream;
+    {
+        TextWriter writer(&stream);
+        printInteractionParameters(&writer, input.ftype, input.iparams);
+    }
+    out << "Function parameters " << stream.toString();
+    out << "Parameters trigger FEP? " << (input.fep ? "true" : "false") << endl;
+    return out;
+}
+
 /*! \brief Utility to fill iatoms struct
  *
  * \param[in]  ftype  Function type
  * \param[out] iatoms Pointer to iatoms struct
  */
-void fillIatoms(int ftype, std::vector<t_iatom> *iatoms)
+void fillIatoms(int ftype, std::vector<t_iatom>iatoms)
 {
-    std::unordered_map<int, std::vector<int> > ia =
-    { { 2, { 0, 0, 1, 0, 1, 2, 0, 2, 3 } },
-      { 3, { 0, 0, 1, 2, 0, 1, 2, 3 } },
-      { 4, { 0, 0, 1, 2, 3 } }};
+    std::unordered_map<int, std::vector<int>> ia = { { 2, { 0, 0, 1, 0, 1, 2, 0, 2, 3 } },
+                                                     { 3, { 0, 0, 1, 2, 0, 1, 2, 3 } },
+                                                     { 4, { 0, 0, 1, 2, 3 } },
+                                                     { 5, { 0, 0, 1, 2, 3, 0 } } };
     EXPECT_TRUE(ftype >= 0 && ftype < F_NRE);
     int nral = interaction_function[ftype].nratoms;
-    for (auto &i : ia[nral])
+    for (autoi : ia[nral])
     {
         iatoms->push_back(i);
     }
 }
 
-class ListedForcesTest : public ::testing::TestWithParam<std::tuple<iListInput, std::vector<gmx::RVec>, int> >
+class ListedForcesTest :
+    public ::testing::TestWithParam<std::tuple<iListInput, PaddedVector<RVec>, PbcType>>
 {
-    protected:
-        matrix                     box_;
-        t_pbc                      pbc_;
-        std::vector<gmx::RVec>     x_;
-        int                        epbc_;
-        iListInput                 input_;
-        test::TestReferenceData    refData_;
-        test::TestReferenceChecker checker_;
-        ListedForcesTest( ) :
-            checker_(refData_.rootChecker())
+protected:
+    matrix                 box_;
+    t_pbc                  pbc_;
+    PaddedVector<RVec>     x_;
+    PbcType                pbcType_;
+    iListInput             input_;
+    TestReferenceData      refData_;
+    TestReferenceChecker   checker_;
+    FloatingPointTolerance shiftForcesTolerance_ = defaultRealTolerance();
+    ListedForcesTest() : checker_(refData_.rootChecker())
+    {
+        input_   = std::get<0>(GetParam());
+        x_       = std::get<1>(GetParam());
+        pbcType_ = std::get<2>(GetParam());
+        clear_mat(box_);
+        box_[0][0] = box_[1][1] = box_[2][2] = 1.5;
+        set_pbc(&pbc_, pbcType_, box_);
+        // We need quite specific tolerances here since angle functions
+        // etc. are not very precise and reproducible.
+        test::FloatingPointTolerance tolerance(test::FloatingPointTolerance(
+                input_.ftoler, input_.dtoler, 1.0e-6, 1.0e-12, 10000, 100, false));
+        checker_.setDefaultTolerance(tolerance);
+        // The SIMD acos() is only accurate to 2-3 ULP, so the angles
+        // computed by it and the non-SIMD code paths (that use
+        // std::acos) differ by enough to require quite large
+        // tolerances for the shift forces in mixed precision.
+        float singleShiftForcesAbsoluteTolerance =
+                ((input_.ftype == F_POLARIZATION) || (input_.ftype == F_ANHARM_POL)
+                                 || (IS_ANGLE(input_.ftype))
+                         ? 5e-3
+                         : 5e-5);
+        // Note that std::numeric_limits isn't required by the standard to
+        // have an implementation for uint64_t(!) but this is likely to
+        // work because that type is likely to be a typedef for one of
+        // the other numerical types that happens to be 64-bits wide.
+        shiftForcesTolerance_ = FloatingPointTolerance(singleShiftForcesAbsoluteTolerance,
+                                                       1e-8,
+                                                       1e-6,
+                                                       1e-12,
+                                                       std::numeric_limits<uint64_t>::max(),
+                                                       std::numeric_limits<uint64_t>::max(),
+                                                       false);
+    }
+    void testOneIfunc(TestReferenceChecker* checker, const std::vector<t_iatom>& iatoms, const real lambda)
+    {
+        SCOPED_TRACE(std::string("Testing PBC type: ") + c_pbcTypeNames[pbcType_]);
+        std::vector<int>  ddgatindex = { 0, 1, 2, 3 };
+        std::vector<real> charge     = { 1.5, -2.0, 1.5, -1.0 };
+        /* Here we run both the standard, plain-C force+shift-forces+energy+free-energy
+         * kernel flavor and the potentially optimized, with SIMD and less output,
+         * force only kernels. Note that we also run the optimized kernel for free-energy
+         * input when lambda=0, as the force output should match the non free-energy case.
+         */
+        std::vector<BondedKernelFlavor> flavors = { BondedKernelFlavor::ForcesAndVirialAndEnergy };
+        if (!input_.fep || lambda == 0)
         {
-            input_ = std::get<0>(GetParam());
-            x_     = std::get<1>(GetParam());
-            epbc_  = std::get<2>(GetParam());
-            clear_mat(box_);
-            box_[0][0] = box_[1][1] = box_[2][2] = 1.5;
-            set_pbc(&pbc_, epbc_, box_);
-
-            // We need quite specific tolerances here since angle functions
-            // etc. are not very precise and reproducible.
-            test::FloatingPointTolerance tolerance(test::FloatingPointTolerance(input_.ftoler, 1.0e-6,
-                                                                                input_.dtoler, 1.0e-12,
-                                                                                10000, 100, false));
-            checker_.setDefaultTolerance(tolerance);
+            flavors.push_back(BondedKernelFlavor::ForcesSimdWhenAvailable);
         }
-        void testOneIfunc(test::TestReferenceChecker *checker,
-                          const std::vector<t_iatom> &iatoms,
-                          const real                  lambda)
+        for (const auto flavor : flavors)
         {
-            SCOPED_TRACE(std::string("Testing PBC ") + epbc_names[epbc_]);
-            std::vector<int> ddgatindex = { 0, 1, 2, 3 };
+            SCOPED_TRACE("Testing bonded kernel flavor: " + c_bondedKernelFlavorStrings[flavor]);
             OutputQuantities output;
-            output.energy = bondedFunction(input_.ftype) (iatoms.size(),
-                                                          iatoms.data(),
-                                                          &input_.iparams,
-                                                          as_rvec_array(x_.data()),
-                                                          output.f, output.fshift,
-                                                          &pbc_,
-                                                          /* const struct t_graph *g */ nullptr,
-                                                          lambda, &output.dvdlambda,
-                                                          /* const struct t_mdatoms *md */ nullptr,
-                                                          /* struct t_fcdata *fcd */ nullptr,
-                                                          ddgatindex.data());
+            output.energy = calculateSimpleBond(input_.ftype,
+                                                iatoms.size(),
+                                                iatoms.data(),
+                                                &input_.iparams,
+                                                as_rvec_array(x_.data()),
+                                                output.f,
+                                                output.fshift,
+                                                &pbc_,
+                                                lambda,
+                                                &output.dvdlambda,
+                                                charge,
+                                                /* struct t_fcdata * */ nullptr,
+                                                nullptr,
+                                                nullptr,
+                                                ddgatindex.data(),
+                                                flavor);
             // Internal consistency test of both test input
             // and bonded functions.
-            EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0)));
-            checkOutput(checker, output);
-        }
-        void testIfunc()
-        {
-            test::TestReferenceChecker thisChecker =
-                checker_.checkCompound("FunctionType",
-                                       interaction_function[input_.ftype].name).checkCompound("FEP", (input_.fep ? "Yes" : "No"));
-            std::vector<t_iatom> iatoms;
-            fillIatoms(input_.ftype, &iatoms);
-            if (input_.fep)
+            EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0))) << "dvdlambda was " << output.dvdlambda;
+            checkOutput(checker, output, flavor);
+            auto shiftForcesChecker = checker->checkCompound("Shift-Forces", "Shift-forces");
+            if (computeVirial(flavor))
             {
-                const int numLambdas = 3;
-                for (int i = 0; i < numLambdas; ++i)
-                {
-                    const real lambda       = i / (numLambdas - 1.0);
-                    auto       valueChecker = thisChecker.checkCompound("Lambda", toString(lambda));
-                    testOneIfunc(&valueChecker, iatoms, lambda);
-                }
+                shiftForcesChecker.setDefaultTolerance(shiftForcesTolerance_);
+                shiftForcesChecker.checkVector(output.fshift[c_centralShiftIndex], "Central");
             }
             else
             {
-                testOneIfunc(&thisChecker, iatoms, 0.0);
+                // Permit omitting to compare shift forces with
+                // reference data when that is useless.
+                shiftForcesChecker.disableUnusedEntriesCheck();
+            }
+        }
+    }
+    void testIfunc()
+    {
+        TestReferenceChecker thisChecker =
+                checker_.checkCompound("FunctionType", interaction_function[input_.ftype].name)
+                        .checkCompound("FEP", (input_.fep ? "Yes" : "No"));
+        std::vector<t_iatom> iatoms;
+        fillIatoms(input_.ftype, &iatoms);
+        if (input_.fep)
+        {
+            const int numLambdas = 3;
+            for (int i = 0; i < numLambdas; ++i)
+            {
+                const real lambda       = i / (numLambdas - 1.0);
+                auto       valueChecker = thisChecker.checkCompound("Lambda", toString(lambda));
+                testOneIfunc(&valueChecker, iatoms, lambda);
             }
         }
+        else
+        {
+            testOneIfunc(&thisChecker, iatoms, 0.0);
+        }
+    }
 };
 
-TEST_P (ListedForcesTest, Ifunc)
+TEST_P(ListedForcesTest, Ifunc)
 {
     testIfunc();
 }
 
 //! Function types for testing bonds. Add new terms at the end.
-std::vector<iListInput> c_InputBonds =
-{
-    { iListInput().setHarmonic(F_BONDS, 0.15, 500.0) },
-    { iListInput(2e-6f, 1e-8).setHarmonic(F_BONDS, 0.15, 500.0, 0.17, 400.0) },
-    { iListInput(1e-4f, 1e-8).setHarmonic(F_G96BONDS, 0.15, 50.0) },
+std::vector<iListInput> c_InputBonds = {
+    { iListInput(2e-6F, 1e-8).setHarmonic(F_BONDS, 0.15, 500.0) },
+    { iListInput(2e-6F, 1e-8).setHarmonic(F_BONDS, 0.15, 500.0, 0.17, 400.0) },
+    { iListInput(1e-4F, 1e-8).setHarmonic(F_G96BONDS, 0.15, 50.0) },
     { iListInput().setHarmonic(F_G96BONDS, 0.15, 50.0, 0.17, 40.0) },
     { iListInput().setCubic(0.16, 50.0, 2.0) },
-    { iListInput(2e-6f, 1e-8).setMorse(0.15, 50.0, 2.0, 0.17, 40.0, 1.6) },
-    { iListInput(2e-6f, 1e-8).setMorse(0.15, 30.0, 2.7) },
+    { iListInput(2e-6F, 1e-8).setMorse(0.15, 50.0, 2.0, 0.17, 40.0, 1.6) },
+    { iListInput(2e-6F, 1e-8).setMorse(0.15, 30.0, 2.7) },
     { iListInput().setFene(0.4, 5.0) }
 };
 
@@ -524,8 +669,7 @@ std::vector<iListInput> c_InputBonds =
 const real cQuarticAngles[5] = { 1.1, 2.3, 4.6, 7.8, 9.2 };
 
 //! Function types for testing angles. Add new terms at the end.
-std::vector<iListInput> c_InputAngles =
-{
+std::vector<iListInput> c_InputAngles = {
     { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 100.0, 50.0) },
     { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 100.15, 50.0, 95.0, 30.0) },
     { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES, 100.0, 50.0) },
@@ -533,10 +677,10 @@ std::vector<iListInput> c_InputAngles =
     { iListInput().setLinearAngle(50.0, 0.4) },
     { iListInput().setLinearAngle(50.0, 0.4, 40.0, 0.6) },
     { iListInput(2e-6, 1e-8).setCrossBondBonds(0.8, 0.7, 45.0) },
-    { iListInput(2e-6, 1e-8).setCrossBondAngles(0.8, 0.7, 0.3, 45.0) },
+    { iListInput(3e-6, 1e-8).setCrossBondAngles(0.8, 0.7, 0.3, 45.0) },
     { iListInput(2e-2, 1e-8).setUreyBradley(950.0, 46.0, 0.3, 5.0) },
     { iListInput(2e-2, 1e-8).setUreyBradley(100.0, 45.0, 0.3, 5.0, 90.0, 47.0, 0.32, 7.0) },
-    { iListInput(8e-4, 1e-8).setQuarticAngles(87.0, cQuarticAngles ) }
+    { iListInput(2e-3, 1e-8).setQuarticAngles(87.0, cQuarticAngles) }
 };
 
 //! Constants for Ryckaert-Bellemans A
@@ -549,33 +693,138 @@ const real rbcB[NR_RBDIHS] = { -6.35, 12.6, 8.1, -10.7, 0.9, 15.4 };
 const real rbc[NR_RBDIHS] = { -7.35, 13.6, 8.4, -16.7, 1.3, 12.4 };
 
 //! Function types for testing dihedrals. Add new terms at the end.
-std::vector<iListInput> c_InputDihs =
-{
-    { iListInput(1e-4, 1e-8).setProperDihedrals(-100.0, 10.0, 2, -80.0, 20.0) },
-    { iListInput(1e-4, 1e-8).setProperDihedrals(-105.0, 15.0, 2) },
+std::vector<iListInput> c_InputDihs = {
+    { iListInput(5e-4, 1e-8).setPDihedrals(F_PDIHS, -100.0, 10.0, 2, -80.0, 20.0) },
+    { iListInput(1e-4, 1e-8).setPDihedrals(F_PDIHS, -105.0, 15.0, 2) },
     { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS, 100.0, 50.0) },
     { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS, 100.15, 50.0, 95.0, 30.0) },
-    { iListInput(3e-4, 1e-8).setRbDihedrals(rbcA, rbcB) },
-    { iListInput(3e-4, 1e-8).setRbDihedrals(rbc) }
+    { iListInput(4e-4, 1e-8).setRbDihedrals(rbcA, rbcB) },
+    { iListInput(4e-4, 1e-8).setRbDihedrals(rbc) }
+};
+
+//! Function types for testing polarization. Add new terms at the end.
+std::vector<iListInput> c_InputPols = {
+    { iListInput(2e-5, 1e-8).setPolarization(0.12) },
+    { iListInput(2e-3, 1e-8).setAnharmPolarization(0.0013, 0.02, 1235.6) },
+    { iListInput(1.4e-3, 1e-8).setTholePolarization(0.26, 0.07, 0.09) },
+    { iListInput(2e-3, 1e-8).setWaterPolarization(0.001, 0.0012, 0.0016, 0.095, 0.15, 0.02) },
+};
+
+//! Function types for testing polarization. Add new terms at the end.
+std::vector<iListInput> c_InputRestraints = {
+    { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRES, -100.0, 10.0, 2, -80.0, 20.0) },
+    { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRES, -105.0, 15.0, 2) },
+    { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRESZ, -100.0, 10.0, 2, -80.0, 20.0) },
+    { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRESZ, -105.0, 15.0, 2) },
+    { iListInput(2e-3, 1e-8).setHarmonic(F_RESTRANGLES, 100.0, 50.0) },
+    { iListInput(2e-3, 1e-8).setHarmonic(F_RESTRANGLES, 100.0, 50.0, 110.0, 45.0) }
 };
 
-//! Coordinates for testing
-std::vector<std::vector<gmx::RVec> > c_coordinatesForTests =
+//! Function types for testing bond with zero length, has zero reference length to make physical sense.
+std::vector<iListInput> c_InputBondsZeroLength = {
+    { iListInput().setHarmonic(F_BONDS, 0.0, 500.0) },
+};
+
+//! Function types for testing angles with zero angle, has zero reference angle to make physical sense.
+std::vector<iListInput> c_InputAnglesZeroAngle = {
+    { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 0.0, 50.0) },
+};
+
+} // namespace
+} // namespace test
+
+//! Print an RVec to \c os
+static void PrintTo(const RVec& value, std::ostream* os)
 {
-    {{  0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.2 }, { 0.005, 0.0, 0.1 }, { -0.001, 0.1, 0.0 }},
-    {{  0.5, 0.0, 0.0 }, { 0.5, 0.0, 0.15 }, { 0.5, 0.07, 0.22 }, { 0.5, 0.18, 0.22 }},
-    {{ -0.1143, -0.0282, 0.0 }, { 0.0, 0.0434, 0.0 }, { 0.1185, -0.0138, 0.0 }, { -0.0195, 0.1498, 0.0 }}
+    *os << value[XX] << " " << value[YY] << " " << value[ZZ] << std::endl;
+}
+
+//! Print a padded vector of RVec to \c os
+static void PrintTo(const PaddedVector<RVec>& vector, std::ostream* os)
+{
+    if (vector.empty())
+    {
+        *os << "Empty vector" << std::endl;
+    }
+    else
+    {
+        *os << "Vector of RVec containing:" << std::endl;
+        std::for_each(vector.begin(), vector.end(), [os](const RVec& v) { PrintTo(v, os); });
+    }
+}
+
+namespace test
+{
+namespace
+{
+
+/*! \brief Coordinates for testing
+ *
+ * Taken from a butane molecule, so we have some
+ * normal-sized bonds and angles to test.
+ *
+ * \todo Test also some weirder values */
+std::vector<PaddedVector<RVec>> c_coordinatesForTests = {
+    { { 1.382, 1.573, 1.482 }, { 1.281, 1.559, 1.596 }, { 1.292, 1.422, 1.663 }, { 1.189, 1.407, 1.775 } }
+};
+
+//! Coordinates for testing bonds with zero length
+std::vector<PaddedVector<RVec>> c_coordinatesForTestsZeroBondLength = {
+    { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.005, 0.0, 0.1 }, { 0.005, 0.0, 0.1 } }
+};
+
+//! Coordinates for testing bonds with zero length
+std::vector<PaddedVector<RVec>> c_coordinatesForTestsZeroAngle = {
+    { { 0.005, 0.0, 0.1 }, { 0.0, 0.0, 0.0 }, { 0.005, 0.0, 0.1 }, { 0.5, 0.18, 0.22 } }
 };
 
 //! PBC values for testing
-std::vector<int> c_pbcForTests = { epbcNONE, epbcXY, epbcXYZ };
+std::vector<PbcType> c_pbcForTests = { PbcType::No, PbcType::XY, PbcType::Xyz };
+
+INSTANTIATE_TEST_SUITE_P(Bond,
+                         ListedForcesTest,
+                         ::testing::Combine(::testing::ValuesIn(c_InputBonds),
+                                            ::testing::ValuesIn(c_coordinatesForTests),
+                                            ::testing::ValuesIn(c_pbcForTests)));
+
+INSTANTIATE_TEST_SUITE_P(Angle,
+                         ListedForcesTest,
+                         ::testing::Combine(::testing::ValuesIn(c_InputAngles),
+                                            ::testing::ValuesIn(c_coordinatesForTests),
+                                            ::testing::ValuesIn(c_pbcForTests)));
+
+INSTANTIATE_TEST_SUITE_P(Dihedral,
+                         ListedForcesTest,
+                         ::testing::Combine(::testing::ValuesIn(c_InputDihs),
+                                            ::testing::ValuesIn(c_coordinatesForTests),
+                                            ::testing::ValuesIn(c_pbcForTests)));
+
+INSTANTIATE_TEST_SUITE_P(Polarize,
+                         ListedForcesTest,
+                         ::testing::Combine(::testing::ValuesIn(c_InputPols),
+                                            ::testing::ValuesIn(c_coordinatesForTests),
+                                            ::testing::ValuesIn(c_pbcForTests)));
+
+INSTANTIATE_TEST_SUITE_P(Restraints,
+                         ListedForcesTest,
+                         ::testing::Combine(::testing::ValuesIn(c_InputRestraints),
+                                            ::testing::ValuesIn(c_coordinatesForTests),
+                                            ::testing::ValuesIn(c_pbcForTests)));
 
-INSTANTIATE_TEST_CASE_P(Bond, ListedForcesTest, ::testing::Combine(::testing::ValuesIn(c_InputBonds), ::testing::ValuesIn(c_coordinatesForTests), ::testing::ValuesIn(c_pbcForTests)));
+INSTANTIATE_TEST_SUITE_P(BondZeroLength,
+                         ListedForcesTest,
+                         ::testing::Combine(::testing::ValuesIn(c_InputBondsZeroLength),
+                                            ::testing::ValuesIn(c_coordinatesForTestsZeroBondLength),
+                                            ::testing::ValuesIn(c_pbcForTests)));
 
-INSTANTIATE_TEST_CASE_P(Angle, ListedForcesTest, ::testing::Combine(::testing::ValuesIn(c_InputAngles), ::testing::ValuesIn(c_coordinatesForTests), ::testing::ValuesIn(c_pbcForTests)));
+INSTANTIATE_TEST_SUITE_P(AngleZero,
+                         ListedForcesTest,
+                         ::testing::Combine(::testing::ValuesIn(c_InputAnglesZeroAngle),
+                                            ::testing::ValuesIn(c_coordinatesForTestsZeroAngle),
+                                            ::testing::ValuesIn(c_pbcForTests)));
 
-INSTANTIATE_TEST_CASE_P(Dihedral, ListedForcesTest, ::testing::Combine(::testing::ValuesIn(c_InputDihs), ::testing::ValuesIn(c_coordinatesForTests), ::testing::ValuesIn(c_pbcForTests)));
+} // namespace
 
-}  // namespace
+} // namespace test
 
-}  // namespace gmx
+} // namespace gmx