Remove the dependence on Constraints object on t_mdatoms

This removes t_mdatoms structure from all function signatures in
Constraints object. Only those data neede for Constraints is
now passed upon setup and stored locally.

diff --git a/src/gromacs/domdec/mdsetup.cpp b/src/gromacs/domdec/mdsetup.cpp
index f8d9172c6b5a7894b261199d6e84558d0ae27ff0..d18168947e04433a0792d86b48d59b391ea37781 100644 (file)
--- a/src/gromacs/domdec/mdsetup.cpp
+++ b/src/gromacs/domdec/mdsetup.cpp
@@ -144,7 +144,8 @@ void mdAlgorithmsSetupAtomData(const t_commrec*        cr,
 
     if (constr)
     {
-        constr->setConstraints(top, *mdatoms);
+        constr->setConstraints(top, mdatoms->nr, mdatoms->homenr, mdatoms->massT, mdatoms->invmass,
+                               mdatoms->nMassPerturbed != 0, mdatoms->lambda, mdatoms->cFREEZE);
     }
 }
 

diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp
index cbc6ba9e3620f5d6f3aea85dfb2b9bec9c680dcb..f90107f9e9cec8bb38f9390bc9f8face15136c59 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -68,7 +68,6 @@
 #include "gromacs/mdtypes/commrec.h"
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/mdtypes/state.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/pulling/pull.h"
@@ -102,7 +101,6 @@ public:
          const t_inputrec&     ir_p,
          pull_t*               pull_work,
          FILE*                 log_p,
-         const t_mdatoms&      md_p,
          const t_commrec*      cr_p,
          const gmx_multisim_t* ms,
          t_nrnb*               nrnb,
@@ -111,7 +109,14 @@ public:
          int                   numConstraints,
          int                   numSettles);
     ~Impl();
-    void setConstraints(gmx_localtop_t* top, const t_mdatoms& md);
+    void setConstraints(gmx_localtop_t* top,
+                        int             numAtoms,
+                        int             numHomeAtoms,
+                        real*           masses,
+                        real*           inverseMasses,
+                        bool            hasMassPerturbedAtoms,
+                        real            lambda,
+                        unsigned short* cFREEZE);
     bool apply(bool                      bLog,
                bool                      bEner,
                int64_t                   step,
@@ -158,8 +163,20 @@ public:
     const gmx_mtop_t& mtop;
     //! Parameters for the interactions in this domain.
     const InteractionDefinitions* idef = nullptr;
-    //! Data about atoms in this domain.
-    const t_mdatoms& md;
+    //! Total number of atoms.
+    int numAtoms_ = 0;
+    //! Number of local atoms.
+    int numHomeAtoms_ = 0;
+    //! Atoms masses.
+    real* masses_;
+    //! Inverse masses.
+    real* inverseMasses_;
+    //! If there are atoms with perturbed mass (for FEP).
+    bool hasMassPerturbedAtoms_;
+    //! FEP lambda value.
+    real lambda_;
+    //! Freeze groups data
+    unsigned short* cFREEZE_;
     //! Whether we need to do pbc for handling bonds.
     bool pbcHandlingRequired_ = false;
 
@@ -356,7 +373,7 @@ bool Constraints::Impl::apply(bool                      bLog,
                               ConstraintVariable        econq)
 {
     bool   bOK, bDump;
-    int    start, homenr;
+    int    start;
     tensor vir_r_m_dr;
     real   scaled_delta_t;
     real   invdt, vir_fac = 0, t;
@@ -378,8 +395,7 @@ bool Constraints::Impl::apply(bool                      bLog,
     bOK   = TRUE;
     bDump = FALSE;
 
-    start  = 0;
-    homenr = md.homenr;
+    start = 0;
 
     scaled_delta_t = step_scaling * ir.delta_t;
 
@@ -457,8 +473,8 @@ bool Constraints::Impl::apply(bool                      bLog,
 
     if (lincsd != nullptr)
     {
-        bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md.invmass, cr, ms, x, xprime,
-                              min_proj, box, pbc_null, md.nMassPerturbed != 0, lambda, dvdlambda,
+        bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, inverseMasses_, cr, ms, x, xprime,
+                              min_proj, box, pbc_null, hasMassPerturbedAtoms_, lambda, dvdlambda,
                               invdt, v.unpaddedArrayRef(), vir != nullptr, vir_r_m_dr, econq, nrnb,
                               maxwarn, &warncount_lincs);
         if (!bOK && maxwarn < INT_MAX)
@@ -474,7 +490,7 @@ bool Constraints::Impl::apply(bool                      bLog,
 
     if (shaked != nullptr)
     {
-        bOK = constrain_shake(log, shaked.get(), md.invmass, *idef, ir, x.unpaddedArrayRef(),
+        bOK = constrain_shake(log, shaked.get(), inverseMasses_, *idef, ir, x.unpaddedArrayRef(),
                               xprime.unpaddedArrayRef(), min_proj, pbc_null, nrnb, lambda,
                               dvdlambda, invdt, v.unpaddedArrayRef(), vir != nullptr, vir_r_m_dr,
                               maxwarn < INT_MAX, econq);
@@ -540,7 +556,7 @@ bool Constraints::Impl::apply(bool                      bLog,
                         }
                         else
                         {
-                            calcvir_atom_end = md.homenr;
+                            calcvir_atom_end = numHomeAtoms_;
                         }
 
                         if (th > 0)
@@ -646,7 +662,7 @@ bool Constraints::Impl::apply(bool                      bLog,
 
     if (bDump)
     {
-        dump_confs(log, step, mtop, start, homenr, cr, x.unpaddedArrayRef(),
+        dump_confs(log, step, mtop, start, numHomeAtoms_, cr, x.unpaddedArrayRef(),
                    xprime.unpaddedArrayRef(), box);
     }
 
@@ -663,7 +679,7 @@ bool Constraints::Impl::apply(bool                      bLog,
                 t = ir.init_t;
             }
             set_pbc(&pbc, ir.pbcType, box);
-            pull_constraint(pull_work, md.massT, &pbc, cr, ir.delta_t, t,
+            pull_constraint(pull_work, masses_, &pbc, cr, ir.delta_t, t,
                             as_rvec_array(x.unpaddedArrayRef().data()),
                             as_rvec_array(xprime.unpaddedArrayRef().data()),
                             as_rvec_array(v.unpaddedArrayRef().data()), *vir);
@@ -677,7 +693,7 @@ bool Constraints::Impl::apply(bool                      bLog,
     }
     wallcycle_stop(wcycle, ewcCONSTR);
 
-    if (!v.empty() && md.cFREEZE)
+    if (!v.empty() && cFREEZE_)
     {
         /* Set the velocities of frozen dimensions to zero */
         ArrayRef<RVec> vRef = v.unpaddedArrayRef();
@@ -685,9 +701,9 @@ bool Constraints::Impl::apply(bool                      bLog,
         int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
 
 #pragma omp parallel for num_threads(numThreads) schedule(static)
-        for (int i = 0; i < md.homenr; i++)
+        for (int i = 0; i < numHomeAtoms_; i++)
         {
-            int freezeGroup = md.cFREEZE[i];
+            int freezeGroup = cFREEZE_[i];
 
             for (int d = 0; d < DIM; d++)
             {
@@ -874,8 +890,23 @@ static std::vector<int> make_at2settle(int natoms, const InteractionList& ilist)
     return at2s;
 }
 
-void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
+void Constraints::Impl::setConstraints(gmx_localtop_t* top,
+                                       int             numAtoms,
+                                       int             numHomeAtoms,
+                                       real*           masses,
+                                       real*           inverseMasses,
+                                       const bool      hasMassPerturbedAtoms,
+                                       const real      lambda,
+                                       unsigned short* cFREEZE)
 {
+    numAtoms_              = numAtoms;
+    numHomeAtoms_          = numHomeAtoms;
+    masses_                = masses;
+    inverseMasses_         = inverseMasses;
+    hasMassPerturbedAtoms_ = hasMassPerturbedAtoms;
+    lambda_                = lambda;
+    cFREEZE_               = cFREEZE;
+
     idef = &top->idef;
 
     if (ncon_tot > 0)
@@ -885,7 +916,7 @@ void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
          */
         if (ir.eConstrAlg == econtLINCS)
         {
-            set_lincs(*idef, md.nr, md.invmass, md.lambda, EI_DYNAMICS(ir.eI), cr, lincsd);
+            set_lincs(*idef, numAtoms_, inverseMasses_, lambda_, EI_DYNAMICS(ir.eI), cr, lincsd);
         }
         if (ir.eConstrAlg == econtSHAKE)
         {
@@ -897,14 +928,14 @@ void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
             }
             else
             {
-                make_shake_sblock_serial(shaked.get(), &top->idef, md.nr);
+                make_shake_sblock_serial(shaked.get(), &top->idef, numAtoms_);
             }
         }
     }
 
     if (settled)
     {
-        settle_set_constraints(settled, idef->il[F_SETTLE], md.homenr, md.massT, md.invmass);
+        settle_set_constraints(settled, idef->il[F_SETTLE], numHomeAtoms_, masses_, inverseMasses_);
     }
 
     /* Make a selection of the local atoms for essential dynamics */
@@ -914,9 +945,17 @@ void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
     }
 }
 
-void Constraints::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
+void Constraints::setConstraints(gmx_localtop_t* top,
+                                 const int       numAtoms,
+                                 const int       numHomeAtoms,
+                                 real*           masses,
+                                 real*           inverseMasses,
+                                 const bool      hasMassPerturbedAtoms,
+                                 const real      lambda,
+                                 unsigned short* cFREEZE)
 {
-    impl_->setConstraints(top, md);
+    impl_->setConstraints(top, numAtoms, numHomeAtoms, masses, inverseMasses, hasMassPerturbedAtoms,
+                          lambda, cFREEZE);
 }
 
 /*! \brief Makes a per-moleculetype container of mappings from atom
@@ -939,7 +978,6 @@ Constraints::Constraints(const gmx_mtop_t&     mtop,
                          const t_inputrec&     ir,
                          pull_t*               pull_work,
                          FILE*                 log,
-                         const t_mdatoms&      md,
                          const t_commrec*      cr,
                          const gmx_multisim_t* ms,
                          t_nrnb*               nrnb,
@@ -947,7 +985,7 @@ Constraints::Constraints(const gmx_mtop_t&     mtop,
                          bool                  pbcHandlingRequired,
                          int                   numConstraints,
                          int                   numSettles) :
-    impl_(new Impl(mtop, ir, pull_work, log, md, cr, ms, nrnb, wcycle, pbcHandlingRequired, numConstraints, numSettles))
+    impl_(new Impl(mtop, ir, pull_work, log, cr, ms, nrnb, wcycle, pbcHandlingRequired, numConstraints, numSettles))
 {
 }
 
@@ -955,7 +993,6 @@ Constraints::Impl::Impl(const gmx_mtop_t&     mtop_p,
                         const t_inputrec&     ir_p,
                         pull_t*               pull_work,
                         FILE*                 log_p,
-                        const t_mdatoms&      md_p,
                         const t_commrec*      cr_p,
                         const gmx_multisim_t* ms_p,
                         t_nrnb*               nrnb_p,
@@ -965,7 +1002,6 @@ Constraints::Impl::Impl(const gmx_mtop_t&     mtop_p,
                         int                   numSettles) :
     ncon_tot(numConstraints),
     mtop(mtop_p),
-    md(md_p),
     pbcHandlingRequired_(pbcHandlingRequired),
     log(log_p),
     cr(cr_p),
@@ -1125,8 +1161,8 @@ ArrayRef<const std::vector<int>> Constraints::atom2settle_moltype() const
 void do_constrain_first(FILE*                     fplog,
                         gmx::Constraints*         constr,
                         const t_inputrec*         ir,
-                        const t_mdatoms*          md,
-                        int                       natoms,
+                        int                       numAtoms,
+                        int                       numHomeAtoms,
                         ArrayRefWithPadding<RVec> x,
                         ArrayRefWithPadding<RVec> v,
                         const matrix              box,
@@ -1137,14 +1173,14 @@ void do_constrain_first(FILE*                     fplog,
     real    dt = ir->delta_t;
     real    dvdl_dum;
 
-    PaddedVector<RVec> savex(natoms);
+    PaddedVector<RVec> savex(numAtoms);
 
     start = 0;
-    end   = md->homenr;
+    end   = numHomeAtoms;
 
     if (debug)
     {
-        fprintf(debug, "vcm: start=%d, homenr=%d, end=%d\n", start, md->homenr, end);
+        fprintf(debug, "vcm: start=%d, homenr=%d, end=%d\n", start, numHomeAtoms, end);
     }
     /* Do a first constrain to reset particles... */
     step = ir->init_step;

diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h
index 680a941acbee8edc4df837da534cdcc5aa4e2908..2d45d410d7b8ddcb91113e36bb1a8bb5f018fc89 100644 (file)
--- a/src/gromacs/mdlib/constr.h
+++ b/src/gromacs/mdlib/constr.h
@@ -67,7 +67,6 @@ struct pull_t;
 struct t_commrec;
 struct t_ilist;
 struct t_inputrec;
-struct t_mdatoms;
 struct t_nrnb;
 struct t_pbc;
 class t_state;
@@ -107,7 +106,6 @@ private:
                 const t_inputrec&     ir,
                 pull_t*               pull_work,
                 FILE*                 log,
-                const t_mdatoms&      md,
                 const t_commrec*      cr,
                 const gmx_multisim_t* ms,
                 t_nrnb*               nrnb,
@@ -134,7 +132,14 @@ public:
      *
      * \todo Make this a callback that is called automatically
      * once a new domain has been made. */
-    void setConstraints(gmx_localtop_t* top, const t_mdatoms& md);
+    void setConstraints(gmx_localtop_t* top,
+                        int             numAtoms,
+                        int             numHomeAtoms,
+                        real*           masses,
+                        real*           inverseMasses,
+                        bool            hasMassPerturbedAtoms,
+                        real            lambda,
+                        unsigned short* cFREEZE);
 
     /*! \brief Applies constraints to coordinates.
      *
@@ -304,8 +309,8 @@ bool inter_charge_group_settles(const gmx_mtop_t& mtop);
 void do_constrain_first(FILE*                     log,
                         gmx::Constraints*         constr,
                         const t_inputrec*         inputrec,
-                        const t_mdatoms*          md,
-                        int                       natoms,
+                        int                       numAtoms,
+                        int                       numHomeAtoms,
                         ArrayRefWithPadding<RVec> x,
                         ArrayRefWithPadding<RVec> v,
                         const matrix              box,

diff --git a/src/gromacs/mdlib/tests/energyoutput.cpp b/src/gromacs/mdlib/tests/energyoutput.cpp
index 06f2bf73280810331d6841ebc9402b34ffe1f74e..e22cb28991f775702f151aac0fecc87308cd53ed 100644 (file)
--- a/src/gromacs/mdlib/tests/energyoutput.cpp
+++ b/src/gromacs/mdlib/tests/energyoutput.cpp
@@ -145,8 +145,6 @@ public:
     t_inputrec inputrec_;
     //! Topology
     gmx_mtop_t mtop_;
-    //! MD atoms
-    t_mdatoms mdatoms_;
     //! Simulation time
     double time_;
     //! Total mass
@@ -378,8 +376,8 @@ public:
         // TODO EnergyOutput should not take Constraints object
         // TODO This object will always return zero as RMSD value.
         //      It is more relevant to have non-zero value for testing.
-        constraints_ = makeConstraints(mtop_, inputrec_, nullptr, false, nullptr, mdatoms_, &cr_,
-                                       nullptr, nullptr, nullptr, false);
+        constraints_ = makeConstraints(mtop_, inputrec_, nullptr, false, nullptr, &cr_, nullptr,
+                                       nullptr, nullptr, false);
 
         // No position/orientation restraints
         fcd_.disres.npair = 0;

diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp
index 06bb8a9211ca777baead3cf21ed89cc0864bb36b..190a77ad3dc347f85d4037e6e8c01b7061e785ab 100644 (file)
--- a/src/gromacs/mdrun/md.cpp
+++ b/src/gromacs/mdrun/md.cpp
@@ -501,8 +501,9 @@ void gmx::LegacySimulator::do_md()
         if (constr)
         {
             /* Constrain the initial coordinates and velocities */
-            do_constrain_first(fplog, constr, ir, mdatoms, state->natoms, state->x.arrayRefWithPadding(),
-                               state->v.arrayRefWithPadding(), state->box, state->lambda[efptBONDED]);
+            do_constrain_first(fplog, constr, ir, mdatoms->nr, mdatoms->homenr,
+                               state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(),
+                               state->box, state->lambda[efptBONDED]);
         }
         if (vsite)
         {

diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp
index 33f3a23b0941a7c777ca9d3cdadf8c064a27c95d..12bf53fbe4ad946136be3a466eb32fc4863d8490 100644 (file)
--- a/src/gromacs/mdrun/runner.cpp
+++ b/src/gromacs/mdrun/runner.cpp
@@ -1579,8 +1579,8 @@ int Mdrunner::mdrunner()
         }
 
         /* Let makeConstraints know whether we have essential dynamics constraints. */
-        auto constr = makeConstraints(mtop, *inputrec, pull_work, doEssentialDynamics, fplog,
-                                      *mdAtoms->mdatoms(), cr, ms, &nrnb, wcycle, fr->bMolPBC);
+        auto constr = makeConstraints(mtop, *inputrec, pull_work, doEssentialDynamics, fplog, cr,
+                                      ms, &nrnb, wcycle, fr->bMolPBC);
 
         /* Energy terms and groups */
         gmx_enerdata_t enerd(mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),

diff --git a/src/gromacs/modularsimulator/constraintelement.cpp b/src/gromacs/modularsimulator/constraintelement.cpp
index 23145242adbda427650e25ebc119b4e05659cb28..9a202f57837cc52fe93dc1a60500948f8b4d17eb 100644 (file)
--- a/src/gromacs/modularsimulator/constraintelement.cpp
+++ b/src/gromacs/modularsimulator/constraintelement.cpp
@@ -90,9 +90,10 @@ void ConstraintsElement<variable>::elementSetup()
                                           ? freeEnergyPerturbationElement_->constLambdaView()[efptBONDED]
                                           : 0;
         // Constrain the initial coordinates and velocities
-        do_constrain_first(fplog_, constr_, inputrec_, mdAtoms_, statePropagatorData_->localNumAtoms(),
-                           statePropagatorData_->positionsView(), statePropagatorData_->velocitiesView(),
-                           statePropagatorData_->box(), lambdaBonded);
+        do_constrain_first(
+                fplog_, constr_, inputrec_, statePropagatorData_->totalNumAtoms(),
+                statePropagatorData_->localNumAtoms(), statePropagatorData_->positionsView(),
+                statePropagatorData_->velocitiesView(), statePropagatorData_->box(), lambdaBonded);
 
         if (isMasterRank_)
         {

diff --git a/src/gromacs/modularsimulator/statepropagatordata.cpp b/src/gromacs/modularsimulator/statepropagatordata.cpp
index 8774aeaf0b244adc368fdd904fc47c4b418e18d6..acb6249455218124142ecf912b1ee1d5a79db412 100644 (file)
--- a/src/gromacs/modularsimulator/statepropagatordata.cpp
+++ b/src/gromacs/modularsimulator/statepropagatordata.cpp
@@ -230,6 +230,11 @@ int StatePropagatorData::localNumAtoms()
     return localNAtoms_;
 }
 
+int StatePropagatorData::totalNumAtoms()
+{
+    return totalNumAtoms_;
+}
+
 std::unique_ptr<t_state> StatePropagatorData::localState()
 {
     auto state   = std::make_unique<t_state>();

diff --git a/src/gromacs/modularsimulator/statepropagatordata.h b/src/gromacs/modularsimulator/statepropagatordata.h
index 3161299b9061244b0fb20e7407db6b5c6c3e68fa..7b72d4dd8e9855c451b6cc0f1929bc2c60995d02 100644 (file)
--- a/src/gromacs/modularsimulator/statepropagatordata.h
+++ b/src/gromacs/modularsimulator/statepropagatordata.h
@@ -151,6 +151,8 @@ public:
     const rvec* constPreviousBox() const;
     //! Get the local number of atoms
     int localNumAtoms();
+    //! Get the total number of atoms
+    int totalNumAtoms();
 
     /*! \brief Register run function for step / time
      *