From a4c37232e2afb96b9c0ee1bb1375ab9c8cb958b2 Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Tue, 25 May 2021 09:41:58 +0000 Subject: [PATCH] Only Pass needed data to nonbonded param setup in forcerec --- .../gmxlib/nonbonded/tests/nb_free_energy.cpp | 5 +- src/gromacs/mdlib/forcerec.cpp | 85 +++++++++---------- src/gromacs/mdlib/forcerec.h | 21 +++-- 3 files changed, 57 insertions(+), 54 deletions(-) diff --git a/src/gromacs/gmxlib/nonbonded/tests/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/tests/nb_free_energy.cpp index 9aa1ad0659..84b7793c61 100644 --- a/src/gromacs/gmxlib/nonbonded/tests/nb_free_energy.cpp +++ b/src/gromacs/gmxlib/nonbonded/tests/nb_free_energy.cpp @@ -238,9 +238,8 @@ public: InteractionModifiers vdwMod) { icHelper_.initInteractionConst(coulType, vdwType, vdwMod); - nbfp_ = makeNonBondedParameterLists(idef, false); - t_forcerec frTmp; - ljPmeC6Grid_ = makeLJPmeC6GridCorrectionParameters(idef, frTmp); + nbfp_ = makeNonBondedParameterLists(idef.atnr, idef.iparams, false); + ljPmeC6Grid_ = makeLJPmeC6GridCorrectionParameters(idef.atnr, idef.iparams, LongRangeVdW::Geom); } void setSoftcoreAlpha(const real scAlpha) { fepVals_.sc_alpha = scAlpha; } diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index 871e7cdc33..61659e59e4 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -85,6 +85,7 @@ #include "gromacs/pbcutil/pbc.h" #include "gromacs/tables/forcetable.h" #include "gromacs/topology/mtop_util.h" +#include "gromacs/topology/idef.h" #include "gromacs/trajectory/trajectoryframe.h" #include "gromacs/utility/cstringutil.h" #include "gromacs/utility/exceptions.h" @@ -112,38 +113,38 @@ void ForceHelperBuffers::resize(int numAtoms) } } -std::vector makeNonBondedParameterLists(const gmx_ffparams_t& forceFieldParams, bool useBuckinghamPotential) +std::vector makeNonBondedParameterLists(const int numAtomTypes, + gmx::ArrayRef iparams, + bool useBuckinghamPotential) { std::vector nbfp; - int atnr; - atnr = forceFieldParams.atnr; if (useBuckinghamPotential) { - nbfp.resize(3 * atnr * atnr); + nbfp.resize(3 * numAtomTypes * numAtomTypes); int k = 0; - for (int i = 0; (i < atnr); i++) + for (int i = 0; (i < numAtomTypes); i++) { - for (int j = 0; (j < atnr); j++, k++) + for (int j = 0; (j < numAtomTypes); j++, k++) { - BHAMA(nbfp, atnr, i, j) = forceFieldParams.iparams[k].bham.a; - BHAMB(nbfp, atnr, i, j) = forceFieldParams.iparams[k].bham.b; + BHAMA(nbfp, numAtomTypes, i, j) = iparams[k].bham.a; + BHAMB(nbfp, numAtomTypes, i, j) = iparams[k].bham.b; /* nbfp now includes the 6.0 derivative prefactor */ - BHAMC(nbfp, atnr, i, j) = forceFieldParams.iparams[k].bham.c * 6.0; + BHAMC(nbfp, numAtomTypes, i, j) = iparams[k].bham.c * 6.0; } } } else { - nbfp.resize(2 * atnr * atnr); + nbfp.resize(2 * numAtomTypes * numAtomTypes); int k = 0; - for (int i = 0; (i < atnr); i++) + for (int i = 0; (i < numAtomTypes); i++) { - for (int j = 0; (j < atnr); j++, k++) + for (int j = 0; (j < numAtomTypes); j++, k++) { /* nbfp now includes the 6.0/12.0 derivative prefactors */ - C6(nbfp, atnr, i, j) = forceFieldParams.iparams[k].lj.c6 * 6.0; - C12(nbfp, atnr, i, j) = forceFieldParams.iparams[k].lj.c12 * 12.0; + C6(nbfp, numAtomTypes, i, j) = iparams[k].lj.c6 * 6.0; + C12(nbfp, numAtomTypes, i, j) = iparams[k].lj.c12 * 12.0; } } } @@ -151,41 +152,39 @@ std::vector makeNonBondedParameterLists(const gmx_ffparams_t& forceFieldPa return nbfp; } -std::vector makeLJPmeC6GridCorrectionParameters(const gmx_ffparams_t& forceFieldParams, - const t_forcerec& forceRec) +std::vector makeLJPmeC6GridCorrectionParameters(const int numAtomTypes, + gmx::ArrayRef iparams, + LongRangeVdW ljpme_combination_rule) { - int i, j, k, atnr; - real c6, c6i, c6j, c12i, c12j, epsi, epsj, sigmai, sigmaj; - /* For LJ-PME simulations, we correct the energies with the reciprocal space * inside of the cut-off. To do this the non-bonded kernels needs to have * access to the C6-values used on the reciprocal grid in pme.c */ - atnr = forceFieldParams.atnr; - std::vector grid(2 * atnr * atnr, 0.0); - for (i = k = 0; (i < atnr); i++) + std::vector grid(2 * numAtomTypes * numAtomTypes, 0.0); + int k = 0; + for (int i = 0; (i < numAtomTypes); i++) { - for (j = 0; (j < atnr); j++, k++) + for (int j = 0; (j < numAtomTypes); j++, k++) { - c6i = forceFieldParams.iparams[i * (atnr + 1)].lj.c6; - c12i = forceFieldParams.iparams[i * (atnr + 1)].lj.c12; - c6j = forceFieldParams.iparams[j * (atnr + 1)].lj.c6; - c12j = forceFieldParams.iparams[j * (atnr + 1)].lj.c12; - c6 = std::sqrt(c6i * c6j); - if (forceRec.ljpme_combination_rule == LongRangeVdW::LB && !gmx_numzero(c6) - && !gmx_numzero(c12i) && !gmx_numzero(c12j)) + real c6i = iparams[i * (numAtomTypes + 1)].lj.c6; + real c12i = iparams[i * (numAtomTypes + 1)].lj.c12; + real c6j = iparams[j * (numAtomTypes + 1)].lj.c6; + real c12j = iparams[j * (numAtomTypes + 1)].lj.c12; + real c6 = std::sqrt(c6i * c6j); + if (ljpme_combination_rule == LongRangeVdW::LB && !gmx_numzero(c6) && !gmx_numzero(c12i) + && !gmx_numzero(c12j)) { - sigmai = gmx::sixthroot(c12i / c6i); - sigmaj = gmx::sixthroot(c12j / c6j); - epsi = c6i * c6i / c12i; - epsj = c6j * c6j / c12j; - c6 = std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj)); + real sigmai = gmx::sixthroot(c12i / c6i); + real sigmaj = gmx::sixthroot(c12j / c6j); + real epsi = c6i * c6i / c12i; + real epsj = c6j * c6j / c12j; + c6 = std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj)); } /* Store the elements at the same relative positions as C6 in nbfp in order * to simplify access in the kernels */ - grid[2 * (atnr * i + j)] = c6 * 6.0; + grid[2 * (numAtomTypes * i + j)] = c6 * 6.0; } } return grid; @@ -891,14 +890,14 @@ void init_forcerec(FILE* fplog, forcerec->shift_vec.resize(gmx::c_numShiftVectors); } - if (forcerec->nbfp.empty()) + GMX_ASSERT(forcerec->nbfp.empty(), "The nonbonded force parameters should not be set up yet."); + forcerec->ntype = mtop.ffparams.atnr; + forcerec->nbfp = makeNonBondedParameterLists( + mtop.ffparams.atnr, mtop.ffparams.iparams, forcerec->haveBuckingham); + if (EVDW_PME(interactionConst->vdwtype)) { - forcerec->ntype = mtop.ffparams.atnr; - forcerec->nbfp = makeNonBondedParameterLists(mtop.ffparams, forcerec->haveBuckingham); - if (EVDW_PME(interactionConst->vdwtype)) - { - forcerec->ljpme_c6grid = makeLJPmeC6GridCorrectionParameters(mtop.ffparams, *forcerec); - } + forcerec->ljpme_c6grid = makeLJPmeC6GridCorrectionParameters( + mtop.ffparams.atnr, mtop.ffparams.iparams, forcerec->ljpme_combination_rule); } /* Copy the energy group exclusions */ diff --git a/src/gromacs/mdlib/forcerec.h b/src/gromacs/mdlib/forcerec.h index 5d8f66d390..88ab95af87 100644 --- a/src/gromacs/mdlib/forcerec.h +++ b/src/gromacs/mdlib/forcerec.h @@ -51,7 +51,8 @@ struct gmx_localtop_t; struct gmx_mtop_t; struct gmx_wallcycle; struct interaction_const_t; -struct gmx_ffparams_t; +union t_iparams; +enum class LongRangeVdW : int; namespace gmx { @@ -61,19 +62,23 @@ class PhysicalNodeCommunicator; /*! \brief Create nonbonded parameter lists * - * \param[in] forceFieldParams The forcefield parameters + * \param[in] numAtomTypes The number of atom types + * \param[in] iparams The LJ parameters * \param[in] useBuckinghamPotential Use Buckingham potential */ -std::vector makeNonBondedParameterLists(const gmx_ffparams_t& forceFieldParams, - bool useBuckinghamPotential); +std::vector makeNonBondedParameterLists(int numAtomTypes, + gmx::ArrayRef iparams, + bool useBuckinghamPotential); /*! \brief Calculate c6 parameters for grid correction * - * \param[in] forceFieldParams The forcefield parameters - * \param[in] forceRec The forcerec + * \param[in] numAtomTypes The number of atom types + * \param[in] iparams The LJ parameters + * \param[in] ljpme_combination_rule How long range LJ is treated */ -std::vector makeLJPmeC6GridCorrectionParameters(const gmx_ffparams_t& forceFieldParams, - const t_forcerec& forceRec); +std::vector makeLJPmeC6GridCorrectionParameters(int numAtomTypes, + gmx::ArrayRef iparams, + LongRangeVdW ljpme_combination_rule); /*! \brief Set the number of charge groups and atoms. * -- 2.22.0