Reimplement constant acceleration groups
[alexxy/gromacs.git] / src / gromacs / mdlib / mdatoms.cpp
index 0d1abd8565da3830a1ad7b112cc5b99bcffb4ccd..e4ee7b7b57abc81466d0d94f3b80790d1e5917a4 100644 (file)
@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
- * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2017,2018,2019,2020,2021, 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.
@@ -148,7 +148,7 @@ std::unique_ptr<MDAtoms> makeMDAtoms(FILE* fp, const gmx_mtop_t& mtop, const t_i
     double totalMassB = 0.0;
 
     md->haveVsites                  = FALSE;
-    gmx_mtop_atomloop_block_t aloop = gmx_mtop_atomloop_block_init(&mtop);
+    gmx_mtop_atomloop_block_t aloop = gmx_mtop_atomloop_block_init(mtop);
     const t_atom*             atom;
     int                       nmol;
     while (gmx_mtop_atomloop_block_next(aloop, &atom, &nmol))
@@ -156,12 +156,12 @@ std::unique_ptr<MDAtoms> makeMDAtoms(FILE* fp, const gmx_mtop_t& mtop, const t_i
         totalMassA += nmol * atom->m;
         totalMassB += nmol * atom->mB;
 
-        if (atom->ptype == eptVSite)
+        if (atom->ptype == ParticleType::VSite)
         {
             md->haveVsites = TRUE;
         }
 
-        if (ir.efep != efepNO && PERTURBED(*atom))
+        if (ir.efep != FreeEnergyPerturbationType::No && PERTURBED(*atom))
         {
             md->nPerturbed++;
             if (atom->mB != atom->m)
@@ -182,10 +182,12 @@ std::unique_ptr<MDAtoms> makeMDAtoms(FILE* fp, const gmx_mtop_t& mtop, const t_i
     md->tmassA = totalMassA;
     md->tmassB = totalMassB;
 
-    if (ir.efep != efepNO && fp)
+    if (ir.efep != FreeEnergyPerturbationType::No && fp)
     {
-        fprintf(fp, "There are %d atoms and %d charges for free energy perturbation\n",
-                md->nPerturbed, md->nChargePerturbed);
+        fprintf(fp,
+                "There are %d atoms and %d charges for free energy perturbation\n",
+                md->nPerturbed,
+                md->nChargePerturbed);
     }
 
     md->havePartiallyFrozenAtoms = FALSE;
@@ -200,15 +202,15 @@ std::unique_ptr<MDAtoms> makeMDAtoms(FILE* fp, const gmx_mtop_t& mtop, const t_i
         }
     }
 
-    md->bOrires = (gmx_mtop_ftype_count(&mtop, F_ORIRES) != 0);
+    md->bOrires = (gmx_mtop_ftype_count(mtop, F_ORIRES) != 0);
 
     return mdAtoms;
 }
 
 } // namespace gmx
 
-void atoms2md(const gmx_mtop_t*  mtop,
-              const t_inputrec*  ir,
+void atoms2md(const gmx_mtop_t&  mtop,
+              const t_inputrec&  inputrec,
               int                nindex,
               gmx::ArrayRef<int> index,
               int                homenr,
@@ -216,15 +218,15 @@ void atoms2md(const gmx_mtop_t*  mtop,
 {
     gmx_bool         bLJPME;
     const t_grpopts* opts;
-    int nthreads gmx_unused;
+    int nthreads     gmx_unused;
 
-    bLJPME = EVDW_PME(ir->vdwtype);
+    bLJPME = EVDW_PME(inputrec.vdwtype);
 
-    opts = &ir->opts;
+    opts = &inputrec.opts;
 
-    const SimulationGroups& groups = mtop->groups;
+    const SimulationGroups& groups = mtop.groups;
 
-    auto md = mdAtoms->mdatoms();
+    auto* md = mdAtoms->mdatoms();
     /* nindex>=0 indicates DD where we use an index */
     if (nindex >= 0)
     {
@@ -232,7 +234,7 @@ void atoms2md(const gmx_mtop_t*  mtop,
     }
     else
     {
-        md->nr = mtop->natoms;
+        md->nr = mtop.natoms;
     }
 
     if (md->nr > md->nalloc)
@@ -287,12 +289,11 @@ void atoms2md(const gmx_mtop_t*  mtop,
             /* We always copy cTC with domain decomposition */
         }
         srenew(md->cENER, md->nalloc);
-        if (opts->ngacc > 1)
+        if (inputrec.useConstantAcceleration)
         {
             srenew(md->cACC, md->nalloc);
         }
-        if (opts->nFreeze
-            && (opts->ngfrz > 1 || opts->nFreeze[0][XX] || opts->nFreeze[0][YY] || opts->nFreeze[0][ZZ]))
+        if (inputrecFrozenAtoms(&inputrec))
         {
             srenew(md->cFREEZE, md->nalloc);
         }
@@ -314,11 +315,11 @@ void atoms2md(const gmx_mtop_t*  mtop,
          * Therefore, when adding code, the user should use something like:
          * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])
          */
-        if (!mtop->groups.groupNumbers[SimulationAtomGroupType::User1].empty())
+        if (!mtop.groups.groupNumbers[SimulationAtomGroupType::User1].empty())
         {
             srenew(md->cU1, md->nalloc);
         }
-        if (!mtop->groups.groupNumbers[SimulationAtomGroupType::User2].empty())
+        if (!mtop.groups.groupNumbers[SimulationAtomGroupType::User2].empty())
         {
             srenew(md->cU2, md->nalloc);
         }
@@ -326,7 +327,7 @@ void atoms2md(const gmx_mtop_t*  mtop,
 
     int molb = 0;
 
-    nthreads = gmx_omp_nthreads_get(emntDefault);
+    nthreads = gmx_omp_nthreads_get(ModuleMultiThread::Default);
 #pragma omp parallel for num_threads(nthreads) schedule(static) firstprivate(molb)
     for (int i = 0; i < md->nr; i++)
     {
@@ -350,13 +351,13 @@ void atoms2md(const gmx_mtop_t*  mtop,
             {
                 md->cFREEZE[i] = getGroupType(groups, SimulationAtomGroupType::Freeze, ag);
             }
-            if (EI_ENERGY_MINIMIZATION(ir->eI))
+            if (EI_ENERGY_MINIMIZATION(inputrec.eI))
             {
                 /* Displacement is proportional to F, masses used for constraints */
                 mA = 1.0;
                 mB = 1.0;
             }
-            else if (ir->eI == eiBD)
+            else if (inputrec.eI == IntegrationAlgorithm::BD)
             {
                 /* With BD the physical masses are irrelevant.
                  * To keep the code simple we use most of the normal MD code path
@@ -368,15 +369,15 @@ void atoms2md(const gmx_mtop_t*  mtop,
                  * Thus with BD v*dt will give the displacement and the reported
                  * temperature can signal bad integration (too large time step).
                  */
-                if (ir->bd_fric > 0)
+                if (inputrec.bd_fric > 0)
                 {
-                    mA = 0.5 * ir->bd_fric * ir->delta_t;
-                    mB = 0.5 * ir->bd_fric * ir->delta_t;
+                    mA = 0.5 * inputrec.bd_fric * inputrec.delta_t;
+                    mB = 0.5 * inputrec.bd_fric * inputrec.delta_t;
                 }
                 else
                 {
                     /* The friction coefficient is mass/tau_t */
-                    fac = ir->delta_t
+                    fac = inputrec.delta_t
                           / opts->tau_t[md->cTC ? groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ag] : 0];
                     mA = 0.5 * atom.m * fac;
                     mB = 0.5 * atom.mB * fac;
@@ -404,8 +405,7 @@ void atoms2md(const gmx_mtop_t*  mtop,
             else if (md->cFREEZE)
             {
                 g = md->cFREEZE[i];
-                GMX_ASSERT(opts->nFreeze != nullptr,
-                           "Must have freeze groups to initialize masses");
+                GMX_ASSERT(opts->nFreeze != nullptr, "Must have freeze groups to initialize masses");
                 if (opts->nFreeze[g][XX] && opts->nFreeze[g][YY] && opts->nFreeze[g][ZZ])
                 {
                     /* Set the mass of completely frozen particles to ALMOST_ZERO
@@ -439,9 +439,9 @@ void atoms2md(const gmx_mtop_t*  mtop,
             md->typeA[i]   = atom.type;
             if (bLJPME)
             {
-                c6  = mtop->ffparams.iparams[atom.type * (mtop->ffparams.atnr + 1)].lj.c6;
-                c12 = mtop->ffparams.iparams[atom.type * (mtop->ffparams.atnr + 1)].lj.c12;
-                md->sqrt_c6A[i] = sqrt(c6);
+                c6  = mtop.ffparams.iparams[atom.type * (mtop.ffparams.atnr + 1)].lj.c6;
+                c12 = mtop.ffparams.iparams[atom.type * (mtop.ffparams.atnr + 1)].lj.c12;
+                md->sqrt_c6A[i] = std::sqrt(c6);
                 if (c6 == 0.0 || c12 == 0)
                 {
                     md->sigmaA[i] = 1.0;
@@ -459,9 +459,9 @@ void atoms2md(const gmx_mtop_t*  mtop,
                 md->typeB[i]      = atom.typeB;
                 if (bLJPME)
                 {
-                    c6  = mtop->ffparams.iparams[atom.typeB * (mtop->ffparams.atnr + 1)].lj.c6;
-                    c12 = mtop->ffparams.iparams[atom.typeB * (mtop->ffparams.atnr + 1)].lj.c12;
-                    md->sqrt_c6B[i] = sqrt(c6);
+                    c6  = mtop.ffparams.iparams[atom.typeB * (mtop.ffparams.atnr + 1)].lj.c6;
+                    c12 = mtop.ffparams.iparams[atom.typeB * (mtop.ffparams.atnr + 1)].lj.c12;
+                    md->sqrt_c6B[i] = std::sqrt(c6);
                     if (c6 == 0.0 || c12 == 0)
                     {
                         md->sigmaB[i] = 1.0;
@@ -528,7 +528,7 @@ void update_mdatoms(t_mdatoms* md, real lambda)
         real L1 = 1 - lambda;
 
         /* Update masses of perturbed atoms for the change in lambda */
-        int gmx_unused nthreads = gmx_omp_nthreads_get(emntDefault);
+        int gmx_unused nthreads = gmx_omp_nthreads_get(ModuleMultiThread::Default);
 #pragma omp parallel for num_threads(nthreads) schedule(static)
         for (int i = 0; i < md->nr; i++)
         {