mcMd/neighbor/CellList.h

#ifndef MCMD_CELL_LIST_H
#define MCMD_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 "Cell.h"
#include "CellTag.h"
#include <mcMd/chemistry/Atom.h>
#include <simp/boundary/Boundary.h>
#include <util/space/IntVector.h>
#include <util/containers/DArray.h>
#include <util/containers/FSArray.h>
#include <util/global.h>

#include <sstream>

class CellListTest;

namespace McMd
{

   using namespace Util;
   using namespace Simp;

   class CellList
   {

   public:

      // Static members

      static const int MaxNeighbor = 27*Cell::MaxAtomCell;

      typedef FSArray<Atom*, MaxNeighbor> NeighborArray;

      // Public member functions

      CellList();

      virtual ~CellList();

      void setAtomCapacity(int atomCapacity);

      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);

      void setup(const Boundary &boundary, double cutoff);

      void clear();

      int cellIndexFromPosition(const Vector &pos) const;

      void addAtom(Atom &atom);

      void deleteAtom(Atom &atom);

      void updateAtomCell(Atom &atom, const Vector &pos);

      void
      getNeighbors(const Vector &pos, NeighborArray &neighbors) const;

      void 
      getCellNeighbors(int ic, NeighborArray &neighbors, int &nInCell) const;

      int gridDimension(int i) const;

      int totCells() const;

      int nAtom() const;

      int atomCapacity() const;

      bool isValid(int nAtom=-1) const;

   private:

      DArray<Cell>    cells_;

      DArray<CellTag> cellTags_;

      Vector lengths_;      

      Vector invCellWidths_;   

      IntVector minDel_;       

      IntVector maxDel_;       

      IntVector numCells_;     

      IntVector minCells_; 

      IntVector maxCells_; 

      int  YZCells_;        
 
      int  totCells_;        

      int  atomCapacity_;        

      const Boundary* boundaryPtr_; 

      void setCellsAxis(int axis, double cutoff);

      int shiftCellCoordAxis(int i, int x) const;

      int cellIndexFromCoord(int cx, int cy, int cz) const;

      void cellCoordFromIndex(int i, int &cx, int &cy, int &cz) const;

      bool isValidAtomId(int atomId)
      { return ( (0 <= atomId) && (atomId < cellTags_.capacity()) ); }


   //friends:

      friend class ::CellListTest;

   }; 


   // Public inline function definitions:

   /*
   * Get total number of cells in this CellList.
   */
   inline int CellList::totCells() const
   { return totCells_; }

   /*
   *  Return index of the cell that contains the position array pos = {x, y, z}.
   */
   inline int CellList::cellIndexFromPosition(const Vector &pos) const
   {
      int cx, cy, cz;

      {
         Vector posG;
         boundaryPtr_->transformCartToGen(pos, posG);
         cx = int(posG[0]*invCellWidths_[0]);
         cy = int(posG[1]*invCellWidths_[1]);
         cz = int(posG[2]*invCellWidths_[2]);
      }

      assert(cx >= 0);
      assert(cx < numCells_[0]);
      assert(cy >= 0);
      assert(cy < numCells_[1]);
      assert(cz >= 0);
      assert(cz < numCells_[2]);

      return cz + cy*numCells_[2] + cx*YZCells_;
   }

   /*
   * Remove an Atom object from its Cell.
   */
   inline void CellList::deleteAtom(Atom &atom)
   {
      int atomId   = atom.id();
      assert(isValidAtomId(atomId));
      CellTag& cellTag = cellTags_[atomId];
      cells_[cellTag.cellId].deleteAtom(cellTag);
   }

   /*
   * Add a Atom to the appropriate cell, based on its position.
   */
   inline void CellList::addAtom(Atom &atom)
   {
      Vector position = atom.position();
      int    cellId   = cellIndexFromPosition(position);
      int    atomId   = atom.id();
      assert(isValidAtomId(atomId));
      cells_[cellId].addAtom(cellTags_[atomId], atom, cellId);
   }

   /*
   * Update CellList to reflect new atom position
   */
   inline void CellList::updateAtomCell(Atom &atom, const Vector &pos)
   {
      int atomId = atom.id();
      assert(isValidAtomId(atomId));
      CellTag& cellTag = cellTags_[atomId];
      int      oldCell = cellTag.cellId;
      int      newCell = cellIndexFromPosition(pos);
      if (oldCell != newCell) {
         cells_[oldCell].deleteAtom(cellTag);
         cells_[newCell].addAtom(cellTag, atom, newCell);
      }
   }

   /*
   * Fill an array with Ids of atoms that are near a specified position.
   */
   inline void
   CellList::getNeighbors(const Vector &pos, NeighborArray &neighbors) const
   {
      int nInCell;
      int ic = cellIndexFromPosition(pos);
      getCellNeighbors(ic, neighbors, nInCell);
   }

   // Private inline member functions

   /*
   * Return value of integer cell coordinate x for axis i shifted to the range
   * minCells_[i] <= x <= maxCells_[i].
   */
   inline int CellList::shiftCellCoordAxis(int i, int x) const
   {
      if (x > maxCells_[i]) return x - numCells_[i];
      if (x < minCells_[i]) return x + numCells_[i];
      return x;
   }

   /*
   * Return integer index of a cell with integer coordinates (cx,cy,cz).
   */
   inline int CellList::cellIndexFromCoord(int cx, int cy, int cz) const
   {
      return cx*YZCells_ + cy*numCells_[2] + cz ;
   }

   /*
   * Return (by reference) int coordinates (cx,cy,cz) of cell with index i.
   */
   inline
   void CellList::cellCoordFromIndex(int i, int &cx, int &cy, int &cz) const
   {
      cx = (int)(i / YZCells_);
      i = i - cx*YZCells_;
      cy = (int)(i / numCells_[2]);
      cz = i - cy*numCells_[2];
   }

   inline int CellList::gridDimension(int i) const
   {  return numCells_[i]; }

   /*
   * Serialize to/from an Archive.
   */
   template <class Archive>
   void CellList::serialize(Archive& ar, const unsigned int version)
   {
      ar & lengths_;      
      ar & invCellWidths_;   
      ar & minDel_;       
      ar & maxDel_;       
      ar & numCells_;     
      ar & minCells_; 
      ar & maxCells_; 
      ar & YZCells_;        
      ar & totCells_;   
      ar & atomCapacity_;
      if (ar.is_loading()) {
         cells_.allocate(totCells_);
         cellTags_.allocate(atomCapacity_);
      }
      clear();
   }

} 
#endif