diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index fc867054c1..3cc1a6fbe3 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -38,6 +38,9 @@ Change Log
   (`#1869 <https://github.com/glotzerlab/hoomd-blue/pull/1869>`__).
 * ``P`` property of ``hoomd.hpmc.compute.SDF``
   (`#1869 <https://github.com/glotzerlab/hoomd-blue/pull/1869>`__).
+* MPCD supports non-cubic collision cells, including both orthorhombic and triclinic cell shapes.
+  The cell size is controlled by the number of cells along each lattice vector that defines the
+  simulation box (`#1950 <https://github.com/glotzerlab/hoomd-blue/pull/1950>`__).
 
 *Changed*
 
diff --git a/hoomd/mpcd/CellList.cc b/hoomd/mpcd/CellList.cc
index 25efd67c28..e4e7b7a702 100644
--- a/hoomd/mpcd/CellList.cc
+++ b/hoomd/mpcd/CellList.cc
@@ -18,22 +18,53 @@
 namespace hoomd
     {
 mpcd::CellList::CellList(std::shared_ptr<SystemDefinition> sysdef, Scalar cell_size, bool shift)
-    : Compute(sysdef), m_mpcd_pdata(m_sysdef->getMPCDParticleData()), m_cell_size(cell_size),
-      m_cell_np_max(4), m_cell_np(m_exec_conf), m_cell_list(m_exec_conf),
-      m_embed_cell_ids(m_exec_conf), m_conditions(m_exec_conf), m_needs_compute_dim(true),
-      m_particles_sorted(false), m_virtual_change(false)
+    : Compute(sysdef), m_mpcd_pdata(m_sysdef->getMPCDParticleData()), m_cell_np_max(4),
+      m_cell_np(m_exec_conf), m_cell_list(m_exec_conf), m_embed_cell_ids(m_exec_conf),
+      m_conditions(m_exec_conf), m_needs_compute_dim(true), m_particles_sorted(false),
+      m_virtual_change(false)
     {
     assert(m_mpcd_pdata);
     m_exec_conf->msg->notice(5) << "Constructing MPCD CellList" << std::endl;
 
-    // by default, grid shifting is initialized to zeroes
+    setCellSize(cell_size);
+    m_origin_idx = make_int3(0, 0, 0);
     m_cell_dim = make_uint3(0, 0, 0);
-    m_global_cell_dim = make_uint3(0, 0, 0);
 
     m_enable_grid_shift = shift;
     m_grid_shift = make_scalar3(0.0, 0.0, 0.0);
-    m_max_grid_shift = 0.5 * m_cell_size;
+
+    resetConditions();
+
+#ifdef ENABLE_MPI
+    m_decomposition = m_pdata->getDomainDecomposition();
+    m_num_extra = 0;
+    m_cover_box = m_pdata->getBox();
+#endif // ENABLE_MPI
+
+    m_mpcd_pdata->getSortSignal().connect<mpcd::CellList, &mpcd::CellList::sort>(this);
+    m_mpcd_pdata->getNumVirtualSignal().connect<mpcd::CellList, &mpcd::CellList::slotNumVirtual>(
+        this);
+    m_pdata->getParticleSortSignal().connect<mpcd::CellList, &mpcd::CellList::slotSorted>(this);
+    m_pdata->getBoxChangeSignal().connect<mpcd::CellList, &mpcd::CellList::slotBoxChanged>(this);
+    }
+
+mpcd::CellList::CellList(std::shared_ptr<SystemDefinition> sysdef,
+                         const uint3& global_cell_dim,
+                         bool shift)
+    : Compute(sysdef), m_mpcd_pdata(m_sysdef->getMPCDParticleData()), m_cell_np_max(4),
+      m_cell_np(m_exec_conf), m_cell_list(m_exec_conf), m_embed_cell_ids(m_exec_conf),
+      m_conditions(m_exec_conf), m_needs_compute_dim(true), m_particles_sorted(false),
+      m_virtual_change(false)
+    {
+    assert(m_mpcd_pdata);
+    m_exec_conf->msg->notice(5) << "Constructing MPCD CellList" << std::endl;
+
+    setGlobalDim(global_cell_dim);
     m_origin_idx = make_int3(0, 0, 0);
+    m_cell_dim = make_uint3(0, 0, 0);
+
+    m_enable_grid_shift = shift;
+    m_grid_shift = make_scalar3(0.0, 0.0, 0.0);
 
     resetConditions();
 
@@ -140,55 +171,11 @@ void mpcd::CellList::reallocate()
     m_cell_list.resize(m_cell_list_indexer.getNumElements());
     }
 
-void mpcd::CellList::updateGlobalBox()
-    {
-    // triclinic boxes are not allowed
-    const BoxDim& global_box = m_pdata->getGlobalBox();
-    if (global_box.getTiltFactorXY() != Scalar(0.0) || global_box.getTiltFactorXZ() != Scalar(0.0)
-        || global_box.getTiltFactorYZ() != Scalar(0.0))
-        {
-        m_exec_conf->msg->error() << "mpcd: box must be orthorhombic" << std::endl;
-        throw std::runtime_error("Box must be orthorhombic");
-        }
-
-    // box must be evenly divisible by cell size
-    const Scalar3 L = global_box.getL();
-    m_global_cell_dim = make_uint3((unsigned int)round(L.x / m_cell_size),
-                                   (unsigned int)round(L.y / m_cell_size),
-                                   (unsigned int)round(L.z / m_cell_size));
-    if (m_sysdef->getNDimensions() == 2)
-        {
-        if (m_global_cell_dim.z > 1)
-            {
-            m_exec_conf->msg->error()
-                << "mpcd: In 2d simulations, box width must be smaller than cell size" << std::endl;
-            throw std::runtime_error("Lz bigger than cell size in 2D!");
-            }
-
-        // force to be only one cell along z in 2d
-        m_global_cell_dim.z = 1;
-        }
-
-    const double eps = 1e-5;
-    if (fabs((double)L.x - m_global_cell_dim.x * (double)m_cell_size) > eps * m_cell_size
-        || fabs((double)L.y - m_global_cell_dim.y * (double)m_cell_size) > eps * m_cell_size
-        || (m_sysdef->getNDimensions() == 3
-            && fabs((double)L.z - m_global_cell_dim.z * (double)m_cell_size) > eps * m_cell_size))
-        {
-        m_exec_conf->msg->error() << "mpcd: Box size must be even multiple of cell size"
-                                  << std::endl;
-        throw std::runtime_error("MPCD cell size must evenly divide box");
-        }
-    }
-
 void mpcd::CellList::computeDimensions()
     {
     if (!m_needs_compute_dim)
         return;
 
-    // first update / validate the global box
-    updateGlobalBox();
-
 #ifdef ENABLE_MPI
     uchar3 communicating = make_uchar3(0, 0, 0);
     if (m_decomposition)
@@ -202,34 +189,39 @@ void mpcd::CellList::computeDimensions()
     // Only do complicated sizing if some direction is being communicated
     if (communicating.x || communicating.y || communicating.z)
         {
-        // Global simulation box for absolute position referencing
-        const BoxDim& global_box = m_pdata->getGlobalBox();
-        const Scalar3 global_lo = global_box.getLo();
-        const Scalar3 global_L = global_box.getL();
-
-        // if global box is valid (no triclinic skew), then local box also has no skew,
-        // so assume orthorhombic in all subsequent calculations
-        const BoxDim& box = m_pdata->getBox();
+        // fractional bounds of domain
+        const auto grid_pos = m_decomposition->getGridPos();
+        Scalar3 fractional_lo = make_scalar3(m_decomposition->getCumulativeFraction(0, grid_pos.x),
+                                             m_decomposition->getCumulativeFraction(1, grid_pos.y),
+                                             m_decomposition->getCumulativeFraction(2, grid_pos.z));
+        const Scalar3 fractional_hi
+            = make_scalar3(m_decomposition->getCumulativeFraction(0, grid_pos.x + 1),
+                           m_decomposition->getCumulativeFraction(1, grid_pos.y + 1),
+                           m_decomposition->getCumulativeFraction(2, grid_pos.z + 1));
 
         // setup lo bin
-        const Scalar3 delta_lo = box.getLo() - global_lo;
-        int3 my_lo_bin = make_int3((int)std::floor((delta_lo.x - m_max_grid_shift) / m_cell_size),
-                                   (int)std::floor((delta_lo.y - m_max_grid_shift) / m_cell_size),
-                                   (int)std::floor((delta_lo.z - m_max_grid_shift) / m_cell_size));
+        const Scalar3 fractional_lo_shifted_down = fractional_lo - m_max_grid_shift;
+        int3 my_lo_bin
+            = make_int3((int)std::floor(fractional_lo_shifted_down.x * m_global_cell_dim.x),
+                        (int)std::floor(fractional_lo_shifted_down.y * m_global_cell_dim.y),
+                        (int)std::floor(fractional_lo_shifted_down.z * m_global_cell_dim.z));
+        const Scalar3 fractional_lo_shifted_up = fractional_lo + m_max_grid_shift;
         int3 lo_neigh_bin
-            = make_int3((int)std::ceil((delta_lo.x + m_max_grid_shift) / m_cell_size),
-                        (int)std::ceil((delta_lo.y + m_max_grid_shift) / m_cell_size),
-                        (int)std::ceil((delta_lo.z + m_max_grid_shift) / m_cell_size));
+            = make_int3((int)std::ceil(fractional_lo_shifted_up.x * m_global_cell_dim.x),
+                        (int)std::ceil(fractional_lo_shifted_up.y * m_global_cell_dim.y),
+                        (int)std::ceil(fractional_lo_shifted_up.z * m_global_cell_dim.z));
 
         // setup hi bin
-        const Scalar3 delta_hi = box.getHi() - global_lo;
-        int3 my_hi_bin = make_int3((int)std::ceil((delta_hi.x + m_max_grid_shift) / m_cell_size),
-                                   (int)std::ceil((delta_hi.y + m_max_grid_shift) / m_cell_size),
-                                   (int)std::ceil((delta_hi.z + m_max_grid_shift) / m_cell_size));
+        const Scalar3 fractional_hi_shifted_up = fractional_hi + m_max_grid_shift;
+        int3 my_hi_bin
+            = make_int3((int)std::ceil(fractional_hi_shifted_up.x * m_global_cell_dim.x),
+                        (int)std::ceil(fractional_hi_shifted_up.y * m_global_cell_dim.y),
+                        (int)std::ceil(fractional_hi_shifted_up.z * m_global_cell_dim.z));
+        const Scalar3 fractional_hi_shifted_down = fractional_hi - m_max_grid_shift;
         int3 hi_neigh_bin
-            = make_int3((int)std::floor((delta_hi.x - m_max_grid_shift) / m_cell_size),
-                        (int)std::floor((delta_hi.y - m_max_grid_shift) / m_cell_size),
-                        (int)std::floor((delta_hi.z - m_max_grid_shift) / m_cell_size));
+            = make_int3((int)std::floor(fractional_hi_shifted_down.x * m_global_cell_dim.x),
+                        (int)std::floor(fractional_hi_shifted_down.y * m_global_cell_dim.y),
+                        (int)std::floor(fractional_hi_shifted_down.z * m_global_cell_dim.z));
 
         // initially size the grid assuming one rank in each direction, and then resize based on
         // communication
@@ -238,9 +230,10 @@ void mpcd::CellList::computeDimensions()
         std::fill(m_num_comm.begin(), m_num_comm.end(), 0);
 
         // Compute size of the box with diffusion layer
-        Scalar3 cover_lo = box.getLo();
-        Scalar3 cover_hi = box.getHi();
-        uchar3 cover_periodic = box.getPeriodic();
+        const BoxDim& global_box = m_pdata->getGlobalBox();
+        uchar3 cover_periodic = global_box.getPeriodic();
+        Scalar3 fractional_cover_lo = fractional_lo;
+        Scalar3 fractional_cover_hi = fractional_hi;
 
         if (communicating.x)
             {
@@ -255,9 +248,9 @@ void mpcd::CellList::computeDimensions()
                 = lo_neigh_bin.x - my_lo_bin.x + m_num_extra;
 
             // "safe" size of the diffusion layer
-            cover_lo.x = m_origin_idx.x * m_cell_size + m_max_grid_shift - 0.5 * global_L.x;
-            cover_hi.x = (m_origin_idx.x + m_cell_dim.x) * m_cell_size - m_max_grid_shift
-                         - 0.5 * global_L.x;
+            fractional_cover_lo.x = m_global_cell_dim_inv.x * m_origin_idx.x + m_max_grid_shift.x;
+            fractional_cover_hi.x
+                = m_global_cell_dim_inv.x * (m_origin_idx.x + m_cell_dim.x) - m_max_grid_shift.x;
             cover_periodic.x = 0;
             }
 
@@ -271,9 +264,9 @@ void mpcd::CellList::computeDimensions()
             m_num_comm[static_cast<unsigned int>(mpcd::detail::face::south)]
                 = lo_neigh_bin.y - my_lo_bin.y + m_num_extra;
 
-            cover_lo.y = m_origin_idx.y * m_cell_size + m_max_grid_shift - 0.5 * global_L.y;
-            cover_hi.y = (m_origin_idx.y + m_cell_dim.y) * m_cell_size - m_max_grid_shift
-                         - 0.5 * global_L.y;
+            fractional_cover_lo.y = m_global_cell_dim_inv.y * m_origin_idx.y + m_max_grid_shift.y;
+            fractional_cover_hi.y
+                = m_global_cell_dim_inv.y * (m_origin_idx.y + m_cell_dim.y) - m_max_grid_shift.y;
             cover_periodic.y = 0;
             }
 
@@ -287,16 +280,22 @@ void mpcd::CellList::computeDimensions()
             m_num_comm[static_cast<unsigned int>(mpcd::detail::face::down)]
                 = lo_neigh_bin.z - my_lo_bin.z + m_num_extra;
 
-            cover_lo.z = m_origin_idx.z * m_cell_size + m_max_grid_shift - 0.5 * global_L.z;
-            cover_hi.z = (m_origin_idx.z + m_cell_dim.z) * m_cell_size - m_max_grid_shift
-                         - 0.5 * global_L.z;
+            fractional_cover_lo.z = m_global_cell_dim_inv.z * m_origin_idx.z + m_max_grid_shift.z;
+            fractional_cover_hi.z
+                = m_global_cell_dim_inv.z * (m_origin_idx.z + m_cell_dim.z) - m_max_grid_shift.z;
             cover_periodic.z = 0;
             }
 
         // set the box covered by this cell list
-        m_cover_box = BoxDim(cover_lo, cover_hi, cover_periodic);
-
-        checkDomainBoundaries();
+        // this requires us to reduce the box back to equivalent cube lengths
+        const Scalar3 global_lo = global_box.getLo();
+        const Scalar3 global_L = global_box.getL();
+        m_cover_box = BoxDim(global_lo + fractional_cover_lo * global_L,
+                             global_lo + fractional_cover_hi * global_L,
+                             cover_periodic);
+        m_cover_box.setTiltFactors(global_box.getTiltFactorXY(),
+                                   global_box.getTiltFactorXZ(),
+                                   global_box.getTiltFactorYZ());
         }
     else
 #endif // ENABLE_MPI
@@ -319,184 +318,6 @@ void mpcd::CellList::computeDimensions()
     }
 
 #ifdef ENABLE_MPI
-void mpcd::CellList::checkDomainBoundaries()
-    {
-    if (!m_decomposition)
-        return;
-
-    MPI_Comm mpi_comm = m_exec_conf->getMPICommunicator();
-
-    for (unsigned int dir = 0; dir < m_num_comm.size(); ++dir)
-        {
-        mpcd::detail::face d = static_cast<mpcd::detail::face>(dir);
-        if (!isCommunicating(d))
-            continue;
-
-        // receive in the opposite direction from which we send
-        unsigned int send_neighbor = m_decomposition->getNeighborRank(dir);
-        unsigned int recv_neighbor;
-        if (dir % 2 == 0)
-            recv_neighbor = m_decomposition->getNeighborRank(dir + 1);
-        else
-            recv_neighbor = m_decomposition->getNeighborRank(dir - 1);
-
-        // first make sure each dimension is sending and receiving the same size data
-        MPI_Request reqs[2];
-        MPI_Status status[2];
-
-        // check that the number received is the same as that being sent from neighbor
-        unsigned int n_send = m_num_comm[dir];
-        unsigned int n_expect_recv;
-        if (dir % 2 == 0)
-            n_expect_recv = m_num_comm[dir + 1];
-        else
-            n_expect_recv = m_num_comm[dir - 1];
-
-        unsigned int n_recv;
-        MPI_Isend(&n_send, 1, MPI_UNSIGNED, send_neighbor, 0, mpi_comm, &reqs[0]);
-        MPI_Irecv(&n_recv, 1, MPI_UNSIGNED, recv_neighbor, 0, mpi_comm, &reqs[1]);
-        MPI_Waitall(2, reqs, status);
-
-        // check if any rank errored out
-        unsigned int recv_error = 0;
-        if (n_expect_recv != n_recv)
-            {
-            recv_error = 1;
-            }
-        unsigned int any_error = 0;
-        MPI_Allreduce(&recv_error, &any_error, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm);
-        if (any_error)
-            {
-            if (recv_error)
-                {
-                m_exec_conf->msg->error() << "mpcd: expected to communicate " << n_expect_recv
-                                          << " cells, but only receiving " << n_recv << std::endl;
-                }
-            throw std::runtime_error("Error setting up MPCD cell list");
-            }
-
-        // check that the same cell ids are communicated
-        std::vector<int> send_cells(n_send), recv_cells(n_recv);
-        for (unsigned int i = 0; i < n_send; ++i)
-            {
-            if (d == mpcd::detail::face::east)
-                {
-                send_cells[i] = m_origin_idx.x + m_cell_dim.x - m_num_extra - n_send + i;
-                }
-            else if (d == mpcd::detail::face::west)
-                {
-                send_cells[i] = m_origin_idx.x + i;
-                }
-            else if (d == mpcd::detail::face::north)
-                {
-                send_cells[i] = m_origin_idx.y + m_cell_dim.y - m_num_extra - n_send + i;
-                }
-            else if (d == mpcd::detail::face::south)
-                {
-                send_cells[i] = m_origin_idx.y + i;
-                }
-            else if (d == mpcd::detail::face::up)
-                {
-                send_cells[i] = m_origin_idx.z + m_cell_dim.z - m_num_extra - n_send + i;
-                }
-            else if (d == mpcd::detail::face::down)
-                {
-                send_cells[i] = m_origin_idx.z + i;
-                }
-            }
-
-        MPI_Isend(&send_cells[0], n_send, MPI_INT, send_neighbor, 1, mpi_comm, &reqs[0]);
-        MPI_Irecv(&recv_cells[0], n_recv, MPI_INT, recv_neighbor, 1, mpi_comm, &reqs[1]);
-        MPI_Waitall(2, reqs, status);
-
-        unsigned int overlap_error = 0;
-        std::array<int, 2> err_pair {0, 0};
-        for (unsigned int i = 0; i < n_recv && !overlap_error; ++i)
-            {
-            // wrap the received cell back into the global box
-            // only two of the entries will be valid, the others are dummies
-            int3 recv_cell = make_int3(0, 0, 0);
-            if (d == mpcd::detail::face::east || d == mpcd::detail::face::west)
-                {
-                recv_cell.x = recv_cells[i];
-                }
-            else if (d == mpcd::detail::face::north || d == mpcd::detail::face::south)
-                {
-                recv_cell.y = recv_cells[i];
-                }
-            else if (d == mpcd::detail::face::up || d == mpcd::detail::face::down)
-                {
-                recv_cell.z = recv_cells[i];
-                }
-            recv_cell = wrapGlobalCell(recv_cell);
-
-            // compute the expected cell to receive, also wrapped
-            int3 expect_recv_cell = make_int3(0, 0, 0);
-            if (d == mpcd::detail::face::east)
-                {
-                expect_recv_cell.x = m_origin_idx.x + i;
-                }
-            else if (d == mpcd::detail::face::west)
-                {
-                expect_recv_cell.x = m_origin_idx.x + m_cell_dim.x - m_num_extra - n_recv + i;
-                }
-            else if (d == mpcd::detail::face::north)
-                {
-                expect_recv_cell.y = m_origin_idx.y + i;
-                }
-            else if (d == mpcd::detail::face::south)
-                {
-                expect_recv_cell.y = m_origin_idx.y + m_cell_dim.y - m_num_extra - n_recv + i;
-                }
-            else if (d == mpcd::detail::face::up)
-                {
-                expect_recv_cell.z = m_origin_idx.z + i;
-                }
-            else if (d == mpcd::detail::face::down)
-                {
-                expect_recv_cell.z = m_origin_idx.z + m_cell_dim.z - m_num_extra - n_recv + i;
-                }
-            expect_recv_cell = wrapGlobalCell(expect_recv_cell);
-
-            if (recv_cell.x != expect_recv_cell.x || recv_cell.y != expect_recv_cell.y
-                || recv_cell.z != expect_recv_cell.z)
-                {
-                overlap_error = i;
-                if (d == mpcd::detail::face::east || d == mpcd::detail::face::west)
-                    {
-                    err_pair[0] = recv_cell.x;
-                    err_pair[1] = expect_recv_cell.x;
-                    }
-                else if (d == mpcd::detail::face::north || d == mpcd::detail::face::south)
-                    {
-                    err_pair[0] = recv_cell.y;
-                    err_pair[1] = expect_recv_cell.y;
-                    }
-                else if (d == mpcd::detail::face::up || d == mpcd::detail::face::down)
-                    {
-                    err_pair[0] = recv_cell.z;
-                    err_pair[1] = expect_recv_cell.z;
-                    }
-                }
-            }
-
-        // check if anyone reported an error, then race to see who gets to write it out
-        any_error = 0;
-        MPI_Allreduce(&overlap_error, &any_error, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm);
-        if (any_error)
-            {
-            if (overlap_error)
-                {
-                m_exec_conf->msg->error()
-                    << "mpcd: communication grid does not overlap. " << "Expected to receive cell "
-                    << err_pair[1] << " from rank " << recv_neighbor << ", but got cell "
-                    << err_pair[0] << "." << std::endl;
-                }
-            throw std::runtime_error("Error setting up MPCD cell list");
-            }
-        }
-    }
-
 /*!
  * \param dir Direction of communication
  * \returns True if communication is occurring along \a dir
@@ -579,7 +400,7 @@ void mpcd::CellList::buildCellList()
         n_global_cells.z += 2 * m_num_extra;
 #endif // ENABLE_MPI
 
-    const Scalar3 global_lo = m_pdata->getGlobalBox().getLo();
+    const BoxDim& global_box = m_pdata->getGlobalBox();
 
     for (unsigned int cur_p = 0; cur_p < N_tot; ++cur_p)
         {
@@ -600,11 +421,11 @@ void mpcd::CellList::buildCellList()
             continue;
             }
 
-        // bin particle assuming orthorhombic box (already validated)
-        const Scalar3 delta = (pos_i - m_grid_shift) - global_lo;
-        int3 global_bin = make_int3((int)std::floor(delta.x / m_cell_size),
-                                    (int)std::floor(delta.y / m_cell_size),
-                                    (int)std::floor(delta.z / m_cell_size));
+        // bin particle
+        const Scalar3 fractional_pos_i = global_box.makeFraction(pos_i) - m_grid_shift;
+        int3 global_bin = make_int3((int)std::floor(fractional_pos_i.x * m_global_cell_dim.x),
+                                    (int)std::floor(fractional_pos_i.y * m_global_cell_dim.y),
+                                    (int)std::floor(fractional_pos_i.z * m_global_cell_dim.z));
 
         // wrap cell back through the boundaries (grid shifting may send +/- 1 outside of range)
         // this is done using periodic from the "local" box, since this will be periodic
@@ -635,6 +456,28 @@ void mpcd::CellList::buildCellList()
         int3 bin = make_int3(global_bin.x - m_origin_idx.x,
                              global_bin.y - m_origin_idx.y,
                              global_bin.z - m_origin_idx.z);
+        // these checks guard against round-off errors with domain decomposition
+        if (!periodic.x)
+            {
+            if (bin.x == -1)
+                bin.x = 0;
+            else if (bin.x == (int)m_cell_dim.x)
+                bin.x = m_cell_dim.x - 1;
+            }
+        if (!periodic.y)
+            {
+            if (bin.y == -1)
+                bin.y = 0;
+            else if (bin.y == (int)m_cell_dim.y)
+                bin.y = m_cell_dim.y - 1;
+            }
+        if (!periodic.z)
+            {
+            if (bin.z == -1)
+                bin.z = 0;
+            else if (bin.z == (int)m_cell_dim.z)
+                bin.z = m_cell_dim.z - 1;
+            }
 
         // validate and make sure no particles blew out of the box
         if ((bin.x < 0 || bin.x >= (int)m_cell_dim.x) || (bin.y < 0 || bin.y >= (int)m_cell_dim.y)
@@ -688,15 +531,16 @@ void mpcd::CellList::sort(uint64_t timestep,
         return;
 
     // if mapping is not valid, signal that we need to force a recompute next time
-    // that the cell list is needed. We don't call forceCompute() directly because this always
-    // runs compute(), and we just want to defer to the next compute() call.
+    // that the cell list is needed. We don't call forceCompute() directly because this
+    // always runs compute(), and we just want to defer to the next compute() call.
     if (rorder.isNull())
         {
         m_force_compute = true;
         return;
         }
 
-    // iterate through particles in cell list, and update their indexes using reverse mapping
+    // iterate through particles in cell list, and update their indexes using reverse
+    // mapping
     ArrayHandle<unsigned int> h_rorder(rorder, access_location::host, access_mode::read);
     ArrayHandle<unsigned int> h_cell_np(m_cell_np, access_location::host, access_mode::read);
     ArrayHandle<unsigned int> h_cell_list(m_cell_list,
@@ -730,9 +574,7 @@ bool mpcd::CellList::needsEmbedMigrate(uint64_t timestep)
     // ensure that the cell list has been sized first
     computeDimensions();
 
-    // coverage box dimensions, assuming orthorhombic
-    const Scalar3 lo = m_cover_box.getLo();
-    const Scalar3 hi = m_cover_box.getHi();
+    // coverage box dimensions
     const uchar3 periodic = m_cover_box.getPeriodic();
     const unsigned int ndim = m_sysdef->getNDimensions();
 
@@ -750,10 +592,11 @@ bool mpcd::CellList::needsEmbedMigrate(uint64_t timestep)
         const unsigned int idx = h_group.data[i];
         const Scalar4 postype = h_pos.data[idx];
         const Scalar3 pos = make_scalar3(postype.x, postype.y, postype.z);
-
-        if ((!periodic.x && (pos.x >= hi.x || pos.x < lo.x))
-            || (!periodic.y && (pos.y >= hi.y || pos.y < lo.y))
-            || (!periodic.z && ndim == 3 && (pos.z >= hi.z || pos.z < lo.z)))
+        const Scalar3 fractional_pos = m_cover_box.makeFraction(pos);
+        if ((!periodic.x && (fractional_pos.x >= Scalar(1.0) || fractional_pos.x < Scalar(0.0)))
+            || (!periodic.y && (fractional_pos.y >= Scalar(1.0) || fractional_pos.y < Scalar(0.0)))
+            || (!periodic.z && ndim == 3
+                && (fractional_pos.z >= Scalar(1.0) || fractional_pos.z < Scalar(0.0))))
             {
             migrate = 1;
             }
@@ -802,6 +645,7 @@ bool mpcd::CellList::checkConditions()
         {
         unsigned int n = conditions.z - 1;
         Scalar4 pos_empty_i;
+        std::string msg;
         if (n < m_mpcd_pdata->getN() + m_mpcd_pdata->getNVirtual())
             {
             ArrayHandle<Scalar4> h_pos(m_mpcd_pdata->getPositions(),
@@ -809,11 +653,9 @@ bool mpcd::CellList::checkConditions()
                                        access_mode::read);
             pos_empty_i = h_pos.data[n];
             if (n < m_mpcd_pdata->getN())
-                m_exec_conf->msg->errorAllRanks()
-                    << "MPCD particle is no longer in the simulation box" << std::endl;
+                msg = "MPCD particle is no longer in the simulation box";
             else
-                m_exec_conf->msg->errorAllRanks()
-                    << "MPCD virtual particle is no longer in the simulation box" << std::endl;
+                msg = "MPCD virtual particle is no longer in the simulation box";
             }
         else
             {
@@ -827,27 +669,17 @@ bool mpcd::CellList::checkConditions()
                 = h_pos_embed
                       .data[h_embed_member_idx
                                 .data[n - (m_mpcd_pdata->getN() + m_mpcd_pdata->getNVirtual())]];
-            m_exec_conf->msg->errorAllRanks()
-                << "Embedded particle is no longer in the simulation box" << std::endl;
+            msg = "Embedded particle is no longer in the simulation box";
             }
 
         Scalar3 pos = make_scalar3(pos_empty_i.x, pos_empty_i.y, pos_empty_i.z);
         m_exec_conf->msg->errorAllRanks()
+            << msg << std::endl
             << "Cartesian coordinates: " << std::endl
             << "x: " << pos.x << " y: " << pos.y << " z: " << pos.z << std::endl
             << "Grid shift: " << std::endl
             << "x: " << m_grid_shift.x << " y: " << m_grid_shift.y << " z: " << m_grid_shift.z
             << std::endl;
-
-        const BoxDim cover_box = getCoverageBox();
-        Scalar3 lo = cover_box.getLo();
-        Scalar3 hi = cover_box.getHi();
-        uchar3 periodic = cover_box.getPeriodic();
-        m_exec_conf->msg->errorAllRanks()
-            << "Covered box lo: (" << lo.x << ", " << lo.y << ", " << lo.z << ")" << std::endl
-            << "            hi: (" << hi.x << ", " << hi.y << ", " << hi.z << ")" << std::endl
-            << "      periodic: (" << ((periodic.x) ? "1" : "0") << " "
-            << ((periodic.y) ? "1" : "0") << " " << ((periodic.z) ? "1" : "0") << ")" << std::endl;
         throw std::runtime_error("Error computing cell list");
         }
 
@@ -873,6 +705,8 @@ void mpcd::CellList::drawGridShift(uint64_t timestep)
     {
     if (m_enable_grid_shift)
         {
+        computeDimensions();
+
         uint16_t seed = m_sysdef->getSeed();
 
         // PRNG using seed and timestep as seeds
@@ -880,29 +714,35 @@ void mpcd::CellList::drawGridShift(uint64_t timestep)
                                    hoomd::Counter());
 
         // draw shift variables from uniform distribution
-        hoomd::UniformDistribution<Scalar> uniform(-m_max_grid_shift, m_max_grid_shift);
-        Scalar3 shift;
-        shift.x = uniform(rng);
-        shift.y = uniform(rng);
-        shift.z = (m_sysdef->getNDimensions() == 3) ? uniform(rng) : Scalar(0.0);
+        hoomd::UniformDistribution<Scalar> uniform(Scalar(-0.5), Scalar(0.5));
+        Scalar3 shift = make_scalar3(
+            uniform(rng) * m_global_cell_dim_inv.x,
+            uniform(rng) * m_global_cell_dim_inv.y,
+            (m_sysdef->getNDimensions() == 3) ? uniform(rng) * m_global_cell_dim_inv.z : 0);
         setGridShift(shift);
         }
     }
 
-void mpcd::CellList::getCellStatistics() const
+void mpcd::CellList::setGlobalDim(const uint3& global_cell_dim)
     {
-    unsigned int min_np(0xffffffff), max_np(0);
-    ArrayHandle<unsigned int> h_cell_np(m_cell_np, access_location::host, access_mode::read);
-    for (unsigned int cur_cell = 0; cur_cell < m_cell_indexer.getNumElements(); ++cur_cell)
+    if (global_cell_dim.x == 0 || global_cell_dim.y == 0)
         {
-        const unsigned int np = h_cell_np.data[cur_cell];
-        if (np < min_np)
-            min_np = np;
-        if (np > max_np)
-            max_np = np;
+        throw std::runtime_error("Global cell dimensions must be at least 1");
         }
-    m_exec_conf->msg->notice(2) << "MPCD cell list stats:" << std::endl;
-    m_exec_conf->msg->notice(2) << "Min: " << min_np << " Max: " << max_np << std::endl;
+
+    m_global_cell_dim = global_cell_dim;
+    if (m_sysdef->getNDimensions() == 2)
+        {
+        m_global_cell_dim.z = 1;
+        }
+
+    m_global_cell_dim_inv = make_scalar3(Scalar(1.0) / global_cell_dim.x,
+                                         Scalar(1.0) / global_cell_dim.y,
+                                         Scalar(1.0) / global_cell_dim.z);
+
+    m_max_grid_shift = Scalar(0.5) * m_global_cell_dim_inv;
+
+    m_needs_compute_dim = true;
     }
 
 /*!
@@ -912,8 +752,10 @@ void mpcd::CellList::getCellStatistics() const
  * \warning The returned cell index may lie outside the local grid. It is the
  *          caller's responsibility to check that the index is valid.
  */
-const int3 mpcd::CellList::getLocalCell(const int3& global) const
+const int3 mpcd::CellList::getLocalCell(const int3& global)
     {
+    computeDimensions();
+
     int3 local = make_int3(global.x - m_origin_idx.x,
                            global.y - m_origin_idx.y,
                            global.z - m_origin_idx.z);
@@ -928,21 +770,65 @@ const int3 mpcd::CellList::getLocalCell(const int3& global) const
  * Local cell coordinates are wrapped around the global box so that a valid global
  * index is computed.
  */
-const int3 mpcd::CellList::getGlobalCell(const int3& local) const
+const int3 mpcd::CellList::getGlobalCell(const int3& local)
     {
+    computeDimensions();
+
     int3 global
         = make_int3(local.x + m_origin_idx.x, local.y + m_origin_idx.y, local.z + m_origin_idx.z);
     return wrapGlobalCell(global);
     }
 
+Scalar3 mpcd::CellList::getCellSize()
+    {
+    computeDimensions();
+
+    const BoxDim& global_box = m_pdata->getGlobalBox();
+    const Scalar3 L = global_box.getL();
+    return make_scalar3(L.x * m_global_cell_dim_inv.x,
+                        L.y * m_global_cell_dim_inv.y,
+                        L.z * m_global_cell_dim_inv.z);
+    }
+
+/*!
+ * \param cell_size Grid spacing
+ * \note Calling forces a resize of the cell list on the next update
+ */
+void mpcd::CellList::setCellSize(Scalar cell_size)
+    {
+    const BoxDim& global_box = m_pdata->getGlobalBox();
+    const Scalar3 L = global_box.getL();
+    uint3 global_cell_dim = make_uint3((unsigned int)round(L.x / cell_size),
+                                       (unsigned int)round(L.y / cell_size),
+                                       (unsigned int)round(L.z / cell_size));
+    if (m_sysdef->getNDimensions() == 2)
+        {
+        global_cell_dim.z = 1;
+        }
+
+    // check that box is a multiple of cell size
+    const double eps = 1e-5;
+    if (fabs((double)L.x - global_cell_dim.x * (double)cell_size) > eps * cell_size
+        || fabs((double)L.y - global_cell_dim.y * (double)cell_size) > eps * cell_size
+        || (m_sysdef->getNDimensions() == 3
+            && fabs((double)L.z - global_cell_dim.z * (double)cell_size) > eps * cell_size))
+        {
+        throw std::runtime_error("MPCD cell size must evenly divide box");
+        }
+
+    setGlobalDim(global_cell_dim);
+    }
+
 /*!
  * \param cell Cell coordinates to wrap back into the global box
  *
  * \warning Only up to one global box size is wrapped. This method is intended
  *          to be used for wrapping cells off by only one or two from the global boundary.
  */
-const int3 mpcd::CellList::wrapGlobalCell(const int3& cell) const
+const int3 mpcd::CellList::wrapGlobalCell(const int3& cell)
     {
+    computeDimensions();
+
     int3 wrap = cell;
 
     if (wrap.x >= (int)m_global_cell_dim.x)
@@ -971,7 +857,19 @@ void export_CellList(pybind11::module& m)
     {
     pybind11::class_<mpcd::CellList, Compute, std::shared_ptr<mpcd::CellList>>(m, "CellList")
         .def(pybind11::init<std::shared_ptr<SystemDefinition>, Scalar, bool>())
-        .def_property("cell_size", &mpcd::CellList::getCellSize, &mpcd::CellList::setCellSize)
+        .def_property(
+            "num_cells",
+            [](const mpcd::CellList& cl)
+            {
+                const auto num_cells = cl.getGlobalDim();
+                return pybind11::make_tuple(num_cells.x, num_cells.y, num_cells.z);
+            },
+            [](mpcd::CellList& cl, const pybind11::tuple& num_cells)
+            {
+                cl.setGlobalDim(make_uint3(pybind11::cast<unsigned int>(num_cells[0]),
+                                           pybind11::cast<unsigned int>(num_cells[1]),
+                                           pybind11::cast<unsigned int>(num_cells[2])));
+            })
         .def_property("shift",
                       &mpcd::CellList::isGridShifting,
                       &mpcd::CellList::enableGridShifting);
diff --git a/hoomd/mpcd/CellList.h b/hoomd/mpcd/CellList.h
index e76ea2ca2b..2ffa165f9e 100644
--- a/hoomd/mpcd/CellList.h
+++ b/hoomd/mpcd/CellList.h
@@ -33,9 +33,12 @@ namespace mpcd
 class PYBIND11_EXPORT CellList : public Compute
     {
     public:
-    //! Constructor
+    //! Constructor by size (deprecated)
     CellList(std::shared_ptr<SystemDefinition> sysdef, Scalar cell_size, bool shift);
 
+    //! Constructor by dimension
+    CellList(std::shared_ptr<SystemDefinition> sysdef, const uint3& global_cell_dim, bool shift);
+
     //! Destructor
     virtual ~CellList();
 
@@ -93,19 +96,21 @@ class PYBIND11_EXPORT CellList : public Compute
         return m_global_cell_dim;
         }
 
+    void setGlobalDim(const uint3& global_cell_dim);
+
     const int3& getOriginIndex() const
         {
         return m_origin_idx;
         }
 
     //! Obtain the local cell index corresponding to a global cell
-    const int3 getLocalCell(const int3& global) const;
+    const int3 getLocalCell(const int3& global);
 
     //! Obtain the global cell corresponding to local cell
-    const int3 getGlobalCell(const int3& local) const;
+    const int3 getGlobalCell(const int3& local);
 
     //! Wrap a cell into a global cell
-    const int3 wrapGlobalCell(const int3& cell) const;
+    const int3 wrapGlobalCell(const int3& cell);
 
     //! Get the maximum number of particles in a cell
     const unsigned int getNmax() const
@@ -113,23 +118,11 @@ class PYBIND11_EXPORT CellList : public Compute
         return m_cell_np_max;
         }
 
-    //! Set the MPCD cell size
-    /*!
-     * \param cell_size Grid spacing
-     * \note Calling forces a resize of the cell list on the next update
-     */
-    void setCellSize(Scalar cell_size)
-        {
-        m_cell_size = cell_size;
-        m_max_grid_shift = 0.5 * m_cell_size;
-        m_needs_compute_dim = true;
-        }
+    //! Get the MPCD cell size (deprecated)
+    Scalar3 getCellSize();
 
-    //! Get the MPCD cell size
-    Scalar getCellSize() const
-        {
-        return m_cell_size;
-        }
+    //! Set the MPCD cell size (deprecated)
+    void setCellSize(Scalar cell_size);
 
     //! Get the box that is covered by the cell list
     /*!
@@ -189,29 +182,27 @@ class PYBIND11_EXPORT CellList : public Compute
             }
         }
 
-    //! Get the maximum permitted grid shift
-    const Scalar getMaxGridShift() const
+    //! Get the maximum permitted grid shift (fractional coordinates)
+    const Scalar3 getMaxGridShift()
         {
+        computeDimensions();
         return m_max_grid_shift;
         }
 
-    // Get the grid shift vector
+    // Get the grid shift vector (fractional coordinates)
     const Scalar3& getGridShift() const
         {
         return m_grid_shift;
         }
 
-    //! Set the grid shift vector
+    //! Set the grid shift vector (fractional coordinates)
     void setGridShift(const Scalar3& shift)
         {
-        if (std::fabs(shift.x) > m_max_grid_shift || std::fabs(shift.y) > m_max_grid_shift
-            || std::fabs(shift.z) > m_max_grid_shift)
+        const Scalar3 max_grid_shift = getMaxGridShift();
+        if (std::fabs(shift.x) > max_grid_shift.x || std::fabs(shift.y) > max_grid_shift.y
+            || std::fabs(shift.z) > max_grid_shift.z)
             {
-            m_exec_conf->msg->error()
-                << "mpcd: Specified cell list grid shift (" << shift.x << ", " << shift.y << ", "
-                << shift.z << ")" << std::endl
-                << "exceeds maximum component magnitude " << m_max_grid_shift << std::endl;
-            throw std::runtime_error("Error setting MPCD grid shift");
+            throw std::runtime_error("MPCD grid shift out of range");
             }
 
         m_grid_shift = shift;
@@ -220,9 +211,6 @@ class PYBIND11_EXPORT CellList : public Compute
     //! Generates the random grid shift vector
     void drawGridShift(uint64_t timestep);
 
-    //! Calculate current cell occupancy statistics
-    virtual void getCellStatistics() const;
-
     //! Gets the group of particles that is coupled to the MPCD solvent through the collision step
     std::shared_ptr<ParticleGroup> getEmbeddedGroup() const
         {
@@ -261,14 +249,14 @@ class PYBIND11_EXPORT CellList : public Compute
 
     bool m_enable_grid_shift; //!< Flag to enable grid shifting
     Scalar3 m_grid_shift;     //!< Amount to shift particle positions when computing cell list
-    Scalar m_max_grid_shift;  //!< Maximum amount grid can be shifted in any direction
+    Scalar3 m_max_grid_shift; //!< Maximum amount grid can be shifted in any direction
 
     //! Allocates internal data arrays
     virtual void reallocate();
 
-    Scalar m_cell_size;            //!< MPCD cell width
     uint3 m_cell_dim;              //!< Number of cells in each direction
     uint3 m_global_cell_dim;       //!< Number of cells in each direction of global simulation box
+    Scalar3 m_global_cell_dim_inv; //!< Inverse of number of cells in each direction of global box
     Index3D m_cell_indexer;        //!< Indexer from 3D into cell list 1D
     Index3D m_global_cell_indexer; //!< Indexer from 3D into 1D for global cell indexes
     Index2D m_cell_list_indexer;   //!< Indexer into cell list members
@@ -337,9 +325,6 @@ class PYBIND11_EXPORT CellList : public Compute
 
 #ifdef ENABLE_MPI
     std::shared_ptr<DomainDecomposition> m_decomposition;
-
-    //! Checks neighboring domains to make sure they are properly overlapping
-    void checkDomainBoundaries();
 #endif // ENABLE_MPI
     };
     } // end namespace mpcd
diff --git a/hoomd/mpcd/CellListGPU.cc b/hoomd/mpcd/CellListGPU.cc
index 3b0cee9a2f..f4d75ec72b 100644
--- a/hoomd/mpcd/CellListGPU.cc
+++ b/hoomd/mpcd/CellListGPU.cc
@@ -35,6 +35,30 @@ mpcd::CellListGPU::CellListGPU(std::shared_ptr<SystemDefinition> sysdef,
 #endif // ENABLE_MPI
     }
 
+mpcd::CellListGPU::CellListGPU(std::shared_ptr<SystemDefinition> sysdef,
+                               const uint3& global_cell_dim,
+                               bool shift)
+    : mpcd::CellList(sysdef, global_cell_dim, shift)
+    {
+    m_tuner_cell.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)},
+                                        m_exec_conf,
+                                        "mpcd_cell"));
+    m_tuner_sort.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)},
+                                        m_exec_conf,
+                                        "mpcd_cell_sort"));
+    m_autotuners.insert(m_autotuners.end(), {m_tuner_cell, m_tuner_sort});
+
+#ifdef ENABLE_MPI
+    m_tuner_embed_migrate.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)},
+                                                 m_exec_conf,
+                                                 "mpcd_cell_embed_migrate"));
+    m_autotuners.push_back(m_tuner_embed_migrate);
+
+    GPUFlags<unsigned int> migrate_flag(m_exec_conf);
+    m_migrate_flag.swap(migrate_flag);
+#endif // ENABLE_MPI
+    }
+
 mpcd::CellListGPU::~CellListGPU() { }
 
 void mpcd::CellListGPU::buildCellList()
@@ -90,9 +114,9 @@ void mpcd::CellListGPU::buildCellList()
                                      m_pdata->getBox().getPeriodic(),
                                      m_origin_idx,
                                      m_grid_shift,
-                                     m_pdata->getGlobalBox().getLo(),
+                                     m_pdata->getGlobalBox(),
                                      n_global_cells,
-                                     m_cell_size,
+                                     m_global_cell_dim,
                                      m_cell_np_max,
                                      m_cell_indexer,
                                      m_cell_list_indexer,
@@ -117,9 +141,9 @@ void mpcd::CellListGPU::buildCellList()
                                      m_pdata->getBox().getPeriodic(),
                                      m_origin_idx,
                                      m_grid_shift,
-                                     m_pdata->getGlobalBox().getLo(),
+                                     m_pdata->getGlobalBox(),
                                      n_global_cells,
-                                     m_cell_size,
+                                     m_global_cell_dim,
                                      m_cell_np_max,
                                      m_cell_indexer,
                                      m_cell_list_indexer,
diff --git a/hoomd/mpcd/CellListGPU.cu b/hoomd/mpcd/CellListGPU.cu
index f1e511f79d..146b5f0483 100644
--- a/hoomd/mpcd/CellListGPU.cu
+++ b/hoomd/mpcd/CellListGPU.cu
@@ -29,9 +29,9 @@ namespace kernel
  * \param periodic Flags if local simulation is periodic
  * \param origin_idx Global origin index for the local box
  * \param grid_shift Random grid shift vector
- * \param global_lo Lower bound of global orthorhombic simulation box
+ * \param global_box Global simulation box
  * \param n_global_cell Global dimensions of the cell list, including padding
- * \param cell_size Cell width
+ * \param global_cell_dim Global cell dimensions, no padding
  * \param cell_np_max Maximum number of particles per cell
  * \param cell_indexer 3D indexer for cell id
  * \param cell_list_indexer 2D indexer for particle position in cell
@@ -56,9 +56,9 @@ __global__ void compute_cell_list(unsigned int* d_cell_np,
                                   const uchar3 periodic,
                                   const int3 origin_idx,
                                   const Scalar3 grid_shift,
-                                  const Scalar3 global_lo,
+                                  const BoxDim global_box,
                                   const uint3 n_global_cell,
-                                  const Scalar cell_size,
+                                  const uint3 global_cell_dim,
                                   const unsigned int cell_np_max,
                                   const Index3D cell_indexer,
                                   const Index2D cell_list_indexer,
@@ -87,11 +87,11 @@ __global__ void compute_cell_list(unsigned int* d_cell_np,
         return;
         }
 
-    // bin particle with grid shift assuming orthorhombic box (already validated)
-    const Scalar3 delta = (pos_i - grid_shift) - global_lo;
-    int3 global_bin = make_int3(std::floor(delta.x / cell_size),
-                                std::floor(delta.y / cell_size),
-                                std::floor(delta.z / cell_size));
+    // bin particle with grid shift
+    const Scalar3 fractional_pos_i = global_box.makeFraction(pos_i) - grid_shift;
+    int3 global_bin = make_int3((int)std::floor(fractional_pos_i.x * global_cell_dim.x),
+                                (int)std::floor(fractional_pos_i.y * global_cell_dim.y),
+                                (int)std::floor(fractional_pos_i.z * global_cell_dim.z));
 
     // wrap cell back through the boundaries (grid shifting may send +/- 1 outside of range)
     // this is done using periodic from the "local" box, since this will be periodic
@@ -122,6 +122,28 @@ __global__ void compute_cell_list(unsigned int* d_cell_np,
     int3 bin = make_int3(global_bin.x - origin_idx.x,
                          global_bin.y - origin_idx.y,
                          global_bin.z - origin_idx.z);
+    // these checks guard against round-off errors with domain decomposition
+    if (!periodic.x)
+        {
+        if (bin.x == -1)
+            bin.x = 0;
+        else if (bin.x == (int)cell_indexer.getW())
+            bin.x = cell_indexer.getW() - 1;
+        }
+    if (!periodic.y)
+        {
+        if (bin.y == -1)
+            bin.y = 0;
+        else if (bin.y == (int)cell_indexer.getH())
+            bin.y = cell_indexer.getH() - 1;
+        }
+    if (!periodic.z)
+        {
+        if (bin.z == -1)
+            bin.z = 0;
+        else if (bin.z == (int)cell_indexer.getD())
+            bin.z = cell_indexer.getD() - 1;
+        }
 
     // validate and make sure no particles blew out of the box
     if ((bin.x < 0 || bin.x >= (int)cell_indexer.getW())
@@ -185,13 +207,12 @@ __global__ void cell_check_migrate_embed(unsigned int* d_migrate_flag,
     const Scalar4 postype = d_pos[idx];
     const Scalar3 pos = make_scalar3(postype.x, postype.y, postype.z);
 
-    const Scalar3 lo = box.getLo();
-    const Scalar3 hi = box.getHi();
     const uchar3 periodic = box.getPeriodic();
-
-    if ((!periodic.x && (pos.x >= hi.x || pos.x < lo.x))
-        || (!periodic.y && (pos.y >= hi.y || pos.y < lo.y))
-        || (!periodic.z && num_dim == 3 && (pos.z >= hi.z || pos.z < lo.z)))
+    const Scalar3 fractional_pos = box.makeFraction(pos);
+    if ((!periodic.x && (fractional_pos.x >= Scalar(1.0) || fractional_pos.x < Scalar(0.0)))
+        || (!periodic.y && (fractional_pos.y >= Scalar(1.0) || fractional_pos.y < Scalar(0.0)))
+        || (!periodic.z && num_dim == 3
+            && (fractional_pos.z >= Scalar(1.0) || fractional_pos.z < Scalar(0.0))))
         {
         atomicMax(d_migrate_flag, 1);
         }
@@ -242,9 +263,9 @@ __global__ void cell_apply_sort(unsigned int* d_cell_list,
  * \param periodic Flags if local simulation is periodic
  * \param origin_idx Global origin index for the local box
  * \param grid_shift Random grid shift vector
- * \param global_lo Lower bound of global orthorhombic simulation box
+ * \param global_box Global simulation box
  * \param n_global_cell Global dimensions of the cell list, including padding
- * \param cell_size Cell width
+ * \param global_cell_dim Global cell dimensions, no padding
  * \param cell_np_max Maximum number of particles per cell
  * \param cell_indexer 3D indexer for cell id
  * \param cell_list_indexer 2D indexer for particle position in cell
@@ -265,9 +286,9 @@ cudaError_t mpcd::gpu::compute_cell_list(unsigned int* d_cell_np,
                                          const uchar3& periodic,
                                          const int3& origin_idx,
                                          const Scalar3& grid_shift,
-                                         const Scalar3& global_lo,
+                                         const BoxDim& global_box,
                                          const uint3& n_global_cell,
-                                         const Scalar cell_size,
+                                         const uint3& global_cell_dim,
                                          const unsigned int cell_np_max,
                                          const Index3D& cell_indexer,
                                          const Index2D& cell_list_indexer,
@@ -299,9 +320,9 @@ cudaError_t mpcd::gpu::compute_cell_list(unsigned int* d_cell_np,
                                                                    periodic,
                                                                    origin_idx,
                                                                    grid_shift,
-                                                                   global_lo,
+                                                                   global_box,
                                                                    n_global_cell,
-                                                                   cell_size,
+                                                                   global_cell_dim,
                                                                    cell_np_max,
                                                                    cell_indexer,
                                                                    cell_list_indexer,
diff --git a/hoomd/mpcd/CellListGPU.cuh b/hoomd/mpcd/CellListGPU.cuh
index dec9bf8a5c..3e84d9eefd 100644
--- a/hoomd/mpcd/CellListGPU.cuh
+++ b/hoomd/mpcd/CellListGPU.cuh
@@ -33,9 +33,9 @@ cudaError_t compute_cell_list(unsigned int* d_cell_np,
                               const uchar3& periodic,
                               const int3& origin_idx,
                               const Scalar3& grid_shift,
-                              const Scalar3& global_lo,
+                              const BoxDim& global_box,
                               const uint3& n_global_cell,
-                              const Scalar cell_size,
+                              const uint3& global_cell_dim,
                               const unsigned int cell_np_max,
                               const Index3D& cell_indexer,
                               const Index2D& cell_list_indexer,
diff --git a/hoomd/mpcd/CellListGPU.h b/hoomd/mpcd/CellListGPU.h
index 24e3f00530..f468988778 100644
--- a/hoomd/mpcd/CellListGPU.h
+++ b/hoomd/mpcd/CellListGPU.h
@@ -24,9 +24,12 @@ namespace mpcd
 class PYBIND11_EXPORT CellListGPU : public mpcd::CellList
     {
     public:
-    //! Constructor
+    //! Constructor by size (deprecated)
     CellListGPU(std::shared_ptr<SystemDefinition> sysdef, Scalar cell_size, bool shift);
 
+    //! Constructor by dimension
+    CellListGPU(std::shared_ptr<SystemDefinition> sysdef, const uint3& global_cell_dim, bool shift);
+
     virtual ~CellListGPU();
 
     protected:
diff --git a/hoomd/mpcd/Communicator.cc b/hoomd/mpcd/Communicator.cc
index 1ccf25f28e..a83f7e852a 100644
--- a/hoomd/mpcd/Communicator.cc
+++ b/hoomd/mpcd/Communicator.cc
@@ -460,27 +460,26 @@ void mpcd::Communicator::setCommFlags(const BoxDim& box)
                                           access_mode::overwrite);
 
     // since box is orthorhombic, just use branching to compute comm flags
-    const Scalar3 lo = box.getLo();
-    const Scalar3 hi = box.getHi();
     for (unsigned int idx = 0; idx < N; ++idx)
         {
         const Scalar4& postype = h_pos.data[idx];
         const Scalar3 pos = make_scalar3(postype.x, postype.y, postype.z);
+        const Scalar3 f = box.makeFraction(pos);
 
         unsigned int flags = 0;
-        if (pos.x >= hi.x)
+        if (f.x >= Scalar(1.0))
             flags |= static_cast<unsigned int>(mpcd::detail::send_mask::east);
-        else if (pos.x < lo.x)
+        else if (f.x < Scalar(0.0))
             flags |= static_cast<unsigned int>(mpcd::detail::send_mask::west);
 
-        if (pos.y >= hi.y)
+        if (f.y >= Scalar(1.0))
             flags |= static_cast<unsigned int>(mpcd::detail::send_mask::north);
-        else if (pos.y < lo.y)
+        else if (f.y < Scalar(0.0))
             flags |= static_cast<unsigned int>(mpcd::detail::send_mask::south);
 
-        if (pos.z >= hi.z)
+        if (f.z >= Scalar(1.0))
             flags |= static_cast<unsigned int>(mpcd::detail::send_mask::up);
-        else if (pos.z < lo.z)
+        else if (f.z < Scalar(0.0))
             flags |= static_cast<unsigned int>(mpcd::detail::send_mask::down);
 
         h_comm_flag.data[idx] = flags;
@@ -654,7 +653,10 @@ BoxDim mpcd::Communicator::getWrapBox(const BoxDim& box)
     periodic.y = (isCommunicating(mpcd::detail::face::north)) ? 1 : 0;
     periodic.z = (isCommunicating(mpcd::detail::face::up)) ? 1 : 0;
 
-    return BoxDim(global_lo + shift, global_hi + shift, periodic);
+    BoxDim wrap_box(global_box);
+    wrap_box.setLoHi(global_lo + shift, global_hi + shift);
+    wrap_box.setPeriodic(periodic);
+    return wrap_box;
     }
 
 void mpcd::Communicator::attachCallbacks()
diff --git a/hoomd/mpcd/CommunicatorGPU.cu b/hoomd/mpcd/CommunicatorGPU.cu
index 1cdc338f0d..91c5ff195f 100644
--- a/hoomd/mpcd/CommunicatorGPU.cu
+++ b/hoomd/mpcd/CommunicatorGPU.cu
@@ -50,21 +50,22 @@ stage_particles(unsigned int* d_comm_flag, const Scalar4* d_pos, unsigned int N,
 
     const Scalar4 postype = d_pos[idx];
     const Scalar3 pos = make_scalar3(postype.x, postype.y, postype.z);
-    const Scalar3 lo = box.getLo();
-    const Scalar3 hi = box.getHi();
+    const Scalar3 f = box.makeFraction(pos);
 
     unsigned int flags = 0;
-    if (pos.x >= hi.x)
+    if (f.x >= Scalar(1.0))
         flags |= static_cast<unsigned int>(mpcd::detail::send_mask::east);
-    else if (pos.x < lo.x)
+    else if (f.x < Scalar(0.0))
         flags |= static_cast<unsigned int>(mpcd::detail::send_mask::west);
-    if (pos.y >= hi.y)
+
+    if (f.y >= Scalar(1.0))
         flags |= static_cast<unsigned int>(mpcd::detail::send_mask::north);
-    else if (pos.y < lo.y)
+    else if (f.y < Scalar(0.0))
         flags |= static_cast<unsigned int>(mpcd::detail::send_mask::south);
-    if (pos.z >= hi.z)
+
+    if (f.z >= Scalar(1.0))
         flags |= static_cast<unsigned int>(mpcd::detail::send_mask::up);
-    else if (pos.z < lo.z)
+    else if (f.z < Scalar(0.0))
         flags |= static_cast<unsigned int>(mpcd::detail::send_mask::down);
 
     d_comm_flag[idx] = flags;
diff --git a/hoomd/mpcd/ParallelPlateGeometryFiller.cc b/hoomd/mpcd/ParallelPlateGeometryFiller.cc
index 149268f9e6..30101e9b55 100644
--- a/hoomd/mpcd/ParallelPlateGeometryFiller.cc
+++ b/hoomd/mpcd/ParallelPlateGeometryFiller.cc
@@ -28,10 +28,21 @@ mpcd::ParallelPlateGeometryFiller::~ParallelPlateGeometryFiller()
     m_exec_conf->msg->notice(5) << "Destroying MPCD ParallelPlateGeometryFiller" << std::endl;
     }
 
+void mpcd::ParallelPlateGeometryFiller::fill(uint64_t timestep)
+    {
+    const BoxDim& box = m_pdata->getBox();
+    if (box.getTiltFactorXY() != Scalar(0.0) || box.getTiltFactorXZ() != Scalar(0.0)
+        || box.getTiltFactorYZ() != Scalar(0.0))
+        {
+        throw std::runtime_error("ParallelPlateGeometryFiller does not work with skewed boxes");
+        }
+    mpcd::ManualVirtualParticleFiller::fill(timestep);
+    }
+
 void mpcd::ParallelPlateGeometryFiller::computeNumFill()
     {
     const BoxDim& global_box = m_pdata->getGlobalBox();
-    const Scalar cell_size = m_cl->getCellSize();
+    const Scalar cell_size = m_cl->getCellSize().y;
 
     // box and slit geometry
     const BoxDim& box = m_pdata->getBox();
@@ -49,7 +60,7 @@ void mpcd::ParallelPlateGeometryFiller::computeNumFill()
      * This is done by round the walls onto the cell grid away from zero, and then including the
      * max shift of this cell edge.
      */
-    const Scalar max_shift = m_cl->getMaxGridShift();
+    const Scalar max_shift = m_cl->getMaxGridShift().y * global_box.getL().y;
     const Scalar global_lo = global_box.getLo().y;
     if (box.getHi().y >= H)
         {
diff --git a/hoomd/mpcd/ParallelPlateGeometryFiller.h b/hoomd/mpcd/ParallelPlateGeometryFiller.h
index c4d45ccabd..0389bd0119 100644
--- a/hoomd/mpcd/ParallelPlateGeometryFiller.h
+++ b/hoomd/mpcd/ParallelPlateGeometryFiller.h
@@ -38,6 +38,8 @@ class PYBIND11_EXPORT ParallelPlateGeometryFiller : public mpcd::ManualVirtualPa
 
     virtual ~ParallelPlateGeometryFiller();
 
+    virtual void fill(uint64_t timestep) override;
+
     std::shared_ptr<const mpcd::ParallelPlateGeometry> getGeometry() const
         {
         return m_geom;
diff --git a/hoomd/mpcd/PlanarPoreGeometryFiller.cc b/hoomd/mpcd/PlanarPoreGeometryFiller.cc
index 31d9210df5..7a836e1b5e 100644
--- a/hoomd/mpcd/PlanarPoreGeometryFiller.cc
+++ b/hoomd/mpcd/PlanarPoreGeometryFiller.cc
@@ -29,7 +29,7 @@ mpcd::PlanarPoreGeometryFiller::PlanarPoreGeometryFiller(
 
     // unphysical values in cache to always force recompute
     m_needs_recompute = true;
-    m_recompute_cache = make_scalar3(-1, -1, -1);
+    m_recompute_cache = make_scalar4(-1, -1, -1, -1);
     m_pdata->getBoxChangeSignal()
         .connect<mpcd::PlanarPoreGeometryFiller, &mpcd::PlanarPoreGeometryFiller::notifyRecompute>(
             this);
@@ -43,23 +43,39 @@ mpcd::PlanarPoreGeometryFiller::~PlanarPoreGeometryFiller()
                     &mpcd::PlanarPoreGeometryFiller::notifyRecompute>(this);
     }
 
+void mpcd::PlanarPoreGeometryFiller::fill(uint64_t timestep)
+    {
+    const BoxDim& box = m_pdata->getBox();
+    if (box.getTiltFactorXY() != Scalar(0.0) || box.getTiltFactorXZ() != Scalar(0.0)
+        || box.getTiltFactorYZ() != Scalar(0.0))
+        {
+        throw std::runtime_error("PlanarPoreGeometryFiller does not work with skewed boxes");
+        }
+    mpcd::ManualVirtualParticleFiller::fill(timestep);
+    }
+
 void mpcd::PlanarPoreGeometryFiller::computeNumFill()
     {
-    const Scalar cell_size = m_cl->getCellSize();
-    const Scalar max_shift = m_cl->getMaxGridShift();
+    const Scalar3 cell_size_vector = m_cl->getCellSize();
+    if (fabs(cell_size_vector.x - cell_size_vector.y) > Scalar(1e-6)
+        || fabs(cell_size_vector.x - cell_size_vector.z) > Scalar(1e-6))
+        {
+        throw std::runtime_error("Cell size must be constant");
+        }
+    const Scalar cell_size = cell_size_vector.y;
+    const BoxDim& global_box = m_pdata->getGlobalBox();
+    const Scalar3 max_shift = m_cl->getMaxGridShift() * global_box.getL();
 
     // check if fill-relevant variables have changed (can't use signal because cell list build may
     // not have triggered yet)
-    m_needs_recompute |= (m_recompute_cache.x != cell_size || m_recompute_cache.y != max_shift
-                          || m_recompute_cache.z != m_density);
+    m_needs_recompute
+        |= (m_recompute_cache.x != cell_size || m_recompute_cache.y != max_shift.x
+            || m_recompute_cache.z != max_shift.y || m_recompute_cache.w != m_density);
 
     // only recompute if needed
     if (!m_needs_recompute)
         return;
 
-    // as a precaution, validate the global box with the current cell list
-    const BoxDim& global_box = m_pdata->getGlobalBox();
-
     // box and slit geometry
     const BoxDim& box = m_pdata->getBox();
     const Scalar3 lo = box.getLo();
@@ -77,13 +93,13 @@ void mpcd::PlanarPoreGeometryFiller::computeNumFill()
     const Scalar3 global_hi = global_box.getHi();
     Scalar2 x_bounds, y_bounds;
     // upper bound on lower wall in x
-    x_bounds.x = cell_size * std::ceil((-L - global_lo.x) / cell_size) + global_lo.x + max_shift;
+    x_bounds.x = cell_size * std::ceil((-L - global_lo.x) / cell_size) + global_lo.x + max_shift.x;
     // lower bound on upper wall in x
-    x_bounds.y = cell_size * std::floor((L - global_lo.x) / cell_size) + global_lo.x - max_shift;
+    x_bounds.y = cell_size * std::floor((L - global_lo.x) / cell_size) + global_lo.x - max_shift.x;
     // lower bound on lower wall in y
-    y_bounds.x = cell_size * std::floor((-H - global_lo.y) / cell_size) + global_lo.y - max_shift;
+    y_bounds.x = cell_size * std::floor((-H - global_lo.y) / cell_size) + global_lo.y - max_shift.y;
     // upper bound on upper wall in y
-    y_bounds.y = cell_size * std::ceil((H - global_lo.y) / cell_size) + global_lo.y + max_shift;
+    y_bounds.y = cell_size * std::ceil((H - global_lo.y) / cell_size) + global_lo.y + max_shift.y;
 
     // define the 6 2D bounding boxes (lo.x,hi.x,lo.y,hi.y) for filling in this geometry (z is
     // infinite) (this is essentially a U shape inside the pore).
@@ -130,7 +146,7 @@ void mpcd::PlanarPoreGeometryFiller::computeNumFill()
 
     // size is now updated, cache the cell dimensions used
     m_needs_recompute = false;
-    m_recompute_cache = make_scalar3(cell_size, max_shift, m_density);
+    m_recompute_cache = make_scalar4(cell_size, max_shift.x, max_shift.y, m_density);
     }
 
 /*!
diff --git a/hoomd/mpcd/PlanarPoreGeometryFiller.h b/hoomd/mpcd/PlanarPoreGeometryFiller.h
index 36ef35a1e9..fa6b96f6d1 100644
--- a/hoomd/mpcd/PlanarPoreGeometryFiller.h
+++ b/hoomd/mpcd/PlanarPoreGeometryFiller.h
@@ -38,6 +38,8 @@ class PYBIND11_EXPORT PlanarPoreGeometryFiller : public mpcd::ManualVirtualParti
 
     virtual ~PlanarPoreGeometryFiller();
 
+    virtual void fill(uint64_t timestep) override;
+
     std::shared_ptr<const mpcd::PlanarPoreGeometry> getGeometry() const
         {
         return m_geom;
@@ -65,7 +67,7 @@ class PYBIND11_EXPORT PlanarPoreGeometryFiller : public mpcd::ManualVirtualParti
 
     private:
     bool m_needs_recompute;
-    Scalar3 m_recompute_cache;
+    Scalar4 m_recompute_cache;
     void notifyRecompute()
         {
         m_needs_recompute = true;
diff --git a/hoomd/mpcd/RejectionVirtualParticleFiller.h b/hoomd/mpcd/RejectionVirtualParticleFiller.h
index 3c9b15b070..fdedb2bf7f 100644
--- a/hoomd/mpcd/RejectionVirtualParticleFiller.h
+++ b/hoomd/mpcd/RejectionVirtualParticleFiller.h
@@ -80,8 +80,14 @@ class PYBIND11_EXPORT RejectionVirtualParticleFiller : public mpcd::VirtualParti
 
 template<class Geometry> void RejectionVirtualParticleFiller<Geometry>::fill(uint64_t timestep)
     {
-    // Number of particles that we need to draw (constant)
     const BoxDim& box = m_pdata->getBox();
+    if (box.getTiltFactorXY() != Scalar(0.0) || box.getTiltFactorXZ() != Scalar(0.0)
+        || box.getTiltFactorYZ() != Scalar(0.0))
+        {
+        throw std::runtime_error("Rejection particle filler does not work with skewed boxes");
+        }
+
+    // Number of particles that we need to draw (constant)
     const Scalar3 lo = box.getLo();
     const Scalar3 hi = box.getHi();
     const unsigned int num_virtual_max
diff --git a/hoomd/mpcd/RejectionVirtualParticleFillerGPU.h b/hoomd/mpcd/RejectionVirtualParticleFillerGPU.h
index e28f3afb67..0fc84c1bdb 100644
--- a/hoomd/mpcd/RejectionVirtualParticleFillerGPU.h
+++ b/hoomd/mpcd/RejectionVirtualParticleFillerGPU.h
@@ -66,8 +66,14 @@ class PYBIND11_EXPORT RejectionVirtualParticleFillerGPU
 
 template<class Geometry> void RejectionVirtualParticleFillerGPU<Geometry>::fill(uint64_t timestep)
     {
-    // Number of particles that we need to draw (constant)
     const BoxDim& box = this->m_pdata->getBox();
+    if (box.getTiltFactorXY() != Scalar(0.0) || box.getTiltFactorXZ() != Scalar(0.0)
+        || box.getTiltFactorYZ() != Scalar(0.0))
+        {
+        throw std::runtime_error("Rejection particle filler does not work with skewed boxes");
+        }
+
+    // Number of particles that we need to draw (constant)
     const Scalar3 lo = box.getLo();
     const Scalar3 hi = box.getHi();
     const unsigned int num_virtual_max
diff --git a/hoomd/mpcd/collide.py b/hoomd/mpcd/collide.py
index 94a1a72f4c..1f8d96ba96 100644
--- a/hoomd/mpcd/collide.py
+++ b/hoomd/mpcd/collide.py
@@ -29,8 +29,10 @@ class CellList(Compute):
     Args:
         shift (bool): When True, randomly shift underlying collision cells.
 
-    The MPCD `CellList` bins particles into cubic cells of edge length 1.0.
-    The simulation box must be orthorhombic, and its edges must be an integer.
+    The MPCD `CellList` bins particles into cells aligned with the lattice
+    vectors that define the simulation box. The default number of cells along
+    each lattice vector is computed from :math:`L_x`, :math:`L_y`, and
+    :math:`L_z` with a length of 1.0.
 
     When the total mean-free path of the MPCD particles is small, the cells
     should be randomly shifted in order to ensure Galilean invariance of the
@@ -85,6 +87,11 @@ def _attach_hook(self):
 
         super()._attach_hook()
 
+    @property
+    def num_cells(self):
+        """tuple[int]: Number of cells along each lattice vector."""
+        return self._cpp_obj.num_cells
+
 
 class CollisionMethod(Operation):
     """Base collision method.
diff --git a/hoomd/mpcd/fill.py b/hoomd/mpcd/fill.py
index 4b6b2438d0..b3b5e97929 100644
--- a/hoomd/mpcd/fill.py
+++ b/hoomd/mpcd/fill.py
@@ -131,6 +131,13 @@ class GeometryFiller(VirtualParticleFiller):
     specified `geometry`. The algorithm for doing the filling depends on the
     specific `geometry`.
 
+    .. rubric:: Limitations:
+
+    This filler **does not** currently support triclinic boxes for any
+    :class:`~hoomd.mpcd.geometry.Geometry`. Additionally, this filler does not
+    support the :class:`~hoomd.mpcd.geometry.PlanarPore` geometry for any
+    non-cubic cell shape. Exceptions will be raised in these cases.
+
     .. rubric:: Example:
 
     Filler for parallel plate geometry.
diff --git a/hoomd/mpcd/pytest/test_collide.py b/hoomd/mpcd/pytest/test_collide.py
index 0db7c21266..04d3a1dddc 100644
--- a/hoomd/mpcd/pytest/test_collide.py
+++ b/hoomd/mpcd/pytest/test_collide.py
@@ -19,6 +19,17 @@ def small_snap():
     return snap
 
 
+def test_cell_list(small_snap, simulation_factory):
+    if small_snap.communicator.rank == 0:
+        small_snap.configuration.box = [20, 30, 40, 0, 0, 0]
+    sim = simulation_factory(small_snap)
+
+    cl = hoomd.mpcd.collide.CellList()
+
+    cl._attach(sim)
+    assert cl.num_cells == (20, 30, 40)
+
+
 @pytest.mark.parametrize(
     "cls, init_args",
     [
diff --git a/hoomd/mpcd/test/cell_communicator_mpi_test.cc b/hoomd/mpcd/test/cell_communicator_mpi_test.cc
index f7e105ad4a..130afd1bb2 100644
--- a/hoomd/mpcd/test/cell_communicator_mpi_test.cc
+++ b/hoomd/mpcd/test/cell_communicator_mpi_test.cc
@@ -155,7 +155,7 @@ void cell_communicator_overdecompose_test(std::shared_ptr<ExecutionConfiguration
     std::shared_ptr<Communicator> pdata_comm(new Communicator(sysdef, decomposition));
     sysdef->setCommunicator(pdata_comm);
 
-    auto cl = std::make_shared<mpcd::CellList>(sysdef, 1.0, false);
+    auto cl = std::make_shared<mpcd::CellList>(sysdef, make_uint3(6, 6, 6), false);
     cl->computeDimensions();
 
     // Don't really care what's in this array, just want to make sure errors get thrown
@@ -191,20 +191,20 @@ void cell_communicator_overdecompose_test(std::shared_ptr<ExecutionConfiguration
     UP_ASSERT_EXCEPTION(std::runtime_error, [&] { comm.communicate(props, pack_op); });
 
     // cut the cell size down, which should make it so that the system can be decomposed again
-    cl->setCellSize(0.5);
+    cl->setGlobalDim(make_uint3(12, 12, 12));
     cl->compute(3);
     props.resize(cl->getNCells());
     comm.communicate(props, pack_op);
 
     // shrink the box size to the minimum that can be decomposed
     cl->setNExtraCells(0);
-    cl->setCellSize(2.0);
+    cl->setGlobalDim(make_uint3(3, 3, 3));
     cl->compute(4);
     props.resize(cl->getNCells());
     comm.communicate(props, pack_op);
 
     // now shrink further to ensure failure
-    cl->setCellSize(3.0);
+    cl->setGlobalDim(make_uint3(2, 2, 2));
     cl->compute(5);
     props.resize(cl->getNCells());
     UP_ASSERT_EXCEPTION(std::runtime_error, [&] { comm.communicate(props, pack_op); });
diff --git a/hoomd/mpcd/test/cell_list_mpi_test.cc b/hoomd/mpcd/test/cell_list_mpi_test.cc
index 3af0996f66..43205a03e9 100644
--- a/hoomd/mpcd/test/cell_list_mpi_test.cc
+++ b/hoomd/mpcd/test/cell_list_mpi_test.cc
@@ -11,10 +11,145 @@
 #include "hoomd/mpcd/Communicator.h"
 #include "hoomd/test/upp11_config.h"
 
+#include "utils.h"
+
 HOOMD_UP_MAIN()
 
 using namespace hoomd;
 
+void checkDomainBoundaries(std::shared_ptr<SystemDefinition> sysdef,
+                           std::shared_ptr<mpcd::CellList> cl)
+    {
+    auto pdata = sysdef->getParticleData();
+    auto exec_conf = pdata->getExecConf();
+    MPI_Comm mpi_comm = exec_conf->getMPICommunicator();
+    auto decomposition = pdata->getDomainDecomposition();
+
+    auto num_comm = cl->getNComm();
+    auto origin_idx = cl->getOriginIndex();
+    auto cell_dim = cl->getDim();
+    auto num_extra = cl->getNExtraCells();
+
+    for (unsigned int dir = 0; dir < num_comm.size(); ++dir)
+        {
+        mpcd::detail::face d = static_cast<mpcd::detail::face>(dir);
+        if (!cl->isCommunicating(d))
+            continue;
+
+        // receive in the opposite direction from which we send
+        unsigned int send_neighbor = decomposition->getNeighborRank(dir);
+        unsigned int recv_neighbor;
+        if (dir % 2 == 0)
+            recv_neighbor = decomposition->getNeighborRank(dir + 1);
+        else
+            recv_neighbor = decomposition->getNeighborRank(dir - 1);
+
+        // first make sure each dimension is sending and receiving the same size data
+        MPI_Request reqs[2];
+        MPI_Status status[2];
+
+        // check that the number received is the same as that being sent from neighbor
+        unsigned int n_send = num_comm[dir];
+        unsigned int n_expect_recv;
+        if (dir % 2 == 0)
+            n_expect_recv = num_comm[dir + 1];
+        else
+            n_expect_recv = num_comm[dir - 1];
+
+        unsigned int n_recv;
+        MPI_Isend(&n_send, 1, MPI_UNSIGNED, send_neighbor, 0, mpi_comm, &reqs[0]);
+        MPI_Irecv(&n_recv, 1, MPI_UNSIGNED, recv_neighbor, 0, mpi_comm, &reqs[1]);
+        MPI_Waitall(2, reqs, status);
+        UP_ASSERT_EQUAL(n_recv, n_expect_recv);
+
+        // check that the same cell ids are communicated
+        std::vector<int> send_cells(n_send), recv_cells(n_recv);
+        for (unsigned int i = 0; i < n_send; ++i)
+            {
+            if (d == mpcd::detail::face::east)
+                {
+                send_cells[i] = origin_idx.x + cell_dim.x - num_extra - n_send + i;
+                }
+            else if (d == mpcd::detail::face::west)
+                {
+                send_cells[i] = origin_idx.x + i;
+                }
+            else if (d == mpcd::detail::face::north)
+                {
+                send_cells[i] = origin_idx.y + cell_dim.y - num_extra - n_send + i;
+                }
+            else if (d == mpcd::detail::face::south)
+                {
+                send_cells[i] = origin_idx.y + i;
+                }
+            else if (d == mpcd::detail::face::up)
+                {
+                send_cells[i] = origin_idx.z + cell_dim.z - num_extra - n_send + i;
+                }
+            else if (d == mpcd::detail::face::down)
+                {
+                send_cells[i] = origin_idx.z + i;
+                }
+            }
+
+        MPI_Isend(&send_cells[0], n_send, MPI_INT, send_neighbor, 1, mpi_comm, &reqs[0]);
+        MPI_Irecv(&recv_cells[0], n_recv, MPI_INT, recv_neighbor, 1, mpi_comm, &reqs[1]);
+        MPI_Waitall(2, reqs, status);
+
+        for (unsigned int i = 0; i < n_recv; ++i)
+            {
+            // wrap the received cell back into the global box
+            // only two of the entries will be valid, the others are dummies
+            int3 recv_cell = make_int3(0, 0, 0);
+            if (d == mpcd::detail::face::east || d == mpcd::detail::face::west)
+                {
+                recv_cell.x = recv_cells[i];
+                }
+            else if (d == mpcd::detail::face::north || d == mpcd::detail::face::south)
+                {
+                recv_cell.y = recv_cells[i];
+                }
+            else if (d == mpcd::detail::face::up || d == mpcd::detail::face::down)
+                {
+                recv_cell.z = recv_cells[i];
+                }
+            recv_cell = cl->wrapGlobalCell(recv_cell);
+
+            // compute the expected cell to receive, also wrapped
+            int3 expect_recv_cell = make_int3(0, 0, 0);
+            if (d == mpcd::detail::face::east)
+                {
+                expect_recv_cell.x = origin_idx.x + i;
+                }
+            else if (d == mpcd::detail::face::west)
+                {
+                expect_recv_cell.x = origin_idx.x + cell_dim.x - num_extra - n_recv + i;
+                }
+            else if (d == mpcd::detail::face::north)
+                {
+                expect_recv_cell.y = origin_idx.y + i;
+                }
+            else if (d == mpcd::detail::face::south)
+                {
+                expect_recv_cell.y = origin_idx.y + cell_dim.y - num_extra - n_recv + i;
+                }
+            else if (d == mpcd::detail::face::up)
+                {
+                expect_recv_cell.z = origin_idx.z + i;
+                }
+            else if (d == mpcd::detail::face::down)
+                {
+                expect_recv_cell.z = origin_idx.z + cell_dim.z - num_extra - n_recv + i;
+                }
+            expect_recv_cell = cl->wrapGlobalCell(expect_recv_cell);
+
+            UP_ASSERT_EQUAL(recv_cell.x, expect_recv_cell.x);
+            UP_ASSERT_EQUAL(recv_cell.y, expect_recv_cell.y);
+            UP_ASSERT_EQUAL(recv_cell.z, expect_recv_cell.z);
+            }
+        }
+    }
+
 //! Test for correct calculation of MPCD grid dimensions
 /*!
  * \param exec_conf Execution configuration
@@ -28,7 +163,8 @@ template<class CL>
 void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
                              bool mpi_x,
                              bool mpi_y,
-                             bool mpi_z)
+                             bool mpi_z,
+                             const Scalar3& tilt)
     {
     // only run tests on first partition
     if (exec_conf->getPartition() != 0)
@@ -36,6 +172,7 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
 
     std::shared_ptr<SnapshotSystemData<Scalar>> snap(new SnapshotSystemData<Scalar>());
     snap->global_box = std::make_shared<BoxDim>(5.0);
+    snap->global_box->setTiltFactors(tilt.x, tilt.y, tilt.z);
     snap->particle_data.type_mapping.push_back("A");
     snap->mpcd_data.resize(1);
     snap->mpcd_data.type_mapping.push_back("A");
@@ -74,11 +211,14 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
 
     // initialize mpcd system
     auto pdata_1 = sysdef->getMPCDParticleData();
-    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, 1.0, false));
+    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, make_uint3(5, 5, 5), false));
 
     // compute the cell list
     cl->computeDimensions();
+    checkDomainBoundaries(sysdef, cl);
 
+    const bool is_orthorhombic = (tilt == make_scalar3(0, 0, 0));
+    if (is_orthorhombic)
         {
         // check domain origins
         const int3 origin = cl->getOriginIndex();
@@ -139,7 +279,14 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
         if (mpi_x)
             {
             // split evenly in x -> both domains are same size
-            UP_ASSERT_EQUAL(dim.x, 4);
+            if (pos.x)
+                {
+                UP_ASSERT_EQUAL(dim.x, 4);
+                }
+            else
+                {
+                UP_ASSERT_EQUAL(dim.x, 4);
+                }
             }
         else
             {
@@ -181,52 +328,18 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
             }
 
         std::array<unsigned int, 6> num_comm = cl->getNComm();
-        if (mpi_x)
-            {
-            if (pos.x)
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)], 2);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)], 1);
-                }
-            else
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)], 1);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)], 2);
-                }
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)], 0);
-            }
-
-        if (mpi_y)
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)], 2);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)], 2);
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)], 0);
-            }
-
-        if (mpi_z)
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)], 2);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)], 2);
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)], 0);
-            }
-
-        // check for grid shifting errors
-        UP_ASSERT_EXCEPTION(std::runtime_error,
-                            [&] { cl->setGridShift(make_scalar3(-0.51, -0.51, -0.51)); });
-        UP_ASSERT_EXCEPTION(std::runtime_error,
-                            [&] { cl->setGridShift(make_scalar3(0.51, 0.51, 0.51)); });
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)],
+                        (mpi_x) ? ((pos.x) ? 2 : 1) : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)],
+                        (mpi_x) ? ((pos.x) ? 1 : 2) : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)],
+                        (mpi_y) ? 2 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)],
+                        (mpi_y) ? 2 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)],
+                        (mpi_z) ? 2 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)],
+                        (mpi_z) ? 2 : 0);
 
         // check coverage box
         const BoxDim coverage = cl->getCoverageBox();
@@ -234,64 +347,66 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
             {
             if (pos.x)
                 {
-                CHECK_CLOSE(coverage.getLo().x, 0.0, tol);
-                CHECK_CLOSE(coverage.getHi().x, 3.0, tol);
+                UP_ASSERT_SMALL(coverage.getLo().x, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().x, 3.0, tol);
                 }
             else
                 {
-                CHECK_CLOSE(coverage.getLo().x, -3.0, tol);
-                CHECK_CLOSE(coverage.getHi().x, 0.0, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().x, -3.0, tol);
+                UP_ASSERT_SMALL(coverage.getHi().x, tol);
                 }
             }
         else
             {
-            CHECK_CLOSE(coverage.getLo().x, -2.5, tol);
-            CHECK_CLOSE(coverage.getHi().x, 2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getLo().x, -2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getHi().x, 2.5, tol);
             }
 
         if (mpi_y)
             {
             if (pos.y)
                 {
-                CHECK_CLOSE(coverage.getLo().y, -1.0, tol);
-                CHECK_CLOSE(coverage.getHi().y, 3.0, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().y, -1.0, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().y, 3.0, tol);
                 }
             else
                 {
-                CHECK_CLOSE(coverage.getLo().y, -3.0, tol);
-                CHECK_CLOSE(coverage.getHi().y, 0.0, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().y, -3.0, tol);
+                UP_ASSERT_SMALL(coverage.getHi().y, tol);
                 }
             }
         else
             {
-            CHECK_CLOSE(coverage.getLo().y, -2.5, tol);
-            CHECK_CLOSE(coverage.getHi().y, 2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getLo().y, -2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getHi().y, 2.5, tol);
             }
 
         if (mpi_z)
             {
             if (pos.z)
                 {
-                CHECK_CLOSE(coverage.getLo().z, 0.0, tol);
-                CHECK_CLOSE(coverage.getHi().z, 3.0, tol);
+                UP_ASSERT_SMALL(coverage.getLo().z, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().z, 3.0, tol);
                 }
             else
                 {
-                CHECK_CLOSE(coverage.getLo().z, -3.0, tol);
-                CHECK_CLOSE(coverage.getHi().z, 1.0, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().z, -3.0, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().z, 1.0, tol);
                 }
             }
         else
             {
-            CHECK_CLOSE(coverage.getLo().z, -2.5, tol);
-            CHECK_CLOSE(coverage.getHi().z, 2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getLo().z, -2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getHi().z, 2.5, tol);
             }
         }
 
     /*******************/
     // Change the cell size, and ensure everything stays up to date
-    cl->setCellSize(0.5);
+    cl->setGlobalDim(make_uint3(10, 10, 10));
     cl->computeDimensions();
+    checkDomainBoundaries(sysdef, cl);
+    if (is_orthorhombic)
         {
         // check domain origins
         const int3 origin = cl->getOriginIndex();
@@ -385,7 +500,9 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
                 }
             else
                 {
-                UP_ASSERT_EQUAL(dim.z, 7);
+                // floating point rounding makes this 8 not 7 (extra cell communicated on this
+                // edge)
+                UP_ASSERT_EQUAL(dim.z, 8);
                 }
             }
         else
@@ -394,124 +511,86 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
             }
 
         std::array<unsigned int, 6> num_comm = cl->getNComm();
-        if (mpi_x)
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)], 2);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)], 2);
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)], 0);
-            }
-
-        if (mpi_y)
-            {
-            if (pos.y)
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)], 2);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)], 1);
-                }
-            else
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)], 1);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)], 2);
-                }
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)], 0);
-            }
-
-        if (mpi_z)
-            {
-            if (pos.z)
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)], 2);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)], 1);
-                }
-            else
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)], 1);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)], 2);
-                }
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)], 0);
-            }
-
-        UP_ASSERT_EXCEPTION(std::runtime_error,
-                            [&] { cl->setGridShift(make_scalar3(-0.3, -0.3, -0.3)); });
-        UP_ASSERT_EXCEPTION(std::runtime_error,
-                            [&] { cl->setGridShift(make_scalar3(0.3, 0.3, 0.3)); });
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)],
+                        (mpi_x) ? 2 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)],
+                        (mpi_x) ? 2 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)],
+                        (mpi_y) ? ((pos.y) ? 2 : 1) : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)],
+                        (mpi_y) ? ((pos.y) ? 1 : 2) : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)],
+                        (mpi_z) ? 2 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)],
+                        (mpi_z) ? 2 : 0);
 
         const BoxDim coverage = cl->getCoverageBox();
         if (mpi_x)
             {
             if (pos.x)
                 {
-                CHECK_CLOSE(coverage.getLo().x, -0.25, tol);
-                CHECK_CLOSE(coverage.getHi().x, 2.75, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().x, -0.25, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().x, 2.75, tol);
                 }
             else
                 {
-                CHECK_CLOSE(coverage.getLo().x, -2.75, tol);
-                CHECK_CLOSE(coverage.getHi().x, 0.25, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().x, -2.75, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().x, 0.25, tol);
                 }
             }
         else
             {
-            CHECK_CLOSE(coverage.getLo().x, -2.5, tol);
-            CHECK_CLOSE(coverage.getHi().x, 2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getLo().x, -2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getHi().x, 2.5, tol);
             }
 
         if (mpi_y)
             {
             if (pos.y)
                 {
-                CHECK_CLOSE(coverage.getLo().y, -0.25, tol);
-                CHECK_CLOSE(coverage.getHi().y, 2.75, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().y, -0.25, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().y, 2.75, tol);
                 }
             else
                 {
-                CHECK_CLOSE(coverage.getLo().y, -2.75, tol);
-                CHECK_CLOSE(coverage.getHi().y, -0.25, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().y, -2.75, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().y, -0.25, tol);
                 }
             }
         else
             {
-            CHECK_CLOSE(coverage.getLo().y, -2.5, tol);
-            CHECK_CLOSE(coverage.getHi().y, 2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getLo().y, -2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getHi().y, 2.5, tol);
             }
 
         if (mpi_z)
             {
             if (pos.z)
                 {
-                CHECK_CLOSE(coverage.getLo().z, 0.25, tol);
-                CHECK_CLOSE(coverage.getHi().z, 2.75, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().z, 0.25, tol);
+                UP_ASSERT_CLOSE(coverage.getHi().z, 2.75, tol);
                 }
             else
                 {
-                CHECK_CLOSE(coverage.getLo().z, -2.75, tol);
-                CHECK_CLOSE(coverage.getHi().z, 0.25, tol);
+                UP_ASSERT_CLOSE(coverage.getLo().z, -2.75, tol);
+                // floating point rounding makes this 0.75 not 0.25
+                UP_ASSERT_CLOSE(coverage.getHi().z, 0.75, tol);
                 }
             }
         else
             {
-            CHECK_CLOSE(coverage.getLo().z, -2.5, tol);
-            CHECK_CLOSE(coverage.getHi().z, 2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getLo().z, -2.5, tol);
+            UP_ASSERT_CLOSE(coverage.getHi().z, 2.5, tol);
             }
         }
 
     /*******************/
-    // Increase the number of communication cells. This will trigger an increase in the size of the
-    // diffusion layer
+    // Increase the number of communication cells. This will trigger an increase in the size of
+    // the diffusion layer
     cl->setNExtraCells(1);
     cl->computeDimensions();
+    checkDomainBoundaries(sysdef, cl);
+    if (is_orthorhombic)
         {
         // all origins should be shifted down by 1 cell
         const int3 origin = cl->getOriginIndex();
@@ -599,7 +678,8 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
                 }
             else
                 {
-                UP_ASSERT_EQUAL(dim.z, 9);
+                // floating point rounding makes thes 10 not 9
+                UP_ASSERT_EQUAL(dim.z, 10);
                 }
             }
         else
@@ -609,59 +689,18 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
 
         // all comms should be increased by 1
         std::array<unsigned int, 6> num_comm = cl->getNComm();
-        if (mpi_x)
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)], 3);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)], 3);
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)], 0);
-            }
-
-        if (mpi_y)
-            {
-            if (pos.y)
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)], 3);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)], 2);
-                }
-            else
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)], 2);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)], 3);
-                }
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)], 0);
-            }
-
-        if (mpi_z)
-            {
-            if (pos.z)
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)], 3);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)], 2);
-                }
-            else
-                {
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)], 2);
-                UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)], 3);
-                }
-            }
-        else
-            {
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)], 0);
-            UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)], 0);
-            }
-
-        UP_ASSERT_EXCEPTION(std::runtime_error,
-                            [&] { cl->setGridShift(make_scalar3(-0.3, -0.3, -0.3)); });
-        UP_ASSERT_EXCEPTION(std::runtime_error,
-                            [&] { cl->setGridShift(make_scalar3(0.3, 0.3, 0.3)); });
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::east)],
+                        (mpi_x) ? 3 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::west)],
+                        (mpi_x) ? 3 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::north)],
+                        (mpi_y) ? ((pos.y) ? 3 : 2) : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::south)],
+                        (mpi_y) ? ((pos.y) ? 2 : 3) : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::up)],
+                        (mpi_z) ? 3 : 0);
+        UP_ASSERT_EQUAL(num_comm[static_cast<unsigned int>(mpcd::detail::face::down)],
+                        (mpi_z) ? 3 : 0);
 
         const BoxDim coverage = cl->getCoverageBox();
         if (mpi_x)
@@ -712,7 +751,8 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
             else
                 {
                 CHECK_CLOSE(coverage.getLo().z, -3.25, tol);
-                CHECK_CLOSE(coverage.getHi().z, 0.75, tol);
+                // floating point rounding makes this 1.25 not 0.75
+                CHECK_CLOSE(coverage.getHi().z, 1.25, tol);
                 }
             }
         else
@@ -724,12 +764,19 @@ void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
     }
 
 //! Test for correct cell listing of a basic system
-template<class CL> void celllist_basic_test(std::shared_ptr<ExecutionConfiguration> exec_conf)
+template<class CL>
+void celllist_basic_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
+                         const Scalar3& L,
+                         const Scalar3& tilt)
     {
     UP_ASSERT_EQUAL(exec_conf->getNRanks(), 8);
 
+    auto ref_box = std::make_shared<BoxDim>(6.0);
+    auto box = std::make_shared<BoxDim>(L);
+    box->setTiltFactors(tilt.x, tilt.y, tilt.z);
+
     std::shared_ptr<SnapshotSystemData<Scalar>> snap(new SnapshotSystemData<Scalar>());
-    snap->global_box = std::make_shared<BoxDim>(6.0);
+    snap->global_box = box;
     snap->particle_data.type_mapping.push_back("A");
     // place each particle in the same cell, but on different ranks
     /*
@@ -746,14 +793,14 @@ template<class CL> void celllist_basic_test(std::shared_ptr<ExecutionConfigurati
      */
     snap->mpcd_data.resize(8);
     snap->mpcd_data.type_mapping.push_back("A");
-    snap->mpcd_data.position[0] = vec3<Scalar>(-0.1, -0.1, -0.1);
-    snap->mpcd_data.position[1] = vec3<Scalar>(0.1, -0.1, -0.1);
-    snap->mpcd_data.position[2] = vec3<Scalar>(-0.1, 0.1, -0.1);
-    snap->mpcd_data.position[3] = vec3<Scalar>(0.1, 0.1, -0.1);
-    snap->mpcd_data.position[4] = vec3<Scalar>(-0.1, -0.1, 0.1);
-    snap->mpcd_data.position[5] = vec3<Scalar>(0.1, -0.1, 0.1);
-    snap->mpcd_data.position[6] = vec3<Scalar>(-0.1, 0.1, 0.1);
-    snap->mpcd_data.position[7] = vec3<Scalar>(0.1, 0.1, 0.1);
+    snap->mpcd_data.position[0] = scale(vec3<Scalar>(-0.1, -0.1, -0.1), ref_box, box);
+    snap->mpcd_data.position[1] = scale(vec3<Scalar>(0.1, -0.1, -0.1), ref_box, box);
+    snap->mpcd_data.position[2] = scale(vec3<Scalar>(-0.1, 0.1, -0.1), ref_box, box);
+    snap->mpcd_data.position[3] = scale(vec3<Scalar>(0.1, 0.1, -0.1), ref_box, box);
+    snap->mpcd_data.position[4] = scale(vec3<Scalar>(-0.1, -0.1, 0.1), ref_box, box);
+    snap->mpcd_data.position[5] = scale(vec3<Scalar>(0.1, -0.1, 0.1), ref_box, box);
+    snap->mpcd_data.position[6] = scale(vec3<Scalar>(-0.1, 0.1, 0.1), ref_box, box);
+    snap->mpcd_data.position[7] = scale(vec3<Scalar>(0.1, 0.1, 0.1), ref_box, box);
 
     std::shared_ptr<DomainDecomposition> decomposition(
         new DomainDecomposition(exec_conf, snap->global_box->getL(), 2, 2, 2));
@@ -763,7 +810,7 @@ template<class CL> void celllist_basic_test(std::shared_ptr<ExecutionConfigurati
 
     // initialize mpcd system
     std::shared_ptr<mpcd::ParticleData> pdata = sysdef->getMPCDParticleData();
-    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, 1.0, false));
+    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, make_uint3(6, 6, 6), false));
     cl->compute(0);
     const unsigned int my_rank = exec_conf->getRank();
         {
@@ -824,7 +871,8 @@ template<class CL> void celllist_basic_test(std::shared_ptr<ExecutionConfigurati
         }
 
     // apply a grid shift so that particles move into the same cell (3,3,3)
-    cl->setGridShift(make_scalar3(-0.5, -0.5, -0.5));
+    const Scalar3 shift = (Scalar(0.5) / 6) * make_scalar3(1, 1, 1);
+    cl->setGridShift(-shift);
     cl->compute(1);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
@@ -884,7 +932,7 @@ template<class CL> void celllist_basic_test(std::shared_ptr<ExecutionConfigurati
         }
 
     // apply a grid shift so that particles move into the same cell (2,2,2)
-    cl->setGridShift(make_scalar3(0.5, 0.5, 0.5));
+    cl->setGridShift(shift);
     cl->compute(2);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
@@ -945,12 +993,20 @@ template<class CL> void celllist_basic_test(std::shared_ptr<ExecutionConfigurati
     }
 
 //! Test for correct cell listing of a system with particles on the edges
-template<class CL> void celllist_edge_test(std::shared_ptr<ExecutionConfiguration> exec_conf)
+template<class CL>
+void celllist_edge_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
+                        const Scalar3& L,
+                        const Scalar3& tilt)
     {
     UP_ASSERT_EQUAL(exec_conf->getNRanks(), 8);
 
+    auto ref_box = std::make_shared<BoxDim>(5.0);
+    auto box = std::make_shared<BoxDim>(L);
+    box->setTiltFactors(tilt.x, tilt.y, tilt.z);
+    bool is_orthorhombic = tilt.x == Scalar(0.0) && tilt.y == Scalar(0.0) && tilt.z == Scalar(0.0);
+
     std::shared_ptr<SnapshotSystemData<Scalar>> snap(new SnapshotSystemData<Scalar>());
-    snap->global_box = std::make_shared<BoxDim>(5.0);
+    snap->global_box = box;
     snap->particle_data.type_mapping.push_back("A");
     // dummy initialize one particle to every domain, we will move them outside the domains for
     // the tests
@@ -968,14 +1024,14 @@ template<class CL> void celllist_edge_test(std::shared_ptr<ExecutionConfiguratio
      */
     snap->mpcd_data.resize(8);
     snap->mpcd_data.type_mapping.push_back("A");
-    snap->mpcd_data.position[0] = vec3<Scalar>(-1.0, -1.0, -1.0);
-    snap->mpcd_data.position[1] = vec3<Scalar>(1.0, -1.0, -1.0);
-    snap->mpcd_data.position[2] = vec3<Scalar>(-1.0, 1.0, -1.0);
-    snap->mpcd_data.position[3] = vec3<Scalar>(1.0, 1.0, -1.0);
-    snap->mpcd_data.position[4] = vec3<Scalar>(-1.0, -1.0, 1.0);
-    snap->mpcd_data.position[5] = vec3<Scalar>(1.0, -1.0, 1.0);
-    snap->mpcd_data.position[6] = vec3<Scalar>(-1.0, 1.0, 1.0);
-    snap->mpcd_data.position[7] = vec3<Scalar>(1.0, 1.0, 1.0);
+    snap->mpcd_data.position[0] = scale(vec3<Scalar>(-1.0, -1.0, -1.0), ref_box, box);
+    snap->mpcd_data.position[1] = scale(vec3<Scalar>(1.0, -1.0, -1.0), ref_box, box);
+    snap->mpcd_data.position[2] = scale(vec3<Scalar>(-1.0, 1.0, -1.0), ref_box, box);
+    snap->mpcd_data.position[3] = scale(vec3<Scalar>(1.0, 1.0, -1.0), ref_box, box);
+    snap->mpcd_data.position[4] = scale(vec3<Scalar>(-1.0, -1.0, 1.0), ref_box, box);
+    snap->mpcd_data.position[5] = scale(vec3<Scalar>(1.0, -1.0, 1.0), ref_box, box);
+    snap->mpcd_data.position[6] = scale(vec3<Scalar>(-1.0, 1.0, 1.0), ref_box, box);
+    snap->mpcd_data.position[7] = scale(vec3<Scalar>(1.0, 1.0, 1.0), ref_box, box);
     std::vector<Scalar> fx {0.5};
     std::vector<Scalar> fy {0.45};
     std::vector<Scalar> fz {0.55};
@@ -986,7 +1042,7 @@ template<class CL> void celllist_edge_test(std::shared_ptr<ExecutionConfiguratio
     sysdef->setCommunicator(pdata_comm);
 
     std::shared_ptr<mpcd::ParticleData> pdata = sysdef->getMPCDParticleData();
-    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, 1.0, false));
+    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, make_uint3(5, 5, 5), false));
 
     // move particles to edges of domains for testing
     const unsigned int my_rank = exec_conf->getRank();
@@ -997,212 +1053,147 @@ template<class CL> void celllist_edge_test(std::shared_ptr<ExecutionConfiguratio
         switch (my_rank)
             {
         case 0:
-            h_pos.data[0] = make_scalar4(-0.01, -0.01, -0.01, __int_as_scalar(0));
+            h_pos.data[0]
+                = scale(make_scalar4(-0.01, -0.01, -0.01, __int_as_scalar(0)), ref_box, box);
             break;
         case 1:
-            h_pos.data[0] = make_scalar4(0.0, -0.01, -0.01, __int_as_scalar(0));
+            h_pos.data[0]
+                = scale(make_scalar4(0.0, -0.01, -0.01, __int_as_scalar(0)), ref_box, box);
             break;
         case 2:
-            h_pos.data[0] = make_scalar4(-0.01, 0.0, -0.01, __int_as_scalar(0));
+            h_pos.data[0]
+                = scale(make_scalar4(-0.01, 0.0, -0.01, __int_as_scalar(0)), ref_box, box);
             break;
         case 3:
-            h_pos.data[0] = make_scalar4(0.0, 0.0, -0.01, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(0.0, 0.0, -0.01, __int_as_scalar(0)), ref_box, box);
             break;
         case 4:
-            h_pos.data[0] = make_scalar4(-0.01, -0.01, 0.0, __int_as_scalar(0));
+            h_pos.data[0]
+                = scale(make_scalar4(-0.01, -0.01, 0.0, __int_as_scalar(0)), ref_box, box);
             break;
         case 5:
-            h_pos.data[0] = make_scalar4(0.0, -0.01, 0.0, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(0.0, -0.01, 0.0, __int_as_scalar(0)), ref_box, box);
             break;
         case 6:
-            h_pos.data[0] = make_scalar4(-0.01, 0.0, 0.0, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(-0.01, 0.0, 0.0, __int_as_scalar(0)), ref_box, box);
             break;
         case 7:
-            h_pos.data[0] = make_scalar4(0.0, 0.0, 0.0, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(0.0, 0.0, 0.0, __int_as_scalar(0)), ref_box, box);
             break;
             };
         }
 
     cl->compute(0);
-
+    checkDomainBoundaries(sysdef, cl);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
                                             access_location::host,
                                             access_mode::read);
-        ArrayHandle<unsigned int> h_cell_list(cl->getCellList(),
-                                              access_location::host,
-                                              access_mode::read);
-        Index3D ci = cl->getCellIndexer();
         ArrayHandle<Scalar4> h_vel(pdata->getVelocities(),
                                    access_location::host,
                                    access_mode::read);
-
-        switch (my_rank)
-            {
-        case 0:
-            // global index is (2,2,2), with origin (-1,-1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(3, 3, 3)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(3, 3, 3));
-            break;
-        case 1:
-            // global index is (2,2,2), with origin (2,-1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 3, 3)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 3, 3));
-            break;
-        case 2:
-            // global index is (2,2,2), with origin (-1,1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(3, 1, 3)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(3, 1, 3));
-            break;
-        case 3:
-            // global index is (2,2,2), with origin (2,1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 1, 3)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 1, 3));
-            break;
-        case 4:
-            // global index is (2,2,2), with origin (-1,-1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(3, 3, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(3, 3, 0));
-            break;
-        case 5:
-            // global index is (2,2,2), with origin (2,-1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 3, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 3, 0));
-            break;
-        case 6:
-            // global index is (2,2,2), with origin (-1,1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(3, 1, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(3, 1, 0));
-            break;
-        case 7:
-            // global index is (2,2,2), with origin (2,1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 1, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 1, 0));
-            break;
-            };
+        const unsigned int local_cell = make_local_cell(cl, 2, 2, 2);
+        UP_ASSERT_EQUAL(h_cell_np.data[local_cell], 1);
+        UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), local_cell);
         }
 
     // apply a grid shift, particles on left internal boundary will move up one cell
     // particles on right internal boundary will stay in the same cell
-    cl->setGridShift(make_scalar3(-0.5, -0.5, -0.5));
-    cl->compute(1);
+    const Scalar3 shift = (Scalar(0.5) / 5) * make_scalar3(1, 1, 1);
         {
+        cl->setGridShift(-shift);
+        cl->compute(1);
+        checkDomainBoundaries(sysdef, cl);
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
                                             access_location::host,
                                             access_mode::read);
-        ArrayHandle<unsigned int> h_cell_list(cl->getCellList(),
-                                              access_location::host,
-                                              access_mode::read);
-        Index3D ci = cl->getCellIndexer();
         ArrayHandle<Scalar4> h_vel(pdata->getVelocities(),
                                    access_location::host,
                                    access_mode::read);
 
+        int3 cell;
         switch (my_rank)
             {
         case 0:
-            // global index is (2,2,2), with origin (-1,-1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(3, 3, 3)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(3, 3, 3));
+            cell = make_int3(2, 2, 2);
             break;
         case 1:
-            // global index is (3,2,2), with origin (2,-1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(1, 3, 3)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(1, 3, 3));
+            cell = make_int3(3, 2, 2);
             break;
         case 2:
-            // global index is (2,3,2), with origin (-1,1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(3, 2, 3)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(3, 2, 3));
+            cell = make_int3(2, 3, 2);
             break;
         case 3:
-            // global index is (3,3,2), with origin (2,1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(1, 2, 3)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(1, 2, 3));
+            cell = make_int3(3, 3, 2);
             break;
         case 4:
-            // global index is (2,2,3), with origin (-1,-1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(3, 3, 1)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(3, 3, 1));
+            cell = make_int3(2, 2, 3);
             break;
         case 5:
-            // global index is (3,2,3), with origin (2,-1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(1, 3, 1)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(1, 3, 1));
+            cell = make_int3(3, 2, 3);
             break;
         case 6:
-            // global index is (2,3,3), with origin (-1,1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(3, 2, 1)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(3, 2, 1));
+            cell = make_int3(2, 3, 3);
             break;
         case 7:
-            // global index is (3,3,3), with origin (2,1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(1, 2, 1)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(1, 2, 1));
+            cell = make_int3(3, 3, 3);
             break;
             };
+
+        const unsigned int local_cell = make_local_cell(cl, cell.x, cell.y, cell.z);
+        UP_ASSERT_EQUAL(h_cell_np.data[local_cell], 1);
+        UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), local_cell);
         }
 
-    // apply a grid shift, particles on left internal boundary will stay in cell
-    // particles on right internal boundary will move down one cell
-    cl->setGridShift(make_scalar3(0.5, 0.5, 0.5));
-    cl->compute(2);
+        // apply a grid shift, particles on left internal boundary will stay in cell
+        // particles on right internal boundary will move down one cell
+        // if (*box == *ref_box)
         {
+        cl->setGridShift(shift);
+        cl->compute(2);
+        checkDomainBoundaries(sysdef, cl);
+
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
                                             access_location::host,
                                             access_mode::read);
-        ArrayHandle<unsigned int> h_cell_list(cl->getCellList(),
-                                              access_location::host,
-                                              access_mode::read);
-        Index3D ci = cl->getCellIndexer();
         ArrayHandle<Scalar4> h_vel(pdata->getVelocities(),
                                    access_location::host,
                                    access_mode::read);
 
+        int3 cell;
         switch (my_rank)
             {
         case 0:
-            // global index is (1,1,1), with origin (-1,-1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(2, 2, 2)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(2, 2, 2));
+            cell = make_int3(1, 1, 1);
             break;
         case 1:
-            // global index is (2,1,1), with origin (2,-1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 2, 2)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 2, 2));
+            cell = make_int3(2, 1, 1);
             break;
         case 2:
-            // global index is (1,2,1), with origin (-1,1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(2, 1, 2)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(2, 1, 2));
+            cell = make_int3(1, 2, 1);
             break;
         case 3:
-            // global index is (2,2,1), with origin (2,1,-1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 1, 2)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 1, 2));
+            cell = make_int3(2, 2, 1);
             break;
         case 4:
-            // global index is (1,1,2), with origin (-1,-1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(2, 2, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(2, 2, 0));
+            cell = make_int3(1, 1, 2);
             break;
         case 5:
-            // global index is (2,1,2), with origin (2,-1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 2, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 2, 0));
+            cell = make_int3(2, 1, 2);
             break;
         case 6:
-            // global index is (1,2,2), with origin (-1,1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(2, 1, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(2, 1, 0));
+            cell = make_int3(1, 2, 2);
             break;
         case 7:
-            // global index is (2,2,2), with origin (2,1,2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 1, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 1, 0));
+            cell = make_int3(2, 2, 2);
             break;
             };
+
+        const unsigned int local_cell = make_local_cell(cl, cell.x, cell.y, cell.z);
+        UP_ASSERT_EQUAL(h_cell_np.data[local_cell], 1);
+        UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), local_cell);
         }
+
     // now we move the particles to their exterior boundaries, and repeat the testing process
     // we are going to pad the cell list with an extra cell just to test that binning now
     cl->setNExtraCells(1);
@@ -1213,28 +1204,28 @@ template<class CL> void celllist_edge_test(std::shared_ptr<ExecutionConfiguratio
         switch (my_rank)
             {
         case 0:
-            h_pos.data[0] = make_scalar4(-4.0, -4.0, -4.0, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(-4.0, -4.0, -4.0, __int_as_scalar(0)), ref_box, box);
             break;
         case 1:
-            h_pos.data[0] = make_scalar4(3.99, -4.0, -4.0, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(3.99, -4.0, -4.0, __int_as_scalar(0)), ref_box, box);
             break;
         case 2:
-            h_pos.data[0] = make_scalar4(-4.0, 3.99, -4.0, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(-4.0, 3.99, -4.0, __int_as_scalar(0)), ref_box, box);
             break;
         case 3:
-            h_pos.data[0] = make_scalar4(3.99, 3.99, -4.0, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(3.99, 3.99, -4.0, __int_as_scalar(0)), ref_box, box);
             break;
         case 4:
-            h_pos.data[0] = make_scalar4(-4.0, -4.0, 3.99, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(-4.0, -4.0, 3.99, __int_as_scalar(0)), ref_box, box);
             break;
         case 5:
-            h_pos.data[0] = make_scalar4(3.99, -4.0, 3.99, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(3.99, -4.0, 3.99, __int_as_scalar(0)), ref_box, box);
             break;
         case 6:
-            h_pos.data[0] = make_scalar4(-4.0, 3.99, 3.99, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(-4.0, 3.99, 3.99, __int_as_scalar(0)), ref_box, box);
             break;
         case 7:
-            h_pos.data[0] = make_scalar4(3.99, 3.99, 3.99, __int_as_scalar(0));
+            h_pos.data[0] = scale(make_scalar4(3.99, 3.99, 3.99, __int_as_scalar(0)), ref_box, box);
             break;
             };
         }
@@ -1242,184 +1233,146 @@ template<class CL> void celllist_edge_test(std::shared_ptr<ExecutionConfiguratio
     // reset the grid shift and recompute
     cl->setGridShift(make_scalar3(0, 0, 0));
     cl->compute(3);
+    checkDomainBoundaries(sysdef, cl);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
                                             access_location::host,
                                             access_mode::read);
-        ArrayHandle<unsigned int> h_cell_list(cl->getCellList(),
-                                              access_location::host,
-                                              access_mode::read);
-        Index3D ci = cl->getCellIndexer();
         ArrayHandle<Scalar4> h_vel(pdata->getVelocities(),
                                    access_location::host,
                                    access_mode::read);
 
+        int3 cell;
         switch (my_rank)
             {
         case 0:
-            // global index is (-2,-2,-2), with origin (-2,-2,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 0, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 0, 0));
+            cell = make_int3(-2, -2, -2);
             break;
         case 1:
-            // global index is (6,-2,-2), with origin (1,-2,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(5, 0, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(5, 0, 0));
+            cell = make_int3(6, -2, -2);
             break;
         case 2:
-            // global index is (-2,6,-2), with origin (-2,0,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 6, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 6, 0));
+            cell = make_int3(-2, 6, -2);
             break;
         case 3:
-            // global index is (6,6,-2), with origin (1,0,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(5, 6, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(5, 6, 0));
+            cell = make_int3(6, 6, -2);
             break;
         case 4:
-            // global index is (-2,-2,6), with origin (-2,-2,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 0, 5)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 0, 5));
+            cell = make_int3(-2, -2, 6);
             break;
         case 5:
-            // global index is (6,-2,6), with origin (1,-2,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(5, 0, 5)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(5, 0, 5));
+            cell = make_int3(6, -2, 6);
             break;
         case 6:
-            // global index is (-2,6,6), with origin (-2,0,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 6, 5)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 6, 5));
+            cell = make_int3(-2, 6, 6);
             break;
         case 7:
-            // global index is (6,6,6), with origin (1,0,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(5, 6, 5)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(5, 6, 5));
+            cell = make_int3(6, 6, 6);
             break;
             };
+
+        const unsigned int local_cell = make_local_cell(cl, cell.x, cell.y, cell.z);
+        UP_ASSERT_EQUAL(h_cell_np.data[local_cell], 1);
+        UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), local_cell);
         }
 
     // shift all particles to the right by 0.5
     // all particles on left bound will move up one cell, all particles on right bound stay
-    cl->setGridShift(make_scalar3(-0.5, -0.5, -0.5));
-    cl->compute(4);
+    // this test has weird round off issues for triclinic, so only run for orthorhombic boxes
+    if (is_orthorhombic)
         {
+        cl->setGridShift(-shift);
+        cl->compute(4);
+        checkDomainBoundaries(sysdef, cl);
+
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
                                             access_location::host,
                                             access_mode::read);
-        ArrayHandle<unsigned int> h_cell_list(cl->getCellList(),
-                                              access_location::host,
-                                              access_mode::read);
-        Index3D ci = cl->getCellIndexer();
         ArrayHandle<Scalar4> h_vel(pdata->getVelocities(),
                                    access_location::host,
                                    access_mode::read);
 
+        int3 cell;
         switch (my_rank)
             {
         case 0:
-            // global index is (-1,-1,-1), with origin (-2,-2,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(1, 1, 1)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(1, 1, 1));
+            cell = make_int3(-2, -2, -2);
             break;
         case 1:
-            // global index is (6,-1,-1), with origin (1,-2,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(5, 1, 1)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(5, 1, 1));
+            cell = make_int3(6, -2, -2);
             break;
         case 2:
-            // global index is (-1,6,-1), with origin (-2,0,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(1, 6, 1)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(1, 6, 1));
+            cell = make_int3(-2, 6, -2);
             break;
         case 3:
-            // global index is (6,6,-1), with origin (1,0,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(5, 6, 1)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(5, 6, 1));
+            cell = make_int3(6, 6, -2);
             break;
         case 4:
-            // global index is (-1,-1,6), with origin (-2,-2,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(1, 1, 5)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(1, 1, 5));
+            cell = make_int3(-2, -2, 6);
             break;
         case 5:
-            // global index is (6,-1,6), with origin (1,-2,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(5, 1, 5)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(5, 1, 5));
+            cell = make_int3(6, -2, 6);
             break;
         case 6:
-            // global index is (-1,6,6), with origin (-2,0,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(1, 6, 5)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(1, 6, 5));
+            cell = make_int3(-2, 6, 6);
             break;
         case 7:
-            // global index is (6,6,6), with origin (1,0,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(5, 6, 5)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(5, 6, 5));
+            cell = make_int3(6, 6, 6);
             break;
             };
+
+        const unsigned int local_cell = make_local_cell(cl, cell.x, cell.y, cell.z);
+        UP_ASSERT_EQUAL(h_cell_np.data[local_cell], 1);
+        UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), local_cell);
         }
 
-    // shift all particles left by 0.5
-    // particles on the lower bound stay in the same bin, and particles on the right bound move down
-    // one
-    cl->setGridShift(make_scalar3(0.5, 0.5, 0.5));
-    cl->compute(5);
+        // shift all particles left by 0.5
+        // particles on the lower bound stay in the same bin, and particles on the right bound move
+        // down one
         {
+        cl->setGridShift(shift);
+        cl->compute(5);
+        checkDomainBoundaries(sysdef, cl);
+
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
                                             access_location::host,
                                             access_mode::read);
-        ArrayHandle<unsigned int> h_cell_list(cl->getCellList(),
-                                              access_location::host,
-                                              access_mode::read);
-        Index3D ci = cl->getCellIndexer();
         ArrayHandle<Scalar4> h_vel(pdata->getVelocities(),
                                    access_location::host,
                                    access_mode::read);
 
+        int3 cell;
         switch (my_rank)
             {
         case 0:
-            // global index is (-2,-2,-2), with origin (-2,-2,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 0, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 0, 0));
+            cell = make_int3(-2, -2, -2);
             break;
         case 1:
-            // global index is (5,-2,-2), with origin (1,-2,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(4, 0, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(4, 0, 0));
+            cell = make_int3(5, -2, -2);
             break;
         case 2:
-            // global index is (-2,5,-2), with origin (-2,0,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 5, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 5, 0));
+            cell = make_int3(-2, 5, -2);
             break;
         case 3:
-            // global index is (5,5,-2), with origin (1,0,-2)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(4, 5, 0)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(4, 5, 0));
+            cell = make_int3(5, 5, -2);
             break;
         case 4:
-            // global index is (-2,-2,5), with origin (-2,-2,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 0, 4)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 0, 4));
+            cell = make_int3(-2, -2, 5);
             break;
         case 5:
-            // global index is (5,-2,5), with origin (1,-2,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(4, 0, 4)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(4, 0, 4));
+            cell = make_int3(5, -2, 5);
             break;
         case 6:
-            // global index is (-2,5,5), with origin (-2,0,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(0, 5, 4)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(0, 5, 4));
+            cell = make_int3(-2, 5, 5);
             break;
         case 7:
-            // global index is (5,5,5), with origin (1,0,1)
-            UP_ASSERT_EQUAL(h_cell_np.data[ci(4, 5, 4)], 1);
-            UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), ci(4, 5, 4));
+            cell = make_int3(5, 5, 5);
             break;
             };
+
+        const unsigned int local_cell = make_local_cell(cl, cell.x, cell.y, cell.z);
+        UP_ASSERT_EQUAL(h_cell_np.data[local_cell], 1);
+        UP_ASSERT_EQUAL(__scalar_as_int(h_vel.data[0].w), local_cell);
         }
     }
 
@@ -1431,40 +1384,159 @@ UP_TEST(mpcd_cell_list_dimensions)
         std::shared_ptr<ExecutionConfiguration> exec_conf(
             new ExecutionConfiguration(ExecutionConfiguration::CPU));
         exec_conf->getMPIConfig()->splitPartitions(2);
-        celllist_dimension_test<mpcd::CellList>(exec_conf, true, false, false);
-        celllist_dimension_test<mpcd::CellList>(exec_conf, false, true, false);
-        celllist_dimension_test<mpcd::CellList>(exec_conf, false, false, true);
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                true,
+                                                false,
+                                                false,
+                                                make_scalar3(0, 0, 0));
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                false,
+                                                true,
+                                                false,
+                                                make_scalar3(0, 0, 0));
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                false,
+                                                false,
+                                                true,
+                                                make_scalar3(0, 0, 0));
         }
         // mpi in 2d
         {
         std::shared_ptr<ExecutionConfiguration> exec_conf(
             new ExecutionConfiguration(ExecutionConfiguration::CPU));
         exec_conf->getMPIConfig()->splitPartitions(4);
-        celllist_dimension_test<mpcd::CellList>(exec_conf, true, true, false);
-        celllist_dimension_test<mpcd::CellList>(exec_conf, true, false, true);
-        celllist_dimension_test<mpcd::CellList>(exec_conf, false, true, true);
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                true,
+                                                true,
+                                                false,
+                                                make_scalar3(0, 0, 0));
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                true,
+                                                false,
+                                                true,
+                                                make_scalar3(0, 0, 0));
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                false,
+                                                true,
+                                                true,
+                                                make_scalar3(0, 0, 0));
         }
         // mpi in 3d
         {
         std::shared_ptr<ExecutionConfiguration> exec_conf(
             new ExecutionConfiguration(ExecutionConfiguration::CPU));
         exec_conf->getMPIConfig()->splitPartitions(8);
-        celllist_dimension_test<mpcd::CellList>(exec_conf, true, true, true);
+        celllist_dimension_test<mpcd::CellList>(exec_conf, true, true, true, make_scalar3(0, 0, 0));
+        }
+    }
+
+//! dimension test case for MPCD CellList class, triclinic
+UP_TEST(mpcd_cell_list_dimensions_triclinic)
+    {
+        // mpi in 1d
+        {
+        auto exec_conf = std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU);
+        exec_conf->getMPIConfig()->splitPartitions(2);
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                true,
+                                                false,
+                                                false,
+                                                make_scalar3(0.5, -0.75, 1.0));
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                false,
+                                                true,
+                                                false,
+                                                make_scalar3(0.5, -0.75, 1.0));
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                false,
+                                                false,
+                                                true,
+                                                make_scalar3(0.5, -0.75, 1.0));
+        }
+        // mpi in 2d
+        {
+        auto exec_conf = std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU);
+        exec_conf->getMPIConfig()->splitPartitions(4);
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                true,
+                                                true,
+                                                false,
+                                                make_scalar3(0.5, -0.75, 1.0));
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                true,
+                                                false,
+                                                true,
+                                                make_scalar3(0.5, -0.75, 1.0));
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                false,
+                                                true,
+                                                true,
+                                                make_scalar3(0.5, -0.75, 1.0));
+        }
+        // mpi in 3d
+        {
+        auto exec_conf = std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU);
+        exec_conf->getMPIConfig()->splitPartitions(8);
+        celllist_dimension_test<mpcd::CellList>(exec_conf,
+                                                true,
+                                                true,
+                                                true,
+                                                make_scalar3(0.5, -0.75, 1.0));
         }
     }
 
 //! basic test case for MPCD CellList class
 UP_TEST(mpcd_cell_list_basic_test)
     {
-    celllist_basic_test<mpcd::CellList>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::CPU)));
+    celllist_basic_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.0, 6.0, 6.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! basic test case for MPCD CellList class, noncubic
+UP_TEST(mpcd_cell_list_basic_test_noncubic)
+    {
+    celllist_basic_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.5, 7.0, 7.5),
+        make_scalar3(0, 0, 0));
+    }
+
+//! basic test case for MPCD CellList class, triclinic
+UP_TEST(mpcd_cell_list_basic_test_triclinic)
+    {
+    celllist_basic_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.0, 6.0, 6.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 
 //! edge test case for MPCD CellList class
 UP_TEST(mpcd_cell_list_edge_test)
     {
-    celllist_edge_test<mpcd::CellList>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::CPU)));
+    celllist_edge_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(5.0, 5.0, 5.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! edge test case for MPCD CellList class, noncubic
+UP_TEST(mpcd_cell_list_edge_test_noncubic)
+    {
+    celllist_edge_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.0, 6.5, 7.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! edge test case for MPCD CellList class, triclinic
+UP_TEST(mpcd_cell_list_edge_test_triclinic)
+    {
+    celllist_edge_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(5.0, 5.0, 5.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 
 #ifdef ENABLE_HIP
@@ -1476,39 +1548,162 @@ UP_TEST(mpcd_cell_list_gpu_dimensions)
         std::shared_ptr<ExecutionConfiguration> exec_conf(
             new ExecutionConfiguration(ExecutionConfiguration::GPU));
         exec_conf->getMPIConfig()->splitPartitions(2);
-        celllist_dimension_test<mpcd::CellListGPU>(exec_conf, true, false, false);
-        celllist_dimension_test<mpcd::CellListGPU>(exec_conf, false, true, false);
-        celllist_dimension_test<mpcd::CellListGPU>(exec_conf, false, false, true);
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   true,
+                                                   false,
+                                                   false,
+                                                   make_scalar3(0, 0, 0));
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   false,
+                                                   true,
+                                                   false,
+                                                   make_scalar3(0, 0, 0));
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   false,
+                                                   false,
+                                                   true,
+                                                   make_scalar3(0, 0, 0));
         }
         // mpi in 2d
         {
         std::shared_ptr<ExecutionConfiguration> exec_conf(
             new ExecutionConfiguration(ExecutionConfiguration::GPU));
         exec_conf->getMPIConfig()->splitPartitions(4);
-        celllist_dimension_test<mpcd::CellListGPU>(exec_conf, true, true, false);
-        celllist_dimension_test<mpcd::CellListGPU>(exec_conf, true, false, true);
-        celllist_dimension_test<mpcd::CellListGPU>(exec_conf, false, true, true);
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   true,
+                                                   true,
+                                                   false,
+                                                   make_scalar3(0, 0, 0));
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   true,
+                                                   false,
+                                                   true,
+                                                   make_scalar3(0, 0, 0));
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   false,
+                                                   true,
+                                                   true,
+                                                   make_scalar3(0, 0, 0));
         }
         // mpi in 3d
         {
         std::shared_ptr<ExecutionConfiguration> exec_conf(
             new ExecutionConfiguration(ExecutionConfiguration::GPU));
         exec_conf->getMPIConfig()->splitPartitions(8);
-        celllist_dimension_test<mpcd::CellListGPU>(exec_conf, true, true, true);
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   true,
+                                                   true,
+                                                   true,
+                                                   make_scalar3(0, 0, 0));
+        }
+    }
+
+//! dimension test case for MPCD CellListGPU class
+UP_TEST(mpcd_cell_list_gpu_dimensions_triclinic)
+    {
+        // mpi in 1d
+        {
+        auto exec_conf = std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU);
+        exec_conf->getMPIConfig()->splitPartitions(2);
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   true,
+                                                   false,
+                                                   false,
+                                                   make_scalar3(0.5, -0.75, 1.0));
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   false,
+                                                   true,
+                                                   false,
+                                                   make_scalar3(0.5, -0.75, 1.0));
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   false,
+                                                   false,
+                                                   true,
+                                                   make_scalar3(0.5, -0.75, 1.0));
+        }
+        // mpi in 2d
+        {
+        auto exec_conf = std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU);
+        exec_conf->getMPIConfig()->splitPartitions(4);
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   true,
+                                                   true,
+                                                   false,
+                                                   make_scalar3(0.5, -0.75, 1.0));
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   true,
+                                                   false,
+                                                   true,
+                                                   make_scalar3(0.5, -0.75, 1.0));
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   false,
+                                                   true,
+                                                   true,
+                                                   make_scalar3(0.5, -0.75, 1.0));
+        }
+        // mpi in 3d
+        {
+        auto exec_conf = std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU);
+        exec_conf->getMPIConfig()->splitPartitions(8);
+        celllist_dimension_test<mpcd::CellListGPU>(exec_conf,
+                                                   true,
+                                                   true,
+                                                   true,
+                                                   make_scalar3(0.5, -0.75, 1.0));
         }
     }
 
 //! basic test case for MPCD CellListGPU class
 UP_TEST(mpcd_cell_list_gpu_basic_test)
     {
-    celllist_basic_test<mpcd::CellListGPU>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::GPU)));
+    celllist_basic_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.0, 6.0, 6.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! basic test case for MPCD CellListGPU class, noncubic
+UP_TEST(mpcd_cell_list_gpu_basic_test_noncubic)
+    {
+    celllist_basic_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.5, 7.0, 7.5),
+        make_scalar3(0, 0, 0));
+    }
+
+//! basic test case for MPCD CellListGPU class, triclinic
+UP_TEST(mpcd_cell_list_gpu_basic_test_triclinic)
+    {
+    celllist_basic_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.0, 6.0, 6.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 
 //! edge test case for MPCD CellListGPU class
 UP_TEST(mpcd_cell_list_gpu_edge_test)
     {
-    celllist_edge_test<mpcd::CellListGPU>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::GPU)));
+    celllist_edge_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(5.0, 5.0, 5.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! edge test case for MPCD CellListGPU class, noncubic
+UP_TEST(mpcd_cell_list_gpu_edge_test_noncubic)
+    {
+    celllist_edge_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.0, 6.5, 7.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! edge test case for MPCD CellListGPU class, triclinic
+UP_TEST(mpcd_cell_list_gpu_edge_test_triclinic)
+    {
+    celllist_edge_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(5.0, 5.0, 5.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 #endif // ENABLE_HIP
diff --git a/hoomd/mpcd/test/cell_list_test.cc b/hoomd/mpcd/test/cell_list_test.cc
index b0c3e22017..0d96bdd458 100644
--- a/hoomd/mpcd/test/cell_list_test.cc
+++ b/hoomd/mpcd/test/cell_list_test.cc
@@ -11,15 +11,21 @@
 #include "hoomd/filter/ParticleFilterType.h"
 #include "hoomd/test/upp11_config.h"
 
+#include "utils.h"
+
 HOOMD_UP_MAIN()
 
 using namespace hoomd;
 
 //! Test for correct calculation of MPCD grid dimensions
-template<class CL> void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf)
+template<class CL>
+void celllist_dimension_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
+                             const Scalar3& L,
+                             const Scalar3& tilt)
     {
     std::shared_ptr<SnapshotSystemData<Scalar>> snap(new SnapshotSystemData<Scalar>());
-    snap->global_box = std::make_shared<BoxDim>(6.0, 8.0, 10.0);
+    snap->global_box = std::make_shared<BoxDim>(L);
+    snap->global_box->setTiltFactors(tilt.x, tilt.y, tilt.z);
     snap->particle_data.type_mapping.push_back("A");
     snap->mpcd_data.resize(1);
     snap->mpcd_data.type_mapping.push_back("A");
@@ -28,7 +34,7 @@ template<class CL> void celllist_dimension_test(std::shared_ptr<ExecutionConfigu
     auto pdata_1 = sysdef->getMPCDParticleData();
 
     // define a system of different edge lengths
-    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, 1.0, false));
+    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, make_uint3(6, 8, 10), false));
 
     // compute the cell list dimensions
     cl->computeDimensions();
@@ -57,7 +63,7 @@ template<class CL> void celllist_dimension_test(std::shared_ptr<ExecutionConfigu
 
     /*******************/
     // Change the cell size, and ensure everything stays up to date
-    cl->setCellSize(2.0);
+    cl->setGlobalDim(make_uint3(3, 4, 5));
     cl->computeDimensions();
 
     dim = cl->getDim();
@@ -76,46 +82,38 @@ template<class CL> void celllist_dimension_test(std::shared_ptr<ExecutionConfigu
 
     // Cell list uses amortized sizing, so must only be at least this big
     UP_ASSERT(cl->getCellList().getNumElements() >= 3 * 4 * 5 * Nmax);
-
-    /*******************/
-    // Change the cell size to something that does not evenly divide a side, and check for an
-    // exception evenly divides no sides
-    cl->setCellSize(4.2);
-    UP_ASSERT_EXCEPTION(std::runtime_error, [&] { cl->computeDimensions(); });
-    // evenly divides z
-    cl->setCellSize(10.0);
-    UP_ASSERT_EXCEPTION(std::runtime_error, [&] { cl->computeDimensions(); });
-    // evenly divides y
-    cl->setCellSize(8.0);
-    UP_ASSERT_EXCEPTION(std::runtime_error, [&] { cl->computeDimensions(); });
-    // evenly divides x
-    cl->setCellSize(6.0);
-    UP_ASSERT_EXCEPTION(std::runtime_error, [&] { cl->computeDimensions(); });
     }
 
 //! Test for correct cell listing of a small system
-template<class CL> void celllist_small_test(std::shared_ptr<ExecutionConfiguration> exec_conf)
+template<class CL>
+void celllist_small_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
+                         const Scalar3& L,
+                         const Scalar3& tilt)
     {
+    auto ref_box = std::make_shared<BoxDim>(2.0);
+    auto box = std::make_shared<BoxDim>(L);
+    box->setTiltFactors(tilt.x, tilt.y, tilt.z);
+
     std::shared_ptr<SnapshotSystemData<Scalar>> snap(new SnapshotSystemData<Scalar>());
-    snap->global_box = std::make_shared<BoxDim>(2.0);
+    snap->global_box = box;
     snap->particle_data.type_mapping.push_back("A");
     // place each particle in a different cell, doubling the first cell
     snap->mpcd_data.resize(9);
     snap->mpcd_data.type_mapping.push_back("A");
-    snap->mpcd_data.position[0] = vec3<Scalar>(-0.5, -0.5, -0.5);
-    snap->mpcd_data.position[1] = vec3<Scalar>(0.5, -0.5, -0.5);
-    snap->mpcd_data.position[2] = vec3<Scalar>(-0.5, 0.5, -0.5);
-    snap->mpcd_data.position[3] = vec3<Scalar>(0.5, 0.5, -0.5);
-    snap->mpcd_data.position[4] = vec3<Scalar>(-0.5, -0.5, 0.5);
-    snap->mpcd_data.position[5] = vec3<Scalar>(0.5, -0.5, 0.5);
-    snap->mpcd_data.position[6] = vec3<Scalar>(-0.5, 0.5, 0.5);
-    snap->mpcd_data.position[7] = vec3<Scalar>(0.5, 0.5, 0.5);
-    snap->mpcd_data.position[8] = vec3<Scalar>(-0.5, -0.5, -0.5);
+    snap->mpcd_data.position[0] = scale(vec3<Scalar>(-0.5, -0.5, -0.5), ref_box, box);
+    snap->mpcd_data.position[1] = scale(vec3<Scalar>(0.5, -0.5, -0.5), ref_box, box);
+    snap->mpcd_data.position[2] = scale(vec3<Scalar>(-0.5, 0.5, -0.5), ref_box, box);
+    snap->mpcd_data.position[3] = scale(vec3<Scalar>(0.5, 0.5, -0.5), ref_box, box);
+    snap->mpcd_data.position[4] = scale(vec3<Scalar>(-0.5, -0.5, 0.5), ref_box, box);
+    snap->mpcd_data.position[5] = scale(vec3<Scalar>(0.5, -0.5, 0.5), ref_box, box);
+    snap->mpcd_data.position[6] = scale(vec3<Scalar>(-0.5, 0.5, 0.5), ref_box, box);
+    snap->mpcd_data.position[7] = scale(vec3<Scalar>(0.5, 0.5, 0.5), ref_box, box);
+    snap->mpcd_data.position[8] = scale(vec3<Scalar>(-0.5, -0.5, -0.5), ref_box, box);
     std::shared_ptr<SystemDefinition> sysdef(new SystemDefinition(snap, exec_conf));
 
     std::shared_ptr<mpcd::ParticleData> pdata_9 = sysdef->getMPCDParticleData();
 
-    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, 1.0, false));
+    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, make_uint3(2, 2, 2), false));
     cl->compute(0);
 
         // check that each particle is in the proper bin (cell list and velocity)
@@ -170,8 +168,8 @@ template<class CL> void celllist_small_test(std::shared_ptr<ExecutionConfigurati
         ArrayHandle<Scalar4> h_pos(pdata_9->getPositions(),
                                    access_location::host,
                                    access_mode::overwrite);
-        h_pos.data[0] = make_scalar4(-0.3, -0.3, -0.3, 0.0);
-        h_pos.data[1] = make_scalar4(0.3, 0.3, 0.3, 0.0);
+        h_pos.data[0] = scale(make_scalar4(-0.3, -0.3, -0.3, __int_as_scalar(0)), ref_box, box);
+        h_pos.data[1] = scale(make_scalar4(0.3, 0.3, 0.3, __int_as_scalar(0)), ref_box, box);
         h_pos.data[2] = h_pos.data[0];
         h_pos.data[3] = h_pos.data[1];
         h_pos.data[4] = h_pos.data[0];
@@ -233,7 +231,7 @@ template<class CL> void celllist_small_test(std::shared_ptr<ExecutionConfigurati
         ArrayHandle<Scalar4> h_pos(pdata_9->getPositions(),
                                    access_location::host,
                                    access_mode::overwrite);
-        h_pos.data[0] = make_scalar4(0.9, -0.4, 0.0, 0.0);
+        h_pos.data[0] = scale(make_scalar4(0.9, -0.4, 0.0, __int_as_scalar(0)), ref_box, box);
         for (unsigned int i = 1; i < 9; ++i)
             h_pos.data[i] = h_pos.data[0];
         }
@@ -260,7 +258,7 @@ template<class CL> void celllist_small_test(std::shared_ptr<ExecutionConfigurati
         ArrayHandle<Scalar4> h_pos(pdata_9->getPositions(),
                                    access_location::host,
                                    access_mode::overwrite);
-        h_pos.data[0] = make_scalar4(2.1, 2.1, 2.1, __int_as_scalar(0));
+        h_pos.data[0] = scale(make_scalar4(2.1, 2.1, 2.1, __int_as_scalar(0)), ref_box, box);
         }
     UP_ASSERT_EXCEPTION(std::runtime_error, [&] { cl->compute(3); });
         // check the other side as well
@@ -268,24 +266,31 @@ template<class CL> void celllist_small_test(std::shared_ptr<ExecutionConfigurati
         ArrayHandle<Scalar4> h_pos(pdata_9->getPositions(),
                                    access_location::host,
                                    access_mode::overwrite);
-        h_pos.data[0] = make_scalar4(-2.1, -2.1, -2.1, __int_as_scalar(0));
+        h_pos.data[0] = scale(make_scalar4(-2.1, -2.1, -2.1, __int_as_scalar(0)), ref_box, box);
         }
     UP_ASSERT_EXCEPTION(std::runtime_error, [&] { cl->compute(4); });
     }
 
 //! Test that particles can be grid shifted correctly
-template<class CL> void celllist_grid_shift_test(std::shared_ptr<ExecutionConfiguration> exec_conf)
+template<class CL>
+void celllist_grid_shift_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
+                              const Scalar3& L,
+                              const Scalar3& tilt)
     {
+    auto ref_box = std::make_shared<BoxDim>(6.0);
+    auto box = std::make_shared<BoxDim>(L);
+    box->setTiltFactors(tilt.x, tilt.y, tilt.z);
+
     std::shared_ptr<SnapshotSystemData<Scalar>> snap(new SnapshotSystemData<Scalar>());
-    snap->global_box = std::make_shared<BoxDim>(6.0);
+    snap->global_box = box;
     snap->particle_data.type_mapping.push_back("A");
     snap->mpcd_data.resize(1);
     snap->mpcd_data.type_mapping.push_back("A");
-    snap->mpcd_data.position[0] = vec3<Scalar>(0.1, 0.1, 0.1);
+    snap->mpcd_data.position[0] = scale(vec3<Scalar>(0.1, 0.1, 0.1), ref_box, box);
     std::shared_ptr<SystemDefinition> sysdef(new SystemDefinition(snap, exec_conf));
 
     std::shared_ptr<mpcd::ParticleData> pdata_1 = sysdef->getMPCDParticleData();
-    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, 1.0, false));
+    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, make_uint3(6, 6, 6), false));
     cl->compute(0);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
@@ -296,7 +301,8 @@ template<class CL> void celllist_grid_shift_test(std::shared_ptr<ExecutionConfig
         }
 
     // shift the grid and see that it falls from (3,3,3) to (2,2,2)
-    cl->setGridShift(make_scalar3(0.5, 0.5, 0.5));
+    const Scalar3 shift = (Scalar(0.5) / 6) * make_scalar3(1, 1, 1);
+    cl->setGridShift(shift);
     cl->compute(1);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
@@ -311,9 +317,9 @@ template<class CL> void celllist_grid_shift_test(std::shared_ptr<ExecutionConfig
         ArrayHandle<Scalar4> h_pos(pdata_1->getPositions(),
                                    access_location::host,
                                    access_mode::overwrite);
-        h_pos.data[0] = make_scalar4(-0.1, -0.1, -0.1, 0.0);
+        h_pos.data[0] = scale(make_scalar4(-0.1, -0.1, -0.1, __int_as_scalar(0)), ref_box, box);
         }
-    cl->setGridShift(make_scalar3(-0.5, -0.5, -0.5));
+    cl->setGridShift(-shift);
     cl->compute(2);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
@@ -328,9 +334,9 @@ template<class CL> void celllist_grid_shift_test(std::shared_ptr<ExecutionConfig
         ArrayHandle<Scalar4> h_pos(pdata_1->getPositions(),
                                    access_location::host,
                                    access_mode::overwrite);
-        h_pos.data[0] = make_scalar4(-2.9, -2.9, -2.9, 0.0);
+        h_pos.data[0] = scale(make_scalar4(-2.9, -2.9, -2.9, __int_as_scalar(0)), ref_box, box);
         }
-    cl->setGridShift(make_scalar3(0.5, 0.5, 0.5));
+    cl->setGridShift(shift);
     cl->compute(3);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
@@ -345,9 +351,9 @@ template<class CL> void celllist_grid_shift_test(std::shared_ptr<ExecutionConfig
         ArrayHandle<Scalar4> h_pos(pdata_1->getPositions(),
                                    access_location::host,
                                    access_mode::overwrite);
-        h_pos.data[0] = make_scalar4(2.9, 2.9, 2.9, 0.0);
+        h_pos.data[0] = scale(make_scalar4(2.9, 2.9, 2.9, __int_as_scalar(0)), ref_box, box);
         }
-    cl->setGridShift(make_scalar3(-0.5, -0.5, -0.5));
+    cl->setGridShift(-shift);
     cl->compute(4);
         {
         ArrayHandle<unsigned int> h_cell_np(cl->getCellSizeArray(),
@@ -358,31 +364,36 @@ template<class CL> void celllist_grid_shift_test(std::shared_ptr<ExecutionConfig
         }
 
     // check for error in grid shifting
-    UP_ASSERT_EXCEPTION(std::runtime_error,
-                        [&] { cl->setGridShift(make_scalar3(-0.51, -0.51, -0.51)); });
-    UP_ASSERT_EXCEPTION(std::runtime_error,
-                        [&] { cl->setGridShift(make_scalar3(0.51, 0.51, 0.51)); });
+    UP_ASSERT_EXCEPTION(std::runtime_error, [&] { cl->setGridShift(Scalar(-1.01) * shift); });
+    UP_ASSERT_EXCEPTION(std::runtime_error, [&] { cl->setGridShift(Scalar(1.01) * shift); });
     }
 
 //! Test that small systems can embed particles
-template<class CL> void celllist_embed_test(std::shared_ptr<ExecutionConfiguration> exec_conf)
+template<class CL>
+void celllist_embed_test(std::shared_ptr<ExecutionConfiguration> exec_conf,
+                         const Scalar3& L,
+                         const Scalar3& tilt)
     {
+    auto ref_box = std::make_shared<BoxDim>(2.0);
+    auto box = std::make_shared<BoxDim>(L);
+    box->setTiltFactors(tilt.x, tilt.y, tilt.z);
+
     // setup a system where both MD and MPCD particles are in each of the cells
     std::shared_ptr<SnapshotSystemData<Scalar>> snap(new SnapshotSystemData<Scalar>());
-    snap->global_box = std::make_shared<BoxDim>(2.0);
+    snap->global_box = box;
         {
         SnapshotParticleData<Scalar>& pdata_snap = snap->particle_data;
         pdata_snap.type_mapping.push_back("A");
         pdata_snap.type_mapping.push_back("B");
         pdata_snap.resize(8);
-        pdata_snap.pos[0] = vec3<Scalar>(-0.5, -0.5, -0.5);
-        pdata_snap.pos[1] = vec3<Scalar>(0.5, -0.5, -0.5);
-        pdata_snap.pos[2] = vec3<Scalar>(-0.5, 0.5, -0.5);
-        pdata_snap.pos[3] = vec3<Scalar>(0.5, 0.5, -0.5);
-        pdata_snap.pos[4] = vec3<Scalar>(-0.5, -0.5, 0.5);
-        pdata_snap.pos[5] = vec3<Scalar>(0.5, -0.5, 0.5);
-        pdata_snap.pos[6] = vec3<Scalar>(-0.5, 0.5, 0.5);
-        pdata_snap.pos[7] = vec3<Scalar>(0.5, 0.5, 0.5);
+        pdata_snap.pos[0] = scale(vec3<Scalar>(-0.5, -0.5, -0.5), ref_box, box);
+        pdata_snap.pos[1] = scale(vec3<Scalar>(0.5, -0.5, -0.5), ref_box, box);
+        pdata_snap.pos[2] = scale(vec3<Scalar>(-0.5, 0.5, -0.5), ref_box, box);
+        pdata_snap.pos[3] = scale(vec3<Scalar>(0.5, 0.5, -0.5), ref_box, box);
+        pdata_snap.pos[4] = scale(vec3<Scalar>(-0.5, -0.5, 0.5), ref_box, box);
+        pdata_snap.pos[5] = scale(vec3<Scalar>(0.5, -0.5, 0.5), ref_box, box);
+        pdata_snap.pos[6] = scale(vec3<Scalar>(-0.5, 0.5, 0.5), ref_box, box);
+        pdata_snap.pos[7] = scale(vec3<Scalar>(0.5, 0.5, 0.5), ref_box, box);
         pdata_snap.type[0] = 0;
         pdata_snap.type[1] = 1;
         pdata_snap.type[2] = 0;
@@ -394,18 +405,18 @@ template<class CL> void celllist_embed_test(std::shared_ptr<ExecutionConfigurati
         }
     snap->mpcd_data.resize(8);
     snap->mpcd_data.type_mapping.push_back("A");
-    snap->mpcd_data.position[0] = vec3<Scalar>(-0.5, -0.5, -0.5);
-    snap->mpcd_data.position[1] = vec3<Scalar>(0.5, -0.5, -0.5);
-    snap->mpcd_data.position[2] = vec3<Scalar>(-0.5, 0.5, -0.5);
-    snap->mpcd_data.position[3] = vec3<Scalar>(0.5, 0.5, -0.5);
-    snap->mpcd_data.position[4] = vec3<Scalar>(-0.5, -0.5, 0.5);
-    snap->mpcd_data.position[5] = vec3<Scalar>(0.5, -0.5, 0.5);
-    snap->mpcd_data.position[6] = vec3<Scalar>(-0.5, 0.5, 0.5);
-    snap->mpcd_data.position[7] = vec3<Scalar>(0.5, 0.5, 0.5);
+    snap->mpcd_data.position[0] = scale(vec3<Scalar>(-0.5, -0.5, -0.5), ref_box, box);
+    snap->mpcd_data.position[1] = scale(vec3<Scalar>(0.5, -0.5, -0.5), ref_box, box);
+    snap->mpcd_data.position[2] = scale(vec3<Scalar>(-0.5, 0.5, -0.5), ref_box, box);
+    snap->mpcd_data.position[3] = scale(vec3<Scalar>(0.5, 0.5, -0.5), ref_box, box);
+    snap->mpcd_data.position[4] = scale(vec3<Scalar>(-0.5, -0.5, 0.5), ref_box, box);
+    snap->mpcd_data.position[5] = scale(vec3<Scalar>(0.5, -0.5, 0.5), ref_box, box);
+    snap->mpcd_data.position[6] = scale(vec3<Scalar>(-0.5, 0.5, 0.5), ref_box, box);
+    snap->mpcd_data.position[7] = scale(vec3<Scalar>(0.5, 0.5, 0.5), ref_box, box);
     std::shared_ptr<SystemDefinition> sysdef(new SystemDefinition(snap, exec_conf));
 
     std::shared_ptr<mpcd::ParticleData> pdata_8 = sysdef->getMPCDParticleData();
-    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, 1.0, false));
+    std::shared_ptr<mpcd::CellList> cl(new CL(sysdef, make_uint3(2, 2, 2), false));
     cl->compute(0);
 
         // at first, there is no embedded particle, so everything should just look like the test
@@ -560,7 +571,7 @@ template<class CL> void celllist_embed_test(std::shared_ptr<ExecutionConfigurati
         ArrayHandle<Scalar4> h_embed_pos(embed_pdata->getPositions(),
                                          access_location::host,
                                          access_mode::overwrite);
-        h_embed_pos.data[1] = make_scalar4(0.5, 0.5, -0.5, __int_as_scalar(1));
+        h_embed_pos.data[1] = scale(make_scalar4(0.5, 0.5, -0.5, __int_as_scalar(1)), ref_box, box);
         }
     cl->compute(2);
         // now there should be a second particle in the cell
@@ -654,57 +665,217 @@ template<class CL> void celllist_embed_test(std::shared_ptr<ExecutionConfigurati
 //! dimension test case for MPCD CellList class
 UP_TEST(mpcd_cell_list_dimensions)
     {
-    celllist_dimension_test<mpcd::CellList>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::CPU)));
+    celllist_dimension_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.0, 8.0, 10.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! dimension test case for MPCD CellList class, noncubic
+UP_TEST(mpcd_cell_list_dimensions_noncubic)
+    {
+    celllist_dimension_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(8.5, 10.1, 5.9),
+        make_scalar3(0, 0, 0));
+    }
+
+//! dimension test case for MPCD CellList class, triclinic
+UP_TEST(mpcd_cell_list_dimensions_triclinic)
+    {
+    celllist_dimension_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.0, 8.0, 10.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 
 //! small system test case for MPCD CellList class
 UP_TEST(mpcd_cell_list_small_test)
     {
-    celllist_small_test<mpcd::CellList>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::CPU)));
+    celllist_small_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(2.0, 2.0, 2.0),
+        make_scalar3(0, 0, 0));
     }
 
-//! small system test case for MPCD CellList class
+//! small system test case for MPCD CellList class, noncubic
+UP_TEST(mpcd_cell_list_small_test_noncubic)
+    {
+    celllist_small_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(3.1, 4.2, 5.3),
+        make_scalar3(0, 0, 0));
+    }
+
+//! small system test case for MPCD CellList class, triclinic
+UP_TEST(mpcd_cell_list_small_test_triclinic)
+    {
+    celllist_small_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(2.0, 2.0, 2.0),
+        make_scalar3(0.5, -0.75, 1.0));
+    }
+
+//! grid shift test case for MPCD CellList class
 UP_TEST(mpcd_cell_list_grid_shift_test)
     {
-    celllist_grid_shift_test<mpcd::CellList>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::CPU)));
+    celllist_grid_shift_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.0, 6.0, 6.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! grid shift test case for MPCD CellList class, noncubic
+UP_TEST(mpcd_cell_list_grid_shift_test_noncubic)
+    {
+    celllist_grid_shift_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.1, 7.2, 8.3),
+        make_scalar3(0, 0, 0));
+    }
+
+//! grid shift test case for MPCD CellList class, triclinic
+UP_TEST(mpcd_cell_list_grid_shift_test_triclinic)
+    {
+    celllist_grid_shift_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(6.0, 6.0, 6.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 
 //! embedded particle test case for MPCD CellList class
 UP_TEST(mpcd_cell_list_embed_test)
     {
-    celllist_embed_test<mpcd::CellList>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::CPU)));
+    celllist_embed_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(2.0, 2.0, 2.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! embedded particle test case for MPCD CellList class, noncubic
+UP_TEST(mpcd_cell_list_embed_test_noncubic)
+    {
+    celllist_embed_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(3.1, 4.2, 5.3),
+        make_scalar3(0, 0, 0));
+    }
+
+//! embedded particle test case for MPCD CellList class, triclinic
+UP_TEST(mpcd_cell_list_embed_test_triclinic)
+    {
+    celllist_embed_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::CPU),
+        make_scalar3(2.0, 2.0, 2.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 
 #ifdef ENABLE_HIP
 //! dimension test case for MPCD CellListGPU class
 UP_TEST(mpcd_cell_list_gpu_dimensions)
     {
-    celllist_dimension_test<mpcd::CellListGPU>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::GPU)));
+    celllist_dimension_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.0, 8.0, 10.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! dimension test case for MPCD CellListGPU class, noncubic
+UP_TEST(mpcd_cell_list_gpu_dimensions_noncubic)
+    {
+    celllist_dimension_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(8.5, 10.1, 5.9),
+        make_scalar3(0, 0, 0));
+    }
+
+//! dimension test case for MPCD CellList class, triclinic
+UP_TEST(mpcd_cell_list_gpu_dimensions_triclinic)
+    {
+    celllist_dimension_test<mpcd::CellList>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.0, 8.0, 10.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 
 //! small system test case for MPCD CellListGPU class
 UP_TEST(mpcd_cell_list_gpu_small_test)
     {
-    celllist_small_test<mpcd::CellListGPU>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::GPU)));
+    celllist_small_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(2.0, 2.0, 2.0),
+        make_scalar3(0, 0, 0));
     }
 
-//! small system test case for MPCD CellListGPU class
+//! small system test case for MPCD CellListGPU class, noncubic
+UP_TEST(mpcd_cell_list_gpu_small_test_noncubic)
+    {
+    celllist_small_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(3.1, 4.2, 5.3),
+        make_scalar3(0, 0, 0));
+    }
+
+//! small system test case for MPCD CellListGPU class, triclinic
+UP_TEST(mpcd_cell_list_gpu_small_test_triclinic)
+    {
+    celllist_small_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(2.0, 2.0, 2.0),
+        make_scalar3(0.5, -0.75, 1.0));
+    }
+
+//! grid shift test case for MPCD CellListGPU class
 UP_TEST(mpcd_cell_list_gpu_grid_shift_test)
     {
-    celllist_grid_shift_test<mpcd::CellListGPU>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::GPU)));
+    celllist_grid_shift_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.0, 6.0, 6.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! grid shift test case for MPCD CellListGPU class, noncubic
+UP_TEST(mpcd_cell_list_gpu_grid_shift_test_noncubic)
+    {
+    celllist_grid_shift_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.1, 7.2, 8.3),
+        make_scalar3(0, 0, 0));
+    }
+
+//! grid shift test case for MPCD CellListGPU class, triclinic
+UP_TEST(mpcd_cell_list_gpu_grid_shift_test_triclinic)
+    {
+    celllist_grid_shift_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(6.0, 6.0, 6.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 
 //! embedded particle test case for MPCD CellListGPU class
 UP_TEST(mpcd_cell_list_gpu_embed_test)
     {
-    celllist_embed_test<mpcd::CellListGPU>(std::shared_ptr<ExecutionConfiguration>(
-        new ExecutionConfiguration(ExecutionConfiguration::GPU)));
+    celllist_embed_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(2.0, 2.0, 2.0),
+        make_scalar3(0, 0, 0));
+    }
+
+//! embedded particle test case for MPCD CellListGPU class, noncubic
+UP_TEST(mpcd_cell_list_gpu_embed_test_noncubic)
+    {
+    celllist_embed_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(3.1, 4.2, 5.3),
+        make_scalar3(0, 0, 0));
+    }
+
+//! embedded particle test case for MPCD CellListGPU class, triclinic
+UP_TEST(mpcd_cell_list_gpu_embed_test_triclinic)
+    {
+    celllist_embed_test<mpcd::CellListGPU>(
+        std::make_shared<ExecutionConfiguration>(ExecutionConfiguration::GPU),
+        make_scalar3(2.0, 2.0, 2.0),
+        make_scalar3(0.5, -0.75, 1.0));
     }
 #endif // ENABLE_HIP
diff --git a/hoomd/mpcd/test/parallel_plate_geometry_filler_test.cc b/hoomd/mpcd/test/parallel_plate_geometry_filler_test.cc
index 42b687c618..b450efef2b 100644
--- a/hoomd/mpcd/test/parallel_plate_geometry_filler_test.cc
+++ b/hoomd/mpcd/test/parallel_plate_geometry_filler_test.cc
@@ -27,7 +27,7 @@ void parallel_plate_fill_basic_test(std::shared_ptr<ExecutionConfiguration> exec
     std::shared_ptr<SystemDefinition> sysdef(new SystemDefinition(snap, exec_conf));
 
     auto pdata = sysdef->getMPCDParticleData();
-    auto cl = std::make_shared<mpcd::CellList>(sysdef, 2.0, false);
+    auto cl = std::make_shared<mpcd::CellList>(sysdef, make_uint3(10, 10, 10), false);
     UP_ASSERT_EQUAL(pdata->getNVirtual(), 0);
 
     // create slit channel with half width 5
@@ -111,7 +111,7 @@ void parallel_plate_fill_basic_test(std::shared_ptr<ExecutionConfiguration> exec
      * Change the cell size so that we lie exactly on a boundary.
      */
     pdata->removeVirtualParticles();
-    cl->setCellSize(1.0);
+    cl->setGlobalDim(make_uint3(20, 20, 20));
     filler->fill(2);
     // volume to fill is now from 5->5.5 on + side, same other parameters
     UP_ASSERT_EQUAL(pdata->getNVirtual(), 2 * (20 * 20 / 2) * 2);
@@ -134,7 +134,7 @@ void parallel_plate_fill_basic_test(std::shared_ptr<ExecutionConfiguration> exec
     /*
      * Test the average properties of the virtual particles.
      */
-    cl->setCellSize(2.0);
+    cl->setGlobalDim(make_uint3(10, 10, 10));
     unsigned int N_lo(0), N_hi(0);
     Scalar3 v_lo = make_scalar3(0, 0, 0);
     Scalar3 v_hi = make_scalar3(0, 0, 0);
diff --git a/hoomd/mpcd/test/planar_pore_geometry_filler_test.cc b/hoomd/mpcd/test/planar_pore_geometry_filler_test.cc
index bb3b9a6aa9..bc027a5e28 100644
--- a/hoomd/mpcd/test/planar_pore_geometry_filler_test.cc
+++ b/hoomd/mpcd/test/planar_pore_geometry_filler_test.cc
@@ -27,7 +27,7 @@ void planar_pore_fill_basic_test(std::shared_ptr<ExecutionConfiguration> exec_co
     std::shared_ptr<SystemDefinition> sysdef(new SystemDefinition(snap, exec_conf));
 
     auto pdata = sysdef->getMPCDParticleData();
-    auto cl = std::make_shared<mpcd::CellList>(sysdef, 2.0, false);
+    auto cl = std::make_shared<mpcd::CellList>(sysdef, make_uint3(10, 10, 10), false);
     UP_ASSERT_EQUAL(pdata->getNVirtual(), 0);
 
     // create slit pore channel with half width 5, half length 8
@@ -125,7 +125,7 @@ void planar_pore_fill_basic_test(std::shared_ptr<ExecutionConfiguration> exec_co
      * Now, all sides of the U have thickness 0.5.
      */
     pdata->removeVirtualParticles();
-    cl->setCellSize(1.0);
+    cl->setGlobalDim(make_uint3(20, 20, 20));
     filler->fill(2);
     UP_ASSERT_EQUAL(pdata->getNVirtual(),
                     (unsigned int)(4 * 2 * (0.5 * 4.5 + 0.5 * 16 + 0.5 * 4.5) * 20));
@@ -152,7 +152,7 @@ void planar_pore_fill_basic_test(std::shared_ptr<ExecutionConfiguration> exec_co
      * Test the average fill properties of the virtual particles.
      */
     filler->setDensity(2.0);
-    cl->setCellSize(2.0);
+    cl->setGlobalDim(make_uint3(10, 10, 10));
     unsigned int N_avg(0);
     Scalar3 v_avg = make_scalar3(0, 0, 0);
     Scalar T_avg(0);
diff --git a/hoomd/mpcd/test/utils.h b/hoomd/mpcd/test/utils.h
index 51ec3e97b0..cf39103485 100644
--- a/hoomd/mpcd/test/utils.h
+++ b/hoomd/mpcd/test/utils.h
@@ -4,6 +4,7 @@
 #ifndef MPCD_TEST_UTILS_H_
 #define MPCD_TEST_UTILS_H_
 
+#include "hoomd/mpcd/CellList.h"
 #include "hoomd/mpcd/CellThermoCompute.h"
 
 namespace hoomd
@@ -50,6 +51,29 @@ class AllThermoRequest
     std::shared_ptr<mpcd::CellThermoCompute> m_thermo;
     };
 
+Scalar3 scale(const Scalar3& ref_pos, std::shared_ptr<BoxDim> ref_box, std::shared_ptr<BoxDim> box)
+    {
+    return box->makeCoordinates(ref_box->makeFraction(ref_pos));
+    }
+
+Scalar4 scale(const Scalar4& ref_pos, std::shared_ptr<BoxDim> ref_box, std::shared_ptr<BoxDim> box)
+    {
+    const Scalar3 pos = scale(make_scalar3(ref_pos.x, ref_pos.y, ref_pos.z), ref_box, box);
+    return make_scalar4(pos.x, pos.y, pos.z, ref_pos.w);
+    }
+
+vec3<Scalar>
+scale(const vec3<Scalar>& ref_pos, std::shared_ptr<BoxDim> ref_box, std::shared_ptr<BoxDim> box)
+    {
+    return vec3<Scalar>(scale(vec_to_scalar3(ref_pos), ref_box, box));
+    }
+
+unsigned int make_local_cell(std::shared_ptr<mpcd::CellList> cl, int ix, int iy, int iz)
+    {
+    const int3 cell = cl->getLocalCell(make_int3(ix, iy, iz));
+    return cl->getCellIndexer()(cell.x, cell.y, cell.z);
+    }
+
     } // end namespace hoomd
 
 #endif // MPCD_TEST_UTILS_H_