Update DD cutoffs for update groups

Due to the distance between atoms to the center of geometry of the
update group they belong to, there is now a difference between
atom-atom, domain-domain and atom-domain cutoffs.
Also adds the cutoff versus box check for update groups.

Change-Id: I474cd1371a1fac22f26b80adcc235bfb10764d34

diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp
index 38de346839096ef2eb89e692196e0aefda377300..1cc507d119eee158a85a50c8dee36689c765a7f6 100644 (file)
--- a/src/gromacs/domdec/domdec.cpp
+++ b/src/gromacs/domdec/domdec.cpp
@@ -2073,6 +2073,7 @@ static bool systemHasConstraintsOrVsites(const gmx_mtop_t &mtop)
 static void setupUpdateGroups(const gmx::MDLogger &mdlog,
                               const gmx_mtop_t    &mtop,
                               const t_inputrec    &inputrec,
+                              real                 cutoffMargin,
                               int                  numMpiRanksTotal,
                               gmx_domdec_comm_t   *comm)
 {
@@ -2110,7 +2111,8 @@ static void setupUpdateGroups(const gmx::MDLogger &mdlog,
         /* To use update groups, the large domain-to-domain cutoff distance
          * should be compatible with the box size.
          */
-        /* TODO: add this cutoff check */
+        comm->useUpdateGroups = (atomToAtomIntoDomainToDomainCutoff(*comm, 0) < cutoffMargin);
+        /* TODO: Enable update groups when all infrastructure is present */
         comm->useUpdateGroups = false;
 
         if (comm->useUpdateGroups)
@@ -2180,7 +2182,8 @@ static void set_dd_limits_and_grid(const gmx::MDLogger &mdlog,
     comm->useUpdateGroups = false;
     if (ir->cutoff_scheme == ecutsVERLET)
     {
-        setupUpdateGroups(mdlog, *mtop, *ir, cr->nnodes, comm);
+        real cutoffMargin = std::sqrt(max_cutoff2(ir->ePBC, box)) - ir->rlist;
+        setupUpdateGroups(mdlog, *mtop, *ir, cutoffMargin, cr->nnodes, comm);
     }
 
     comm->bInterCGBondeds = ((ncg_mtop(mtop) > gmx_mtop_num_molecules(*mtop)) ||
@@ -2194,8 +2197,16 @@ static void set_dd_limits_and_grid(const gmx::MDLogger &mdlog,
         comm->bInterCGMultiBody = FALSE;
     }
 
-    dd->splitConstraints = gmx::inter_charge_group_constraints(*mtop);
-    dd->splitSettles     = gmx::inter_charge_group_settles(*mtop);
+    if (comm->useUpdateGroups)
+    {
+        dd->splitConstraints = false;
+        dd->splitSettles     = false;
+    }
+    else
+    {
+        dd->splitConstraints = gmx::inter_charge_group_constraints(*mtop);
+        dd->splitSettles     = gmx::inter_charge_group_settles(*mtop);
+    }
 
     if (ir->rlist == 0)
     {
@@ -2207,7 +2218,7 @@ static void set_dd_limits_and_grid(const gmx::MDLogger &mdlog,
     }
     else
     {
-        comm->cutoff   = ir->rlist;
+        comm->cutoff   = atomToAtomIntoDomainToDomainCutoff(*comm, ir->rlist);
     }
     comm->cutoff_mbody = 0;
 
@@ -2245,7 +2256,7 @@ static void set_dd_limits_and_grid(const gmx::MDLogger &mdlog,
     {
         if (options.minimumCommunicationRange > 0)
         {
-            comm->cutoff_mbody = options.minimumCommunicationRange;
+            comm->cutoff_mbody = atomToAtomIntoDomainToDomainCutoff(*comm, options.minimumCommunicationRange);
             if (options.useBondedCommunication)
             {
                 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);

diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp
index 3e98a16f5261f07113f5cd7641126e5d6806de6c..a46518dfa4632e8ef53898423cf13a166e636f56 100644 (file)
--- a/src/gromacs/domdec/domdec_topology.cpp
+++ b/src/gromacs/domdec/domdec_topology.cpp
@@ -1278,33 +1278,12 @@ check_assign_interactions_atom(int i, int i_gl,
     j = ind_start;
     while (j < ind_end)
     {
-        int            ftype;
-        int            nral;
-        t_iatom        tiatoms[1 + MAXATOMLIST];
-
-        ftype  = rtil[j++];
-        gmx::ArrayRef<const t_iatom> iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j);
-        nral   = NRAL(ftype);
-        if (ftype == F_SETTLE)
-        {
-            /* Settles are only in the reverse top when they
-             * operate within a charge group. So we can assign
-             * them without checks. We do this only for performance
-             * reasons; it could be handled by the code below.
-             */
-            if (iz == 0)
-            {
-                /* Home zone: add this settle to the local topology */
-                tiatoms[0] = iatoms[0];
-                tiatoms[1] = i;
-                tiatoms[2] = i + iatoms[2] - iatoms[1];
-                tiatoms[3] = i + iatoms[3] - iatoms[1];
-                add_ifunc(nral, tiatoms, &idef->il[ftype]);
-                (*nbonded_local)++;
-            }
-            j += 1 + nral;
-        }
-        else if (interaction_function[ftype].flags & IF_VSITE)
+        t_iatom   tiatoms[1 + MAXATOMLIST];
+
+        const int ftype  = rtil[j++];
+        auto      iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j);
+        const int nral   = NRAL(ftype);
+        if (interaction_function[ftype].flags & IF_VSITE)
         {
             assert(!bInterMolInteractions);
             /* The vsite construction goes where the vsite itself is */

diff --git a/src/gromacs/domdec/partition.cpp b/src/gromacs/domdec/partition.cpp
index 36582eebdf1d6dad723bca64107c305400c5f7a2..6d4218cde5aabee79987ccc98d289d11588db2ed 100644 (file)
--- a/src/gromacs/domdec/partition.cpp
+++ b/src/gromacs/domdec/partition.cpp
@@ -1892,7 +1892,6 @@ static void setup_dd_communication(gmx_domdec_t *dd,
     gmx_domdec_comm_dim_t *cd;
     cginfo_mb_t           *cginfo_mb;
     gmx_bool               bBondComm, bDist2B, bDistMB, bDistBonded;
-    real                   r_comm2, r_bcomm2;
     dd_corners_t           corners;
     rvec                  *cg_cm, *normal, *v_d, *v_0 = nullptr, *v_1 = nullptr;
     real                   skew_fac2_d, skew_fac_01;
@@ -1935,8 +1934,8 @@ static void setup_dd_communication(gmx_domdec_t *dd,
     /* Do we need to determine extra distances for only two-body bondeds? */
     bDist2B = (bBondComm && !bDistMB);
 
-    r_comm2  = gmx::square(comm->cutoff);
-    r_bcomm2 = gmx::square(comm->cutoff_mbody);
+    const real r_comm2  = gmx::square(domainToDomainIntoAtomToDomainCutoff(*comm, comm->cutoff));
+    const real r_bcomm2 = gmx::square(domainToDomainIntoAtomToDomainCutoff(*comm, comm->cutoff_mbody));
 
     if (debug)
     {

diff --git a/src/gromacs/domdec/utility.h b/src/gromacs/domdec/utility.h
index dc0a6654a52a3d2912762c06ccb65aaab8ab5f72..cab76565649d315deec601f5d57425a0d0667c00 100644 (file)
--- a/src/gromacs/domdec/utility.h
+++ b/src/gromacs/domdec/utility.h
@@ -100,4 +100,34 @@ void dd_check_alloc_ncg(t_forcerec       *fr,
                         PaddedRVecVector *f,
                         int               numChargeGroups);
 
+/*! \brief Returns a domain-to-domain cutoff distance given an atom-to-atom cutoff */
+static inline real
+atomToAtomIntoDomainToDomainCutoff(const gmx_domdec_comm_t &comm,
+                                   real                     cutoff)
+{
+    if (comm.useUpdateGroups)
+    {
+        GMX_ASSERT(comm.updateGroupsCog, "updateGroupsCog should be initialized here");
+
+        cutoff += 2*comm.updateGroupsCog->maxUpdateGroupRadius();
+    }
+
+    return cutoff;
+}
+
+/*! \brief Returns an atom-to-domain cutoff distance given a domain-to-domain cutoff */
+static inline real
+domainToDomainIntoAtomToDomainCutoff(const gmx_domdec_comm_t &comm,
+                                     real                     cutoff)
+{
+    if (comm.useUpdateGroups)
+    {
+        GMX_ASSERT(comm.updateGroupsCog, "updateGroupsCog should be initialized here");
+
+        cutoff -= comm.updateGroupsCog->maxUpdateGroupRadius();
+    }
+
+    return cutoff;
+}
+
 #endif