From bb4acb6f88487ff70e543b3a3ea06eb52131a532 Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Mon, 7 Jun 2021 23:37:24 +0000 Subject: [PATCH] Add hasPerturbedChargeIn14Interaction to atomInfo bits --- src/gromacs/domdec/localtopology.cpp | 4 +- src/gromacs/domdec/localtopology.h | 2 +- src/gromacs/domdec/mdsetup.cpp | 2 +- src/gromacs/mdlib/forcerec.cpp | 11 +++- src/gromacs/mdtypes/atominfo.h | 15 +++--- src/gromacs/topology/mtop_util.cpp | 51 +++++++++++++++--- src/gromacs/topology/mtop_util.h | 14 +++++ src/gromacs/topology/tests/mtop.cpp | 74 ++++++++++++++++++++++++-- src/gromacs/topology/tests/topsort.cpp | 30 +++++------ src/gromacs/topology/topsort.cpp | 23 ++++---- src/gromacs/topology/topsort.h | 10 +++- 11 files changed, 179 insertions(+), 57 deletions(-) diff --git a/src/gromacs/domdec/localtopology.cpp b/src/gromacs/domdec/localtopology.cpp index ea4a68c725..dfe8f85703 100644 --- a/src/gromacs/domdec/localtopology.cpp +++ b/src/gromacs/domdec/localtopology.cpp @@ -1019,11 +1019,11 @@ int dd_make_local_top(gmx_domdec_t* dd, return numBondedInteractionsToReduce; } -void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop) +void dd_sort_local_top(const gmx_domdec_t& dd, gmx::ArrayRef atomInfo, gmx_localtop_t* ltop) { if (dd.reverse_top->doSorting()) { - gmx_sort_ilist_fe(<op->idef, mdatoms->chargeA, mdatoms->chargeB); + gmx_sort_ilist_fe(<op->idef, atomInfo); } else { diff --git a/src/gromacs/domdec/localtopology.h b/src/gromacs/domdec/localtopology.h index 6adfc26554..8a66ea43ed 100644 --- a/src/gromacs/domdec/localtopology.h +++ b/src/gromacs/domdec/localtopology.h @@ -74,6 +74,6 @@ int dd_make_local_top(struct gmx_domdec_t* dd, gmx_localtop_t* ltop); /*! \brief Sort ltop->ilist when we are doing free energy. */ -void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop); +void dd_sort_local_top(const gmx_domdec_t& dd, gmx::ArrayRef atomInfo, gmx_localtop_t* ltop); #endif diff --git a/src/gromacs/domdec/mdsetup.cpp b/src/gromacs/domdec/mdsetup.cpp index 631d167227..1176ddb4e1 100644 --- a/src/gromacs/domdec/mdsetup.cpp +++ b/src/gromacs/domdec/mdsetup.cpp @@ -110,7 +110,7 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, t_mdatoms* mdatoms = mdAtoms->mdatoms(); if (usingDomDec) { - dd_sort_local_top(*cr->dd, mdatoms, top); + dd_sort_local_top(*cr->dd, fr->atomInfo, top); } else { diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index 3cbb1124e7..f5e784134a 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -330,9 +330,16 @@ makeAtomInfoForEachMoleculeBlock(const gmx_mtop_t& mtop, const t_forcerec* fr) { atomInfo |= gmx::sc_atomInfo_HasCharge; } - if (fr->efep != FreeEnergyPerturbationType::No && PERTURBED(atom)) + if (fr->efep != FreeEnergyPerturbationType::No) { - atomInfo |= gmx::sc_atomInfo_FreeEnergyPerturbation; + if (PERTURBED(atom)) + { + atomInfo |= gmx::sc_atomInfo_FreeEnergyPerturbation; + } + if (atomHasPerturbedChargeIn14Interaction(a, molt)) + { + atomInfo |= gmx::sc_atomInfo_HasPerturbedChargeIn14Interaction; + } } } } diff --git a/src/gromacs/mdtypes/atominfo.h b/src/gromacs/mdtypes/atominfo.h index 8697dc2501..104b1f072e 100644 --- a/src/gromacs/mdtypes/atominfo.h +++ b/src/gromacs/mdtypes/atominfo.h @@ -59,13 +59,14 @@ namespace gmx * here, reserving bits 0-7 for the energy-group ID. */ //! \{ -static constexpr int64_t sc_atomInfo_FreeEnergyPerturbation = 1 << 15; -static constexpr int64_t sc_atomInfo_Exclusion = 1 << 17; -static constexpr int64_t sc_atomInfo_Constraint = 1 << 20; -static constexpr int64_t sc_atomInfo_Settle = 1 << 21; -static constexpr int64_t sc_atomInfo_BondCommunication = 1 << 22; -static constexpr int64_t sc_atomInfo_HasVdw = 1 << 23; -static constexpr int64_t sc_atomInfo_HasCharge = 1 << 24; +static constexpr int64_t sc_atomInfo_FreeEnergyPerturbation = 1 << 15; +static constexpr int64_t sc_atomInfo_HasPerturbedChargeIn14Interaction = 1 << 16; +static constexpr int64_t sc_atomInfo_Exclusion = 1 << 17; +static constexpr int64_t sc_atomInfo_Constraint = 1 << 20; +static constexpr int64_t sc_atomInfo_Settle = 1 << 21; +static constexpr int64_t sc_atomInfo_BondCommunication = 1 << 22; +static constexpr int64_t sc_atomInfo_HasVdw = 1 << 23; +static constexpr int64_t sc_atomInfo_HasCharge = 1 << 24; //! \} //! The first 8 bits are reserved for energy-group ID static constexpr int64_t sc_atomInfo_EnergyGroupIdMask = 0b11111111; diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp index 1e92998ad7..30e1b96147 100644 --- a/src/gromacs/topology/mtop_util.cpp +++ b/src/gromacs/topology/mtop_util.cpp @@ -45,6 +45,7 @@ #include #include "gromacs/math/vectypes.h" +#include "gromacs/mdtypes/atominfo.h" #include "gromacs/topology/atoms.h" #include "gromacs/topology/block.h" #include "gromacs/topology/exclusionblocks.h" @@ -828,18 +829,52 @@ static void addMimicExclusions(gmx::ListOfLists* excls, const gmx::ArrayRef gmx::mergeExclusions(excls, qmexcl2); } +bool atomHasPerturbedChargeIn14Interaction(const int atomIndex, const gmx_moltype_t& molt) +{ + if (molt.atoms.nr > 0) + { + // Is the charge perturbed at all? + const t_atom& atom = molt.atoms.atom[atomIndex]; + if (atom.q != atom.qB) + { + // Loop over 1-4 interactions + const InteractionList& ilist = molt.ilist[F_LJ14]; + gmx::ArrayRef iatoms = ilist.iatoms; + const int nral = NRAL(F_LJ14); + for (int i = 0; i < ilist.size(); i += nral + 1) + { + // Compare the atom indices in this 1-4 interaction to + // atomIndex. + int ia1 = iatoms[i + 1]; + int ia2 = iatoms[i + 2]; + if ((atomIndex == ia1) || (atomIndex == ia2)) + { + return true; + } + } + } + } + return false; +} + static void sortFreeEnergyInteractionsAtEnd(const gmx_mtop_t& mtop, InteractionDefinitions* idef) { - std::vector qA(mtop.natoms); - std::vector qB(mtop.natoms); - for (const AtomProxy atomP : AtomRange(mtop)) + std::vector atomInfo(mtop.natoms, 0); + + for (size_t mb = 0; mb < mtop.molblock.size(); mb++) { - const t_atom& local = atomP.atom(); - int index = atomP.globalAtomNumber(); - qA[index] = local.q; - qB[index] = local.qB; + const gmx_molblock_t& molb = mtop.molblock[mb]; + const gmx_moltype_t& molt = mtop.moltype[molb.type]; + for (int a = 0; a < molt.atoms.nr; a++) + { + if (atomHasPerturbedChargeIn14Interaction(a, molt)) + { + atomInfo[a] |= gmx::sc_atomInfo_HasPerturbedChargeIn14Interaction; + break; + } + } } - gmx_sort_ilist_fe(idef, qA.data(), qB.data()); + gmx_sort_ilist_fe(idef, atomInfo); } static void gen_local_top(const gmx_mtop_t& mtop, diff --git a/src/gromacs/topology/mtop_util.h b/src/gromacs/topology/mtop_util.h index 4d0e877946..22c90e0cdb 100644 --- a/src/gromacs/topology/mtop_util.h +++ b/src/gromacs/topology/mtop_util.h @@ -262,6 +262,20 @@ gmx::EnumerationArray gmx_mtop_particletype_count(const gmx_m /* Returns a single t_atoms struct for the whole system */ t_atoms gmx_mtop_global_atoms(const gmx_mtop_t& mtop); +/*! \brief Return whether the atom with the given index is within a + * 1-4 interaction within the given molecule type and its charge is + * perturbed. + * + * Some 1-4 interactions like F_COUL14 have the charges stored in the + * iparams list, but others do not. The commonly used F_LJ14 gets its + * charges from the topology, so we need more detailed checks for it, + * when FEP is active. + * + * \param[in] atomIndex Index within range [0, molt.nr) + * \param[in] molt Moleculetype to consider + * \return Whether this atom is in a 1-4 interaction and charge is perturbed + * */ +bool atomHasPerturbedChargeIn14Interaction(int atomIndex, const gmx_moltype_t& molt); /*! \brief * Populate a 'local' topology for the whole system. diff --git a/src/gromacs/topology/tests/mtop.cpp b/src/gromacs/topology/tests/mtop.cpp index 08713b39eb..a52f080133 100644 --- a/src/gromacs/topology/tests/mtop.cpp +++ b/src/gromacs/topology/tests/mtop.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 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. @@ -36,15 +36,16 @@ * \brief * Implements test of mtop routines * - * \author Roland Schulz + * \author Mark Abraham * \ingroup module_topology */ #include "gmxpre.h" +#include "gromacs/topology/mtop_util.h" + #include -#include "gromacs/topology/mtop_util.h" -#include "gromacs/utility/basedefinitions.h" +#include "gromacs/topology/ifunc.h" namespace gmx { @@ -163,6 +164,71 @@ TEST(MtopTest, CanFindResidueStartAndEndAtoms) } } +TEST(MtopTest, AtomHasPerturbedChargeIn14Interaction) +{ + gmx_moltype_t molt; + molt.atoms.nr = 0; + molt.atoms.atom = nullptr; + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(0, molt)) << "empty moltype has none"; + + { + SCOPED_TRACE("Use a moltype that has no perturbed charges at all"); + const int numAtoms = 4; + init_t_atoms(&molt.atoms, numAtoms, false); + for (int i = 0; i != numAtoms; ++i) + { + molt.atoms.atom[i].q = 1.0; + molt.atoms.atom[i].qB = 1.0; + } + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(0, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(1, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(2, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(3, molt)); + } + + { + SCOPED_TRACE("Use a moltype that has a perturbed charge, but no interactions"); + molt.atoms.atom[1].q = 1.0; + molt.atoms.atom[1].qB = -1.0; + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(0, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(1, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(2, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(3, molt)); + } + + { + SCOPED_TRACE("Use a moltype that has a perturbed charge but no 1-4 interactions"); + std::array bondAtoms = { 0, 1 }; + molt.ilist[F_BONDS].push_back(0, bondAtoms); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(0, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(1, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(2, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(3, molt)); + } + + { + SCOPED_TRACE( + "Use a moltype that has a perturbed charge, but not on an atom with a 1-4 " + "interaction"); + std::array pairAtoms = { 2, 3 }; + molt.ilist[F_LJ14].push_back(0, pairAtoms); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(0, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(1, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(2, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(3, molt)); + } + + { + SCOPED_TRACE("Use a moltype that has a perturbed charge on an atom with a 1-4 interaction"); + std::array perturbedPairAtoms = { 1, 2 }; + molt.ilist[F_LJ14].push_back(0, perturbedPairAtoms); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(0, molt)); + EXPECT_TRUE(atomHasPerturbedChargeIn14Interaction(1, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(2, molt)); + EXPECT_FALSE(atomHasPerturbedChargeIn14Interaction(3, molt)); + } +} + } // namespace } // namespace gmx diff --git a/src/gromacs/topology/tests/topsort.cpp b/src/gromacs/topology/tests/topsort.cpp index 577f154b20..f9efe2bfc7 100644 --- a/src/gromacs/topology/tests/topsort.cpp +++ b/src/gromacs/topology/tests/topsort.cpp @@ -63,12 +63,11 @@ using testing::Pointwise; TEST(TopSortTest, WorksOnEmptyIdef) { - gmx_ffparams_t forcefieldParams; - InteractionDefinitions idef(forcefieldParams); - real* emptyChargeA = nullptr; - real* emptyChargeB = nullptr; + gmx_ffparams_t forcefieldParams; + InteractionDefinitions idef(forcefieldParams); + ArrayRef emptyAtomInfo; - gmx_sort_ilist_fe(&idef, emptyChargeA, emptyChargeB); + gmx_sort_ilist_fe(&idef, emptyAtomInfo); EXPECT_EQ(0, idef.numNonperturbedInteractions[F_BONDS]) << "empty idef has no perturbed bonds"; EXPECT_THAT(idef.il[F_BONDS].iatoms, IsEmpty()); @@ -120,10 +119,9 @@ TEST(TopSortTest, WorksOnIdefWithNoPerturbedInteraction) // F_LJ14 forcefieldParams.iparams.push_back(makeUnperturbedLJ14Params(100, 10000)); idef.il[F_LJ14].push_back(forcefieldParams.iparams.size() - 1, bondAtoms); - std::vector chargeA{ -1, 2 }; - std::vector chargeB{ -1, 2 }; + std::vector atomInfo{ 0, 0 }; - gmx_sort_ilist_fe(&idef, chargeA.data(), chargeB.data()); + gmx_sort_ilist_fe(&idef, atomInfo); EXPECT_EQ(1, idef.numNonperturbedInteractions[F_BONDS] / (NRAL(F_BONDS) + 1)) << "idef with no perturbed bonds has no perturbed bonds"; @@ -155,10 +153,9 @@ TEST(TopSortTest, WorksOnIdefWithPerturbedInteractions) forcefieldParams.iparams.push_back(makePerturbedLJ14Params(100, 10000, 200, 20000)); idef.il[F_LJ14].push_back(forcefieldParams.iparams.size() - 1, perturbedBondAtoms); // Perturb the charge of atom 2, affecting the non-perturbed LJ14 above - std::vector chargeA{ -1, 2, -2, 0 }; - std::vector chargeB{ -1, 2, 4, 0 }; + std::vector atomInfo{ 0, 0, sc_atomInfo_HasPerturbedChargeIn14Interaction, 0 }; - gmx_sort_ilist_fe(&idef, chargeA.data(), chargeB.data()); + gmx_sort_ilist_fe(&idef, atomInfo); EXPECT_EQ(1, idef.numNonperturbedInteractions[F_BONDS] / (NRAL(F_BONDS) + 1)) << "idef with a perturbed bond has a perturbed bond"; @@ -191,10 +188,9 @@ TEST(TopSortTest, SortsIdefWithPerturbedInteractions) forcefieldParams.iparams.push_back(makeUnperturbedLJ14Params(100, 10000)); idef.il[F_LJ14].push_back(forcefieldParams.iparams.size() - 1, bondAtoms); // Perturb the charge of atom 2, affecting the non-perturbed LJ14 above - std::vector chargeA{ -1, 2, -2, 0 }; - std::vector chargeB{ -1, 2, 4, 0 }; + std::vector atomInfo{ 0, 0, sc_atomInfo_HasPerturbedChargeIn14Interaction, 0 }; - gmx_sort_ilist_fe(&idef, chargeA.data(), chargeB.data()); + gmx_sort_ilist_fe(&idef, atomInfo); EXPECT_EQ(1, idef.numNonperturbedInteractions[F_BONDS] / (NRAL(F_BONDS) + 1)) << "idef with a perturbed bond has a perturbed bond"; @@ -236,10 +232,10 @@ TEST(TopSortTest, SortsMoreComplexIdefWithPerturbedInteractions) forcefieldParams.iparams.push_back(makeUnperturbedLJ14Params(100, 10000)); idef.il[F_LJ14].push_back(forcefieldParams.iparams.size() - 1, moreBondAtoms); // Perturb the charge of atom 2, affecting the non-perturbed LJ14 above - std::vector chargeA{ -1, 2, -2, 0, 1, -2, 3, -3 }; - std::vector chargeB{ -1, 2, 4, 0, 1, -2, 3, -3 }; + std::vector atomInfo{ 0, 0, sc_atomInfo_HasPerturbedChargeIn14Interaction, 0, 0, 0, + 0, 0 }; - gmx_sort_ilist_fe(&idef, chargeA.data(), chargeB.data()); + gmx_sort_ilist_fe(&idef, atomInfo); EXPECT_EQ(2, idef.numNonperturbedInteractions[F_BONDS] / (NRAL(F_BONDS) + 1)) << "idef with some perturbed bonds has some perturbed bonds"; diff --git a/src/gromacs/topology/topsort.cpp b/src/gromacs/topology/topsort.cpp index 8270ca0307..98048319e6 100644 --- a/src/gromacs/topology/topsort.cpp +++ b/src/gromacs/topology/topsort.cpp @@ -41,6 +41,7 @@ #include +#include "gromacs/mdtypes/atominfo.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/topology.h" #include "gromacs/utility/arrayref.h" @@ -131,13 +132,12 @@ static gmx_bool ip_pert(int ftype, const t_iparams* ip) return bPert; } -static gmx_bool ip_q_pert(int ftype, const t_iatom* ia, const t_iparams* ip, const real* qA, const real* qB) + +//! Return whether the two atom indices in a 1-4 interaction have perturbed charges per \c atomInfo +static bool hasPerturbedChargeIn14Interaction(int atom0, int atom1, gmx::ArrayRef atomInfo) { - /* 1-4 interactions do not have the charges stored in the iparams list, - * so we need a separate check for those. - */ - return (ip_pert(ftype, ip + ia[0]) - || (ftype == F_LJ14 && (qA[ia[1]] != qB[ia[1]] || qA[ia[2]] != qB[ia[2]]))); + return (bool(atomInfo[atom0] & gmx::sc_atomInfo_HasPerturbedChargeIn14Interaction) + || bool(atomInfo[atom1] & gmx::sc_atomInfo_HasPerturbedChargeIn14Interaction)); } gmx_bool gmx_mtop_bondeds_free_energy(const gmx_mtop_t* mtop) @@ -176,13 +176,8 @@ gmx_bool gmx_mtop_bondeds_free_energy(const gmx_mtop_t* mtop) return bPert; } -void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* qB) +void gmx_sort_ilist_fe(InteractionDefinitions* idef, gmx::ArrayRef atomInfo) { - if (qB == nullptr) - { - qB = qA; - } - bool havePerturbedInteractions = false; int iabuf_nalloc = 0; @@ -201,7 +196,9 @@ void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* while (i < ilist->size()) { /* Check if this interaction is perturbed */ - if (ip_q_pert(ftype, iatoms + i, idef->iparams.data(), qA, qB)) + if (ip_pert(ftype, idef->iparams.data() + iatoms[i]) + || (ftype == F_LJ14 + && hasPerturbedChargeIn14Interaction(iatoms[i + 1], iatoms[i + 2], atomInfo))) { /* Copy to the perturbed buffer */ if (ib + 1 + nral > iabuf_nalloc) diff --git a/src/gromacs/topology/topsort.h b/src/gromacs/topology/topsort.h index 2237b65903..fd471a95a8 100644 --- a/src/gromacs/topology/topsort.h +++ b/src/gromacs/topology/topsort.h @@ -2,7 +2,7 @@ * This file is part of the GROMACS molecular simulation package. * * Copyright (c) 2008,2009,2010,2013,2014 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, 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. @@ -42,12 +42,18 @@ struct gmx_mtop_t; class InteractionDefinitions; +namespace gmx +{ +template +class ArrayRef; +} // namespace gmx + /* Returns if there are perturbed bonded interactions */ gmx_bool gmx_mtop_bondeds_free_energy(const struct gmx_mtop_t* mtop); /* Sort all the bonded ilists in idef to have the perturbed ones at the end * and set nr_nr_nonperturbed in ilist. */ -void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* qB); +void gmx_sort_ilist_fe(InteractionDefinitions* idef, gmx::ArrayRef atomInfo); #endif -- 2.22.0