From: Berk Hess Date: Thu, 10 Dec 2020 13:34:59 +0000 (+0000) Subject: A check for perturbed listed pairs beyond rlist X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=da9026923650548e3b421eb144e7ec02fced5ee1;p=alexxy%2Fgromacs.git A check for perturbed listed pairs beyond rlist Excluded atoms pairs should be within the pairlist cut-off, as otherwise long-range corrections are not applied. This is problematic for free-energy calculations where intra-molecular interactions are not coupled and replaced by pair interactions. Now a fatal error is generated when perturbed pair interactions are beyond rlist. Also added PME grid correction for LJ excluded pairs beyond rvdw, which was missing from #3403. Added missing release note for the fix of #3403. Fixes #3403 Fixes #3809 --- diff --git a/api/nblib/gmxsetup.cpp b/api/nblib/gmxsetup.cpp index 8501ff062a..f6b3209500 100644 --- a/api/nblib/gmxsetup.cpp +++ b/api/nblib/gmxsetup.cpp @@ -280,7 +280,7 @@ void NbvSetupUtil::setupInteractionConst(const NBKernelOptions& options) } gmxForceCalculator_->interactionConst_->coulombEwaldTables = std::make_unique(); - init_interaction_const_tables(nullptr, gmxForceCalculator_->interactionConst_.get(), 0); + init_interaction_const_tables(nullptr, gmxForceCalculator_->interactionConst_.get(), 0, 0); } } diff --git a/docs/release-notes/2021/major/bugs-fixed.rst b/docs/release-notes/2021/major/bugs-fixed.rst index a38d2f26ba..5fec2401a9 100644 --- a/docs/release-notes/2021/major/bugs-fixed.rst +++ b/docs/release-notes/2021/major/bugs-fixed.rst @@ -25,6 +25,18 @@ H atoms in particular. :issue:`3469` +Correct excluded perturbed interactions beyond the non-bonded cut-off distance +"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +With free-energy calculations without coupling of intermolecular interactions, +non-bonded pair interactions at distance longer than the cut-off distance can +be excluded. These interactions would still have PME long-range contributions. +The contributions are now removed. In addition, mdrun will stop with a fatal +error when interactions beyond the pair-list cut-off are present. + +:issue:`3403` +:issue:`3808` + Corrected AWH initial histogram size """""""""""""""""""""""""""""""""""" diff --git a/src/gromacs/ewald/pme_load_balancing.cpp b/src/gromacs/ewald/pme_load_balancing.cpp index 73b40cc43a..542c730156 100644 --- a/src/gromacs/ewald/pme_load_balancing.cpp +++ b/src/gromacs/ewald/pme_load_balancing.cpp @@ -833,7 +833,7 @@ static void pme_load_balance(pme_load_balancing_t* pme_lb, } /* We always re-initialize the tables whether they are used or not */ - init_interaction_const_tables(nullptr, ic, ir.tabext); + init_interaction_const_tables(nullptr, ic, set->rlistOuter, ir.tabext); Nbnxm::gpu_pme_loadbal_update_param(nbv, ic); diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp index 53f3dbf9d7..6573583c7e 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp @@ -375,6 +375,10 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, real* gmx_restrict f = &(forceWithShiftForces->force()[0][0]); real* gmx_restrict fshift = &(forceWithShiftForces->shiftForces()[0][0]); + const real rlistSquared = gmx::square(fr->rlist); + + int numExcludedPairsBeyondRlist = 0; + for (int n = 0; n < nri; n++) { int npair_within_cutoff = 0; @@ -431,6 +435,11 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, } npair_within_cutoff++; + if (rsq > rlistSquared) + { + numExcludedPairsBeyondRlist++; + } + if (rsq > 0) { /* Note that unlike in the nbnxn kernels, we do not need @@ -722,7 +731,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, } } - if (vdwInteractionTypeIsEwald && r < rvdw) + if (vdwInteractionTypeIsEwald && (r < rvdw || !bPairIncluded)) { /* See comment in the preamble. When using LJ-Ewald interactions * (unless we use a switch modifier) we subtract the reciprocal-space @@ -839,6 +848,19 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, */ #pragma omp atomic inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri * 12 + nlist->jindex[nri] * 150); + + if (numExcludedPairsBeyondRlist > 0) + { + gmx_fatal(FARGS, + "There are %d perturbed non-bonded pair interactions beyond the pair-list cutoff " + "of %g nm, which is not supported. This can happen because the system is " + "unstable or because intra-molecular interactions at long distances are " + "excluded. If the " + "latter is the case, you can try to increase nstlist or rlist to avoid this." + "The error is likely triggered by the use of couple-intramol=no " + "and the maximal distance in the decoupled molecule exceeding rlist.", + numExcludedPairsBeyondRlist, fr->rlist); + } } typedef void (*KernelFunction)(const t_nblist* gmx_restrict nlist, diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index 0a9f0dc89f..8d7eb27108 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -724,7 +724,8 @@ static void initVdwEwaldParameters(FILE* fp, const t_inputrec* ir, interaction_c * both accuracy requirements, when relevant. */ static void init_ewald_f_table(const interaction_const_t& ic, - const real tableExtensionLength, + const real rlist, + const real tabext, EwaldCorrectionTables* coulombTables, EwaldCorrectionTables* vdwTables) { @@ -739,7 +740,7 @@ static void init_ewald_f_table(const interaction_const_t& ic, const bool havePerturbedNonbondeds = (ic.softCoreParameters != nullptr); real tableLen = ic.rcoulomb; - if (useCoulombTable && havePerturbedNonbondeds && tableExtensionLength > 0.0) + if ((useCoulombTable || useVdwTable) && havePerturbedNonbondeds && rlist + tabext > 0.0) { /* TODO: Ideally this should also check if couple-intramol == no, but that isn't * stored in ir. Grompp puts that info into an opts structure that doesn't make it into the tpr. @@ -747,7 +748,7 @@ static void init_ewald_f_table(const interaction_const_t& ic, * couple-intramol == no. Meanwhile, always having larger tables should only affect * memory consumption, not speed (barring cache issues). */ - tableLen = ic.rcoulomb + tableExtensionLength; + tableLen = rlist + tabext; } const int tableSize = static_cast(tableLen * tableScale) + 2; @@ -763,11 +764,11 @@ static void init_ewald_f_table(const interaction_const_t& ic, } } -void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real tableExtensionLength) +void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real rlist, const real tableExtensionLength) { if (EEL_PME_EWALD(ic->eeltype) || EVDW_PME(ic->vdwtype)) { - init_ewald_f_table(*ic, tableExtensionLength, ic->coulombEwaldTables.get(), + init_ewald_f_table(*ic, rlist, tableExtensionLength, ic->coulombEwaldTables.get(), ic->vdwEwaldTables.get()); if (fp != nullptr) { @@ -1081,7 +1082,7 @@ void init_forcerec(FILE* fp, /* fr->ic is used both by verlet and group kernels (to some extent) now */ init_interaction_const(fp, &fr->ic, ir, mtop, systemHasNetCharge); - init_interaction_const_tables(fp, fr->ic, ir->tabext); + init_interaction_const_tables(fp, fr->ic, fr->rlist, ir->tabext); const interaction_const_t* ic = fr->ic; diff --git a/src/gromacs/mdlib/forcerec.h b/src/gromacs/mdlib/forcerec.h index aed645f186..1a01ff1413 100644 --- a/src/gromacs/mdlib/forcerec.h +++ b/src/gromacs/mdlib/forcerec.h @@ -82,9 +82,10 @@ void forcerec_set_ranges(t_forcerec* fr, int natoms_force, int natoms_force_cons * Initializes the tables in the interaction constant data structure. * \param[in] fp File for debugging output * \param[in] ic Structure holding the table constant + * \param[in] rlist Length of the neighbour list * \param[in] tableExtensionLength Length by which to extend the tables. Taken from the input record. */ -void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, real tableExtensionLength); +void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, real rlist, real tableExtensionLength); /*! \brief Initialize forcerec structure. * diff --git a/src/gromacs/nbnxm/benchmark/bench_setup.cpp b/src/gromacs/nbnxm/benchmark/bench_setup.cpp index 9fffb257ea..af0b292699 100644 --- a/src/gromacs/nbnxm/benchmark/bench_setup.cpp +++ b/src/gromacs/nbnxm/benchmark/bench_setup.cpp @@ -165,7 +165,7 @@ static interaction_const_t setupInteractionConst(const KernelBenchOptions& optio GMX_RELEASE_ASSERT(options.ewaldcoeff_q > 0, "Ewald coefficient should be > 0"); ic.ewaldcoeff_q = options.ewaldcoeff_q; ic.coulombEwaldTables = std::make_unique(); - init_interaction_const_tables(nullptr, &ic, 0); + init_interaction_const_tables(nullptr, &ic, 0, 0); } return ic;