#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
-#include "gromacs/mdlib/forcerec_threading.h"
+#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/enerdata.h"
#include "gromacs/mdtypes/forceoutput.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/simulation_workload.h"
-#include "gromacs/pbcutil/ishift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
using gmx::ArrayRef;
using gmx::RVec;
}
}
-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 RVec> coordinates,
- gmx::ForceWithVirial* forceWithVirial,
- gmx_enerdata_t* enerd,
- const matrix box,
- gmx::ArrayRef<const real> lambda,
- gmx::ArrayRef<const gmx::RVec> mu_tot,
- const gmx::StepWorkload& stepWork,
- const DDBalanceRegionHandler& ddBalanceRegionHandler)
+CpuPpLongRangeNonbondeds::CpuPpLongRangeNonbondeds(int numberOfTestPaticles,
+ real ewaldCoeffQ,
+ real epsilonR,
+ gmx::ArrayRef<const double> chargeC6Sum,
+ CoulombInteractionType eeltype,
+ VanDerWaalsType vdwtype,
+ const t_inputrec& inputrec,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle,
+ FILE* fplog) :
+ numTpiAtoms_(numberOfTestPaticles),
+ ewaldCoeffQ_(ewaldCoeffQ),
+ epsilonR_(epsilonR),
+ chargeC6Sum_(chargeC6Sum),
+ coulombInteractionType_(eeltype),
+ vanDerWaalsType_(vdwtype),
+ ewaldGeometry_(inputrec.ewald_geometry),
+ epsilonSurface_(inputrec.epsilon_surface),
+ haveEwaldSurfaceTerm_(haveEwaldSurfaceContribution(inputrec)),
+ wallEwaldZfac_(inputrec.wall_ewald_zfac),
+ havePbcXY2Walls_(inputrecPbcXY2Walls(&inputrec)),
+ freeEnergyPerturbationType_(inputrec.efep),
+ nrnb_(nrnb),
+ wcycle_(wcycle)
{
- const bool computePmeOnCpu = (EEL_PME(fr->ic->eeltype) || EVDW_PME(fr->ic->vdwtype))
- && thisRankHasDuty(cr, DUTY_PME)
- && (pme_run_mode(fr->pmedata) == PmeRunMode::CPU);
+ GMX_ASSERT(epsilonR == inputrec.epsilon_r,
+ "Forcerec and inputrec relative dielectrics don't match");
- const bool haveEwaldSurfaceTerm = haveEwaldSurfaceContribution(ir);
+ // Use the thread count for the bonded module since reducing CPU-side
+ // non-bonded contributions does not currently have its own thread count.
+ outputPerThread_.resize(gmx_omp_nthreads_get(ModuleMultiThread::Bonded));
+
+ if (inputrec.coulombtype == CoulombInteractionType::Ewald)
+ {
+ ewaldTable_ = std::make_unique<gmx_ewald_tab_t>(inputrec, fplog);
+ }
+}
+
+CpuPpLongRangeNonbondeds::~CpuPpLongRangeNonbondeds() = default;
+
+void CpuPpLongRangeNonbondeds::updateAfterPartition(const t_mdatoms& md)
+{
+ homenr_ = md.homenr;
+ havePerturbed_ = md.nChargePerturbed != 0;
+ chargeA_ = md.chargeA ? gmx::constArrayRefFromArray(md.chargeA, md.nr) : gmx::ArrayRef<const real>{};
+ chargeB_ = md.chargeB ? gmx::constArrayRefFromArray(md.chargeB, md.nr) : gmx::ArrayRef<const real>{};
+ sqrt_c6A_ = md.sqrt_c6A ? gmx::constArrayRefFromArray(md.sqrt_c6A, md.nr)
+ : gmx::ArrayRef<const real>{};
+ sqrt_c6B_ = md.sqrt_c6B ? gmx::constArrayRefFromArray(md.sqrt_c6B, md.nr)
+ : gmx::ArrayRef<const real>{};
+ sigmaA_ = md.sigmaA ? gmx::constArrayRefFromArray(md.sigmaA, md.nr) : gmx::ArrayRef<const real>{};
+ sigmaB_ = md.sigmaB ? gmx::constArrayRefFromArray(md.sigmaB, md.nr) : gmx::ArrayRef<const real>{};
+}
+
+void CpuPpLongRangeNonbondeds::calculate(gmx_pme_t* pmedata,
+ const t_commrec* commrec,
+ gmx::ArrayRef<const RVec> coordinates,
+ gmx::ForceWithVirial* forceWithVirial,
+ gmx_enerdata_t* enerd,
+ const matrix box,
+ gmx::ArrayRef<const real> lambda,
+ gmx::ArrayRef<const gmx::RVec> mu_tot,
+ const gmx::StepWorkload& stepWork,
+ const DDBalanceRegionHandler& ddBalanceRegionHandler)
+{
+ const bool computePmeOnCpu = (EEL_PME(coulombInteractionType_) || EVDW_PME(vanDerWaalsType_))
+ && thisRankHasDuty(commrec, DUTY_PME)
+ && (pme_run_mode(pmedata) == PmeRunMode::CPU);
/* Do long-range electrostatics and/or LJ-PME
* and compute PME surface terms when necessary.
*/
- if ((computePmeOnCpu || fr->ic->eeltype == CoulombInteractionType::Ewald || haveEwaldSurfaceTerm)
+ if ((computePmeOnCpu || coulombInteractionType_ == CoulombInteractionType::Ewald || haveEwaldSurfaceTerm_)
&& stepWork.computeNonbondedForces)
{
- int status = 0;
real Vlr_q = 0, Vlr_lj = 0;
-
/* We reduce all virial, dV/dlambda and energy contributions, except
* for the reciprocal energies (Vlr_q, Vlr_lj) into the same struct.
*/
- ewald_corr_thread_t& ewaldOutput = fr->ewc_t[0];
+ ewald_corr_thread_t& ewaldOutput = outputPerThread_[0];
clearEwaldThreadOutput(&ewaldOutput);
- if (EEL_PME_EWALD(fr->ic->eeltype) || EVDW_PME(fr->ic->vdwtype))
+ if (EEL_PME_EWALD(coulombInteractionType_) || EVDW_PME(vanDerWaalsType_))
{
/* Calculate the Ewald surface force and energy contributions, when necessary */
- if (haveEwaldSurfaceTerm)
+ if (haveEwaldSurfaceTerm_)
{
- wallcycle_sub_start(wcycle, WallCycleSubCounter::EwaldCorrection);
+ wallcycle_sub_start(wcycle_, WallCycleSubCounter::EwaldCorrection);
- int nthreads = fr->nthread_ewc;
+ int nthreads = gmx::ssize(outputPerThread_);
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (int t = 0; t < nthreads; t++)
{
try
{
- ewald_corr_thread_t& ewc_t = fr->ewc_t[t];
+ ewald_corr_thread_t& ewc_t = outputPerThread_[t];
if (t > 0)
{
clearEwaldThreadOutput(&ewc_t);
* the forces in the normal, single forceWithVirial->force_ array.
*/
ewald_LRcorrection(
- md->homenr,
- cr,
+ homenr_,
+ commrec,
nthreads,
t,
- fr->ic->epsilon_r,
- fr->qsum,
- ir.ewald_geometry,
- ir.epsilon_surface,
- inputrecPbcXY2Walls(&ir),
- ir.wall_ewald_zfac,
- md->chargeA ? gmx::constArrayRefFromArray(md->chargeA, md->nr)
- : gmx::ArrayRef<const real>{},
- md->chargeB ? gmx::constArrayRefFromArray(md->chargeB, md->nr)
- : gmx::ArrayRef<const real>{},
- (md->nChargePerturbed != 0),
+ epsilonR_,
+ chargeC6Sum_,
+ ewaldGeometry_,
+ epsilonSurface_,
+ havePbcXY2Walls_,
+ wallEwaldZfac_,
+ chargeA_,
+ chargeB_,
+ havePerturbed_,
coordinates,
box,
mu_tot,
}
if (nthreads > 1)
{
- reduceEwaldThreadOuput(nthreads, fr->ewc_t);
+ reduceEwaldThreadOuput(nthreads, outputPerThread_);
}
- wallcycle_sub_stop(wcycle, WallCycleSubCounter::EwaldCorrection);
+ wallcycle_sub_stop(wcycle_, WallCycleSubCounter::EwaldCorrection);
}
- if (EEL_PME_EWALD(fr->ic->eeltype) && fr->n_tpi == 0)
+ if (EEL_PME_EWALD(coulombInteractionType_) && numTpiAtoms_ == 0)
{
/* This is not in a subcounter because it takes a
negligible and constant-sized amount of time */
ewaldOutput.Vcorr_q += ewald_charge_correction(
- cr,
- fr->ic->epsilon_r,
- fr->ic->ewaldcoeff_q,
- fr->qsum,
+ commrec,
+ epsilonR_,
+ ewaldCoeffQ_,
+ chargeC6Sum_,
lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)],
box,
&ewaldOutput.dvdl[FreeEnergyPerturbationCouplingType::Coul],
if (computePmeOnCpu)
{
/* Do reciprocal PME for Coulomb and/or LJ. */
- assert(fr->n_tpi >= 0);
- if (fr->n_tpi == 0 || stepWork.stateChanged)
+ assert(numTpiAtoms_ >= 0);
+ if (numTpiAtoms_ == 0 || stepWork.stateChanged)
{
/* With domain decomposition we close the CPU side load
* balancing region here, because PME does global
*/
ddBalanceRegionHandler.closeAfterForceComputationCpu();
- wallcycle_start(wcycle, WallCycleCounter::PmeMesh);
- status = gmx_pme_do(
- fr->pmedata,
- gmx::constArrayRefFromArray(coordinates.data(), md->homenr - fr->n_tpi),
+ wallcycle_start(wcycle_, WallCycleCounter::PmeMesh);
+ int status = gmx_pme_do(
+ pmedata,
+ gmx::constArrayRefFromArray(coordinates.data(), homenr_ - numTpiAtoms_),
forceWithVirial->force_,
- md->chargeA ? gmx::constArrayRefFromArray(md->chargeA, md->nr)
- : gmx::ArrayRef<const real>{},
- md->chargeB ? gmx::constArrayRefFromArray(md->chargeB, md->nr)
- : gmx::ArrayRef<const real>{},
- md->sqrt_c6A ? gmx::constArrayRefFromArray(md->sqrt_c6A, md->nr)
- : gmx::ArrayRef<const real>{},
- md->sqrt_c6B ? gmx::constArrayRefFromArray(md->sqrt_c6B, md->nr)
- : gmx::ArrayRef<const real>{},
- md->sigmaA ? gmx::constArrayRefFromArray(md->sigmaA, md->nr)
- : gmx::ArrayRef<const real>{},
- md->sigmaB ? gmx::constArrayRefFromArray(md->sigmaB, md->nr)
- : gmx::ArrayRef<const real>{},
+ chargeA_,
+ chargeB_,
+ sqrt_c6A_,
+ sqrt_c6B_,
+ sigmaA_,
+ sigmaB_,
box,
- cr,
- DOMAINDECOMP(cr) ? dd_pme_maxshift_x(*cr->dd) : 0,
- DOMAINDECOMP(cr) ? dd_pme_maxshift_y(*cr->dd) : 0,
- nrnb,
- wcycle,
+ commrec,
+ DOMAINDECOMP(commrec) ? dd_pme_maxshift_x(*commrec->dd) : 0,
+ DOMAINDECOMP(commrec) ? dd_pme_maxshift_y(*commrec->dd) : 0,
+ nrnb_,
+ wcycle_,
ewaldOutput.vir_q,
ewaldOutput.vir_lj,
&Vlr_q,
&ewaldOutput.dvdl[FreeEnergyPerturbationCouplingType::Coul],
&ewaldOutput.dvdl[FreeEnergyPerturbationCouplingType::Vdw],
stepWork);
- wallcycle_stop(wcycle, WallCycleCounter::PmeMesh);
+ wallcycle_stop(wcycle_, WallCycleCounter::PmeMesh);
if (status != 0)
{
gmx_fatal(FARGS, "Error %d in reciprocal PME routine", status);
* of the force call (without PME).
*/
}
- if (fr->n_tpi > 0)
+ if (numTpiAtoms_ > 0)
{
/* Determine the PME grid energy of the test molecule
* with the PME grid potential of the other charges.
*/
Vlr_q = gmx_pme_calc_energy(
- fr->pmedata,
- coordinates.subArray(md->homenr - fr->n_tpi, fr->n_tpi),
- gmx::arrayRefFromArray(md->chargeA + md->homenr - fr->n_tpi, fr->n_tpi));
+ pmedata,
+ coordinates.subArray(homenr_ - numTpiAtoms_, numTpiAtoms_),
+ chargeA_.subArray(homenr_ - numTpiAtoms_, numTpiAtoms_));
}
}
}
- if (fr->ic->eeltype == CoulombInteractionType::Ewald)
+ if (coulombInteractionType_ == CoulombInteractionType::Ewald)
{
- Vlr_q = do_ewald(inputrecPbcXY2Walls(&ir),
- ir.wall_ewald_zfac,
- ir.epsilon_r,
- ir.efep,
+ Vlr_q = do_ewald(havePbcXY2Walls_,
+ wallEwaldZfac_,
+ epsilonR_,
+ freeEnergyPerturbationType_,
coordinates,
forceWithVirial->force_,
- gmx::arrayRefFromArray(md->chargeA, md->nr),
- gmx::arrayRefFromArray(md->chargeB, md->nr),
+ chargeA_,
+ chargeB_,
box,
- cr,
- md->homenr,
+ commrec,
+ homenr_,
ewaldOutput.vir_q,
- fr->ic->ewaldcoeff_q,
+ ewaldCoeffQ_,
lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)],
&ewaldOutput.dvdl[FreeEnergyPerturbationCouplingType::Coul],
- fr->ewald_table.get());
+ ewaldTable_.get());
}
/* Note that with separate PME nodes we get the real energies later */
if (debug)
{
- print_nrnb(debug, nrnb);
+ print_nrnb(debug, nrnb_);
}
}
#include <cstdio>
+#include <memory>
+
#include "gromacs/math/vectypes.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/enumerationhelpers.h"
class DDBalanceRegionHandler;
struct gmx_edsam;
struct gmx_localtop_t;
struct gmx_multisim_t;
struct gmx_wallcycle;
+struct gmx_pme_t;
class history_t;
class InteractionDefinitions;
struct pull_t;
struct t_lambda;
struct t_mdatoms;
struct t_nrnb;
+struct gmx_ewald_tab_t;
+class CpuPpLongRangeNonbondeds;
namespace gmx
{
template<typename>
-class ArrayRef;
-template<typename>
class ArrayRefWithPadding;
class Awh;
class ForceBuffersView;
class VirtualSitesHandler;
} // namespace gmx
+struct ewald_corr_thread_t
+{
+ real Vcorr_q;
+ real Vcorr_lj;
+ gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> dvdl;
+ tensor vir_q;
+ tensor vir_lj;
+};
+
/* Perform the force and, if requested, energy computation
*
* Without multiple time stepping the force is returned in force->force().
rvec mu_tot,
double t,
gmx_edsam* ed,
+ CpuPpLongRangeNonbondeds* longRangeNonbondeds,
int legacyFlags,
const DDBalanceRegionHandler& ddBalanceRegionHandler);
*/
-/* Calculate CPU Ewald or PME-mesh forces when done on this rank and Ewald corrections, when used
- *
- * 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 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,
- gmx::ArrayRef<const real> lambda,
- gmx::ArrayRef<const gmx::RVec> mu_tot,
- const gmx::StepWorkload& stepWork,
- const DDBalanceRegionHandler& ddBalanceRegionHandler);
+class CpuPpLongRangeNonbondeds
+{
+public:
+ /* \brief Constructor
+ *
+ * Should be called after init_forcerec if params come from a populated forcerec
+ */
+ CpuPpLongRangeNonbondeds(int numberOfTestPaticles,
+ real ewaldCoeffQ,
+ real epsilonR,
+ gmx::ArrayRef<const double> chargeC6Sum,
+ CoulombInteractionType eeltype,
+ VanDerWaalsType vdwtype,
+ const t_inputrec& inputrec,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle,
+ FILE* fplog);
+
+ ~CpuPpLongRangeNonbondeds();
+
+ void updateAfterPartition(const t_mdatoms& md);
+
+ /* Calculate CPU Ewald or PME-mesh forces when done on this rank and Ewald corrections, when used
+ *
+ * Note that Ewald dipole and net charge corrections are always computed here, independently
+ * of whether the PME-mesh contribution is computed on a separate PME rank or on a GPU.
+ */
+ void calculate(gmx_pme_t* pmedata,
+ const t_commrec* commrec,
+ gmx::ArrayRef<const gmx::RVec> coordinates,
+ gmx::ForceWithVirial* forceWithVirial,
+ gmx_enerdata_t* enerd,
+ const matrix box,
+ gmx::ArrayRef<const real> lambda,
+ gmx::ArrayRef<const gmx::RVec> mu_tot,
+ const gmx::StepWorkload& stepWork,
+ const DDBalanceRegionHandler& ddBalanceRegionHandler);
+
+private:
+ //! Number of particles for test particle insertion
+ int numTpiAtoms_;
+ //! Ewald charge coefficient
+ real ewaldCoeffQ_;
+ //! Dielectric constant
+ real epsilonR_;
+ //! [0]: sum of charges; [1]: sum of C6's
+ gmx::ArrayRef<const double> chargeC6Sum_;
+ //! Cut-off treatment for Coulomb
+ CoulombInteractionType coulombInteractionType_;
+ //! Van der Waals interaction treatment
+ VanDerWaalsType vanDerWaalsType_;
+ //! Ewald geometry
+ EwaldGeometry ewaldGeometry_;
+ //! Epsilon for PME dipole correction
+ real epsilonSurface_;
+ //! Whether a long range correction is used
+ bool haveEwaldSurfaceTerm_;
+ //! Scaling factor for the box for Ewald
+ real wallEwaldZfac_;
+ //! Whether the simulation is 2D periodic with two walls
+ bool havePbcXY2Walls_;
+ //! Free energy perturbation type
+ FreeEnergyPerturbationType freeEnergyPerturbationType_;
+ //! Number of atoms on this node
+ int homenr_;
+ //! Whether there are perturbed interactions
+ bool havePerturbed_;
+ //! State A charge
+ gmx::ArrayRef<const real> chargeA_;
+ //! State B charge
+ gmx::ArrayRef<const real> chargeB_;
+ //! State A LJ c6
+ gmx::ArrayRef<const real> sqrt_c6A_;
+ //! State B LJ c6
+ gmx::ArrayRef<const real> sqrt_c6B_;
+ //! State A LJ sigma
+ gmx::ArrayRef<const real> sigmaA_;
+ //! State B LJ sigma
+ gmx::ArrayRef<const real> sigmaB_;
+ //! Ewald correction thread local virial and energy data
+ std::vector<ewald_corr_thread_t> outputPerThread_;
+ //! Ewald table
+ std::unique_ptr<gmx_ewald_tab_t> ewaldTable_;
+ //! Non bonded kernel flop counters
+ t_nrnb* nrnb_;
+ //! Wall cycle counters
+ gmx_wallcycle* wcycle_;
+};
#endif
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/dispersioncorrection.h"
#include "gromacs/mdlib/force.h"
-#include "gromacs/mdlib/forcerec_threading.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/wall.h"
const interaction_const_t* interactionConst = forcerec->ic.get();
- /* TODO: Replace this Ewald table or move it into interaction_const_t */
- if (inputrec.coulombtype == CoulombInteractionType::Ewald)
- {
- forcerec->ewald_table = std::make_unique<gmx_ewald_tab_t>(inputrec, fplog);
- }
-
/* Electrostatics: Translate from interaction-setting-in-mdp-file to kernel interaction format */
switch (interactionConst->eeltype)
{
forcerec->print_force = print_force;
- forcerec->nthread_ewc = gmx_omp_nthreads_get(ModuleMultiThread::Bonded);
- forcerec->ewc_t.resize(forcerec->nthread_ewc);
-
if (inputrec.eDispCorr != DispersionCorrectionType::No)
{
forcerec->dispersionCorrection = std::make_unique<DispersionCorrection>(
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2014,2015,2018,2019,2021, 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.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-
-#ifndef GMX_MDLIB_FORCEREC_THREADING_H
-#define GMX_MDLIB_FORCEREC_THREADING_H
-
-#include "gromacs/math/vectypes.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/utility/enumerationhelpers.h"
-
-struct ewald_corr_thread_t
-{
- real Vcorr_q;
- real Vcorr_lj;
- gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> dvdl;
- tensor vir_q;
- tensor vir_lj;
-};
-
-#endif
rvec muTotal,
double t,
gmx_edsam* ed,
+ CpuPpLongRangeNonbondeds* longRangeNonbondeds,
int legacyFlags,
const DDBalanceRegionHandler& ddBalanceRegionHandler)
{
if (stepWork.computeSlowForces)
{
- calculateLongRangeNonbondeds(fr,
- inputrec,
- cr,
- nrnb,
- wcycle,
- mdatoms,
- x.unpaddedConstArrayRef(),
- &forceOutMtsLevel1->forceWithVirial(),
- enerd,
- box,
- lambda,
- dipoleData.muStateAB,
- stepWork,
- ddBalanceRegionHandler);
+ longRangeNonbondeds->calculate(fr->pmedata,
+ cr,
+ x.unpaddedConstArrayRef(),
+ &forceOutMtsLevel1->forceWithVirial(),
+ enerd,
+ box,
+ lambda,
+ dipoleData.muStateAB,
+ stepWork,
+ ddBalanceRegionHandler);
}
wallcycle_stop(wcycle, WallCycleCounter::Force);
: gmx::ArrayRef<const unsigned short>(),
md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
: gmx::ArrayRef<const unsigned short>());
+ fr->longRangeNonbondeds->updateAfterPartition(*md);
}
else
{
: gmx::ArrayRef<const unsigned short>(),
md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
: gmx::ArrayRef<const unsigned short>());
+ fr->longRangeNonbondeds->updateAfterPartition(*md);
}
std::unique_ptr<UpdateConstrainGpu> integrator;
: gmx::ArrayRef<const unsigned short>(),
md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
: gmx::ArrayRef<const unsigned short>());
+ fr->longRangeNonbondeds->updateAfterPartition(*md);
}
}
&f.view(),
force_vir,
*md,
+ fr->longRangeNonbondeds.get(),
nrnb,
wcycle,
shellfc,
mu_tot,
t,
ed ? ed->getLegacyED() : nullptr,
+ fr->longRangeNonbondeds.get(),
(bNS ? GMX_FORCE_NS : 0) | force_flags,
ddBalanceRegionHandler);
}
: gmx::ArrayRef<const unsigned short>(),
md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
: gmx::ArrayRef<const unsigned short>());
+ fr->longRangeNonbondeds->updateAfterPartition(*md);
}
bFirstStep = FALSE;
// (Global topology should persist.)
update_mdatoms(mdatoms, state->lambda[FreeEnergyPerturbationCouplingType::Mass]);
+ fr->longRangeNonbondeds->updateAfterPartition(*mdatoms);
if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl != 0)
{
update_mdatoms(mdatoms, state->lambda[FreeEnergyPerturbationCouplingType::Mass]);
}
+ fr->longRangeNonbondeds->updateAfterPartition(*mdatoms);
+
force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES
| GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 is solved
GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0));
&f.view(),
force_vir,
*mdatoms,
+ fr->longRangeNonbondeds.get(),
nrnb,
wcycle,
shellfc,
mu_tot,
t,
ed,
+ fr->longRangeNonbondeds.get(),
GMX_FORCE_NS | force_flags,
ddBalanceRegionHandler);
}
setCoordinates(&pairSearchCoordinates, localCoordinates);
}
+ fr->longRangeNonbondeds->updateAfterPartition(*mdAtoms->mdatoms());
+
/* Calc force & energy on new trial position */
/* 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
mu_tot,
t,
nullptr,
+ fr->longRangeNonbondeds.get(),
GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES | GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY
| (bNS ? GMX_FORCE_NS : 0),
DDBalanceRegionHandler(cr));
em_state_t state_work{};
+ fr->longRangeNonbondeds->updateAfterPartition(*mdAtoms->mdatoms());
ObservablesReducer observablesReducer = observablesReducerBuilder->build();
/* Init em and store the local state in state_minimum */
&state_work.f.view(),
vir,
*mdatoms,
+ fr->longRangeNonbondeds.get(),
nrnb,
wcycle,
shellfc,
}
auto* mdatoms = mdAtoms->mdatoms();
+ fr->longRangeNonbondeds->updateAfterPartition(*mdatoms);
// NOTE: The global state is no longer used at this point.
// But state_global is still used as temporary storage space for writing
update_mdatoms(mdatoms, state->lambda[FreeEnergyPerturbationCouplingType::Mass]);
}
+ fr->longRangeNonbondeds->updateAfterPartition(*mdatoms);
+
force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES
| GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 and #3400 are solved
GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0));
&f.view(),
force_vir,
*mdatoms,
+ fr->longRangeNonbondeds.get(),
nrnb,
wcycle,
shellfc,
mu_tot,
t,
ed,
+ fr->longRangeNonbondeds.get(),
GMX_FORCE_NS | force_flags,
ddBalanceRegionHandler);
}
deviceStreamManager->bondedStream(havePPDomainDecomposition(cr)),
wcycle.get());
}
+ fr->longRangeNonbondeds = std::make_unique<CpuPpLongRangeNonbondeds>(fr->n_tpi,
+ fr->ic->ewaldcoeff_q,
+ fr->ic->epsilon_r,
+ fr->qsum,
+ fr->ic->eeltype,
+ fr->ic->vdwtype,
+ *inputrec,
+ &nrnb,
+ wcycle.get(),
+ fplog);
/* Initialize the mdAtoms structure.
* mdAtoms is not filled with atom data,
gmx::ForceBuffersView* f,
tensor force_vir,
const t_mdatoms& md,
+ CpuPpLongRangeNonbondeds* longRangeNonbondeds,
t_nrnb* nrnb,
gmx_wallcycle* wcycle,
gmx_shellfc_t* shfc,
mu_tot,
t,
nullptr,
+ longRangeNonbondeds,
(bDoNS ? GMX_FORCE_NS : 0) | shellfc_flags,
ddBalanceRegionHandler);
mu_tot,
t,
nullptr,
+ longRangeNonbondeds,
shellfc_flags,
ddBalanceRegionHandler);
accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals.get());
struct t_mdatoms;
struct t_nrnb;
class t_state;
+class CpuPpLongRangeNonbondeds;
namespace gmx
{
gmx::ForceBuffersView* f,
tensor force_vir,
const t_mdatoms& md,
+ CpuPpLongRangeNonbondeds* longRangeNonbondeds,
t_nrnb* nrnb,
gmx_wallcycle* wcycle,
gmx_shellfc_t* shfc,
/* Note: NonLocal refers to the inserted molecule */
fr->nbv->convertCoordinates(AtomLocality::NonLocal, x);
+ fr->longRangeNonbondeds->updateAfterPartition(*mdatoms);
/* Clear some matrix variables */
clear_mat(force_vir);
mu_tot,
t,
nullptr,
+ fr->longRangeNonbondeds.get(),
GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY | (bStateChanged ? GMX_FORCE_STATECHANGED : 0),
DDBalanceRegionHandler(nullptr));
std::feclearexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
class DeviceContext;
class DispersionCorrection;
class ListedForces;
+class CpuPpLongRangeNonbondeds;
struct t_fcdata;
struct t_forcetable;
struct interaction_const_t;
/* Forward declaration of type for managing Ewald tables */
struct gmx_ewald_tab_t;
-struct ewald_corr_thread_t;
-
/*! \brief Helper force buffers for ForceOutputs
*
* This class stores intermediate force buffers that are used
gmx_pme_t* pmedata = nullptr;
LongRangeVdW ljpme_combination_rule = LongRangeVdW::Geom;
- /* PME/Ewald stuff */
- std::unique_ptr<gmx_ewald_tab_t> ewald_table;
-
/* Non bonded Parameter lists */
int ntype = 0; /* Number of atom types */
bool haveBuckingham = false;
// The listed forces calculation data for GPU
std::unique_ptr<gmx::ListedForcesGpu> listedForcesGpu;
- /* Ewald correction thread local virial and energy data */
- int nthread_ewc = 0;
- std::vector<ewald_corr_thread_t> ewc_t;
+ // The long range non-bonded forces
+ std::unique_ptr<CpuPpLongRangeNonbondeds> longRangeNonbondeds;
gmx::ForceProviders* forceProviders = nullptr;
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdrun/shellfc.h"
#include "gromacs/mdtypes/forcebuffers.h"
+#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/mdrunoptions.h"
#include "gromacs/mdtypes/simulation_workload.h"
const MDAtoms* mdAtoms,
t_nrnb* nrnb,
t_forcerec* fr,
-
- gmx_wallcycle* wcycle,
- MdrunScheduleWorkload* runScheduleWork,
- VirtualSitesHandler* vsite,
- ImdSession* imdSession,
- pull_t* pull_work,
- Constraints* constr,
- const gmx_mtop_t& globalTopology,
- gmx_enfrot* enforcedRotation) :
+ gmx_wallcycle* wcycle,
+ MdrunScheduleWorkload* runScheduleWork,
+ VirtualSitesHandler* vsite,
+ ImdSession* imdSession,
+ pull_t* pull_work,
+ Constraints* constr,
+ const gmx_mtop_t& globalTopology,
+ gmx_enfrot* enforcedRotation) :
shellfc_(init_shell_flexcon(fplog,
globalTopology,
constr ? constr->numFlexibleConstraints() : 0,
isVerbose_(isVerbose),
nShellRelaxationSteps_(0),
ddBalanceRegionHandler_(cr),
+ longRangeNonbondeds_(std::make_unique<CpuPpLongRangeNonbondeds>(fr->n_tpi,
+ fr->ic->ewaldcoeff_q,
+ fr->ic->epsilon_r,
+ fr->qsum,
+ fr->ic->eeltype,
+ fr->ic->vdwtype,
+ *inputrec,
+ nrnb,
+ wcycle,
+ fplog)),
lambda_(),
fplog_(fplog),
cr_(cr),
ArrayRef<real> lambda =
freeEnergyPerturbationData_ ? freeEnergyPerturbationData_->lambdaView() : lambda_;
+ longRangeNonbondeds_->updateAfterPartition(*mdAtoms_->mdatoms());
+
if (doShellFC)
{
auto v = statePropagatorData_->velocitiesView();
&forces,
force_vir,
*mdAtoms_->mdatoms(),
+ longRangeNonbondeds_.get(),
nrnb_,
wcycle_,
shellfc_,
energyData_->muTot(),
time,
ed,
+ longRangeNonbondeds_.get(),
static_cast<int>(flags),
ddBalanceRegionHandler_);
}
return std::nullopt;
}
+DomDecCallback ForceElement::registerDomDecCallback()
+{
+ return [this]() { longRangeNonbondeds_->updateAfterPartition(*mdAtoms_->mdatoms()); };
+}
+
ISimulatorElement*
ForceElement::getElementPointerImpl(LegacySimulatorData* legacySimulatorData,
ModularSimulatorAlgorithmBuilderHelper* builderHelper,
struct gmx_enfrot;
struct gmx_shellfc_t;
struct gmx_wallcycle;
+class CpuPpLongRangeNonbondeds;
struct pull_t;
struct t_nrnb;
public ISimulatorElement,
public ITopologyHolderClient,
public INeighborSearchSignallerClient,
- public IEnergySignallerClient
+ public IEnergySignallerClient,
+ public IDomDecHelperClient
{
public:
//! Constructor
GlobalCommunicationHelper* globalCommunicationHelper,
ObservablesReducer* observablesReducer);
+ //! Callback on domain decomposition repartitioning
+ DomDecCallback registerDomDecCallback() override;
+
private:
//! ITopologyHolderClient implementation
void setTopology(const gmx_localtop_t* top) override;
//! DD / DLB helper object
const DDBalanceRegionHandler ddBalanceRegionHandler_;
+ //! Long range force calculator
+ std::unique_ptr<CpuPpLongRangeNonbondeds> longRangeNonbondeds_;
/* \brief The FEP lambda vector
*
simulationsShareState);
registerWithInfrastructureAndSignallers(trajectoryElement.get());
- // Build domdec helper
+ // Build domdec helper (free energy element is a client, so keep this after it is built)
if (DOMAINDECOMP(legacySimulatorData_->cr))
{
algorithm.domDecHelper_ =
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff