84b7793c61b276a9198fdb7880f8605f0f44105a
[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / tests / nb_free_energy.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
36 /*! \internal \file
37  * \brief Implements test of nonbonded fep kernel
38  *
39  * Implements the test logic from the bonded interactions also for the
40  * nonbonded fep kernel. This requires setting up some more input
41  * structures that in the bonded case.
42  *
43  * The test setup consists of an atom pair that is evaluated in an fep setting
44  * (vanishing charge and lennard-jones parameters of atom #2) with and without
45  * softcore Potentials.
46  *
47  * \author Sebastian Kehl <sebastian.kehl@mpcdf.mpg.de>
48  * \ingroup module_gmxlib_nonbonded
49  */
50 #include "gmxpre.h"
51
52 #include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
53 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
54
55 #include <cmath>
56
57 #include <gtest/gtest.h>
58
59 #include "gromacs/math/paddedvector.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/math/vectypes.h"
63 #include "gromacs/math/arrayrefwithpadding.h"
64 #include "gromacs/mdtypes/mdatom.h"
65 #include "gromacs/mdtypes/enerdata.h"
66 #include "gromacs/mdtypes/forcerec.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/interaction_const.h"
69 #include "gromacs/mdtypes/nblist.h"
70 #include "gromacs/mdtypes/forceoutput.h"
71 #include "gromacs/mdlib/forcerec.h"
72 #include "gromacs/tables/forcetable.h"
73 #include "gromacs/pbcutil/ishift.h"
74 #include "gromacs/pbcutil/pbc.h"
75 #include "gromacs/topology/idef.h"
76 #include "gromacs/topology/forcefieldparameters.h"
77 #include "gromacs/utility/enumerationhelpers.h"
78 #include "gromacs/utility/strconvert.h"
79 #include "gromacs/utility/stringstream.h"
80 #include "gromacs/utility/textwriter.h"
81 #include "gromacs/utility/arrayref.h"
82 #include "gromacs/utility/smalloc.h"
83 #include "gromacs/ewald/ewald_utils.h"
84 #include "gromacs/gmxlib/nrnb.h"
85
86 #include "testutils/refdata.h"
87 #include "testutils/testasserts.h"
88
89 namespace gmx
90 {
91 namespace test
92 {
93 namespace
94 {
95
96 //! Number of atoms used in these tests.
97 constexpr int c_numAtoms     = 4;
98 constexpr int c_numAtomTypes = 3;
99
100 /*! \brief Output from nonbonded fep kernel
101  *
102  */
103 struct OutputQuantities
104 {
105     OutputQuantities() :
106         energy(static_cast<int>(NonBondedEnergyTerms::Count)),
107         dvdLambda(static_cast<int>(FreeEnergyPerturbationCouplingType::Count), 0.0),
108         fShift(1, { 0.0, 0.0, 0.0 }),
109         f(c_numAtoms, { 0.0, 0.0, 0.0 })
110     {
111     }
112
113     //! Energies of this interaction (size EgNR)
114     gmx_grppairener_t energy;
115     //! Derivative with respect to lambda (size efptNR)
116     std::vector<real> dvdLambda;
117     //! Shift force vectors (size N_IVEC but in this test only 1)
118     std::vector<RVec> fShift;
119     //! Forces (size c_numAtoms)
120     PaddedVector<RVec> f;
121 };
122
123 /*! \brief Utility to check the output from nonbonded test
124  *
125  * \param[in] checker Reference checker
126  * \param[in] output  The output from the test to check
127  */
128 void checkOutput(TestReferenceChecker* checker, const OutputQuantities& output)
129 {
130     checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJSR][0], "EVdw ");
131     checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR][0],
132                        "ECoul ");
133     checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)],
134                        "dVdlCoul ");
135     checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)],
136                        "dVdlVdw ");
137
138     checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
139
140     auto shiftForcesChecker = checker->checkCompound("Shift-Forces", "Shift-forces");
141     shiftForcesChecker.checkVector(output.fShift[0], "Central");
142 }
143
144 class InteractionConstHelper
145 {
146 public:
147     InteractionConstHelper() {}
148
149     //! init data to construct interaction_const
150     void initInteractionConst(CoulombInteractionType coulType, VanDerWaalsType vdwType, InteractionModifiers vdwMod)
151     {
152         coulType_ = coulType;
153         vdwType_  = vdwType;
154         vdwMod_   = vdwMod;
155
156         // initialize correction tables
157         interaction_const_t tmp;
158         tmp.ewaldcoeff_q       = calc_ewaldcoeff_q(1.0, 1.0e-5);
159         tmp.ewaldcoeff_lj      = calc_ewaldcoeff_lj(1.0, 1.0e-5);
160         tmp.eeltype            = coulType;
161         tmp.vdwtype            = vdwType;
162         tmp.coulombEwaldTables = std::make_unique<EwaldCorrectionTables>();
163         tmp.vdwEwaldTables     = std::make_unique<EwaldCorrectionTables>();
164
165         init_interaction_const_tables(nullptr, &tmp, 1.0, 0.0);
166         coulombTables_ = *tmp.coulombEwaldTables;
167         vdwTables_     = *tmp.vdwEwaldTables;
168     }
169
170     /*! \brief Setup interaction_const_t
171      *
172      * \param[in]  fepVals t_lambda struct of fep values
173      * \param[out] ic      interaction_const_t pointer with data
174      */
175     void getInteractionConst(const t_lambda& fepVals, interaction_const_t* ic)
176     {
177         ic->softCoreParameters = std::unique_ptr<interaction_const_t::SoftCoreParameters>(
178                 new interaction_const_t::SoftCoreParameters(fepVals));
179
180         ic->coulombEwaldTables  = std::unique_ptr<EwaldCorrectionTables>(new EwaldCorrectionTables);
181         *ic->coulombEwaldTables = coulombTables_;
182
183         ic->vdwEwaldTables  = std::unique_ptr<EwaldCorrectionTables>(new EwaldCorrectionTables);
184         *ic->vdwEwaldTables = vdwTables_;
185
186         // set coulomb and vdw types
187         ic->eeltype      = coulType_;
188         ic->vdwtype      = vdwType_;
189         ic->vdw_modifier = vdwMod_;
190
191         // some non default parameters used in this testcase
192         ic->epsfac                   = gmx::c_one4PiEps0 * 0.25;
193         ic->reactionFieldCoefficient = 0.0; // former k_rf
194         ic->reactionFieldShift       = 1.0; // former c_rf
195         ic->sh_ewald                 = 1.0e-5;
196         ic->sh_lj_ewald              = -1.0;
197         ic->dispersion_shift.cpot    = -1.0;
198         ic->repulsion_shift.cpot     = -1.0;
199     }
200
201 private:
202     //! correction tables
203     EwaldCorrectionTables coulombTables_;
204     EwaldCorrectionTables vdwTables_;
205
206     //! coulomb and vdw type specifiers
207     CoulombInteractionType coulType_;
208     VanDerWaalsType        vdwType_;
209     InteractionModifiers   vdwMod_;
210 };
211
212
213 /* \brief Utility class to setup forcerec
214  *
215  * This helper takes care of handling the neccessary data that are kept
216  * at various places and in various forms in the forcerec hierarchy such
217  * that this class can safely be used.
218  *
219  * Data is only initialized as necessary for the nonbonded kernel to work!
220  */
221 class ForcerecHelper
222 {
223 public:
224     ForcerecHelper()
225     {
226         fepVals_.sc_alpha     = 0.3;
227         fepVals_.sc_power     = 1;
228         fepVals_.sc_r_power   = 6.0;
229         fepVals_.sc_sigma     = 0.3;
230         fepVals_.sc_sigma_min = 0.3;
231         fepVals_.bScCoul      = true;
232     }
233
234     //! initialize data structure to construct forcerec
235     void initForcerec(const gmx_ffparams_t&  idef,
236                       CoulombInteractionType coulType,
237                       VanDerWaalsType        vdwType,
238                       InteractionModifiers   vdwMod)
239     {
240         icHelper_.initInteractionConst(coulType, vdwType, vdwMod);
241         nbfp_        = makeNonBondedParameterLists(idef.atnr, idef.iparams, false);
242         ljPmeC6Grid_ = makeLJPmeC6GridCorrectionParameters(idef.atnr, idef.iparams, LongRangeVdW::Geom);
243     }
244
245     void setSoftcoreAlpha(const real scAlpha) { fepVals_.sc_alpha = scAlpha; }
246     void setSoftcoreCoulomb(const bool scCoulomb) { fepVals_.bScCoul = scCoulomb; }
247
248     //! get forcerec data as wanted by the nonbonded kernel
249     void getForcerec(t_forcerec* fr)
250     {
251         fr->ic = std::make_unique<interaction_const_t>();
252
253         // set data in ic
254         icHelper_.getInteractionConst(fepVals_, fr->ic.get());
255
256         // set data in fr
257         fr->ljpme_c6grid = ljPmeC6Grid_;
258         fr->nbfp         = nbfp_;
259         fr->shift_vec    = { { 0.0, 0.0, 0.0 } };
260         fr->rlist        = rlist_;
261         fr->ntype        = c_numAtomTypes;
262
263         // simd
264         fr->use_simd_kernels = GMX_USE_SIMD_KERNELS;
265     }
266
267 private:
268     InteractionConstHelper icHelper_;
269     std::vector<real>      ljPmeC6Grid_;
270     std::vector<real>      nbfp_;
271     t_lambda               fepVals_;
272     real                   rlist_ = 1.0;
273 };
274
275 /*! \brief Utility structure to hold atoms data
276  *
277  * A system having 4 interactions with the following perturbation pattern:
278  *  - no perturbation
279  *  - vdw- and coulomb-perturbation
280  *  - coulomb-perturbation only
281  *  - vdw-perturbation only
282  *
283  * This is realized by defining 3 different atom types that control
284  * the vdw-perturbation. The coulomb-perturbation is controlled by directly
285  * setting the charge of the atoms at the lambda states A/B.
286  */
287 struct AtomData
288 {
289     AtomData()
290     {
291         idef.atnr = c_numAtomTypes;
292         idef.iparams.resize(2 * c_numAtomTypes * c_numAtomTypes);
293
294         // set interaction parameters for different combinations of types
295         idef.iparams[0].lj = { 0.001458, 1.0062882e-6 }; // 0-0
296         idef.iparams[1].lj = { 0.0, 0.0 };               // 0-1
297         idef.iparams[2].lj = { 0.001458, 1.0062882e-6 }; // 0-2
298         idef.iparams[3].lj = { 0.0, 0.0 };               // 1-0
299         idef.iparams[4].lj = { 0.0, 0.0 };               // 1-1
300         idef.iparams[5].lj = { 0.0, 0.0 };               // 1-2
301         idef.iparams[6].lj = { 0.001458, 1.0062882e-6 }; // 2-0
302         idef.iparams[7].lj = { 0.0, 0.0 };               // 2-1
303         idef.iparams[8].lj = { 0.001458, 1.0062882e-6 }; // 2-2
304
305         GMX_ASSERT(chargeA.size() == c_numAtoms, "This test wants 4 atoms");
306         GMX_ASSERT(chargeB.size() == c_numAtoms, "This test wants 4 atoms");
307         GMX_ASSERT(typeA.size() == c_numAtoms, "This test wants 4 atoms");
308         GMX_ASSERT(typeB.size() == c_numAtoms, "This test wants 4 atoms");
309     }
310
311     // forcefield parameters
312     gmx_ffparams_t idef;
313
314     // atom data
315     std::vector<real> chargeA = { 1.0, -1.0, -1.0, 1.0 };
316     std::vector<real> chargeB = { 1.0, 0.0, 0.0, 1.0 };
317     std::vector<int>  typeA   = { 0, 0, 0, 0 };
318     std::vector<int>  typeB   = { 0, 1, 2, 1 };
319     // perturbation pattern: {no-pert, vdw- and coul-pert, coul-pert, vdw-pert}
320
321     // neighbourhood information
322     std::vector<int> iAtoms  = { 0 };
323     std::vector<int> jAtoms  = { 0, 1, 2, 3 };
324     std::vector<int> jIndex  = { 0, 4 };
325     std::vector<int> shift   = { 0 };
326     std::vector<int> gid     = { 0 };
327     std::vector<int> exclFep = { 0, 1, 1, 1 };
328
329     // construct t_nblist
330     t_nblist getNbList()
331     {
332         t_nblist nbl;
333         nbl.nri      = 1;
334         nbl.nrj      = 4;
335         nbl.iinr     = iAtoms;
336         nbl.jindex   = jIndex;
337         nbl.jjnr     = jAtoms;
338         nbl.shift    = shift;
339         nbl.gid      = gid;
340         nbl.excl_fep = exclFep;
341         return nbl;
342     }
343 };
344
345 /*! \brief Input structure for nonbonded fep kernel
346  */
347 struct ListInput
348 {
349 public:
350     //! Function type
351     int fType = F_LJ;
352     //! Tolerance for float evaluation
353     float floatToler = 1e-6;
354     //! Tolerance for double evaluation
355     double doubleToler = 1e-8;
356     //! atom parameters
357     AtomData atoms;
358     //! forcerec helper
359     ForcerecHelper frHelper;
360
361     //! Constructor
362     ListInput() {}
363
364     /*! \brief Constructor with tolerance
365      *
366      * \param[in] ftol Single precision tolerance
367      * \param[in] dtol Double precision tolerance
368      */
369     ListInput(float ftol, double dtol)
370     {
371         floatToler  = ftol;
372         doubleToler = dtol;
373     }
374
375     /*! \brief Set parameters for nonbonded interaction
376      *
377      * \param[in] coulType coulomb type
378      * \param[in] vdwType  vdw type
379      * \param[in] vdwMod   vdw potential modifier
380      */
381     ListInput setInteraction(CoulombInteractionType coulType, VanDerWaalsType vdwType, InteractionModifiers vdwMod)
382     {
383         frHelper.initForcerec(atoms.idef, coulType, vdwType, vdwMod);
384         return *this;
385     }
386 };
387
388 class NonbondedFepTest :
389     public ::testing::TestWithParam<std::tuple<ListInput, PaddedVector<RVec>, real, real, bool>>
390 {
391 protected:
392     PaddedVector<RVec>   x_;
393     ListInput            input_;
394     real                 lambda_;
395     real                 softcoreAlpha_;
396     bool                 softcoreCoulomb_;
397     TestReferenceData    refData_;
398     TestReferenceChecker checker_;
399
400     NonbondedFepTest() : checker_(refData_.rootChecker())
401     {
402         input_           = std::get<0>(GetParam());
403         x_               = std::get<1>(GetParam());
404         lambda_          = std::get<2>(GetParam());
405         softcoreAlpha_   = std::get<3>(GetParam());
406         softcoreCoulomb_ = std::get<4>(GetParam());
407
408         test::FloatingPointTolerance tolerance(
409                 input_.floatToler, input_.doubleToler, 1.0e-6, 1.0e-12, 10000, 100, false);
410         checker_.setDefaultTolerance(tolerance);
411     }
412
413     void testKernel()
414     {
415         input_.frHelper.setSoftcoreAlpha(softcoreAlpha_);
416         input_.frHelper.setSoftcoreCoulomb(softcoreCoulomb_);
417
418         // get forcerec and interaction_const
419         t_forcerec fr;
420         input_.frHelper.getForcerec(&fr);
421
422         // t_nblist
423         t_nblist nbl = input_.atoms.getNbList();
424
425         // output buffers
426         OutputQuantities output;
427
428         // lambda vector
429         int numFepCouplingTerms = static_cast<int>(FreeEnergyPerturbationCouplingType::Count);
430         std::vector<real> lambdas(numFepCouplingTerms, lambda_);
431
432         // fep kernel data
433         int doNBFlags = 0;
434         doNBFlags |= GMX_NONBONDED_DO_FORCE;
435         doNBFlags |= GMX_NONBONDED_DO_SHIFTFORCE;
436         doNBFlags |= GMX_NONBONDED_DO_POTENTIAL;
437
438         // force buffers
439         bool                      unusedBool = true; // this bool has no effect in the kernel
440         gmx::ForceWithShiftForces forces(output.f.arrayRefWithPadding(), unusedBool, output.fShift);
441
442         // dummy counter
443         t_nrnb nrnb;
444
445         // run fep kernel
446         gmx_nb_free_energy_kernel(nbl,
447                                   x_.arrayRefWithPadding().unpaddedArrayRef(),
448                                   &forces,
449                                   fr.use_simd_kernels,
450                                   fr.ntype,
451                                   fr.rlist,
452                                   *fr.ic,
453                                   fr.shift_vec,
454                                   fr.nbfp,
455                                   fr.ljpme_c6grid,
456                                   input_.atoms.chargeA,
457                                   input_.atoms.chargeB,
458                                   input_.atoms.typeA,
459                                   input_.atoms.typeB,
460                                   doNBFlags,
461                                   lambdas,
462                                   output.dvdLambda,
463                                   output.energy.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR],
464                                   output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJSR],
465                                   &nrnb);
466
467         checkOutput(&checker_, output);
468     }
469 };
470
471 TEST_P(NonbondedFepTest, testKernel)
472 {
473     testKernel();
474 }
475
476 //! configurations to test
477 std::vector<ListInput> c_interaction = {
478     { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Cut, VanDerWaalsType::Cut, InteractionModifiers::None) },
479     { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Cut, VanDerWaalsType::Cut, InteractionModifiers::PotSwitch) },
480     { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Pme, VanDerWaalsType::Pme, InteractionModifiers::None) }
481 };
482
483 //! test parameters
484 std::vector<real> c_fepLambdas      = { 0.0, 0.5, 1.0 };
485 std::vector<real> c_softcoreAlphas  = { 0.0, 0.3 };
486 std::vector<bool> c_softcoreCoulomb = { true, false };
487
488 //! Coordinates for testing
489 std::vector<PaddedVector<RVec>> c_coordinates = {
490     { { 1.0, 1.0, 1.0 }, { 1.1, 1.15, 1.2 }, { 0.9, 0.85, 0.8 }, { 1.1, 1.15, 0.8 } }
491 };
492
493 INSTANTIATE_TEST_CASE_P(NBInteraction,
494                         NonbondedFepTest,
495                         ::testing::Combine(::testing::ValuesIn(c_interaction),
496                                            ::testing::ValuesIn(c_coordinates),
497                                            ::testing::ValuesIn(c_fepLambdas),
498                                            ::testing::ValuesIn(c_softcoreAlphas),
499                                            ::testing::ValuesIn(c_softcoreCoulomb)));
500
501 } // namespace
502
503 } // namespace test
504
505 } // namespace gmx