From: Joe Jordan <ejjordan12@gmail.com>
Date: Mon, 12 Apr 2021 10:24:29 +0000 (+0000)
Subject: Replace defines with constexpr in ishift
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=54e7cd7f3abf36752aa1ba96e89476e23c01cb87;p=alexxy%2Fgromacs.git

Replace defines with constexpr in ishift
---

diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp
index d941620c12..8bec7e2409 100644
--- a/src/gromacs/domdec/domdec.cpp
+++ b/src/gromacs/domdec/domdec.cpp
@@ -371,7 +371,7 @@ void dd_move_f(gmx_domdec_t* dd, gmx::ForceWithShiftForces* forceWithShiftForces
         /* Determine which shift vector we need */
         ivec vis        = { 0, 0, 0 };
         vis[dd->dim[d]] = 1;
-        const int is    = IVEC2IS(vis);
+        const int is    = gmx::ivecToShiftIndex(vis);
 
         /* Loop over the pulses */
         const gmx_domdec_comm_dim_t& cd = comm.cd[d];
diff --git a/src/gromacs/domdec/domdec_specatomcomm.cpp b/src/gromacs/domdec/domdec_specatomcomm.cpp
index 5a87c7d339..948041d7d7 100644
--- a/src/gromacs/domdec/domdec_specatomcomm.cpp
+++ b/src/gromacs/domdec/domdec_specatomcomm.cpp
@@ -112,7 +112,7 @@ void dd_move_f_specat(const gmx_domdec_t* dd, gmx_domdec_specat_comm_t* spac, rv
                 {
                     clear_ivec(vis);
                     vis[dim] = (dir == 0 ? 1 : -1);
-                    int is   = IVEC2IS(vis);
+                    int is   = gmx::ivecToShiftIndex(vis);
                     if (!bScrew)
                     {
                         /* Sum and add to shift forces */
diff --git a/src/gromacs/gmxana/gmx_disre.cpp b/src/gromacs/gmxana/gmx_disre.cpp
index 54865bc612..3490334d18 100644
--- a/src/gromacs/gmxana/gmx_disre.cpp
+++ b/src/gromacs/gmxana/gmx_disre.cpp
@@ -236,7 +236,7 @@ static void check_viol(FILE*                          log,
         dr[clust_id].aver_3[ndr] += drt;
         dr[clust_id].aver_6[ndr] += disresdata->Rt_6[label];
 
-        snew(fshift, SHIFTS);
+        snew(fshift, gmx::c_numShiftVectors);
         ta_disres(n, &forceatoms[i], forceparams.data(), x, f, fshift, pbc, lam, &dvdl, {}, nullptr, disresdata, nullptr, nullptr);
         sfree(fshift);
         viol = disresdata->sumviol;
diff --git a/src/gromacs/listed_forces/bonded.cpp b/src/gromacs/listed_forces/bonded.cpp
index fbddaf5621..87c06ee74f 100644
--- a/src/gromacs/listed_forces/bonded.cpp
+++ b/src/gromacs/listed_forces/bonded.cpp
@@ -135,7 +135,7 @@ int pbc_rvec_sub(const t_pbc* pbc, const rvec xi, const rvec xj, rvec dx)
     else
     {
         rvec_sub(xi, xj, dx);
-        return CENTRAL;
+        return c_centralShiftIndex;
     }
 }
 
@@ -231,7 +231,7 @@ inline void spreadBondForces(const real bondForce,
         if (computeVirial(flavor))
         {
             fshift[shiftIndex][m] += fij;
-            fshift[CENTRAL][m] -= fij;
+            fshift[c_centralShiftIndex][m] -= fij;
         }
     }
 }
@@ -915,7 +915,7 @@ real water_pol(int             nbonds,
                 if (computeVirial(flavor))
                 {
                     fshift[ki][m] += fij;
-                    fshift[CENTRAL][m] -= fij;
+                    fshift[c_centralShiftIndex][m] -= fij;
                 }
             }
         }
@@ -949,7 +949,7 @@ do_1_thole(const rvec xi, const rvec xj, rvec fi, rvec fj, const t_pbc* pbc, rea
         if (computeVirial(flavor))
         {
             fshift[t][m] += fff;
-            fshift[CENTRAL][m] -= fff;
+            fshift[c_centralShiftIndex][m] -= fff;
         }
     } /* 15 */
 
@@ -1084,7 +1084,7 @@ angles(int             nbonds,
             if (computeVirial(flavor))
             {
                 rvec_inc(fshift[t1], f_i);
-                rvec_inc(fshift[CENTRAL], f_j);
+                rvec_inc(fshift[c_centralShiftIndex], f_j);
                 rvec_inc(fshift[t2], f_k);
             }
         } /* 161 TOTAL	*/
@@ -1328,7 +1328,7 @@ real linear_angles(int             nbonds,
         if (computeVirial(flavor))
         {
             rvec_inc(fshift[t1], f_i);
-            rvec_inc(fshift[CENTRAL], f_j);
+            rvec_inc(fshift[c_centralShiftIndex], f_j);
             rvec_inc(fshift[t2], f_k);
         }
     } /* 57 TOTAL	*/
@@ -1414,7 +1414,7 @@ urey_bradley(int             nbonds,
             if (computeVirial(flavor))
             {
                 rvec_inc(fshift[t1], f_i);
-                rvec_inc(fshift[CENTRAL], f_j);
+                rvec_inc(fshift[c_centralShiftIndex], f_j);
                 rvec_inc(fshift[t2], f_k);
             }
         } /* 161 TOTAL	*/
@@ -1435,7 +1435,7 @@ urey_bradley(int             nbonds,
             if (computeVirial(flavor))
             {
                 fshift[ki][m] += fik;
-                fshift[CENTRAL][m] -= fik;
+                fshift[c_centralShiftIndex][m] -= fik;
             }
         }
     }
@@ -1669,7 +1669,7 @@ real quartic_angles(int             nbonds,
             if (computeVirial(flavor))
             {
                 rvec_inc(fshift[t1], f_i);
-                rvec_inc(fshift[CENTRAL], f_j);
+                rvec_inc(fshift[c_centralShiftIndex], f_j);
                 rvec_inc(fshift[t2], f_k);
             }
         } /* 153 TOTAL	*/
@@ -1858,11 +1858,11 @@ void do_dih_fup(int          i,
             }
             else
             {
-                t3 = CENTRAL;
+                t3 = c_centralShiftIndex;
             }
 
             rvec_inc(fshift[t1], f_i);
-            rvec_dec(fshift[CENTRAL], f_j);
+            rvec_dec(fshift[c_centralShiftIndex], f_j);
             rvec_dec(fshift[t2], f_k);
             rvec_inc(fshift[t3], f_l);
         }
@@ -2421,11 +2421,11 @@ real low_angres(int             nbonds,
             if (computeVirial(flavor))
             {
                 rvec_inc(fshift[t1], f_i);
-                rvec_dec(fshift[CENTRAL], f_i);
+                rvec_dec(fshift[c_centralShiftIndex], f_i);
                 if (!bZAxis)
                 {
                     rvec_inc(fshift[t2], f_k);
-                    rvec_dec(fshift[CENTRAL], f_k);
+                    rvec_dec(fshift[c_centralShiftIndex], f_k);
                 }
             }
         }
@@ -2682,7 +2682,7 @@ real restrangles(int             nbonds,
         if (computeVirial(flavor))
         {
             rvec_inc(fshift[t1], f_i);
-            rvec_inc(fshift[CENTRAL], f_j);
+            rvec_inc(fshift[c_centralShiftIndex], f_j);
             rvec_inc(fshift[t2], f_k);
         }
     }
@@ -2800,11 +2800,11 @@ real restrdihs(int             nbonds,
             }
             else
             {
-                t3 = CENTRAL;
+                t3 = c_centralShiftIndex;
             }
 
             rvec_inc(fshift[t1], f_i);
-            rvec_inc(fshift[CENTRAL], f_j);
+            rvec_inc(fshift[c_centralShiftIndex], f_j);
             rvec_inc(fshift[t2], f_k);
             rvec_inc(fshift[t3], f_l);
         }
@@ -2917,11 +2917,11 @@ real cbtdihs(int             nbonds,
             }
             else
             {
-                t3 = CENTRAL;
+                t3 = c_centralShiftIndex;
             }
 
             rvec_inc(fshift[t1], f_i);
-            rvec_inc(fshift[CENTRAL], f_j);
+            rvec_inc(fshift[c_centralShiftIndex], f_j);
             rvec_inc(fshift[t2], f_k);
             rvec_inc(fshift[t3], f_l);
         }
@@ -3437,17 +3437,17 @@ real cmap_dihs(int                 nbonds,
             }
             else
             {
-                t31 = CENTRAL;
-                t32 = CENTRAL;
+                t31 = c_centralShiftIndex;
+                t32 = c_centralShiftIndex;
             }
 
             rvec_inc(fshift[t11], f1_i);
-            rvec_inc(fshift[CENTRAL], f1_j);
+            rvec_inc(fshift[c_centralShiftIndex], f1_j);
             rvec_inc(fshift[t21], f1_k);
             rvec_inc(fshift[t31], f1_l);
 
             rvec_inc(fshift[t12], f2_i);
-            rvec_inc(fshift[CENTRAL], f2_j);
+            rvec_inc(fshift[c_centralShiftIndex], f2_j);
             rvec_inc(fshift[t22], f2_k);
             rvec_inc(fshift[t32], f2_l);
         }
