1 // Copyright (C) 2008-2011 Jakob Schiotz and Center for Individual
2 // Nanoparticle Functionality, Department of Physics, Technical
3 // University of Denmark. Email: schiotz@fysik.dtu.dk
4 //
5 // This file is part of Asap version 3.
6 // Asap is released under the GNU Lesser Public License (LGPL) version 3.
7 // However, the parts of Asap distributed within the OpenKIM project
8 // (including this file) are also released under the Common Development
9 // and Distribution License (CDDL) version 1.0.
10 //
11 // This program is free software: you can redistribute it and/or
12 // modify it under the terms of the GNU Lesser General Public License
13 // version 3 as published by the Free Software Foundation. Permission
14 // to use other versions of the GNU Lesser General Public License may
15 // granted by Jakob Schiotz or the head of department of the
16 // Department of Physics, Technical University of Denmark, as
17 // described in section 14 of the GNU General Public License.
18 //
19 // This program is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // and the GNU Lesser Public License along with this program. If not,
26 // see <http://www.gnu.org/licenses/>.
27
28 // Asap: Copyright (C) 2008 CINF/CAMD and Jakob Schiotz
29
30 #include "NeighborCellLocator.h"
31 #include "Matrix3x3.h"
32 #include "Atoms.h"
33 #include "Exception.h"
34 #include "Timing.h"
35 //#define ASAPDEBUG
36 #include "Debug.h"
37 #include <cstdio>
38 #include <iostream>
39 #include <string.h>
40 using std::cerr;
41 using std::endl;
42 using std::flush;
43
44
45 // Choose between 3x3x3 and 5x5x5 search patterns.
46 #undef PATTERN5
47
48 #define TESTOPTIM
49
NeighborCellLocator(Atoms * a,double rCut,double driftfactor)50 NeighborCellLocator::NeighborCellLocator(Atoms *a, double rCut,
51 double driftfactor)
52 {
53 CONSTRUCTOR;
54 DEBUGPRINT;
55 if (a != NULL)
56 {
57 atoms = a;
58 AsapAtoms_INCREF(atoms);
59 }
60 else
61 atoms = new NormalAtoms();
62 this->rCut = rCut;
63 includeghostneighbors = false; // Only true in SecondaryNeighborLocator
64 rCut2 = rCut * rCut;
65 #ifdef PATTERN5
66 minboxsize = (1 + 2.0 * driftfactor) * rCut / 2.0;
67 #else
68 minboxsize = (1 + 2.0 * driftfactor) * rCut;
69 #endif
70 nCells[0] = nCells[1] = nCells[2] = 0;
71 nTotalCells[0] = nTotalCells[1] = nTotalCells[2] = nTotalCells[3] = 0;
72 nCellsTrue[0] = nCellsTrue[1] = nCellsTrue[2] = 0;
73 nCellsGapStart[0] = nCellsGapStart[1] = nCellsGapStart[2] = 0;
74 nCellsGapSize[0] = nCellsGapSize[1] = nCellsGapSize[2] = 0;
75 invalid = true;
76 scaledPositionsValid = false;
77 wrappedPositionsValid = false;
78 supercell_counter = 0;
79 MakeTranslationTable();
80 DEBUGPRINT;
81 }
82
~NeighborCellLocator()83 NeighborCellLocator::~NeighborCellLocator()
84 {
85 DESTRUCTOR;
86 DEBUGPRINT;
87 for (int i = 0; i < nbCells_onthefly.size(); ++i)
88 delete nbCells_onthefly[i];
89 nbCells_onthefly.clear();
90 DEBUGPRINT;
91 AsapAtoms_DECREF(atoms);
92 DEBUGPRINT;
93 }
94
EnableNeighborsOfGhosts()95 void NeighborCellLocator::EnableNeighborsOfGhosts()
96 {
97 invalid = true;
98 includeghostneighbors = true;
99 }
100
MakeTranslationTable()101 void NeighborCellLocator::MakeTranslationTable()
102 {
103 DEBUGPRINT;
104 translationTable.resize(27);
105 for (int i = 0; i < 3; i++)
106 for (int j = 0; j < 3; j++)
107 for (int k = 0; k < 3; k++)
108 translationTable[i + 3*j + 9*k] = IVec(i == 2 ? -1 : i,
109 j == 2 ? -1 : j,
110 k == 2 ? -1 : k);
111 DEBUGPRINT;
112 }
113
MakeList()114 void NeighborCellLocator::MakeList()
115 {
116 RETURNIFASAPERROR;
117 DEBUGPRINT;
118 MEMORY;
119 USETIMER("NeighborCellLocator::MakeList");
120 if (verbose >= 1)
121 cerr << " NeighborCellLocator-Update ";
122 nAtoms = atoms->GetNumberOfAtoms();
123 nAllAtoms = nAtoms + atoms->GetNumberOfGhostAtoms();
124 if (includeghostneighbors)
125 nAtoms = nAllAtoms;
126 /// Update positions and scaledpositions
127 ScaleAndNormalizePositions();
128 ASSERT(scaledPositions.size() == nAllAtoms);
129 MEMORY;
130
131 const double *cellHeights = atoms->GetCellHeights();
132 const bool *atomsperiodic = atoms->GetBoundaryConditions();
133 for (int i = 0; i < 3; i++)
134 periodic[i] = atomsperiodic[i];
135 // Check the size. As this locator does not use minimum image
136 // convention, we can allow height down to rCut rather than the
137 // usual 2*rCut.
138 for (int i = 0; i < 3; i++)
139 if (periodic[i] && cellHeights[i] < rCut)
140 THROW( AsapError("NeighborCellLocator: The height of the cell (")
141 << cellHeights[i] << ") must be larger than " << rCut );
142
143 const vector<Vec> &UNUSED(positions) = GetWrappedPositions();
144 const vector<Vec> &scaledpositions = GetScaledPositions();
145
146 MEMORY;
147 // Find the system size in scaled space coordinates. Use this to
148 // determine the number of cells. If this processor only handles
149 // part of space in a periodic direction, we must search for a range
150 // of coordinates not present (it will be found on the outermost
151 // processors where ghosts from the other side of the simulation
152 // cell appear). This is used to massively cut down on the memory
153 // consumption.
154 IVec cpulayout = atoms->GetNumberOfCells();
155 bool cellLayoutChanged = false;
156 for (int i = 0; i < 3; i++)
157 {
158 bool lookforgap = periodic[i] && cpulayout[i] > 1;
159 bool gapvalid = false;
160 double gaplow;
161 double gaphigh;
162 if (periodic[i] && cpulayout[i] == 1)
163 {
164 // Periodic boundary conditions, use whole axis
165 size[i] = 1.0;
166 minimum[i] = 0.0;
167 }
168 else if (lookforgap)
169 {
170 // Periodic boundary conditions, use part of the axis (parallel sim).
171 double gap1start = 0.25;
172 double gap2start = 0.75;
173 double gap1low = -1.0;
174 double gap1hi = 2.0;
175 double gap2low = -1.0;
176 double gap2hi = 2.0;
177 bool gap1below = false;
178 bool gap1above = false;
179 bool gap2below = false;
180 bool gap2above = false;
181 double min, max;
182 min = max = scaledpositions[0][i];
183 for (int a = 1; a < nAllAtoms; a++)
184 {
185 double x = scaledpositions[a][i];
186 if (x > max)
187 max = x;
188 else if (x < min)
189 min = x;
190 if (x < gap1start && x > gap1low)
191 {
192 gap1low = x;
193 gap1below = true;
194 }
195 else if (x > gap1start && x < gap1hi)
196 {
197 gap1hi = x;
198 gap1above = true;
199 }
200 if (x < gap2start && x > gap2low)
201 {
202 gap2low = x;
203 gap2below = true;
204 }
205 else if (x > gap2start && x < gap2hi)
206 {
207 gap2hi = x;
208 gap2above = true;
209 }
210 }
211 minimum[i] = min;
212 size[i] = max - min;
213 #if 0
214 cerr << "GAP: " << gap1below << " " << gap1above << " "
215 << gap1low << " " << gap1hi << " --- "
216 << gap2below << " " << gap2above << " "
217 << gap2low << " " << gap2hi << endl;
218 #endif
219 gapvalid = false;
220 if (gap1below && gap1above)
221 {
222 gapvalid = true;
223 gaplow = gap1low;
224 gaphigh = gap1hi;
225 }
226 if (gap2below && gap2above &&
227 (!gapvalid || (gap2hi - gap2low > gap1hi - gap1low)))
228 {
229 gapvalid = true;
230 gaplow = gap2low;
231 gaphigh = gap2hi;
232 }
233 }
234 else
235 {
236 // Free boundary conditions.
237 // Could be merged with section above, but if statements in
238 // loops hurt performance.
239 double min, max;
240 min = max = scaledpositions[0][i];
241 for (int a = 1; a < nAllAtoms; a++)
242 {
243 double x = scaledpositions[a][i];
244 if (x > max)
245 max = x;
246 else if (x < min)
247 min = x;
248 }
249 minimum[i] = min;
250 size[i] = max - min;
251 }
252 minimum[i] -= 1e-4;
253 size[i] += 2e-4;
254 int nc = int(size[i] * cellHeights[i] / minboxsize);
255 int gapstart, gapsize;
256 if (gapvalid)
257 {
258 gapstart = int(ceil(nc * (gaplow - minimum[i]) / size[i])) + 2;
259 int gapend = int(floor(nc * (gaphigh - minimum[i]) / size[i])) - 2;
260 gapsize = gapend - gapstart;
261 if (gapsize < 3)
262 {
263 gapsize = 0;
264 gapstart = nc + 1;
265 }
266 }
267 else
268 {
269 gapsize = 0;
270 gapstart = nc + 1;
271 }
272 if (nc == 0)
273 {
274 nc = 1;
275 minimum[i] -= 0.5 * (minboxsize / cellHeights[i] - size[i]);
276 size[i] = minboxsize / cellHeights[i];
277 }
278 if (nCells[i] != nc - gapsize || nCellsTrue[i] != nc
279 || nCellsGapStart[i] != gapstart || nCellsGapSize[i] != gapsize)
280 {
281 nCells[i] = nc - gapsize;
282 nCellsTrue[i] = nc;
283 nCellsGapStart[i] = gapstart;
284 nCellsGapSize[i] = gapsize;
285 cellLayoutChanged = true;
286 #if 0
287 cerr << "Axis " << i << "(pbc:" << periodic[i] << ") nCells="
288 << nCells[i] << " (gap: " << gapstart << ", " << gapsize
289 << ")" << " minimum: " << minimum[i] << " size: "
290 << size[i] << endl;
291 #endif
292 }
293 }
294 nTotalCells[0] = 1;
295 for (int i = 0; i < 3; i++)
296 nTotalCells[i + 1] = nTotalCells[i] * nCells[i];
297
298 if (verbose >= 2)
299 cerr << endl << "nNeighborCellLocator: " << nAtoms << " atoms, "
300 << nCells[0] << "*" << nCells[1] << "*" << nCells[2] << " = "
301 << nTotalCells[3] << " cells, "
302 << nAtoms / (double) nTotalCells[3] << " atoms/cell" << endl;
303 // Now check if the boundary conditions have changed
304 for (int i = 0; i < 3; i++)
305 {
306 if (periodic[i] != oldperiodic[i])
307 cellLayoutChanged = true;
308 oldperiodic[i] = periodic[i];
309 }
310
311 MEMORY;
312 // Now reallocate the cell list, if its size has changed
313 if (cellLayoutChanged || (cells.size() != nTotalCells[3]))
314 {
315 // Clear the old cells, create new ones.
316 cells.clear();
317 cells.resize(nTotalCells[3]);
318 MEMORY;
319 MakeNeighboringCellLists();
320 }
321 else
322 {
323 // Empty the old cells, but keep them to minimize reallocations
324 ASSERT(cells.size() == nTotalCells[3]);
325 for (int i = 0; i < cells.size(); i++)
326 cells[i].clear();
327 }
328 MEMORY;
329 cellIndices.resize(nAllAtoms);
330 MEMORY;
331
332 // Now we are ready to put the atoms into the cells!
333 for (int a = 0; a < nAllAtoms; a++)
334 {
335 int index = 0;
336 for (int i = 0; i < 3; i++)
337 {
338 int k =
339 int((scaledpositions[a][i] - minimum[i]) / size[i] * nCellsTrue[i]);
340 ASSERT(k >= 0);
341 if (k > nCellsGapStart[i])
342 {
343 ASSERT (k > nCellsGapStart[i] + nCellsGapSize[i]);
344 k -= nCellsGapSize[i];
345 }
346 if (k == nCells[i])
347 k--;
348 ASSERT(k < nCells[i]);
349 index += nTotalCells[i] * k;
350 }
351 cells[index].push_back(a);
352 cellIndices[a] = index;
353 }
354
355 // update reference positions for neighbor list:
356 atoms->GetPositions(referencePositions);
357 MEMORY;
358
359 // Find the largest cell size, and use it to calculate the max length
360 maxLength = 0;
361 for (int i = 0; i < cells.size(); i++)
362 {
363 int x = cells[i].size();
364 if (x > maxLength)
365 maxLength = x;
366 }
367 #ifdef PATTERN5
368 maxLength *= 5*5*5;
369 #else
370 maxLength *= 3*3*3;
371 #endif
372 invalid = false;
373 DEBUGPRINT;
374 MEMORY;
375 }
376
MakeNeighboringCellLists()377 void NeighborCellLocator::MakeNeighboringCellLists()
378 {
379 DEBUGPRINT;
380
381 #ifdef PATTERN5
382 int dx = 2;
383 #else
384 int dx = 1;
385 #endif
386 neighborCellOffsets.clear();
387 for (int i = -dx; i <= dx; i++)
388 for (int j = -dx; j <= dx; j++)
389 for (int k = -dx; k <= dx; k++)
390 neighborCellOffsets.push_back(IVec(i,j,k));
391
392 nbCells_inside.clear();
393 nbCells_left.clear();
394 nbCells_right.clear();
395 nbCells_top.clear();
396 nbCells_bottom.clear();
397 nbCells_front.clear();
398 nbCells_back.clear();
399 for(vector<IVec>::const_iterator nb = neighborCellOffsets.begin();
400 nb != neighborCellOffsets.end(); ++nb)
401 {
402 IVec idx = *nb;
403 IVec magic(1,3,9);
404 // inside
405 pair<int, translationsidx_t> data0;
406 pair<int, translationsidx_t> data;
407 data0.first = idx[0] + idx[1] * nTotalCells[1]
408 + idx[2] * nTotalCells[2];
409 data0.second = 0;
410 nbCells_inside.push_back(data0);
411 // left side
412 data = data0;
413 if (idx[0] != -1)
414 nbCells_left.push_back(data);
415 else if (periodic[0])
416 {
417 data.first += nCells[0];
418 data.second = 1 * magic[0];
419 nbCells_left.push_back(data);
420 }
421 // right side
422 data = data0;
423 if (idx[0] != 1)
424 nbCells_right.push_back(data);
425 else if (periodic[0])
426 {
427 data.first -= nCells[0];
428 data.second = 2 * magic[0];
429 nbCells_right.push_back(data);
430 }
431 // bottom
432 data = data0;
433 if (idx[1] != -1)
434 nbCells_bottom.push_back(data);
435 else if (periodic[1])
436 {
437 data.first += nTotalCells[2];
438 data.second = 1 * magic[1];
439 nbCells_bottom.push_back(data);
440 }
441 // top
442 data = data0;
443 if (idx[1] != 1)
444 nbCells_top.push_back(data);
445 else if (periodic[1])
446 {
447 data.first -= nTotalCells[2];
448 data.second = 2 * magic[1];
449 nbCells_top.push_back(data);
450 }
451 // front
452 data = data0;
453 if (idx[2] != -1)
454 nbCells_front.push_back(data);
455 else if (periodic[2])
456 {
457 data.first += nTotalCells[3];
458 data.second = 1 * magic[2];
459 nbCells_front.push_back(data);
460 }
461 // back
462 data = data0;
463 if (idx[2] != 1)
464 nbCells_back.push_back(data);
465 else if (periodic[2])
466 {
467 data.first -= nTotalCells[3];
468 data.second = 2 * magic[2];
469 nbCells_back.push_back(data);
470 }
471 }
472 DEBUGPRINT;
473 //std::cerr << "************** cells.size(): " << cells.size() << std::endl;
474 // Make the corner cases (literally!), and build the map
475 nbCells_all.clear();
476 for (int i = 0; i < nbCells_onthefly.size(); ++i)
477 delete nbCells_onthefly[i];
478 nbCells_onthefly.clear();
479 for (int i = 0; i < cells.size(); i++)
480 makeNbCells(i);
481 ASSERT(nbCells_all.size() == cells.size());
482 }
483
GetNeighbors(int a1,int * neighbors,Vec * diffs,double * diffs2,int & size,double r) const484 int NeighborCellLocator::GetNeighbors(int a1, int *neighbors, Vec *diffs,
485 double *diffs2, int& size,
486 double r) const
487 {
488 return CommonGetNeighbors(a1, neighbors, diffs, diffs2, size, r, false);
489 }
490
GetNeighbors(int a1,vector<int> & neighbors) const491 void NeighborCellLocator::GetNeighbors(int a1, vector<int> &neighbors) const
492 {
493 CommonGetNeighbors(a1, neighbors, false);
494 }
495
GetFullNeighbors(int a1,int * neighbors,Vec * diffs,double * diffs2,int & size,double r) const496 int NeighborCellLocator::GetFullNeighbors(int a1, int *neighbors, Vec *diffs,
497 double *diffs2, int& size,
498 double r) const
499 {
500 return CommonGetNeighbors(a1, neighbors, diffs, diffs2, size, r, true);
501 }
502
503
GetFullNeighbors(int a1,vector<int> & neighbors) const504 void NeighborCellLocator::GetFullNeighbors(int a1, vector<int> &neighbors) const
505 {
506 CommonGetNeighbors(a1, neighbors, true);
507 }
508
makeNbCells(int thiscell)509 void NeighborCellLocator::makeNbCells(int thiscell)
510 {
511 IVec cellidx(thiscell % nTotalCells[1],
512 (thiscell % nTotalCells[2]) / nTotalCells[1],
513 thiscell / nTotalCells[2]);
514 ASSERT(thiscell == (cellidx[0] * nTotalCells[0] +
515 cellidx[1] * nTotalCells[1] +
516 cellidx[2] * nTotalCells[2]));
517 int celltype = (cellidx[0] == 0) + 2 * (cellidx[0] == nCells[0]-1)
518 + 4 * (cellidx[1] == 0) + 8 * (cellidx[1] == nCells[1]-1)
519 + 16 * (cellidx[2] == 0) + 32 * (cellidx[2] == nCells[2]-1);
520
521 if (celltype == 0)
522 nbCells_all[thiscell] = &nbCells_inside;
523 else if (celltype == 1)
524 nbCells_all[thiscell] = &nbCells_left;
525 else if (celltype == 2)
526 nbCells_all[thiscell] = &nbCells_right;
527 else if (celltype == 4)
528 nbCells_all[thiscell] = &nbCells_bottom;
529 else if (celltype == 8)
530 nbCells_all[thiscell] = &nbCells_top;
531 else if (celltype == 16)
532 nbCells_all[thiscell] = &nbCells_front;
533 else if (celltype == 32)
534 nbCells_all[thiscell] = &nbCells_back;
535 else
536 {
537 // Edge or corner cell - make it on the fly.
538 nbcell_t *NbCells = new nbcell_t;
539 nbCells_onthefly.push_back(NbCells);
540 NbCells->clear();
541 nbCells_all[thiscell] = NbCells;
542
543 for (vector<IVec>::const_iterator i = neighborCellOffsets.begin();
544 i != neighborCellOffsets.end(); ++i)
545 {
546 IVec othercell = cellidx + *i;
547 translationsidx_t xlat = 0;
548 IVec xlatvec(0,0,0);
549 IVec magic(1,3,9);
550 bool ok = true;
551 for (int j = 0; j < 3; j++)
552 {
553 if (othercell[j] < 0)
554 {
555 if (!periodic[j])
556 {
557 ok = false;
558 break;
559 }
560 othercell[j] += nCells[j];
561 //std::cerr << "*** Magic j=" << j << ": " << magic[j] << std::endl;
562 xlat += magic[j] * 1;
563 xlatvec[j] = 1;
564 }
565 else if (othercell[j] >= nCells[j])
566 {
567 if (!periodic[j])
568 {
569 ok = false;
570 break;
571 }
572 othercell[j] -= nCells[j];
573 //std::cerr << "*** Magic j=" << j << ": " << magic[j] << std::endl;
574 xlat += magic[j] * 2;
575 xlatvec[j] = -1;
576 }
577 //std::cerr << xlat << " ";
578 }
579 if (ok) {
580 //std::cerr << xlatvec << " ==?== " << xlat << ": " << translationTable[xlat] << std::endl;
581 ASSERT(xlatvec == translationTable.at(xlat));
582 pair<int, translationsidx_t> data;
583 data.first = (othercell[0] * nTotalCells[0] +
584 othercell[1] * nTotalCells[1] +
585 othercell[2] * nTotalCells[2]) - thiscell;
586 data.second = xlat;
587 NbCells->push_back(data);
588 }
589 }
590 }
591 }
592
CommonGetNeighbors(int a1,int * RESTRICT neighbors,Vec * RESTRICT diffs,double * RESTRICT diffs2,int & size,double r,bool wantfull) const593 int NeighborCellLocator::CommonGetNeighbors(int a1, int *RESTRICT neighbors, Vec *RESTRICT diffs,
594 double *RESTRICT diffs2, int& size,
595 double r, bool wantfull) const
596 {
597 if (invalid)
598 throw AsapError("NeighborCellLocator has been invalidated, possibly by another NeighborList using the same atoms.");
599
600 const vector<Vec> &RESTRICT positions = GetWrappedPositions();
601 // Need to use GET_CELL instead of GetCell as the atoms are not open
602 // when called from the Python interface.
603 const Vec *RESTRICT superCell = atoms->GET_CELL();
604
605 // Find the cell of this atom
606 int thiscell = cellIndices[a1];
607 double rC2;
608 if (r > 0.0)
609 rC2 = r*r;
610 else
611 rC2 = rCut2;
612
613 int nNeighbors = 0;
614
615 // Ghost atoms have no neighbors - but they are neighbors to real atoms.
616 if (a1 < nAtoms)
617 {
618 const nbcell_t *neighborCells = nbCells_all.at(thiscell);
619
620 // Loop over all neighboring cells, including this one
621 for (vector< pair<int, translationsidx_t> >::const_iterator i =
622 neighborCells->begin(); i < neighborCells->end(); ++i)
623 {
624 int othcell = i->first + thiscell;
625 IVec celltranslation = translationTable[i->second];
626 Vec pos1 = positions[a1] + celltranslation[0] * superCell[0] +
627 celltranslation[1] * superCell[1] + celltranslation[2] * superCell[2];
628 // Loop over all atoms in the cell.
629 vector<int>::const_iterator aend = cells[othcell].end();
630
631 for (vector<int>::const_iterator a2 = cells[othcell].begin();
632 a2 < aend; ++a2)
633 {
634 diffs[nNeighbors] = positions[*a2] - pos1;
635 diffs2[nNeighbors] = Length2(diffs[nNeighbors]);
636 neighbors[nNeighbors] = *a2;
637 nNeighbors++;
638 }
639 }
640 }
641 // Remove neighbors too far away
642 int j = 0;
643 for (int i = 0; i < nNeighbors; i++)
644 {
645 int a2 = neighbors[i];
646 if (i != j)
647 {
648 diffs[j] = diffs[i];
649 diffs2[j] = diffs2[i];
650 neighbors[j] = a2;
651 }
652 if ((diffs2[i] < rC2) && ((a2 > a1) || (wantfull && (a2 != a1))))
653 j++;
654 }
655 nNeighbors = j;
656 size -= nNeighbors;
657 ASSERT(size >= 0);
658 return nNeighbors;
659 }
660
CommonGetNeighbors(int a1,vector<int> & neighbors,bool wantfull) const661 void NeighborCellLocator::CommonGetNeighbors(int a1, vector<int> &neighbors,
662 bool wantfull) const
663 {
664 if (invalid)
665 throw AsapError("NeighborCellLocator has been invalidated, possibly by another NeighborList using the same atoms.");
666
667 const vector<Vec> &positions_x = GetWrappedPositions();
668 const Vec *RESTRICT superCell = atoms->GET_CELL();
669
670 // Find the cell of this atom
671 int thiscell = cellIndices[a1];
672 neighbors.resize(maxLength);
673 // vector<double> d2(maxLength);
674 double *RESTRICT d2 = new double[maxLength];
675 int n = 0;
676
677 const Vec *RESTRICT positions = &positions_x[0];
678 int *RESTRICT neighbors_p = &neighbors[0];
679
680 const IVec *translationTable = &(this->translationTable)[0];
681 // Ghost atoms have no neighbors - but they are neighbors to real atoms.
682 if (a1 < nAtoms)
683 {
684 const nbcell_t *neighborCells = nbCells_all.at(thiscell);
685
686 // Loop over all neighboring cells, including this one
687 vector< pair<int, translationsidx_t> >::const_iterator end = neighborCells->end();
688 for (vector< pair<int, translationsidx_t> >::const_iterator i =
689 neighborCells->begin(); i < end; ++i)
690 {
691 int othcell = i->first + thiscell;
692 IVec celltranslation = translationTable[i->second];
693 Vec pos1 = positions[a1] + celltranslation[0] * superCell[0] +
694 celltranslation[1] * superCell[1] + celltranslation[2] * superCell[2];
695 // Loop over all atoms in the cell.
696 vector<int>::const_iterator aend = cells[othcell].end();
697 for (vector<int>::const_iterator a2 = cells[othcell].begin();
698 a2 < aend; ++a2)
699 {
700 Vec diff = positions[*a2] - pos1;
701 d2[n] = Length2(diff);
702 neighbors_p[n] = *a2;
703 n++;
704 }
705 }
706 }
707 // Remove neighbors too far away
708 int j = 0;
709 for (int i = 0; i < n; i++)
710 {
711 int a2 = neighbors[i];
712 if (i != j)
713 neighbors[j] = a2;
714 if ((d2[i] < rCut2) && ((a2 > a1) || (wantfull && (a2 != a1))))
715 j++;
716 }
717 neighbors.resize(j);
718 delete[] d2;
719 }
720
GetListAndTranslations(int a1,vector<neighboritem_t> & neighbors) const721 int NeighborCellLocator::GetListAndTranslations(int a1,
722 vector<neighboritem_t> &neighbors) const
723 {
724 if (invalid)
725 THROW( AsapError("NeighborCellLocator has been invalidated, possibly by another NeighborList using the same atoms.") );
726 RETURNIFASAPERROR2(0);
727
728 const vector<Vec> &positions = GetWrappedPositions();
729 const Vec *superCell = atoms->GetCell();
730
731 // Find the cell of this atom
732 int thiscell = cellIndices[a1];
733 double rC2 = rCut2;
734
735 neighbors.clear();
736 // Ghost atoms have no neighbors - but they are neighbors to real atoms.
737 if (a1 < nAtoms)
738 {
739 const nbcell_t *neighborCells = nbCells_all.at(thiscell);
740
741 // Loop over all neighboring cells, including this one
742 for (vector< pair<int, translationsidx_t> >::const_iterator i =
743 neighborCells->begin(); i < neighborCells->end(); ++i)
744 {
745 int othcell = i->first + thiscell;
746 IVec celltranslation = translationTable[i->second];
747 Vec pos1 = positions[a1] + celltranslation[0] * superCell[0] +
748 celltranslation[1] * superCell[1] + celltranslation[2] * superCell[2];
749 // Loop over all atoms in the cell.
750 vector<int>::const_iterator aend = cells[othcell].end();
751 for (vector<int>::const_iterator a2 = cells[othcell].begin();
752 a2 < aend; ++a2)
753 {
754 double l2;
755 if ((*a2 > a1) && ((l2 = Length2(positions[*a2] - pos1)) < rC2))
756 {
757 #ifdef SLOWASSERT
758 if (l2 < 1e-6)
759 THROW( AsapError("XX Collision between atoms ") << a1 << " and " << *a2 );
760 #endif
761 //neighboritem_t data(*a2, i->second);
762 neighboritem_t NEIGHBORITEM_PACK(data, *a2, i->second);
763 neighbors.push_back(data);
764 }
765 }
766 }
767 }
768 return neighbors.size();
769 }
770
GetComplementaryListAndTranslations(int a1,vector<neighboritem_t> & neighbors) const771 int NeighborCellLocator::GetComplementaryListAndTranslations(int a1,
772 vector<neighboritem_t> &neighbors) const
773 {
774 RETURNIFASAPERROR2(0);
775 if (invalid)
776 THROW( AsapError("NeighborCellLocator has been invalidated, possibly by another NeighborList using the same atoms.") );
777
778 const vector<Vec> &positions = GetWrappedPositions();
779 const Vec *superCell = atoms->GetCell();
780
781 // Find the cell of this atom
782 int thiscell = cellIndices[a1];
783 double rC2 = rCut2;
784
785 neighbors.clear();
786
787 // Ghost atoms have no neighbors - but they are neighbors to real atoms.
788 if (a1 < nAtoms)
789 {
790 const nbcell_t *neighborCells = nbCells_all.at(thiscell);
791
792 // Loop over all neighboring cells, including this one
793 for (vector< pair<int, translationsidx_t> >::const_iterator i =
794 neighborCells->begin(); i < neighborCells->end(); ++i)
795 {
796 int othcell = i->first + thiscell;
797 IVec celltranslation = translationTable[i->second];
798 Vec pos1 = positions[a1] + celltranslation[0] * superCell[0] +
799 celltranslation[1] * superCell[1] + celltranslation[2] * superCell[2];
800 // Loop over all atoms in the cell.
801 vector<int>::const_iterator aend = cells[othcell].end();
802 for (vector<int>::const_iterator a2 = cells[othcell].begin();
803 a2 < aend; ++a2)
804 {
805 if ((*a2 < a1) && (Length2(positions[*a2] - pos1) < rC2))
806 {
807 //neighboritem_t data(*a2, i->second);
808 neighboritem_t NEIGHBORITEM_PACK(data, *a2, i->second);
809 neighbors.push_back(data);
810 }
811 }
812 }
813 }
814 return neighbors.size();
815 }
816
get_drift() const817 double NeighborCellLocator::get_drift() const
818 {
819 // Find the max allowed drift of an atom.
820 const double *superCellHeight = atoms->GetCellHeights();
821 double height = superCellHeight[0]/nCellsTrue[0];
822 for (int i = 1; i < 3; i++)
823 {
824 double h = superCellHeight[i]/nCellsTrue[i];
825 if (h < height)
826 height = h;
827 }
828 #ifdef PATTERN5
829 return 0.5 * (2 * height - rCut);
830 #else
831 return 0.5 * (height - rCut);
832 #endif
833 }
834
CheckNeighborList()835 bool NeighborCellLocator::CheckNeighborList()
836 {
837 DEBUGPRINT;
838 USETIMER("NeighborCellLocator::CheckNeighborList");
839
840 const bool *newpbc = atoms->GetBoundaryConditions();
841 if (nAtoms != atoms->GetNumberOfAtoms() || oldperiodic[0] != newpbc[0]
842 || oldperiodic[1] != newpbc[1] || oldperiodic[2] != newpbc[2])
843 invalid = true;
844
845 if (invalid)
846 return true; // An invalid neighbor list always need an update!
847
848 RenormalizePositions();
849
850 double drift = get_drift();
851 double drift2 = drift * drift;
852 bool updateRequired = invalid;
853 const Vec *positions = atoms->GetPositions();
854
855 if (!updateRequired)
856 for (int n = 0; n < nAtoms; n++)
857 if (Length2(positions[n] - referencePositions[n]) > drift2)
858 {
859 updateRequired = true;
860 //break; // prevents vectorization
861 }
862 // If we are NOT going to update the list, we need to update wrapped positions
863 DEBUGPRINT;
864 return updateRequired;
865 }
866
UpdateNeighborList()867 void NeighborCellLocator::UpdateNeighborList()
868 {
869 RETURNIFASAPERROR;
870 DEBUGPRINT;
871 USETIMER("NeighborCellLocator::UpdateNeighborList");
872 if (invalid && verbose)
873 cerr << "NeighborCellLocator::UpdateNeighborList: NBList has been marked invalid." << endl;
874
875 MakeList();
876 DEBUGPRINT;
877 }
878
CheckAndUpdateNeighborList()879 bool NeighborCellLocator::CheckAndUpdateNeighborList()
880 {
881 DEBUGPRINT;
882 bool update = CheckNeighborList();
883 if (update)
884 UpdateNeighborList();
885 DEBUGPRINT;
886 return update;
887 }
888
CheckAndUpdateNeighborList(PyObject * atoms_obj)889 bool NeighborCellLocator::CheckAndUpdateNeighborList(PyObject *atoms_obj)
890 {
891 atoms->Begin(atoms_obj);
892 CHECKNOASAPERROR;
893 bool res = CheckAndUpdateNeighborList();
894 PROPAGATEASAPERROR;
895 atoms->End();
896 return res;
897 }
898
GetTranslationTable(vector<IVec> & table) const899 void NeighborCellLocator::GetTranslationTable(vector<IVec> &table) const
900 {
901 DEBUGPRINT;
902 table.clear();
903 table.insert(table.begin(), translationTable.begin(), translationTable.end());
904 DEBUGPRINT;
905 }
906
907 #if 0
908 void NeighborCellLocator::RenormalizeReferencePositions(const vector<Vec>
909 &oldpos)
910 {
911 if (!invalid)
912 {
913 int n = atoms->GetNumberOfAtoms();
914 ASSERT(referencePositions.size() == n);
915 const Vec *pos = atoms->GetPositions();
916 for (int i = 0; i < n; ++i)
917 referencePositions[i] += pos[i] - oldpos[i];
918 }
919 }
920 #endif
921
GetScaledPositions() const922 const vector<Vec> &NeighborCellLocator::GetScaledPositions() const
923 {
924 DEBUGPRINT;
925 ASSERT(scaledPositionsValid);
926 return scaledPositions;
927 }
928
ScaleAndNormalizePositions()929 void NeighborCellLocator::ScaleAndNormalizePositions()
930 {
931 DEBUGPRINT;
932 atoms->GetScaledPositions(scaledPositions, true); // Get also ghosts
933 ASSERT(scaledPositions.size() == nAllAtoms);
934 const bool *pbc = atoms->GetBoundaryConditions();
935 // Special-case fully periodic and fully free boundaries
936 if (pbc[0] && pbc[1] && pbc[2])
937 {
938 // Fully periodic boundaries
939 int spsz = scaledPositions.size();
940 if (wrappedPositions.capacity() < spsz)
941 wrappedPositions.reserve(spsz + spsz / 25);
942 wrappedPositions.resize(spsz);
943 if (offsetPositions.capacity() < spsz)
944 offsetPositions.reserve(spsz + spsz / 25);
945 offsetPositions.resize(scaledPositions.size());
946 scaledOffsetPositions.clear(); // Not used
947 const Vec *pos = atoms->GetPositions();
948 const Vec *cell = atoms->GetCell();
949 int n = scaledPositions.size();
950 Vec *RESTRICT scaledPositions = &(this->scaledPositions)[0]; // Helps vectorization
951 Vec *RESTRICT wrappedPositions = &(this->wrappedPositions)[0];
952 Vec *RESTRICT offsetPositions = &(this->offsetPositions)[0];
953 for (int i = 0; i < n; i++)
954 {
955 scaledPositions[i].x -= floor(scaledPositions[i].x);
956 scaledPositions[i].y -= floor(scaledPositions[i].y);
957 scaledPositions[i].z -= floor(scaledPositions[i].z);
958 wrappedPositions[i] = cell[0] * scaledPositions[i].x
959 + cell[1] * scaledPositions[i].y
960 + cell[2] * scaledPositions[i].z;
961 offsetPositions[i] = wrappedPositions[i] - pos[i];
962 }
963 scaledPositionsValid = wrappedPositionsValid = true;
964 }
965 else if (!pbc[0] && !pbc[1] && !pbc[2])
966 {
967 // Fully free boundaries: wrappedPositions are the positions
968 atoms->GetPositions(wrappedPositions, true);
969 offsetPositions.clear();
970 scaledOffsetPositions.clear();
971 scaledPositionsValid = wrappedPositionsValid = true;
972 }
973 else
974 {
975 // Mixed boundary conditions
976 double xpbc0 = (double) pbc[0];
977 double xpbc1 = (double) pbc[1];
978 double xpbc2 = (double) pbc[2];
979 int spsz = scaledPositions.size();
980 if (wrappedPositions.capacity() < spsz)
981 wrappedPositions.reserve(spsz + spsz / 25);
982 wrappedPositions.resize(spsz);
983 if (scaledOffsetPositions.capacity() < spsz)
984 scaledOffsetPositions.reserve(spsz + spsz / 25);
985 scaledOffsetPositions.resize(scaledPositions.size());
986 offsetPositions.clear(); // Not used
987 const Vec *cell = atoms->GetCell();
988 int n = scaledPositions.size();
989 Vec *RESTRICT scaledPositions = &(this->scaledPositions)[0]; // Helps vectorization
990 Vec *RESTRICT wrappedPositions = &(this->wrappedPositions)[0];
991 Vec *RESTRICT scaledOffsetPositions = &(this->scaledOffsetPositions)[0];
992 for (int i = 0; i < n; i++)
993 {
994 scaledOffsetPositions[i].x = -floor(scaledPositions[i].x) * xpbc0;
995 scaledPositions[i].x += scaledOffsetPositions[i].x;
996 scaledOffsetPositions[i].y = -floor(scaledPositions[i].y) * xpbc1;
997 scaledPositions[i].y += scaledOffsetPositions[i].y;
998 scaledOffsetPositions[i].z = -floor(scaledPositions[i].z) * xpbc2;
999 scaledPositions[i].z += scaledOffsetPositions[i].z;
1000 wrappedPositions[i] = cell[0] * scaledPositions[i].x
1001 + cell[1] * scaledPositions[i].y
1002 + cell[2] * scaledPositions[i].z;
1003 }
1004 scaledPositionsValid = wrappedPositionsValid = true;
1005 }
1006 // We need the inverse cell later.
1007 memcpy(old_inverse_cell, atoms->GetInverseCell(), 3*sizeof(Vec));
1008 supercell_counter = atoms->GetCellCounter();
1009 DEBUGPRINT;
1010 }
1011
ScaleAndNormalizePositions(const set<int> & modified,vector<Vec> & scaledpos)1012 void NeighborCellLocator::ScaleAndNormalizePositions(const set<int> &modified,
1013 vector<Vec> &scaledpos)
1014 {
1015 DEBUGPRINT;
1016 ASSERT(modified.size() == scaledpos.size());
1017 atoms->GetScaledPositions(scaledpos, modified);
1018 const bool *pbc = atoms->GetBoundaryConditions();
1019 // Special-case fully periodic and fully free boundaries, as they
1020 // are special-cased in the normal version.
1021 if (pbc[0] && pbc[1] && pbc[2])
1022 {
1023 const Vec *pos = atoms->GetPositions();
1024 const Vec *cell = atoms->GetCell();
1025 vector<Vec>::iterator sp = scaledpos.begin();
1026 for (set<int>::const_iterator ii = modified.begin();
1027 ii != modified.end(); ++ii)
1028 {
1029 int i = *ii;
1030 scaledPositions[i] = *sp;
1031 scaledPositions[i][0] -= floor(scaledPositions[i][0]);
1032 scaledPositions[i][1] -= floor(scaledPositions[i][1]);
1033 scaledPositions[i][2] -= floor(scaledPositions[i][2]);
1034 *(sp++) = scaledPositions[i];
1035 wrappedPositions[i] = cell[0] * scaledPositions[i][0]
1036 + cell[1] * scaledPositions[i][1]
1037 + cell[2] * scaledPositions[i][2];
1038 offsetPositions[i] = wrappedPositions[i] - pos[i];
1039 }
1040 scaledPositionsValid = wrappedPositionsValid = true;
1041 }
1042 else if (!pbc[0] && !pbc[1] && !pbc[2])
1043 {
1044 // Fully free boundaries: wrappedPositions are the positions
1045 const Vec *r = atoms->GetPositions();
1046 vector<Vec>::iterator sp = scaledpos.begin();
1047 for (set<int>::const_iterator ii = modified.begin();
1048 ii != modified.end(); ++ii)
1049 {
1050 int i = *ii;
1051 scaledPositions[i] = *(sp++);
1052 wrappedPositions[i] = r[i];
1053 }
1054 scaledPositionsValid = wrappedPositionsValid = true;
1055 }
1056 else
1057 {
1058 // Mixed boundary conditions
1059 int xpbc0 = (int) pbc[0];
1060 int xpbc1 = (int) pbc[1];
1061 int xpbc2 = (int) pbc[2];
1062 const Vec *cell = atoms->GetCell();
1063 vector<Vec>::iterator sp = scaledpos.begin();
1064 for (set<int>::const_iterator ii = modified.begin();
1065 ii != modified.end(); ++ii)
1066 {
1067 int i = *ii;
1068 scaledPositions[i] = *sp;
1069 scaledOffsetPositions[i][0] = -floor(scaledPositions[i][0]) * xpbc0;
1070 scaledPositions[i][0] += scaledOffsetPositions[i][0];
1071 scaledOffsetPositions[i][1] = -floor(scaledPositions[i][1]) * xpbc1;
1072 scaledPositions[i][1] += scaledOffsetPositions[i][1];
1073 scaledOffsetPositions[i][2] = -floor(scaledPositions[i][2]) * xpbc2;
1074 scaledPositions[i][2] += scaledOffsetPositions[i][2];
1075 *(sp++) = scaledPositions[i];
1076 wrappedPositions[i] = cell[0] * scaledPositions[i][0]
1077 + cell[1] * scaledPositions[i][1]
1078 + cell[2] * scaledPositions[i][2];
1079 }
1080 scaledPositionsValid = wrappedPositionsValid = true;
1081 }
1082 }
1083
RenormalizePositions()1084 void NeighborCellLocator::RenormalizePositions()
1085 {
1086 DEBUGPRINT;
1087 scaledPositionsValid = false;
1088 int nAtoms = this->nAtoms; // Helps vectorization.
1089 int nAllAtoms = this->nAllAtoms;
1090 const bool *pbc = atoms->GetBoundaryConditions();
1091 if (pbc[0] && pbc[1] && pbc[2])
1092 {
1093 // Full periodic boundaries
1094 if (atoms->GetCellCounter() != supercell_counter)
1095 {
1096 // The unit cell has changed, and we have full period
1097 // boundary conditions. Assume that atoms have moved along
1098 // the rescaling.
1099 Vec transformation[3];
1100 matrixMultiply3x3(transformation, old_inverse_cell, atoms->GetCell());
1101 memcpy(old_inverse_cell, atoms->GetInverseCell(), 3*sizeof(Vec));
1102 supercell_counter = atoms->GetCellCounter();
1103 //cerr << endl << __FUNCTION__ << ": " << referencePositions.size()
1104 // << " " << offsetPositions.size() << endl;
1105 ASSERT(referencePositions.size() == nAtoms);
1106 ASSERT(offsetPositions.size() == nAllAtoms);
1107 vector<Vec>::iterator rp = referencePositions.begin();
1108 vector<Vec>::iterator op = offsetPositions.begin();
1109 for(int i = 0; i < nAtoms; ++i, ++rp, ++op)
1110 {
1111 *op = transformation[0] * (*op)[0] + transformation[1] * (*op)[1]
1112 +transformation[2] * (*op)[2];
1113 *rp = transformation[0] * (*rp)[0] + transformation[1] * (*rp)[1]
1114 +transformation[2] * (*rp)[2];
1115 }
1116 ASSERT(rp == referencePositions.end());
1117 for(int i = nAtoms; i < nAllAtoms; ++i, ++op)
1118 *op = transformation[0] * (*op)[0] + transformation[1] * (*op)[1]
1119 +transformation[2] * (*op)[2];
1120 ASSERT(op == offsetPositions.end());
1121 }
1122 // We now need to rewrap
1123 ASSERT(wrappedPositions.size() == nAllAtoms);
1124 const Vec *pos = atoms->GetPositions();
1125 vector<Vec>::const_iterator op = offsetPositions.begin();
1126 vector<Vec>::iterator end = wrappedPositions.end();
1127 for (vector<Vec>::iterator wp = wrappedPositions.begin(); wp != end;
1128 ++wp, ++pos, ++op)
1129 *wp = *pos + *op;
1130 wrappedPositionsValid = true;
1131 }
1132 else
1133 {
1134 // Free or mixed boundary conditions
1135 if (atoms->GetCellCounter() != supercell_counter)
1136 {
1137 // The unit cell has changed. Assume that atoms have moved
1138 // along the rescaling.
1139 Vec transformation[3];
1140 matrixMultiply3x3(transformation, old_inverse_cell, atoms->GetCell());
1141 memcpy(old_inverse_cell, atoms->GetInverseCell(), 3*sizeof(Vec));
1142 supercell_counter = atoms->GetCellCounter();
1143 vector<Vec>::iterator end = referencePositions.end();
1144 for (vector<Vec>::iterator rp = referencePositions.begin();
1145 rp < end; ++rp)
1146 {
1147 *rp = transformation[0] * (*rp)[0] + transformation[1] * (*rp)[1]
1148 +transformation[2] * (*rp)[2];
1149 }
1150 }
1151 if (!pbc[0] && !pbc[1] && !pbc[2])
1152 {
1153 // Fully free boundary conditions.
1154 atoms->GetPositions(wrappedPositions, true); // Include ghosts
1155 wrappedPositionsValid = true;
1156 }
1157 else
1158 {
1159 // Mixed boundary conditions
1160 const Vec *cell = atoms->GetCell();
1161 atoms->GetScaledPositions(scaledPositions, true);
1162 ASSERT(scaledPositions.size() == scaledOffsetPositions.size());
1163 ASSERT(wrappedPositions.size() == scaledOffsetPositions.size());
1164 vector<Vec>::const_iterator sop = scaledOffsetPositions.begin();
1165 vector<Vec>::iterator sp = scaledPositions.begin();
1166 vector<Vec>::iterator end = wrappedPositions.end();
1167 for (vector<Vec>::iterator wp = wrappedPositions.begin();
1168 wp != end; ++wp, ++sop, ++sp)
1169 {
1170 *sp += *sop;
1171 *wp = cell[0] * (*sp)[0] + cell[1] * (*sp)[1]
1172 + cell[2] * (*sp)[2];
1173 }
1174 wrappedPositionsValid = true;
1175 }
1176 }
1177 DEBUGPRINT;
1178 }
1179
RemakeLists_Simple(const set<int> & modified)1180 void NeighborCellLocator::RemakeLists_Simple(const set<int> &modified)
1181 {
1182 DEBUGPRINT;
1183 ASSERT(modified.size() > 0);
1184 // Transform positions to scaled space
1185 vector<Vec> scaledpos(modified.size());
1186 ScaleAndNormalizePositions(modified, scaledpos);
1187
1188 const vector<Vec> &positions = GetWrappedPositions();
1189
1190 // Assign to boxes. Caveat: atom may be outsize the expected
1191 // interval if there are free boundary conditions.
1192 vector<Vec>::const_iterator sp = scaledpos.begin();
1193 for (set<int>::const_iterator a = modified.begin(); a != modified.end();
1194 ++a, ++sp) // Loop over modified and scaledpos
1195 {
1196 // Find the new index of the atom
1197 int index = 0;
1198 for (int i = 0; i < 3; i++)
1199 {
1200 double p = (*sp)[i];
1201 if (p < minimum[i])
1202 p = minimum[i];
1203 if (p > minimum[i] + size[i])
1204 p = minimum[i] + size[i];
1205 int k = int((p - minimum[i]) / size[i] * nCellsTrue[i]);
1206 if (k > nCellsGapStart[i])
1207 {
1208 #if 0
1209 ASSERT (k > nCellsGapStart[i] + nCellsGapSize[i]);
1210 #endif
1211 k -= nCellsGapSize[i];
1212 }
1213 #if 0
1214 ASSERT(k >= 0);
1215 #endif
1216 if (k == nCells[i])
1217 k--;
1218 #if 0
1219 ASSERT(k < nCells[i]);
1220 #endif
1221 index += nTotalCells[i] * k;
1222 }
1223 if (index != cellIndices[*a])
1224 {
1225 // Remove atom from old cell
1226 vector<int>::iterator i = cells[cellIndices[*a]].begin();
1227 vector<int>::iterator terminate = cells[cellIndices[*a]].end();
1228 while ((*i != *a) && (i != terminate))
1229 ++i;
1230 ASSERT(*i == *a);
1231 cells[cellIndices[*a]].erase(i);
1232 // Add to new cell
1233 cells[index].push_back(*a);
1234 cellIndices[*a] = index;
1235 }
1236 // We know that the atom is now in the right box
1237 referencePositions[*a] = positions[*a];
1238 }
1239 DEBUGPRINT;
1240 }
1241
1242 /// Update reference positions of some atoms
UpdateReferencePositions(const set<int> & modified)1243 void NeighborCellLocator::UpdateReferencePositions(const set<int> &modified)
1244 {
1245 const Vec *pos = atoms->GetPositions();
1246 for (set<int>::const_iterator i = modified.begin(); i != modified.end();
1247 ++i)
1248 referencePositions[*i] = pos[*i];
1249 }
1250
print_info(int n)1251 void NeighborCellLocator::print_info(int n)
1252 {
1253 cerr << "NeighborCellLocator info on atom " << n << ":" << endl;
1254 if (referencePositions.size() > n)
1255 cerr << "referencePositions: " << referencePositions[n] << endl;
1256 if (wrappedPositions.size() > n)
1257 cerr << "wrappedPositions: " << wrappedPositions[n] << endl;
1258 if (scaledPositions.size() > n)
1259 cerr << "scaledPositions: " << scaledPositions[n] << endl;
1260 if (offsetPositions.size() > n)
1261 cerr << "offsetPositions: " << offsetPositions[n] << endl;
1262 if (scaledOffsetPositions.size() > n)
1263 cerr << "scaledOffsetPositions: " << scaledOffsetPositions[n] << endl;
1264 cerr << "cellIndex: " << cellIndices[n] << endl;
1265 }
1266
PrintMemory() const1267 long NeighborCellLocator::PrintMemory() const
1268 {
1269 long mem1 = 0;
1270 long mem2 = 0;
1271 long mem3 = 0;
1272 long needed = 0;
1273 mem1 += referencePositions.capacity() * sizeof(Vec);
1274 mem1 += wrappedPositions.capacity() * sizeof(Vec);
1275 mem1 += scaledPositions.capacity() * sizeof(Vec);
1276 mem1 += offsetPositions.capacity() * sizeof(Vec);
1277 mem1 += scaledOffsetPositions.capacity() * sizeof(Vec);
1278 mem2 += cellIndices.capacity() * sizeof(int);
1279 mem2 += cells.capacity() * sizeof(vector<int>);
1280 needed += referencePositions.size() * sizeof(Vec);
1281 needed += wrappedPositions.size() * sizeof(Vec);
1282 needed += scaledPositions.size() * sizeof(Vec);
1283 needed += offsetPositions.size() * sizeof(Vec);
1284 needed += scaledOffsetPositions.size() * sizeof(Vec);
1285 needed += cellIndices.size() * sizeof(int);
1286 needed += cells.size() * sizeof(vector<int>);
1287 int longest = 0;
1288 int empty = 0;
1289 for (vector< vector<int> >::const_iterator i = cells.begin();
1290 i != cells.end(); ++i)
1291 {
1292 mem2 += i->capacity() * sizeof(int);
1293 needed += i->size() * sizeof(int);
1294 if (i->size() > longest)
1295 longest = i->size();
1296 if (i->size() == 0)
1297 empty++;
1298 }
1299 long mem = (mem1 + mem2 + mem3 + 512*1024) / (1024*1024);
1300 mem1 = (mem1 + 512*1024) / (1024*1024);
1301 mem2 = (mem2 + 512*1024) / (1024*1024);
1302 mem3 = (mem3 + 512*1024) / (1024*1024);
1303 long overhead = mem - (needed + 512*1024) / (1024*1024);
1304 char buffer[500];
1305 snprintf(buffer, 500,
1306 "*MEM* NeighborCellLocator %ld MB. [ cells: %ld MB (longest: %d, empty: %d/%d), other: %ld MB, overhead: %ld MB ]",
1307 mem, mem2, longest, empty, (int) cells.size(), mem1, overhead);
1308 cerr << buffer << endl;
1309 return mem;
1310 }
1311
1312