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