@@ -3609,7 +3609,7 @@ real g96angles(int             nbonds,
         if (computeVirial(flavor))
         {
             rvec_inc(fshift[t1], f_i);
-            rvec_inc(fshift[CENTRAL], f_j);
+            rvec_inc(fshift[c_centralShiftIndex], f_j);
             rvec_inc(fshift[t2], f_k); /* 9 */
         }
         /* 163 TOTAL	*/
@@ -3683,7 +3683,7 @@ real cross_bond_bond(int             nbonds,
         if (computeVirial(flavor))
         {
             rvec_inc(fshift[t1], f_i);
-            rvec_inc(fshift[CENTRAL], f_j);
+            rvec_inc(fshift[c_centralShiftIndex], f_j);
             rvec_inc(fshift[t2], f_k); /* 9 */
         }
         /* 163 TOTAL	*/
@@ -3767,7 +3767,7 @@ real cross_bond_angle(int             nbonds,
         if (computeVirial(flavor))
         {
             rvec_inc(fshift[t1], f_i);
-            rvec_inc(fshift[CENTRAL], f_j);
+            rvec_inc(fshift[c_centralShiftIndex], f_j);
             rvec_inc(fshift[t2], f_k); /* 9 */
         }
         /* 163 TOTAL	*/
@@ -3960,7 +3960,7 @@ real tab_angles(int             nbonds,
             if (computeVirial(flavor))
             {
                 rvec_inc(fshift[t1], f_i);
-                rvec_inc(fshift[CENTRAL], f_j);
+                rvec_inc(fshift[c_centralShiftIndex], f_j);
                 rvec_inc(fshift[t2], f_k);
             }
         } /* 169 TOTAL	*/
diff --git a/src/gromacs/listed_forces/disre.cpp b/src/gromacs/listed_forces/disre.cpp
index de1928792a..981cbefa47 100644
--- a/src/gromacs/listed_forces/disre.cpp
+++ b/src/gromacs/listed_forces/disre.cpp
@@ -585,7 +585,7 @@ real ta_disres(int              nfa,
             int pair = (faOffset + fa) / 3;
             int ai   = forceatoms[fa + 1];
             int aj   = forceatoms[fa + 2];
-            int ki   = CENTRAL;
+            int ki   = gmx::c_centralShiftIndex;
             if (pbc)
             {
                 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
@@ -624,7 +624,7 @@ real ta_disres(int              nfa,
                 if (fshift)
                 {
                     fshift[ki][m] += fij;
-                    fshift[CENTRAL][m] -= fij;
+                    fshift[gmx::c_centralShiftIndex][m] -= fij;
                 }
             }
         }
diff --git a/src/gromacs/listed_forces/gpubonded_impl.cu b/src/gromacs/listed_forces/gpubonded_impl.cu
index 0689950b56..4a8bbe41ce 100644
--- a/src/gromacs/listed_forces/gpubonded_impl.cu
+++ b/src/gromacs/listed_forces/gpubonded_impl.cu
@@ -78,9 +78,10 @@ GpuBonded::Impl::Impl(const gmx_ffparams_t& ffparams,
     GMX_RELEASE_ASSERT(deviceStream.isValid(),
                        "Can't run GPU version of bonded forces in stream that is not valid.");
 
-    static_assert(c_threadsPerBlock >= SHIFTS,
-                  "Threads per block in GPU bonded must be >= SHIFTS for the virial kernel "
-                  "(calcVir=true)");
+    static_assert(
+            c_threadsPerBlock >= c_numShiftVectors,
+            "Threads per block in GPU bonded must be >= c_numShiftVectors for the virial kernel "
+            "(calcVir=true)");
 
     wcycle_ = wcycle;
 
@@ -121,7 +122,7 @@ GpuBonded::Impl::Impl(const gmx_ffparams_t& ffparams,
     kernelLaunchConfig_.gridSize[1]  = 1;
     kernelLaunchConfig_.gridSize[2]  = 1;
     kernelLaunchConfig_.sharedMemorySize =
-            SHIFTS * sizeof(float3) + (c_threadsPerBlock / warp_size) * 3 * sizeof(float);
+            c_numShiftVectors * sizeof(float3) + (c_threadsPerBlock / warp_size) * 3 * sizeof(float);
 }
 
 GpuBonded::Impl::~Impl()
diff --git a/src/gromacs/listed_forces/gpubondedkernels.cu b/src/gromacs/listed_forces/gpubondedkernels.cu
index 42d24e40e1..57fbbc1be4 100644
--- a/src/gromacs/listed_forces/gpubondedkernels.cu
+++ b/src/gromacs/listed_forces/gpubondedkernels.cu
@@ -132,10 +132,10 @@ __device__ void bonds_gpu(const int       i,
             float3 fij = fbond * dx;
             atomicAdd(&gm_f[ai], fij);
             atomicAdd(&gm_f[aj], -fij);
-            if (calcVir && ki != CENTRAL)
+            if (calcVir && ki != gmx::c_centralShiftIndex)
             {
                 atomicAdd(&sm_fShiftLoc[ki], fij);
-                atomicAdd(&sm_fShiftLoc[CENTRAL], -fij);
+                atomicAdd(&sm_fShiftLoc[gmx::c_centralShiftIndex], -fij);
             }
         }
     }
@@ -228,7 +228,7 @@ __device__ void angles_gpu(const int       i,
             if (calcVir)
             {
                 atomicAdd(&sm_fShiftLoc[t1], f_i);
-                atomicAdd(&sm_fShiftLoc[CENTRAL], f_j);
+                atomicAdd(&sm_fShiftLoc[gmx::c_centralShiftIndex], f_j);
                 atomicAdd(&sm_fShiftLoc[t2], f_k);
             }
         }
@@ -310,7 +310,7 @@ __device__ void urey_bradley_gpu(const int       i,
             if (calcVir)
             {
                 atomicAdd(&sm_fShiftLoc[t1], f_i);
-                atomicAdd(&sm_fShiftLoc[CENTRAL], f_j);
+                atomicAdd(&sm_fShiftLoc[gmx::c_centralShiftIndex], f_j);
                 atomicAdd(&sm_fShiftLoc[t2], f_k);
             }
         }
@@ -329,10 +329,10 @@ __device__ void urey_bradley_gpu(const int       i,
             atomicAdd(&gm_f[ai], fik);
             atomicAdd(&gm_f[ak], -fik);
 
-            if (calcVir && ki != CENTRAL)
+            if (calcVir && ki != gmx::c_centralShiftIndex)
             {
                 atomicAdd(&sm_fShiftLoc[ki], fik);
-                atomicAdd(&sm_fShiftLoc[CENTRAL], -fik);
+                atomicAdd(&sm_fShiftLoc[gmx::c_centralShiftIndex], -fik);
             }
         }
     }
@@ -432,7 +432,7 @@ __device__ static void do_dih_fup_gpu(const int      i,
             int    t3 = pbcDxAiuc<calcVir>(pbcAiuc, gm_xq[l], gm_xq[j], dx_jl);
 
             atomicAdd(&sm_fShiftLoc[t1], f_i);
-            atomicAdd(&sm_fShiftLoc[CENTRAL], -f_j);
+            atomicAdd(&sm_fShiftLoc[gmx::c_centralShiftIndex], -f_j);
             atomicAdd(&sm_fShiftLoc[t2], -f_k);
             atomicAdd(&sm_fShiftLoc[t3], f_l);
         }
@@ -709,10 +709,10 @@ __device__ void pairs_gpu(const int       i,
         /* Add the forces */
         atomicAdd(&gm_f[ai], f);
         atomicAdd(&gm_f[aj], -f);
-        if (calcVir && fshift_index != CENTRAL)
+        if (calcVir && fshift_index != gmx::c_centralShiftIndex)
         {
             atomicAdd(&sm_fShiftLoc[fshift_index], f);
-            atomicAdd(&sm_fShiftLoc[CENTRAL], -f);
+            atomicAdd(&sm_fShiftLoc[gmx::c_centralShiftIndex], -f);
         }
 
         if (calcEner)
@@ -738,11 +738,11 @@ __global__ void exec_kernel_gpu(BondedCudaKernelParameters kernelParams)
     extern __shared__ char sm_dynamicShmem[];
     char*                  sm_nextSlotPtr = sm_dynamicShmem;
     float3*                sm_fShiftLoc   = (float3*)sm_nextSlotPtr;
-    sm_nextSlotPtr += SHIFTS * sizeof(float3);
+    sm_nextSlotPtr += c_numShiftVectors * sizeof(float3);
 
     if (calcVir)
     {
-        if (threadIdx.x < SHIFTS)
+        if (threadIdx.x < c_numShiftVectors)
         {
             sm_fShiftLoc[threadIdx.x] = make_float3(0.0f, 0.0f, 0.0f);
         }
@@ -893,11 +893,11 @@ __global__ void exec_kernel_gpu(BondedCudaKernelParameters kernelParams)
             atomicAdd(vtotElec, sm_vTotElec[warpId]);
         }
     }
-    /* Accumulate shift vectors from shared memory to global memory on the first SHIFTS threads of the block. */
+    /* Accumulate shift vectors from shared memory to global memory on the first c_numShiftVectors threads of the block. */
     if (calcVir)
     {
         __syncthreads();
-        if (threadIdx.x < SHIFTS)
+        if (threadIdx.x < c_numShiftVectors)
         {
             atomicAdd(kernelParams.d_fShift[threadIdx.x], sm_fShiftLoc[threadIdx.x]);
         }
diff --git a/src/gromacs/listed_forces/listed_forces.cpp b/src/gromacs/listed_forces/listed_forces.cpp
index 82b93fa3e2..3a226ed961 100644
--- a/src/gromacs/listed_forces/listed_forces.cpp
+++ b/src/gromacs/listed_forces/listed_forces.cpp
@@ -179,7 +179,7 @@ void ListedForces::setup(const InteractionDefinitions& domainIdef, const int num
     if (idef_->ilsort == ilsortFE_SORTED)
     {
         forceBufferLambda_.resize(numAtomsForce * sizeof(rvec4) / sizeof(real));
-        shiftForceBufferLambda_.resize(SHIFTS);
+        shiftForceBufferLambda_.resize(gmx::c_numShiftVectors);
     }
 }
 
@@ -214,7 +214,7 @@ void zero_thread_output(f_thread_t* f_t)
         }
     }
 
-    for (int i = 0; i < SHIFTS; i++)
+    for (int i = 0; i < gmx::c_numShiftVectors; i++)
     {
         clear_rvec(f_t->fshift[i]);
     }
@@ -322,7 +322,7 @@ void reduce_thread_output(gmx::ForceWithShiftForces* forceWithShiftForces,
 
         if (stepWork.computeVirial)
         {
-            for (int i = 0; i < SHIFTS; i++)
+            for (int i = 0; i < gmx::c_numShiftVectors; i++)
             {
                 for (int t = 1; t < bt->nthreads; t++)
                 {
diff --git a/src/gromacs/listed_forces/listed_internal.h b/src/gromacs/listed_forces/listed_internal.h
index 5013a6917a..9366efcce1 100644
--- a/src/gromacs/listed_forces/listed_internal.h
+++ b/src/gromacs/listed_forces/listed_internal.h
@@ -111,7 +111,7 @@ struct f_thread_t
     //! Index to touched blocks
     std::vector<int> block_index;
 
-    //! Shift force array, size SHIFTS
+    //! Shift force array, size c_numShiftVectors
     std::vector<gmx::RVec> fshift;
     //! Energy array
     real ener[F_NRE];
diff --git a/src/gromacs/listed_forces/manage_threading.cpp b/src/gromacs/listed_forces/manage_threading.cpp
index db95e75c03..0eeac68641 100644
--- a/src/gromacs/listed_forces/manage_threading.cpp
+++ b/src/gromacs/listed_forces/manage_threading.cpp
@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 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.
@@ -490,7 +490,9 @@ void setup_bonded_threading(bonded_threading_t*           bt,
     }
 }
 
-f_thread_t::f_thread_t(int numEnergyGroups) : fshift(SHIFTS), grpp(numEnergyGroups) {}
+f_thread_t::f_thread_t(int numEnergyGroups) : fshift(gmx::c_numShiftVectors), grpp(numEnergyGroups)
+{
+}
 
 bonded_threading_t::bonded_threading_t(const int numThreads, const int numEnergyGroups, FILE* fplog) :
     nthreads(numThreads),
diff --git a/src/gromacs/listed_forces/orires.cpp b/src/gromacs/listed_forces/orires.cpp
index 6b41b649ea..8adb3d20a5 100644
--- a/src/gromacs/listed_forces/orires.cpp
+++ b/src/gromacs/listed_forces/orires.cpp
@@ -663,7 +663,7 @@ real orires(int             nfa,
             t_oriresdata*            oriresdata,
             int gmx_unused* global_atom_index)
 {
-    int      ex, power, ki = CENTRAL;
+    int      ex, power, ki = gmx::c_centralShiftIndex;
     real     r2, invr, invr2, fc, smooth_fc, dev, devins, pfac;
     rvec     r, Sr, fij;
     real     vtot;
@@ -745,7 +745,7 @@ real orires(int             nfa,
                 if (fshift)
                 {
                     fshift[ki][i] += fij[i];
-                    fshift[CENTRAL][i] -= fij[i];
+                    fshift[gmx::c_centralShiftIndex][i] -= fij[i];
                 }
             }
         }
diff --git a/src/gromacs/listed_forces/pairs.cpp b/src/gromacs/listed_forces/pairs.cpp
index 43d292b6fc..0f4bef2714 100644
--- a/src/gromacs/listed_forces/pairs.cpp
+++ b/src/gromacs/listed_forces/pairs.cpp
@@ -483,7 +483,7 @@ static real do_pairs_general(int                 ftype,
         }
         else
         {
-            fshift_index = CENTRAL;
+            fshift_index = c_centralShiftIndex;
             rvec_sub(x[ai], x[aj], dx);
         }
         r2 = norm2(dx);
@@ -553,10 +553,10 @@ static real do_pairs_general(int                 ftype,
 
         if (computeVirial(flavor))
         {
-            if (fshift_index != CENTRAL)
+            if (fshift_index != c_centralShiftIndex)
             {
                 rvec_inc(fshift[fshift_index], dx);
-                rvec_dec(fshift[CENTRAL], dx);
+                rvec_dec(fshift[c_centralShiftIndex], dx);
             }
         }
     }
diff --git a/src/gromacs/listed_forces/tests/bonded.cpp b/src/gromacs/listed_forces/tests/bonded.cpp
index 6299f3384a..76aa9e97c7 100644
--- a/src/gromacs/listed_forces/tests/bonded.cpp
+++ b/src/gromacs/listed_forces/tests/bonded.cpp
@@ -92,7 +92,7 @@ struct OutputQuantities
     //! Derivative with respect to lambda
     real dvdlambda = 0;
     //! Shift vectors
-    rvec fshift[N_IVEC] = { { 0 } };
+    rvec fshift[c_numShiftVectors] = { { 0 } };
     //! Forces
     alignas(GMX_REAL_MAX_SIMD_WIDTH * sizeof(real)) rvec4 f[c_numAtoms] = { { 0 } };
 };
@@ -616,7 +616,7 @@ protected:
             if (computeVirial(flavor))
             {
                 shiftForcesChecker.setDefaultTolerance(shiftForcesTolerance_);
-                shiftForcesChecker.checkVector(output.fshift[CENTRAL], "Central");
+                shiftForcesChecker.checkVector(output.fshift[c_centralShiftIndex], "Central");
             }
             else
             {
diff --git a/src/gromacs/listed_forces/tests/pairs.cpp b/src/gromacs/listed_forces/tests/pairs.cpp
index d9087d9400..5b038f1240 100644
--- a/src/gromacs/listed_forces/tests/pairs.cpp
+++ b/src/gromacs/listed_forces/tests/pairs.cpp
@@ -111,7 +111,7 @@ struct OutputQuantities
     //! Derivative with respect to lambda
     std::vector<real> dvdLambda;
     //! Shift vectors
-    rvec fShift[N_IVEC] = { { 0 } };
+    rvec fShift[gmx::detail::c_numIvecs] = { { 0 } };
     //! Forces
     alignas(GMX_REAL_MAX_SIMD_WIDTH * sizeof(real)) rvec4 f[c_numAtoms] = { { 0 } };
 };
@@ -393,7 +393,7 @@ protected:
 
             if (computeVirial(flavor))
             {
-                shiftForcesChecker.checkVector(output.fShift[CENTRAL], "Central");
+                shiftForcesChecker.checkVector(output.fShift[gmx::c_centralShiftIndex], "Central");
             }
             else
             {
diff --git a/src/gromacs/mdlib/calcvir.cpp b/src/gromacs/mdlib/calcvir.cpp
index 0cf7c41952..42076ab292 100644
--- a/src/gromacs/mdlib/calcvir.cpp
+++ b/src/gromacs/mdlib/calcvir.cpp
@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2013,2014,2015,2016,2018 by the GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 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.
@@ -84,7 +84,7 @@ static void calc_x_times_f(int nxf, const rvec x[], const rvec f[], gmx_bool bSc
 
         if (bScrewPBC)
         {
-            int isx = IS2X(i);
+            int isx = gmx::shiftIndexToXDim(i);
             /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
             if (isx == 1 || isx == -1)
             {
diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp
index d2a3bd8692..7628ce37d1 100644
--- a/src/gromacs/mdlib/forcerec.cpp
+++ b/src/gromacs/mdlib/forcerec.cpp
@@ -105,7 +105,7 @@
 ForceHelperBuffers::ForceHelperBuffers(bool haveDirectVirialContributions) :
     haveDirectVirialContributions_(haveDirectVirialContributions)
 {
-    shiftForces_.resize(SHIFTS);
+    shiftForces_.resize(gmx::c_numShiftVectors);
 }
 
 void ForceHelperBuffers::resize(int numAtoms)
@@ -1217,7 +1217,7 @@ void init_forcerec(FILE*                            fplog,
 
     if (forcerec->shift_vec.empty())
     {
-        forcerec->shift_vec.resize(SHIFTS);
+        forcerec->shift_vec.resize(gmx::c_numShiftVectors);
     }
 
     if (forcerec->nbfp.empty())
diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp
index 9f6711386b..60a046fec2 100644
--- a/src/gromacs/mdlib/sim_util.cpp
+++ b/src/gromacs/mdlib/sim_util.cpp
@@ -167,8 +167,8 @@ static void calc_virial(int                              start,
     /* The short-range virial from surrounding boxes */
     const rvec* fshift          = as_rvec_array(forceWithShiftForces.shiftForces().data());
     const rvec* shiftVecPointer = as_rvec_array(fr->shift_vec.data());
-    calc_vir(SHIFTS, shiftVecPointer, fshift, vir_part, pbcType == PbcType::Screw, box);
-    inc_nrnb(nrnb, eNR_VIRIAL, SHIFTS);
+    calc_vir(gmx::c_numShiftVectors, shiftVecPointer, fshift, vir_part, pbcType == PbcType::Screw, box);
+    inc_nrnb(nrnb, eNR_VIRIAL, gmx::c_numShiftVectors);
 
     /* Calculate partial virial, for local atoms only, based on short range.
      * Total virial is computed in global_stat, called from do_md
diff --git a/src/gromacs/mdlib/vsite.cpp b/src/gromacs/mdlib/vsite.cpp
index adcf0ea8d6..6dd3c6f104 100644
--- a/src/gromacs/mdlib/vsite.cpp
+++ b/src/gromacs/mdlib/vsite.cpp
@@ -165,7 +165,7 @@ struct VsiteThread
     //! The interaction lists, only vsite entries are used
     std::array<InteractionList, F_NRE> ilist;
     //! Local fshift accumulation buffer
-    std::array<RVec, SHIFTS> fshift;
+    std::array<RVec, c_numShiftVectors> fshift;
     //! Local virial dx*df accumulation buffer
     matrix dxdf;
     //! Tells if interdependent task idTask should be used (in addition to the rest of this task), this bool has the same value on all threads
@@ -313,7 +313,7 @@ static int pbc_rvec_sub(const t_pbc* pbc, const rvec xi, const rvec xj, rvec dx)
     else
     {
         rvec_sub(xi, xj, dx);
-        return CENTRAL;
+        return c_centralShiftIndex;
     }
 }
 
@@ -1119,7 +1119,7 @@ static void construct_vsites_thread(ArrayRef<RVec>                  x,
                     /* Keep the vsite in the same periodic image as before */
                     rvec dx;
                     int  ishift = pbc_dx_aiuc(pbc_null, x[avsite], xv, dx);
-                    if (ishift != CENTRAL)
+                    if (ishift != c_centralShiftIndex)
                     {
                         rvec_add(xv, dx, x[avsite]);
                     }
@@ -1306,14 +1306,14 @@ static void spread_vsite2(const t_iatom        ia[],
         }
         else
         {
-            siv = CENTRAL;
-            sij = CENTRAL;
+            siv = c_centralShiftIndex;
+            sij = c_centralShiftIndex;
         }
 
-        if (siv != CENTRAL || sij != CENTRAL)
+        if (siv != c_centralShiftIndex || sij != c_centralShiftIndex)
         {
             rvec_inc(fshift[siv], f[av]);
-            rvec_dec(fshift[CENTRAL], fi);
+            rvec_dec(fshift[c_centralShiftIndex], fi);
             rvec_dec(fshift[sij], fj);
         }
     }
@@ -1396,15 +1396,15 @@ static void spread_vsite2FD(const t_iatom        ia[],
         }
         else
         {
-            svi = CENTRAL;
+            svi = c_centralShiftIndex;
         }
 
-        if (svi != CENTRAL || sji != CENTRAL)
+        if (svi != c_centralShiftIndex || sji != c_centralShiftIndex)
         {
             rvec_dec(fshift[svi], fv);
-            fshift[CENTRAL][XX] += fv[XX] - fj[XX];
-            fshift[CENTRAL][YY] += fv[YY] - fj[YY];
-            fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ];
+            fshift[c_centralShiftIndex][XX] += fv[XX] - fj[XX];
+            fshift[c_centralShiftIndex][YY] += fv[YY] - fj[YY];
+            fshift[c_centralShiftIndex][ZZ] += fv[ZZ] - fj[ZZ];
             fshift[sji][XX] += fj[XX];
             fshift[sji][YY] += fj[YY];
             fshift[sji][ZZ] += fj[ZZ];
@@ -1478,15 +1478,15 @@ static void spread_vsite3(const t_iatom        ia[],
         }
         else
         {
-            siv = CENTRAL;
-            sij = CENTRAL;
-            sik = CENTRAL;
+            siv = c_centralShiftIndex;
+            sij = c_centralShiftIndex;
+            sik = c_centralShiftIndex;
         }
 
-        if (siv != CENTRAL || sij != CENTRAL || sik != CENTRAL)
+        if (siv != c_centralShiftIndex || sij != c_centralShiftIndex || sik != c_centralShiftIndex)
         {
             rvec_inc(fshift[siv], f[av]);
-            rvec_dec(fshift[CENTRAL], fi);
+            rvec_dec(fshift[c_centralShiftIndex], fi);
             rvec_dec(fshift[sij], fj);
             rvec_dec(fshift[sik], fk);
         }
@@ -1561,15 +1561,15 @@ static void spread_vsite3FD(const t_iatom        ia[],
         }
         else
         {
-            svi = CENTRAL;
+            svi = c_centralShiftIndex;
         }
 
-        if (svi != CENTRAL || sji != CENTRAL || skj != CENTRAL)
+        if (svi != c_centralShiftIndex || sji != c_centralShiftIndex || skj != c_centralShiftIndex)
         {
             rvec_dec(fshift[svi], fv);
-            fshift[CENTRAL][XX] += fv[XX] - (1 + a) * temp[XX];
-            fshift[CENTRAL][YY] += fv[YY] - (1 + a) * temp[YY];
-            fshift[CENTRAL][ZZ] += fv[ZZ] - (1 + a) * temp[ZZ];
+            fshift[c_centralShiftIndex][XX] += fv[XX] - (1 + a) * temp[XX];
+            fshift[c_centralShiftIndex][YY] += fv[YY] - (1 + a) * temp[YY];
+            fshift[c_centralShiftIndex][ZZ] += fv[ZZ] - (1 + a) * temp[ZZ];
             fshift[sji][XX] += temp[XX];
             fshift[sji][YY] += temp[YY];
             fshift[sji][ZZ] += temp[ZZ];
@@ -1683,15 +1683,15 @@ static void spread_vsite3FAD(const t_iatom        ia[],
         }
         else
         {
-            svi = CENTRAL;
+            svi = c_centralShiftIndex;
         }
 
-        if (svi != CENTRAL || sji != CENTRAL || skj != CENTRAL)
+        if (svi != c_centralShiftIndex || sji != c_centralShiftIndex || skj != c_centralShiftIndex)
         {
             rvec_dec(fshift[svi], fv);
-            fshift[CENTRAL][XX] += fv[XX] - f1[XX] - (1 - c1) * f2[XX] + f3[XX];
-            fshift[CENTRAL][YY] += fv[YY] - f1[YY] - (1 - c1) * f2[YY] + f3[YY];
-            fshift[CENTRAL][ZZ] += fv[ZZ] - f1[ZZ] - (1 - c1) * f2[ZZ] + f3[ZZ];
+            fshift[c_centralShiftIndex][XX] += fv[XX] - f1[XX] - (1 - c1) * f2[XX] + f3[XX];
+            fshift[c_centralShiftIndex][YY] += fv[YY] - f1[YY] - (1 - c1) * f2[YY] + f3[YY];
+            fshift[c_centralShiftIndex][ZZ] += fv[ZZ] - f1[ZZ] - (1 - c1) * f2[ZZ] + f3[ZZ];
             fshift[sji][XX] += f1[XX] - c1 * f2[XX] - f3[XX];
             fshift[sji][YY] += f1[YY] - c1 * f2[YY] - f3[YY];
             fshift[sji][ZZ] += f1[ZZ] - c1 * f2[ZZ] - f3[ZZ];
@@ -1777,15 +1777,15 @@ static void spread_vsite3OUT(const t_iatom        ia[],
         }
         else
         {
-            svi = CENTRAL;
+            svi = c_centralShiftIndex;
         }
 
-        if (svi != CENTRAL || sji != CENTRAL || ski != CENTRAL)
+        if (svi != c_centralShiftIndex || sji != c_centralShiftIndex || ski != c_centralShiftIndex)
         {
             rvec_dec(fshift[svi], fv);
-            fshift[CENTRAL][XX] += fv[XX] - fj[XX] - fk[XX];
-            fshift[CENTRAL][YY] += fv[YY] - fj[YY] - fk[YY];
-            fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ];
+            fshift[c_centralShiftIndex][XX] += fv[XX] - fj[XX] - fk[XX];
+            fshift[c_centralShiftIndex][YY] += fv[YY] - fj[YY] - fk[YY];
+            fshift[c_centralShiftIndex][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ];
             rvec_inc(fshift[sji], fj);
             rvec_inc(fshift[ski], fk);
         }
@@ -1879,15 +1879,16 @@ static void spread_vsite4FD(const t_iatom        ia[],
         }
         else
         {
-            svi = CENTRAL;
+            svi = c_centralShiftIndex;
         }
 
-        if (svi != CENTRAL || sji != CENTRAL || skj != CENTRAL || slj != CENTRAL)
+        if (svi != c_centralShiftIndex || sji != c_centralShiftIndex || skj != c_centralShiftIndex
+            || slj != c_centralShiftIndex)
         {
             rvec_dec(fshift[svi], fv);
             for (m = 0; m < DIM; m++)
             {
-                fshift[CENTRAL][m] += fv[m] - (1 + a + b) * temp[m];
+                fshift[c_centralShiftIndex][m] += fv[m] - (1 + a + b) * temp[m];
                 fshift[sji][m] += temp[m];
                 fshift[skj][m] += a * temp[m];
                 fshift[slj][m] += b * temp[m];
@@ -2038,15 +2039,16 @@ static void spread_vsite4FDN(const t_iatom        ia[],
         }
         else
         {
-            svi = CENTRAL;
+            svi = c_centralShiftIndex;
         }
 
-        if (svi != CENTRAL || sij != CENTRAL || sik != CENTRAL || sil != CENTRAL)
+        if (svi != c_centralShiftIndex || sij != c_centralShiftIndex || sik != c_centralShiftIndex
+            || sil != c_centralShiftIndex)
         {
             rvec_dec(fshift[svi], fv);
-            fshift[CENTRAL][XX] += fv[XX] - fj[XX] - fk[XX] - fl[XX];
-            fshift[CENTRAL][YY] += fv[YY] - fj[YY] - fk[YY] - fl[YY];
-            fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ] - fl[ZZ];
+            fshift[c_centralShiftIndex][XX] += fv[XX] - fj[XX] - fk[XX] - fl[XX];
+            fshift[c_centralShiftIndex][YY] += fv[YY] - fj[YY] - fk[YY] - fl[YY];
+            fshift[c_centralShiftIndex][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ] - fl[ZZ];
             rvec_inc(fshift[sij], fj);
             rvec_inc(fshift[sik], fk);
             rvec_inc(fshift[sil], fl);
@@ -2098,16 +2100,16 @@ static int spread_vsiten(const t_iatom             ia[],
         }
         else
         {
-            siv = CENTRAL;
+            siv = c_centralShiftIndex;
         }
         a = ip[ia[i]].vsiten.a;
         svmul(a, f[av], fi);
         rvec_inc(f[ai], fi);
 
-        if (virialHandling == VirialHandling::Pbc && siv != CENTRAL)
+        if (virialHandling == VirialHandling::Pbc && siv != c_centralShiftIndex)
         {
             rvec_inc(fshift[siv], fi);
-            rvec_dec(fshift[CENTRAL], fi);
+            rvec_dec(fshift[c_centralShiftIndex], fi);
         }
         /* 6 Flops */
     }
@@ -2355,7 +2357,7 @@ void VirtualSitesHandler::Impl::spreadForces(ArrayRef<const RVec> x,
                     {
                         fshift_t = tData.fshift;
 
-                        for (int i = 0; i < SHIFTS; i++)
+                        for (int i = 0; i < c_numShiftVectors; i++)
                         {
                             clear_rvec(fshift_t[i]);
                         }
@@ -2434,7 +2436,7 @@ void VirtualSitesHandler::Impl::spreadForces(ArrayRef<const RVec> x,
         {
             for (int th = 1; th < numThreads; th++)
             {
-                for (int i = 0; i < SHIFTS; i++)
+                for (int i = 0; i < c_numShiftVectors; i++)
                 {
                     rvec_inc(fshift[i], threadingInfo_.threadData(th).fshift[i]);
                 }
diff --git a/src/gromacs/mdtypes/forceoutput.h b/src/gromacs/mdtypes/forceoutput.h
index c31bcf05ba..a9a9ae6b2c 100644
--- a/src/gromacs/mdtypes/forceoutput.h
+++ b/src/gromacs/mdtypes/forceoutput.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 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.
@@ -74,7 +74,7 @@ public:
      *
      * \param[in] force          A force buffer that will be used for storing forces
      * \param[in] computeVirial  True when algorithms are required to provide their virial contribution (for the current force evaluation)
-     * \param[in] shiftForces    A shift forces buffer of size SHIFTS, only used with \p computeVirial = true
+     * \param[in] shiftForces    A shift forces buffer of size c_numShiftVectors, only used with \p computeVirial = true
      */
     ForceWithShiftForces(const gmx::ArrayRefWithPadding<gmx::RVec>& force,
                          const bool                                 computeVirial,
@@ -110,7 +110,7 @@ private:
     gmx::ArrayRefWithPadding<gmx::RVec> force_;
     //! True when virial computation is requested
     bool computeVirial_;
-    //! A buffer for storing the shift forces, size SHIFTS
+    //! A buffer for storing the shift forces, size c_numShiftVectors
     gmx::ArrayRef<gmx::RVec> shiftForces_;
     //! Tells whether we have spread the vsite forces
     bool haveSpreadVsiteForces_ = false;
diff --git a/src/gromacs/mdtypes/forcerec.h b/src/gromacs/mdtypes/forcerec.h
index de1c67659f..414e533a6f 100644
--- a/src/gromacs/mdtypes/forcerec.h
+++ b/src/gromacs/mdtypes/forcerec.h
@@ -149,7 +149,7 @@ public:
         return forceBufferForDirectVirialContributions_;
     }
 
-    //! Returns the buffer for shift forces, size SHIFTS
+    //! Returns the buffer for shift forces, size c_numShiftVectors
     gmx::ArrayRef<gmx::RVec> shiftForces() { return shiftForces_; }
 
     //! Resizes the direct virial contribution buffer, when present
@@ -160,7 +160,7 @@ private:
     bool haveDirectVirialContributions_ = false;
     //! Force buffer for force computation with direct virial contributions
     std::vector<gmx::RVec> forceBufferForDirectVirialContributions_;
-    //! Shift force array for computing the virial, size SHIFTS
+    //! Shift force array for computing the virial, size c_numShiftVectors
     std::vector<gmx::RVec> shiftForces_;
 };
 // NOLINTNEXTLINE (clang-analyzer-optin.performance.Padding)
diff --git a/src/gromacs/nbnxm/atomdata.cpp b/src/gromacs/nbnxm/atomdata.cpp
index 1fa91c6f03..d2e497b6da 100644
--- a/src/gromacs/nbnxm/atomdata.cpp
+++ b/src/gromacs/nbnxm/atomdata.cpp
@@ -107,7 +107,7 @@ nbnxn_atomdata_output_t::nbnxn_atomdata_output_t(Nbnxm::KernelType  kernelType,
     Vvdw({}, { pinningPolicy }),
     Vc({}, { pinningPolicy })
 {
-    fshift.resize(SHIFTS * DIM);
+    fshift.resize(gmx::c_numShiftVectors * DIM);
     Vvdw.resize(numEnergyGroups * numEnergyGroups);
     Vc.resize(numEnergyGroups * numEnergyGroups);
 
@@ -666,7 +666,7 @@ void nbnxn_atomdata_init(const gmx::MDLogger&    mdlog,
         nbat->FFormat = nbatXYZ;
     }
 
-    nbat->shift_vec.resize(SHIFTS);
+    nbat->shift_vec.resize(gmx::c_numShiftVectors);
 
     nbat->xstride = (nbat->XFormat == nbatXYZQ ? STRIDE_XYZQ : DIM);
     nbat->fstride = (nbat->FFormat == nbatXYZQ ? STRIDE_XYZQ : DIM);
@@ -1276,7 +1276,7 @@ void nbnxn_atomdata_add_nbat_fshift_to_fshift(const nbnxn_atomdata_t& nbat, gmx:
 {
     gmx::ArrayRef<const nbnxn_atomdata_output_t> outputBuffers = nbat.out;
 
-    for (int s = 0; s < SHIFTS; s++)
+    for (int s = 0; s < gmx::c_numShiftVectors; s++)
     {
         rvec sum;
         clear_rvec(sum);
diff --git a/src/gromacs/nbnxm/atomdata.h b/src/gromacs/nbnxm/atomdata.h
index d676ae1f4c..cb878aa779 100644
--- a/src/gromacs/nbnxm/atomdata.h
+++ b/src/gromacs/nbnxm/atomdata.h
@@ -123,7 +123,7 @@ struct nbnxn_atomdata_output_t
 
     //! f, size natoms*fstride
     gmx::HostVector<real> f;
-    //! Shift force array, size SHIFTS*DIM
+    //! Shift force array, size c_numShiftVectors*DIM
     gmx::HostVector<real> fshift;
     //! Temporary Van der Waals group energy storage
     gmx::HostVector<real> Vvdw;
diff --git a/src/gromacs/nbnxm/benchmark/bench_system.cpp b/src/gromacs/nbnxm/benchmark/bench_system.cpp
index 3076939ce5..1ba1c01d90 100644
--- a/src/gromacs/nbnxm/benchmark/bench_system.cpp
+++ b/src/gromacs/nbnxm/benchmark/bench_system.cpp
@@ -200,7 +200,7 @@ BenchmarkSystem::BenchmarkSystem(const int multiplicationFactor, const std::stri
 
     forceRec.ntype = numAtomTypes;
     forceRec.nbfp  = nonbondedParameters;
-    forceRec.shift_vec.resize(SHIFTS);
+    forceRec.shift_vec.resize(gmx::c_numShiftVectors);
     calc_shifts(box, forceRec.shift_vec);
     if (!outputFile.empty())
     {
diff --git a/src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh b/src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh
index 49297a2f0a..7df9231d51 100644
--- a/src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh
+++ b/src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh
@@ -326,7 +326,8 @@ __launch_bounds__(THREADS_PER_BLOCK)
     E_el         = 0.0f;
 
 #        ifdef EXCLUSION_FORCES /* Ewald or RF */
-    if (nb_sci.shift == CENTRAL && pl_cj4[cij4_start].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
+    if (nb_sci.shift == gmx::c_centralShiftIndex
+        && pl_cj4[cij4_start].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
     {
         /* we have the diagonal: add the charge and LJ self interaction energy term */
         for (i = 0; i < c_nbnxnGpuNumClusterPerSupercluster; i++)
@@ -364,7 +365,7 @@ __launch_bounds__(THREADS_PER_BLOCK)
 #    endif /* CALC_ENERGIES */
 
 #    ifdef EXCLUSION_FORCES
-    const int nonSelfInteraction = !(nb_sci.shift == CENTRAL & tidxj <= tidxi);
+    const int nonSelfInteraction = !(nb_sci.shift == gmx::c_centralShiftIndex & tidxj <= tidxi);
 #    endif
 
     /* loop over the j clusters = seen by any of the atoms in the current super-cluster;
@@ -627,7 +628,7 @@ __launch_bounds__(THREADS_PER_BLOCK)
     }
 
     /* skip central shifts when summing shift forces */
-    if (nb_sci.shift == CENTRAL)
+    if (nb_sci.shift == gmx::c_centralShiftIndex)
     {
         bCalcFshift = false;
     }
diff --git a/src/gromacs/nbnxm/gpu_common.h b/src/gromacs/nbnxm/gpu_common.h
index e357433f0f..cdc35d093a 100644
--- a/src/gromacs/nbnxm/gpu_common.h
+++ b/src/gromacs/nbnxm/gpu_common.h
@@ -153,7 +153,7 @@ static inline void gpu_reduce_staged_outputs(const NBStagingData&      nbst,
 
         if (reduceFshift)
         {
-            for (int i = 0; i < SHIFTS; i++)
+            for (int i = 0; i < gmx::c_numShiftVectors; i++)
             {
                 rvec_inc(fshift[i], nbst.fShift[i]);
             }
diff --git a/src/gromacs/nbnxm/kernel_common.cpp b/src/gromacs/nbnxm/kernel_common.cpp
index b28b9a0168..bcc60fc1d4 100644
--- a/src/gromacs/nbnxm/kernel_common.cpp
+++ b/src/gromacs/nbnxm/kernel_common.cpp
@@ -98,7 +98,7 @@ void clearForceBuffer(nbnxn_atomdata_t* nbat, int outputIndex)
 
 void clear_fshift(real* fshift)
 {
-    for (int i = 0; i < SHIFTS * DIM; i++)
+    for (int i = 0; i < gmx::c_numShiftVectors * DIM; i++)
     {
         fshift[i] = 0;
     }
diff --git a/src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp b/src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp
index 6002eb5be7..32f70f2cea 100644
--- a/src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp
+++ b/src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp
@@ -119,7 +119,8 @@ void nbnxn_kernel_gpu_ref(const NbnxnPairlistGpu*        nbl,
         real       vctot    = 0;
         real       Vvdwtot  = 0;
 
-        if (nbln.shift == CENTRAL && nbl->cj4[cj4_ind0].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
+        if (nbln.shift == gmx::c_centralShiftIndex
+            && nbl->cj4[cj4_ind0].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
         {
             /* we have the diagonal:
              * add the charge self interaction energy term
@@ -186,7 +187,7 @@ void nbnxn_kernel_gpu_ref(const NbnxnPairlistGpu*        nbl,
                             {
                                 const int ja = cj * c_clSize + jc;
 
-                                if (nbln.shift == CENTRAL && ci == cj && ja <= ia)
+                                if (nbln.shift == gmx::c_centralShiftIndex && ci == cj && ja <= ia)
                                 {
                                     continue;
                                 }
diff --git a/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h b/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h
index 1203b8e4ba..a74d579ec5 100644
--- a/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h
+++ b/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h
@@ -191,7 +191,7 @@ void
         const int cjind1 = ciEntry.cj_ind_end;
         /* Currently only works super-cells equal to sub-cells */
         const int ci    = ciEntry.ci;
-        const int ci_sh = (ish == CENTRAL ? ci : -1);
+        const int ci_sh = (ish == gmx::c_centralShiftIndex ? ci : -1);
 
         /* We have 5 LJ/C combinations, but use only three inner loops,
          * as the other combinations are unlikely and/or not much faster:
diff --git a/src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h b/src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h
index f45fc6f843..be27518f93 100644
--- a/src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h
+++ b/src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h
@@ -340,7 +340,7 @@
         const int cjind0 = ciEntry.cj_ind_start;
         const int cjind1 = ciEntry.cj_ind_end;
         const int ci     = ciEntry.ci;
-        const int ci_sh  = (ish == CENTRAL ? ci : -1);
+        const int ci_sh  = (ish == gmx::c_centralShiftIndex ? ci : -1);
 
         shX_S = SimdReal(shiftvec[ish3]);
         shY_S = SimdReal(shiftvec[ish3 + 1]);
diff --git a/src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h b/src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h
index ea0ec1fbac..1946355f26 100644
--- a/src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h
+++ b/src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h
@@ -361,7 +361,7 @@
         const int cjind0 = ciEntry.cj_ind_start;
         const int cjind1 = ciEntry.cj_ind_end;
         const int ci     = ciEntry.ci;
-        const int ci_sh  = (ish == CENTRAL ? ci : -1);
+        const int ci_sh  = (ish == gmx::c_centralShiftIndex ? ci : -1);
 
         shX_S = SimdReal(shiftvec[ish3]);
         shY_S = SimdReal(shiftvec[ish3 + 1]);
diff --git a/src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp b/src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp
index e1df9d0a41..ce93e61e01 100644
--- a/src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp
+++ b/src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp
@@ -249,14 +249,14 @@ static inline void initAtomdataFirst(NBAtomDataGpu*       atomdata,
                                      const DeviceStream&  localStream)
 {
     atomdata->numTypes = numTypes;
-    allocateDeviceBuffer(&atomdata->shiftVec, SHIFTS, deviceContext);
+    allocateDeviceBuffer(&atomdata->shiftVec, gmx::c_numShiftVectors, deviceContext);
     atomdata->shiftVecUploaded = false;
 
-    allocateDeviceBuffer(&atomdata->fShift, SHIFTS, deviceContext);
+    allocateDeviceBuffer(&atomdata->fShift, gmx::c_numShiftVectors, deviceContext);
     allocateDeviceBuffer(&atomdata->eLJ, 1, deviceContext);
     allocateDeviceBuffer(&atomdata->eElec, 1, deviceContext);
 
-    clearDeviceBufferAsync(&atomdata->fShift, 0, SHIFTS, localStream);
+    clearDeviceBufferAsync(&atomdata->fShift, 0, gmx::c_numShiftVectors, localStream);
     clearDeviceBufferAsync(&atomdata->eElec, 0, 1, localStream);
     clearDeviceBufferAsync(&atomdata->eLJ, 0, 1, localStream);
 
@@ -452,7 +452,7 @@ NbnxmGpu* gpu_init(const gmx::DeviceStreamManager& deviceStreamManager,
     /* init nbst */
     pmalloc(reinterpret_cast<void**>(&nb->nbst.eLJ), sizeof(*nb->nbst.eLJ));
     pmalloc(reinterpret_cast<void**>(&nb->nbst.eElec), sizeof(*nb->nbst.eElec));
-    pmalloc(reinterpret_cast<void**>(&nb->nbst.fShift), SHIFTS * sizeof(*nb->nbst.fShift));
+    pmalloc(reinterpret_cast<void**>(&nb->nbst.fShift), gmx::c_numShiftVectors * sizeof(*nb->nbst.fShift));
 
     init_plist(nb->plist[InteractionLocality::Local]);
 
@@ -518,7 +518,7 @@ void gpu_upload_shiftvec(NbnxmGpu* nb, const nbnxn_atomdata_t* nbatom)
         copyToDeviceBuffer(&adat->shiftVec,
                            gmx::asGenericFloat3Pointer(nbatom->shift_vec),
                            0,
-                           SHIFTS,
+                           gmx::c_numShiftVectors,
                            localStream,
                            GpuApiCallBehavior::Async,
                            nullptr);
@@ -718,7 +718,7 @@ void gpu_clear_outputs(NbnxmGpu* nb, bool computeVirial)
     // Clear shift force array and energies if the outputs were used in the current step
     if (computeVirial)
     {
-        clearDeviceBufferAsync(&adat->fShift, 0, SHIFTS, localStream);
+        clearDeviceBufferAsync(&adat->fShift, 0, gmx::c_numShiftVectors, localStream);
         clearDeviceBufferAsync(&adat->eLJ, 0, 1, localStream);
         clearDeviceBufferAsync(&adat->eElec, 0, 1, localStream);
     }
@@ -859,7 +859,7 @@ void gpu_launch_cpyback(NbnxmGpu*                nb,
             copyFromDeviceBuffer(nb->nbst.fShift,
                                  &adat->fShift,
                                  0,
-                                 SHIFTS,
+                                 gmx::c_numShiftVectors,
                                  deviceStream,
                                  GpuApiCallBehavior::Async,
                                  bDoTime ? timers->xf[atomLocality].nb_d2h.fetchNextEvent() : nullptr);
diff --git a/src/gromacs/nbnxm/opencl/CMakeLists.txt b/src/gromacs/nbnxm/opencl/CMakeLists.txt
index dab11de278..e4e875c2af 100644
--- a/src/gromacs/nbnxm/opencl/CMakeLists.txt
+++ b/src/gromacs/nbnxm/opencl/CMakeLists.txt
@@ -2,7 +2,7 @@
 # This file is part of the GROMACS molecular simulation package.
 #
 # Copyright (c) 2012,2013,2014,2015,2018 by the GROMACS development team.
-# Copyright (c) 2019,2020, by the GROMACS development team, led by
+# Copyright (c) 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.
@@ -82,6 +82,7 @@ foreach(ELEC_DEF IN LISTS ELEC_DEFS)
                 -DNBNXM_MIN_DISTANCE_SQUARED_VALUE_FLOAT=3.82e-07
                 -Dc_nbnxnGpuNumClusterPerSupercluster=8
                 -Dc_nbnxnGpuJgroupSize=4
+                -Dc_centralShiftIndex=22
                 -DIATYPE_SHMEM
                 -c -I ${CMAKE_SOURCE_DIR}/src -std=cl1.2
                 -Weverything  -Wno-conversion -Wno-missing-variable-declarations -Wno-used-but-marked-unused
diff --git a/src/gromacs/nbnxm/opencl/nbnxm_ocl_data_mgmt.cpp b/src/gromacs/nbnxm/opencl/nbnxm_ocl_data_mgmt.cpp
index e9ae6dfc75..ea901eb715 100644
--- a/src/gromacs/nbnxm/opencl/nbnxm_ocl_data_mgmt.cpp
+++ b/src/gromacs/nbnxm/opencl/nbnxm_ocl_data_mgmt.cpp
@@ -66,7 +66,6 @@
 #include "gromacs/nbnxm/nbnxm_gpu.h"
 #include "gromacs/nbnxm/nbnxm_gpu_data_mgmt.h"
 #include "gromacs/nbnxm/pairlistsets.h"
-#include "gromacs/pbcutil/ishift.h"
 #include "gromacs/timing/gpu_timing.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
diff --git a/src/gromacs/nbnxm/opencl/nbnxm_ocl_jit_support.cpp b/src/gromacs/nbnxm/opencl/nbnxm_ocl_jit_support.cpp
index 6a494e0821..dfc2c5a6e2 100644
--- a/src/gromacs/nbnxm/opencl/nbnxm_ocl_jit_support.cpp
+++ b/src/gromacs/nbnxm/opencl/nbnxm_ocl_jit_support.cpp
@@ -2,7 +2,7 @@
  * This file is part of the GROMACS molecular simulation package.
  *
  * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 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.
@@ -199,11 +199,13 @@ void nbnxn_gpu_compile_kernels(NbnxmGpu* nb)
                 " -DNBNXM_MIN_DISTANCE_SQUARED_VALUE_FLOAT=%g"
                 " -Dc_nbnxnGpuNumClusterPerSupercluster=%d"
                 " -Dc_nbnxnGpuJgroupSize=%d"
+                " -Dc_centralShiftIndex=%d"
                 "%s",
                 c_nbnxnGpuClusterSize,
                 c_nbnxnMinDistanceSquared,
                 c_nbnxnGpuNumClusterPerSupercluster,
                 c_nbnxnGpuJgroupSize,
+                gmx::c_centralShiftIndex,
                 (nb->bPrefetchLjParam) ? " -DIATYPE_SHMEM" : "");
         try
         {
diff --git a/src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel.clh b/src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel.clh
index e51957eae6..68c1a9df7a 100644
--- a/src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel.clh
+++ b/src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel.clh
@@ -276,7 +276,8 @@ __kernel void NB_KERNEL_FUNC_NAME(nbnxn_kernel, _F_opencl)
     float E_el = 0.0F;
 
 #    if defined EXCLUSION_FORCES /* Ewald or RF */
-    if (nb_sci.shift == CENTRAL && pl_cj4[cij4_start].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
+    if (nb_sci.shift == c_centralShiftIndex
+        && pl_cj4[cij4_start].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
     {
         /* we have the diagonal: add the charge and LJ self interaction energy term */
         for (int i = 0; i < c_nbnxnGpuNumClusterPerSupercluster; i++)
@@ -312,7 +313,7 @@ __kernel void NB_KERNEL_FUNC_NAME(nbnxn_kernel, _F_opencl)
 #endif /* CALC_ENERGIES */
 
 #ifdef EXCLUSION_FORCES
-    const int nonSelfInteraction = !(nb_sci.shift == CENTRAL & tidxj <= tidxi);
+    const int nonSelfInteraction = !(nb_sci.shift == c_centralShiftIndex & tidxj <= tidxi);
 #endif
 
     /* loop over the j clusters = seen by any of the atoms in the current super-cluster */
@@ -604,7 +605,7 @@ __kernel void NB_KERNEL_FUNC_NAME(nbnxn_kernel, _F_opencl)
     }
 
     /* skip central shifts when summing shift forces */
-    if (nb_sci.shift == CENTRAL)
+    if (nb_sci.shift == c_centralShiftIndex)
     {
         bCalcFshift = 0;
     }
diff --git a/src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel_utils.clh b/src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel_utils.clh
index 6fafa99fb4..4e2e208d40 100644
--- a/src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel_utils.clh
+++ b/src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel_utils.clh
@@ -47,7 +47,6 @@
 
 #include "gromacs/gpu_utils/device_utils.clh"
 #include "gromacs/gpu_utils/vectype_ops.clh"
-#include "gromacs/pbcutil/ishift.h"
 
 #include "nbnxm_ocl_consts.h"
 
diff --git a/src/gromacs/nbnxm/pairlist.cpp b/src/gromacs/nbnxm/pairlist.cpp
index 0f220d3d91..5f9ea05056 100644
--- a/src/gromacs/nbnxm/pairlist.cpp
+++ b/src/gromacs/nbnxm/pairlist.cpp
@@ -743,8 +743,8 @@ static void print_nblist_statistics(FILE*                   fp,
             "nbl average j cell list length %.1f\n",
             0.25 * nbl.ncjInUse / std::max(static_cast<double>(nbl.ci.size()), 1.0));
 
-    int cs[SHIFTS] = { 0 };
-    int npexcl     = 0;
+    int cs[gmx::c_numShiftVectors] = { 0 };
+    int npexcl                     = 0;
     for (const nbnxn_ci_t& ciEntry : nbl.ci)
     {
         cs[ciEntry.shift & NBNXN_CI_SHIFT] += ciEntry.cj_ind_end - ciEntry.cj_ind_start;
@@ -761,7 +761,7 @@ static void print_nblist_statistics(FILE*                   fp,
             nbl.cj.size(),
             npexcl,
             100 * npexcl / std::max(static_cast<double>(nbl.cj.size()), 1.0));
-    for (int s = 0; s < SHIFTS; s++)
+    for (int s = 0; s < gmx::c_numShiftVectors; s++)
     {
         if (cs[s] > 0)
         {
@@ -3261,11 +3261,11 @@ static void nbnxn_make_pairlist_part(const Nbnxm::GridSet&   gridSet,
 
                 for (int tx = -shp[XX]; tx <= shp[XX]; tx++)
                 {
-                    const int shift = XYZ2IS(tx, ty, tz);
+                    const int shift = xyzToShiftIndex(tx, ty, tz);
 
-                    const bool excludeSubDiagonal = (isIntraGridList && shift == CENTRAL);
+                    const bool excludeSubDiagonal = (isIntraGridList && shift == gmx::c_centralShiftIndex);
 
-                    if (c_pbcShiftBackward && isIntraGridList && shift > CENTRAL)
+                    if (c_pbcShiftBackward && isIntraGridList && shift > gmx::c_centralShiftIndex)
                     {
                         continue;
                     }
@@ -3327,10 +3327,10 @@ static void nbnxn_make_pairlist_part(const Nbnxm::GridSet&   gridSet,
                         /* When true, leave the pairs with i > j.
                          * Skip half of y when i and j have the same x.
                          */
-                        const bool skipHalfY =
-                                (isIntraGridList && cx == 0
-                                 && (!c_pbcShiftBackward || shift == CENTRAL) && cyf < ci_y);
-                        const int cyf_x = skipHalfY ? ci_y : cyf;
+                        const bool skipHalfY = (isIntraGridList && cx == 0
+                                                && (!c_pbcShiftBackward || shift == gmx::c_centralShiftIndex)
+                                                && cyf < ci_y);
+                        const int  cyf_x     = skipHalfY ? ci_y : cyf;
 
                         for (int cy = cyf_x; cy <= cyl; cy++)
                         {
@@ -3436,7 +3436,7 @@ static void nbnxn_make_pairlist_part(const Nbnxm::GridSet&   gridSet,
                                     /* We want each atom/cell pair only once,
                                      * only use cj >= ci.
                                      */
-                                    if (!c_pbcShiftBackward || shift == CENTRAL)
+                                    if (!c_pbcShiftBackward || shift == gmx::c_centralShiftIndex)
                                     {
                                         firstCell = std::max(firstCell, ci);
                                     }
diff --git a/src/gromacs/nbnxm/sycl/nbnxm_sycl_kernel.cpp b/src/gromacs/nbnxm/sycl/nbnxm_sycl_kernel.cpp
index 1916010e23..e2d4e151fa 100644
--- a/src/gromacs/nbnxm/sycl/nbnxm_sycl_kernel.cpp
+++ b/src/gromacs/nbnxm/sycl/nbnxm_sycl_kernel.cpp
@@ -612,7 +612,8 @@ auto nbnxmKernel(cl::sycl::handler&                                   cgh,
         }
         if constexpr (doCalcEnergies && doExclusionForces)
         {
-            if (nbSci.shift == CENTRAL && a_plistCJ4[cij4Start].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
+            if (nbSci.shift == gmx::c_centralShiftIndex
+                && a_plistCJ4[cij4Start].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
             {
                 // we have the diagonal: add the charge and LJ self interaction energy term
                 for (int i = 0; i < c_nbnxnGpuNumClusterPerSupercluster; i++)
@@ -648,11 +649,11 @@ auto nbnxmKernel(cl::sycl::handler&                                   cgh,
                     energyElec /= epsFac * c_clSize;
                     energyElec *= -ewaldBeta * c_OneOverSqrtPi; /* last factor 1/sqrt(pi) */
                 }
-            } // (nbSci.shift == CENTRAL && a_plistCJ4[cij4Start].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
+            } // (nbSci.shift == gmx::c_centralShiftIndex && a_plistCJ4[cij4Start].cj[0] == sci * c_nbnxnGpuNumClusterPerSupercluster)
         }     // (doCalcEnergies && doExclusionForces)
 
         // Only needed if (doExclusionForces)
-        const bool nonSelfInteraction = !(nbSci.shift == CENTRAL & tidxj <= tidxi);
+        const bool nonSelfInteraction = !(nbSci.shift == gmx::c_centralShiftIndex & tidxj <= tidxi);
 
         // loop over the j clusters = seen by any of the atoms in the current super-cluster
         for (int j4 = cij4Start + tidxz; j4 < cij4End; j4 += 1)
@@ -915,7 +916,7 @@ auto nbnxmKernel(cl::sycl::handler&                                   cgh,
         } // for (int j4 = cij4Start; j4 < cij4End; j4 += 1)
 
         /* skip central shifts when summing shift forces */
-        const bool doCalcShift = (calcShift && !(nbSci.shift == CENTRAL));
+        const bool doCalcShift = (calcShift && !(nbSci.shift == gmx::c_centralShiftIndex));
 
         reduceForceIAndFShift(
                 sm_reductionBuffer, fCiBuf, doCalcShift, itemIdx, tidxi, tidxj, sci, nbSci.shift, a_f, a_fShift);
diff --git a/src/gromacs/pbcutil/ishift.h b/src/gromacs/pbcutil/ishift.h
index 8374337437..7a25fd39c2 100644
--- a/src/gromacs/pbcutil/ishift.h
+++ b/src/gromacs/pbcutil/ishift.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2010,2014,2019, by the GROMACS development team, led by
+ * Copyright (c) 2010,2014,2019,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.
@@ -37,20 +37,42 @@
 #ifndef GMX_PBCUTIL_ISHIFT_H
 #define GMX_PBCUTIL_ISHIFT_H
 
-#define D_BOX_Z 1
-#define D_BOX_Y 1
-#define D_BOX_X 2
-#define N_BOX_Z (2 * D_BOX_Z + 1)
-#define N_BOX_Y (2 * D_BOX_Y + 1)
-#define N_BOX_X (2 * D_BOX_X + 1)
-#define N_IVEC (N_BOX_Z * N_BOX_Y * N_BOX_X)
-#define CENTRAL (N_IVEC / 2)
-#define SHIFTS N_IVEC
+namespace gmx
+{
+//! Maximum dimensions of grid expressing shifts across PBC
+//! \{
+constexpr int c_dBoxZ = 1;
+constexpr int c_dBoxY = 1;
+constexpr int c_dBoxX = 2;
+//! \}
+namespace detail
+{
+constexpr int c_nBoxZ    = 2 * gmx::c_dBoxZ + 1;
+constexpr int c_nBoxY    = 2 * gmx::c_dBoxY + 1;
+constexpr int c_nBoxX    = 2 * gmx::c_dBoxX + 1;
+constexpr int c_numIvecs = detail::c_nBoxZ * detail::c_nBoxY * detail::c_nBoxX;
+} // namespace detail
 
-#define XYZ2IS(x, y, z) (N_BOX_X * (N_BOX_Y * ((z) + D_BOX_Z) + (y) + D_BOX_Y) + (x) + D_BOX_X)
-#define IVEC2IS(iv) (XYZ2IS((iv)[XX], (iv)[YY], (iv)[ZZ]))
-#define IS2X(iv) (((iv) % N_BOX_X) - D_BOX_X)
-#define IS2Y(iv) ((((iv) / N_BOX_X) % N_BOX_Y) - D_BOX_Y)
-#define IS2Z(iv) ((iv) / (N_BOX_X * N_BOX_Y) - D_BOX_Z)
+constexpr int c_centralShiftIndex = detail::c_numIvecs / 2;
+constexpr int c_numShiftVectors   = detail::c_numIvecs;
 
+//! Convert grid coordinates to shift index
+static inline int xyzToShiftIndex(int x, int y, int z)
+{
+    return (detail::c_nBoxX * (detail::c_nBoxY * ((z) + gmx::c_dBoxZ) + (y) + gmx::c_dBoxY) + (x)
+            + gmx::c_dBoxX);
+}
+
+//! Convert grid coordinates to shift index
+static inline int ivecToShiftIndex(ivec iv)
+{
+    return (xyzToShiftIndex((iv)[XX], (iv)[YY], (iv)[ZZ]));
+}
+
+//! Return the shift in the X dimension of grid space corresponding to \c iv
+static inline int shiftIndexToXDim(int iv)
+{
+    return (((iv) % detail::c_nBoxX) - gmx::c_dBoxX);
+}
+} // namespace gmx
 #endif
diff --git a/src/gromacs/pbcutil/pbc.cpp b/src/gromacs/pbcutil/pbc.cpp
index 900663ad0c..e490aaedb1 100644
--- a/src/gromacs/pbcutil/pbc.cpp
+++ b/src/gromacs/pbcutil/pbc.cpp
@@ -1066,10 +1066,10 @@ int pbc_dx_aiuc(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx)
                       "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
     }
 
-    is = IVEC2IS(ishift);
+    is = gmx::ivecToShiftIndex(ishift);
     if (debug)
     {
-        range_check_mesg(is, 0, SHIFTS, "PBC shift vector index range check.");
+        range_check_mesg(is, 0, gmx::c_numShiftVectors, "PBC shift vector index range check.");
     }
 
     return is;
@@ -1198,11 +1198,11 @@ void pbc_dx_d(const t_pbc* pbc, const dvec x1, const dvec x2, dvec dx)
 
 void calc_shifts(const matrix box, gmx::ArrayRef<gmx::RVec> shift_vec)
 {
-    for (int n = 0, m = -D_BOX_Z; m <= D_BOX_Z; m++)
+    for (int n = 0, m = -gmx::c_dBoxZ; m <= gmx::c_dBoxZ; m++)
     {
-        for (int l = -D_BOX_Y; l <= D_BOX_Y; l++)
+        for (int l = -gmx::c_dBoxY; l <= gmx::c_dBoxY; l++)
         {
-            for (int k = -D_BOX_X; k <= D_BOX_X; k++, n++)
+            for (int k = -gmx::c_dBoxX; k <= gmx::c_dBoxX; k++, n++)
             {
                 for (int d = 0; d < DIM; d++)
                 {
diff --git a/src/gromacs/pbcutil/pbc.h b/src/gromacs/pbcutil/pbc.h
index fce60b814f..acda69c65e 100644
--- a/src/gromacs/pbcutil/pbc.h
+++ b/src/gromacs/pbcutil/pbc.h
@@ -260,7 +260,7 @@ void pbc_dx(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx);
  * \param[in]    x2  Coordinates for particle 2
  * \param[out]   dx  Distance vector
  * \return the ishift required to shift x1 at closest distance to x2;
- * i.e. if 0<=ishift<SHIFTS then x1 - x2 + shift_vec[ishift] = dx
+ * i.e. if 0<=ishift<c_numShiftVectors then x1 - x2 + shift_vec[ishift] = dx
  * (see calc_shifts below on how to obtain shift_vec)
  */
 int pbc_dx_aiuc(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx);
diff --git a/src/gromacs/pbcutil/pbc_aiuc_cuda.cuh b/src/gromacs/pbcutil/pbc_aiuc_cuda.cuh
index 96661613f6..4911850f4f 100644
--- a/src/gromacs/pbcutil/pbc_aiuc_cuda.cuh
+++ b/src/gromacs/pbcutil/pbc_aiuc_cuda.cuh
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 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.
@@ -54,8 +54,20 @@
 #define GMX_PBCUTIL_PBC_AIUC_CUDA_CUH
 
 #include "gromacs/gpu_utils/vectype_ops.cuh"
+#include "gromacs/pbcutil/ishift.h"
 #include "gromacs/pbcutil/pbc_aiuc.h"
 
+static inline __device__ int xyzToShiftIndex(int x, int y, int z)
+{
+    return (gmx::detail::c_nBoxX * (gmx::detail::c_nBoxY * ((z) + gmx::c_dBoxZ) + (y) + gmx::c_dBoxY)
+            + (x) + gmx::c_dBoxX);
+}
+
+static inline __device__ int int3ToShiftIndex(int3 iv)
+{
+    return (xyzToShiftIndex(iv.x, iv.y, iv.z));
+}
+
 /*! \brief Computes the vector between two points taking PBC into account.
  *
  * Computes the vector dr between points r2 and r1, taking into account the
@@ -102,13 +114,13 @@ static __forceinline__ __device__ int
 
     if (returnShift)
     {
-        ivec ishift;
+        int3 ishift;
 
-        ishift[XX] = -__float2int_rn(shx);
-        ishift[YY] = -__float2int_rn(shy);
-        ishift[ZZ] = -__float2int_rn(shz);
+        ishift.x = -__float2int_rn(shx);
+        ishift.y = -__float2int_rn(shy);
+        ishift.z = -__float2int_rn(shz);
 
-        return IVEC2IS(ishift);
+        return int3ToShiftIndex(ishift);
     }
     else
     {
diff --git a/src/gromacs/pbcutil/tests/pbc.cpp b/src/gromacs/pbcutil/tests/pbc.cpp
index 6ad87561fa..79c126c41b 100644
--- a/src/gromacs/pbcutil/tests/pbc.cpp
+++ b/src/gromacs/pbcutil/tests/pbc.cpp
@@ -63,7 +63,7 @@ TEST(PbcTest, CalcShiftsWorks)
     // shift vector values when the largest box shift in any dimension
     // is two.
     const matrix           box = { { 0.01, 1, -100 }, { 300, -0.03, 3 }, { -6, -600, 0.06 } };
-    std::vector<gmx::RVec> shiftVectors(SHIFTS);
+    std::vector<gmx::RVec> shiftVectors(c_numShiftVectors);
 
     calc_shifts(box, shiftVectors);