#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"
#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"
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))
{
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;
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());)
} 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)
{
{
if (index[j] == forceparams[type].disres.label)
{
- vvindex[j] = gmx::invsixthroot(fcd->disres.Rt_6[label]);
+ vvindex[j] = gmx::invsixthroot(disresdata->Rt_6[label]);
}
}
}
};
FILE * out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
- t_fcdata fcd;
int i, j, kkk;
t_trxstatus* status;
real t;
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)
}
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
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;
}
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)
{
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);
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<const t_functype> functype = top->idef.functype;
- gmx::ArrayRef<const t_iparams> iparams = top->idef.iparams;
+ gmx::ArrayRef<const t_functype> functype = idef.functype;
+ gmx::ArrayRef<const t_iparams> 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");
/* Fill the index array */
label1 = -1;
- const InteractionList& disres = top->idef.il[F_DISRES];
+ const InteractionList& disres = idef.il[F_DISRES];
gmx::ArrayRef<const int> iatom = disres.iatoms;
for (i = j = k = 0; (i < disres.size());)
{
top = std::make_unique<gmx_localtop_t>(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);
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<int>(rt);
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;
/* 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))
{
* 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");
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;
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.
gmx_sum(2 * dd->nres, dd->Rt_6, cr);
}
- if (fcd->disres.nsystems > 1)
+ if (dd->nsystems > 1)
{
real invn = 1.0 / dd->nsystems;
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;
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];
}
}
}
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;
};
/*! \brief
- * Initiates *fcd data.
+ * Initiates *disresdata.
*
* Must be called once, nbonds is the number
* of iatoms in the ilist of the idef struct.
NumRanks numRanks,
MPI_Comm communicator,
const gmx_multisim_t* ms,
- t_fcdata* fcd,
+ t_disresdata* disresdata,
t_state* state,
gmx_bool bIsREMD);
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.
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
#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"
#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<bonded_threading_t>(numThreads, numEnergyGroups, fplog)),
+ fcdata_(std::make_unique<t_fcdata>())
+{
+}
+
+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
{
#define MAX_BONDED_THREADS 256
/*! \brief Reduce thread-local force buffers */
-void reduce_thread_forces(int n, gmx::ArrayRef<gmx::RVec> force, const bonded_threading_t* bt, int nthreads)
+void reduce_thread_forces(gmx::ArrayRef<gmx::RVec> force, const bonded_threading_t* bt, int nthreads)
{
if (nthreads > MAX_BONDED_THREADS)
{
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)
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++)
}
/*! \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,
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());
/*! \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,
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++)
{
}
else
{
- fshift = threadBuffers.fshift;
+ fshift = as_rvec_array(threadBuffers.fshift.data());
epot = threadBuffers.ener;
grpp = &threadBuffers.grpp;
dvdlt = threadBuffers.dvdl;
{
ArrayRef<const int> 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;
}
}
}
}
-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
/*! \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,
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();
/* 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)
{
/* 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;
}
}
* 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<real> forceBufferLambda,
+ gmx::ArrayRef<gmx::RVec> shiftForceBufferLambda,
gmx_grppairener_t* grpp,
real* epot,
gmx::ArrayRef<real> dvdl,
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;
}
/* 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<rvec4*>(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++)
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<const gmx::RVec> 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<const gmx::RVec> 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())
{
}
/* 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.
{
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++)
and reduction of output data across threads.
*
* \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \author Berk Hess <hess@kth.se>
*
*/
/*! \libinternal \file
* should call functions declared here.
*
* \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \author Berk Hess <hess@kth.se>
*
* \inlibraryapi
* \ingroup module_listed_forces
#ifndef GMX_LISTED_FORCES_LISTED_FORCES_H
#define GMX_LISTED_FORCES_LISTED_FORCES_H
+#include <memory>
+
#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;
//! 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<const gmx::RVec> 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<const gmx::RVec> 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<bonded_threading_t> threading_;
+ //! Data for bonded tables and NMR restraining, needs to be refactored
+ std::unique_ptr<t_fcdata> fcdata_;
+ //! Force buffer for free-energy forces
+ std::vector<real> forceBufferLambda_;
+ //! Shift force buffer for free-energy forces
+ std::vector<gmx::RVec> shiftForceBufferLambda_;
+
+ GMX_DISALLOW_COPY_MOVE_AND_ASSIGN(ListedForces);
+};
#endif
/*
* 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.
* internally by the module.
*
* \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \author Berk Hess <hess@kth.se>
* \ingroup module_listed_forces
*/
#ifndef GMX_LISTED_FORCES_LISTED_INTERNAL_H
#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
//! 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<real, gmx::AlignedAllocator<real>> fBuffer;
+ //! Mask for marking which parts of f are filled, working array for constructing mask in bonded_threading_t
+ std::vector<gmx_bitmask_t> mask;
+ //! Number of blocks touched by our thread
+ int nblock_used = 0;
+ //! Index to touched blocks
+ std::vector<int> 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<gmx::RVec> 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<std::unique_ptr<f_thread_t>> 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<int> 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<gmx_bitmask_t> 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);
};
* interactions.
*
* \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \author Berk Hess <hess@kth.se>
* \ingroup module_listed_forces
*/
#include "gmxpre.h"
#include <string>
#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"
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<rvec4*>(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<gmx_bitmask_t> mask = f_thread->mask;
+
for (int ftype = 0; ftype < F_NRE; ftype++)
{
if (ftype_is_bonded_potential(ftype))
}
void setup_bonded_threading(bonded_threading_t* bt,
- int numAtoms,
+ int numAtomsForce,
bool useGpuForBondeds,
const InteractionDefinitions& idef)
{
assert(bt->nthreads >= 1);
+ bt->numAtomsForce = numAtomsForce;
+
/* Divide the bonded interaction over the threads */
divide_bondeds_over_threads(bt, useGpuForBondeds, idef);
{
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
}
/* 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<size_t>(nblock_tot) > bt->block_index.size())
{
{
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<double>(numAtoms),
+ ctot * reduction_block_size / static_cast<double>(numAtomsForce),
ctot / static_cast<double>(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++)
}
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;
}
* 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 <mark.j.abraham@gmail.com>
- * \inlibraryapi
* \ingroup module_listed_forces
*/
#ifndef GMX_LISTED_FORCES_MANAGE_THREADING_H
* 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
ArrayRef<const RVec> 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)
{
gmx_bool bTAV;
vtot = 0;
- od = &(fcd->orires);
+ od = fcd->orires;
if (od->fc != 0)
{
/* 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];
}
}
}
struct t_inputrec;
struct t_pbc;
struct t_commrec;
-struct t_fcdata;
struct t_oriresdata;
class t_state;
gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
const rvec x[],
const t_pbc* pbc,
- t_fcdata* fcd,
+ t_oriresdata* oriresdata,
history_t* hist);
/*! \brief
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
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];
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;
}
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 */
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);
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,
history_t* hist,
gmx::ForceOutputs* forceOutputs,
gmx_enerdata_t* enerd,
- t_fcdata* fcd,
const matrix box,
const real* lambda,
const rvec* mu_tot,
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,
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))
class InteractionDefinitions;
struct pull_t;
struct t_commrec;
-struct t_fcdata;
struct t_forcerec;
struct t_inputrec;
struct t_lambda;
tensor vir_force,
const t_mdatoms* mdatoms,
gmx_enerdata_t* enerd,
- t_fcdata* fcd,
gmx::ArrayRef<real> lambda,
t_forcerec* fr,
gmx::MdrunScheduleWorkload* runScheduleWork,
*/
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,
history_t* hist,
gmx::ForceOutputs* forceOutputs,
gmx_enerdata_t* enerd,
- t_fcdata* fcd,
const matrix box,
const real* lambda,
const rvec* mu_tot,
#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"
*
* 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<const std::string> tabbfnm,
- const char* tabext)
+static std::vector<bondedtable_t> make_bonded_tables(FILE* fplog,
+ int ftype1,
+ int ftype2,
+ const gmx_mtop_t* mtop,
+ gmx::ArrayRef<const std::string> tabbfnm,
+ const char* tabext)
{
- int ncount, *count;
- bondedtable_t* tab;
+ std::vector<bondedtable_t> 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?
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,
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<ListedForces>(
+ 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
}
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);
/* Note: This code will disappear when types are converted to C++ */
sfree(shift_vec);
sfree(ewc_t);
- tear_down_bonded_threading(bondedThreading);
}
struct gmx_hw_info_t;
struct t_commrec;
-struct t_fcdata;
struct t_forcerec;
struct t_filenm;
struct t_inputrec;
* \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
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,
#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"
/*! \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)
{
// 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
tensor vir_force,
const t_mdatoms* mdatoms,
gmx_enerdata_t* enerd,
- t_fcdata* fcd,
gmx::ArrayRef<real> lambda,
t_forcerec* fr,
gmx::MdrunScheduleWorkload* runScheduleWork,
// 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);
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);
std::vector<char*> 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)
// 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
{
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);
#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"
const t_mdatoms* md,
t_state* state,
const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
- const t_fcdata* fcd,
+ const t_fcdata& fcdata,
const gmx_ekindata_t* ekind,
const matrix M,
int UpdatePart,
const t_mdatoms* md,
t_state* state,
const gmx::ArrayRefWithPadding<const gmx::RVec>& 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);
}
const t_mdatoms* md,
t_state* state,
const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
- const t_fcdata* fcd,
+ const t_fcdata& fcdata,
const gmx_ekindata_t* ekind,
const matrix M,
int updatePart,
/* 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 ######### */
* \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.
const t_mdatoms* md,
t_state* state,
const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
- const t_fcdata* fcd,
+ const t_fcdata& fcdata,
const gmx_ekindata_t* ekind,
const matrix M,
int updatePart,
struct pull_t;
struct ReplicaExchangeParameters;
struct t_commrec;
-struct t_fcdata;
struct t_forcerec;
struct t_filenm;
struct t_inputrec;
pull_t* pull_work,
t_swap* swap,
gmx_mtop_t* top_global,
- t_fcdata* fcd,
t_state* state_global,
ObservablesHistory* observablesHistory,
MDAtoms* mdAtoms,
pull_work(pull_work),
swap(swap),
top_global(top_global),
- fcd(fcd),
state_global(state_global),
observablesHistory(observablesHistory),
mdAtoms(mdAtoms),
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.
#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"
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
"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
{
/* 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,
*/
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);
}
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);
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);
}
}
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);
/* 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);
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)
#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"
#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"
{
/* 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,
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
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))
{
#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"
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.
*/
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));
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
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 */
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. */
{
/* 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);
}
}
* 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))
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.
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. */
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... */
{
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
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);
}
}
/* 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;
/* 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,
#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"
#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"
{
/* 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,
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
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))
{
#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"
{
matrix box;
t_forcerec* fr = nullptr;
- t_fcdata* fcd = nullptr;
real ewaldcoeff_q = 0;
real ewaldcoeff_lj = 0;
int nChargePerturbed = -1, nTypePerturbed = 0;
* 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,
/* 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.
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();
/* 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)
{
}
free_gpu(deviceInfo);
- sfree(fcd);
if (doMembed)
{
const gmx_localtop_t* top,
gmx::Constraints* constr,
gmx_enerdata_t* enerd,
- t_fcdata* fcd,
int natoms,
ArrayRefWithPadding<RVec> xPadded,
ArrayRefWithPadding<RVec> vPadded,
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;
/* 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)
class history_t;
struct pull_t;
struct t_forcerec;
-struct t_fcdata;
struct t_inputrec;
struct t_mdatoms;
struct t_nrnb;
const gmx_localtop_t* top,
gmx::Constraints* constr,
gmx_enerdata_t* enerd,
- t_fcdata* fcd,
int natoms,
gmx::ArrayRefWithPadding<gmx::RVec> x,
gmx::ArrayRefWithPadding<gmx::RVec> v,
struct pull_t;
struct ReplicaExchangeParameters;
struct t_commrec;
-struct t_fcdata;
struct t_forcerec;
struct t_filenm;
struct t_inputrec;
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));
#ifndef GMX_MDTYPES_FCDATA_H
#define GMX_MDTYPES_FCDATA_H
+#include <vector>
+
#include "gromacs/math/vectypes.h"
#include "gromacs/topology/idef.h"
#include "gromacs/utility/basedefinitions.h"
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<real> data; /* the actual table data, per point there are 4 numbers */
+};
/*
* Data struct used in the force calculation routines
* (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<bondedtable_t> bondtab;
+ std::vector<bondedtable_t> angletab;
+ std::vector<bondedtable_t> 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
struct bonded_threading_t;
class DeviceContext;
class DispersionCorrection;
+class ListedForces;
struct t_forcetable;
struct t_QMMMrec;
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> listedForces;
/* TODO: Replace the pointer by an object once we got rid of C */
gmx::GpuBonded* gpuBonded = nullptr;
const MDAtoms* mdAtoms,
t_nrnb* nrnb,
t_forcerec* fr,
- t_fcdata* fcd,
gmx_wallcycle* wcycle,
MdrunScheduleWorkload* runScheduleWork,
VirtualSitesHandler* vsite,
vsite_(vsite),
imdSession_(imdSession),
pull_work_(pull_work),
- fcd_(fcd),
runScheduleWork_(runScheduleWork),
constr_(constr),
enforcedRotation_(enforcedRotation)
relax_shell_flexcon(
fplog_, cr_, ms, isVerbose_, enforcedRotation_, step, inputrec_, imdSession_,
pull_work_, step == nextNSStep_, static_cast<int>(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_++;
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<int>(flags), ddBalanceRegionHandler_);
}
energyElement_->addToForceVirial(force_vir, step);
struct gmx_shellfc_t;
struct gmx_wallcycle;
struct pull_t;
-struct t_fcdata;
struct t_nrnb;
namespace gmx
const MDAtoms* mdAtoms,
t_nrnb* nrnb,
t_forcerec* fr,
- t_fcdata* fcd,
gmx_wallcycle* wcycle,
MdrunScheduleWorkload* runScheduleWork,
VirtualSitesHandler* vsite,
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.
#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"
auto energyElement = std::make_unique<EnergyElement>(
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());
/*
auto forceElement = std::make_unique<ForceElement>(
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()));
int numEnsembleRestraintSystems;
if (fcd)
{
- numEnsembleRestraintSystems = fcd->disres.nsystems;
+ numEnsembleRestraintSystems = fcd->disres->nsystems;
}
else
{
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,
#include "pmeloadbalancehelper.h"
#include "topologyholder.h"
+struct t_fcdata;
+
namespace gmx
{
class DomDecHelper;
}
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;
qB = qA;
}
+ bool havePerturbedInteractions = false;
+
iabuf_nalloc = 0;
iabuf = nullptr;
{
iabuf[ib++] = iatoms[i++];
}
+
+ havePerturbedInteractions = true;
}
else
{
sfree(iabuf);
- idef->ilsort = ilsortFE_SORTED;
+ idef->ilsort = (havePerturbedInteractions ? ilsortFE_SORTED : ilsortNO_FE);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff