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<const int64_t> 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
{
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<const int64_t> atomInfo, gmx_localtop_t* ltop);
#endif
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
{
{
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;
+ }
}
}
}
* 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;
#include <cstring>
#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"
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<const int> 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<real> qA(mtop.natoms);
- std::vector<real> qB(mtop.natoms);
- for (const AtomProxy atomP : AtomRange(mtop))
+ std::vector<int64_t> 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,
/* 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.
/*
* 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.
* \brief
* Implements test of mtop routines
*
- * \author Roland Schulz <roland.schulz@intel.com>
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
* \ingroup module_topology
*/
#include "gmxpre.h"
+#include "gromacs/topology/mtop_util.h"
+
#include <gtest/gtest.h>
-#include "gromacs/topology/mtop_util.h"
-#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/topology/ifunc.h"
namespace gmx
{
}
}
+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<int, 2> 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<int, 2> 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<int, 2> 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
TEST(TopSortTest, WorksOnEmptyIdef)
{
- gmx_ffparams_t forcefieldParams;
- InteractionDefinitions idef(forcefieldParams);
- real* emptyChargeA = nullptr;
- real* emptyChargeB = nullptr;
+ gmx_ffparams_t forcefieldParams;
+ InteractionDefinitions idef(forcefieldParams);
+ ArrayRef<const int64_t> 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());
// F_LJ14
forcefieldParams.iparams.push_back(makeUnperturbedLJ14Params(100, 10000));
idef.il[F_LJ14].push_back(forcefieldParams.iparams.size() - 1, bondAtoms);
- std::vector<real> chargeA{ -1, 2 };
- std::vector<real> chargeB{ -1, 2 };
+ std::vector<int64_t> 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";
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<real> chargeA{ -1, 2, -2, 0 };
- std::vector<real> chargeB{ -1, 2, 4, 0 };
+ std::vector<int64_t> 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";
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<real> chargeA{ -1, 2, -2, 0 };
- std::vector<real> chargeB{ -1, 2, 4, 0 };
+ std::vector<int64_t> 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";
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<real> chargeA{ -1, 2, -2, 0, 1, -2, 3, -3 };
- std::vector<real> chargeB{ -1, 2, 4, 0, 1, -2, 3, -3 };
+ std::vector<int64_t> 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";
#include <cstdio>
+#include "gromacs/mdtypes/atominfo.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/arrayref.h"
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<const int64_t> 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)
return bPert;
}
-void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* qB)
+void gmx_sort_ilist_fe(InteractionDefinitions* idef, gmx::ArrayRef<const int64_t> atomInfo)
{
- if (qB == nullptr)
- {
- qB = qA;
- }
-
bool havePerturbedInteractions = false;
int iabuf_nalloc = 0;
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)
* 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.
struct gmx_mtop_t;
class InteractionDefinitions;
+namespace gmx
+{
+template<typename>
+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<const int64_t> atomInfo);
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff