Convert gmx_update_t to C++
authorKevin Boyd <kevin.boyd@uconn.edu>
Sat, 22 Dec 2018 15:21:07 +0000 (10:21 -0500)
committerDavid van der Spoel <spoel@xray.bmc.uu.se>
Mon, 21 Jan 2019 10:42:24 +0000 (11:42 +0100)
gmx_update_t replaced with Update class with pImpled
implementation class

Fixes an end of mdrun leak

Put Update on the stack in mdrun

Moved Update construction out of init_md

Removed unused deform parameters

Change-Id: I1fc98ac1ab7630461a75a8aee1a4af07afa58794

src/gromacs/mdlib/coupling.cpp
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdlib/sim_util.h
src/gromacs/mdlib/update.cpp
src/gromacs/mdlib/update.h
src/gromacs/mdrun/md.cpp

index 538be1dffc65472e15a034db5bba59e0ef7b70e1..d6689b89053e1843560d87f8dd51b3a810137367 100644 (file)
@@ -1619,7 +1619,7 @@ void rescale_velocities(const gmx_ekindata_t *ekind, const t_mdatoms *mdatoms,
 
 
 /* set target temperatures if we are annealing */
-void update_annealing_target_temp(t_inputrec *ir, real t, gmx_update_t *upd)
+void update_annealing_target_temp(t_inputrec *ir, real t, gmx::Update *upd)
 {
     int  i, j, n, npoints;
     real pert, thist = 0, x;
@@ -1679,5 +1679,5 @@ void update_annealing_target_temp(t_inputrec *ir, real t, gmx_update_t *upd)
         }
     }
 
-    update_temperature_constants(upd, ir);
+    update_temperature_constants(upd->sd(), ir);
 }
index 9c381d555932002f634af6a17c96457b9395cd67..c18287aad0f430c9f22bca64d15ba4ee7c64372c 100644 (file)
@@ -2938,8 +2938,7 @@ void init_md(FILE *fplog,
              double *t, double *t0,
              t_state *globalState, double *lam0,
              t_nrnb *nrnb, gmx_mtop_t *mtop,
-             gmx_update_t **upd,
-             gmx::BoxDeformation *deform,
+             gmx::Update *upd,
              int nfile, const t_filenm fnm[],
              gmx_mdoutf_t *outf, t_mdebin **mdebin,
              tensor force_vir, tensor shake_vir,
@@ -2978,14 +2977,9 @@ void init_md(FILE *fplog,
         initialize_lambdas(fplog, ir, &tmpFepState, tmpLambda, lam0);
     }
 
-    // TODO upd is never NULL in practice, but the analysers don't know that
-    if (upd)
-    {
-        *upd = init_update(ir, deform);
-    }
     if (*bSimAnn)
     {
-        update_annealing_target_temp(ir, ir->init_t, upd ? *upd : nullptr);
+        update_annealing_target_temp(ir, ir->init_t, upd);
     }
 
     if (EI_DYNAMICS(ir->eI) && !mdrunOptions.continuationOptions.appendFiles)
index 4ca378b5157f87c0ca625157fc27ccb5b5859600..f908a3184d458dcf42be2c111aacef7f10d672c1 100644 (file)
@@ -47,7 +47,6 @@
 
 struct gmx_output_env_t;
 struct gmx_pme_t;
-struct gmx_update_t;
 struct MdrunOptions;
 struct nonbonded_verlet_t;
 struct t_forcerec;
@@ -60,6 +59,7 @@ class BoxDeformation;
 class Constraints;
 class IMDOutputProvider;
 class MDLogger;
+class Update;
 }
 
 typedef struct gmx_global_stat *gmx_global_stat_t;
@@ -149,8 +149,7 @@ void init_md(FILE *fplog,
              double *t, double *t0,
              t_state *globalState, double *lam0,
              t_nrnb *nrnb, gmx_mtop_t *mtop,
-             gmx_update_t **upd,
-             gmx::BoxDeformation *deform,
+             gmx::Update *upd,
              int nfile, const t_filenm fnm[],
              gmx_mdoutf_t *outf, t_mdebin **mdebin,
              tensor force_vir, tensor shake_vir,
index c34571000f2094eadb316703d0f86844166d3785..e6f39df9bc17b0f77db5e4b1dc377f60747cd98c 100644 (file)
 
 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
 
-typedef struct {
+struct gmx_sd_const_t {
     double em;
-} gmx_sd_const_t;
+};
 
-typedef struct {
+struct gmx_sd_sigma_t {
     real V;
-} gmx_sd_sigma_t;
+};
 
 struct gmx_stochd_t
 {
@@ -104,22 +104,46 @@ struct gmx_stochd_t
     std::vector<bool>           randomize_group;
     std::vector<real>           boltzfac;
 
-    gmx_stochd_t(const t_inputrec *ir);
+    explicit gmx_stochd_t(const t_inputrec *ir);
 };
 
-struct gmx_update_t
+//! pImpled implementation for Update
+class Update::Impl
 {
-    std::unique_ptr<gmx_stochd_t> sd;
-    /* xprime for constraint algorithms */
-    PaddedVector<gmx::RVec>       xp;
+    public:
+        //! Constructor
+        Impl(const t_inputrec    *ir, BoxDeformation *boxDeformation);
+        //! Destructor
+        ~Impl() = default;
+        //! stochastic dynamics struct
+        std::unique_ptr<gmx_stochd_t> sd;
+        //! xprime for constraint algorithms
+        PaddedVector<RVec>            xp;
+        //! Box deformation handler (or nullptr if inactive).
+        BoxDeformation               *deform = nullptr;
+};
 
-    /* Variables for the deform algorithm */
-    int64_t           deformref_step;
-    matrix            deformref_box;
+Update::Update(const t_inputrec    *ir, BoxDeformation *boxDeformation)
+    : impl_(new Impl(ir, boxDeformation))
+{};
 
-    //! Box deformation handler (or nullptr if inactive).
-    gmx::BoxDeformation *deform;
-};
+Update::~Update()
+{};
+
+gmx_stochd_t* Update::sd() const
+{
+    return impl_->sd.get();
+}
+
+PaddedVector<RVec> * Update::xp()
+{
+    return &impl_->xp;
+}
+
+BoxDeformation * Update::deform() const
+{
+    return impl_->deform;
+}
 
 static bool isTemperatureCouplingStep(int64_t step, const t_inputrec *ir)
 {
@@ -819,7 +843,7 @@ gmx_stochd_t::gmx_stochd_t(const t_inputrec *ir)
     }
 }
 
-void update_temperature_constants(gmx_update_t *upd, const t_inputrec *ir)
+void update_temperature_constants(gmx_stochd_t *sd, const t_inputrec *ir)
 {
     if (ir->eI == eiBD)
     {
@@ -827,14 +851,14 @@ void update_temperature_constants(gmx_update_t *upd, const t_inputrec *ir)
         {
             for (int gt = 0; gt < ir->opts.ngtc; gt++)
             {
-                upd->sd->bd_rf[gt] = std::sqrt(2.0*BOLTZ*ir->opts.ref_t[gt]/(ir->bd_fric*ir->delta_t));
+                sd->bd_rf[gt] = std::sqrt(2.0*BOLTZ*ir->opts.ref_t[gt]/(ir->bd_fric*ir->delta_t));
             }
         }
         else
         {
             for (int gt = 0; gt < ir->opts.ngtc; gt++)
             {
-                upd->sd->bd_rf[gt] = std::sqrt(2.0*BOLTZ*ir->opts.ref_t[gt]);
+                sd->bd_rf[gt] = std::sqrt(2.0*BOLTZ*ir->opts.ref_t[gt]);
             }
         }
     }
@@ -844,32 +868,23 @@ void update_temperature_constants(gmx_update_t *upd, const t_inputrec *ir)
         {
             real kT = BOLTZ*ir->opts.ref_t[gt];
             /* The mass is accounted for later, since this differs per atom */
-            upd->sd->sdsig[gt].V  = std::sqrt(kT*(1 - upd->sd->sdc[gt].em*upd->sd->sdc[gt].em));
+            sd->sdsig[gt].V  = std::sqrt(kT*(1 - sd->sdc[gt].em * sd->sdc[gt].em));
         }
     }
 }
 
-gmx_update_t *init_update(const t_inputrec    *ir,
-                          gmx::BoxDeformation *deform)
+Update::Impl::Impl(const t_inputrec    *ir, BoxDeformation *boxDeformation)
 {
-    gmx_update_t *upd = new(gmx_update_t);
-
-    upd->sd    = gmx::compat::make_unique<gmx_stochd_t>(ir);
-
-    update_temperature_constants(upd, ir);
-
-    upd->xp.resizeWithPadding(0);
-
-    upd->deform = deform;
-
-    return upd;
+    sd = gmx::compat::make_unique<gmx_stochd_t>(ir);
+    update_temperature_constants(sd.get(), ir);
+    xp.resizeWithPadding(0);
+    deform = boxDeformation;
 }
 
-void update_realloc(gmx_update_t *upd, int natoms)
+void Update::setNumAtoms(int nAtoms)
 {
-    GMX_ASSERT(upd, "upd must be allocated before its fields can be reallocated");
 
-    upd->xp.resizeWithPadding(natoms);
+    impl_->xp.resizeWithPadding(nAtoms);
 }
 
 /*! \brief Sets the SD update type */
@@ -1478,15 +1493,15 @@ void constrain_velocities(int64_t                        step,
     }
 }
 
-void constrain_coordinates(int64_t                        step,
-                           real                          *dvdlambda, /* the contribution to be added to the bonded interactions */
-                           t_state                       *state,
-                           tensor                         vir_part,
-                           gmx_update_t                  *upd,
-                           gmx::Constraints              *constr,
-                           gmx_bool                       bCalcVir,
-                           bool                           do_log,
-                           bool                           do_ene)
+void constrain_coordinates(int64_t           step,
+                           real             *dvdlambda, /* the contribution to be added to the bonded interactions */
+                           t_state          *state,
+                           tensor            vir_part,
+                           Update           *upd,
+                           gmx::Constraints *constr,
+                           gmx_bool          bCalcVir,
+                           bool              do_log,
+                           bool              do_ene)
 {
     if (!constr)
     {
@@ -1502,7 +1517,7 @@ void constrain_coordinates(int64_t                        step,
         /* Constrain the coordinates upd->xp */
         constr->apply(do_log, do_ene,
                       step, 1, 1.0,
-                      state->x.rvec_array(), upd->xp.rvec_array(), nullptr,
+                      state->x.rvec_array(), upd->xp()->rvec_array(), nullptr,
                       state->box,
                       state->lambda[efptBONDED], dvdlambda,
                       as_rvec_array(state->v.data()), bCalcVir ? &vir_con : nullptr, ConstraintVariable::Positions);
@@ -1515,18 +1530,18 @@ void constrain_coordinates(int64_t                        step,
 }
 
 void
-update_sd_second_half(int64_t                        step,
-                      real                          *dvdlambda,   /* the contribution to be added to the bonded interactions */
-                      const t_inputrec              *inputrec,    /* input record and box stuff        */
-                      const t_mdatoms               *md,
-                      t_state                       *state,
-                      const t_commrec               *cr,
-                      t_nrnb                        *nrnb,
-                      gmx_wallcycle_t                wcycle,
-                      gmx_update_t                  *upd,
-                      gmx::Constraints              *constr,
-                      bool                           do_log,
-                      bool                           do_ene)
+update_sd_second_half(int64_t           step,
+                      real             *dvdlambda, /* the contribution to be added to the bonded interactions */
+                      const t_inputrec *inputrec,  /* input record and box stuff       */
+                      const t_mdatoms  *md,
+                      t_state          *state,
+                      const t_commrec  *cr,
+                      t_nrnb           *nrnb,
+                      gmx_wallcycle_t   wcycle,
+                      Update           *upd,
+                      gmx::Constraints *constr,
+                      bool              do_log,
+                      bool              do_ene)
 {
     if (!constr)
     {
@@ -1559,12 +1574,12 @@ update_sd_second_half(int64_t                        step,
                 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
 
                 doSDUpdateGeneral<SDUpdate::FrictionAndNoiseOnly>
-                    (*upd->sd,
+                    (*upd->sd(),
                     start_th, end_th, dt,
                     inputrec->opts.acc, inputrec->opts.nFreeze,
                     md->invmass, md->ptype,
                     md->cFREEZE, nullptr, md->cTC,
-                    state->x.rvec_array(), upd->xp.rvec_array(),
+                    state->x.rvec_array(), upd->xp()->rvec_array(),
                     state->v.rvec_array(), nullptr,
                     step, inputrec->ld_seed,
                     DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
@@ -1577,21 +1592,21 @@ update_sd_second_half(int64_t                        step,
         /* Constrain the coordinates upd->xp for half a time step */
         constr->apply(do_log, do_ene,
                       step, 1, 0.5,
-                      state->x.rvec_array(), upd->xp.rvec_array(), nullptr,
+                      state->x.rvec_array(), upd->xp()->rvec_array(), nullptr,
                       state->box,
                       state->lambda[efptBONDED], dvdlambda,
                       as_rvec_array(state->v.data()), nullptr, ConstraintVariable::Positions);
     }
 }
 
-void finish_update(const t_inputrec              *inputrec,  /* input record and box stuff     */
-                   const t_mdatoms               *md,
-                   t_state                       *state,
-                   const t_graph                 *graph,
-                   t_nrnb                        *nrnb,
-                   gmx_wallcycle_t                wcycle,
-                   gmx_update_t                  *upd,
-                   const gmx::Constraints        *constr)
+void finish_update(const t_inputrec       *inputrec, /* input record and box stuff     */
+                   const t_mdatoms        *md,
+                   t_state                *state,
+                   const t_graph          *graph,
+                   t_nrnb                 *nrnb,
+                   gmx_wallcycle_t         wcycle,
+                   Update                 *upd,
+                   const gmx::Constraints *constr)
 {
     int homenr = md->homenr;
 
@@ -1624,7 +1639,7 @@ void finish_update(const t_inputrec              *inputrec,  /* input record and
             }
             if (partialFreezeAndConstraints)
             {
-                auto xp = makeArrayRef(upd->xp).subArray(0, homenr);
+                auto xp = makeArrayRef(*upd->xp()).subArray(0, homenr);
                 auto x  = makeConstArrayRef(state->x).subArray(0, homenr);
                 for (int i = 0; i < homenr; i++)
                 {
@@ -1643,7 +1658,7 @@ void finish_update(const t_inputrec              *inputrec,  /* input record and
 
         if (graph && (graph->nnodes > 0))
         {
-            unshift_x(graph, state->box, state->x.rvec_array(), upd->xp.rvec_array());
+            unshift_x(graph, state->box, state->x.rvec_array(), upd->xp()->rvec_array());
             if (TRICLINIC(state->box))
             {
                 inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);
@@ -1655,7 +1670,7 @@ void finish_update(const t_inputrec              *inputrec,  /* input record and
         }
         else
         {
-            auto           xp = makeConstArrayRef(upd->xp).subArray(0, homenr);
+            auto           xp = makeConstArrayRef(*upd->xp()).subArray(0, homenr);
             auto           x  = makeArrayRef(state->x).subArray(0, homenr);
 #ifndef __clang_analyzer__
             int gmx_unused nth = gmx_omp_nthreads_get(emntUpdate);
@@ -1682,7 +1697,7 @@ void update_pcouple_after_coordinates(FILE             *fplog,
                                       const matrix      parrinellorahmanMu,
                                       t_state          *state,
                                       t_nrnb           *nrnb,
-                                      gmx_update_t     *upd)
+                                      Update           *upd)
 {
     int  start  = 0;
     int  homenr = md->homenr;
@@ -1764,10 +1779,10 @@ void update_pcouple_after_coordinates(FILE             *fplog,
             break;
     }
 
-    if (upd->deform)
+    if (upd->deform())
     {
         auto localX = makeArrayRef(state->x).subArray(start, homenr);
-        upd->deform->apply(localX, state->box, step);
+        upd->deform()->apply(localX, state->box, step);
     }
 }
 
@@ -1779,7 +1794,7 @@ void update_coords(int64_t                             step,
                    const t_fcdata                     *fcd,
                    const gmx_ekindata_t               *ekind,
                    const matrix                        M,
-                   gmx_update_t                       *upd,
+                   Update                             *upd,
                    int                                 UpdatePart,
                    const t_commrec                    *cr, /* these shouldn't be here -- need to think about it */
                    const gmx::Constraints             *constr)
@@ -1820,7 +1835,7 @@ void update_coords(int64_t                             step,
             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
 
             const rvec *x_rvec  = state->x.rvec_array();
-            rvec       *xp_rvec = upd->xp.rvec_array();
+            rvec       *xp_rvec = upd->xp()->rvec_array();
             rvec       *v_rvec  = state->v.rvec_array();
             const rvec *f_rvec  = as_rvec_array(f.unpaddedArrayRef().data());
 
@@ -1837,7 +1852,7 @@ void update_coords(int64_t                             step,
                     {
                         // With constraints, the SD update is done in 2 parts
                         doSDUpdateGeneral<SDUpdate::ForcesOnly>
-                            (*upd->sd,
+                            (*upd->sd(),
                             start_th, end_th, dt,
                             inputrec->opts.acc, inputrec->opts.nFreeze,
                             md->invmass, md->ptype,
@@ -1848,7 +1863,7 @@ void update_coords(int64_t                             step,
                     else
                     {
                         doSDUpdateGeneral<SDUpdate::Combined>
-                            (*upd->sd,
+                            (*upd->sd(),
                             start_th, end_th, dt,
                             inputrec->opts.acc, inputrec->opts.nFreeze,
                             md->invmass, md->ptype,
@@ -1864,7 +1879,7 @@ void update_coords(int64_t                             step,
                                  md->cFREEZE, md->cTC,
                                  x_rvec, xp_rvec, v_rvec, f_rvec,
                                  inputrec->bd_fric,
-                                 upd->sd->bd_rf.data(),
+                                 upd->sd()->bd_rf.data(),
                                  step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
                     break;
                 case (eiVV):
@@ -1909,7 +1924,7 @@ void update_coords(int64_t                             step,
 extern gmx_bool update_randomize_velocities(const t_inputrec *ir, int64_t step, const t_commrec *cr,
                                             const t_mdatoms *md,
                                             gmx::ArrayRef<gmx::RVec> v,
-                                            const gmx_update_t *upd,
+                                            const Update *upd,
                                             const gmx::Constraints *constr)
 {
 
@@ -1931,8 +1946,8 @@ extern gmx_bool update_randomize_velocities(const t_inputrec *ir, int64_t step,
     if ((ir->etc == etcANDERSEN) || do_per_step(step, roundToInt(1.0/rate)))
     {
         andersen_tcoupl(ir, step, cr, md, v, rate,
-                        upd->sd->randomize_group,
-                        upd->sd->boltzfac);
+                        upd->sd()->randomize_group,
+                        upd->sd()->boltzfac);
         return TRUE;
     }
     return FALSE;
index 366d0621f687d377914f50f74b9e170fcbe3b911..4d24ee32f54ce451b620b0b51c7a9bb886198d4f 100644 (file)
@@ -43,6 +43,7 @@
 #include "gromacs/timing/wallcycle.h"
 #include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/classhelpers.h"
 #include "gromacs/utility/real.h"
 
 class ekinstate_t;
@@ -58,33 +59,48 @@ struct t_nrnb;
 class t_state;
 
 /* Abstract type for update */
-struct gmx_update_t;
+struct gmx_stochd_t;
 
 namespace gmx
 {
 class BoxDeformation;
 class Constraints;
-}
 
-/* Initialize the stochastic dynamics struct */
-gmx_update_t *init_update(const t_inputrec    *ir,
-                          gmx::BoxDeformation *deform);
+
+/*! \libinternal
+ * \brief Contains data for update phase */
+class Update
+{
+    public:
+        //! Constructor
+        Update(const t_inputrec    *ir, BoxDeformation *boxDeformation);
+        ~Update();
+        // TODO Get rid of getters when more free functions are incorporated as member methods
+        //! Returns handle to stochd_t struct
+        gmx_stochd_t * sd() const;
+        //! Returns pointer to PaddedVector xp
+        PaddedVector<gmx::RVec> * xp();
+        //! Returns handle to box deformation class
+        BoxDeformation * deform() const;
+        //! Resizes xp
+        void setNumAtoms(int nAtoms);
+    private:
+        //! Implementation type.
+        class Impl;
+        //! Implementation object.
+        PrivateImplPointer<Impl> impl_;
+};
+
+}; // namespace gmx
 
 /* Update pre-computed constants that depend on the reference
  * temperature for coupling.
  *
  * This could change e.g. in simulated annealing. */
-void update_temperature_constants(gmx_update_t *upd, const t_inputrec *ir);
+void update_temperature_constants(gmx_stochd_t *sd, const t_inputrec *ir);
 
 /* Update the size of per-atom arrays (e.g. after DD re-partitioning,
    which might increase the number of home atoms). */
-void update_realloc(gmx_update_t *upd, int natoms);
-
-/* Store the box at step step
- * as a reference state for simulations with box deformation.
- */
-void set_deform_reference_box(gmx_update_t *upd,
-                              int64_t step, matrix box);
 
 void update_tcouple(int64_t           step,
                     const t_inputrec *inputrec,
@@ -118,27 +134,27 @@ void update_pcouple_after_coordinates(FILE             *fplog,
                                       const matrix      parrinellorahmanMu,
                                       t_state          *state,
                                       t_nrnb           *nrnb,
-                                      gmx_update_t     *upd);
-
-void update_coords(int64_t                              step,
-                   const t_inputrec                    *inputrec, /* input record and box stuff        */
-                   const t_mdatoms                     *md,
-                   t_state                             *state,
-                   gmx::ArrayRefWithPadding<gmx::RVec>  f, /* forces on home particles */
-                   const t_fcdata                      *fcd,
-                   const gmx_ekindata_t                *ekind,
-                   const matrix                         M,
-                   gmx_update_t                        *upd,
-                   int                                  bUpdatePart,
-                   const t_commrec                     *cr, /* these shouldn't be here -- need to think about it */
-                   const gmx::Constraints              *constr);
+                                      gmx::Update      *upd);
+
+void update_coords(int64_t                             step,
+                   const t_inputrec                   *inputrec, /* input record and box stuff */
+                   const t_mdatoms                    *md,
+                   t_state                            *state,
+                   gmx::ArrayRefWithPadding<gmx::RVec> f, /* forces on home particles */
+                   const t_fcdata                     *fcd,
+                   const gmx_ekindata_t               *ekind,
+                   const matrix                        M,
+                   gmx::Update                        *upd,
+                   int                                 bUpdatePart,
+                   const t_commrec                    *cr, /* these shouldn't be here -- need to think about it */
+                   const gmx::Constraints             *constr);
 
 /* Return TRUE if OK, FALSE in case of Shake Error */
 
 extern gmx_bool update_randomize_velocities(const t_inputrec *ir, int64_t step, const t_commrec *cr,
                                             const t_mdatoms *md,
                                             gmx::ArrayRef<gmx::RVec> v,
-                                            const gmx_update_t *upd,
+                                            const gmx::Update *upd,
                                             const gmx::Constraints *constr);
 
 void constrain_velocities(int64_t                        step,
@@ -150,37 +166,37 @@ void constrain_velocities(int64_t                        step,
                           bool                           do_log,
                           bool                           do_ene);
 
-void constrain_coordinates(int64_t                        step,
-                           real                          *dvdlambda, /* the contribution to be added to the bonded interactions */
-                           t_state                       *state,
-                           tensor                         vir_part,
-                           gmx_update_t                  *upd,
-                           gmx::Constraints              *constr,
-                           gmx_bool                       bCalcVir,
-                           bool                           do_log,
-                           bool                           do_ene);
-
-void update_sd_second_half(int64_t                        step,
-                           real                          *dvdlambda, /* the contribution to be added to the bonded interactions */
-                           const t_inputrec              *inputrec,  /* input record and box stuff */
-                           const t_mdatoms               *md,
-                           t_state                       *state,
-                           const t_commrec               *cr,
-                           t_nrnb                        *nrnb,
-                           gmx_wallcycle_t                wcycle,
-                           gmx_update_t                  *upd,
-                           gmx::Constraints              *constr,
-                           bool                           do_log,
-                           bool                           do_ene);
-
-void finish_update(const t_inputrec              *inputrec,
-                   const t_mdatoms               *md,
-                   t_state                       *state,
-                   const t_graph                 *graph,
-                   t_nrnb                        *nrnb,
-                   gmx_wallcycle_t                wcycle,
-                   gmx_update_t                  *upd,
-                   const gmx::Constraints        *constr);
+void constrain_coordinates(int64_t           step,
+                           real             *dvdlambda, /* the contribution to be added to the bonded interactions */
+                           t_state          *state,
+                           tensor            vir_part,
+                           gmx::Update      *upd,
+                           gmx::Constraints *constr,
+                           gmx_bool          bCalcVir,
+                           bool              do_log,
+                           bool              do_ene);
+
+void update_sd_second_half(int64_t           step,
+                           real             *dvdlambda, /* the contribution to be added to the bonded interactions */
+                           const t_inputrec *inputrec,  /* input record and box stuff */
+                           const t_mdatoms  *md,
+                           t_state          *state,
+                           const t_commrec  *cr,
+                           t_nrnb           *nrnb,
+                           gmx_wallcycle_t   wcycle,
+                           gmx::Update      *upd,
+                           gmx::Constraints *constr,
+                           bool              do_log,
+                           bool              do_ene);
+
+void finish_update(const t_inputrec       *inputrec,
+                   const t_mdatoms        *md,
+                   t_state                *state,
+                   const t_graph          *graph,
+                   t_nrnb                 *nrnb,
+                   gmx_wallcycle_t         wcycle,
+                   gmx::Update            *upd,
+                   const gmx::Constraints *constr);
 
 /* Return TRUE if OK, FALSE in case of Shake Error */
 
@@ -248,7 +264,7 @@ void rescale_velocities(const gmx_ekindata_t *ekind, const t_mdatoms *mdatoms,
 /* Rescale the velocities with the scaling factor in ekind */
 
 // TODO: This is the only function in update.h altering the inputrec
-void update_annealing_target_temp(t_inputrec *ir, real t, gmx_update_t *upd);
+void update_annealing_target_temp(t_inputrec *ir, real t, gmx::Update *upd);
 /* Set reference temp for simulated annealing at time t*/
 
 real calc_temp(real ekin, real nrdf);
index 1aed281feb0abbea5f6bf05a5af14d7cfa93bbee..d4bd7daa17fd171d343d88ae8f3a062fdd1e9b68 100644 (file)
@@ -173,7 +173,6 @@ void gmx::Integrator::do_md()
     gmx_enerdata_t         *enerd;
     PaddedVector<gmx::RVec> f {};
     gmx_global_stat_t       gstat;
-    gmx_update_t           *upd   = nullptr;
     t_graph                *graph = nullptr;
     gmx_groups_t           *groups;
     gmx_shellfc_t          *shellfc;
@@ -248,10 +247,11 @@ void gmx::Integrator::do_md()
                         oenv, mdrunOptions.continuationOptions.appendFiles);
     }
 
+    Update upd(ir, deform);
     /* Initial values */
     init_md(fplog, cr, outputProvider, ir, oenv, mdrunOptions,
             &t, &t0, state_global, lam0,
-            nrnb, top_global, &upd, deform,
+            nrnb, top_global, &upd,
             nfile, fnm, &outf, &mdebin,
             force_vir, shake_vir, total_vir, pres, mu_tot, &bSimAnn, wcycle);
 
@@ -308,7 +308,7 @@ void gmx::Integrator::do_md()
                             vsite, constr,
                             nrnb, nullptr, FALSE);
         shouldCheckNumberOfBondedInteractions = true;
-        update_realloc(upd, state->natoms);
+        upd.setNumAtoms(state->natoms);
     }
     else
     {
@@ -321,7 +321,7 @@ void gmx::Integrator::do_md()
         mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr,
                                   &graph, mdAtoms, constr, vsite, shellfc);
 
-        update_realloc(upd, state->natoms);
+        upd.setNumAtoms(state->natoms);
     }
 
     auto mdatoms = mdAtoms->mdatoms();
@@ -723,7 +723,7 @@ void gmx::Integrator::do_md()
 
         if (bSimAnn)
         {
-            update_annealing_target_temp(ir, t, upd);
+            update_annealing_target_temp(ir, t, &upd);
         }
 
         /* Stop Center of Mass motion */
@@ -771,7 +771,7 @@ void gmx::Integrator::do_md()
                                     nrnb, wcycle,
                                     do_verbose && !bPMETunePrinting);
                 shouldCheckNumberOfBondedInteractions = true;
-                update_realloc(upd, state->natoms);
+                upd.setNumAtoms(state->natoms);
             }
         }
 
@@ -909,7 +909,7 @@ void gmx::Integrator::do_md()
             }
 
             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
-                          ekind, M, upd, etrtVELOCITY1,
+                          ekind, M, &upd, etrtVELOCITY1,
                           cr, constr);
 
             wallcycle_stop(wcycle, ewcUPDATE);
@@ -1089,7 +1089,7 @@ void gmx::Integrator::do_md()
         if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
         {
             gmx_bool bIfRandomize;
-            bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, upd, constr);
+            bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
             /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
             if (constr && bIfRandomize)
             {
@@ -1132,7 +1132,7 @@ void gmx::Integrator::do_md()
         {
             /* velocity half-step update */
             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
-                          ekind, M, upd, etrtVELOCITY2,
+                          ekind, M, &upd, etrtVELOCITY2,
                           cr, constr);
         }
 
@@ -1153,18 +1153,18 @@ void gmx::Integrator::do_md()
         }
 
         update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
-                      ekind, M, upd, etrtPOSITION, cr, constr);
+                      ekind, M, &upd, etrtPOSITION, cr, constr);
         wallcycle_stop(wcycle, ewcUPDATE);
 
         constrain_coordinates(step, &dvdl_constr, state,
                               shake_vir,
-                              upd, constr,
+                              &upd, constr,
                               bCalcVir, do_log, do_ene);
         update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
-                              cr, nrnb, wcycle, upd, constr, do_log, do_ene);
+                              cr, nrnb, wcycle, &upd, constr, do_log, do_ene);
         finish_update(ir, mdatoms,
                       state, graph,
-                      nrnb, wcycle, upd, constr);
+                      nrnb, wcycle, &upd, constr);
 
         if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
         {
@@ -1187,7 +1187,7 @@ void gmx::Integrator::do_md()
             copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
 
             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
-                          ekind, M, upd, etrtPOSITION, cr, constr);
+                          ekind, M, &upd, etrtPOSITION, cr, constr);
             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? */
@@ -1197,7 +1197,7 @@ void gmx::Integrator::do_md()
              * For now, will call without actually constraining, constr=NULL*/
             finish_update(ir, mdatoms,
                           state, graph,
-                          nrnb, wcycle, upd, nullptr);
+                          nrnb, wcycle, &upd, nullptr);
         }
         if (EI_VV(ir->eI))
         {
@@ -1298,7 +1298,7 @@ void gmx::Integrator::do_md()
         update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
                                          pres, force_vir, shake_vir,
                                          parrinellorahmanMu,
-                                         state, nrnb, upd);
+                                         state, nrnb, &upd);
 
         /* ################# END UPDATE STEP 2 ################# */
         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
@@ -1433,7 +1433,7 @@ void gmx::Integrator::do_md()
                                 vsite, constr,
                                 nrnb, wcycle, FALSE);
             shouldCheckNumberOfBondedInteractions = true;
-            update_realloc(upd, state->natoms);
+            upd.setNumAtoms(state->natoms);
         }
 
         bFirstStep             = FALSE;