Merge pull request #1950 from glotzerlab/feature-mpcd-triclinic

Allow triclinic boxes in MPCD

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*
 

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,43 +96,33 @@ 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
         {
         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

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,

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,

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,

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: