src/gromacs/listed_forces/tests/pairs.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, 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 Implements test of 1-4 interactions
  37  *
  38  * This test is copied from the bonded interactions test and slightly
  39  * modified since 'do_pairs' takes a different set of arguments than
  40  * 'calculateSimpleBond'. To keep the test setup uncluttered this test is
  41  * therefore not merged into the bonded test but implemented standalone.
  42  *
  43  * The test setup consists of 2 atom pairs that are tested in an fep setting
  44  * (vanishing charge and lennard-jones parameters of one atom) and without
  45  * fep. Placement of the atoms in the box is such that shift-forces and pbc
  46  * paths in do_pairs are covered.
  47  *
  48  * \author Sebastian Kehl <sebastian.kehl@mpcdf.mpg.de>
  49  * \ingroup module_listed_forces
  50  */
  51 #include "gmxpre.h"
  52
  53 #include "gromacs/listed_forces/bonded.h"
  54
  55 #include <cmath>
  56
  57 #include <memory>
  58 #include <unordered_map>
  59
  60 #include <gtest/gtest.h>
  61
  62 #include "gromacs/listed_forces/listed_forces.h"
  63 #include "gromacs/listed_forces/pairs.h"
  64 #include "gromacs/math/paddedvector.h"
  65 #include "gromacs/math/units.h"
  66 #include "gromacs/math/vec.h"
  67 #include "gromacs/math/vectypes.h"
  68 #include "gromacs/mdtypes/mdatom.h"
  69 #include "gromacs/mdtypes/simulation_workload.h"
  70 #include "gromacs/mdtypes/enerdata.h"
  71 #include "gromacs/mdtypes/forcerec.h"
  72 #include "gromacs/mdtypes/inputrec.h"
  73 #include "gromacs/mdtypes/interaction_const.h"
  74 #include "gromacs/mdtypes/nblist.h"
  75 #include "gromacs/tables/forcetable.h"
  76 #include "gromacs/pbcutil/ishift.h"
  77 #include "gromacs/pbcutil/pbc.h"
  78 #include "gromacs/topology/idef.h"
  79 #include "gromacs/utility/enumerationhelpers.h"
  80 #include "gromacs/utility/strconvert.h"
  81 #include "gromacs/utility/stringstream.h"
  82 #include "gromacs/utility/textwriter.h"
  83 #include "gromacs/utility/fatalerror.h"
  84
  85 #include "testutils/refdata.h"
  86 #include "testutils/testasserts.h"
  87
  88 namespace gmx
  89 {
  90 namespace test
  91 {
  92 namespace
  93 {
  94
  95 //! Number of atoms used in these tests.
  96 constexpr int c_numAtoms = 3;
  97
  98 /*! \brief Output from pairs kernels
  99  *
 100  */
 101 struct OutputQuantities
 102 {
 103     OutputQuantities(int energyGroup) :
 104         energy(energyGroup), dvdLambda(static_cast<int>(FreeEnergyPerturbationCouplingType::Count), 0.0)
 105     {
 106     }
 107
 108     //! Energy of this interaction
 109     gmx_grppairener_t energy;
 110     //! Derivative with respect to lambda
 111     std::vector<real> dvdLambda;
 112     //! Shift vectors
 113     rvec fShift[gmx::detail::c_numIvecs] = { { 0 } };
 114     //! Forces
 115     alignas(GMX_REAL_MAX_SIMD_WIDTH * sizeof(real)) rvec4 f[c_numAtoms] = { { 0 } };
 116 };
 117
 118 /*! \brief Utility to check the output from pairs tests
 119  *
 120  * \param[in] checker Reference checker
 121  * \param[in] output  The output from the test to check
 122  * \param[in] bondedKernelFlavor  Flavor for determining what output to check
 123  * \param[in] functionType type of the interaction
 124  */
 125 void checkOutput(TestReferenceChecker*    checker,
 126                  const OutputQuantities&  output,
 127                  const BondedKernelFlavor bondedKernelFlavor,
 128                  const int                functionType)
 129 {
 130     if (computeEnergy(bondedKernelFlavor))
 131     {
 132         switch (functionType)
 133         {
 134             case F_LJ14:
 135             case F_LJC14_Q:
 136                 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14][0],
 137                                    "Epot Coulomb14");
 138                 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJ14][0],
 139                                    "Epot LJ14");
 140                 break;
 141             case F_LJC_PAIRS_NB:
 142                 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR][0],
 143                                    "Epot Coulomb14");
 144                 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJSR][0],
 145                                    "Epot LJ14");
 146                 break;
 147             default: gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", functionType);
 148         }
 149         checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)],
 150                            "dVdlCoul ");
 151         checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)],
 152                            "dVdlVdw ");
 153     }
 154     checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
 155 }
 156
 157 /* \brief Utility class to setup forcerec and interaction parameters
 158  *
 159  * Data is only initialized as necessary for the 1-4 interactions to work!
 160  */
 161 class ForcerecHelper
 162 {
 163 public:
 164     ForcerecHelper()
 165     {
 166         fepVals_.sc_alpha                = 0.3;
 167         fepVals_.sc_power                = 1;
 168         fepVals_.sc_r_power              = 6.0;
 169         fepVals_.sc_sigma                = 0.3;
 170         fepVals_.sc_sigma_min            = 0.3;
 171         fepVals_.bScCoul                 = true;
 172         fepVals_.scGapsysScaleLinpointLJ = 0.85;
 173         fepVals_.scGapsysScaleLinpointQ  = 0.3;
 174         fepVals_.scGapsysSigmaLJ         = 0.3;
 175         fepVals_.softcoreFunction        = SoftcoreType::Beutler;
 176
 177         fr_.ic = std::make_unique<interaction_const_t>();
 178         // set data in ic
 179         fr_.ic->softCoreParameters = std::make_unique<interaction_const_t::SoftCoreParameters>(fepVals_);
 180
 181         // set data in fr
 182         real tableRange = 2.9;
 183         fr_.pairsTable = make_tables(nullptr, fr_.ic.get(), nullptr, tableRange, GMX_MAKETABLES_14ONLY);
 184         fr_.efep       = haveFep_;
 185         fr_.fudgeQQ          = 0.5;
 186         fr_.use_simd_kernels = useSimd_;
 187         fr_.bMolPBC          = haveMolPBC_;
 188     }
 189
 190     t_forcerec* get() { return &fr_; }
 191
 192 private:
 193     FreeEnergyPerturbationType haveFep_    = FreeEnergyPerturbationType::No;
 194     bool                       useSimd_    = false;
 195     bool                       haveMolPBC_ = false;
 196
 197     t_lambda   fepVals_;
 198     t_forcerec fr_;
 199 };
 200
 201 ForcerecHelper frHelper;
 202
 203
 204 /*! \brief Input structure for listed forces tests
 205  */
 206 struct ListInput
 207 {
 208 public:
 209     //! Function type
 210     int fType = -1;
 211     //! do fep
 212     bool fep = false;
 213     //! Tolerance for float evaluation
 214     float floatToler = 1e-6;
 215     //! Tolerance for double evaluation
 216     double doubleToler = 1e-8;
 217     //! Interaction parameters
 218     t_iparams iparams = { { 0 } };
 219
 220     //! Constructor
 221     ListInput() {}
 222
 223     /*! \brief Constructor with tolerance
 224      *
 225      * \param[in] ftol Single precision tolerance
 226      * \param[in] dtol Double precision tolerance
 227      */
 228     ListInput(float ftol, double dtol)
 229     {
 230         floatToler  = ftol;
 231         doubleToler = dtol;
 232     }
 233
 234     /*! \brief Set parameters for lj14 interaction
 235      *
 236      * Fep is used if either c6A != c6B or c12A != c12B.
 237      *
 238      * \param[in] c6A  lj-c6 of state A
 239      * \param[in] c12A lj-c12 of state A
 240      * \param[in] c6B  lj-c6 of state B
 241      * \param[in] c12B lj-c12 of state B
 242      */
 243     ListInput setLj14Interaction(real c6A, real c12A, real c6B, real c12B)
 244     {
 245         fType             = F_LJ14;
 246         fep               = (c6A != c6B || c12A != c12B);
 247         iparams.lj14.c6A  = c6A;
 248         iparams.lj14.c12A = c12A;
 249         iparams.lj14.c6B  = c6B;
 250         iparams.lj14.c12B = c12B;
 251
 252         return *this;
 253     }
 254
 255     /*! \brief Set parameters for ljc14 interaction
 256      *
 257      * \param[in] qi     charge i
 258      * \param[in] qj     charge j
 259      * \param[in] c6     lj-c6
 260      * \param[in] c12    lj-c12
 261      * \param[in] fudgeQ fudge factor
 262      */
 263     ListInput setLjc14Interaction(real qi, real qj, real c6, real c12, real fudgeQ)
 264     {
 265         fType             = F_LJC14_Q;
 266         iparams.ljc14.qi  = qi;
 267         iparams.ljc14.qj  = qj;
 268         iparams.ljc14.c6  = c6;
 269         iparams.ljc14.c12 = c12;
 270         iparams.ljc14.fqq = fudgeQ;
 271
 272         return *this;
 273     }
 274
 275     /*! \brief Set parameters for ljcnb interaction
 276      *
 277      * \param[in] qi  charge i
 278      * \param[in] qj  charge j
 279      * \param[in] c6  lj-c6
 280      * \param[in] c12 lj-c12
 281      */
 282     ListInput setLjcnbInteraction(real qi, real qj, real c6, real c12)
 283     {
 284         fType             = F_LJC_PAIRS_NB;
 285         iparams.ljcnb.qi  = qi;
 286         iparams.ljcnb.qj  = qj;
 287         iparams.ljcnb.c6  = c6;
 288         iparams.ljcnb.c12 = c12;
 289
 290         return *this;
 291     }
 292 };
 293
 294 class ListedForcesPairsTest :
 295     public ::testing::TestWithParam<std::tuple<ListInput, PaddedVector<RVec>, PbcType>>
 296 {
 297 protected:
 298     matrix               box_;
 299     t_pbc                pbc_;
 300     PaddedVector<RVec>   x_;
 301     PbcType              pbcType_;
 302     ListInput            input_;
 303     TestReferenceData    refData_;
 304     TestReferenceChecker checker_;
 305
 306     ListedForcesPairsTest() : checker_(refData_.rootChecker())
 307     {
 308         input_   = std::get<0>(GetParam());
 309         x_       = std::get<1>(GetParam());
 310         pbcType_ = std::get<2>(GetParam());
 311         clear_mat(box_);
 312         box_[0][0] = box_[1][1] = box_[2][2] = 1.0;
 313         set_pbc(&pbc_, pbcType_, box_);
 314
 315         FloatingPointTolerance tolerance = relativeToleranceAsPrecisionDependentFloatingPoint(
 316                 1.0, input_.floatToler, input_.doubleToler);
 317
 318         checker_.setDefaultTolerance(tolerance);
 319     }
 320
 321     void testOneIfunc(TestReferenceChecker* checker, const real lambda, const SoftcoreType softcoreType)
 322     {
 323         SCOPED_TRACE(std::string("Testing PBC type: ") + c_pbcTypeNames[pbcType_]);
 324
 325         // 'definition of pairs' is a concatenation of #npairs (here 2)
 326         // 'nAtomsPerPair+1'-tuples (fType a_0 a_i ... a_nAtomsPerPair)
 327         std::vector<t_iatom> iatoms = { 0, 1, 2, 0, 0, 2 };
 328
 329         std::vector<int>            ddgatindex = { 0, 1, 2 };
 330         std::vector<real>           chargeA    = { 1.0, -0.5, -0.5 };
 331         std::vector<real>           chargeB    = { 0.0, 0.0, 0.0 };
 332         std::vector<unsigned short> egrp       = { 0, 0, 0 };
 333         t_mdatoms                   mdatoms    = { 0 };
 334
 335         mdatoms.chargeA = chargeA.data();
 336         mdatoms.chargeB = chargeB.data();
 337         mdatoms.cENER   = egrp.data();
 338         // nPerturbed is not decisive for fep to be used; it is overruled by
 339         // other conditions in do_pairs_general; just here to not segfault
 340         // upon query
 341         mdatoms.nPerturbed = 0;
 342
 343         t_forcerec* fr = frHelper.get();
 344         fr->efep = input_.fep ? FreeEnergyPerturbationType::Yes : FreeEnergyPerturbationType::No;
 345         fr->ic->softCoreParameters->softcoreType = softcoreType;
 346         if (pbcType_ != PbcType::No)
 347         {
 348             fr->bMolPBC = true;
 349         }
 350
 351         std::vector<BondedKernelFlavor> flavors = { BondedKernelFlavor::ForcesAndVirialAndEnergy };
 352
 353         if (!input_.fep || lambda == 0)
 354         {
 355             fr->use_simd_kernels = true;
 356             flavors.push_back(BondedKernelFlavor::ForcesSimdWhenAvailable);
 357         }
 358
 359         for (const auto flavor : flavors)
 360         {
 361             SCOPED_TRACE("Testing bonded kernel flavor: " + c_bondedKernelFlavorStrings[flavor]);
 362
 363             StepWorkload stepWork;
 364             stepWork.computeVirial = computeVirial(flavor);
 365             stepWork.computeEnergy = computeEnergy(flavor);
 366
 367             bool havePerturbedInteractions = input_.fep;
 368             if (flavor == BondedKernelFlavor::ForcesSimdWhenAvailable)
 369             {
 370                 havePerturbedInteractions = false;
 371             }
 372
 373             int numEnergyTerms      = static_cast<int>(NonBondedEnergyTerms::Count);
 374             int numFepCouplingTerms = static_cast<int>(FreeEnergyPerturbationCouplingType::Count);
 375             OutputQuantities  output(numEnergyTerms);
 376             std::vector<real> lambdas(numFepCouplingTerms, lambda);
 377
 378             do_pairs(input_.fType,
 379                      iatoms.size(),
 380                      iatoms.data(),
 381                      &input_.iparams,
 382                      as_rvec_array(x_.data()),
 383                      output.f,
 384                      output.fShift,
 385                      &pbc_,
 386                      lambdas.data(),
 387                      output.dvdLambda.data(),
 388                      gmx::arrayRefFromArray(mdatoms.chargeA, mdatoms.nr),
 389                      gmx::arrayRefFromArray(mdatoms.chargeB, mdatoms.nr),
 390                      gmx::arrayRefFromArray(mdatoms.bPerturbed, mdatoms.nr),
 391                      gmx::arrayRefFromArray(mdatoms.cENER, mdatoms.nr),
 392                      mdatoms.nPerturbed,
 393                      fr,
 394                      havePerturbedInteractions,
 395                      stepWork,
 396                      &output.energy,
 397                      ddgatindex.data());
 398
 399             checkOutput(checker, output, flavor, input_.fType);
 400             auto shiftForcesChecker = checker->checkCompound("Shift-Forces", "Shift-forces");
 401
 402             if (computeVirial(flavor))
 403             {
 404                 shiftForcesChecker.checkVector(output.fShift[gmx::c_centralShiftIndex], "Central");
 405             }
 406             else
 407             {
 408                 // Permit omitting to compare shift forces with
 409                 // reference data when that is useless.
 410                 shiftForcesChecker.disableUnusedEntriesCheck();
 411             }
 412         }
 413     }
 414
 415     void testIfunc()
 416     {
 417         TestReferenceChecker thisChecker =
 418                 checker_.checkCompound("FunctionType", interaction_function[input_.fType].name)
 419                         .checkCompound("FEP", (input_.fep ? "Yes" : "No"));
 420
 421         if (input_.fep)
 422         {
 423             const int numLambdas = 3;
 424             for (int i = 0; i < numLambdas; ++i)
 425             {
 426                 const real lambda = i / (numLambdas - 1.0);
 427                 for (SoftcoreType c : EnumerationWrapper<SoftcoreType>{})
 428                 {
 429                     auto lambdaChecker = thisChecker.checkCompound("Lambda", toString(lambda));
 430                     auto softcoreChecker = lambdaChecker.checkCompound("Sofcore", enumValueToString(c));
 431                     testOneIfunc(&softcoreChecker, lambda, c);
 432                 }
 433             }
 434         }
 435         else
 436         {
 437             testOneIfunc(&thisChecker, 0.0, SoftcoreType::Beutler);
 438         }
 439     }
 440 };
 441
 442 TEST_P(ListedForcesPairsTest, Ifunc)
 443 {
 444     testIfunc();
 445 }
 446
 447 //! Function types for testing 1-4 interaction. Add new terms at the end.
 448 std::vector<ListInput> c_14Interaction = {
 449     { ListInput(1e-5, 1e-7).setLj14Interaction(0.001458, 1.0062882e-6, 0.0, 0.0) },
 450     { ListInput(1e-5, 1e-7).setLj14Interaction(0.001458, 1.0062882e-6, 0.001458, 1.0062882e-6) },
 451     { ListInput(1e-5, 1e-7).setLjc14Interaction(1.0, -1.0, 0.001458, 1.0062882e-6, 0.5) },
 452     { ListInput(1e-5, 1e-7).setLjcnbInteraction(1.0, -1.0, 0.001458, 1.0062882e-6) }
 453 };
 454
 455 //! PBC values for testing
 456 std::vector<PbcType> c_pbcForTests = { PbcType::No, PbcType::XY, PbcType::Xyz };
 457
 458 /*! \brief Coordinates for testing 1-4 interaction
 459  *
 460  * Define coordinates for 3 atoms here, which will be used in 2 interactions.
 461  */
 462 std::vector<PaddedVector<RVec>> c_coordinatesFor14Interaction = {
 463     { { 0.0, 0.0, 0.0 }, { 1.0, 1.0, 1.0 }, { 1.1, 1.2, 1.3 } }
 464 };
 465
 466 INSTANTIATE_TEST_SUITE_P(14Interaction,
 467                          ListedForcesPairsTest,
 468                          ::testing::Combine(::testing::ValuesIn(c_14Interaction),
 469                                             ::testing::ValuesIn(c_coordinatesFor14Interaction),
 470                                             ::testing::ValuesIn(c_pbcForTests)));
 471
 472 } // namespace
 473
 474 } // namespace test
 475
 476 } // namespace gmx