Add ForceBuffers class
authorBerk Hess <hess@kth.se>
Thu, 27 Aug 2020 06:09:17 +0000 (06:09 +0000)
committerBerk Hess <hess@kth.se>
Thu, 27 Aug 2020 06:09:17 +0000 (06:09 +0000)
This change replaces the type of the force buffer in mdrun from
PaddedHostVector<RVec> 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.

27 files changed:
src/gromacs/domdec/mdsetup.cpp
src/gromacs/domdec/mdsetup.h
src/gromacs/domdec/partition.cpp
src/gromacs/domdec/partition.h
src/gromacs/mdlib/force.h
src/gromacs/mdlib/mdoutf.cpp
src/gromacs/mdlib/mdoutf.h
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdlib/trajectory_writing.cpp
src/gromacs/mdlib/trajectory_writing.h
src/gromacs/mdrun/md.cpp
src/gromacs/mdrun/mimic.cpp
src/gromacs/mdrun/minimize.cpp
src/gromacs/mdrun/rerun.cpp
src/gromacs/mdrun/shellfc.cpp
src/gromacs/mdrun/shellfc.h
src/gromacs/mdrun/tpi.cpp
src/gromacs/mdtypes/CMakeLists.txt
src/gromacs/mdtypes/forcebuffers.cpp [new file with mode: 0644]
src/gromacs/mdtypes/forcebuffers.h [new file with mode: 0644]
src/gromacs/mdtypes/tests/CMakeLists.txt [new file with mode: 0644]
src/gromacs/mdtypes/tests/forcebuffers.cpp [new file with mode: 0644]
src/gromacs/modularsimulator/domdechelper.cpp
src/gromacs/modularsimulator/forceelement.cpp
src/gromacs/modularsimulator/propagator.cpp
src/gromacs/modularsimulator/statepropagatordata.cpp
src/gromacs/modularsimulator/statepropagatordata.h

index 9c71dd410cb2fbcb8a47857763ba6194448062b0..6a8559217f633a7300ebf1a9162cba2a4b3b87d6 100644 (file)
@@ -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<RVec>* 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,
index d406b09f7fb824d9f11659c9202ea26fe558d36b..ab0da8fcf16f398beeace4c119e42ea610c462a9 100644 (file)
@@ -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<RVec>* 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
 
index 2d5d6cb0b82c978d8224574f0aed4fccd7dd3378..c7a01ae473b785fee4bc134ef612d1ab82d9a8ce 100644 (file)
@@ -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<gmx::RVec>* 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;
index e06138a1b7c01b9e0d938d9ab6957793804671b3..de05c4a73847b7cebcb1f9cd78892d2074085e0d 100644 (file)
@@ -45,7 +45,8 @@
 #ifndef GMX_DOMDEC_PARTITION_H
 #define GMX_DOMDEC_PARTITION_H
 
-#include "gromacs/gpu_utils/hostallocator.h"
+#include <cstdio>
+
 #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<typename>
+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<gmx::RVec>* 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
  *
index 3db3f96db42a71ae9fb0d666f9a73de7c95b64cd..711ce856b88ae2b7ad1d434bed282707885bccd5 100644 (file)
@@ -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<gmx::RVec> coordinates,
               history_t*                          hist,
-              gmx::ArrayRefWithPadding<gmx::RVec> force,
+              gmx::ForceBuffersView*              force,
               tensor                              vir_force,
               const t_mdatoms*                    mdatoms,
               gmx_enerdata_t*                     enerd,
index 795aae408ba957e2a0f9ffa064e25cc457c07289..5b94dbd1d4c356c91d034bd46799ecd931da4e46 100644 (file)
@@ -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<gmx::RVec> 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<const gmx::RVec> 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<gmx::RVec*>(f_global), f_local.size()));
+            dd_collect_vec(
+                    cr->dd, state_local, f_local,
+                    gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(of->f_global), f_local.size()));
         }
     }
     else
index 4e704efecfce1da0c2eb82a7af7a6ed0ba829431..fe534fab60e1ab2149a1b7f131a93b9cfed09d4a 100644 (file)
@@ -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<gmx::RVec> 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<const gmx::RVec> f_local);
 
 /*! \brief Get the output interval of box size of uncompressed TNG output.
  * Returns 0 if no uncompressed TNG file is open.
index 3f2cd68dd8156f42dc4a4dc06cc21f7317336c68..dda01b711fd0cf8420644e157c3240f5fd4351bf 100644 (file)
@@ -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<gmx::RVec> x,
               history_t*                          hist,
-              gmx::ArrayRefWithPadding<gmx::RVec> 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");
index 9c87f46c4c237a84cbf8c13ef33ae17b91915c51..6d250ceeb2b30a8eb2c03a90596a92603e81301c 100644 (file)
 #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<gmx::RVec> 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<const gmx::RVec> f,
+                              gmx_bool                       bCPT,
+                              gmx_bool                       bRerunMD,
+                              gmx_bool                       bLastStep,
+                              gmx_bool                       bDoConfOut,
+                              gmx_bool                       bSumEkinhOld)
 {
     int   mdof_flags;
     rvec* x_for_confout = nullptr;
index bd17e22d26376839fcf8fe5c44489f9b251383ec..a74c707e778abf205d4f755850cfc15cda2eca03 100644 (file)
@@ -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<gmx::RVec> 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<const gmx::RVec> f,
+                              gmx_bool                       bCPT,
+                              gmx_bool                       bRerunMD,
+                              gmx_bool                       bLastStep,
+                              gmx_bool                       bDoConfOut,
+                              gmx_bool                       bSumEkinhOld);
 
 
 #endif
index 6e683a8dd89141cd0cd34cf978adb5729a6ebacc..ac08b56afc0c24752e833e6f335b619a6110c8a0 100644 (file)
 #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<gmx::RVec> 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<RVec>           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<RVec> 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<UpdateConstrainGpu> 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<t_state>();
@@ -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<UpdateConstrainGpu> 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<RVec>(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<RVec>(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? */
index 84bea6856b6f6215fbb32a8918198917834c0c7d..29d846b8515ac559ec732ed34b35034e8a2adc79 100644 (file)
 #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<gmx::RVec> 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<gmx::RVec>     fglobal(top_global->natoms);
             gmx::ArrayRef<gmx::RVec>       ftemp;
