tools/neighbor/CellList.h

#ifndef TOOLS_CELL_LIST_H
#define TOOLS_CELL_LIST_H

/*
* Simpatico - Simulation Package for Polymeric and Molecular Liquids
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <tools/neighbor/Cell.h>
#include <tools/neighbor/CellAtom.h>
#include <tools/chemistry/Atom.h>
#include <simp/boundary/Boundary.h>
#include <util/space/Grid.h>
#include <util/containers/DArray.h>
#include <util/containers/GArray.h>
#include <util/containers/FArray.h>
#include <util/global.h>

namespace Tools
{

   using namespace Util;
   using namespace Simp;

   class CellList
   {

   public:

      // Public methods

      CellList();

      virtual ~CellList();

      void allocate(int atomCapacity, const Vector& lower, const Vector& upper, 
                    const Vector& cutoffs);

      void allocate(int atomCapacity, const Vector& lower, const Vector& upper, 
                    double cutoff);

      void 
      makeGrid(const Vector& lower, const Vector& upper, const Vector& cutoffs);

      void placeAtom(Atom &atom);

      void build();

      void update();

      void clear();

      const Grid& grid() const;

      double cellLength(int i) const;

      int cellIndexFromPosition(const Vector& position) const;

      const Cell* begin() const;

      const Cell& cell(int i) const;

      int nAtom() const;

      int nReject() const;

      #ifdef UTIL_DEBUG

      int maxNAtomCell() const;
      #endif

      int atomCapacity() const;

      int cellCapacity() const;

      bool isAllocated() const;

      bool isBuilt() const;

      bool isValid() const;

   private:

      /*
      * Temporary storage for atom pointers, before copying to cells. 
      */
      struct Tag {
         Atom* ptr;
         int cellRank;
      };

      Grid grid_;

      DArray<Tag> tags_;

      DArray<CellAtom> atoms_;

      GArray<Cell> cells_;

      Vector lower_;

      Vector upper_;

      Vector cellLengths_;

      Cell* begin_;

      int nAtom_;

      int nReject_;

      #ifdef UTIL_DEBUG
      int maxNAtomCell_;
      #endif

      bool isBuilt_;

      int gridOffsets_ [3];

      void setGridDimensions(const Vector& lower, const Vector& upper,
                             const Vector& cutoffs);

      bool isValidAtomId(int atomId);

      int calculateAxisOffset(int cellId, int i, int x) const;
   };

   // Private inline method definitions:

   /*
   * Return value of offset for cell a distance x along axis i shifted to the range
   * 0 <= x <= grid_.Dimension(i)
   */
   inline int CellList::calculateAxisOffset(int cellId, int i, int x) const
   {
      // mod the shift so it's within the grid dimension
      x %= grid_.dimension(i);

      // calculate cell's position within the grid
      int p[3]; int j;
      for (j = 0; j < Dimension -1; ++j) {
         p[j] = cellId/gridOffsets_[j];
         cellId -= p[j]*gridOffsets_[j];
      }
      p[j] = cellId;

      if (p[i] + x >= grid_.dimension(i)) return (x - grid_.dimension(i))*gridOffsets_[i];
      if (p[i] + x < 0) return (grid_.dimension(i) + x)*gridOffsets_[i];
      return x*gridOffsets_[i];
   }

   // Public inline method definitions:

   /*
   * Identify the cell within which a position lies.
   */
   inline int CellList::cellIndexFromPosition(const Vector& position) const
   {
      IntVector r;
      for (int i = 0; i < Dimension; ++i) {
         if (position[i] <= lower_[i]) {
            return -1;
         }
         if (position[i] >= upper_[i]) {
            return -1;
         }
         r[i] = int( (position[i] - lower_[i]) / cellLengths_[i] );
         assert(r[i] < grid_.dimension(i));
      }
      return grid_.rank(r);
   }

   /*
   * Add an atom to the appropriate cell, based on its position.
   */
   inline void CellList::placeAtom(Atom &atom)
   {
      // Preconditon
      assert(nAtom_ < tags_.capacity());

      int rank = cellIndexFromPosition(atom.position);
      if (rank >= 0) {
         tags_[nAtom_].cellRank = rank;
         tags_[nAtom_].ptr = &atom;
         cells_[rank].incrementCapacity();
         ++nAtom_;
      } else {
         ++nReject_;
      }
   }

   /*
   * Return true iff atomId is valid, i.e., if 0 <= 0 < atomCapacity.
   */
   inline bool CellList::isValidAtomId(int atomId)
   {  return ( (0 <= atomId) && (atomId < tags_.capacity()) ); }

   /*
   * Return associated Grid object.
   */
   inline const Grid& CellList::grid() const
   {  return grid_; }

   /*
   * Return length of each cell in direction i.
   */
   inline double CellList::cellLength(int i) const
   {  return cellLengths_[i]; }

   /*
   * Return reference to cell number i.
   */
   inline const Cell& CellList::cell(int i) const
   {  return cells_[i]; }

   /*
   * Return pointer to first Cell.
   */
   inline const Cell* CellList::begin() const
   {  return begin_; }

   /*
   * Is this CellList allocated?
   */
   inline bool CellList::isAllocated() const
   {  return (cells_.capacity() > 0); }

}
#endif