From: Berk Hess Date: Thu, 27 Aug 2020 06:09:17 +0000 (+0000) Subject: Add ForceBuffers class X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=8601caf88e4fb4047c987830c14572ede4c994c6;p=alexxy%2Fgromacs.git Add ForceBuffers class This change replaces the type of the force buffer in mdrun from PaddedHostVector to a new class ForceBuffers. This new class currently only holds one force buffer, but this is in preparation for multiple time stepping where a second force buffer is needed. --- diff --git a/src/gromacs/domdec/mdsetup.cpp b/src/gromacs/domdec/mdsetup.cpp index 9c71dd410c..6a8559217f 100644 --- a/src/gromacs/domdec/mdsetup.cpp +++ b/src/gromacs/domdec/mdsetup.cpp @@ -44,6 +44,7 @@ #include "gromacs/mdlib/mdatoms.h" #include "gromacs/mdlib/vsite.h" #include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/interaction_const.h" @@ -62,16 +63,16 @@ namespace gmx * The final solution should be an MD algorithm base class with methods * for initialization and atom-data setup. */ -void mdAlgorithmsSetupAtomData(const t_commrec* cr, - const t_inputrec* ir, - const gmx_mtop_t& top_global, - gmx_localtop_t* top, - t_forcerec* fr, - PaddedHostVector* force, - MDAtoms* mdAtoms, - Constraints* constr, - VirtualSitesHandler* vsite, - gmx_shellfc_t* shellfc) +void mdAlgorithmsSetupAtomData(const t_commrec* cr, + const t_inputrec* ir, + const gmx_mtop_t& top_global, + gmx_localtop_t* top, + t_forcerec* fr, + ForceBuffers* force, + MDAtoms* mdAtoms, + Constraints* constr, + VirtualSitesHandler* vsite, + gmx_shellfc_t* shellfc) { bool usingDomDec = DOMAINDECOMP(cr); @@ -94,10 +95,7 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, if (force != nullptr) { - /* We need to allocate one element extra, since we might use - * (unaligned) 4-wide SIMD loads to access rvec entries. - */ - force->resizeWithPadding(numTotalAtoms); + force->resize(numTotalAtoms); } atoms2md(&top_global, ir, numAtomIndex, diff --git a/src/gromacs/domdec/mdsetup.h b/src/gromacs/domdec/mdsetup.h index d406b09f7f..ab0da8fcf1 100644 --- a/src/gromacs/domdec/mdsetup.h +++ b/src/gromacs/domdec/mdsetup.h @@ -43,7 +43,6 @@ #ifndef GMX_DOMDEC_MDSETUP_H #define GMX_DOMDEC_MDSETUP_H -#include "gromacs/gpu_utils/hostallocator.h" #include "gromacs/math/vectypes.h" struct bonded_threading_t; @@ -58,6 +57,7 @@ struct t_mdatoms; namespace gmx { class Constraints; +class ForceBuffers; class MDAtoms; class VirtualSitesHandler; @@ -87,16 +87,16 @@ void make_local_shells(const t_commrec* cr, const t_mdatoms* md, gmx_shellfc_t* * \param[in,out] vsite The virtual site data, can be NULL * \param[in,out] shellfc The shell/flexible-constraint data, can be NULL */ -void mdAlgorithmsSetupAtomData(const t_commrec* cr, - const t_inputrec* ir, - const gmx_mtop_t& top_global, - gmx_localtop_t* top, - t_forcerec* fr, - PaddedHostVector* force, - MDAtoms* mdAtoms, - Constraints* constr, - VirtualSitesHandler* vsite, - gmx_shellfc_t* shellfc); +void mdAlgorithmsSetupAtomData(const t_commrec* cr, + const t_inputrec* ir, + const gmx_mtop_t& top_global, + gmx_localtop_t* top, + t_forcerec* fr, + ForceBuffers* force, + MDAtoms* mdAtoms, + Constraints* constr, + VirtualSitesHandler* vsite, + gmx_shellfc_t* shellfc); } // namespace gmx diff --git a/src/gromacs/domdec/partition.cpp b/src/gromacs/domdec/partition.cpp index 2d5d6cb0b8..c7a01ae473 100644 --- a/src/gromacs/domdec/partition.cpp +++ b/src/gromacs/domdec/partition.cpp @@ -2624,27 +2624,27 @@ void print_dd_statistics(const t_commrec* cr, const t_inputrec* ir, FILE* fplog) } //!\brief TODO Remove fplog when group scheme and charge groups are gone -void dd_partition_system(FILE* fplog, - const gmx::MDLogger& mdlog, - int64_t step, - const t_commrec* cr, - gmx_bool bMasterState, - int nstglobalcomm, - t_state* state_global, - const gmx_mtop_t& top_global, - const t_inputrec* ir, - gmx::ImdSession* imdSession, - pull_t* pull_work, - t_state* state_local, - PaddedHostVector* f, - gmx::MDAtoms* mdAtoms, - gmx_localtop_t* top_local, - t_forcerec* fr, - gmx::VirtualSitesHandler* vsite, - gmx::Constraints* constr, - t_nrnb* nrnb, - gmx_wallcycle* wcycle, - gmx_bool bVerbose) +void dd_partition_system(FILE* fplog, + const gmx::MDLogger& mdlog, + int64_t step, + const t_commrec* cr, + gmx_bool bMasterState, + int nstglobalcomm, + t_state* state_global, + const gmx_mtop_t& top_global, + const t_inputrec* ir, + gmx::ImdSession* imdSession, + pull_t* pull_work, + t_state* state_local, + gmx::ForceBuffers* f, + gmx::MDAtoms* mdAtoms, + gmx_localtop_t* top_local, + t_forcerec* fr, + gmx::VirtualSitesHandler* vsite, + gmx::Constraints* constr, + t_nrnb* nrnb, + gmx_wallcycle* wcycle, + gmx_bool bVerbose) { gmx_domdec_t* dd; gmx_domdec_comm_t* comm; diff --git a/src/gromacs/domdec/partition.h b/src/gromacs/domdec/partition.h index e06138a1b7..de05c4a738 100644 --- a/src/gromacs/domdec/partition.h +++ b/src/gromacs/domdec/partition.h @@ -45,7 +45,8 @@ #ifndef GMX_DOMDEC_PARTITION_H #define GMX_DOMDEC_PARTITION_H -#include "gromacs/gpu_utils/hostallocator.h" +#include + #include "gromacs/math/vectypes.h" #include "gromacs/utility/basedefinitions.h" #include "gromacs/utility/real.h" @@ -64,7 +65,10 @@ class t_state; namespace gmx { +template +class ArrayRef; class Constraints; +class ForceBuffers; class ImdSession; class MDAtoms; class MDLogger; @@ -97,7 +101,7 @@ void print_dd_statistics(const t_commrec* cr, const t_inputrec* ir, FILE* fplog) * \param[in] pull_work Pulling data * \param[in] state_local Local state * \param[in] f Force buffer - * \param[in] mdatoms MD atoms + * \param[in] mdAtoms MD atoms * \param[in] top_local Local topology * \param[in] fr Force record * \param[in] vsite Virtual sites handler @@ -106,27 +110,27 @@ void print_dd_statistics(const t_commrec* cr, const t_inputrec* ir, FILE* fplog) * \param[in] wcycle Timers * \param[in] bVerbose Be verbose */ -void dd_partition_system(FILE* fplog, - const gmx::MDLogger& mdlog, - int64_t step, - const t_commrec* cr, - gmx_bool bMasterState, - int nstglobalcomm, - t_state* state_global, - const gmx_mtop_t& top_global, - const t_inputrec* ir, - gmx::ImdSession* imdSession, - pull_t* pull_work, - t_state* state_local, - gmx::PaddedHostVector* f, - gmx::MDAtoms* mdatoms, - gmx_localtop_t* top_local, - t_forcerec* fr, - gmx::VirtualSitesHandler* vsite, - gmx::Constraints* constr, - t_nrnb* nrnb, - gmx_wallcycle* wcycle, - gmx_bool bVerbose); +void dd_partition_system(FILE* fplog, + const gmx::MDLogger& mdlog, + int64_t step, + const t_commrec* cr, + gmx_bool bMasterState, + int nstglobalcomm, + t_state* state_global, + const gmx_mtop_t& top_global, + const t_inputrec* ir, + gmx::ImdSession* imdSession, + pull_t* pull_work, + t_state* state_local, + gmx::ForceBuffers* f, + gmx::MDAtoms* mdAtoms, + gmx_localtop_t* top_local, + t_forcerec* fr, + gmx::VirtualSitesHandler* vsite, + gmx::Constraints* constr, + t_nrnb* nrnb, + gmx_wallcycle* wcycle, + gmx_bool bVerbose); /*! \brief Check whether bonded interactions are missing, if appropriate * diff --git a/src/gromacs/mdlib/force.h b/src/gromacs/mdlib/force.h index 3db3f96db4..711ce856b8 100644 --- a/src/gromacs/mdlib/force.h +++ b/src/gromacs/mdlib/force.h @@ -63,6 +63,7 @@ struct t_nrnb; namespace gmx { class Awh; +class ForceBuffersView; class ForceOutputs; class ForceWithVirial; class ImdSession; @@ -87,7 +88,7 @@ void do_force(FILE* log, const matrix box, gmx::ArrayRefWithPadding coordinates, history_t* hist, - gmx::ArrayRefWithPadding force, + gmx::ForceBuffersView* force, tensor vir_force, const t_mdatoms* mdatoms, gmx_enerdata_t* enerd, diff --git a/src/gromacs/mdlib/mdoutf.cpp b/src/gromacs/mdlib/mdoutf.cpp index 795aae408b..5b94dbd1d4 100644 --- a/src/gromacs/mdlib/mdoutf.cpp +++ b/src/gromacs/mdlib/mdoutf.cpp @@ -476,19 +476,19 @@ static void write_checkpoint(const char* fn, #endif /* end GMX_FAHCORE block */ } -void mdoutf_write_to_trajectory_files(FILE* fplog, - const t_commrec* cr, - gmx_mdoutf_t of, - int mdof_flags, - int natoms, - int64_t step, - double t, - t_state* state_local, - t_state* state_global, - ObservablesHistory* observablesHistory, - gmx::ArrayRef f_local) +void mdoutf_write_to_trajectory_files(FILE* fplog, + const t_commrec* cr, + gmx_mdoutf_t of, + int mdof_flags, + int natoms, + int64_t step, + double t, + t_state* state_local, + t_state* state_global, + ObservablesHistory* observablesHistory, + gmx::ArrayRef f_local) { - rvec* f_global; + const rvec* f_global; if (DOMAINDECOMP(cr)) { @@ -512,8 +512,9 @@ void mdoutf_write_to_trajectory_files(FILE* fplog, f_global = of->f_global; if (mdof_flags & MDOF_F) { - dd_collect_vec(cr->dd, state_local, f_local, - gmx::arrayRefFromArray(reinterpret_cast(f_global), f_local.size())); + dd_collect_vec( + cr->dd, state_local, f_local, + gmx::arrayRefFromArray(reinterpret_cast(of->f_global), f_local.size())); } } else diff --git a/src/gromacs/mdlib/mdoutf.h b/src/gromacs/mdlib/mdoutf.h index 4e704efecf..fe534fab60 100644 --- a/src/gromacs/mdlib/mdoutf.h +++ b/src/gromacs/mdlib/mdoutf.h @@ -121,17 +121,17 @@ void done_mdoutf(gmx_mdoutf_t of); * \param[in] observablesHistory Pointer to the ObservableHistory object. * \param[in] f_local The local forces. */ -void mdoutf_write_to_trajectory_files(FILE* fplog, - const t_commrec* cr, - gmx_mdoutf_t of, - int mdof_flags, - int natoms, - int64_t step, - double t, - t_state* state_local, - t_state* state_global, - ObservablesHistory* observablesHistory, - gmx::ArrayRef f_local); +void mdoutf_write_to_trajectory_files(FILE* fplog, + const t_commrec* cr, + gmx_mdoutf_t of, + int mdof_flags, + int natoms, + int64_t step, + double t, + t_state* state_local, + t_state* state_global, + ObservablesHistory* observablesHistory, + gmx::ArrayRef f_local); /*! \brief Get the output interval of box size of uncompressed TNG output. * Returns 0 if no uncompressed TNG file is open. diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 3f2cd68dd8..dda01b711f 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -84,6 +84,7 @@ #include "gromacs/mdlib/wholemoleculetransform.h" #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forceoutput.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/iforceprovider.h" @@ -1015,7 +1016,7 @@ void do_force(FILE* fplog, const matrix box, gmx::ArrayRefWithPadding x, history_t* hist, - gmx::ArrayRefWithPadding force, + gmx::ForceBuffersView* forceView, tensor vir_force, const t_mdatoms* mdatoms, gmx_enerdata_t* enerd, @@ -1029,6 +1030,7 @@ void do_force(FILE* fplog, int legacyFlags, const DDBalanceRegionHandler& ddBalanceRegionHandler) { + auto force = forceView->forceWithPadding(); GMX_ASSERT(force.unpaddedArrayRef().ssize() >= fr->natoms_force_constr, "The size of the force buffer should be at least the number of atoms to compute " "forces for"); diff --git a/src/gromacs/mdlib/trajectory_writing.cpp b/src/gromacs/mdlib/trajectory_writing.cpp index 9c87f46c4c..6d250ceeb2 100644 --- a/src/gromacs/mdlib/trajectory_writing.cpp +++ b/src/gromacs/mdlib/trajectory_writing.cpp @@ -54,28 +54,28 @@ #include "gromacs/topology/topology.h" #include "gromacs/utility/smalloc.h" -void do_md_trajectory_writing(FILE* fplog, - t_commrec* cr, - int nfile, - const t_filenm fnm[], - int64_t step, - int64_t step_rel, - double t, - t_inputrec* ir, - t_state* state, - t_state* state_global, - ObservablesHistory* observablesHistory, - const gmx_mtop_t* top_global, - t_forcerec* fr, - gmx_mdoutf_t outf, - const gmx::EnergyOutput& energyOutput, - gmx_ekindata_t* ekind, - gmx::ArrayRef f, - gmx_bool bCPT, - gmx_bool bRerunMD, - gmx_bool bLastStep, - gmx_bool bDoConfOut, - gmx_bool bSumEkinhOld) +void do_md_trajectory_writing(FILE* fplog, + t_commrec* cr, + int nfile, + const t_filenm fnm[], + int64_t step, + int64_t step_rel, + double t, + t_inputrec* ir, + t_state* state, + t_state* state_global, + ObservablesHistory* observablesHistory, + const gmx_mtop_t* top_global, + t_forcerec* fr, + gmx_mdoutf_t outf, + const gmx::EnergyOutput& energyOutput, + gmx_ekindata_t* ekind, + gmx::ArrayRef f, + gmx_bool bCPT, + gmx_bool bRerunMD, + gmx_bool bLastStep, + gmx_bool bDoConfOut, + gmx_bool bSumEkinhOld) { int mdof_flags; rvec* x_for_confout = nullptr; diff --git a/src/gromacs/mdlib/trajectory_writing.h b/src/gromacs/mdlib/trajectory_writing.h index bd17e22d26..a74c707e77 100644 --- a/src/gromacs/mdlib/trajectory_writing.h +++ b/src/gromacs/mdlib/trajectory_writing.h @@ -60,28 +60,28 @@ class EnergyOutput; * * This routine does communication (e.g. collecting distributed coordinates) */ -void do_md_trajectory_writing(FILE* fplog, - struct t_commrec* cr, - int nfile, - const t_filenm fnm[], - int64_t step, - int64_t step_rel, - double t, - t_inputrec* ir, - t_state* state, - t_state* state_global, - ObservablesHistory* observablesHistory, - const gmx_mtop_t* top_global, - t_forcerec* fr, - gmx_mdoutf_t outf, - const gmx::EnergyOutput& energyOutput, - gmx_ekindata_t* ekind, - gmx::ArrayRef f, - gmx_bool bCPT, - gmx_bool bRerunMD, - gmx_bool bLastStep, - gmx_bool bDoConfOut, - gmx_bool bSumEkinhOld); +void do_md_trajectory_writing(FILE* fplog, + struct t_commrec* cr, + int nfile, + const t_filenm fnm[], + int64_t step, + int64_t step_rel, + double t, + t_inputrec* ir, + t_state* state, + t_state* state_global, + ObservablesHistory* observablesHistory, + const gmx_mtop_t* top_global, + t_forcerec* fr, + gmx_mdoutf_t outf, + const gmx::EnergyOutput& energyOutput, + gmx_ekindata_t* ekind, + gmx::ArrayRef f, + gmx_bool bCPT, + gmx_bool bRerunMD, + gmx_bool bLastStep, + gmx_bool bDoConfOut, + gmx_bool bSumEkinhOld); #endif diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 6e683a8dd8..ac08b56afc 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -110,6 +110,7 @@ #include "gromacs/mdtypes/df_history.h" #include "gromacs/mdtypes/energyhistory.h" #include "gromacs/mdtypes/fcdata.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" @@ -174,16 +175,15 @@ void gmx::LegacySimulator::do_md() int i, m; rvec mu_tot; matrix pressureCouplingMu, M; - gmx_repl_ex_t repl_ex = nullptr; - PaddedHostVector f{}; - gmx_global_stat_t gstat; - gmx_shellfc_t* shellfc; - gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition; - gmx_bool bTemp, bPres, bTrotter; - real dvdl_constr; - std::vector cbuf; - matrix lastbox; - int lamnew = 0; + gmx_repl_ex_t repl_ex = nullptr; + gmx_global_stat_t gstat; + gmx_shellfc_t* shellfc; + gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition; + gmx_bool bTemp, bPres, bTrotter; + real dvdl_constr; + std::vector cbuf; + matrix lastbox; + int lamnew = 0; /* for FEP */ int nstfep = 0; double cycles; @@ -322,8 +322,15 @@ void gmx::LegacySimulator::do_md() auto mdatoms = mdAtoms->mdatoms(); - std::unique_ptr integrator; + const auto& simulationWork = runScheduleWork->simulationWork; + const bool useGpuForPme = simulationWork.useGpuPme; + const bool useGpuForNonbonded = simulationWork.useGpuNonbonded; + const bool useGpuForBufferOps = simulationWork.useGpuBufferOps; + const bool useGpuForUpdate = simulationWork.useGpuUpdate; + ForceBuffers f(((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate) + ? PinningPolicy::PinnedIfSupported + : PinningPolicy::CannotBePinned); if (DOMAINDECOMP(cr)) { stateInstance = std::make_unique(); @@ -349,11 +356,7 @@ void gmx::LegacySimulator::do_md() upd.setNumAtoms(state->natoms); } - const auto& simulationWork = runScheduleWork->simulationWork; - const bool useGpuForPme = simulationWork.useGpuPme; - const bool useGpuForNonbonded = simulationWork.useGpuNonbonded; - const bool useGpuForBufferOps = simulationWork.useGpuBufferOps; - const bool useGpuForUpdate = simulationWork.useGpuUpdate; + std::unique_ptr integrator; StatePropagatorDataGpu* stateGpu = fr->stateGpu; @@ -421,10 +424,6 @@ void gmx::LegacySimulator::do_md() { changePinningPolicy(&state->x, PinningPolicy::PinnedIfSupported); } - if ((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate) - { - changePinningPolicy(&f, PinningPolicy::PinnedIfSupported); - } if (useGpuForUpdate) { changePinningPolicy(&state->v, PinningPolicy::PinnedIfSupported); @@ -925,8 +924,8 @@ void gmx::LegacySimulator::do_md() relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir, 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, + state->v.arrayRefWithPadding(), state->box, state->lambda, + &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc, fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); } else @@ -951,8 +950,8 @@ void gmx::LegacySimulator::do_md() */ do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, state->lambda, fr, - runScheduleWork, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr, + &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, + vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr, (bNS ? GMX_FORCE_NS : 0) | force_flags, ddBalanceRegionHandler); } @@ -985,8 +984,8 @@ void gmx::LegacySimulator::do_md() trotter_seq, ettTSEQ1); } - upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcdata, ekind, M, - etrtVELOCITY1, cr, constr != nullptr); + upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind, + M, etrtVELOCITY1, cr, constr != nullptr); wallcycle_stop(wcycle, ewcUPDATE); constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir); @@ -1155,7 +1154,7 @@ void gmx::LegacySimulator::do_md() if (runScheduleWork->stepWork.useGpuFBufferOps && (simulationWork.useGpuUpdate && !vsite) && do_per_step(step, ir->nstfout)) { - stateGpu->copyForcesFromGpu(ArrayRef(f), AtomLocality::Local); + stateGpu->copyForcesFromGpu(f.view().force(), AtomLocality::Local); stateGpu->waitForcesReadyOnHost(AtomLocality::Local); } /* Now we have the energies and forces corresponding to the @@ -1163,9 +1162,9 @@ void gmx::LegacySimulator::do_md() * the update. */ do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, state_global, - observablesHistory, top_global, fr, outf, energyOutput, ekind, f, - checkpointHandler->isCheckpointingStep(), bRerunMD, bLastStep, - mdrunOptions.writeConfout, bSumEkinhOld); + observablesHistory, top_global, fr, outf, energyOutput, ekind, + f.view().force(), checkpointHandler->isCheckpointingStep(), + bRerunMD, bLastStep, mdrunOptions.writeConfout, bSumEkinhOld); /* Check if IMD step and do IMD communication, if bIMD is TRUE. */ bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x.rvec_array(), t); @@ -1233,8 +1232,8 @@ void gmx::LegacySimulator::do_md() if (EI_VV(ir->eI)) { /* velocity half-step update */ - upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcdata, ekind, M, - etrtVELOCITY2, cr, constr != nullptr); + upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind, + M, etrtVELOCITY2, cr, constr != nullptr); } /* Above, initialize just copies ekinh into ekin, @@ -1276,7 +1275,7 @@ void gmx::LegacySimulator::do_md() // If the buffer ops were not offloaded this step, the forces are on the host and have to be copied if (!runScheduleWork->stepWork.useGpuFBufferOps) { - stateGpu->copyForcesToGpu(ArrayRef(f), AtomLocality::Local); + stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::Local); } const bool doTemperatureScaling = @@ -1299,8 +1298,8 @@ void gmx::LegacySimulator::do_md() } else { - upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcdata, ekind, M, - etrtPOSITION, cr, constr != nullptr); + upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind, + M, etrtPOSITION, cr, constr != nullptr); wallcycle_stop(wcycle, ewcUPDATE); @@ -1330,8 +1329,8 @@ void gmx::LegacySimulator::do_md() /* now we know the scaling, we can compute the positions again */ std::copy(cbuf.begin(), cbuf.end(), state->x.begin()); - upd.update_coords(*ir, step, mdatoms, state, f.arrayRefWithPadding(), fcdata, ekind, M, - etrtPOSITION, cr, constr != nullptr); + upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind, + M, etrtPOSITION, cr, constr != nullptr); wallcycle_stop(wcycle, ewcUPDATE); /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */ diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index 84bea6856b..29d846b851 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -104,6 +104,7 @@ #include "gromacs/mdtypes/df_history.h" #include "gromacs/mdtypes/enerdata.h" #include "gromacs/mdtypes/energyhistory.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" @@ -139,17 +140,17 @@ using gmx::SimulationSignaller; void gmx::LegacySimulator::do_mimic() { - t_inputrec* ir = inputrec; - int64_t step, step_rel; - double t, lam0[efptNR]; - bool isLastStep = false; - bool doFreeEnergyPerturbation = false; - unsigned int force_flags; - tensor force_vir, shake_vir, total_vir, pres; - rvec mu_tot; - PaddedHostVector f{}; - gmx_global_stat_t gstat; - gmx_shellfc_t* shellfc; + t_inputrec* ir = inputrec; + int64_t step, step_rel; + double t, lam0[efptNR]; + bool isLastStep = false; + bool doFreeEnergyPerturbation = false; + unsigned int force_flags; + tensor force_vir, shake_vir, total_vir, pres; + rvec mu_tot; + ForceBuffers f; + gmx_global_stat_t gstat; + gmx_shellfc_t* shellfc; double cycles; @@ -415,8 +416,8 @@ void gmx::LegacySimulator::do_mimic() relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir, 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, + state->v.arrayRefWithPadding(), state->box, state->lambda, + &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc, fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); } else @@ -430,7 +431,7 @@ void gmx::LegacySimulator::do_mimic() gmx_edsam* ed = nullptr; do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, + &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler); } @@ -443,8 +444,8 @@ void gmx::LegacySimulator::do_mimic() const bool bSumEkinhOld = false; do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, state_global, observablesHistory, top_global, fr, outf, - energyOutput, ekind, f, isCheckpointingStep, doRerun, - isLastStep, mdrunOptions.writeConfout, bSumEkinhOld); + energyOutput, ekind, f.view().force(), isCheckpointingStep, + doRerun, isLastStep, mdrunOptions.writeConfout, bSumEkinhOld); } stopHandler->setSignal(); @@ -478,7 +479,7 @@ void gmx::LegacySimulator::do_mimic() { gmx::HostVector fglobal(top_global->natoms); gmx::ArrayRef ftemp; - gmx::ArrayRef flocal = gmx::makeArrayRef(f); + gmx::ArrayRef flocal = f.view().force(); if (DOMAINDECOMP(cr)) { ftemp = gmx::makeArrayRef(fglobal); @@ -486,7 +487,7 @@ void gmx::LegacySimulator::do_mimic() } else { - ftemp = gmx::makeArrayRef(f); + ftemp = f.view().force(); } if (MASTER(cr)) diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index c96ef39acd..7e6b030e3c 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -91,6 +91,7 @@ #include "gromacs/mdrunutility/handlerestart.h" #include "gromacs/mdrunutility/printtime.h" #include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/interaction_const.h" @@ -123,7 +124,7 @@ typedef struct em_state //! Copy of the global state t_state s; //! Force array - PaddedHostVector f; + gmx::ForceBuffers f; //! Potential energy real epot; //! Norm of the force @@ -255,13 +256,13 @@ static void print_converged(FILE* fp, } //! Compute the norm and max of the force array in parallel -static void get_f_norm_max(const t_commrec* cr, - t_grpopts* opts, - t_mdatoms* mdatoms, - const rvec* f, - real* fnorm, - real* fmax, - int* a_fmax) +static void get_f_norm_max(const t_commrec* cr, + t_grpopts* opts, + t_mdatoms* mdatoms, + gmx::ArrayRef f, + real* fnorm, + real* fmax, + int* a_fmax) { double fnorm2, *sum; real fmax2, fam; @@ -355,7 +356,7 @@ static void get_f_norm_max(const t_commrec* cr, //! Compute the norm of the force static void get_state_f_norm_max(const t_commrec* cr, t_grpopts* opts, t_mdatoms* mdatoms, em_state_t* ems) { - get_f_norm_max(cr, opts, mdatoms, ems->f.rvec_array(), &ems->fnorm, &ems->fmax, &ems->a_fmax); + get_f_norm_max(cr, opts, mdatoms, ems->f.view().force(), &ems->fnorm, &ems->fmax, &ems->a_fmax); } //! Initialize the energy minimization @@ -533,7 +534,7 @@ static void write_em_traj(FILE* fplog, mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags, top_global->natoms, step, static_cast(step), &state->s, state_global, - observablesHistory, state->f); + observablesHistory, state->f.view().force()); if (confout != nullptr) { @@ -569,15 +570,15 @@ static void write_em_traj(FILE* fplog, //! \brief Do one minimization step // // \returns true when the step succeeded, false when a constraint error occurred -static bool do_em_step(const t_commrec* cr, - t_inputrec* ir, - t_mdatoms* md, - em_state_t* ems1, - real a, - const PaddedHostVector* force, - em_state_t* ems2, - gmx::Constraints* constr, - int64_t count) +static bool do_em_step(const t_commrec* cr, + t_inputrec* ir, + t_mdatoms* md, + em_state_t* ems1, + real a, + gmx::ArrayRefWithPadding force, + em_state_t* ems2, + gmx::Constraints* constr, + int64_t count) { t_state *s1, *s2; @@ -600,7 +601,7 @@ static bool do_em_step(const t_commrec* cr, if (s2->natoms != s1->natoms) { state_change_natoms(s2, s1->natoms); - ems2->f.resizeWithPadding(s2->natoms); + ems2->f.resize(s2->natoms); } if (DOMAINDECOMP(cr) && s2->cg_gl.size() != s1->cg_gl.size()) { @@ -620,7 +621,7 @@ static bool do_em_step(const t_commrec* cr, { const rvec* x1 = s1->x.rvec_array(); rvec* x2 = s2->x.rvec_array(); - const rvec* f = force->rvec_array(); + const rvec* f = as_rvec_array(force.unpaddedArrayRef().data()); int gf = 0; #pragma omp for schedule(static) nowait @@ -840,9 +841,8 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres, * We do not unshift, so molecules are always whole in congrad.c */ 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, ems->s.lambda, fr, - runScheduleWork, vsite, mu_tot, t, nullptr, + top, ems->s.box, ems->s.x.arrayRefWithPadding(), &ems->s.hist, &ems->f.view(), 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)); @@ -887,7 +887,7 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres, bool computeEnergy = false; bool computeVirial = true; dvdl_constr = 0; - auto f = ems->f.arrayRefWithPadding(); + auto f = ems->f.view().forceWithPadding(); constr->apply(needsLogging, computeEnergy, count, 0, 1.0, ems->s.x.arrayRefWithPadding(), f, f.unpaddedArrayRef(), ems->s.box, ems->s.lambda[efptBONDED], &dvdl_constr, gmx::ArrayRefWithPadding(), computeVirial, shake_vir, @@ -920,16 +920,16 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres, static double reorder_partsum(const t_commrec* cr, t_grpopts* opts, const gmx_mtop_t* top_global, - em_state_t* s_min, - em_state_t* s_b) + const em_state_t* s_min, + const em_state_t* s_b) { if (debug) { fprintf(debug, "Doing reorder_partsum\n"); } - const rvec* fm = s_min->f.rvec_array(); - const rvec* fb = s_b->f.rvec_array(); + auto fm = s_min->f.view().force(); + auto fb = s_b->f.view().force(); /* Collect fm in a global vector fmg. * This conflicts with the spirit of domain decomposition, @@ -982,8 +982,8 @@ static real pr_beta(const t_commrec* cr, t_grpopts* opts, t_mdatoms* mdatoms, const gmx_mtop_t* top_global, - em_state_t* s_min, - em_state_t* s_b) + const em_state_t* s_min, + const em_state_t* s_b) { double sum; @@ -995,10 +995,10 @@ static real pr_beta(const t_commrec* cr, if (!DOMAINDECOMP(cr) || (s_min->s.ddp_count == cr->dd->ddp_count && s_b->s.ddp_count == cr->dd->ddp_count)) { - const rvec* fm = s_min->f.rvec_array(); - const rvec* fb = s_b->f.rvec_array(); - sum = 0; - int gf = 0; + auto fm = s_min->f.view().force(); + auto fb = s_b->f.view().force(); + sum = 0; + int gf = 0; /* This part of code can be incorrect with DD, * since the atom ordering in s_b and s_min might differ. */ @@ -1159,10 +1159,10 @@ void LegacySimulator::do_cg() */ /* Calculate the new direction in p, and the gradient in this direction, gpa */ - rvec* pm = s_min->s.cg_p.rvec_array(); - const rvec* sfm = s_min->f.rvec_array(); - double gpa = 0; - int gf = 0; + gmx::ArrayRef pm = s_min->s.cg_p; + gmx::ArrayRef sfm = s_min->f.view().force(); + double gpa = 0; + int gf = 0; for (int i = 0; i < mdatoms->homenr; i++) { if (mdatoms->cFREEZE) @@ -1279,16 +1279,17 @@ void LegacySimulator::do_cg() } /* Take a trial step (new coords in s_c) */ - do_em_step(cr, inputrec, mdatoms, s_min, c, &s_min->s.cg_p, s_c, constr, -1); + do_em_step(cr, inputrec, mdatoms, s_min, c, s_min->s.cg_p.constArrayRefWithPadding(), s_c, + constr, -1); neval++; /* Calculate energy for the trial step */ energyEvaluator.run(s_c, mu_tot, vir, pres, -1, FALSE); /* Calc derivative along line */ - const rvec* pc = s_c->s.cg_p.rvec_array(); - const rvec* sfc = s_c->f.rvec_array(); - double gpc = 0; + const rvec* pc = s_c->s.cg_p.rvec_array(); + gmx::ArrayRef sfc = s_c->f.view().force(); + double gpc = 0; for (int i = 0; i < mdatoms->homenr; i++) { for (m = 0; m < DIM; m++) @@ -1380,7 +1381,8 @@ void LegacySimulator::do_cg() } /* Take a trial step to this new point - new coords in s_b */ - do_em_step(cr, inputrec, mdatoms, s_min, b, &s_min->s.cg_p, s_b, constr, -1); + do_em_step(cr, inputrec, mdatoms, s_min, b, + s_min->s.cg_p.constArrayRefWithPadding(), s_b, constr, -1); neval++; /* Calculate energy for the trial step */ @@ -1389,9 +1391,9 @@ void LegacySimulator::do_cg() /* p does not change within a step, but since the domain decomposition * might change, we have to use cg_p of s_b here. */ - const rvec* pb = s_b->s.cg_p.rvec_array(); - const rvec* sfb = s_b->f.rvec_array(); - gpb = 0; + const rvec* pb = s_b->s.cg_p.rvec_array(); + gmx::ArrayRef sfb = s_b->f.view().force(); + gpb = 0; for (int i = 0; i < mdatoms->homenr; i++) { for (m = 0; m < DIM; m++) @@ -1803,7 +1805,7 @@ void LegacySimulator::do_lbfgs() point = 0; // Set initial search direction to the force (-gradient), or 0 for frozen particles. - real* fInit = static_cast(ems.f.rvec_array()[0]); + real* fInit = static_cast(ems.f.view().force().data()[0]); for (i = 0; i < n; i++) { if (!frozen[i]) @@ -1856,7 +1858,7 @@ void LegacySimulator::do_lbfgs() mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags, top_global->natoms, step, static_cast(step), &ems.s, state_global, - observablesHistory, ems.f); + observablesHistory, ems.f.view().force()); /* Do the linesearching in the direction dx[point][0..(n-1)] */ @@ -1864,7 +1866,7 @@ void LegacySimulator::do_lbfgs() s = dx[point]; real* xx = static_cast(ems.s.x.rvec_array()[0]); - real* ff = static_cast(ems.f.rvec_array()[0]); + real* ff = static_cast(ems.f.view().force().data()[0]); // calculate line gradient in position A for (gpa = 0, i = 0; i < n; i++) @@ -1896,7 +1898,7 @@ void LegacySimulator::do_lbfgs() // Before taking any steps along the line, store the old position *last = ems; real* lastx = static_cast(last->s.x.data()[0]); - real* lastf = static_cast(last->f.data()[0]); + real* lastf = static_cast(last->f.view().force().data()[0]); Epot0 = ems.epot; *sa = ems; @@ -1968,7 +1970,7 @@ void LegacySimulator::do_lbfgs() energyEvaluator.run(sc, mu_tot, vir, pres, step, FALSE); // Calc line gradient in position C - real* fc = static_cast(sc->f.rvec_array()[0]); + real* fc = static_cast(sc->f.view().force()[0]); for (gpc = 0, i = 0; i < n; i++) { gpc -= s[i] * fc[i]; /* f is negative gradient, thus the sign */ @@ -2050,7 +2052,7 @@ void LegacySimulator::do_lbfgs() fnorm = sb->fnorm; // Calculate gradient in point B - real* fb = static_cast(sb->f.rvec_array()[0]); + real* fb = static_cast(sb->f.view().force()[0]); for (gpb = 0, i = 0; i < n; i++) { gpb -= s[i] * fb[i]; /* f is negative gradient, thus the sign */ @@ -2428,7 +2430,8 @@ void LegacySimulator::do_steep() bool validStep = true; if (count > 0) { - validStep = do_em_step(cr, inputrec, mdatoms, s_min, stepsize, &s_min->f, s_try, constr, count); + validStep = do_em_step(cr, inputrec, mdatoms, s_min, stepsize, + s_min->f.view().forceWithPadding(), s_try, constr, count); } if (validStep) @@ -2730,7 +2733,7 @@ void LegacySimulator::do_nm() /* Steps are divided one by one over the nodes */ bool bNS = true; auto state_work_x = makeArrayRef(state_work.s.x); - auto state_work_f = makeArrayRef(state_work.f); + auto state_work_f = state_work.f.view().force(); for (index aid = cr->nodeid; aid < ssize(atom_index); aid += nnodes) { size_t atom = atom_index[aid]; @@ -2758,13 +2761,13 @@ void LegacySimulator::do_nm() if (shellfc) { /* 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, 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, - fr, runScheduleWork, t, mu_tot, vsite, DDBalanceRegionHandler(nullptr)); + relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, nullptr, step, inputrec, + imdSession, 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.view(), + vir, mdatoms, nrnb, wcycle, shellfc, fr, runScheduleWork, t, + mu_tot, vsite, DDBalanceRegionHandler(nullptr)); bNS = false; step++; } diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index 3cbd2d91bc..8f85eb8227 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -102,6 +102,7 @@ #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/df_history.h" #include "gromacs/mdtypes/energyhistory.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" @@ -171,20 +172,20 @@ void gmx::LegacySimulator::do_rerun() // alias to avoid a large ripple of nearly useless changes. // t_inputrec is being replaced by IMdpOptionsProvider, so this // will go away eventually. - t_inputrec* ir = inputrec; - int64_t step, step_rel; - double t, lam0[efptNR]; - bool isLastStep = false; - bool doFreeEnergyPerturbation = false; - unsigned int force_flags; - tensor force_vir, shake_vir, total_vir, pres; - t_trxstatus* status = nullptr; - rvec mu_tot; - t_trxframe rerun_fr; - gmx_localtop_t top(top_global->ffparams); - PaddedHostVector f{}; - gmx_global_stat_t gstat; - gmx_shellfc_t* shellfc; + t_inputrec* ir = inputrec; + int64_t step, step_rel; + double t, lam0[efptNR]; + bool isLastStep = false; + bool doFreeEnergyPerturbation = false; + unsigned int force_flags; + tensor force_vir, shake_vir, total_vir, pres; + t_trxstatus* status = nullptr; + rvec mu_tot; + t_trxframe rerun_fr; + gmx_localtop_t top(top_global->ffparams); + ForceBuffers f; + gmx_global_stat_t gstat; + gmx_shellfc_t* shellfc; double cycles; @@ -524,8 +525,8 @@ void gmx::LegacySimulator::do_rerun() relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir, 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, + state->v.arrayRefWithPadding(), state->box, state->lambda, + &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc, fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); } else @@ -539,7 +540,7 @@ void gmx::LegacySimulator::do_rerun() gmx_edsam* ed = nullptr; do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, + &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler); } @@ -552,8 +553,8 @@ void gmx::LegacySimulator::do_rerun() const bool bSumEkinhOld = false; do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, state_global, observablesHistory, top_global, fr, outf, - energyOutput, ekind, f, isCheckpointingStep, doRerun, - isLastStep, mdrunOptions.writeConfout, bSumEkinhOld); + energyOutput, ekind, f.view().force(), isCheckpointingStep, + doRerun, isLastStep, mdrunOptions.writeConfout, bSumEkinhOld); } stopHandler->setSignal(); diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index c33a68d3c2..d08e7f8660 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -65,6 +65,7 @@ #include "gromacs/mdlib/vsite.h" #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" @@ -916,7 +917,7 @@ void relax_shell_flexcon(FILE* fplog, const matrix box, ArrayRef lambda, history_t* hist, - ArrayRefWithPadding f, + gmx::ForceBuffersView* f, tensor force_vir, const t_mdatoms* md, t_nrnb* nrnb, @@ -1018,10 +1019,11 @@ void relax_shell_flexcon(FILE* fplog, { pr_rvecs(debug, 0, "x b4 do_force", as_rvec_array(x.data()), homenr); } - int shellfc_flags = force_flags | (bVerbose ? GMX_FORCE_ENERGY : 0); + int shellfc_flags = force_flags | (bVerbose ? GMX_FORCE_ENERGY : 0); + gmx::ForceBuffersView forceViewInit = gmx::ForceBuffersView(forceWithPadding[Min]); do_force(fplog, cr, ms, inputrec, nullptr, enforcedRotation, imdSession, pull_work, mdstep, - nrnb, wcycle, top, box, xPadded, hist, forceWithPadding[Min], force_vir, md, enerd, - lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, + nrnb, wcycle, top, box, xPadded, hist, &forceViewInit, force_vir, md, enerd, lambda, + fr, runScheduleWork, vsite, mu_tot, t, nullptr, (bDoNS ? GMX_FORCE_NS : 0) | shellfc_flags, ddBalanceRegionHandler); sf_dir = 0; @@ -1105,10 +1107,10 @@ void relax_shell_flexcon(FILE* fplog, pr_rvecs(debug, 0, "RELAX: pos[Try] ", as_rvec_array(pos[Try].data()), homenr); } /* 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, lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, shellfc_flags, - ddBalanceRegionHandler); + gmx::ForceBuffersView forceViewTry = gmx::ForceBuffersView(forceWithPadding[Try]); + do_force(fplog, cr, ms, inputrec, nullptr, enforcedRotation, imdSession, pull_work, 1, nrnb, wcycle, + top, box, posWithPadding[Try], hist, &forceViewTry, force_vir, md, enerd, lambda, fr, + runScheduleWork, vsite, mu_tot, t, nullptr, shellfc_flags, ddBalanceRegionHandler); accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals); if (gmx_debug_at) { @@ -1202,7 +1204,7 @@ void relax_shell_flexcon(FILE* fplog, /* Copy back the coordinates and the forces */ std::copy(pos[Min].begin(), pos[Min].end(), x.data()); - std::copy(force[Min].begin(), force[Min].end(), f.unpaddedArrayRef().begin()); + std::copy(force[Min].begin(), force[Min].end(), f->force().begin()); } void done_shellfc(FILE* fplog, gmx_shellfc_t* shfc, int64_t numSteps) diff --git a/src/gromacs/mdrun/shellfc.h b/src/gromacs/mdrun/shellfc.h index 1a6248334b..d14646b646 100644 --- a/src/gromacs/mdrun/shellfc.h +++ b/src/gromacs/mdrun/shellfc.h @@ -66,6 +66,7 @@ class ArrayRef; template class ArrayRefWithPadding; class Constraints; +class ForceBuffersView; class ImdSession; class MdrunScheduleWorkload; class VirtualSitesHandler; @@ -101,7 +102,7 @@ void relax_shell_flexcon(FILE* log, const matrix box, gmx::ArrayRef lambda, history_t* hist, - gmx::ArrayRefWithPadding f, + gmx::ForceBuffersView* f, tensor force_vir, const t_mdatoms* md, t_nrnb* nrnb, diff --git a/src/gromacs/mdrun/tpi.cpp b/src/gromacs/mdrun/tpi.cpp index 3e780cade8..7536de08eb 100644 --- a/src/gromacs/mdrun/tpi.cpp +++ b/src/gromacs/mdrun/tpi.cpp @@ -76,6 +76,7 @@ #include "gromacs/mdlib/vsite.h" #include "gromacs/mdrunutility/printtime.h" #include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" @@ -164,33 +165,33 @@ void LegacySimulator::do_tpi() { GMX_RELEASE_ASSERT(gmx_omp_nthreads_get(emntDefault) == 1, "TPI does not support OpenMP"); - gmx_localtop_t top(top_global->ffparams); - PaddedHostVector f{}; - real lambda, t, temp, beta, drmax, epot; - double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all; - t_trxstatus* status; - t_trxframe rerun_fr; - gmx_bool bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS; - tensor force_vir, shake_vir, vir, pres; - int a_tp0, a_tp1, ngid, gid_tp, nener, e; - rvec* x_mol; - rvec mu_tot, x_init, dx; - int nnodes, frame; - int64_t frame_step_prev, frame_step; - int64_t nsteps, stepblocksize = 0, step; - int64_t seed; - int i; - FILE* fp_tpi = nullptr; - char * ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN]; - double dbl, dump_ener; - gmx_bool bCavity; - int nat_cavity = 0, d; - real * mass_cavity = nullptr, mass_tot; - int nbin; - double invbinw, *bin, refvolshift, logV, bUlogV; - gmx_bool bEnergyOutOfBounds; - const char* tpid_leg[2] = { "direct", "reweighted" }; - auto mdatoms = mdAtoms->mdatoms(); + gmx_localtop_t top(top_global->ffparams); + gmx::ForceBuffers f; + real lambda, t, temp, beta, drmax, epot; + double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all; + t_trxstatus* status; + t_trxframe rerun_fr; + gmx_bool bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS; + tensor force_vir, shake_vir, vir, pres; + int a_tp0, a_tp1, ngid, gid_tp, nener, e; + rvec* x_mol; + rvec mu_tot, x_init, dx; + int nnodes, frame; + int64_t frame_step_prev, frame_step; + int64_t nsteps, stepblocksize = 0, step; + int64_t seed; + int i; + FILE* fp_tpi = nullptr; + char * ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN]; + double dbl, dump_ener; + gmx_bool bCavity; + int nat_cavity = 0, d; + real * mass_cavity = nullptr, mass_tot; + int nbin; + double invbinw, *bin, refvolshift, logV, bUlogV; + gmx_bool bEnergyOutOfBounds; + const char* tpid_leg[2] = { "direct", "reweighted" }; + auto mdatoms = mdAtoms->mdatoms(); GMX_UNUSED_VALUE(outputProvider); @@ -300,7 +301,7 @@ void LegacySimulator::do_tpi() atoms2md(top_global, inputrec, -1, {}, top_global->natoms, mdAtoms); update_mdatoms(mdatoms, inputrec->fepvals->init_lambda); - f.resizeWithPadding(top_global->natoms); + f.resize(top_global->natoms); /* Print to log file */ walltime_accounting_start_time(walltime_accounting); @@ -754,7 +755,7 @@ void LegacySimulator::do_tpi() std::feholdexcept(&floatingPointEnvironment); do_force(fplog, cr, ms, inputrec, nullptr, nullptr, imdSession, pull_work, step, nrnb, wcycle, &top, state_global->box, state_global->x.arrayRefWithPadding(), - &state_global->hist, f.arrayRefWithPadding(), force_vir, mdatoms, enerd, + &state_global->hist, &f.view(), force_vir, mdatoms, enerd, state_global->lambda, fr, runScheduleWork, nullptr, mu_tot, t, nullptr, GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY | (bStateChanged ? GMX_FORCE_STATECHANGED : 0), DDBalanceRegionHandler(nullptr)); diff --git a/src/gromacs/mdtypes/CMakeLists.txt b/src/gromacs/mdtypes/CMakeLists.txt index f13e63c8d8..0c4ecb5ab3 100644 --- a/src/gromacs/mdtypes/CMakeLists.txt +++ b/src/gromacs/mdtypes/CMakeLists.txt @@ -34,6 +34,7 @@ file(GLOB MDTYPES_SOURCES df_history.cpp + forcebuffers.cpp group.cpp iforceprovider.cpp inputrec.cpp @@ -52,6 +53,9 @@ else() ) endif() +if (BUILD_TESTING) + add_subdirectory(tests) +endif() set(LIBGROMACS_SOURCES ${LIBGROMACS_SOURCES} ${MDTYPES_SOURCES} PARENT_SCOPE) diff --git a/src/gromacs/mdtypes/forcebuffers.cpp b/src/gromacs/mdtypes/forcebuffers.cpp new file mode 100644 index 0000000000..e76c02ad43 --- /dev/null +++ b/src/gromacs/mdtypes/forcebuffers.cpp @@ -0,0 +1,78 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 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. + * + * 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. + */ +/*! \internal \file + * \brief + * Implements the ForceBuffers class + * + * \author Berk Hess + * \ingroup module_mdtypes + */ + +#include "gmxpre.h" + +#include "forcebuffers.h" + +#include "gromacs/gpu_utils/hostallocator.h" + +namespace gmx +{ + +ForceBuffers::ForceBuffers(PinningPolicy pinningPolicy) : force_({}, { pinningPolicy }), view_({}) +{ +} + +ForceBuffers::~ForceBuffers() = default; + +ForceBuffers& ForceBuffers::operator=(ForceBuffers const& o) +{ + auto oForce = o.view().force(); + resize(oForce.size()); + std::copy(oForce.begin(), oForce.end(), view().force().begin()); + + return *this; +} + +PinningPolicy ForceBuffers::pinningPolicy() const +{ + return force_.get_allocator().pinningPolicy(); +} + +void ForceBuffers::resize(int numAtoms) +{ + force_.resizeWithPadding(numAtoms); + view_ = ForceBuffersView(force_.arrayRefWithPadding()); +} + +} // namespace gmx diff --git a/src/gromacs/mdtypes/forcebuffers.h b/src/gromacs/mdtypes/forcebuffers.h new file mode 100644 index 0000000000..2d20ff2a2b --- /dev/null +++ b/src/gromacs/mdtypes/forcebuffers.h @@ -0,0 +1,150 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 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. + * + * 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. + */ + +/*! \libinternal \file + * + * \brief + * This file contains the definition of a container for force buffers. + * + * \author Berk Hess + * + * \inlibraryapi + * \ingroup module_mdtypes + */ + +#ifndef GMX_MDTYPES_FORCEBUFFERS_H +#define GMX_MDTYPES_FORCEBUFFERS_H + +#include "gromacs/gpu_utils/hostallocator.h" +#include "gromacs/math/arrayrefwithpadding.h" +#include "gromacs/math/vectypes.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/classhelpers.h" + +namespace gmx +{ + +enum class PinningPolicy : int; + +/*! \libinternal \brief A view of the force buffer + * + * This is used for read and/or write access to the force buffer. + */ +class ForceBuffersView +{ +public: + //! Constructor, creates a view to \p force + ForceBuffersView(const ArrayRefWithPadding& force) : force_(force) {} + + //! Copy constructor deleted to avoid creating non-const from const + ForceBuffersView(const ForceBuffersView& o) = delete; + + //! Move constructor deleted, but could be implemented + ForceBuffersView(ForceBuffersView&& o) = delete; + + ~ForceBuffersView() = default; + + //! Copy assignment deleted to avoid creating non-const from const + ForceBuffersView& operator=(const ForceBuffersView& v) = delete; + + //! Move assignment, moves the view + ForceBuffersView& operator=(ForceBuffersView&& v) = default; + + //! Returns a const arrayref to the force buffer without padding + ArrayRef force() const { return force_.unpaddedConstArrayRef(); } + + //! Returns an arrayref to the force buffer without padding + ArrayRef force() { return force_.unpaddedArrayRef(); } + + //! Returns an ArrayRefWithPadding to the force buffer + ArrayRefWithPadding forceWithPadding() { return force_; } + +private: + //! The force buffer + ArrayRefWithPadding force_; +}; + +/*! \libinternal \brief Object that holds the force buffers + * + * More buffers can be added when needed. Those should also be added + * to ForceBuffersView. + * Can be pinned for efficient transfer to/from GPUs. + * All access happens through the ForceBuffersView object. + */ +class ForceBuffers +{ +public: + //! Constructor, creates an empty force buffer with pinning not active + ForceBuffers(PinningPolicy pinningPolicy = PinningPolicy::CannotBePinned); + + //! Copy constructor deleted, but could be implemented + ForceBuffers(const ForceBuffers& o) = delete; + + //! Move constructor deleted, but could be implemented + ForceBuffers(ForceBuffers&& o) = delete; + + ~ForceBuffers(); + + //! Copy assignment operator, sets the pinning policy to CannotBePinned + ForceBuffers& operator=(ForceBuffers const& o); + + //! Move assignment operator, deleted but could be implemented + ForceBuffers& operator=(ForceBuffers&& o) = delete; + + //! Returns a const view to the force buffers + const ForceBuffersView& view() const { return view_; } + + //! Returns a view to the force buffer + ForceBuffersView& view() { return view_; } + + //! Resize the force buffer + void resize(int numAtoms); + + /*! \brief Returns the active pinning policy. + * + * Does not throw. + */ + PinningPolicy pinningPolicy() const; + +private: + //! The force buffer + PaddedHostVector force_; + //! The view to the force buffer + ForceBuffersView view_; +}; + +} // namespace gmx + +#endif diff --git a/src/gromacs/mdtypes/tests/CMakeLists.txt b/src/gromacs/mdtypes/tests/CMakeLists.txt new file mode 100644 index 0000000000..151938839a --- /dev/null +++ b/src/gromacs/mdtypes/tests/CMakeLists.txt @@ -0,0 +1,38 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright (c) 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. +# +# 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. + +gmx_add_unit_test(MdtypesUnitTest mdtypes-test + CPP_SOURCE_FILES + forcebuffers.cpp + ) diff --git a/src/gromacs/mdtypes/tests/forcebuffers.cpp b/src/gromacs/mdtypes/tests/forcebuffers.cpp new file mode 100644 index 0000000000..d23b91023c --- /dev/null +++ b/src/gromacs/mdtypes/tests/forcebuffers.cpp @@ -0,0 +1,133 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 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. + * + * 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. + */ +/*! \internal \file + * \brief + * Tests for the ForceBuffers and ForceBuffersView classes. + * + * \author berk Hess + * \ingroup module_mdtypes + */ +#include "gmxpre.h" + +#include "gromacs/mdtypes/forcebuffers.h" + +#include + +#include +#include + +#include "testutils/testasserts.h" + +namespace gmx +{ + +const std::array c_forces = { { { 0.5, 0.1, 1.2 }, { -2.1, 0.2, 0.3 } } }; + +TEST(ForceBuffers, ConstructsUnpinned) +{ + ForceBuffers forceBuffers; + + EXPECT_EQ(forceBuffers.pinningPolicy(), PinningPolicy::CannotBePinned); +} + +TEST(ForceBuffers, ConstructsPinned) +{ + ForceBuffers forceBuffers(PinningPolicy::PinnedIfSupported); + + EXPECT_EQ(forceBuffers.pinningPolicy(), PinningPolicy::PinnedIfSupported); +} + +TEST(ForceBuffers, ConstructsEmpty) +{ + ForceBuffers forceBuffers; + + EXPECT_EQ(forceBuffers.view().force().size(), 0); +} + +TEST(ForceBuffers, ResizeWorks) +{ + ForceBuffers forceBuffers; + + forceBuffers.resize(2); + EXPECT_EQ(forceBuffers.view().force().size(), 2); +} + +TEST(ForceBuffers, PaddingWorks) +{ + ForceBuffers forceBuffers; + + forceBuffers.resize(2); + auto paddedRef = forceBuffers.view().forceWithPadding(); + EXPECT_EQ(paddedRef.unpaddedArrayRef().size(), 2); + EXPECT_GT(paddedRef.size(), 2); +} + + +TEST(ForceBuffers, CopyWorks) +{ + ForceBuffers forceBuffers; + + forceBuffers.resize(2); + auto force = forceBuffers.view().force(); + index i = 0; + for (RVec& v : force) + { + v = c_forces[i]; + i++; + } + + ForceBuffers forceBuffersCopy; + forceBuffersCopy = forceBuffers; + auto forceCopy = forceBuffersCopy.view().force(); + EXPECT_EQ(forceBuffersCopy.view().force().size(), 2); + for (index i = 0; i < ssize(forceCopy); i++) + { + for (int d = 0; d < DIM; d++) + { + EXPECT_EQ(forceCopy[i][d], force[i][d]); + } + } +} + +TEST(ForceBuffers, CopyDoesNotPin) +{ + ForceBuffers forceBuffers(PinningPolicy::PinnedIfSupported); + + ForceBuffers forceBuffersCopy; + forceBuffersCopy = forceBuffers; + EXPECT_EQ(forceBuffersCopy.pinningPolicy(), PinningPolicy::CannotBePinned); +} + +} // namespace gmx diff --git a/src/gromacs/modularsimulator/domdechelper.cpp b/src/gromacs/modularsimulator/domdechelper.cpp index 6aa03d656c..4d21b62d10 100644 --- a/src/gromacs/modularsimulator/domdechelper.cpp +++ b/src/gromacs/modularsimulator/domdechelper.cpp @@ -100,7 +100,7 @@ void DomDecHelper::setup() { std::unique_ptr localState = statePropagatorData_->localState(); t_state* globalState = statePropagatorData_->globalState(); - PaddedHostVector* forcePointer = statePropagatorData_->forcePointer(); + ForceBuffers* forcePointer = statePropagatorData_->forcePointer(); // constant choices for this call to dd_partition_system const bool verbose = false; @@ -127,7 +127,7 @@ void DomDecHelper::run(Step step, Time gmx_unused time) } std::unique_ptr localState = statePropagatorData_->localState(); t_state* globalState = statePropagatorData_->globalState(); - PaddedHostVector* forcePointer = statePropagatorData_->forcePointer(); + ForceBuffers* forcePointer = statePropagatorData_->forcePointer(); // constant choices for this call to dd_partition_system const bool verbose = isVerbose_ && (step % verbosePrintInterval_ == 0 || step == inputrec_->init_step); diff --git a/src/gromacs/modularsimulator/forceelement.cpp b/src/gromacs/modularsimulator/forceelement.cpp index aa97fe72e8..4d43376706 100644 --- a/src/gromacs/modularsimulator/forceelement.cpp +++ b/src/gromacs/modularsimulator/forceelement.cpp @@ -50,6 +50,7 @@ #include "gromacs/mdlib/force_flags.h" #include "gromacs/mdlib/mdatoms.h" #include "gromacs/mdrun/shellfc.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/mdatom.h" #include "gromacs/mdtypes/mdrunoptions.h" @@ -178,7 +179,7 @@ void ForceElement::run(Step step, Time time, unsigned int flags) * Check comments in sim_util.c */ auto x = statePropagatorData_->positionsView(); - auto forces = statePropagatorData_->forcesView(); + auto& forces = statePropagatorData_->forcesView(); auto box = statePropagatorData_->constBox(); history_t* hist = nullptr; // disabled @@ -195,7 +196,7 @@ void ForceElement::run(Step step, Time time, unsigned int flags) fplog_, cr_, ms, isVerbose_, enforcedRotation_, step, inputrec_, imdSession_, pull_work_, step == nextNSStep_, static_cast(flags), localTopology_, constr_, energyData_->enerdata(), statePropagatorData_->localNumAtoms(), x, v, box, lambda, - hist, forces, force_vir, mdAtoms_->mdatoms(), nrnb_, wcycle_, shellfc_, fr_, + hist, &forces, force_vir, mdAtoms_->mdatoms(), nrnb_, wcycle_, shellfc_, fr_, runScheduleWork_, time, energyData_->muTot(), vsite_, ddBalanceRegionHandler_); nShellRelaxationSteps_++; } @@ -206,7 +207,7 @@ void ForceElement::run(Step step, Time time, unsigned int flags) gmx_edsam* ed = nullptr; do_force(fplog_, cr_, ms, inputrec_, awh, enforcedRotation_, imdSession_, pull_work_, step, - nrnb_, wcycle_, localTopology_, box, x, hist, forces, force_vir, + nrnb_, wcycle_, localTopology_, box, x, hist, &forces, force_vir, mdAtoms_->mdatoms(), energyData_->enerdata(), lambda, fr_, runScheduleWork_, vsite_, energyData_->muTot(), time, ed, static_cast(flags), ddBalanceRegionHandler_); } diff --git a/src/gromacs/modularsimulator/propagator.cpp b/src/gromacs/modularsimulator/propagator.cpp index 238bc6ca0b..fd626e1cbc 100644 --- a/src/gromacs/modularsimulator/propagator.cpp +++ b/src/gromacs/modularsimulator/propagator.cpp @@ -156,7 +156,7 @@ void Propagator::run() wallcycle_start(wcycle_, ewcUPDATE); auto v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data()); - auto f = as_rvec_array(statePropagatorData_->constForcesView().paddedArrayRef().data()); + auto f = as_rvec_array(statePropagatorData_->constForcesView().force().data()); auto invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim; const real lambda = @@ -253,7 +253,7 @@ void Propagator::run() auto xp = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data()); auto x = as_rvec_array(statePropagatorData_->constPreviousPositionsView().paddedArrayRef().data()); auto v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data()); - auto f = as_rvec_array(statePropagatorData_->constForcesView().paddedArrayRef().data()); + auto f = as_rvec_array(statePropagatorData_->constForcesView().force().data()); auto invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim; const real lambda = @@ -351,7 +351,7 @@ void Propagator::run() auto xp = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data()); auto x = as_rvec_array(statePropagatorData_->constPreviousPositionsView().paddedArrayRef().data()); auto v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data()); - auto f = as_rvec_array(statePropagatorData_->constForcesView().paddedArrayRef().data()); + auto f = as_rvec_array(statePropagatorData_->constForcesView().force().data()); auto invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim; int nth = gmx_omp_nthreads_get(emntUpdate); diff --git a/src/gromacs/modularsimulator/statepropagatordata.cpp b/src/gromacs/modularsimulator/statepropagatordata.cpp index de554c1884..5ea269c272 100644 --- a/src/gromacs/modularsimulator/statepropagatordata.cpp +++ b/src/gromacs/modularsimulator/statepropagatordata.cpp @@ -54,6 +54,7 @@ #include "gromacs/mdlib/stat.h" #include "gromacs/mdlib/update.h" #include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/mdatom.h" @@ -115,7 +116,7 @@ StatePropagatorData::StatePropagatorData(int numAtoms, else { state_change_natoms(globalState, globalState->natoms); - f_.resizeWithPadding(globalState->natoms); + f_.resize(globalState->natoms); localNAtoms_ = globalState->natoms; x_ = globalState->x; v_ = globalState->v; @@ -204,14 +205,14 @@ ArrayRefWithPadding StatePropagatorData::constVelocitiesView() const return v_.constArrayRefWithPadding(); } -ArrayRefWithPadding StatePropagatorData::forcesView() +ForceBuffersView& StatePropagatorData::forcesView() { - return f_.arrayRefWithPadding(); + return f_.view(); } -ArrayRefWithPadding StatePropagatorData::constForcesView() const +const ForceBuffersView& StatePropagatorData::constForcesView() const { - return f_.constArrayRefWithPadding(); + return f_.view(); } rvec* StatePropagatorData::box() @@ -286,7 +287,7 @@ t_state* StatePropagatorData::globalState() return globalState_; } -PaddedHostVector* StatePropagatorData::forcePointer() +ForceBuffers* StatePropagatorData::forcePointer() { return &f_; } @@ -426,7 +427,7 @@ void StatePropagatorData::Element::write(gmx_mdoutf_t outf, Step currentStep, Ti mdoutf_write_to_trajectory_files(fplog_, cr_, outf, static_cast(mdof_flags), statePropagatorData_->totalNumAtoms_, currentStep, currentTime, localStateBackup_.get(), statePropagatorData_->globalState_, - observablesHistory, statePropagatorData_->f_); + observablesHistory, statePropagatorData_->f_.view().force()); if (currentStep != lastStep_ || !isRegularSimulationEnd_) { diff --git a/src/gromacs/modularsimulator/statepropagatordata.h b/src/gromacs/modularsimulator/statepropagatordata.h index f702fa1608..38713f2629 100644 --- a/src/gromacs/modularsimulator/statepropagatordata.h +++ b/src/gromacs/modularsimulator/statepropagatordata.h @@ -47,6 +47,7 @@ #include "gromacs/gpu_utils/hostallocator.h" #include "gromacs/math/paddedvector.h" #include "gromacs/math/vectypes.h" +#include "gromacs/mdtypes/forcebuffers.h" #include "modularsimulatorinterfaces.h" #include "topologyholder.h" @@ -121,9 +122,9 @@ public: //! Get read access to velocity vector ArrayRefWithPadding constVelocitiesView() const; //! Get write access to force vector - ArrayRefWithPadding forcesView(); + ForceBuffersView& forcesView(); //! Get read access to force vector - ArrayRefWithPadding constForcesView() const; + const ForceBuffersView& constForcesView() const; //! Get pointer to box rvec* box(); //! Get const pointer to box @@ -162,7 +163,7 @@ private: //! The velocity vector PaddedHostVector v_; //! The force vector - PaddedHostVector f_; + ForceBuffers f_; //! The box matrix matrix box_; //! The box matrix of the previous step @@ -186,7 +187,7 @@ private: //! Get a pointer to the global state t_state* globalState(); //! Get a force pointer - PaddedHostVector* forcePointer(); + ForceBuffers* forcePointer(); //! Whether we're doing VV and need to reset velocities after the first half step bool vvResetVelocities_;