1 // -*- C++ -*- 2 // NeighborCellLocator.h: Cell-based algorithm for finding neighbors. 3 // 4 // Copyright (C) 2008-2011 Jakob Schiotz and Center for Individual 5 // Nanoparticle Functionality, Department of Physics, Technical 6 // University of Denmark. Email: schiotz@fysik.dtu.dk 7 // 8 // This file is part of Asap version 3. 9 // Asap is released under the GNU Lesser Public License (LGPL) version 3. 10 // However, the parts of Asap distributed within the OpenKIM project 11 // (including this file) are also released under the Common Development 12 // and Distribution License (CDDL) version 1.0. 13 // 14 // This program is free software: you can redistribute it and/or 15 // modify it under the terms of the GNU Lesser General Public License 16 // version 3 as published by the Free Software Foundation. Permission 17 // to use other versions of the GNU Lesser General Public License may 18 // granted by Jakob Schiotz or the head of department of the 19 // Department of Physics, Technical University of Denmark, as 20 // described in section 14 of the GNU General Public License. 21 // 22 // This program is distributed in the hope that it will be useful, 23 // but WITHOUT ANY WARRANTY; without even the implied warranty of 24 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 25 // GNU General Public License for more details. 26 // 27 // You should have received a copy of the GNU General Public License 28 // and the GNU Lesser Public License along with this program. If not, 29 // see <http://www.gnu.org/licenses/>. 30 31 32 #ifndef NEIGHBORCELLLOCATOR 33 #define NEIGHBORCELLLOCATOR 34 35 #include "AsapPython.h" 36 #include "Asap.h" 37 #include "Vec.h" 38 #include "NeighborLocator.h" 39 #include "IVec.h" 40 #include "Templates.h" 41 #include <vector> 42 using std::vector; 43 #include <map> 44 #include <utility> 45 using std::pair; 46 #include <math.h> 47 #include <stdint.h> 48 49 #define COMPACT_NBLIST 50 51 namespace ASAPSPACE { 52 53 typedef int translationsidx_t; // Could also be unsigned char 54 55 // The datatpe neighboritem_t is mainly used in NeigbhorList.cpp 56 // It is a combination of an index into the atoms and an index 57 // into a list of translation vectors. 58 #ifdef COMPACT_NBLIST 59 // A neighborlist item is a single 32-bit unsigned integer, 60 // the first 5 bits are the tranlation vector index, the last 27 are 61 // the actual neighborlist index. 62 typedef uint32_t neighboritem_t; 63 #define NEIGHBOR_ITEM_MASK ((1 << 27) - 1) 64 #define NEIGHBOR_XLAT_SHIFT 27 65 #define NEIGHBOR_INDEX(x) ((x) & NEIGHBOR_ITEM_MASK) 66 #define NEIGHBOR_XLAT(x) ((x) >> NEIGHBOR_XLAT_SHIFT) 67 #define NEIGHBORITEM_PACK(x,a,b) x = a | (b << NEIGHBOR_XLAT_SHIFT) 68 #else // COMPACT_NBLIST 69 // A neighborlist item is a pair of an int (the actual neighborlist 70 // index) and a translationsidx_t, the translation vector index. 71 typedef pair<int,translationsidx_t> neighboritem_t; 72 #define NEIGHBOR_INDEX(x) (x).first 73 #define NEIGHBOR_XLAT(x) (x).second 74 #define NEIGHBORITEM_PACK(x,a,b) x; (x).first = a; (x).second = b; 75 #endif // COMPACT_NBLIST 76 77 PyAsap_NeighborLocatorObject *PyAsap_NewNeighborCellLocator(Atoms *a, 78 double rCut, double driftfactor = 0.05); 79 80 class NeighborCellLocator : public NeighborLocator 81 { 82 protected: 83 /// Generate a neighbor list for atoms a with cutoff rCut. 84 85 /// The neighbor list will contain all neighbors within the distance 86 /// rCut. The neighborlist can be reused until an atom has moved 87 /// more than rCut*driftfactor. 88 NeighborCellLocator(Atoms *a, double rCut, double driftfactor = 0.05); 89 virtual ~NeighborCellLocator(); 90 91 friend PyAsap_NeighborLocatorObject *PyAsap_NewNeighborCellLocator( 92 Atoms *a, double rCut, double driftfactor); 93 94 friend void PyAsap_Dealloc<PyAsap_NeighborLocatorObject>(PyObject *self); 95 96 public: 97 /// Enable neighbors of ghost atoms by calling this just after the constructor 98 void EnableNeighborsOfGhosts(); 99 100 /// Check if the neighbor list can still be reused, update if not. 101 bool CheckAndUpdateNeighborList(); 102 103 /// Check if the neighbor list can still be reused, update if not. 104 /// 105 /// This version is used when called from Python 106 virtual bool CheckAndUpdateNeighborList(PyObject *atoms); 107 108 /// Check the neighbor list. 109 /// 110 /// Check if the neighbor list can still be reused, return true if 111 /// it should be updated. 112 virtual bool CheckNeighborList(); 113 114 /// Update neighbor list 115 virtual void UpdateNeighborList(); 116 117 /// Get wrapped positions of all the atoms GetWrappedPositions()118 const vector<Vec> &GetWrappedPositions() const {ASSERT(wrappedPositionsValid); 119 return wrappedPositions;} 120 GetWrappedPositions(vector<Vec> & wp)121 void GetWrappedPositions(vector<Vec> &wp) const {ASSERT(wrappedPositionsValid); 122 wp.insert(wp.begin(), wrappedPositions.begin(), wrappedPositions.end());} 123 124 /// Get scaled positions of all the atoms 125 const vector<Vec> &GetScaledPositions() const; 126 127 /// Get reference positions (used by NeighborList) GetReferencePositions()128 const Vec *GetReferencePositions() const {return &referencePositions[0];} 129 130 /// Get information about the neighbors of atom n ("half" neighbor list) 131 /// 132 /// Input values: n is the number of the atom. r (optional) is a 133 /// cutoff, must be less than rCut in the constructor (not 134 /// checked!). 135 /// 136 /// In-out values: size contains the maximum space in the arrays. 137 /// It is decremented by the number of neighbors placed in the 138 /// arrays. It is an error to call GetNeighbors with too small a 139 /// value of size. 140 /// 141 /// Out values: neighbors[] contains the numbers of the atoms, 142 /// diffs[] contains the \em relative positions of the atoms, 143 /// diffs2[] contains the norms of the diffs vectors. 144 /// 145 /// Return value: The number of neighbors. 146 int GetNeighbors(int n, int *neighbors, Vec *diffs, double *diffs2, 147 int& size, double r = -1.0) const; 148 149 150 /// Get information about the neighbors of atom n ("half" neighbor list) 151 /// 152 /// This version of GetNeighbors only returns the numbers of the neighbors. 153 /// It is intended for the Python interface. 154 void GetNeighbors(int n, vector<int> &neighbors) const; 155 156 int GetFullNeighbors(int n, int *neighbors, Vec *diffs, double *diffs2, 157 int& size, double r = -1.0) const; 158 159 /// Get information about the neighbors of atom n (full neighbor list) 160 /// 161 /// This version of GetNeighbors only returns the numbers of the neighbors. 162 /// It is intended for the Python interface. 163 void GetFullNeighbors(int n, vector<int> &neighbors) const; 164 165 166 /// Return the neighbors and the corresponding translations. 167 /// 168 /// The vectors are cleared before data is put into them. 169 /// 170 /// Return value: The number of neighbors. 171 int GetListAndTranslations(int n, vector<neighboritem_t> &neighbors) const; 172 173 int GetComplementaryListAndTranslations(int n, 174 vector<neighboritem_t> &neighbors) const; 175 176 /// Remake neighbor lists when a few atoms have been modified. 177 /// 178 /// This version, unlike NeighborList::RemakeList does 179 /// not report back which other atoms have been affected. 180 void RemakeLists_Simple(const set<int> &modified); 181 182 /// Return the guaranteed maximal length of a single atom's NB list. 183 184 /// Call this before using GetNeighbors() to make sure the arrays 185 /// are big enough. The value may change when the neighbor list is 186 /// updated. MaxNeighborListLength()187 int MaxNeighborListLength() const {return maxLength;} 188 189 /// Get the number of atoms in the corresponding list of atoms. GetNumberOfAtoms()190 int GetNumberOfAtoms() const {return nAtoms;} // Used from swig. 191 192 /// Return the cutoff distance (rCut) specified when creating this nblist. GetCutoffRadius()193 double GetCutoffRadius() const {return rCut;} 194 195 /// Return the cutoff distance including twice the drift. 196 /// 197 /// For a NeighborCellLocator, the drift cannot be predicted and is ignored. GetCutoffRadiusWithDrift()198 double GetCutoffRadiusWithDrift() const {return rCut;} 199 200 /// Get a copy of the table of translations (27 entries) 201 void GetTranslationTable(vector<IVec> &table) const; 202 203 /// Normalize the positions and calculate scaled space version 204 /// 205 /// This is used when a neighbor list is updated 206 void ScaleAndNormalizePositions(); 207 208 /// Normalize some positions and calculate scaled space version 209 /// 210 /// The first argument is a set of atoms to be normalized, the 211 /// corresponding scaled positions are placed in scaledpos. 212 213 void ScaleAndNormalizePositions(const set<int> &modified, 214 vector<Vec> &scaledpos); 215 216 /// Normalizing new positions using the normalization already calculated. 217 /// 218 /// This is used when checking a neighborlist. 219 void RenormalizePositions(); 220 221 #if 0 222 /// Renormalize the old positions using old and new positions of the atoms. 223 /// 224 /// This is called by a master neighbor list when it has normalized 225 /// the positions of the atoms. The argument is the positions of 226 /// the atoms *before* renormalization, this method uses them 227 /// together with the current positions to update its own list of 228 /// old positions. 229 void RenormalizeReferencePositions(const vector<Vec> &oldpos); 230 #endif 231 232 /// Update reference positions of some atoms 233 void UpdateReferencePositions(const set<int> &modified); 234 235 /// Return the atoms access object. Used by a few tool functions. GetAtoms()236 virtual Atoms *GetAtoms() const {return atoms;} 237 GetName()238 string GetName() const {return "NeighborCellLocator";} 239 240 /// Print internal info about an atom 241 virtual void print_info(int n); 242 243 /// Print memory usage 244 virtual long PrintMemory() const; 245 246 protected: 247 /// Generate a new neighbor list. 248 virtual void MakeList(); 249 250 /// Make the lists of neighboring cells. 251 void MakeNeighboringCellLists(); 252 253 /// Make translation table 254 void MakeTranslationTable(); 255 256 typedef vector< pair<int, translationsidx_t> > nbcell_t; 257 258 void makeNbCells(int thiscell); 259 260 int CommonGetNeighbors(int n, int *neighbors, Vec *diffs, double *diffs2, 261 int& size, double r, bool wantfull) const; 262 263 void CommonGetNeighbors(int a1, vector<int> &neighbors, 264 bool wantfull) const; 265 266 double get_drift() const; 267 268 protected: 269 Atoms *atoms; ///< A pointer to the atoms. 270 int nAtoms; ///< The number of atoms excluding ghosts. 271 int nAllAtoms; ///< The number of atoms including ghosts. 272 double rCut; ///< The cutoff radius. 273 double rCut2; ///< The square of the cutoff radius. 274 double minboxsize; ///< The minimal box size 275 bool periodic[3]; ///< The boundary conditions. 276 bool oldperiodic[3]; ///< The boundary conditions at the last update. 277 int maxLength; ///< The guaranteed max length of a neighbor list. 278 int nCells[3]; ///< Number of cells 279 int nTotalCells[4]; 280 int nCellsTrue[3]; 281 int nCellsGapStart[3]; 282 int nCellsGapSize[3]; 283 bool includeghostneighbors; ///< Include neighbors to ghost atoms. 284 double size[3]; ///< Size of system in scaled space 285 double minimum[3]; ///< Offset of system in scaled space 286 vector<Vec> referencePositions; ///< Positions at last update. 287 vector<Vec> wrappedPositions; ///< Wrapped positions. 288 vector<Vec> scaledPositions; ///< Scaled positions. 289 vector<Vec> offsetPositions; ///< wrappedPositions - positions. 290 vector<Vec> scaledOffsetPositions; 291 bool scaledPositionsValid; 292 bool wrappedPositionsValid; 293 294 Vec old_inverse_cell[3]; ///< Inverse unit cell of last renormalization. 295 int supercell_counter; ///< When was old_inverse_cell last updated? 296 297 /// The list of cells containing the atoms 298 vector< vector<int> > cells; 299 300 /// The number of the cell to which an atom belongs 301 vector<int> cellIndices; 302 303 /// For each cell, a list of the neighboring cells, and their offset 304 /// across the periodic boundaries. 305 vector<IVec> neighborCellOffsets; 306 307 /// List of neighboring cells, valid for ... 308 nbcell_t nbCells_inside; // ... a cell not touching the boundary. 309 nbcell_t nbCells_left; // ... a cell touching a single boundary. 310 nbcell_t nbCells_right; 311 nbcell_t nbCells_top; 312 nbcell_t nbCells_bottom; 313 nbcell_t nbCells_front; 314 nbcell_t nbCells_back; 315 316 /// Lists of neighboring cells for all cells. Center and sides are 317 /// pointers to the objects above, edges and corners are unique. 318 std::map<int, nbcell_t*> nbCells_all; 319 vector<nbcell_t*> nbCells_onthefly; // On-the-fly elements in nbCells_all 320 // (for memory management). 321 322 /// Table of possible translation vectors 323 vector<IVec> translationTable; 324 }; 325 326 } // end namespace 327 328 #endif // NEIGHBORCELLLOCATOR 329