1 /* -*- c++ -*- */
2 #ifndef CELLLISTENUMERATOR_STANDARD_H
3 #define CELLLISTENUMERATOR_STANDARD_H
4 
5 #include "CellListEnumerator.h"
6 
7 #include "CubicCellManager.h"
8 #include "VacuumBoundaryConditions.h"
9 #include "Topology.h"
10 #include <vector>
11 
12 //#define DEBUG_CELLLISTENUMERATOR_STANDARD
13 namespace ProtoMol {
14   //_________________________________________________________________ CellListEnumerator_standard
15 
16   /**
17    * Specialization of the cell enumerator for vacuum and cubic cell manager
18    */
19   template<>
20   class CellListEnumerator<VacuumBoundaryConditions,CubicCellManager> {
21   public:
22     //typedef pair<int,int> CellPair; // The first atom from each cell list in the pair
23     struct CellPair {int first; int second;};
24 
25   public:
CellListEnumerator()26     CellListEnumerator():myCutoff(-1.0),
27 			 myCellSize(Vector3D(-1.0,-1.0,-1.0)),
28 			 myMax(CubicCellManager::Cell(-1,-1,-1)){}
29 
initialize(const Topology<VacuumBoundaryConditions,CubicCellManager> * topo,Real cutoff)30     void initialize(const Topology<VacuumBoundaryConditions,CubicCellManager>* topo,
31 		    Real cutoff) {
32 
33       the_beginning=topo->cellLists.begin();
34       the_end=topo->cellLists.end();
35 
36       i=the_beginning;
37       j=i;
38 
39       the_first = 0;
40       the_counter = the_first;
41       myCellListStruct = &(topo->cellLists);
42 
43       if(myCutoff != cutoff ||
44 	 topo->cellManager.getRealCellSize() != myCellSize ||
45 	 !(myMax == topo->cellManager.findCell(topo->max-topo->min))){
46 
47 
48 	myCellSize = topo->cellManager.getRealCellSize();
49 	myCutoff=cutoff;
50 
51 	myMax = topo->cellManager.findCell(topo->max-topo->min);
52 
53 // 	Report::report <<"New: "<<myMax.x<<","<< myMax.y<<","<< myMax.z<<Report::endr;
54 	Real cutoff2 = cutoff*cutoff;
55 
56 	CubicCellManager::Cell zero(0,0,0);
57 	myDeltaList.clear();
58 	int nx = (int)(cutoff/myCellSize.x+1.0+Constant::EPSILON);
59 	int ny = (int)(cutoff/myCellSize.y+1.0+Constant::EPSILON);
60 	int nz = (int)(cutoff/myCellSize.z+1.0+Constant::EPSILON);
61 	// Do not consider deltas bigger than the dimesion of the simulation box
62 	int n0 = std::min(nx,topo->cellLists.getDimX()-1);
63 	int n1 = std::min(ny,topo->cellLists.getDimY()-1);
64 	int n2 = std::min(nz,topo->cellLists.getDimZ()-1);
65 	Real xx = myCellSize.x*myCellSize.x;
66 	Real yy = myCellSize.y*myCellSize.y;
67 	Real zz = myCellSize.z*myCellSize.z;
68 	for(int k =-n0;k<=n0;k++){
69 	  int x=abs(k)-1; if(x<0) x=0;
70 	  Real d0 = x*x*xx;
71 	  for(int l =-n1;l<=n1;l++){
72 	    int y=abs(l)-1; if(y<0) y=0;
73 	    Real d1 = d0+y*y*yy;
74 	    for(int m =-n2;m<=n2;m++){
75 	      int z=abs(m)-1; if(z<0) z=0;
76 	      if(d1+z*z*zz<cutoff2) {
77 		CubicCellManager::Cell delta(k,l,m);
78 		if(zero < delta)
79 		  myDeltaList.push_back(delta);
80 	      }
81 	    }
82 	  }
83 	}
84 	std::sort(myDeltaList.begin(),myDeltaList.end());
85 // 	for(unsigned int i=0;i<myDeltaList.size();i++){
86 // 	  Report::report <<"Delta["<<i<<"] "<< myDeltaList[i].x<<","<< myDeltaList[i].y<<","<< myDeltaList[i].z<<Report::endr;
87 // 	}
88 
89 	//Report::report << Report::debug(2) << "Cell list algorithm: "<< topo->cellLists.getDimX()*topo->cellLists.getDimY()*topo->cellLists.getDimZ()<<", "<<myDeltaList.size()<<", "<< topo->cellLists.getDimX()*topo->cellLists.getDimY()*topo->cellLists.getDimZ()*myCellSize.x*myCellSize.y*myCellSize.z << ", "<< myDeltaList.size()*myCellSize.x*myCellSize.y*myCellSize.z <<Report::endr;
90       }
91       the_last = myDeltaList.size();
92     }
93 
94     /// retrieve the current cell pair
get(CellPair & cp)95     void get(CellPair &cp) {cp.first=i->second; cp.second=j->second; }
96     /// if the cells of the current cell pair are the same
notSameCell()97     bool notSameCell(){return (i != j);}
98     /// reached the end of the list of cell pairs
done()99     bool done() {return (i==the_end);}
100     /// goto the end of pair of cells with the same first cell
gotoEndPair()101     void gotoEndPair() {j=the_end;};
102 
103     /// advance by inc in the cell list in respect to first cell
nextNewPair(int inc)104     void nextNewPair(int inc){
105       if(inc < 1)
106 	return;
107       while(inc > 0 && i != the_end){
108 	++i;
109 	--inc;
110 	while(i != the_end && i->second < 0)
111 	  ++i;
112       }
113       j=i;
114     }
115 
116     /// get next pair
next()117     void next(){
118       if(i!=the_end){
119 	j=the_end;
120 	while(the_counter != the_last &&
121 	      the_end == (j=myCellListStruct->find(i->first+myDeltaList[the_counter++])));
122 
123 	if(j==the_end){
124 	  ++i;
125 	  while(i != the_end && i->second < 0)
126 	    ++i;
127 	  j=i;
128 	  the_counter = the_first;
129 	}
130       }
131     }
132 
133     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
134     // My data members
135     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
136   private:
137     CubicCellManager::CellListStructure::const_iterator i,j,the_beginning,the_end;
138     Real myCutoff;
139     std::vector<CubicCellManager::Cell> myDeltaList;
140     int the_first, the_last, the_counter;
141     const CubicCellManager::CellListStructure *myCellListStruct;
142     Vector3D myCellSize;
143     CubicCellManager::Cell myMax;
144   };
145 }
146 #endif /* CELLLISTENUMERATOR_STANDARD_H */
147