-            gmx::ArrayRef<const gmx::RVec> flocal = gmx::makeArrayRef(f);
+            gmx::ArrayRef<const gmx::RVec> 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))
index c96ef39acd12aa0c80e94a612e5d1b26d77fc633..7e6b030e3c9b177fa57966e13247a53da7d64017 100644 (file)
@@ -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<gmx::RVec> 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<const gmx::RVec> 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<double>(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<gmx::RVec>* 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<const gmx::RVec> 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<RVec>(), 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<gmx::RVec>       pm  = s_min->s.cg_p;
+        gmx::ArrayRef<const gmx::RVec> 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<const gmx::RVec> 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<const gmx::RVec> 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<real*>(ems.f.rvec_array()[0]);
+    real* fInit = static_cast<real*>(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<real>(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<real*>(ems.s.x.rvec_array()[0]);
-        real* ff = static_cast<real*>(ems.f.rvec_array()[0]);
+        real* ff = static_cast<real*>(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<real*>(last->s.x.data()[0]);
-        real* lastf = static_cast<real*>(last->f.data()[0]);
+        real* lastf = static_cast<real*>(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<real*>(sc->f.rvec_array()[0]);
+        real* fc = static_cast<real*>(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<real*>(sb->f.rvec_array()[0]);
+                real* fb = static_cast<real*>(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++;
                 }
index 3cbd2d91bc42e75eb513d4e3c0b613de5fae1b60..8f85eb82278b10afc3f9900f4d97e1786cb9c038 100644 (file)
 #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<gmx::RVec> 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();
index c33a68d3c26227a88f651039aca9ffdae26aac37..d08e7f86604549cb302847c9ce85516076c7d667 100644 (file)
@@ -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<real>                lambda,
                          history_t*                    hist,
-                         ArrayRefWithPadding<RVec>     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)
index 1a6248334b7f3ca5e3e994aa26218de8b5aa610a..d14646b6469593ddf8838bd10095551beb8948c4 100644 (file)
@@ -66,6 +66,7 @@ class ArrayRef;
 template<typename>
 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<real>                 lambda,
                          history_t*                          hist,
-                         gmx::ArrayRefWithPadding<gmx::RVec> f,
+                         gmx::ForceBuffersView*              f,
                          tensor                              force_vir,
                          const t_mdatoms*                    md,
                          t_nrnb*                             nrnb,
index 3e780cade8480300c8e0d5309d281ae602cafc19..7536de08ebf5c8a641f6f670309453d292ec4f64 100644 (file)
@@ -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<gmx::RVec> 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));
index f13e63c8d8b402d12c50eab52bf3ddf22bc94eab..0c4ecb5ab3feffdf9274b54563e65f7672b45f8a 100644 (file)
@@ -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 (file)
index 0000000..e76c02a
--- /dev/null
@@ -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 <hess@kth.se>
+ * \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 (file)
index 0000000..2d20ff2
--- /dev/null
@@ -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<RVec>& 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<const RVec> force() const { return force_.unpaddedConstArrayRef(); }
+
+    //! Returns an arrayref to the force buffer without padding
+    ArrayRef<RVec> force() { return force_.unpaddedArrayRef(); }
+
+    //! Returns an ArrayRefWithPadding to the force buffer
+    ArrayRefWithPadding<RVec> forceWithPadding() { return force_; }
+
+private:
+    //! The force buffer
+    ArrayRefWithPadding<RVec> 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<RVec> 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 (file)
index 0000000..1519388
--- /dev/null
@@ -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 (file)
index 0000000..d23b910
--- /dev/null
@@ -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 <hess@kth.se>
+ * \ingroup module_mdtypes
+ */
+#include "gmxpre.h"
+
+#include "gromacs/mdtypes/forcebuffers.h"
+
+#include <array>
+
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+#include "testutils/testasserts.h"
+
+namespace gmx
+{
+
+const std::array<RVec, 2> 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
index 6aa03d656c781c5ddb9c535abd7cb910054c232b..4d21b62d1000036d0dda8790071ace5cfe2390b3 100644 (file)
@@ -100,7 +100,7 @@ void DomDecHelper::setup()
 {
     std::unique_ptr<t_state> localState   = statePropagatorData_->localState();
     t_state*                 globalState  = statePropagatorData_->globalState();
-    PaddedHostVector<RVec>*  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<t_state> localState   = statePropagatorData_->localState();
     t_state*                 globalState  = statePropagatorData_->globalState();
-    PaddedHostVector<RVec>*  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);
index aa97fe72e85c56976897408a9ca30b78ab20842d..4d433767069cea38af23b0b684e9511a5f9bad55 100644 (file)
@@ -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<int>(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<int>(flags), ddBalanceRegionHandler_);
     }
index 238bc6ca0b844634a71d22e8aa8e926d35c882ce..fd626e1cbcf69932f45cdbb2cc5f229a19059f2d 100644 (file)
@@ -156,7 +156,7 @@ void Propagator<IntegrationStep::VelocitiesOnly>::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<IntegrationStep::LeapFrog>::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<IntegrationStep::VelocityVerletPositionsAndVelocities>::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);
index de554c1884997cd684f9657f4ee2224150426c9b..5ea269c272f211f67a7523917c710df3b3ffe48e 100644 (file)
@@ -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<const RVec> StatePropagatorData::constVelocitiesView() const
     return v_.constArrayRefWithPadding();
 }
 
-ArrayRefWithPadding<RVec> StatePropagatorData::forcesView()
+ForceBuffersView& StatePropagatorData::forcesView()
 {
-    return f_.arrayRefWithPadding();
+    return f_.view();
 }
 
-ArrayRefWithPadding<const RVec> 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<RVec>* 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<int>(mdof_flags),
                                      statePropagatorData_->totalNumAtoms_, currentStep, currentTime,
                                      localStateBackup_.get(), statePropagatorData_->globalState_,
-                                     observablesHistory, statePropagatorData_->f_);
+                                     observablesHistory, statePropagatorData_->f_.view().force());
 
     if (currentStep != lastStep_ || !isRegularSimulationEnd_)
     {
index f702fa1608f96dec6fce0f2a7cfd225d83b83631..38713f2629d118117f35b8b0a2fe58e0225cf2d1 100644 (file)
@@ -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<const RVec> constVelocitiesView() const;
     //! Get write access to force vector
-    ArrayRefWithPadding<RVec> forcesView();
+    ForceBuffersView& forcesView();
     //! Get read access to force vector
-    ArrayRefWithPadding<const RVec> 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<RVec> v_;
     //! The force vector
-    PaddedHostVector<RVec> 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<gmx::RVec>* forcePointer();
+    ForceBuffers* forcePointer();
 
     //! Whether we're doing VV and need to reset velocities after the first half step
     bool vvResetVelocities_;