*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013-2020, by the GROMACS development team, led by
+ * Copyright (c) 2013-2019,2020, 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.
#include "gromacs/mdlib/forcerec_threading.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/md_support.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/rf_util.h"
#include "gromacs/mdlib/wall.h"
+#include "gromacs/mdlib/wholemoleculetransform.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/group.h"
#include "gromacs/mdtypes/iforceprovider.h"
#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/nbnxm/gpu_data_mgmt.h"
#include "gromacs/nbnxm/nbnxm.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/strconvert.h"
-/*! \brief environment variable to enable GPU P2P communication */
-static const bool c_enableGpuPmePpComms =
- (getenv("GMX_GPU_PME_PP_COMMS") != nullptr) && GMX_THREAD_MPI && (GMX_GPU == GMX_GPU_CUDA);
-
-static real* mk_nbfp(const gmx_ffparams_t* idef, gmx_bool bBHAM)
+static std::vector<real> mk_nbfp(const gmx_ffparams_t* idef, gmx_bool bBHAM)
{
- real* nbfp;
- int i, j, k, atnr;
+ std::vector<real> nbfp;
+ int atnr;
atnr = idef->atnr;
if (bBHAM)
{
- snew(nbfp, 3 * atnr * atnr);
- for (i = k = 0; (i < atnr); i++)
+ nbfp.resize(3 * atnr * atnr);
+ int k = 0;
+ for (int i = 0; (i < atnr); i++)
{
- for (j = 0; (j < atnr); j++, k++)
+ for (int j = 0; (j < atnr); j++, k++)
{
BHAMA(nbfp, atnr, i, j) = idef->iparams[k].bham.a;
BHAMB(nbfp, atnr, i, j) = idef->iparams[k].bham.b;
}
else
{
- snew(nbfp, 2 * atnr * atnr);
- for (i = k = 0; (i < atnr); i++)
+ nbfp.resize(2 * atnr * atnr);
+ int k = 0;
+ for (int i = 0; (i < atnr); i++)
{
- for (j = 0; (j < atnr); j++, k++)
+ for (int j = 0; (j < atnr); j++, k++)
{
/* nbfp now includes the 6.0/12.0 derivative prefactors */
C6(nbfp, atnr, i, j) = idef->iparams[k].lj.c6 * 6.0;
acSETTLE
};
-static cginfo_mb_t* init_cginfo_mb(const gmx_mtop_t* mtop, const t_forcerec* fr, gmx_bool* bFEP_NonBonded)
+static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t* mtop, const t_forcerec* fr)
{
- cginfo_mb_t* cginfo_mb;
- gmx_bool* type_VDW;
- int* cginfo;
- int* a_con;
-
- snew(cginfo_mb, mtop->molblock.size());
+ gmx_bool* type_VDW;
+ int* a_con;
snew(type_VDW, fr->ntype);
for (int ai = 0; ai < fr->ntype; ai++)
}
}
- *bFEP_NonBonded = FALSE;
-
- int a_offset = 0;
+ std::vector<cginfo_mb_t> cginfoPerMolblock;
+ int a_offset = 0;
for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
{
const gmx_molblock_t& molb = mtop->molblock[mb];
const gmx_moltype_t& molt = mtop->moltype[molb.type];
- const t_blocka& excl = molt.excls;
+ const auto& excl = molt.excls;
/* Check if the cginfo is identical for all molecules in this block.
* If so, we only need an array of the size of one molecule.
}
}
- cginfo_mb[mb].cg_start = a_offset;
- cginfo_mb[mb].cg_end = a_offset + molb.nmol * molt.atoms.nr;
- cginfo_mb[mb].cg_mod = (bId ? 1 : molb.nmol) * molt.atoms.nr;
- snew(cginfo_mb[mb].cginfo, cginfo_mb[mb].cg_mod);
- cginfo = cginfo_mb[mb].cginfo;
+ cginfo_mb_t cginfo_mb;
+ cginfo_mb.cg_start = a_offset;
+ cginfo_mb.cg_end = a_offset + molb.nmol * molt.atoms.nr;
+ cginfo_mb.cg_mod = (bId ? 1 : molb.nmol) * molt.atoms.nr;
+ cginfo_mb.cginfo.resize(cginfo_mb.cg_mod);
+ gmx::ArrayRef<int> cginfo = cginfo_mb.cginfo;
/* Set constraints flags for constrained atoms */
snew(a_con, molt.atoms.nr);
bool haveExclusions = false;
/* Loop over all the exclusions of atom ai */
- for (int j = excl.index[a]; j < excl.index[a + 1]; j++)
+ for (const int j : excl[a])
{
- if (excl.a[j] != a)
+ if (j != a)
{
haveExclusions = true;
break;
if (fr->efep != efepNO && PERTURBED(atom))
{
SET_CGINFO_FEP(atomInfo);
- *bFEP_NonBonded = TRUE;
}
}
}
sfree(a_con);
+ cginfoPerMolblock.push_back(cginfo_mb);
+
a_offset += molb.nmol * molt.atoms.nr;
}
sfree(type_VDW);
- return cginfo_mb;
+ return cginfoPerMolblock;
}
-static std::vector<int> cginfo_expand(const int nmb, const cginfo_mb_t* cgi_mb)
+static std::vector<int> cginfo_expand(const int nmb, gmx::ArrayRef<const cginfo_mb_t> cgi_mb)
{
const int ncg = cgi_mb[nmb - 1].cg_end;
return cginfo;
}
-static void done_cginfo_mb(cginfo_mb_t* cginfo_mb, int numMolBlocks)
-{
- if (cginfo_mb == nullptr)
- {
- return;
- }
- for (int mb = 0; mb < numMolBlocks; ++mb)
- {
- sfree(cginfo_mb[mb].cginfo);
- }
- sfree(cginfo_mb);
-}
-
/* Sets the sum of charges (squared) and C6 in the system in fr.
* Returns whether the system has a net charge.
*/
fr->natoms_force = natoms_force;
fr->natoms_force_constr = natoms_force_constr;
- if (fr->natoms_force_constr > fr->nalloc_force)
- {
- fr->nalloc_force = over_alloc_dd(fr->natoms_force_constr);
- }
-
if (fr->haveDirectVirialContributions)
{
fr->forceBufferForDirectVirialContributions.resize(natoms_f_novirsum);
* both accuracy requirements, when relevant.
*/
static void init_ewald_f_table(const interaction_const_t& ic,
+ const real tableExtensionLength,
EwaldCorrectionTables* coulombTables,
EwaldCorrectionTables* vdwTables)
{
*/
const real tableScale = ewald_spline3_table_scale(ic, useCoulombTable, useVdwTable);
- const int tableSize = static_cast<int>(ic.rcoulomb * tableScale) + 2;
+ real tableLen = ic.rcoulomb;
+ if (useCoulombTable && tableExtensionLength > 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.
+ * The alternative is to look through all the exclusions and check if they come from
+ * couple-intramol == no. Meanwhile, always having larger tables should only affect
+ * memory consumption, not speed (barring cache issues).
+ */
+ tableLen = ic.rcoulomb + tableExtensionLength;
+ }
+ const int tableSize = static_cast<int>(tableLen * tableScale) + 2;
if (useCoulombTable)
{
}
}
-void init_interaction_const_tables(FILE* fp, interaction_const_t* ic)
+void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real tableExtensionLength)
{
if (EEL_PME_EWALD(ic->eeltype) || EVDW_PME(ic->vdwtype))
{
- init_ewald_f_table(*ic, ic->coulombEwaldTables.get(), ic->vdwEwaldTables.get());
+ init_ewald_f_table(*ic, tableExtensionLength, ic->coulombEwaldTables.get(),
+ ic->vdwEwaldTables.get());
if (fp != nullptr)
{
fprintf(fp, "Initialized non-bonded Ewald tables, spacing: %.2e size: %zu\n\n",
*interaction_const = ic;
}
-bool areMoleculesDistributedOverPbc(const t_inputrec& ir, const gmx_mtop_t& mtop, const gmx::MDLogger& mdlog)
-{
- bool areMoleculesDistributedOverPbc = false;
- const bool useEwaldSurfaceCorrection = (EEL_PME_EWALD(ir.coulombtype) && ir.epsilon_surface != 0);
-
- const bool bSHAKE =
- (ir.eConstrAlg == econtSHAKE
- && (gmx_mtop_ftype_count(mtop, F_CONSTR) > 0 || gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0));
-
- /* The group cut-off scheme and SHAKE assume charge groups
- * are whole, but not using molpbc is faster in most cases.
- * With intermolecular interactions we need PBC for calculating
- * distances between atoms in different molecules.
- */
- if (bSHAKE && !mtop.bIntermolecularInteractions)
- {
- areMoleculesDistributedOverPbc = ir.bPeriodicMols;
-
- if (areMoleculesDistributedOverPbc)
- {
- gmx_fatal(FARGS, "SHAKE is not supported with periodic molecules");
- }
- }
- else
- {
- /* Not making molecules whole is faster in most cases,
- * but with orientation restraints or non-tinfoil boundary
- * conditions we need whole molecules.
- */
- areMoleculesDistributedOverPbc =
- (gmx_mtop_ftype_count(mtop, F_ORIRES) == 0 && !useEwaldSurfaceCorrection);
-
- if (getenv("GMX_USE_GRAPH") != nullptr)
- {
- areMoleculesDistributedOverPbc = false;
-
- GMX_LOG(mdlog.warning)
- .asParagraph()
- .appendText(
- "GMX_USE_GRAPH is set, using the graph for bonded "
- "interactions");
-
- if (mtop.bIntermolecularInteractions)
- {
- GMX_LOG(mdlog.warning)
- .asParagraph()
- .appendText(
- "WARNING: Molecules linked by intermolecular interactions "
- "have to reside in the same periodic image, otherwise "
- "artifacts will occur!");
- }
- }
-
- GMX_RELEASE_ASSERT(areMoleculesDistributedOverPbc || !mtop.bIntermolecularInteractions,
- "We need to use PBC within molecules with inter-molecular interactions");
-
- if (bSHAKE && areMoleculesDistributedOverPbc)
- {
- gmx_fatal(FARGS,
- "SHAKE is not properly supported with intermolecular interactions. "
- "For short simulations where linked molecules remain in the same "
- "periodic image, the environment variable GMX_USE_GRAPH can be used "
- "to override this check.\n");
- }
- }
-
- return areMoleculesDistributedOverPbc;
-}
-
void init_forcerec(FILE* fp,
const gmx::MDLogger& mdlog,
t_forcerec* fr,
const char* tabfn,
const char* tabpfn,
gmx::ArrayRef<const std::string> tabbfnm,
- const gmx_hw_info_t& hardwareInfo,
- const gmx_device_info_t* deviceInfo,
- const bool useGpuForBonded,
- const bool pmeOnlyRankUsesGpu,
- real print_force,
- gmx_wallcycle* wcycle)
+ real print_force)
{
- real rtab;
- char* env;
- double dbl;
- gmx_bool bFEP_NonBonded;
+ /* The CMake default turns SIMD kernels on, but it might be turned off further down... */
+ fr->use_simd_kernels = GMX_USE_SIMD_KERNELS;
- /* By default we turn SIMD kernels on, but it might be turned off further down... */
- fr->use_simd_kernels = TRUE;
-
- if (check_box(ir->ePBC, box))
+ if (check_box(ir->pbcType, box))
{
- gmx_fatal(FARGS, "%s", check_box(ir->ePBC, box));
+ gmx_fatal(FARGS, "%s", check_box(ir->pbcType, box));
}
/* Test particle insertion ? */
fr->sc_r_power = ir->fepvals->sc_r_power;
fr->sc_sigma6_def = gmx::power6(ir->fepvals->sc_sigma);
- env = getenv("GMX_SCSIGMA_MIN");
+ char* env = getenv("GMX_SCSIGMA_MIN");
if (env != nullptr)
{
- dbl = 0;
+ double dbl = 0;
sscanf(env, "%20lf", &dbl);
fr->sc_sigma6_min = gmx::power6(dbl);
if (fp)
/* Neighbour searching stuff */
fr->cutoff_scheme = ir->cutoff_scheme;
- fr->ePBC = ir->ePBC;
+ fr->pbcType = ir->pbcType;
/* Determine if we will do PBC for distances in bonded interactions */
- if (fr->ePBC == epbcNONE)
+ if (fr->pbcType == PbcType::No)
{
fr->bMolPBC = FALSE;
}
{
const bool useEwaldSurfaceCorrection =
(EEL_PME_EWALD(ir->coulombtype) && ir->epsilon_surface != 0);
+ const bool haveOrientationRestraints = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
if (!DOMAINDECOMP(cr))
{
- fr->bMolPBC = areMoleculesDistributedOverPbc(*ir, *mtop, mdlog);
- // The assert below is equivalent to fcd->orires.nr != gmx_mtop_ftype_count(mtop, F_ORIRES)
- GMX_RELEASE_ASSERT(!fr->bMolPBC || fcd->orires.nr == 0,
- "Molecules broken over PBC exist in a simulation including "
- "orientation restraints. "
- "This likely means that the global topology and the force constant "
- "data have gotten out of sync.");
- if (useEwaldSurfaceCorrection)
+ fr->bMolPBC = true;
+
+ if (useEwaldSurfaceCorrection || haveOrientationRestraints)
{
- gmx_fatal(FARGS,
- "In GROMACS 2020, Ewald dipole correction is disabled when not "
- "using domain decomposition. With domain decomposition, it only works "
- "when each molecule consists of a single update group (e.g. water). "
- "This will be fixed in GROMACS 2021.");
+ fr->wholeMoleculeTransform =
+ std::make_unique<gmx::WholeMoleculeTransform>(*mtop, ir->pbcType);
}
}
else
{
- fr->bMolPBC = dd_bonded_molpbc(cr->dd, fr->ePBC);
+ fr->bMolPBC = dd_bonded_molpbc(cr->dd, fr->pbcType);
- if (useEwaldSurfaceCorrection && !dd_moleculesAreAlwaysWhole(*cr->dd))
+ /* With Ewald surface correction it is useful to support e.g. running water
+ * in parallel with update groups.
+ * With orientation restraints there is no sensible use case supported with DD.
+ */
+ if ((useEwaldSurfaceCorrection && !dd_moleculesAreAlwaysWhole(*cr->dd)) || haveOrientationRestraints)
{
gmx_fatal(FARGS,
- "You requested dipole correction (epsilon_surface > 0), but molecules "
- "are broken "
+ "You requested Ewald surface correction or orientation restraints, "
+ "but molecules are broken "
"over periodic boundary conditions by the domain decomposition. "
"Run without domain decomposition instead.");
}
if (useEwaldSurfaceCorrection)
{
- GMX_RELEASE_ASSERT((!DOMAINDECOMP(cr) && !fr->bMolPBC)
- || (DOMAINDECOMP(cr) && dd_moleculesAreAlwaysWhole(*cr->dd)),
+ GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr) || dd_moleculesAreAlwaysWhole(*cr->dd),
"Molecules can not be broken by PBC with epsilon_surface > 0");
}
}
/* 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);
+ init_interaction_const_tables(fp, fr->ic, ir->tabext);
const interaction_const_t* ic = fr->ic;
case eelSWITCH:
case eelSHIFT:
case eelUSER:
- case eelENCADSHIFT:
case eelPMESWITCH:
case eelPMEUSER:
case eelPMEUSERSWITCH:
case evdwSWITCH:
case evdwSHIFT:
- case evdwUSER:
- case evdwENCADSHIFT:
- fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_CUBICSPLINETABLE;
- break;
+ case evdwUSER: fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_CUBICSPLINETABLE; break;
default: gmx_fatal(FARGS, "Unsupported vdw interaction: %s", evdw_names[ic->vdwtype]);
}
fr->shiftForces.resize(SHIFTS);
- if (fr->nbfp == nullptr)
+ if (fr->nbfp.empty())
{
fr->ntype = mtop->ffparams.atnr;
fr->nbfp = mk_nbfp(&mtop->ffparams, fr->bBHAM);
* in that case grompp should already have checked that we do not need
* normal tables and we only generate tables for 1-4 interactions.
*/
- rtab = ir->rlist + ir->tabext;
+ real rtab = ir->rlist + ir->tabext;
/* We want to use unmodified tables for 1-4 coulombic
* interactions, so we must in general have an extra set of
}
// QM/MM initialization if requested
- fr->bQMMM = ir->bQMMM;
- if (fr->bQMMM)
+ if (ir->bQMMM)
{
- // Initialize QM/MM if supported
- if (GMX_QMMM)
- {
- GMX_LOG(mdlog.info)
- .asParagraph()
- .appendText(
- "Large parts of the QM/MM support is deprecated, and may be removed in "
- "a future "
- "version. Please get in touch with the developers if you find the "
- "support useful, "
- "as help is needed if the functionality is to continue to be "
- "available.");
- fr->qr = mk_QMMMrec();
- init_QMMMrec(cr, mtop, ir, fr);
- }
- else
- {
- gmx_incons(
- "QM/MM was requested, but is only available when GROMACS "
- "is configured with QM/MM support");
- }
+ gmx_incons("QM/MM was requested, but is no longer available in GROMACS");
}
/* Set all the static charge group info */
- fr->cginfo_mb = init_cginfo_mb(mtop, fr, &bFEP_NonBonded);
+ fr->cginfo_mb = init_cginfo_mb(mtop, fr);
if (!DOMAINDECOMP(cr))
{
fr->cginfo = cginfo_expand(mtop->molblock.size(), fr->cginfo_mb);
fr->nthread_ewc = gmx_omp_nthreads_get(emntBonded);
snew(fr->ewc_t, fr->nthread_ewc);
- if (fr->cutoff_scheme == ecutsVERLET)
- {
- // We checked the cut-offs in grompp, but double-check here.
- // We have PME+LJcutoff kernels for rcoulomb>rvdw.
- if (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == eelCUT)
- {
- GMX_RELEASE_ASSERT(ir->rcoulomb >= ir->rvdw,
- "With Verlet lists and PME we should have rcoulomb>=rvdw");
- }
- else
- {
- GMX_RELEASE_ASSERT(
- ir->rcoulomb == ir->rvdw,
- "With Verlet lists and no PME rcoulomb and rvdw should be identical");
- }
-
- fr->nbv = Nbnxm::init_nb_verlet(mdlog, bFEP_NonBonded, ir, fr, cr, hardwareInfo, deviceInfo,
- mtop, box, wcycle);
-
- if (useGpuForBonded)
- {
- auto stream = havePPDomainDecomposition(cr)
- ? Nbnxm::gpu_get_command_stream(
- fr->nbv->gpu_nbv, gmx::InteractionLocality::NonLocal)
- : Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv,
- gmx::InteractionLocality::Local);
- // TODO the heap allocation is only needed while
- // t_forcerec lacks a constructor.
- fr->gpuBonded = new gmx::GpuBonded(mtop->ffparams, stream, wcycle);
- }
- }
-
if (ir->eDispCorr != edispcNO)
{
fr->dispersionCorrection = std::make_unique<DispersionCorrection>(
- *mtop, *ir, fr->bBHAM, fr->ntype,
- gmx::arrayRefFromArray(fr->nbfp, fr->ntype * fr->ntype * 2), *fr->ic, tabfn);
+ *mtop, *ir, fr->bBHAM, fr->ntype, fr->nbfp, *fr->ic, tabfn);
fr->dispersionCorrection->print(mdlog);
}
*/
fprintf(fp, "\n");
}
-
- if (pmeOnlyRankUsesGpu && c_enableGpuPmePpComms)
- {
- fr->pmePpCommGpu = std::make_unique<gmx::PmePpCommGpu>(cr->mpi_comm_mysim, cr->dd->pme_nodeid);
- }
}
t_forcerec::t_forcerec() = default;
-t_forcerec::~t_forcerec() = default;
-
-/* Frees GPU memory and sets a tMPI node barrier.
- *
- * Note that this function needs to be called even if GPUs are not used
- * in this run because the PME ranks have no knowledge of whether GPUs
- * are used or not, but all ranks need to enter the barrier below.
- * \todo Remove physical node barrier from this function after making sure
- * that it's not needed anymore (with a shared GPU run).
- */
-void free_gpu_resources(t_forcerec* fr,
- const gmx::PhysicalNodeCommunicator& physicalNodeCommunicator,
- const gmx_gpu_info_t& gpu_info)
-{
- bool isPPrankUsingGPU = (fr != nullptr) && (fr->nbv != nullptr) && fr->nbv->useGpu();
-
- /* stop the GPU profiler (only CUDA) */
- if (gpu_info.n_dev > 0)
- {
- stopGpuProfiler();
- }
-
- if (isPPrankUsingGPU)
- {
- /* Free data in GPU memory and pinned memory before destroying the GPU context */
- fr->nbv.reset();
-
- delete fr->gpuBonded;
- fr->gpuBonded = nullptr;
- }
-
- /* With tMPI we need to wait for all ranks to finish deallocation before
- * destroying the CUDA context in free_gpu() as some tMPI ranks may be sharing
- * GPU and context.
- *
- * This is not a concern in OpenCL where we use one context per rank which
- * is freed in nbnxn_gpu_free().
- *
- * Note: it is safe to not call the barrier on the ranks which do not use GPU,
- * but it is easier and more futureproof to call it on the whole node.
- */
- if (GMX_THREAD_MPI)
- {
- physicalNodeCommunicator.barrier();
- }
-}
-
-void done_forcerec(t_forcerec* fr, int numMolBlocks)
+t_forcerec::~t_forcerec()
{
- if (fr == nullptr)
- {
- // PME-only ranks don't have a forcerec
- return;
- }
- done_cginfo_mb(fr->cginfo_mb, numMolBlocks);
- sfree(fr->nbfp);
- delete fr->ic;
- sfree(fr->shift_vec);
- sfree(fr->ewc_t);
- tear_down_bonded_threading(fr->bondedThreading);
- GMX_RELEASE_ASSERT(fr->gpuBonded == nullptr, "Should have been deleted earlier, when used");
- fr->bondedThreading = nullptr;
- delete fr;
+ /* Note: This code will disappear when types are converted to C++ */
+ sfree(shift_vec);
+ sfree(ewc_t);
+ tear_down_bonded_threading(bondedThreading);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff