Remove dependence of constraints on t_mdatoms
authorArtem Zhmurov <zhmurov@gmail.com>
Fri, 22 May 2020 07:40:50 +0000 (07:40 +0000)
committerMark Abraham <mark.j.abraham@gmail.com>
Fri, 22 May 2020 07:40:50 +0000 (07:40 +0000)
This changes the signatures of the constrants set/compute functions
by removing the t_mdatoms from the set of arguments and replacing
it with the parts that are actually used.

Refs. #3535

12 files changed:
src/gromacs/mdlib/constr.cpp
src/gromacs/mdlib/lincs.cpp
src/gromacs/mdlib/lincs.h
src/gromacs/mdlib/lincs_gpu.cu
src/gromacs/mdlib/lincs_gpu.cuh
src/gromacs/mdlib/shake.cpp
src/gromacs/mdlib/shake.h
src/gromacs/mdlib/tests/constrtestdata.cpp
src/gromacs/mdlib/tests/constrtestdata.h
src/gromacs/mdlib/tests/constrtestrunners.cpp
src/gromacs/mdlib/tests/constrtestrunners.cu
src/gromacs/mdlib/update_constrain_gpu_impl.cu

index 16b5c6eba070c08518a1a1f34ac34631597fe9a4..a158fd6c050ae4206f6e759e4e7dea54387b64f4 100644 (file)
@@ -457,9 +457,10 @@ bool Constraints::Impl::apply(bool                      bLog,
 
     if (lincsd != nullptr)
     {
-        bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md, cr, ms, x, xprime, min_proj, box,
-                              pbc_null, lambda, dvdlambda, invdt, v.unpaddedArrayRef(),
-                              vir != nullptr, vir_r_m_dr, econq, nrnb, maxwarn, &warncount_lincs);
+        bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md.invmass, cr, ms, x, xprime,
+                              min_proj, box, pbc_null, md.nMassPerturbed != 0, lambda, dvdlambda,
+                              invdt, v.unpaddedArrayRef(), vir != nullptr, vir_r_m_dr, econq, nrnb,
+                              maxwarn, &warncount_lincs);
         if (!bOK && maxwarn < INT_MAX)
         {
             if (log != nullptr)
@@ -699,7 +700,7 @@ bool Constraints::Impl::apply(bool                      bLog,
     }
 
     return bOK;
-}
+} // namespace gmx
 
 ArrayRef<real> Constraints::rmsdData() const
 {
@@ -884,7 +885,7 @@ void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
          */
         if (ir.eConstrAlg == econtLINCS)
         {
-            set_lincs(*idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
+            set_lincs(*idef, md.nr, md.invmass, md.lambda, EI_DYNAMICS(ir.eI), cr, lincsd);
         }
         if (ir.eConstrAlg == econtSHAKE)
         {
@@ -896,7 +897,7 @@ void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
             }
             else
             {
-                make_shake_sblock_serial(shaked.get(), &top->idef, md);
+                make_shake_sblock_serial(shaked.get(), &top->idef, md.nr);
             }
         }
     }
index 8b13772cdfc58b36ad022021305846b0495a022e..b9fff7fe051271059fc5dffedf9446129b82f2db 100644 (file)
@@ -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/pbcutil/pbc.h"
 #include "gromacs/pbcutil/pbc_simd.h"
 #include "gromacs/simd/simd.h"
@@ -575,7 +574,7 @@ static void do_lincsp(ArrayRefWithPadding<const RVec> xPadded,
                       t_pbc*                          pbc,
                       Lincs*                          lincsd,
                       int                             th,
-                      real*                           invmass,
+                      const real*                     invmass,
                       ConstraintVariable              econq,
                       bool                            bCalcDHDL,
                       bool                            bCalcVir,
@@ -1302,7 +1301,7 @@ static void set_lincs_matrix_task(Lincs* li, Task* li_task, const real* invmass,
 }
 
 /*! \brief Sets the elements in the LINCS matrix. */
-static void set_lincs_matrix(Lincs* li, real* invmass, real lambda)
+static void set_lincs_matrix(Lincs* li, const real* invmass, real lambda)
 {
     const real invsqrt2 = 0.7071067811865475244;
 
@@ -1850,7 +1849,13 @@ static void set_matrix_indices(Lincs* li, const Task& li_task, const ListOfLists
     }
 }
 
-void set_lincs(const InteractionDefinitions& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li)
+void set_lincs(const InteractionDefinitions& idef,
+               const int                     numAtoms,
+               const real*                   invmass,
+               const real                    lambda,
+               bool                          bDynamics,
+               const t_commrec*              cr,
+               Lincs*                        li)
 {
     li->nc_real = 0;
     li->nc      = 0;
@@ -1899,7 +1904,7 @@ void set_lincs(const InteractionDefinitions& idef, const t_mdatoms& md, bool bDy
     }
     else
     {
-        natoms = md.homenr;
+        natoms = numAtoms;
     }
 
     const ListOfLists<int> at2con =
@@ -2123,10 +2128,10 @@ void set_lincs(const InteractionDefinitions& idef, const t_mdatoms& md, bool bDy
 
     if (li->ntask > 1)
     {
-        lincs_thread_setup(li, md.nr);
+        lincs_thread_setup(li, numAtoms);
     }
 
-    set_lincs_matrix(li, md.invmass, md.lambda);
+    set_lincs_matrix(li, invmass, lambda);
 }
 
 //! Issues a warning when LINCS constraints cannot be satisfied.
@@ -2252,7 +2257,7 @@ bool constrain_lincs(bool                            computeRmsd,
                      const t_inputrec&               ir,
                      int64_t                         step,
                      Lincs*                          lincsd,
-                     const t_mdatoms&                md,
+                     const real*                     invmass,
                      const t_commrec*                cr,
                      const gmx_multisim_t*           ms,
                      ArrayRefWithPadding<const RVec> xPadded,
@@ -2260,6 +2265,7 @@ bool constrain_lincs(bool                            computeRmsd,
                      ArrayRef<RVec>                  min_proj,
                      const matrix                    box,
                      t_pbc*                          pbc,
+                     const bool                      hasMassPerturbed,
                      real                            lambda,
                      real*                           dvdlambda,
                      real                            invdt,
@@ -2299,9 +2305,9 @@ bool constrain_lincs(bool                            computeRmsd,
          */
         if (ir.efep != efepNO)
         {
-            if (md.nMassPerturbed && lincsd->matlam != md.lambda)
+            if (hasMassPerturbed && lincsd->matlam != lambda)
             {
-                set_lincs_matrix(lincsd, md.invmass, md.lambda);
+                set_lincs_matrix(lincsd, invmass, lambda);
             }
 
             for (int i = 0; i < lincsd->nc; i++)
@@ -2365,7 +2371,7 @@ bool constrain_lincs(bool                            computeRmsd,
 
                 clear_mat(lincsd->task[th].vir_r_m_dr);
 
-                do_lincs(xPadded, xprimePadded, box, pbc, lincsd, th, md.invmass, cr, bCalcDHDL,
+                do_lincs(xPadded, xprimePadded, box, pbc, lincsd, th, invmass, cr, bCalcDHDL,
                          ir.LincsWarnAngle, &bWarn, invdt, v, bCalcVir,
                          th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
             }
@@ -2444,7 +2450,7 @@ bool constrain_lincs(bool                            computeRmsd,
             {
                 int th = gmx_omp_get_thread_num();
 
-                do_lincsp(xPadded, xprimePadded, min_proj, pbc, lincsd, th, md.invmass, econq,
+                do_lincsp(xPadded, xprimePadded, min_proj, pbc, lincsd, th, invmass, econq,
                           bCalcDHDL, bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
             }
             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
index 2ab8ee8062d321eef6ec1c9082fd3e535a85ebc9..0b00370c018719e946f92c5f214b8d946579c0c2 100644 (file)
@@ -56,7 +56,6 @@ struct gmx_multisim_t;
 class InteractionDefinitions;
 struct t_commrec;
 struct t_inputrec;
-struct t_mdatoms;
 struct t_nrnb;
 struct t_pbc;
 
@@ -91,7 +90,9 @@ void done_lincs(Lincs* li);
 
 /*! \brief Initialize lincs stuff */
 void set_lincs(const InteractionDefinitions& idef,
-               const t_mdatoms&              md,
+               int                           numAtoms,
+               const real*                   invmass,
+               real                          lambda,
                bool                          bDynamics,
                const t_commrec*              cr,
                Lincs*                        li);
@@ -103,7 +104,7 @@ bool constrain_lincs(bool                            computeRmsd,
                      const t_inputrec&               ir,
                      int64_t                         step,
                      Lincs*                          lincsd,
-                     const t_mdatoms&                md,
+                     const real*                     invmass,
                      const t_commrec*                cr,
                      const gmx_multisim_t*           ms,
                      ArrayRefWithPadding<const RVec> x,
@@ -111,6 +112,7 @@ bool constrain_lincs(bool                            computeRmsd,
                      ArrayRef<RVec>                  min_proj,
                      const matrix                    box,
                      t_pbc*                          pbc,
+                     bool                            hasMassPerturbed,
                      real                            lambda,
                      real*                           dvdlambda,
                      real                            invdt,
index 1d6bf1ca7158e88633b58aac5cb549d5304d9b5d..7b3c067d46559184a9d0d725b2355e3f57e4215f 100644 (file)
@@ -725,9 +725,8 @@ bool LincsGpu::isNumCoupledConstraintsSupported(const gmx_mtop_t& mtop)
     return true;
 }
 
-void LincsGpu::set(const InteractionDefinitions& idef, const t_mdatoms& md)
+void LincsGpu::set(const InteractionDefinitions& idef, const int numAtoms, const real* invmass)
 {
-    int numAtoms = md.nr;
     // List of constrained atoms (CPU memory)
     std::vector<int2> constraintsHost;
     // Equilibrium distances for the constraints (CPU)
@@ -861,10 +860,10 @@ void LincsGpu::set(const InteractionDefinitions& idef, const t_mdatoms& md)
 
                 int center = c1a1;
 
-                float sqrtmu1 = 1.0 / sqrt(md.invmass[c1a1] + md.invmass[c1a2]);
-                float sqrtmu2 = 1.0 / sqrt(md.invmass[c2a1] + md.invmass[c2a2]);
+                float sqrtmu1 = 1.0 / sqrt(invmass[c1a1] + invmass[c1a2]);
+                float sqrtmu2 = 1.0 / sqrt(invmass[c2a1] + invmass[c2a2]);
 
-                massFactorsHost.at(index) = -sign * md.invmass[center] * sqrtmu1 * sqrtmu2;
+                massFactorsHost.at(index) = -sign * invmass[center] * sqrtmu1 * sqrtmu2;
 
                 coupledConstraintsCountsHost.at(splitMap.at(c1))++;
             }
@@ -888,10 +887,10 @@ void LincsGpu::set(const InteractionDefinitions& idef, const t_mdatoms& md)
 
                 int center = c1a2;
 
-                float sqrtmu1 = 1.0 / sqrt(md.invmass[c1a1] + md.invmass[c1a2]);
-                float sqrtmu2 = 1.0 / sqrt(md.invmass[c2a1] + md.invmass[c2a2]);
+                float sqrtmu1 = 1.0 / sqrt(invmass[c1a1] + invmass[c1a2]);
+                float sqrtmu2 = 1.0 / sqrt(invmass[c2a1] + invmass[c2a2]);
 
-                massFactorsHost.at(index) = sign * md.invmass[center] * sqrtmu1 * sqrtmu2;
+                massFactorsHost.at(index) = sign * invmass[center] * sqrtmu1 * sqrtmu2;
 
                 coupledConstraintsCountsHost.at(splitMap.at(c1))++;
             }
@@ -958,8 +957,8 @@ void LincsGpu::set(const InteractionDefinitions& idef, const t_mdatoms& md)
                        maxCoupledConstraints * kernelParams_.numConstraintsThreads, deviceStream_,
                        GpuApiCallBehavior::Sync, nullptr);
 
-    GMX_RELEASE_ASSERT(md.invmass != nullptr, "Masses of atoms should be specified.\n");
-    copyToDeviceBuffer(&kernelParams_.d_inverseMasses, md.invmass, 0, numAtoms, deviceStream_,
+    GMX_RELEASE_ASSERT(invmass != nullptr, "Masses of atoms should be specified.\n");
+    copyToDeviceBuffer(&kernelParams_.d_inverseMasses, invmass, 0, numAtoms, deviceStream_,
                        GpuApiCallBehavior::Sync, nullptr);
 }
 
index ef035164318a50bb6397f0daecef1c14e7aaca06..15a962bac1dd05920dfb8a2144265a80d3c4c835 100644 (file)
@@ -47,7 +47,6 @@
 #include "gromacs/gpu_utils/device_context.h"
 #include "gromacs/gpu_utils/gputraits.cuh"
 #include "gromacs/mdlib/constr.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/pbcutil/pbc_aiuc.h"
 #include "gromacs/utility/classhelpers.h"
 
@@ -155,13 +154,12 @@ public:
      * Information about constraints is taken from:
      *     idef.il[F_CONSTR].iatoms  --- type (T) of constraint and two atom indexes (i1, i2)
      *     idef.iparams[T].constr.dA --- target length for constraint of type T
-     * From t_mdatom, the code takes:
-     *     md.invmass  --- array of inverse square root of masses for each atom in the system.
      *
-     * \param[in] idef  Local topology data to get information on constraints from.
-     * \param[in] md    Atoms data to get atom masses from.
+     * \param[in] idef      Local topology data to get information on constraints from.
+     * \param[in] numAtoms  Number of atoms.
+     * \param[in] invmass   Inverse masses of atoms.
      */
-    void set(const InteractionDefinitions& idef, const t_mdatoms& md);
+    void set(const InteractionDefinitions& idef, int numAtoms, const real* invmass);
 
     /*! \brief
      * Returns whether the maximum number of coupled constraints is supported
index 62e351dfe9643bd40427a57d1bd2bb75cbf13270..fd72911b7745479b673c31ecf46c730088353de0 100644 (file)
@@ -58,7 +58,6 @@
 #include "gromacs/mdlib/splitter.h"
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/invblock.h"
 #include "gromacs/utility/fatalerror.h"
@@ -122,7 +121,7 @@ static void resizeLagrangianData(shakedata* shaked, int ncons)
     shaked->scaled_lagrange_multiplier.resize(ncons);
 }
 
-void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md)
+void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const int numAtoms)
 {
     int          i, m, ncons;
     int          bstart, bnr;
@@ -137,7 +136,7 @@ void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, c
 
     init_blocka(&sblocks);
     sfree(sblocks.index); // To solve memory leak
-    gen_sblocks(nullptr, md.homenr, *idef, &sblocks, FALSE);
+    gen_sblocks(nullptr, numAtoms, *idef, &sblocks, FALSE);
 
     /*
        bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
@@ -150,7 +149,7 @@ void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, c
     }
 
     /* Calculate block number for each atom */
-    inv_sblock = make_invblocka(&sblocks, md.nr);
+    inv_sblock = make_invblocka(&sblocks, numAtoms);
 
     done_blocka(&sblocks);
 
index 34db77879c33f0fa80c55f186b9c32e77c3b220c..8559c57553855291095dedf5eadb2247512dc50c 100644 (file)
@@ -52,7 +52,6 @@
 struct InteractionList;
 class InteractionDefinitions;
 struct t_inputrec;
-struct t_mdatoms;
 struct t_nrnb;
 struct t_pbc;
 
@@ -96,7 +95,7 @@ struct shakedata
 };
 
 //! Make SHAKE blocks when not using DD.
-void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md);
+void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, int numAtoms);
 
 //! Make SHAKE blocks when using DD.
 void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon);
index b410504f8d4ff846b2cb551159d0dc0ff407d822..3de0c11b896b37adb798905aded2347f098e6e56 100644 (file)
@@ -101,13 +101,6 @@ ConstraintsTestData::ConstraintsTestData(const std::string&       title,
     ir_.delta_t = timestep;
     ir_.eI      = 0;
 
-    // MD atoms data
-    md_.nMassPerturbed = 0;
-    md_.lambda         = 0.0;
-    md_.invmass        = invmass_.data();
-    md_.nr             = numAtoms;
-    md_.homenr         = numAtoms;
-
     // Virial evaluation
     computeVirial_ = computeVirial;
     if (computeVirial)
@@ -122,8 +115,9 @@ ConstraintsTestData::ConstraintsTestData(const std::string&       title,
         }
     }
 
-
     // Free energy evaluation
+    hasMassPerturbed_  = false;
+    lambda_            = 0.0;
     compute_dHdLambda_ = compute_dHdLambda;
     dHdLambda_         = 0;
     if (compute_dHdLambda_)
index 7374fa19243eb448696be7744910d5bede2ce679..65d6e03f56ce20c3908cc3e4fe0bb636ceb93c7d 100644 (file)
@@ -58,7 +58,6 @@
 #include "gromacs/mdlib/lincs.h"
 #include "gromacs/mdlib/shake.h"
 #include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/idef.h"
 #include "gromacs/topology/ifunc.h"
@@ -92,8 +91,6 @@ public:
     t_inputrec ir_;
     //! Local topology
     std::unique_ptr<InteractionDefinitions> idef_;
-    //! MD atoms
-    t_mdatoms md_;
     //! Computational time array (normally used to benchmark performance)
     t_nrnb nrnb_;
 
@@ -109,6 +106,10 @@ public:
     tensor virialScaledRef_;
     //! If the free energy is computed
     bool compute_dHdLambda_;
+    //! If there are atoms with perturbed mass
+    bool hasMassPerturbed_ = false;
+    //! Lambda value
+    real lambda_ = 0.0;
     //! For free energy computation
     real dHdLambda_;
     //! For free energy computation (reference value)
index 9ef51bdee7536aa588fde40c743b56613abe28e9..1700e6e90b4495a14a9a6bebeefd02221392c896 100644 (file)
@@ -66,7 +66,6 @@
 #include "gromacs/mdrunutility/multisim.h"
 #include "gromacs/mdtypes/commrec.h"
 #include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/idef.h"
 #include "gromacs/topology/ifunc.h"
@@ -90,10 +89,10 @@ namespace test
 void applyShake(ConstraintsTestData* testData, t_pbc gmx_unused pbc)
 {
     shakedata shaked;
-    make_shake_sblock_serial(&shaked, testData->idef_.get(), testData->md_);
+    make_shake_sblock_serial(&shaked, testData->idef_.get(), testData->numAtoms_);
     bool success = constrain_shake(
             nullptr, &shaked, testData->invmass_.data(), *testData->idef_, testData->ir_, testData->x_,
-            testData->xPrime_, testData->xPrime2_, nullptr, &testData->nrnb_, testData->md_.lambda,
+            testData->xPrime_, testData->xPrime2_, nullptr, &testData->nrnb_, testData->lambda_,
             &testData->dHdLambda_, testData->invdt_, testData->v_, testData->computeVirial_,
             testData->virialScaled_, false, gmx::ConstraintVariable::Positions);
     EXPECT_TRUE(success) << "Test failed with a false return value in SHAKE.";
@@ -133,14 +132,15 @@ void applyLincs(ConstraintsTestData* testData, t_pbc pbc)
     // Initialize LINCS
     lincsd = init_lincs(nullptr, testData->mtop_, testData->nflexcon_, at2con_mt, false,
                         testData->ir_.nLincsIter, testData->ir_.nProjOrder);
-    set_lincs(*testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &cr, lincsd);
+    set_lincs(*testData->idef_, testData->numAtoms_, testData->invmass_.data(), testData->lambda_,
+              EI_DYNAMICS(testData->ir_.eI), &cr, lincsd);
 
     // Evaluate constraints
     bool success = constrain_lincs(
-            false, testData->ir_, 0, lincsd, testData->md_, &cr, &ms,
+            false, testData->ir_, 0, lincsd, testData->invmass_.data(), &cr, &ms,
             testData->x_.arrayRefWithPadding(), testData->xPrime_.arrayRefWithPadding(),
             testData->xPrime2_.arrayRefWithPadding().unpaddedArrayRef(), pbc.box, &pbc,
-            testData->md_.lambda, &testData->dHdLambda_, testData->invdt_,
+            testData->hasMassPerturbed_, testData->lambda_, &testData->dHdLambda_, testData->invdt_,
             testData->v_.arrayRefWithPadding().unpaddedArrayRef(), testData->computeVirial_,
             testData->virialScaled_, gmx::ConstraintVariable::Positions, &testData->nrnb_, maxwarn,
             &warncount_lincs);
index 8c2385daabfa75ff78fec6a939a2056f5cdcb1e0..161c63b5d525c5e13948cbce974a370f8f1ad813 100644 (file)
@@ -81,7 +81,7 @@ void applyLincsGpu(ConstraintsTestData* testData, t_pbc pbc)
     int     numAtoms         = testData->numAtoms_;
     float3 *d_x, *d_xp, *d_v;
 
-    lincsGpu->set(*testData->idef_, testData->md_);
+    lincsGpu->set(*testData->idef_, testData->numAtoms_, testData->invmass_.data());
     PbcAiuc pbcAiuc;
     setPbcAiuc(pbc.ndim_ePBC, pbc.box, &pbcAiuc);
 
index 76899bbd82c21ed484d8f7f4977b524d62f90295..b40978bf90bae69f4c1d4305dcf6058415388dd3 100644 (file)
@@ -214,7 +214,7 @@ void UpdateConstrainGpu::Impl::set(DeviceBuffer<RVec>            d_x,
 
     // Integrator should also update something, but it does not even have a method yet
     integrator_->set(md, numTempScaleValues, md.cTC);
-    lincsGpu_->set(idef, md);
+    lincsGpu_->set(idef, md.nr, md.invmass);
     settleGpu_->set(idef, md);
 
     coordinateScalingKernelLaunchConfig_.gridSize[0] =