Renamed do_force_lowlevel() to calculateLongRangeElectrostatics().
Also replace C pointer for coordinates by ArrayRefWithPadding
in calculate() of ListedForces and do_walls().
} // namespace
-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,
- t_fcdata* fcdata,
- 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)
+void ListedForces::calculate(struct gmx_wallcycle* wcycle,
+ const matrix box,
+ const t_lambda* fepvals,
+ const t_commrec* cr,
+ const gmx_multisim_t* ms,
+ gmx::ArrayRefWithPadding<const gmx::RVec> coordinates,
+ gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
+ t_fcdata* fcdata,
+ 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 (interactionSelection_.none() || !stepWork.computeListedForces)
{
const InteractionDefinitions& idef = *idef_;
+ // Todo: replace all rvec use here with ArrayRefWithPadding
+ const rvec* x = as_rvec_array(coordinates.paddedArrayRef().data());
+
t_pbc pbc_full; /* Full PBC is needed for position restraints */
if (haveRestraints(*fcdata))
{
#define GMX_LISTED_FORCES_LISTED_FORCES_H
#include <memory>
+#include <vector>
#include <bitset>
#include "gromacs/math/vectypes.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/ifunc.h"
-#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/classhelpers.h"
class StepWorkload;
template<typename>
class ArrayRef;
+template<typename>
+class ArrayRefWithPadding;
} // namespace gmx
//! Type of CPU function to compute a bonded interaction.
* 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,
- t_fcdata* fcdata,
- 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);
+ void calculate(struct gmx_wallcycle* wcycle,
+ const matrix box,
+ const t_lambda* fepvals,
+ const t_commrec* cr,
+ const gmx_multisim_t* ms,
+ gmx::ArrayRefWithPadding<const gmx::RVec> coordinates,
+ gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
+ t_fcdata* fcdata,
+ 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;
#include "gromacs/ewald/pme.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/math/vecdump.h"
#include "gromacs/mdlib/forcerec_threading.h"
-#include "gromacs/mdlib/rf_util.h"
-#include "gromacs/mdlib/wall.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/enerdata.h"
#include "gromacs/mdtypes/forceoutput.h"
}
}
-void do_force_lowlevel(t_forcerec* fr,
- const t_inputrec* ir,
- const t_commrec* cr,
- const gmx_multisim_t* ms,
- t_nrnb* nrnb,
- gmx_wallcycle_t wcycle,
- const t_mdatoms* md,
- gmx::ArrayRefWithPadding<const RVec> coordinates,
- ArrayRef<const RVec> xWholeMolecules,
- history_t* hist,
- gmx::ForceOutputs* forceOutputs,
- gmx_enerdata_t* enerd,
- const matrix box,
- const real* lambda,
- const rvec* mu_tot,
- const gmx::StepWorkload& stepWork,
- const DDBalanceRegionHandler& ddBalanceRegionHandler)
+void calculateLongRangeNonbondeds(t_forcerec* fr,
+ const t_inputrec* ir,
+ const t_commrec* cr,
+ t_nrnb* nrnb,
+ gmx_wallcycle_t wcycle,
+ const t_mdatoms* md,
+ gmx::ArrayRef<const RVec> coordinates,
+ gmx::ForceWithVirial* forceWithVirial,
+ gmx_enerdata_t* enerd,
+ const matrix box,
+ const real* lambda,
+ const rvec* mu_tot,
+ const gmx::StepWorkload& stepWork,
+ const DDBalanceRegionHandler& ddBalanceRegionHandler)
{
- // TODO: Replace all uses of x by const coordinates
- const rvec* x = as_rvec_array(coordinates.paddedArrayRef().data());
-
- auto& forceWithVirial = forceOutputs->forceWithVirial();
-
- /* Call the short range functions all in one go. */
-
- if (ir->nwall)
- {
- /* foreign lambda component for walls */
- real dvdl_walls = do_walls(*ir, *fr, box, *md, x, &forceWithVirial, lambda[efptVDW],
- enerd->grpp.ener[egLJSR].data(), nrnb);
- enerd->dvdl_lin[efptVDW] += dvdl_walls;
- }
-
- {
- t_pbc pbc;
-
- /* Check whether we need to take into account PBC in listed interactions. */
- ListedForces& listedForces = fr->listedForces[0];
- const auto needPbcForListedForces = fr->bMolPBC && stepWork.computeListedForces
- && listedForces.haveCpuListedForces(*fr->fcdata);
- if (needPbcForListedForces)
- {
- /* Since all atoms are in the rectangular or triclinic unit-cell,
- * only single box vector shifts (2 in x) are required.
- */
- set_pbc_dd(&pbc, fr->pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : nullptr, TRUE, box);
- }
-
- listedForces.calculate(wcycle, box, ir->fepvals, cr, ms, x, xWholeMolecules, fr->fcdata.get(),
- 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))
&& thisRankHasDuty(cr, DUTY_PME)
&& (pme_run_mode(fr->pmedata) == PmeRunMode::CPU);
* exclusion forces) are calculated, so we can store
* the forces in the normal, single forceWithVirial->force_ array.
*/
+ const rvec* x = as_rvec_array(coordinates.data());
ewald_LRcorrection(md->homenr, cr, nthreads, t, *fr, *ir, md->chargeA,
md->chargeB, (md->nChargePerturbed != 0), x, box, mu_tot,
- as_rvec_array(forceWithVirial.force_.data()),
+ as_rvec_array(forceWithVirial->force_.data()),
&ewc_t.Vcorr_q, lambda[efptCOUL], &ewc_t.dvdl[efptCOUL]);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
wallcycle_start(wcycle, ewcPMEMESH);
status = gmx_pme_do(
fr->pmedata,
- gmx::constArrayRefFromArray(coordinates.unpaddedConstArrayRef().data(),
- md->homenr - fr->n_tpi),
- forceWithVirial.force_, md->chargeA, md->chargeB, md->sqrt_c6A,
+ gmx::constArrayRefFromArray(coordinates.data(), md->homenr - fr->n_tpi),
+ forceWithVirial->force_, md->chargeA, md->chargeB, md->sqrt_c6A,
md->sqrt_c6B, md->sigmaA, md->sigmaB, box, cr,
DOMAINDECOMP(cr) ? dd_pme_maxshift_x(cr->dd) : 0,
DOMAINDECOMP(cr) ? dd_pme_maxshift_y(cr->dd) : 0, nrnb, wcycle,
* with the PME grid potential of the other charges.
*/
gmx_pme_calc_energy(
- fr->pmedata,
- coordinates.unpaddedConstArrayRef().subArray(md->homenr - fr->n_tpi, fr->n_tpi),
+ fr->pmedata, coordinates.subArray(md->homenr - fr->n_tpi, fr->n_tpi),
gmx::arrayRefFromArray(md->chargeA + md->homenr - fr->n_tpi, fr->n_tpi),
&Vlr_q);
}
if (fr->ic->eeltype == eelEWALD)
{
- Vlr_q = do_ewald(ir, x, as_rvec_array(forceWithVirial.force_.data()), md->chargeA,
+ const rvec* x = as_rvec_array(coordinates.data());
+ Vlr_q = do_ewald(ir, x, as_rvec_array(forceWithVirial->force_.data()), md->chargeA,
md->chargeB, box, cr, md->homenr, ewaldOutput.vir_q, fr->ic->ewaldcoeff_q,
lambda[efptCOUL], &ewaldOutput.dvdl[efptCOUL], fr->ewald_table);
}
/* Note that with separate PME nodes we get the real energies later */
// TODO it would be simpler if we just accumulated a single
// long-range virial contribution.
- forceWithVirial.addVirialContribution(ewaldOutput.vir_q);
- forceWithVirial.addVirialContribution(ewaldOutput.vir_lj);
+ forceWithVirial->addVirialContribution(ewaldOutput.vir_q);
+ forceWithVirial->addVirialContribution(ewaldOutput.vir_lj);
enerd->dvdl_lin[efptCOUL] += ewaldOutput.dvdl[efptCOUL];
enerd->dvdl_lin[efptVDW] += ewaldOutput.dvdl[efptVDW];
enerd->term[F_COUL_RECIP] = Vlr_q + ewaldOutput.Vcorr_q;
fprintf(debug, "Vlr_q = %g, Vcorr_q = %g, Vlr_corr_q = %g\n", Vlr_q,
ewaldOutput.Vcorr_q, enerd->term[F_COUL_RECIP]);
pr_rvecs(debug, 0, "vir_el_recip after corr", ewaldOutput.vir_q, DIM);
- rvec* fshift = as_rvec_array(forceOutputs->forceWithShiftForces().shiftForces().data());
- pr_rvecs(debug, 0, "fshift after LR Corrections", fshift, SHIFTS);
fprintf(debug, "Vlr_lj: %g, Vcorr_lj = %g, Vlr_corr_lj = %g\n", Vlr_lj,
ewaldOutput.Vcorr_lj, enerd->term[F_LJ_RECIP]);
pr_rvecs(debug, 0, "vir_lj_recip after corr", ewaldOutput.vir_lj, DIM);
{
print_nrnb(debug, nrnb);
}
-
- if (debug)
- {
- rvec* fshift = as_rvec_array(forceOutputs->forceWithShiftForces().shiftForces().data());
- pr_rvecs(debug, 0, "fshift after bondeds", fshift, SHIFTS);
- }
}
#ifndef GMX_MDLIB_FORCE_H
#define GMX_MDLIB_FORCE_H
-#include "gromacs/math/arrayrefwithpadding.h"
+#include <cstdio>
+
#include "gromacs/math/vectypes.h"
-#include "gromacs/utility/arrayref.h"
class DDBalanceRegionHandler;
struct gmx_edsam;
namespace gmx
{
+template<typename>
+class ArrayRef;
+template<typename>
+class ArrayRefWithPadding;
class Awh;
class ForceBuffersView;
-class ForceOutputs;
class ForceWithVirial;
class ImdSession;
class MdrunScheduleWorkload;
*/
-/* Compute listed forces, Ewald, PME corrections add when (when used).
+/* Calculate CPU Ewald or PME-mesh forces when done on this rank and Ewald corrections, when used
*
- * xWholeMolecules only needs to contain whole molecules when orientation
- * restraints need to be computed and can be empty otherwise.
+ * Note that Ewald dipole and net charge corrections are always computed here, independently
+ * on whether the PME-mesh contribution is computed on a separate PME rank or on a GPU.
*/
-void do_force_lowlevel(t_forcerec* fr,
- const t_inputrec* ir,
- const t_commrec* cr,
- const gmx_multisim_t* ms,
- t_nrnb* nrnb,
- gmx_wallcycle* wcycle,
- const t_mdatoms* md,
- gmx::ArrayRefWithPadding<const gmx::RVec> coordinates,
- gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
- history_t* hist,
- gmx::ForceOutputs* forceOutputs,
- gmx_enerdata_t* enerd,
- const matrix box,
- const real* lambda,
- const rvec* mu_tot,
- const gmx::StepWorkload& stepWork,
- const DDBalanceRegionHandler& ddBalanceRegionHandler);
-/* Call all the force routines */
+void calculateLongRangeNonbondeds(t_forcerec* fr,
+ const t_inputrec* ir,
+ const t_commrec* cr,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle,
+ const t_mdatoms* md,
+ gmx::ArrayRef<const gmx::RVec> coordinates,
+ gmx::ForceWithVirial* forceWithVirial,
+ gmx_enerdata_t* enerd,
+ const matrix box,
+ const real* lambda,
+ const rvec* mu_tot,
+ const gmx::StepWorkload& stepWork,
+ const DDBalanceRegionHandler& ddBalanceRegionHandler);
#endif
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/update.h"
#include "gromacs/mdlib/vsite.h"
+#include "gromacs/mdlib/wall.h"
#include "gromacs/mdlib/wholemoleculetransform.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/enerdata.h"
/* Wait for non-local coordinate data to be copied from device */
stateGpu->waitCoordinatesReadyOnHost(AtomLocality::NonLocal);
}
- /* Compute the bonded and non-bonded energies and optionally forces */
- 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);
+
+ // Compute wall interactions, when present.
+ // Note: should be moved to special forces.
+ if (inputrec->nwall && stepWork.computeNonbondedForces)
+ {
+ /* foreign lambda component for walls */
+ real dvdl_walls = do_walls(*inputrec, *fr, box, *mdatoms, x.unpaddedConstArrayRef(),
+ &forceOut.forceWithVirial(), lambda[efptVDW],
+ enerd->grpp.ener[egLJSR].data(), nrnb);
+ enerd->dvdl_lin[efptVDW] += dvdl_walls;
+ }
+
+ if (stepWork.computeListedForces)
+ {
+ /* Check whether we need to take into account PBC in listed interactions */
+ bool needMolPbc = false;
+ for (const auto& listedForces : fr->listedForces)
+ {
+ if (listedForces.haveCpuListedForces(*fr->fcdata))
+ {
+ needMolPbc = fr->bMolPBC;
+ }
+ }
+
+ t_pbc pbc;
+
+ if (needMolPbc)
+ {
+ /* Since all atoms are in the rectangular or triclinic unit-cell,
+ * only single box vector shifts (2 in x) are required.
+ */
+ set_pbc_dd(&pbc, fr->pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : nullptr, TRUE, box);
+ }
+
+ for (auto& listedForces : fr->listedForces)
+ {
+ listedForces.calculate(
+ wcycle, box, inputrec->fepvals, cr, ms, x, xWholeMolecules, fr->fcdata.get(),
+ hist, &forceOut, fr, &pbc, enerd, nrnb, lambda.data(), mdatoms,
+ DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr, stepWork);
+ }
+ }
+
+ calculateLongRangeNonbondeds(fr, inputrec, cr, nrnb, wcycle, mdatoms, x.unpaddedConstArrayRef(),
+ &forceOut.forceWithVirial(), enerd, box, lambda.data(),
+ as_rvec_array(dipoleData.muStateAB), stepWork, ddBalanceRegionHandler);
wallcycle_stop(wcycle, ewcFORCE);
#include "gromacs/mdtypes/nblist.h"
#include "gromacs/tables/forcetable.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
}
}
-[[noreturn]] static void wall_error(int a, const rvec* x, real r)
+[[noreturn]] static void wall_error(int a, gmx::ArrayRef<const gmx::RVec> x, real r)
{
gmx_fatal(FARGS,
"An atom is beyond the wall: coordinates %f %f %f, distance %f\n"
}
}
-real do_walls(const t_inputrec& ir,
- const t_forcerec& fr,
- const matrix box,
- const t_mdatoms& md,
- const rvec* x,
- gmx::ForceWithVirial* forceWithVirial,
- real lambda,
- real Vlj[],
- t_nrnb* nrnb)
+real do_walls(const t_inputrec& ir,
+ const t_forcerec& fr,
+ const matrix box,
+ const t_mdatoms& md,
+ gmx::ArrayRef<const gmx::RVec> x,
+ gmx::ForceWithVirial* forceWithVirial,
+ real lambda,
+ real Vlj[],
+ t_nrnb* nrnb)
{
constexpr real sixth = 1.0 / 6.0;
constexpr real twelfth = 1.0 / 12.0;
namespace gmx
{
+template<typename>
+class ArrayRef;
class ForceWithVirial;
-}
+} // namespace gmx
void make_wall_tables(FILE* fplog,
const t_inputrec* ir,
const SimulationGroups* groups,
t_forcerec* fr);
-real do_walls(const t_inputrec& ir,
- const t_forcerec& fr,
- const matrix box,
- const t_mdatoms& md,
- const rvec x[],
- gmx::ForceWithVirial* forceWithVirial,
- real lambda,
- real Vlj[],
- t_nrnb* nrnb);
+real do_walls(const t_inputrec& ir,
+ const t_forcerec& fr,
+ const matrix box,
+ const t_mdatoms& md,
+ gmx::ArrayRef<const gmx::RVec> x,
+ gmx::ForceWithVirial* forceWithVirial,
+ real lambda,
+ real Vlj[],
+ t_nrnb* nrnb);
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff