From: Berk Hess Date: Wed, 24 Jun 2020 17:35:16 +0000 (+0000) Subject: Add class ListedForces X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=c444d582d185b361cf3c10700a7b642579a5ebd2;p=alexxy%2Fgromacs.git Add class ListedForces Gather all data that is only used for the listed force calculation in a new class called ListedForces. This is a first step in refactoring, more parameters to the listed forces calculation can be moved into class. Also converted bonded_threading_t, t_fcdata and bondedtable_t to C++. Made listed_threading.h internal to the listed_forces module. This change is only refactoring. --- diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index 5170ca8744..3e4996d3b6 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -67,7 +67,6 @@ #include "gromacs/gpu_utils/device_stream_manager.h" #include "gromacs/gpu_utils/gpu_utils.h" #include "gromacs/hardware/hw_info.h" -#include "gromacs/listed_forces/manage_threading.h" #include "gromacs/math/vec.h" #include "gromacs/math/vectypes.h" #include "gromacs/mdlib/calc_verletbuf.h" diff --git a/src/gromacs/domdec/mdsetup.cpp b/src/gromacs/domdec/mdsetup.cpp index d18168947e..3027d3b342 100644 --- a/src/gromacs/domdec/mdsetup.cpp +++ b/src/gromacs/domdec/mdsetup.cpp @@ -39,7 +39,7 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_struct.h" #include "gromacs/ewald/pme.h" -#include "gromacs/listed_forces/manage_threading.h" +#include "gromacs/listed_forces/listed_forces.h" #include "gromacs/mdlib/constr.h" #include "gromacs/mdlib/mdatoms.h" #include "gromacs/mdlib/vsite.h" @@ -131,7 +131,7 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, make_local_shells(cr, mdatoms, shellfc); } - setup_bonded_threading(fr->bondedThreading, fr->natoms_force, fr->gpuBonded != nullptr, top->idef); + fr->listedForces->setup(top->idef, fr->natoms_force, fr->gpuBonded != nullptr); if (EEL_PME(fr->ic->eeltype) && (cr->duty & DUTY_PME)) { diff --git a/src/gromacs/gmxana/gmx_disre.cpp b/src/gromacs/gmxana/gmx_disre.cpp index 6a84cc8a81..3fd41a98dd 100644 --- a/src/gromacs/gmxana/gmx_disre.cpp +++ b/src/gromacs/gmxana/gmx_disre.cpp @@ -158,7 +158,7 @@ static void check_viol(FILE* log, int isize, const int index[], real vvindex[], - t_fcdata* fcd) + t_disresdata* disresdata) { int i, j, nat, n, type, nviol, ndr, label; real rt, mviol, tviol, viol, lam, dvdl, drt; @@ -201,6 +201,11 @@ static void check_viol(FILE* log, i / nat, label, label_old, label_old + 1); } } + + // Set up t_fcdata, only needed for calling ta_disres() + t_fcdata fcd; + fcd.disres = disresdata; + // Get offset for label index label_old = forceparams[forceatoms[0]].disres.label; for (i = 0; (i < disres.size());) @@ -218,25 +223,25 @@ static void check_viol(FILE* log, } while (((i + n) < disres.size()) && (forceparams[forceatoms[i + n]].disres.label == label + label_old)); - calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i], x, pbc, fcd, nullptr); + calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i], x, pbc, disresdata, nullptr); - if (fcd->disres.Rt_6[label] <= 0) + if (disresdata->Rt_6[label] <= 0) { - gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[label]); + gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, disresdata->Rt_6[label]); } - rt = gmx::invsixthroot(fcd->disres.Rt_6[label]); + rt = gmx::invsixthroot(disresdata->Rt_6[label]); dr[clust_id].aver1[ndr] += rt; dr[clust_id].aver2[ndr] += gmx::square(rt); drt = 1.0 / gmx::power3(rt); dr[clust_id].aver_3[ndr] += drt; - dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label]; + dr[clust_id].aver_6[ndr] += disresdata->Rt_6[label]; snew(fshift, SHIFTS); ta_disres(n, &forceatoms[i], forceparams.data(), x, f, fshift, pbc, lam, &dvdl, nullptr, - fcd, nullptr); + &fcd, nullptr); sfree(fshift); - viol = fcd->disres.sumviol; + viol = disresdata->sumviol; if (viol > 0) { @@ -254,7 +259,7 @@ static void check_viol(FILE* log, { if (index[j] == forceparams[type].disres.label) { - vvindex[j] = gmx::invsixthroot(fcd->disres.Rt_6[label]); + vvindex[j] = gmx::invsixthroot(disresdata->Rt_6[label]); } } } @@ -699,7 +704,6 @@ int gmx_disre(int argc, char* argv[]) }; FILE * out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr; - t_fcdata fcd; int i, j, kkk; t_trxstatus* status; real t; @@ -771,6 +775,7 @@ int gmx_disre(int argc, char* argv[]) gmx_localtop_t top(topInfo.mtop()->ffparams); gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO); + const InteractionDefinitions& idef = top.idef; pbc_null = nullptr; if (ir->pbcType != PbcType::No) @@ -800,20 +805,21 @@ int gmx_disre(int argc, char* argv[]) } ir->dr_tau = 0.0; + t_disresdata disresdata; init_disres(fplog, topInfo.mtop(), ir, DisResRunMode::AnalysisTool, DDRole::Master, - NumRanks::Single, MPI_COMM_NULL, nullptr, &fcd, nullptr, FALSE); + NumRanks::Single, MPI_COMM_NULL, nullptr, &disresdata, nullptr, FALSE); int natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box); snew(f, 5 * natoms); - init_dr_res(&dr, fcd.disres.nres); + init_dr_res(&dr, disresdata.nres); if (opt2bSet("-c", NFILE, fnm)) { clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm)); snew(dr_clust, clust->clust->nr + 1); for (i = 0; (i <= clust->clust->nr); i++) { - init_dr_res(&dr_clust[i], fcd.disres.nres); + init_dr_res(&dr_clust[i], disresdata.nres); } } else @@ -829,7 +835,7 @@ int gmx_disre(int argc, char* argv[]) update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda); if (ir->pbcType != PbcType::No) { - gpbc = gmx_rmpbc_init(top.idef, ir->pbcType, natoms); + gpbc = gmx_rmpbc_init(idef, ir->pbcType, natoms); } j = 0; @@ -858,13 +864,13 @@ int gmx_disre(int argc, char* argv[]) } my_clust = clust->inv_clust[j]; range_check(my_clust, 0, clust->clust->nr); - check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, dr_clust, - my_clust, isize, index, vvindex, &fcd); + check_viol(fplog, idef.il[F_DISRES], idef.iparams, x, f, pbc_null, dr_clust, my_clust, + isize, index, vvindex, &disresdata); } else { - check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, &dr, 0, - isize, index, vvindex, &fcd); + check_viol(fplog, idef.il[F_DISRES], idef.iparams, x, f, pbc_null, &dr, 0, isize, index, + vvindex, &disresdata); } if (bPDB) { @@ -907,19 +913,19 @@ int gmx_disre(int argc, char* argv[]) if (clust) { - dump_clust_stats(fplog, fcd.disres, top.idef.il[F_DISRES], top.idef.iparams, clust->clust, - dr_clust, clust->grpname, isize, index); + dump_clust_stats(fplog, disresdata, idef.il[F_DISRES], idef.iparams, clust->clust, dr_clust, + clust->grpname, isize, index); } else { - dump_stats(fplog, j, fcd.disres, top.idef.il[F_DISRES], top.idef.iparams, &dr, isize, index, + dump_stats(fplog, j, disresdata, idef.il[F_DISRES], idef.iparams, &dr, isize, index, bPDB ? atoms.get() : nullptr); if (bPDB) { write_sto_conf(opt2fn("-q", NFILE, fnm), "Coloured by average violation in Angstrom", atoms.get(), xav, nullptr, ir->pbcType, box); } - dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres, j, top.idef, + dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, disresdata.nres, j, idef, topInfo.mtop(), max_dr, nlevels, bThird); xvgrclose(out); xvgrclose(aver); diff --git a/src/gromacs/gmxana/gmx_nmr.cpp b/src/gromacs/gmxana/gmx_nmr.cpp index bc20833ec6..8d597542cd 100644 --- a/src/gromacs/gmxana/gmx_nmr.cpp +++ b/src/gromacs/gmxana/gmx_nmr.cpp @@ -198,18 +198,18 @@ static void get_orires_parms(const char* topnm, t_inputrec* ir, int* nor, int* n done_top_mtop(&top, &mtop); } -static int get_bounds(real** bounds, int** index, int** dr_pair, int* npairs, gmx_localtop_t* top) +static int get_bounds(real** bounds, int** index, int** dr_pair, int* npairs, const InteractionDefinitions& idef) { int i, j, k, type, ftype, natom; real* b; int * ind, *pair; int nb, label1; - gmx::ArrayRef functype = top->idef.functype; - gmx::ArrayRef iparams = top->idef.iparams; + gmx::ArrayRef functype = idef.functype; + gmx::ArrayRef iparams = idef.iparams; /* Count how many distance restraint there are... */ - nb = top->idef.il[F_DISRES].size(); + nb = idef.il[F_DISRES].size(); if (nb == 0) { gmx_fatal(FARGS, "No distance restraints in topology!\n"); @@ -238,7 +238,7 @@ static int get_bounds(real** bounds, int** index, int** dr_pair, int* npairs, gm /* Fill the index array */ label1 = -1; - const InteractionList& disres = top->idef.il[F_DISRES]; + const InteractionList& disres = idef.il[F_DISRES]; gmx::ArrayRef iatom = disres.iatoms; for (i = j = k = 0; (i < disres.size());) { @@ -617,7 +617,7 @@ int gmx_nmr(int argc, char* argv[]) top = std::make_unique(topInfo.mtop()->ffparams); gmx_mtop_generate_local_top(*topInfo.mtop(), top.get(), ir->efep != efepNO); } - nbounds = get_bounds(&bounds, &index, &pair, &npairs, top.get()); + nbounds = get_bounds(&bounds, &index, &pair, &npairs, top->idef); snew(violaver, npairs); out_disre = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv); xvgr_legend(out_disre, 2, drleg, oenv); diff --git a/src/gromacs/listed_forces/bonded.cpp b/src/gromacs/listed_forces/bonded.cpp index 75fd2c8e8f..60b42ca235 100644 --- a/src/gromacs/listed_forces/bonded.cpp +++ b/src/gromacs/listed_forces/bonded.cpp @@ -3542,14 +3542,14 @@ real bonded_tab(const char* type, real* V, real* F) { - real k, tabscale, *VFtab, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF; + real k, tabscale, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF; int n0, nnn; real dvdlambda; k = (1.0 - lambda) * kA + lambda * kB; - tabscale = table->scale; - VFtab = table->data; + tabscale = table->scale; + const real* VFtab = table->data.data(); rt = r * tabscale; n0 = static_cast(rt); diff --git a/src/gromacs/listed_forces/disre.cpp b/src/gromacs/listed_forces/disre.cpp index 5025cc010b..1c38611eb3 100644 --- a/src/gromacs/listed_forces/disre.cpp +++ b/src/gromacs/listed_forces/disre.cpp @@ -75,19 +75,16 @@ void init_disres(FILE* fplog, NumRanks numRanks, MPI_Comm communicator, const gmx_multisim_t* ms, - t_fcdata* fcd, + t_disresdata* dd, t_state* state, gmx_bool bIsREMD) { int fa, nmol, npair, np; - t_disresdata* dd; history_t* hist; gmx_mtop_ilistloop_t iloop; char* ptr; int type_min, type_max; - dd = &(fcd->disres); - if (gmx_mtop_ftype_count(mtop, F_DISRES) == 0) { dd->nres = 0; @@ -241,7 +238,7 @@ void init_disres(FILE* fplog, /* We use to allow any value of nsystems which was a divisor * of ms->nsim. But this required an extra communicator which - * was stored in t_fcdata. This pulled in mpi.h in nearly all C files. + * pulled in mpi.h in nearly all C files. */ if (!(ms->nsim == 1 || ms->nsim == dd->nsystems)) { @@ -292,7 +289,7 @@ void init_disres(FILE* fplog, * succeed...) */ if ((disResRunMode == DisResRunMode::MDRun) && ms && dd->nsystems > 1 && (ddRole == DDRole::Master)) { - check_multi_int(fplog, ms, fcd->disres.nres, "the number of distance restraints", FALSE); + check_multi_int(fplog, ms, dd->nres, "the number of distance restraints", FALSE); } please_cite(fplog, "Tropp80a"); please_cite(fplog, "Torda89a"); @@ -305,16 +302,14 @@ void calc_disres_R_6(const t_commrec* cr, const t_iatom forceatoms[], const rvec x[], const t_pbc* pbc, - t_fcdata* fcd, + t_disresdata* dd, history_t* hist) { - rvec dx; - real * rt, *rm3tav, *Rtl_6, *Rt_6, *Rtav_6; - t_disresdata* dd; - real ETerm, ETerm1, cf1 = 0, cf2 = 0; - gmx_bool bTav; + rvec dx; + real * rt, *rm3tav, *Rtl_6, *Rt_6, *Rtav_6; + real ETerm, ETerm1, cf1 = 0, cf2 = 0; + gmx_bool bTav; - dd = &(fcd->disres); bTav = (dd->dr_tau != 0); ETerm = dd->ETerm; ETerm1 = dd->ETerm1; @@ -365,7 +360,7 @@ void calc_disres_R_6(const t_commrec* cr, rt[pair] = rt2 * rt_1; if (bTav) { - /* Here we update rm3tav in t_fcdata using the data + /* Here we update rm3tav in t_disresdata using the data * in history_t. * Thus the results stay correct when this routine * is called multiple times. @@ -391,7 +386,7 @@ void calc_disres_R_6(const t_commrec* cr, gmx_sum(2 * dd->nres, dd->Rt_6, cr); } - if (fcd->disres.nsystems > 1) + if (dd->nsystems > 1) { real invn = 1.0 / dd->nsystems; @@ -447,7 +442,7 @@ real ta_disres(int nfa, int dr_weighting; gmx_bool dr_bMixed; - dd = &(fcd->disres); + dd = fcd->disres; dr_weighting = dd->dr_weighting; dr_bMixed = dd->dr_bMixed; Rtl_6 = dd->Rtl_6; @@ -643,21 +638,17 @@ real ta_disres(int nfa, return vtot; } -void update_disres_history(const t_fcdata* fcd, history_t* hist) +void update_disres_history(const t_disresdata& dd, history_t* hist) { - const t_disresdata* dd; - int pair; - - dd = &(fcd->disres); - if (dd->dr_tau != 0) + if (dd.dr_tau != 0) { /* Copy the new time averages that have been calculated * in calc_disres_R_6. */ - hist->disre_initf = dd->exp_min_t_tau; - for (pair = 0; pair < dd->npair; pair++) + hist->disre_initf = dd.exp_min_t_tau; + for (int pair = 0; pair < dd.npair; pair++) { - hist->disre_rm3tav[pair] = dd->rm3tav[pair]; + hist->disre_rm3tav[pair] = dd.rm3tav[pair]; } } } diff --git a/src/gromacs/listed_forces/disre.h b/src/gromacs/listed_forces/disre.h index 43983bfb82..7144dd0f12 100644 --- a/src/gromacs/listed_forces/disre.h +++ b/src/gromacs/listed_forces/disre.h @@ -55,6 +55,8 @@ struct gmx_mtop_t; struct gmx_multisim_t; class history_t; struct t_commrec; +struct t_disresdata; +struct t_fcdata; struct t_inputrec; struct t_pbc; class t_state; @@ -69,7 +71,7 @@ enum class DisResRunMode }; /*! \brief - * Initiates *fcd data. + * Initiates *disresdata. * * Must be called once, nbonds is the number * of iatoms in the ilist of the idef struct. @@ -86,7 +88,7 @@ void init_disres(FILE* fplog, NumRanks numRanks, MPI_Comm communicator, const gmx_multisim_t* ms, - t_fcdata* fcd, + t_disresdata* disresdata, t_state* state, gmx_bool bIsREMD); @@ -100,7 +102,7 @@ void calc_disres_R_6(const t_commrec* cr, const t_iatom* fa, const rvec* x, const t_pbc* pbc, - t_fcdata* fcd, + t_disresdata* disresdata, history_t* hist); //! Calculates the distance restraint forces, return the potential. @@ -114,10 +116,10 @@ real ta_disres(int nfa, real lambda, real* dvdlambda, const t_mdatoms* md, - t_fcdata* fcd, + t_fcdata* fcdata, int* global_atom_index); //! Copies the new time averages that have been calculated in calc_disres_R_6. -void update_disres_history(const t_fcdata* fcd, history_t* hist); +void update_disres_history(const t_disresdata& disresdata, history_t* hist); #endif diff --git a/src/gromacs/listed_forces/listed_forces.cpp b/src/gromacs/listed_forces/listed_forces.cpp index 10738a759c..ea2c9a0a9e 100644 --- a/src/gromacs/listed_forces/listed_forces.cpp +++ b/src/gromacs/listed_forces/listed_forces.cpp @@ -61,6 +61,7 @@ #include "gromacs/math/vec.h" #include "gromacs/mdlib/enerdata_utils.h" #include "gromacs/mdlib/force.h" +#include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/fcdata.h" #include "gromacs/mdtypes/forceoutput.h" @@ -74,11 +75,32 @@ #include "gromacs/topology/topology.h" #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/smalloc.h" #include "listed_internal.h" +#include "manage_threading.h" #include "utilities.h" +ListedForces::ListedForces(const int numEnergyGroups, const int numThreads, FILE* fplog) : + threading_(std::make_unique(numThreads, numEnergyGroups, fplog)), + fcdata_(std::make_unique()) +{ +} + +ListedForces::~ListedForces() = default; + +void ListedForces::setup(const InteractionDefinitions& idef, const int numAtomsForce, const bool useGpu) +{ + idef_ = &idef; + + setup_bonded_threading(threading_.get(), numAtomsForce, useGpu, *idef_); + + if (idef_->ilsort == ilsortFE_SORTED) + { + forceBufferLambda_.resize(numAtomsForce * sizeof(rvec4) / sizeof(real)); + shiftForceBufferLambda_.resize(SHIFTS); + } +} + namespace { @@ -137,7 +159,7 @@ void zero_thread_output(f_thread_t* f_t) #define MAX_BONDED_THREADS 256 /*! \brief Reduce thread-local force buffers */ -void reduce_thread_forces(int n, gmx::ArrayRef force, const bonded_threading_t* bt, int nthreads) +void reduce_thread_forces(gmx::ArrayRef force, const bonded_threading_t* bt, int nthreads) { if (nthreads > MAX_BONDED_THREADS) { @@ -146,6 +168,8 @@ void reduce_thread_forces(int n, gmx::ArrayRef force, const bonded_th rvec* gmx_restrict f = as_rvec_array(force.data()); + const int numAtomsForce = bt->numAtomsForce; + /* This reduction can run on any number of threads, * independently of bt->nthreads. * But if nthreads matches bt->nthreads (which it currently does) @@ -177,7 +201,7 @@ void reduce_thread_forces(int n, gmx::ArrayRef force, const bonded_th int a0 = ind * reduction_block_size; int a1 = (ind + 1) * reduction_block_size; /* It would be nice if we could pad f to avoid this min */ - a1 = std::min(a1, n); + a1 = std::min(a1, numAtomsForce); for (int a = a0; a < a1; a++) { for (int fb = 0; fb < nfb; fb++) @@ -192,8 +216,7 @@ void reduce_thread_forces(int n, gmx::ArrayRef force, const bonded_th } /*! \brief Reduce thread-local forces, shift forces and energies */ -void reduce_thread_output(int n, - gmx::ForceWithShiftForces* forceWithShiftForces, +void reduce_thread_output(gmx::ForceWithShiftForces* forceWithShiftForces, real* ener, gmx_grppairener_t* grpp, real* dvdl, @@ -205,7 +228,7 @@ void reduce_thread_output(int n, if (bt->nblock_used > 0) { /* Reduce the bonded force buffer */ - reduce_thread_forces(n, forceWithShiftForces->force(), bt, bt->nthreads); + reduce_thread_forces(forceWithShiftForces->force(), bt, bt->nthreads); } rvec* gmx_restrict fshift = as_rvec_array(forceWithShiftForces->shiftForces().data()); @@ -388,6 +411,7 @@ real calc_one_bond(int thread, /*! \brief Compute the bonded part of the listed forces, parallelized over threads */ static void calcBondedForces(const InteractionDefinitions& idef, + bonded_threading_t* bt, const rvec x[], const t_forcerec* fr, const t_pbc* pbc_null, @@ -401,8 +425,6 @@ static void calcBondedForces(const InteractionDefinitions& idef, const gmx::StepWorkload& stepWork, int* global_atom_index) { - bonded_threading_t* bt = fr->bondedThreading; - #pragma omp parallel for num_threads(bt->nthreads) schedule(static) for (int thread = 0; thread < bt->nthreads; thread++) { @@ -432,7 +454,7 @@ static void calcBondedForces(const InteractionDefinitions& idef, } else { - fshift = threadBuffers.fshift; + fshift = as_rvec_array(threadBuffers.fshift.data()); epot = threadBuffers.ener; grpp = &threadBuffers.grpp; dvdlt = threadBuffers.dvdl; @@ -445,8 +467,8 @@ static void calcBondedForces(const InteractionDefinitions& idef, { ArrayRef iatoms = gmx::makeConstArrayRef(ilist.iatoms); v = calc_one_bond(thread, ftype, idef, iatoms, idef.numNonperturbedInteractions[ftype], - fr->bondedThreading->workDivision, x, ft, fshift, fr, pbc_null, - grpp, nrnb, lambda, dvdlt, md, fcd, stepWork, global_atom_index); + bt->workDivision, x, ft, fshift, fr, pbc_null, grpp, nrnb, + lambda, dvdlt, md, fcd, stepWork, global_atom_index); epot[ftype] += v; } } @@ -455,20 +477,23 @@ static void calcBondedForces(const InteractionDefinitions& idef, } } -bool haveRestraints(const InteractionDefinitions& idef, const t_fcdata& fcd) +bool ListedForces::haveRestraints() const { - return (!idef.il[F_POSRES].empty() || !idef.il[F_FBPOSRES].empty() || fcd.orires.nr > 0 - || fcd.disres.nres > 0); + GMX_ASSERT(fcdata_, "Need valid fcdata"); + GMX_ASSERT(fcdata_->orires && fcdata_->disres, "NMR restraints objects should be set up"); + + return (!idef_->il[F_POSRES].empty() || !idef_->il[F_FBPOSRES].empty() + || fcdata_->orires->nr > 0 || fcdata_->disres->nres > 0); } -bool haveCpuBondeds(const t_forcerec& fr) +bool ListedForces::haveCpuBondeds() const { - return fr.bondedThreading->haveBondeds; + return threading_->haveBondeds; } -bool haveCpuListedForces(const t_forcerec& fr, const InteractionDefinitions& idef, const t_fcdata& fcd) +bool ListedForces::haveCpuListedForces() const { - return haveCpuBondeds(fr) || haveRestraints(idef, fcd); + return haveCpuBondeds() || haveRestraints(); } namespace @@ -477,6 +502,7 @@ namespace /*! \brief Calculates all listed force interactions. */ void calc_listed(struct gmx_wallcycle* wcycle, const InteractionDefinitions& idef, + bonded_threading_t* bt, const rvec x[], gmx::ForceOutputs* forceOutputs, const t_forcerec* fr, @@ -489,9 +515,7 @@ void calc_listed(struct gmx_wallcycle* wcycle, int* global_atom_index, const gmx::StepWorkload& stepWork) { - bonded_threading_t* bt = fr->bondedThreading; - - if (haveCpuBondeds(*fr)) + if (bt->haveBondeds) { gmx::ForceWithShiftForces& forceWithShiftForces = forceOutputs->forceWithShiftForces(); @@ -499,14 +523,13 @@ void calc_listed(struct gmx_wallcycle* wcycle, /* The dummy array is to have a place to store the dhdl at other values of lambda, which will be thrown away in the end */ real dvdl[efptNR] = { 0 }; - calcBondedForces(idef, x, fr, fr->bMolPBC ? pbc : nullptr, + calcBondedForces(idef, bt, x, fr, fr->bMolPBC ? pbc : nullptr, as_rvec_array(forceWithShiftForces.shiftForces().data()), enerd, nrnb, lambda, dvdl, md, fcd, stepWork, global_atom_index); wallcycle_sub_stop(wcycle, ewcsLISTED); wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS); - reduce_thread_output(fr->natoms_force, &forceWithShiftForces, enerd->term, &enerd->grpp, - dvdl, bt, stepWork); + reduce_thread_output(&forceWithShiftForces, enerd->term, &enerd->grpp, dvdl, bt, stepWork); if (stepWork.computeDhdl) { @@ -521,7 +544,7 @@ void calc_listed(struct gmx_wallcycle* wcycle, /* Copy the sum of violations for the distance restraints from fcd */ if (fcd) { - enerd->term[F_DISRESVIOL] = fcd->disres.sumviol; + enerd->term[F_DISRESVIOL] = fcd->disres->sumviol; } } @@ -531,9 +554,12 @@ void calc_listed(struct gmx_wallcycle* wcycle, * The shift forces in fr are not affected. */ void calc_listed_lambda(const InteractionDefinitions& idef, + bonded_threading_t* bt, const rvec x[], const t_forcerec* fr, const struct t_pbc* pbc, + gmx::ArrayRef forceBufferLambda, + gmx::ArrayRef shiftForceBufferLambda, gmx_grppairener_t* grpp, real* epot, gmx::ArrayRef dvdl, @@ -543,12 +569,9 @@ void calc_listed_lambda(const InteractionDefinitions& idef, t_fcdata* fcd, int* global_atom_index) { - real v; - rvec4* f; - rvec* fshift; - const t_pbc* pbc_null; - WorkDivision& workDivision = fr->bondedThreading->foreignLambdaWorkDivision; + WorkDivision& workDivision = bt->foreignLambdaWorkDivision; + const t_pbc* pbc_null; if (fr->bMolPBC) { pbc_null = pbc; @@ -559,9 +582,11 @@ void calc_listed_lambda(const InteractionDefinitions& idef, } /* We already have the forces, so we use temp buffers here */ - // TODO: Get rid of these allocations by using permanent force buffers - snew(f, fr->natoms_force); - snew(fshift, SHIFTS); + std::fill(forceBufferLambda.begin(), forceBufferLambda.end(), 0.0_real); + std::fill(shiftForceBufferLambda.begin(), shiftForceBufferLambda.end(), + gmx::RVec{ 0.0_real, 0.0_real, 0.0_real }); + rvec4* f = reinterpret_cast(forceBufferLambda.data()); + rvec* fshift = as_rvec_array(shiftForceBufferLambda.data()); /* Loop over all bonded force types to calculate the bonded energies */ for (int ftype = 0; (ftype < F_NRE); ftype++) @@ -581,47 +606,45 @@ void calc_listed_lambda(const InteractionDefinitions& idef, gmx::StepWorkload tempFlags; tempFlags.computeEnergy = true; - v = calc_one_bond(0, ftype, idef, iatomsPerturbed, iatomsPerturbed.ssize(), - workDivision, x, f, fshift, fr, pbc_null, grpp, nrnb, lambda, - dvdl.data(), md, fcd, tempFlags, global_atom_index); + real v = calc_one_bond(0, ftype, idef, iatomsPerturbed, iatomsPerturbed.ssize(), + workDivision, x, f, fshift, fr, pbc_null, grpp, nrnb, lambda, + dvdl.data(), md, fcd, tempFlags, global_atom_index); epot[ftype] += v; } } } - - sfree(fshift); - sfree(f); } } // namespace -void do_force_listed(struct gmx_wallcycle* wcycle, - const matrix box, - const t_lambda* fepvals, - const t_commrec* cr, - const gmx_multisim_t* ms, - const InteractionDefinitions& idef, - const rvec x[], - gmx::ArrayRef xWholeMolecules, - history_t* hist, - gmx::ForceOutputs* forceOutputs, - const t_forcerec* fr, - const struct t_pbc* pbc, - gmx_enerdata_t* enerd, - t_nrnb* nrnb, - const real* lambda, - const t_mdatoms* md, - t_fcdata* fcd, - int* global_atom_index, - const gmx::StepWorkload& stepWork) +void ListedForces::calculate(struct gmx_wallcycle* wcycle, + const matrix box, + const t_lambda* fepvals, + const t_commrec* cr, + const gmx_multisim_t* ms, + const rvec x[], + gmx::ArrayRef xWholeMolecules, + history_t* hist, + gmx::ForceOutputs* forceOutputs, + const t_forcerec* fr, + const struct t_pbc* pbc, + gmx_enerdata_t* enerd, + t_nrnb* nrnb, + const real* lambda, + const t_mdatoms* md, + int* global_atom_index, + const gmx::StepWorkload& stepWork) { if (!stepWork.computeListedForces) { return; } + const InteractionDefinitions& idef = *idef_; + t_fcdata& fcdata = *fcdata_; + t_pbc pbc_full; /* Full PBC is needed for position restraints */ - if (haveRestraints(idef, *fcd)) + if (haveRestraints()) { if (!idef.il[F_POSRES].empty() || !idef.il[F_FBPOSRES].empty()) { @@ -647,24 +670,24 @@ void do_force_listed(struct gmx_wallcycle* wcycle, } /* Do pre force calculation stuff which might require communication */ - if (fcd->orires.nr > 0) + if (fcdata.orires->nr > 0) { GMX_ASSERT(!xWholeMolecules.empty(), "Need whole molecules for orienation restraints"); enerd->term[F_ORIRESDEV] = calc_orires_dev( ms, idef.il[F_ORIRES].size(), idef.il[F_ORIRES].iatoms.data(), idef.iparams.data(), - md, xWholeMolecules, x, fr->bMolPBC ? pbc : nullptr, fcd, hist); + md, xWholeMolecules, x, fr->bMolPBC ? pbc : nullptr, fcdata.orires, hist); } - if (fcd->disres.nres > 0) + if (fcdata.disres->nres > 0) { calc_disres_R_6(cr, ms, idef.il[F_DISRES].size(), idef.il[F_DISRES].iatoms.data(), x, - fr->bMolPBC ? pbc : nullptr, fcd, hist); + fr->bMolPBC ? pbc : nullptr, fcdata.disres, hist); } wallcycle_sub_stop(wcycle, ewcsRESTRAINTS); } - calc_listed(wcycle, idef, x, forceOutputs, fr, pbc, enerd, nrnb, lambda, md, fcd, - global_atom_index, stepWork); + calc_listed(wcycle, idef, threading_.get(), x, forceOutputs, fr, pbc, enerd, nrnb, lambda, md, + &fcdata, global_atom_index, stepWork); /* Check if we have to determine energy differences * at foreign lambda's. @@ -692,8 +715,9 @@ void do_force_listed(struct gmx_wallcycle* wcycle, { lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i - 1]); } - calc_listed_lambda(idef, x, fr, pbc, &(enerd->foreign_grpp), enerd->foreign_term, - dvdl, nrnb, lam_i, md, fcd, global_atom_index); + calc_listed_lambda(idef, threading_.get(), x, fr, pbc, forceBufferLambda_, + shiftForceBufferLambda_, &(enerd->foreign_grpp), enerd->foreign_term, + dvdl, nrnb, lam_i, md, &fcdata, global_atom_index); sum_epot(&(enerd->foreign_grpp), enerd->foreign_term); enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT]; for (int j = 0; j < efptNR; j++) diff --git a/src/gromacs/listed_forces/listed_forces.h b/src/gromacs/listed_forces/listed_forces.h index 67c078430b..d821973fc3 100644 --- a/src/gromacs/listed_forces/listed_forces.h +++ b/src/gromacs/listed_forces/listed_forces.h @@ -49,6 +49,7 @@ and reduction of output data across threads. * * \author Mark Abraham + * \author Berk Hess * */ /*! \libinternal \file @@ -60,6 +61,7 @@ * should call functions declared here. * * \author Mark Abraham + * \author Berk Hess * * \inlibraryapi * \ingroup module_listed_forces @@ -67,12 +69,19 @@ #ifndef GMX_LISTED_FORCES_LISTED_FORCES_H #define GMX_LISTED_FORCES_LISTED_FORCES_H +#include + #include "gromacs/math/vectypes.h" #include "gromacs/topology/ifunc.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/classhelpers.h" +struct bonded_threading_t; struct gmx_enerdata_t; +struct gmx_ffparams_t; struct gmx_grppairener_t; +struct gmx_localtop_t; struct gmx_multisim_t; class history_t; class InteractionDefinitions; @@ -109,43 +118,87 @@ using BondedFunction = real (*)(int nbonds, //! Getter for finding a callable CPU function to compute an \c ftype interaction. BondedFunction bondedFunction(int ftype); -/*! \brief Do all aspects of energy and force calculations for mdrun - * on the set of listed interactions - * - * xWholeMolecules only needs to contain whole molecules when orientation - * restraints need to be computed and can be empty otherwise. - */ -void do_force_listed(struct gmx_wallcycle* wcycle, - const matrix box, - const t_lambda* fepvals, - const t_commrec* cr, - const gmx_multisim_t* ms, - const InteractionDefinitions& idef, - const rvec x[], - gmx::ArrayRef xWholeMolecules, - history_t* hist, - gmx::ForceOutputs* forceOutputs, - const t_forcerec* fr, - const struct t_pbc* pbc, - gmx_enerdata_t* enerd, - t_nrnb* nrnb, - const real* lambda, - const t_mdatoms* md, - struct t_fcdata* fcd, - int* global_atom_index, - const gmx::StepWorkload& stepWork); - -/*! \brief Returns true if there are position, distance or orientation restraints. */ -bool haveRestraints(const InteractionDefinitions& idef, const t_fcdata& fcd); - -/*! \brief Returns true if there are CPU (i.e. not GPU-offloaded) bonded interactions to compute. */ -bool haveCpuBondeds(const t_forcerec& fr); - -/*! \brief Returns true if there are listed interactions to compute. - * - * NOTE: the current implementation returns true if there are position restraints - * or any bonded interactions computed on the CPU. +/*! \libinternal + * \brief Class for calculating listed interactions, uses OpenMP parallelization */ -bool haveCpuListedForces(const t_forcerec& fr, const InteractionDefinitions& idef, const t_fcdata& fcd); +class ListedForces +{ +public: + /*! \brief Constructor + * + * \param[in] numEnergyGroups The number of energy groups, used for storage of pair energies + * \param[in] numThreads The number of threads used for computed listed interactions + * \param[in] fplog Log file for printing env.var. override, can be nullptr + */ + ListedForces(int numEnergyGroups, int numThreads, FILE* fplog); + + //! Destructor which is actually default but in the source file to hide implementation classes + ~ListedForces(); + + /*! \brief Copy the listed interactions from \p idef and set up the thread parallelization + * + * \param[in] idef The idef with all listed interactions to be computed on this rank + * \param[in] numAtomsForce Force are, potentially, computed for atoms 0 to \p numAtomsForce + * \param[in] useGpu Whether a GPU is used to compute (part of) the listed interactions + */ + void setup(const InteractionDefinitions& idef, int numAtomsForce, bool useGpu); + + /*! \brief Do all aspects of energy and force calculations for mdrun + * on the set of listed interactions + * + * xWholeMolecules only needs to contain whole molecules when orientation + * restraints need to be computed and can be empty otherwise. + */ + void calculate(struct gmx_wallcycle* wcycle, + const matrix box, + const t_lambda* fepvals, + const t_commrec* cr, + const gmx_multisim_t* ms, + const rvec x[], + gmx::ArrayRef xWholeMolecules, + history_t* hist, + gmx::ForceOutputs* forceOutputs, + const t_forcerec* fr, + const struct t_pbc* pbc, + gmx_enerdata_t* enerd, + t_nrnb* nrnb, + const real* lambda, + const t_mdatoms* md, + int* global_atom_index, + const gmx::StepWorkload& stepWork); + + //! Returns whether bonded interactions are assigned to the CPU + bool haveCpuBondeds() const; + + /*! \brief Returns whether listed forces are computed on the CPU + * + * NOTE: the current implementation returns true if there are position restraints + * or any bonded interactions computed on the CPU. + */ + bool haveCpuListedForces() const; + + //! Returns true if there are position, distance or orientation restraints + bool haveRestraints() const; + + //! Returns a reference to the force calculation data + const t_fcdata& fcdata() const { return *fcdata_; } + + //! Returns a reference to the force calculation data + t_fcdata& fcdata() { return *fcdata_; } + +private: + //! The interaction definitions + InteractionDefinitions const* idef_ = nullptr; + //! Thread parallelization setup, unique_ptr to avoid declaring bonded_threading_t + std::unique_ptr threading_; + //! Data for bonded tables and NMR restraining, needs to be refactored + std::unique_ptr fcdata_; + //! Force buffer for free-energy forces + std::vector forceBufferLambda_; + //! Shift force buffer for free-energy forces + std::vector shiftForceBufferLambda_; + + GMX_DISALLOW_COPY_MOVE_AND_ASSIGN(ListedForces); +}; #endif diff --git a/src/gromacs/listed_forces/listed_internal.h b/src/gromacs/listed_forces/listed_internal.h index a9e698deb8..2e130dcb83 100644 --- a/src/gromacs/listed_forces/listed_internal.h +++ b/src/gromacs/listed_forces/listed_internal.h @@ -1,7 +1,8 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2014,2015,2016,2018,2019, by the GROMACS development team, led by + * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team. + * Copyright (c) 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. @@ -38,6 +39,7 @@ * internally by the module. * * \author Mark Abraham + * \author Berk Hess * \ingroup module_listed_forces */ #ifndef GMX_LISTED_FORCES_LISTED_INTERNAL_H @@ -49,7 +51,9 @@ #include "gromacs/mdtypes/enerdata.h" #include "gromacs/topology/idef.h" #include "gromacs/topology/ifunc.h" +#include "gromacs/utility/alignedallocator.h" #include "gromacs/utility/bitmask.h" +#include "gromacs/utility/classhelpers.h" /* We reduce the force array in blocks of 32 atoms. This is large enough * to not cause overhead and 32*sizeof(rvec) is a multiple of the cache-line @@ -93,52 +97,65 @@ struct f_thread_t //! Constructor f_thread_t(int numEnergyGroups); - ~f_thread_t(); + ~f_thread_t() = default; - rvec4* f = nullptr; /**< Force array */ - int f_nalloc = 0; /**< Allocation size of f */ - gmx_bitmask_t* mask = - nullptr; /**< Mask for marking which parts of f are filled, working array for constructing mask in bonded_threading_t */ - int nblock_used; /**< Number of blocks touched by our thread */ - int* block_index = nullptr; /**< Index to touched blocks, size nblock_used */ - int block_nalloc = 0; /**< Allocation size of f (*reduction_block_size), mask_index, mask */ + //! Force array pointer, equals fBuffer.data(), needed because rvec4 is not a C++ type + rvec4* f = nullptr; + //! Force array buffer + std::vector> fBuffer; + //! Mask for marking which parts of f are filled, working array for constructing mask in bonded_threading_t + std::vector mask; + //! Number of blocks touched by our thread + int nblock_used = 0; + //! Index to touched blocks + std::vector block_index; - rvec* fshift; /**< Shift force array, size SHIFTS */ - real ener[F_NRE]; /**< Energy array */ - gmx_grppairener_t grpp; /**< Group pair energy data for pairs */ - real dvdl[efptNR]; /**< Free-energy dV/dl output */ + //! Shift force array, size SHIFTS + std::vector fshift; + //! Energy array + real ener[F_NRE]; + //! Group pair energy data for pairs + gmx_grppairener_t grpp; + //! Free-energy dV/dl output + real dvdl[efptNR]; + + GMX_DISALLOW_COPY_MOVE_AND_ASSIGN(f_thread_t); }; /*! \internal \brief struct contain all data for bonded force threading */ struct bonded_threading_t { //! Constructor - bonded_threading_t(int numThreads, int numEnergyGroups); + bonded_threading_t(int numThreads, int numEnergyGroups, FILE* fplog); //! Number of threads to be used for bondeds - int nthreads; + int nthreads = 0; //! Force/energy data per thread, size nthreads, stored in unique_ptr to allow thread local allocation std::vector> f_t; //! The number of force blocks to reduce - int nblock_used; + int nblock_used = 0; //! Index of size nblock_used into mask std::vector block_index; //! Mask array, one element corresponds to a block of reduction_block_size atoms of the force array, bit corresponding to thread indices set if a thread writes to that block std::vector mask; //! true if we have and thus need to reduce bonded forces - bool haveBondeds; + bool haveBondeds = false; + //! The number of atoms forces are computed for + int numAtomsForce = 0; /* There are two different ways to distribute the bonded force calculation * over the threads. We dedice which to use based on the number of threads. */ //! Maximum thread count for uniform distribution of bondeds over threads - int max_nthread_uniform; + int max_nthread_uniform = 0; //! The division of work in the t_list over threads. WorkDivision workDivision; //! Work division for free-energy foreign lambda calculations, always uses 1 thread WorkDivision foreignLambdaWorkDivision; + + GMX_DISALLOW_COPY_MOVE_AND_ASSIGN(bonded_threading_t); }; diff --git a/src/gromacs/listed_forces/manage_threading.cpp b/src/gromacs/listed_forces/manage_threading.cpp index f663cbebf3..f666b5a012 100644 --- a/src/gromacs/listed_forces/manage_threading.cpp +++ b/src/gromacs/listed_forces/manage_threading.cpp @@ -40,6 +40,7 @@ * interactions. * * \author Mark Abraham + * \author Berk Hess * \ingroup module_listed_forces */ #include "gmxpre.h" @@ -57,13 +58,11 @@ #include #include "gromacs/listed_forces/gpubonded.h" -#include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/pbcutil/ishift.h" #include "gromacs/topology/ifunc.h" #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/gmxassert.h" -#include "gromacs/utility/smalloc.h" #include "listed_internal.h" #include "utilities.h" @@ -348,25 +347,21 @@ static void calc_bonded_reduction_mask(int natoms, GMX_ASSERT(bondedThreading.nthreads <= BITMASK_SIZE, "We need at least nthreads bits in the mask"); - int nblock = (natoms + reduction_block_size - 1) >> reduction_block_bits; + const int nblock = (natoms + reduction_block_size - 1) >> reduction_block_bits; - if (nblock > f_thread->block_nalloc) - { - f_thread->block_nalloc = over_alloc_large(nblock); - srenew(f_thread->mask, f_thread->block_nalloc); - srenew(f_thread->block_index, f_thread->block_nalloc); - // NOTE: It seems f_thread->f does not need to be aligned - sfree_aligned(f_thread->f); - snew_aligned(f_thread->f, f_thread->block_nalloc * reduction_block_size, 128); - } - - gmx_bitmask_t* mask = f_thread->mask; + f_thread->mask.resize(nblock); + f_thread->block_index.resize(nblock); + // NOTE: It seems f_thread->f does not need to be aligned + f_thread->fBuffer.resize(nblock * reduction_block_size * sizeof(rvec4) / sizeof(real)); + f_thread->f = reinterpret_cast(f_thread->fBuffer.data()); - for (int b = 0; b < nblock; b++) + for (gmx_bitmask_t& mask : f_thread->mask) { - bitmask_clear(&mask[b]); + bitmask_clear(&mask); } + gmx::ArrayRef mask = f_thread->mask; + for (int ftype = 0; ftype < F_NRE; ftype++) { if (ftype_is_bonded_potential(ftype)) @@ -404,7 +399,7 @@ static void calc_bonded_reduction_mask(int natoms, } void setup_bonded_threading(bonded_threading_t* bt, - int numAtoms, + int numAtomsForce, bool useGpuForBondeds, const InteractionDefinitions& idef) { @@ -412,6 +407,8 @@ void setup_bonded_threading(bonded_threading_t* bt, assert(bt->nthreads >= 1); + bt->numAtomsForce = numAtomsForce; + /* Divide the bonded interaction over the threads */ divide_bondeds_over_threads(bt, useGpuForBondeds, idef); @@ -429,7 +426,7 @@ void setup_bonded_threading(bonded_threading_t* bt, { try { - calc_bonded_reduction_mask(numAtoms, bt->f_t[t].get(), idef, t, *bt); + calc_bonded_reduction_mask(numAtomsForce, bt->f_t[t].get(), idef, t, *bt); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } @@ -437,7 +434,7 @@ void setup_bonded_threading(bonded_threading_t* bt, /* Reduce the masks over the threads and determine which blocks * we need to reduce over. */ - int nblock_tot = (numAtoms + reduction_block_size - 1) >> reduction_block_bits; + int nblock_tot = (numAtomsForce + reduction_block_size - 1) >> reduction_block_bits; /* Ensure we have sufficient space for all blocks */ if (static_cast(nblock_tot) > bt->block_index.size()) { @@ -485,36 +482,29 @@ void setup_bonded_threading(bonded_threading_t* bt, { fprintf(debug, "Number of %d atom blocks to reduce: %d\n", reduction_block_size, bt->nblock_used); fprintf(debug, "Reduction density %.2f for touched blocks only %.2f\n", - ctot * reduction_block_size / static_cast(numAtoms), + ctot * reduction_block_size / static_cast(numAtomsForce), ctot / static_cast(bt->nblock_used)); } } -void tear_down_bonded_threading(bonded_threading_t* bt) -{ - delete bt; -} - -f_thread_t::f_thread_t(int numEnergyGroups) : grpp(numEnergyGroups) -{ - snew(fshift, SHIFTS); -} - -f_thread_t::~f_thread_t() -{ - sfree(mask); - sfree(fshift); - sfree(block_index); - sfree_aligned(f); -} +f_thread_t::f_thread_t(int numEnergyGroups) : fshift(SHIFTS), grpp(numEnergyGroups) {} -bonded_threading_t::bonded_threading_t(const int numThreads, const int numEnergyGroups) : +bonded_threading_t::bonded_threading_t(const int numThreads, const int numEnergyGroups, FILE* fplog) : nthreads(numThreads), nblock_used(0), haveBondeds(false), workDivision(nthreads), foreignLambdaWorkDivision(1) { + /* These thread local data structures are used for bondeds only. + * + * Note that we also use there structures when running single-threaded. + * This is because the bonded force buffer uses type rvec4, whereas + * the normal force buffer is uses type rvec. This leads to a little + * reduction overhead, but the speed gain in the bonded calculations + * of doing transposeScatterIncr/DecrU with aligment 4 instead of 3 + * is much larger than the reduction overhead. + */ f_t.resize(numThreads); #pragma omp parallel for num_threads(nthreads) schedule(static) for (int t = 0; t < nthreads; t++) @@ -528,41 +518,25 @@ bonded_threading_t::bonded_threading_t(const int numThreads, const int numEnergy } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } -} - -bonded_threading_t* init_bonded_threading(FILE* fplog, const int nenergrp) -{ - /* These thread local data structures are used for bondeds only. - * - * Note that we also use there structures when running single-threaded. - * This is because the bonded force buffer uses type rvec4, whereas - * the normal force buffer is uses type rvec. This leads to a little - * reduction overhead, but the speed gain in the bonded calculations - * of doing transposeScatterIncr/DecrU with aligment 4 instead of 3 - * is much larger than the reduction overhead. - */ - bonded_threading_t* bt = new bonded_threading_t(gmx_omp_nthreads_get(emntBonded), nenergrp); /* The optimal value after which to switch from uniform to localized * bonded interaction distribution is 3, 4 or 5 depending on the system * and hardware. */ - const int max_nthread_uniform = 4; + const int max_nthread_uniform_default = 4; char* ptr; if ((ptr = getenv("GMX_BONDED_NTHREAD_UNIFORM")) != nullptr) { - sscanf(ptr, "%d", &bt->max_nthread_uniform); + sscanf(ptr, "%d", &max_nthread_uniform); if (fplog != nullptr) { fprintf(fplog, "\nMax threads for uniform bonded distribution set to %d by env.var.\n", - bt->max_nthread_uniform); + max_nthread_uniform); } } else { - bt->max_nthread_uniform = max_nthread_uniform; + max_nthread_uniform = max_nthread_uniform_default; } - - return bt; } diff --git a/src/gromacs/listed_forces/manage_threading.h b/src/gromacs/listed_forces/manage_threading.h index c2b27d493b..9f369fb35d 100644 --- a/src/gromacs/listed_forces/manage_threading.h +++ b/src/gromacs/listed_forces/manage_threading.h @@ -35,11 +35,10 @@ * To help us fund GROMACS development, we humbly ask that you cite * the research papers on the package. Check out http://www.gromacs.org. */ -/*! \libinternal \file +/*! \internal \file * \brief Declares functions for managing threading of listed forces * * \author Mark Abraham - * \inlibraryapi * \ingroup module_listed_forces */ #ifndef GMX_LISTED_FORCES_MANAGE_THREADING_H @@ -58,20 +57,8 @@ class InteractionDefinitions; * i.e. at start-up without domain decomposition and at DD. */ void setup_bonded_threading(bonded_threading_t* bt, - int numAtoms, - bool useGpuForBondes, + int numAtomsForce, + bool useGpuForBondeds, const InteractionDefinitions& idef); -//! Destructor. -void tear_down_bonded_threading(bonded_threading_t* bt); - -/*! \brief Initialize the bonded threading data structures - * - * Allocates and initializes a bonded threading data structure. - * A pointer to this struct is returned as \p *bb_ptr. - * - * \todo Avoid explicit pointers by using Impl - */ -bonded_threading_t* init_bonded_threading(FILE* fplog, int nenergrp); - #endif diff --git a/src/gromacs/listed_forces/orires.cpp b/src/gromacs/listed_forces/orires.cpp index 8e257559f4..4b35af3192 100644 --- a/src/gromacs/listed_forces/orires.cpp +++ b/src/gromacs/listed_forces/orires.cpp @@ -388,20 +388,17 @@ real calc_orires_dev(const gmx_multisim_t* ms, ArrayRef xWholeMolecules, const rvec x[], const t_pbc* pbc, - t_fcdata* fcd, + t_oriresdata* od, history_t* hist) { - int nref; - real edt, edt_1, invn, pfac, r2, invr, corrfac, wsv2, sw, dev; - OriresMatEq* matEq; - real* mref; - double mtot; - rvec * xref, *xtmp, com, r_unrot, r; - t_oriresdata* od; - gmx_bool bTAV; - const real two_thr = 2.0 / 3.0; - - od = &(fcd->orires); + int nref; + real edt, edt_1, invn, pfac, r2, invr, corrfac, wsv2, sw, dev; + OriresMatEq* matEq; + real* mref; + double mtot; + rvec * xref, *xtmp, com, r_unrot, r; + gmx_bool bTAV; + const real two_thr = 2.0 / 3.0; if (od->nr == 0) { @@ -666,7 +663,7 @@ real orires(int nfa, gmx_bool bTAV; vtot = 0; - od = &(fcd->orires); + od = fcd->orires; if (od->fc != 0) { @@ -753,21 +750,19 @@ real orires(int nfa, /* Approx. 80*nfa/3 flops */ } -void update_orires_history(const t_fcdata* fcd, history_t* hist) +void update_orires_history(const t_oriresdata& od, history_t* hist) { - const t_oriresdata* od = &(fcd->orires); - - if (od->edt != 0) + if (od.edt != 0) { /* Copy the new time averages that have been calculated * in calc_orires_dev. */ - hist->orire_initf = od->exp_min_t_tau; - for (int pair = 0; pair < od->nr; pair++) + hist->orire_initf = od.exp_min_t_tau; + for (int pair = 0; pair < od.nr; pair++) { for (int i = 0; i < 5; i++) { - hist->orire_Dtav[pair * 5 + i] = od->Dtav[pair][i]; + hist->orire_Dtav[pair * 5 + i] = od.Dtav[pair][i]; } } } diff --git a/src/gromacs/listed_forces/orires.h b/src/gromacs/listed_forces/orires.h index 2891ec2dbb..37065e8be5 100644 --- a/src/gromacs/listed_forces/orires.h +++ b/src/gromacs/listed_forces/orires.h @@ -55,7 +55,6 @@ class history_t; struct t_inputrec; struct t_pbc; struct t_commrec; -struct t_fcdata; struct t_oriresdata; class t_state; @@ -94,7 +93,7 @@ real calc_orires_dev(const gmx_multisim_t* ms, gmx::ArrayRef xWholeMolecules, const rvec x[], const t_pbc* pbc, - t_fcdata* fcd, + t_oriresdata* oriresdata, history_t* hist); /*! \brief @@ -123,6 +122,6 @@ real orires(int nfa, int* global_atom_index); //! Copies the new time averages that have been calculated in calc_orires_dev. -void update_orires_history(const t_fcdata* fcd, history_t* hist); +void update_orires_history(const t_oriresdata& oriresdata, history_t* hist); #endif diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp index 567af7d880..589d88fce3 100644 --- a/src/gromacs/mdlib/energyoutput.cpp +++ b/src/gromacs/mdlib/energyoutput.cpp @@ -1176,7 +1176,7 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene, fr.nsum = ebin_->nsum; fr.nre = (bEne) ? ebin_->nener : 0; fr.ener = ebin_->e; - int ndisre = bDR ? fcd->disres.npair : 0; + int ndisre = bDR ? fcd->disres->npair : 0; /* these are for the old-style blocks (1 subblock, only reals), because there can be only one per ID for these */ int nr[enxNR]; @@ -1188,17 +1188,18 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene, nr[i] = 0; } - if (bOR && fcd->orires.nr > 0) + if (bOR && fcd->orires->nr > 0) { - diagonalize_orires_tensors(&(fcd->orires)); - nr[enxOR] = fcd->orires.nr; - block[enxOR] = fcd->orires.otav; + t_oriresdata& orires = *fcd->orires; + diagonalize_orires_tensors(&orires); + nr[enxOR] = orires.nr; + block[enxOR] = orires.otav; id[enxOR] = enxOR; - nr[enxORI] = (fcd->orires.oinsl != fcd->orires.otav) ? fcd->orires.nr : 0; - block[enxORI] = fcd->orires.oinsl; + nr[enxORI] = (orires.oinsl != orires.otav) ? orires.nr : 0; + block[enxORI] = orires.oinsl; id[enxORI] = enxORI; - nr[enxORT] = fcd->orires.nex * 12; - block[enxORT] = fcd->orires.eig; + nr[enxORT] = orires.nex * 12; + block[enxORT] = orires.eig; id[enxORT] = enxORT; } @@ -1237,19 +1238,20 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene, add_blocks_enxframe(&fr, fr.nblock); add_subblocks_enxblock(&(fr.block[db]), 2); - fr.block[db].id = enxDISRE; - fr.block[db].sub[0].nr = ndisre; - fr.block[db].sub[1].nr = ndisre; + const t_disresdata& disres = *fcd->disres; + fr.block[db].id = enxDISRE; + fr.block[db].sub[0].nr = ndisre; + fr.block[db].sub[1].nr = ndisre; #if !GMX_DOUBLE fr.block[db].sub[0].type = xdr_datatype_float; fr.block[db].sub[1].type = xdr_datatype_float; - fr.block[db].sub[0].fval = fcd->disres.rt; - fr.block[db].sub[1].fval = fcd->disres.rm3tav; + fr.block[db].sub[0].fval = disres.rt; + fr.block[db].sub[1].fval = disres.rm3tav; #else fr.block[db].sub[0].type = xdr_datatype_double; fr.block[db].sub[1].type = xdr_datatype_double; - fr.block[db].sub[0].dval = fcd->disres.rt; - fr.block[db].sub[1].dval = fcd->disres.rm3tav; + fr.block[db].sub[0].dval = disres.rt; + fr.block[db].sub[1].dval = disres.rm3tav; #endif } /* here we can put new-style blocks */ @@ -1283,9 +1285,9 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene, free_enxframe(&fr); if (log) { - if (bOR && fcd->orires.nr > 0) + if (bOR && fcd->orires->nr > 0) { - print_orires_log(log, &(fcd->orires)); + print_orires_log(log, fcd->orires); } fprintf(log, " Energies (%s)\n", unit_energy); diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp index 591e64f661..dd5887ced2 100644 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@ -103,7 +103,6 @@ static void reduceEwaldThreadOuput(int nthreads, ewald_corr_thread_t* ewc_t) void do_force_lowlevel(t_forcerec* fr, const t_inputrec* ir, - const InteractionDefinitions& idef, const t_commrec* cr, const gmx_multisim_t* ms, t_nrnb* nrnb, @@ -114,7 +113,6 @@ void do_force_lowlevel(t_forcerec* fr, history_t* hist, gmx::ForceOutputs* forceOutputs, gmx_enerdata_t* enerd, - t_fcdata* fcd, const matrix box, const real* lambda, const rvec* mu_tot, @@ -145,8 +143,9 @@ void do_force_lowlevel(t_forcerec* fr, t_pbc pbc; /* Check whether we need to take into account PBC in listed interactions. */ - const auto needPbcForListedForces = - fr->bMolPBC && stepWork.computeListedForces && haveCpuListedForces(*fr, idef, *fcd); + ListedForces& listedForces = *fr->listedForces; + const auto needPbcForListedForces = + fr->bMolPBC && stepWork.computeListedForces && listedForces.haveCpuListedForces(); if (needPbcForListedForces) { /* Since all atoms are in the rectangular or triclinic unit-cell, @@ -155,9 +154,9 @@ void do_force_lowlevel(t_forcerec* fr, set_pbc_dd(&pbc, fr->pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : nullptr, TRUE, box); } - do_force_listed(wcycle, box, ir->fepvals, cr, ms, idef, x, xWholeMolecules, hist, - forceOutputs, fr, &pbc, enerd, nrnb, lambda, md, fcd, - DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr, stepWork); + listedForces.calculate(wcycle, box, ir->fepvals, cr, ms, x, xWholeMolecules, hist, + forceOutputs, fr, &pbc, enerd, nrnb, lambda, md, + DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr, stepWork); } const bool computePmeOnCpu = (EEL_PME(fr->ic->eeltype) || EVDW_PME(fr->ic->vdwtype)) diff --git a/src/gromacs/mdlib/force.h b/src/gromacs/mdlib/force.h index 387f1e1148..81312671b9 100644 --- a/src/gromacs/mdlib/force.h +++ b/src/gromacs/mdlib/force.h @@ -54,7 +54,6 @@ class history_t; class InteractionDefinitions; struct pull_t; struct t_commrec; -struct t_fcdata; struct t_forcerec; struct t_inputrec; struct t_lambda; @@ -92,7 +91,6 @@ void do_force(FILE* log, tensor vir_force, const t_mdatoms* mdatoms, gmx_enerdata_t* enerd, - t_fcdata* fcd, gmx::ArrayRef lambda, t_forcerec* fr, gmx::MdrunScheduleWorkload* runScheduleWork, @@ -120,7 +118,6 @@ void do_force(FILE* log, */ void do_force_lowlevel(t_forcerec* fr, const t_inputrec* ir, - const InteractionDefinitions& idef, const t_commrec* cr, const gmx_multisim_t* ms, t_nrnb* nrnb, @@ -131,7 +128,6 @@ void do_force_lowlevel(t_forcerec* fr, history_t* hist, gmx::ForceOutputs* forceOutputs, gmx_enerdata_t* enerd, - t_fcdata* fcd, const matrix box, const real* lambda, const rvec* mu_tot, diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index dad37af31c..7e7ab9702d 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -60,7 +60,7 @@ #include "gromacs/gpu_utils/gpu_utils.h" #include "gromacs/hardware/hw_info.h" #include "gromacs/listed_forces/gpubonded.h" -#include "gromacs/listed_forces/manage_threading.h" +#include "gromacs/listed_forces/listed_forces.h" #include "gromacs/listed_forces/pairs.h" #include "gromacs/math/functions.h" #include "gromacs/math/units.h" @@ -552,26 +552,23 @@ static void count_tables(int ftype1, int ftype2, const gmx_mtop_t* mtop, int* nc * * A fatal error occurs if no matching filename is found. */ -static bondedtable_t* make_bonded_tables(FILE* fplog, - int ftype1, - int ftype2, - const gmx_mtop_t* mtop, - gmx::ArrayRef tabbfnm, - const char* tabext) +static std::vector make_bonded_tables(FILE* fplog, + int ftype1, + int ftype2, + const gmx_mtop_t* mtop, + gmx::ArrayRef tabbfnm, + const char* tabext) { - int ncount, *count; - bondedtable_t* tab; + std::vector tab; - tab = nullptr; - - ncount = 0; - count = nullptr; + int ncount = 0; + int* count = nullptr; count_tables(ftype1, ftype2, mtop, &ncount, &count); // Are there any relevant tabulated bond interactions? if (ncount > 0) { - snew(tab, ncount); + tab.resize(ncount); for (int i = 0; i < ncount; i++) { // Do any interactions exist that requires this table? @@ -949,7 +946,6 @@ static void init_interaction_const(FILE* fp, void init_forcerec(FILE* fp, const gmx::MDLogger& mdlog, t_forcerec* fr, - t_fcdata* fcd, const t_inputrec* ir, const gmx_mtop_t* mtop, const t_commrec* cr, @@ -1263,15 +1259,21 @@ void init_forcerec(FILE* fp, make_wall_tables(fp, ir, tabfn, &mtop->groups, fr); } - if (fcd && !tabbfnm.empty()) + /* Initialize the thread working data for bonded interactions */ + fr->listedForces = std::make_unique( + mtop->groups.groups[SimulationAtomGroupType::EnergyOutput].size(), + gmx_omp_nthreads_get(emntBonded), fp); + + if (!tabbfnm.empty()) { + t_fcdata& fcdata = fr->listedForces->fcdata(); // Need to catch std::bad_alloc // TODO Don't need to catch this here, when merging with master branch try { - fcd->bondtab = make_bonded_tables(fp, F_TABBONDS, F_TABBONDSNC, mtop, tabbfnm, "b"); - fcd->angletab = make_bonded_tables(fp, F_TABANGLES, -1, mtop, tabbfnm, "a"); - fcd->dihtab = make_bonded_tables(fp, F_TABDIHS, -1, mtop, tabbfnm, "d"); + fcdata.bondtab = make_bonded_tables(fp, F_TABBONDS, F_TABBONDSNC, mtop, tabbfnm, "b"); + fcdata.angletab = make_bonded_tables(fp, F_TABANGLES, -1, mtop, tabbfnm, "a"); + fcdata.dihtab = make_bonded_tables(fp, F_TABDIHS, -1, mtop, tabbfnm, "d"); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } @@ -1305,10 +1307,6 @@ void init_forcerec(FILE* fp, fr->print_force = print_force; - /* Initialize the thread working data for bonded interactions */ - fr->bondedThreading = init_bonded_threading( - fp, mtop->groups.groups[SimulationAtomGroupType::EnergyOutput].size()); - fr->nthread_ewc = gmx_omp_nthreads_get(emntBonded); snew(fr->ewc_t, fr->nthread_ewc); @@ -1336,5 +1334,4 @@ t_forcerec::~t_forcerec() /* Note: This code will disappear when types are converted to C++ */ sfree(shift_vec); sfree(ewc_t); - tear_down_bonded_threading(bondedThreading); } diff --git a/src/gromacs/mdlib/forcerec.h b/src/gromacs/mdlib/forcerec.h index a3b5f1f915..4fbbee3ab2 100644 --- a/src/gromacs/mdlib/forcerec.h +++ b/src/gromacs/mdlib/forcerec.h @@ -44,7 +44,6 @@ struct gmx_hw_info_t; struct t_commrec; -struct t_fcdata; struct t_forcerec; struct t_filenm; struct t_inputrec; @@ -93,7 +92,6 @@ void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, real table * \param[in] fplog File for printing * \param[in] mdlog File for printing * \param[out] fr The forcerec - * \param[in] fcd Force constant data * \param[in] ir Inputrec structure * \param[in] mtop Molecular topology * \param[in] cr Communication structures @@ -106,7 +104,6 @@ void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, real table void init_forcerec(FILE* fplog, const gmx::MDLogger& mdlog, t_forcerec* fr, - t_fcdata* fcd, const t_inputrec* ir, const gmx_mtop_t* mtop, const t_commrec* cr, diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 5c07b3a507..725cdd16ce 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -64,7 +64,6 @@ #include "gromacs/listed_forces/disre.h" #include "gromacs/listed_forces/gpubonded.h" #include "gromacs/listed_forces/listed_forces.h" -#include "gromacs/listed_forces/manage_threading.h" #include "gromacs/listed_forces/orires.h" #include "gromacs/math/arrayrefwithpadding.h" #include "gromacs/math/functions.h" @@ -830,13 +829,11 @@ static ForceOutputs setupForceOutputs(ForceHelperBuffers* forceH /*! \brief Set up flags that have the lifetime of the domain indicating what type of work is there to compute. */ -static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec& inputrec, - const t_forcerec& fr, - const pull_t* pull_work, - const gmx_edsam* ed, - const InteractionDefinitions& idef, - const t_fcdata& fcd, - const t_mdatoms& mdatoms, +static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec& inputrec, + const t_forcerec& fr, + const pull_t* pull_work, + const gmx_edsam* ed, + const t_mdatoms& mdatoms, const SimulationWorkload& simulationWork, const StepWorkload& stepWork) { @@ -844,10 +841,10 @@ static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec& // Note that haveSpecialForces is constant over the whole run domainWork.haveSpecialForces = haveSpecialForces(inputrec, *fr.forceProviders, pull_work, stepWork.computeForces, ed); - domainWork.haveCpuBondedWork = haveCpuBondeds(fr); + domainWork.haveCpuBondedWork = fr.listedForces->haveCpuBondeds(); domainWork.haveGpuBondedWork = ((fr.gpuBonded != nullptr) && fr.gpuBonded->haveInteractions()); - domainWork.haveRestraintsWork = haveRestraints(idef, fcd); - domainWork.haveCpuListedForceWork = haveCpuListedForces(fr, idef, fcd); + domainWork.haveRestraintsWork = fr.listedForces->haveRestraints(); + domainWork.haveCpuListedForceWork = fr.listedForces->haveCpuListedForces(); // Note that haveFreeEnergyWork is constant over the whole run domainWork.haveFreeEnergyWork = (fr.efep != efepNO && mdatoms.nPerturbed != 0); // We assume we have local force work if there are CPU @@ -1013,7 +1010,6 @@ void do_force(FILE* fplog, tensor vir_force, const t_mdatoms* mdatoms, gmx_enerdata_t* enerd, - t_fcdata* fcd, gmx::ArrayRef lambda, t_forcerec* fr, gmx::MdrunScheduleWorkload* runScheduleWork, @@ -1248,7 +1244,7 @@ void do_force(FILE* fplog, // Need to run after the GPU-offload bonded interaction lists // are set up to be able to determine whether there is bonded work. runScheduleWork->domainWork = setupDomainLifetimeWorkload( - *inputrec, *fr, pull_work, ed, top->idef, *fcd, *mdatoms, simulationWork, stepWork); + *inputrec, *fr, pull_work, ed, *mdatoms, simulationWork, stepWork); wallcycle_start_nocount(wcycle, ewcNS); wallcycle_sub_start(wcycle, ewcsNBS_SEARCH_LOCAL); @@ -1575,9 +1571,9 @@ void do_force(FILE* fplog, stateGpu->waitCoordinatesReadyOnHost(AtomLocality::NonLocal); } /* Compute the bonded and non-bonded energies and optionally forces */ - do_force_lowlevel(fr, inputrec, top->idef, cr, ms, nrnb, wcycle, mdatoms, x, xWholeMolecules, - hist, &forceOut, enerd, fcd, box, lambda.data(), - as_rvec_array(dipoleData.muStateAB), stepWork, ddBalanceRegionHandler); + do_force_lowlevel(fr, inputrec, cr, ms, nrnb, wcycle, mdatoms, x, xWholeMolecules, hist, + &forceOut, enerd, box, lambda.data(), as_rvec_array(dipoleData.muStateAB), + stepWork, ddBalanceRegionHandler); wallcycle_stop(wcycle, ewcFORCE); diff --git a/src/gromacs/mdlib/tests/energyoutput.cpp b/src/gromacs/mdlib/tests/energyoutput.cpp index e22cb28991..422213bb85 100644 --- a/src/gromacs/mdlib/tests/energyoutput.cpp +++ b/src/gromacs/mdlib/tests/energyoutput.cpp @@ -173,8 +173,6 @@ public: std::vector groupNameHandles_; //! Total dipole momentum rvec muTotal_; - //! Distance and orientation restraints data - t_fcdata fcd_; //! Communication record t_commrec cr_; //! Constraints object (for constraints RMSD output in case of LINCS) @@ -378,10 +376,6 @@ public: // It is more relevant to have non-zero value for testing. constraints_ = makeConstraints(mtop_, inputrec_, nullptr, false, nullptr, &cr_, nullptr, nullptr, nullptr, false); - - // No position/orientation restraints - fcd_.disres.npair = 0; - fcd_.orires.nr = 0; } /*! \brief Helper function to generate synthetic data to output diff --git a/src/gromacs/mdlib/tests/leapfrogtestrunners.cpp b/src/gromacs/mdlib/tests/leapfrogtestrunners.cpp index 85f5052a67..939ac5c168 100644 --- a/src/gromacs/mdlib/tests/leapfrogtestrunners.cpp +++ b/src/gromacs/mdlib/tests/leapfrogtestrunners.cpp @@ -86,7 +86,7 @@ void integrateLeapFrogSimple(LeapFrogTestData* testData, int numSteps) { testData->update_->update_coords( testData->inputRecord_, step, &testData->mdAtoms_, &testData->state_, testData->f_, - &testData->forceCalculationData_, &testData->kineticEnergyData_, + testData->forceCalculationData_, &testData->kineticEnergyData_, testData->velocityScalingMatrix_, etrtNONE, nullptr, false); testData->update_->finish_update(testData->inputRecord_, &testData->mdAtoms_, &testData->state_, nullptr, false); diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp index 9b47cf8ca3..0cc8ca8e5b 100644 --- a/src/gromacs/mdlib/update.cpp +++ b/src/gromacs/mdlib/update.cpp @@ -64,6 +64,7 @@ #include "gromacs/mdlib/stat.h" #include "gromacs/mdlib/tgroup.h" #include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/fcdata.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" @@ -124,7 +125,7 @@ public: const t_mdatoms* md, t_state* state, const gmx::ArrayRefWithPadding& f, - const t_fcdata* fcd, + const t_fcdata& fcdata, const gmx_ekindata_t* ekind, const matrix M, int UpdatePart, @@ -198,14 +199,14 @@ void Update::update_coords(const t_inputrec& inpu const t_mdatoms* md, t_state* state, const gmx::ArrayRefWithPadding& f, - const t_fcdata* fcd, + const t_fcdata& fcdata, const gmx_ekindata_t* ekind, const matrix M, int updatePart, const t_commrec* cr, const bool haveConstraints) { - return impl_->update_coords(inputRecord, step, md, state, f, fcd, ekind, M, updatePart, cr, + return impl_->update_coords(inputRecord, step, md, state, f, fcdata, ekind, M, updatePart, cr, haveConstraints); } @@ -1399,7 +1400,7 @@ void Update::Impl::update_coords(const t_inputrec& const t_mdatoms* md, t_state* state, const gmx::ArrayRefWithPadding& f, - const t_fcdata* fcd, + const t_fcdata& fcdata, const gmx_ekindata_t* ekind, const matrix M, int updatePart, @@ -1420,11 +1421,11 @@ void Update::Impl::update_coords(const t_inputrec& /* We need to update the NMR restraint history when time averaging is used */ if (state->flags & (1 << estDISRE_RM3TAV)) { - update_disres_history(fcd, &state->hist); + update_disres_history(*fcdata.disres, &state->hist); } if (state->flags & (1 << estORIRE_DTAV)) { - update_orires_history(fcd, &state->hist); + update_orires_history(*fcdata.orires, &state->hist); } /* ############# START The update of velocities and positions ######### */ diff --git a/src/gromacs/mdlib/update.h b/src/gromacs/mdlib/update.h index b7853f3902..3f122b59f7 100644 --- a/src/gromacs/mdlib/update.h +++ b/src/gromacs/mdlib/update.h @@ -107,7 +107,7 @@ public: * \param[in] md MD atoms data. * \param[in] state System state object. * \param[in] f Buffer with atomic forces for home particles. - * \param[in] fcd Force calculation data to update distance and orientation restraints. + * \param[in] fcdata Force calculation data to update distance and orientation restraints. * \param[in] ekind Kinetic energy data (for temperature coupling, energy groups, etc.). * \param[in] M Parrinello-Rahman velocity scaling matrix. * \param[in] updatePart What should be updated, coordinates or velocities. This enum only used in VV integrator. @@ -119,7 +119,7 @@ public: const t_mdatoms* md, t_state* state, const gmx::ArrayRefWithPadding& f, - const t_fcdata* fcd, + const t_fcdata& fcdata, const gmx_ekindata_t* ekind, const matrix M, int updatePart, diff --git a/src/gromacs/mdrun/isimulator.h b/src/gromacs/mdrun/isimulator.h index c51aae29a4..c6a74e848a 100644 --- a/src/gromacs/mdrun/isimulator.h +++ b/src/gromacs/mdrun/isimulator.h @@ -57,7 +57,6 @@ struct ObservablesHistory; struct pull_t; struct ReplicaExchangeParameters; struct t_commrec; -struct t_fcdata; struct t_forcerec; struct t_filenm; struct t_inputrec; @@ -119,7 +118,6 @@ public: pull_t* pull_work, t_swap* swap, gmx_mtop_t* top_global, - t_fcdata* fcd, t_state* state_global, ObservablesHistory* observablesHistory, MDAtoms* mdAtoms, @@ -154,7 +152,6 @@ public: pull_work(pull_work), swap(swap), top_global(top_global), - fcd(fcd), state_global(state_global), observablesHistory(observablesHistory), mdAtoms(mdAtoms), @@ -213,8 +210,6 @@ protected: t_swap* swap; //! Full system topology. const gmx_mtop_t* top_global; - //! Helper struct for force calculations. - t_fcdata* fcd; //! Full simulation state (only non-nullptr on master rank). t_state* state_global; //! History of simulation observables. diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index b2eb8fa378..4ace36f893 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -70,7 +70,7 @@ #include "gromacs/gpu_utils/device_stream_manager.h" #include "gromacs/gpu_utils/gpu_utils.h" #include "gromacs/imd/imd.h" -#include "gromacs/listed_forces/manage_threading.h" +#include "gromacs/listed_forces/listed_forces.h" #include "gromacs/math/functions.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" @@ -258,13 +258,15 @@ void gmx::LegacySimulator::do_md() const bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd); const bool useReplicaExchange = (replExParams.exchangeInterval > 0); + const t_fcdata& fcdata = fr->listedForces->fcdata(); + bool simulationsShareState = false; int nstSignalComm = nstglobalcomm; { // TODO This implementation of ensemble orientation restraints is nasty because // a user can't just do multi-sim with single-sim orientation restraints. bool usingEnsembleRestraints = - (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0)); + (fcdata.disres->nsystems > 1) || ((ms != nullptr) && (fcdata.orires->nr != 0)); bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr)); // Replica exchange, ensemble restraints and AWH need all @@ -380,7 +382,7 @@ void gmx::LegacySimulator::do_md() "Essential dynamics is not supported with the GPU update.\n"); GMX_RELEASE_ASSERT(!ir->bPull || !pull_have_constraint(ir->pull), "Constraints pulling is not supported with the GPU update.\n"); - GMX_RELEASE_ASSERT(fcd->orires.nr == 0, + GMX_RELEASE_ASSERT(fcdata.orires->nr == 0, "Orientation restraints are not supported with the GPU update.\n"); GMX_RELEASE_ASSERT( ir->efep == efepNO @@ -944,7 +946,7 @@ void gmx::LegacySimulator::do_md() { /* Now is the time to relax the shells */ relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir, - imdSession, pull_work, bNS, force_flags, &top, constr, enerd, fcd, + imdSession, pull_work, bNS, force_flags, &top, constr, enerd, state->natoms, state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(), state->box, state->lambda, &state->hist, f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, shellfc, @@ -972,7 +974,7 @@ void gmx::LegacySimulator::do_md() */ do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, fr, + f.arrayRefWithPadding(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr, (bNS ? GMX_FORCE_NS : 0) | force_flags, ddBalanceRegionHandler); } @@ -1006,7 +1008,7 @@ void gmx::LegacySimulator::do_md() trotter_seq, ettTSEQ1); } - upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, + upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcdata, ekind, M, etrtVELOCITY1, cr, constr != nullptr); wallcycle_stop(wcycle, ewcUPDATE); @@ -1256,7 +1258,7 @@ void gmx::LegacySimulator::do_md() if (EI_VV(ir->eI)) { /* velocity half-step update */ - upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, + upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcdata, ekind, M, etrtVELOCITY2, cr, constr != nullptr); } @@ -1322,7 +1324,7 @@ void gmx::LegacySimulator::do_md() } else { - upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, + upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcdata, ekind, M, etrtPOSITION, cr, constr != nullptr); wallcycle_stop(wcycle, ewcUPDATE); @@ -1354,7 +1356,7 @@ void gmx::LegacySimulator::do_md() /* now we know the scaling, we can compute the positions again */ std::copy(cbuf.begin(), cbuf.end(), state->x.begin()); - upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, + upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcdata, ekind, M, etrtPOSITION, cr, constr != nullptr); wallcycle_stop(wcycle, ewcUPDATE); @@ -1555,7 +1557,8 @@ void gmx::LegacySimulator::do_md() if (do_log || do_ene || do_dr || do_or) { energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, - do_log ? fplog : nullptr, step, t, fcd, awh.get()); + do_log ? fplog : nullptr, step, t, + &fr->listedForces->fcdata(), awh.get()); } if (ir->bPull) diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index 0b7e667ff4..ef5a04f3d5 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -66,7 +66,7 @@ #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nrnb.h" #include "gromacs/gpu_utils/gpu_utils.h" -#include "gromacs/listed_forces/manage_threading.h" +#include "gromacs/listed_forces/listed_forces.h" #include "gromacs/math/functions.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" @@ -104,7 +104,6 @@ #include "gromacs/mdtypes/df_history.h" #include "gromacs/mdtypes/enerdata.h" #include "gromacs/mdtypes/energyhistory.h" -#include "gromacs/mdtypes/fcdata.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" @@ -414,7 +413,7 @@ void gmx::LegacySimulator::do_mimic() { /* Now is the time to relax the shells */ relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir, - imdSession, pull_work, bNS, force_flags, &top, constr, enerd, fcd, + imdSession, pull_work, bNS, force_flags, &top, constr, enerd, state->natoms, state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(), state->box, state->lambda, &state->hist, f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, shellfc, @@ -431,9 +430,8 @@ void gmx::LegacySimulator::do_mimic() gmx_edsam* ed = nullptr; do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, fr, - runScheduleWork, vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, - ddBalanceRegionHandler); + f.arrayRefWithPadding(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, + vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler); } /* Now we have the energies and forces corresponding to the @@ -528,7 +526,8 @@ void gmx::LegacySimulator::do_mimic() EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts)); energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, - do_log ? fplog : nullptr, step, t, fcd, awh); + do_log ? fplog : nullptr, step, t, + &fr->listedForces->fcdata(), awh); if (do_per_step(step, ir->nstlog)) { diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index ddd0870ba6..64882158e2 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -68,7 +68,7 @@ #include "gromacs/gmxlib/nrnb.h" #include "gromacs/imd/imd.h" #include "gromacs/linearalgebra/sparsematrix.h" -#include "gromacs/listed_forces/manage_threading.h" +#include "gromacs/listed_forces/listed_forces.h" #include "gromacs/math/functions.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/constr.h" @@ -786,8 +786,6 @@ public: VirtualSitesHandler* vsite; //! Handles constraints. gmx::Constraints* constr; - //! Handles strange things. - t_fcdata* fcd; //! Per-atom data for this domain. gmx::MDAtoms* mdAtoms; //! Handles how to calculate the forces. @@ -841,8 +839,8 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres, */ do_force(fplog, cr, ms, inputrec, nullptr, nullptr, imdSession, pull_work, count, nrnb, wcycle, top, ems->s.box, ems->s.x.arrayRefWithPadding(), &ems->s.hist, - ems->f.arrayRefWithPadding(), force_vir, mdAtoms->mdatoms(), enerd, fcd, ems->s.lambda, - fr, runScheduleWork, vsite, mu_tot, t, nullptr, + ems->f.arrayRefWithPadding(), force_vir, mdAtoms->mdatoms(), enerd, ems->s.lambda, fr, + runScheduleWork, vsite, mu_tot, t, nullptr, GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES | GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY | (bNS ? GMX_FORCE_NS : 0), DDBalanceRegionHandler(cr)); @@ -1106,11 +1104,9 @@ void LegacySimulator::do_cg() sp_header(fplog, CG, inputrec->em_tol, number_steps); } - EnergyEvaluator energyEvaluator{ - fplog, mdlog, cr, ms, top_global, &top, inputrec, imdSession, pull_work, - nrnb, wcycle, gstat, vsite, constr, fcd, mdAtoms, fr, runScheduleWork, - enerd - }; + EnergyEvaluator energyEvaluator{ fplog, mdlog, cr, ms, top_global, &top, + inputrec, imdSession, pull_work, nrnb, wcycle, gstat, + vsite, constr, mdAtoms, fr, runScheduleWork, enerd }; /* Call the force routine and some auxiliary (neighboursearching etc.) */ /* do_force always puts the charge groups in the box and shifts again * We do not unshift, so molecules are always whole in congrad.c @@ -1127,7 +1123,7 @@ void LegacySimulator::do_cg() EnergyOutput::printHeader(fplog, step, step); energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, - step, fcd, nullptr); + step, &fr->listedForces->fcdata(), nullptr); } /* Estimate/guess the initial stepsize */ @@ -1542,7 +1538,8 @@ void LegacySimulator::do_cg() EnergyOutput::printHeader(fplog, step, step); } energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE, - do_log ? fplog : nullptr, step, step, fcd, nullptr); + do_log ? fplog : nullptr, step, step, + &fr->listedForces->fcdata(), nullptr); } /* Send energies and positions to the IMD client if bIMD is TRUE. */ @@ -1585,7 +1582,8 @@ void LegacySimulator::do_cg() { /* Write final energy file entries */ energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE, - !do_log ? fplog : nullptr, step, step, fcd, nullptr); + !do_log ? fplog : nullptr, step, step, + &fr->listedForces->fcdata(), nullptr); } } @@ -1758,11 +1756,9 @@ void LegacySimulator::do_lbfgs() * We do not unshift, so molecules are always whole */ neval++; - EnergyEvaluator energyEvaluator{ - fplog, mdlog, cr, ms, top_global, &top, inputrec, imdSession, pull_work, - nrnb, wcycle, gstat, vsite, constr, fcd, mdAtoms, fr, runScheduleWork, - enerd - }; + EnergyEvaluator energyEvaluator{ fplog, mdlog, cr, ms, top_global, &top, + inputrec, imdSession, pull_work, nrnb, wcycle, gstat, + vsite, constr, mdAtoms, fr, runScheduleWork, enerd }; energyEvaluator.run(&ems, mu_tot, vir, pres, -1, TRUE); if (MASTER(cr)) @@ -1775,7 +1771,7 @@ void LegacySimulator::do_lbfgs() EnergyOutput::printHeader(fplog, step, step); energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, - step, fcd, nullptr); + step, &fr->listedForces->fcdata(), nullptr); } /* Set the initial step. @@ -2267,7 +2263,8 @@ void LegacySimulator::do_lbfgs() EnergyOutput::printHeader(fplog, step, step); } energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE, - do_log ? fplog : nullptr, step, step, fcd, nullptr); + do_log ? fplog : nullptr, step, step, + &fr->listedForces->fcdata(), nullptr); } /* Send x and E to IMD client, if bIMD is TRUE. */ @@ -2309,7 +2306,8 @@ void LegacySimulator::do_lbfgs() if (!do_ene || !do_log) /* Write final energy file entries */ { energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE, - !do_log ? fplog : nullptr, step, step, fcd, nullptr); + !do_log ? fplog : nullptr, step, step, + &fr->listedForces->fcdata(), nullptr); } /* Print some stuff... */ @@ -2404,11 +2402,9 @@ void LegacySimulator::do_steep() { sp_header(fplog, SD, inputrec->em_tol, nsteps); } - EnergyEvaluator energyEvaluator{ - fplog, mdlog, cr, ms, top_global, &top, inputrec, imdSession, pull_work, - nrnb, wcycle, gstat, vsite, constr, fcd, mdAtoms, fr, runScheduleWork, - enerd - }; + EnergyEvaluator energyEvaluator{ fplog, mdlog, cr, ms, top_global, &top, + inputrec, imdSession, pull_work, nrnb, wcycle, gstat, + vsite, constr, mdAtoms, fr, runScheduleWork, enerd }; /**** HERE STARTS THE LOOP **** * count is the counter for the number of steps @@ -2473,8 +2469,8 @@ void LegacySimulator::do_steep() const bool do_dr = do_per_step(steps_accepted, inputrec->nstdisreout); const bool do_or = do_per_step(steps_accepted, inputrec->nstorireout); - energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, do_dr, do_or, - fplog, count, count, fcd, nullptr); + energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, do_dr, do_or, fplog, + count, count, &fr->listedForces->fcdata(), nullptr); fflush(fplog); } } @@ -2697,11 +2693,9 @@ void LegacySimulator::do_nm() /* Make evaluate_energy do a single node force calculation */ cr->nnodes = 1; - EnergyEvaluator energyEvaluator{ - fplog, mdlog, cr, ms, top_global, &top, inputrec, imdSession, pull_work, - nrnb, wcycle, gstat, vsite, constr, fcd, mdAtoms, fr, runScheduleWork, - enerd - }; + EnergyEvaluator energyEvaluator{ fplog, mdlog, cr, ms, top_global, &top, + inputrec, imdSession, pull_work, nrnb, wcycle, gstat, + vsite, constr, mdAtoms, fr, runScheduleWork, enerd }; energyEvaluator.run(&state_work, mu_tot, vir, pres, -1, TRUE); cr->nnodes = nnodes; @@ -2761,7 +2755,7 @@ void LegacySimulator::do_nm() /* Now is the time to relax the shells */ relax_shell_flexcon( fplog, cr, ms, mdrunOptions.verbose, nullptr, step, inputrec, imdSession, - pull_work, bNS, force_flags, &top, constr, enerd, fcd, state_work.s.natoms, + pull_work, bNS, force_flags, &top, constr, enerd, state_work.s.natoms, state_work.s.x.arrayRefWithPadding(), state_work.s.v.arrayRefWithPadding(), state_work.s.box, state_work.s.lambda, &state_work.s.hist, state_work.f.arrayRefWithPadding(), vir, mdatoms, nrnb, wcycle, shellfc, diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index 5a66c06a4d..164c508238 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -65,7 +65,7 @@ #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nrnb.h" #include "gromacs/gpu_utils/gpu_utils.h" -#include "gromacs/listed_forces/manage_threading.h" +#include "gromacs/listed_forces/listed_forces.h" #include "gromacs/math/functions.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" @@ -102,7 +102,6 @@ #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/df_history.h" #include "gromacs/mdtypes/energyhistory.h" -#include "gromacs/mdtypes/fcdata.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" @@ -522,7 +521,7 @@ void gmx::LegacySimulator::do_rerun() { /* Now is the time to relax the shells */ relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir, - imdSession, pull_work, bNS, force_flags, &top, constr, enerd, fcd, + imdSession, pull_work, bNS, force_flags, &top, constr, enerd, state->natoms, state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(), state->box, state->lambda, &state->hist, f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, shellfc, @@ -539,9 +538,8 @@ void gmx::LegacySimulator::do_rerun() gmx_edsam* ed = nullptr; do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, fr, - runScheduleWork, vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, - ddBalanceRegionHandler); + f.arrayRefWithPadding(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, + vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler); } /* Now we have the energies and forces corresponding to the @@ -614,7 +612,8 @@ void gmx::LegacySimulator::do_rerun() EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts)); energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, - do_log ? fplog : nullptr, step, t, fcd, awh); + do_log ? fplog : nullptr, step, t, + &fr->listedForces->fcdata(), awh); if (do_per_step(step, ir->nstlog)) { diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index 928a249dc3..646bfba264 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -82,6 +82,7 @@ #include "gromacs/imd/imd.h" #include "gromacs/listed_forces/disre.h" #include "gromacs/listed_forces/gpubonded.h" +#include "gromacs/listed_forces/listed_forces.h" #include "gromacs/listed_forces/orires.h" #include "gromacs/math/functions.h" #include "gromacs/math/utilities.h" @@ -704,7 +705,6 @@ int Mdrunner::mdrunner() { matrix box; t_forcerec* fr = nullptr; - t_fcdata* fcd = nullptr; real ewaldcoeff_q = 0; real ewaldcoeff_lj = 0; int nChargePerturbed = -1, nTypePerturbed = 0; @@ -1034,14 +1034,17 @@ int Mdrunner::mdrunner() * So the PME-only nodes (if present) will also initialize * the distance restraints. */ - snew(fcd, 1); /* This needs to be called before read_checkpoint to extend the state */ + t_disresdata* disresdata; + snew(disresdata, 1); init_disres(fplog, &mtop, inputrec, DisResRunMode::MDRun, MASTER(cr) ? DDRole::Master : DDRole::Agent, - PAR(cr) ? NumRanks::Multiple : NumRanks::Single, cr->mpi_comm_mysim, ms, fcd, + PAR(cr) ? NumRanks::Multiple : NumRanks::Single, cr->mpi_comm_mysim, ms, disresdata, globalState.get(), replExParams.exchangeInterval > 0); - init_orires(fplog, &mtop, inputrec, cr, ms, globalState.get(), &(fcd->orires)); + t_oriresdata* oriresdata; + snew(oriresdata, 1); + init_orires(fplog, &mtop, inputrec, cr, ms, globalState.get(), oriresdata); auto deform = prepareBoxDeformation(globalState->box, MASTER(cr) ? DDRole::Master : DDRole::Agent, PAR(cr) ? NumRanks::Multiple : NumRanks::Single, @@ -1364,10 +1367,13 @@ int Mdrunner::mdrunner() /* Initiate forcerecord */ fr = new t_forcerec; fr->forceProviders = mdModules_->initForceProviders(); - init_forcerec(fplog, mdlog, fr, fcd, inputrec, &mtop, cr, box, + init_forcerec(fplog, mdlog, fr, inputrec, &mtop, cr, box, opt2fn("-table", filenames.size(), filenames.data()), opt2fn("-tablep", filenames.size(), filenames.data()), opt2fns("-tableb", filenames.size(), filenames.data()), pforce); + // Dirty hack, for fixing disres and orires should be made mdmodules + fr->listedForces->fcdata().disres = disresdata; + fr->listedForces->fcdata().orires = oriresdata; // Save a handle to device stream manager to use elsewhere in the code // TODO: Forcerec is not a correct place to store it. @@ -1640,7 +1646,7 @@ int Mdrunner::mdrunner() filenames.data(), oenv, mdrunOptions, startingBehavior, vsite.get(), constr.get(), enforcedRotation ? enforcedRotation->getLegacyEnfrot() : nullptr, deform.get(), mdModules_->outputProvider(), mdModules_->notifier(), inputrec, imdSession.get(), - pull_work, swap, &mtop, fcd, globalState.get(), &observablesHistory, mdAtoms.get(), + pull_work, swap, &mtop, globalState.get(), &observablesHistory, mdAtoms.get(), &nrnb, wcycle, fr, &enerd, &ekind, &runScheduleWork, replExParams, membed, walltime_accounting, std::move(stopHandlerBuilder_), doRerun); simulator->run(); @@ -1696,6 +1702,9 @@ int Mdrunner::mdrunner() /* Free pinned buffers in *fr */ delete fr; fr = nullptr; + // TODO convert to C++ so we can get rid of these frees + sfree(disresdata); + sfree(oriresdata); if (hwinfo->gpu_info.n_dev > 0) { @@ -1725,7 +1734,6 @@ int Mdrunner::mdrunner() } free_gpu(deviceInfo); - sfree(fcd); if (doMembed) { diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index 9d79a88962..01ddc14ceb 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -910,7 +910,6 @@ void relax_shell_flexcon(FILE* fplog, const gmx_localtop_t* top, gmx::Constraints* constr, gmx_enerdata_t* enerd, - t_fcdata* fcd, int natoms, ArrayRefWithPadding xPadded, ArrayRefWithPadding vPadded, @@ -1022,7 +1021,7 @@ void relax_shell_flexcon(FILE* fplog, int shellfc_flags = force_flags | (bVerbose ? GMX_FORCE_ENERGY : 0); do_force(fplog, cr, ms, inputrec, nullptr, enforcedRotation, imdSession, pull_work, mdstep, nrnb, wcycle, top, box, xPadded, hist, forceWithPadding[Min], force_vir, md, enerd, - fcd, lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, + lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, (bDoNS ? GMX_FORCE_NS : 0) | shellfc_flags, ddBalanceRegionHandler); sf_dir = 0; @@ -1108,7 +1107,7 @@ void relax_shell_flexcon(FILE* fplog, /* Try the new positions */ do_force(fplog, cr, ms, inputrec, nullptr, enforcedRotation, imdSession, pull_work, 1, nrnb, wcycle, top, box, posWithPadding[Try], hist, forceWithPadding[Try], force_vir, md, - enerd, fcd, lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, shellfc_flags, + enerd, lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, shellfc_flags, ddBalanceRegionHandler); sum_epot(&(enerd->grpp), enerd->term); if (gmx_debug_at) diff --git a/src/gromacs/mdrun/shellfc.h b/src/gromacs/mdrun/shellfc.h index 7b5be6592d..1a6248334b 100644 --- a/src/gromacs/mdrun/shellfc.h +++ b/src/gromacs/mdrun/shellfc.h @@ -54,7 +54,6 @@ struct gmx_mtop_t; class history_t; struct pull_t; struct t_forcerec; -struct t_fcdata; struct t_inputrec; struct t_mdatoms; struct t_nrnb; @@ -96,7 +95,6 @@ void relax_shell_flexcon(FILE* log, const gmx_localtop_t* top, gmx::Constraints* constr, gmx_enerdata_t* enerd, - t_fcdata* fcd, int natoms, gmx::ArrayRefWithPadding x, gmx::ArrayRefWithPadding v, diff --git a/src/gromacs/mdrun/simulatorbuilder.h b/src/gromacs/mdrun/simulatorbuilder.h index 9b938b7224..18251ebf7d 100644 --- a/src/gromacs/mdrun/simulatorbuilder.h +++ b/src/gromacs/mdrun/simulatorbuilder.h @@ -61,7 +61,6 @@ struct ObservablesHistory; struct pull_t; struct ReplicaExchangeParameters; struct t_commrec; -struct t_fcdata; struct t_forcerec; struct t_filenm; struct t_inputrec; diff --git a/src/gromacs/mdrun/tpi.cpp b/src/gromacs/mdrun/tpi.cpp index 280d1f22cb..6239d06ba0 100644 --- a/src/gromacs/mdrun/tpi.cpp +++ b/src/gromacs/mdrun/tpi.cpp @@ -754,7 +754,7 @@ void LegacySimulator::do_tpi() std::feholdexcept(&floatingPointEnvironment); do_force(fplog, cr, ms, inputrec, nullptr, nullptr, imdSession, pull_work, step, nrnb, wcycle, &top, state_global->box, state_global->x.arrayRefWithPadding(), - &state_global->hist, f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, + &state_global->hist, f.arrayRefWithPadding(), force_vir, mdatoms, enerd, state_global->lambda, fr, runScheduleWork, nullptr, mu_tot, t, nullptr, GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY | (bStateChanged ? GMX_FORCE_STATECHANGED : 0), DDBalanceRegionHandler(nullptr)); diff --git a/src/gromacs/mdtypes/fcdata.h b/src/gromacs/mdtypes/fcdata.h index 927db7ac89..09dce8f273 100644 --- a/src/gromacs/mdtypes/fcdata.h +++ b/src/gromacs/mdtypes/fcdata.h @@ -38,6 +38,8 @@ #ifndef GMX_MDTYPES_FCDATA_H #define GMX_MDTYPES_FCDATA_H +#include + #include "gromacs/math/vectypes.h" #include "gromacs/topology/idef.h" #include "gromacs/utility/basedefinitions.h" @@ -112,12 +114,13 @@ typedef struct t_oriresdata double** v; } t_oriresdata; -typedef struct bondedtable_t +/* Cubic spline table for tabulated bonded interactions */ +struct bondedtable_t { - int n; /* n+1 is the number of points */ - real scale; /* distance between two points */ - real* data; /* the actual table data, per point there are 4 numbers */ -} bondedtable_t; + int n; /* n+1 is the number of points */ + real scale; /* distance between two points */ + std::vector data; /* the actual table data, per point there are 4 numbers */ +}; /* * Data struct used in the force calculation routines @@ -126,14 +129,15 @@ typedef struct bondedtable_t * (for instance for time averaging in distance retraints) * or for storing output, since force routines only return the potential. */ -typedef struct t_fcdata +struct t_fcdata { - bondedtable_t* bondtab; - bondedtable_t* angletab; - bondedtable_t* dihtab; + std::vector bondtab; + std::vector angletab; + std::vector dihtab; - t_disresdata disres; - t_oriresdata orires; -} t_fcdata; + // TODO: Convert to C++ and unique_ptr (currently this data is not freed) + t_disresdata* disres = nullptr; + t_oriresdata* orires = nullptr; +}; #endif diff --git a/src/gromacs/mdtypes/forcerec.h b/src/gromacs/mdtypes/forcerec.h index 044f8d332c..ea81170069 100644 --- a/src/gromacs/mdtypes/forcerec.h +++ b/src/gromacs/mdtypes/forcerec.h @@ -55,6 +55,7 @@ struct nonbonded_verlet_t; struct bonded_threading_t; class DeviceContext; class DispersionCorrection; +class ListedForces; struct t_forcetable; struct t_QMMMrec; @@ -297,8 +298,8 @@ struct t_forcerec real userreal3 = 0; real userreal4 = 0; - /* Pointer to struct for managing threading of bonded force calculation */ - struct bonded_threading_t* bondedThreading = nullptr; + /* The listed forces calculation data */ + std::unique_ptr listedForces; /* TODO: Replace the pointer by an object once we got rid of C */ gmx::GpuBonded* gpuBonded = nullptr; diff --git a/src/gromacs/modularsimulator/forceelement.cpp b/src/gromacs/modularsimulator/forceelement.cpp index d2ca77e8b2..86122ccf12 100644 --- a/src/gromacs/modularsimulator/forceelement.cpp +++ b/src/gromacs/modularsimulator/forceelement.cpp @@ -76,7 +76,6 @@ ForceElement::ForceElement(StatePropagatorData* statePropagatorData, const MDAtoms* mdAtoms, t_nrnb* nrnb, t_forcerec* fr, - t_fcdata* fcd, gmx_wallcycle* wcycle, MdrunScheduleWorkload* runScheduleWork, VirtualSitesHandler* vsite, @@ -114,7 +113,6 @@ ForceElement::ForceElement(StatePropagatorData* statePropagatorData, vsite_(vsite), imdSession_(imdSession), pull_work_(pull_work), - fcd_(fcd), runScheduleWork_(runScheduleWork), constr_(constr), enforcedRotation_(enforcedRotation) @@ -192,7 +190,7 @@ void ForceElement::run(Step step, Time time, unsigned int flags) relax_shell_flexcon( fplog_, cr_, ms, isVerbose_, enforcedRotation_, step, inputrec_, imdSession_, pull_work_, step == nextNSStep_, static_cast(flags), localTopology_, constr_, - energyElement_->enerdata(), fcd_, statePropagatorData_->localNumAtoms(), x, v, box, + energyElement_->enerdata(), statePropagatorData_->localNumAtoms(), x, v, box, lambda, hist, forces, force_vir, mdAtoms_->mdatoms(), nrnb_, wcycle_, shellfc_, fr_, runScheduleWork_, time, energyElement_->muTot(), vsite_, ddBalanceRegionHandler_); nShellRelaxationSteps_++; @@ -205,7 +203,7 @@ void ForceElement::run(Step step, Time time, unsigned int flags) do_force(fplog_, cr_, ms, inputrec_, awh, enforcedRotation_, imdSession_, pull_work_, step, nrnb_, wcycle_, localTopology_, box, x, hist, forces, force_vir, mdAtoms_->mdatoms(), - energyElement_->enerdata(), fcd_, lambda, fr_, runScheduleWork_, vsite_, + energyElement_->enerdata(), lambda, fr_, runScheduleWork_, vsite_, energyElement_->muTot(), time, ed, static_cast(flags), ddBalanceRegionHandler_); } energyElement_->addToForceVirial(force_vir, step); diff --git a/src/gromacs/modularsimulator/forceelement.h b/src/gromacs/modularsimulator/forceelement.h index 9abbe21ec7..26e7661d03 100644 --- a/src/gromacs/modularsimulator/forceelement.h +++ b/src/gromacs/modularsimulator/forceelement.h @@ -58,7 +58,6 @@ struct gmx_enfrot; struct gmx_shellfc_t; struct gmx_wallcycle; struct pull_t; -struct t_fcdata; struct t_nrnb; namespace gmx @@ -98,7 +97,6 @@ public: const MDAtoms* mdAtoms, t_nrnb* nrnb, t_forcerec* fr, - t_fcdata* fcd, gmx_wallcycle* wcycle, MdrunScheduleWorkload* runScheduleWork, VirtualSitesHandler* vsite, @@ -194,8 +192,6 @@ private: ImdSession* imdSession_; //! The pull work object. pull_t* pull_work_; - //! Helper struct for force calculations. - t_fcdata* fcd_; //! Schedule of work for each MD step for this task. MdrunScheduleWorkload* runScheduleWork_; //! Handles constraints. diff --git a/src/gromacs/modularsimulator/modularsimulator.cpp b/src/gromacs/modularsimulator/modularsimulator.cpp index a64b56e03a..896d7945cf 100644 --- a/src/gromacs/modularsimulator/modularsimulator.cpp +++ b/src/gromacs/modularsimulator/modularsimulator.cpp @@ -50,6 +50,7 @@ #include "gromacs/ewald/pme_pp.h" #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nrnb.h" +#include "gromacs/listed_forces/listed_forces.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/checkpointhandler.h" #include "gromacs/mdlib/constr.h" @@ -381,8 +382,8 @@ void ModularSimulator::constructElementsAndSignallers() auto energyElement = std::make_unique( statePropagatorDataPtr, freeEnergyPerturbationElementPtr, top_global, inputrec, mdAtoms, - enerd, ekind, constr, fplog, fcd, mdModulesNotifier, MASTER(cr), observablesHistory, - startingBehavior); + enerd, ekind, constr, fplog, &fr->listedForces->fcdata(), mdModulesNotifier, MASTER(cr), + observablesHistory, startingBehavior); auto energyElementPtr = compat::make_not_null(energyElement.get()); /* @@ -552,7 +553,7 @@ ModularSimulator::buildForces(SignallerBuilder* neighbo auto forceElement = std::make_unique( statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElement, isVerbose, - isDynamicBox, fplog, cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle, runScheduleWork, vsite, + isDynamicBox, fplog, cr, inputrec, mdAtoms, nrnb, fr, wcycle, runScheduleWork, vsite, imdSession, pull_work, constr, &topologyHolder_->globalTopology(), enforcedRotation); topologyHolder_->registerClient(forceElement.get()); neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get())); @@ -881,7 +882,7 @@ bool ModularSimulator::isInputCompatible(bool exitOn int numEnsembleRestraintSystems; if (fcd) { - numEnsembleRestraintSystems = fcd->disres.nsystems; + numEnsembleRestraintSystems = fcd->disres->nsystems; } else { @@ -944,8 +945,8 @@ bool ModularSimulator::isInputCompatible(bool exitOn void ModularSimulator::checkInputForDisabledFunctionality() { - isInputCompatible(true, inputrec, doRerun, *top_global, ms, replExParams, fcd, - opt2bSet("-ei", nfile, fnm), membed != nullptr); + isInputCompatible(true, inputrec, doRerun, *top_global, ms, replExParams, + &fr->listedForces->fcdata(), opt2bSet("-ei", nfile, fnm), membed != nullptr); if (observablesHistory->edsamHistory) { gmx_fatal(FARGS, diff --git a/src/gromacs/modularsimulator/modularsimulator.h b/src/gromacs/modularsimulator/modularsimulator.h index f5f6bb4e99..0724501364 100644 --- a/src/gromacs/modularsimulator/modularsimulator.h +++ b/src/gromacs/modularsimulator/modularsimulator.h @@ -56,6 +56,8 @@ #include "pmeloadbalancehelper.h" #include "topologyholder.h" +struct t_fcdata; + namespace gmx { class DomDecHelper; diff --git a/src/gromacs/tables/forcetable.cpp b/src/gromacs/tables/forcetable.cpp index af3382855a..97b824ce78 100644 --- a/src/gromacs/tables/forcetable.cpp +++ b/src/gromacs/tables/forcetable.cpp @@ -1396,8 +1396,8 @@ bondedtable_t make_bonded_table(FILE* fplog, const char* fn, int angle) } tab.n = td.nx; tab.scale = td.tabscale; - snew(tab.data, tab.n * stride); - copy2table(tab.n, 0, stride, td.x, td.v, td.f, 1.0, tab.data); + tab.data.resize(tab.n * stride); + copy2table(tab.n, 0, stride, td.x, td.v, td.f, 1.0, tab.data.data()); done_tabledata(&td); return tab; diff --git a/src/gromacs/topology/topsort.cpp b/src/gromacs/topology/topsort.cpp index 3ce7318e75..5971a45094 100644 --- a/src/gromacs/topology/topsort.cpp +++ b/src/gromacs/topology/topsort.cpp @@ -187,6 +187,8 @@ void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* qB = qA; } + bool havePerturbedInteractions = false; + iabuf_nalloc = 0; iabuf = nullptr; @@ -215,6 +217,8 @@ void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* { iabuf[ib++] = iatoms[i++]; } + + havePerturbedInteractions = true; } else { @@ -245,5 +249,5 @@ void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* sfree(iabuf); - idef->ilsort = ilsortFE_SORTED; + idef->ilsort = (havePerturbedInteractions ? ilsortFE_SORTED : ilsortNO_FE); }