/*
* 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/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
{
//! Derivative with respect to lambda
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");
}
//! Interaction parameters
t_iparams iparams = { { 0 } };
+ friend std::ostream& operator<<(std::ostream& out, const iListInput& input);
+
//! Constructor
iListInput() {}
* \param[in] a Thole factor
* \param[in] alpha1 Polarizability 1 (nm^3)
* \param[in] alpha2 Polarizability 2 (nm^3)
- * \param[in] rfac Distance factor
* \return The structure itself.
*/
- iListInput setTholePolarization(real a, real alpha1, real alpha2, real rfac)
+ 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;
- iparams.thole.rfac = rfac;
return *this;
}
/*! \brief Set parameters for Water Polarization
}
};
+//! 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
}
class ListedForcesTest :
- public ::testing::TestWithParam<std::tuple<iListInput, std::vector<gmx::RVec>, int>>
+ 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_;
+ 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());
- epbc_ = std::get<2>(GetParam());
+ 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_, epbc_, box_);
+ 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, 1.0e-6, input_.dtoler, 1.0e-12, 10000, 100, false));
+ 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(test::TestReferenceChecker* checker, const std::vector<t_iatom>& iatoms, const real lambda)
+ void testOneIfunc(TestReferenceChecker* checker, const std::vector<t_iatom>& iatoms, const real lambda)
{
- SCOPED_TRACE(std::string("Testing PBC ") + epbc_names[epbc_]);
+ SCOPED_TRACE(std::string("Testing PBC type: ") + c_pbcTypeNames[pbcType_]);
std::vector<int> ddgatindex = { 0, 1, 2, 3 };
- std::vector<real> chargeA = { 1.5, -2.0, 1.5, -1.0 };
- t_mdatoms mdatoms = { 0 };
- mdatoms.chargeA = chargeA.data();
- OutputQuantities output;
- output.energy = calculateSimpleBond(
- 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, &mdatoms,
- /* struct t_fcdata * */ nullptr, ddgatindex.data(),
- BondedKernelFlavor::ForcesAndVirialAndEnergy);
- // Internal consistency test of both test input
- // and bonded functions.
- EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0)));
- checkOutput(checker, output);
+ 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)
+ {
+ flavors.push_back(BondedKernelFlavor::ForcesSimdWhenAvailable);
+ }
+ for (const auto flavor : flavors)
+ {
+ SCOPED_TRACE("Testing bonded kernel flavor: " + c_bondedKernelFlavorStrings[flavor]);
+ OutputQuantities output;
+ 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))) << "dvdlambda was " << output.dvdlambda;
+ checkOutput(checker, output, flavor);
+ auto shiftForcesChecker = checker->checkCompound("Shift-Forces", "Shift-forces");
+ if (computeVirial(flavor))
+ {
+ shiftForcesChecker.setDefaultTolerance(shiftForcesTolerance_);
+ shiftForcesChecker.checkVector(output.fshift[c_centralShiftIndex], "Central");
+ }
+ else
+ {
+ // Permit omitting to compare shift forces with
+ // reference data when that is useless.
+ shiftForcesChecker.disableUnusedEntriesCheck();
+ }
+ }
}
void testIfunc()
{
- test::TestReferenceChecker thisChecker =
+ TestReferenceChecker thisChecker =
checker_.checkCompound("FunctionType", interaction_function[input_.ftype].name)
.checkCompound("FEP", (input_.fep ? "Yes" : "No"));
std::vector<t_iatom> iatoms;
//! 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) },
{ 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) },
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, 1.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) },
};
{ 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 = {
- { { 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 } }
+//! 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)
+{
+ *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 };
-
-// Those tests give errors with the intel compiler and nothing else, so we disable them only there.
-#ifndef __INTEL_COMPILER
-INSTANTIATE_TEST_CASE_P(Bond,
- ListedForcesTest,
- ::testing::Combine(::testing::ValuesIn(c_InputBonds),
- ::testing::ValuesIn(c_coordinatesForTests),
- ::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_CASE_P(Dihedral,
- ListedForcesTest,
- ::testing::Combine(::testing::ValuesIn(c_InputDihs),
- ::testing::ValuesIn(c_coordinatesForTests),
- ::testing::ValuesIn(c_pbcForTests)));
-
-INSTANTIATE_TEST_CASE_P(Polarize,
- ListedForcesTest,
- ::testing::Combine(::testing::ValuesIn(c_InputPols),
- ::testing::ValuesIn(c_coordinatesForTests),
- ::testing::ValuesIn(c_pbcForTests)));
-
-INSTANTIATE_TEST_CASE_P(Restraints,
- ListedForcesTest,
- ::testing::Combine(::testing::ValuesIn(c_InputRestraints),
- ::testing::ValuesIn(c_coordinatesForTests),
- ::testing::ValuesIn(c_pbcForTests)));
-#endif
+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_SUITE_P(BondZeroLength,
+ ListedForcesTest,
+ ::testing::Combine(::testing::ValuesIn(c_InputBondsZeroLength),
+ ::testing::ValuesIn(c_coordinatesForTestsZeroBondLength),
+ ::testing::ValuesIn(c_pbcForTests)));
+
+INSTANTIATE_TEST_SUITE_P(AngleZero,
+ ListedForcesTest,
+ ::testing::Combine(::testing::ValuesIn(c_InputAnglesZeroAngle),
+ ::testing::ValuesIn(c_coordinatesForTestsZeroAngle),
+ ::testing::ValuesIn(c_pbcForTests)));
+
} // namespace
+} // namespace test
+
} // namespace gmx
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff