| 32 | } |
| 33 | |
| 34 | void CellList::init() { |
| 35 | // Find cell limits |
| 36 | this->xmin = this->xmax = this->positions(0, 0); |
| 37 | this->ymin = this->ymax = this->positions(0, 1); |
| 38 | this->zmin = this->zmax = this->positions(0, 2); |
| 39 | for (int i = 0; i < this->positions.shape(0); i++) { |
| 40 | double x = this->positions(i, 0); |
| 41 | double y = this->positions(i, 1); |
| 42 | double z = this->positions(i, 2); |
| 43 | if (x < this->xmin) { |
| 44 | this->xmin = x; |
| 45 | }; |
| 46 | if (x > this->xmax) { |
| 47 | this->xmax = x; |
| 48 | }; |
| 49 | if (y < this->ymin) { |
| 50 | this->ymin = y; |
| 51 | }; |
| 52 | if (y > this->ymax) { |
| 53 | this->ymax = y; |
| 54 | }; |
| 55 | if (z < this->zmin) { |
| 56 | this->zmin = z; |
| 57 | }; |
| 58 | if (z > this->zmax) { |
| 59 | this->zmax = z; |
| 60 | }; |
| 61 | }; |
| 62 | |
| 63 | // Add small padding to avoid floating point precision problems at the |
| 64 | // boundary |
| 65 | double padding = 0.0001; |
| 66 | this->xmin -= padding; |
| 67 | this->xmax += padding; |
| 68 | this->ymin -= padding; |
| 69 | this->ymax += padding; |
| 70 | this->zmin -= padding; |
| 71 | this->zmax += padding; |
| 72 | |
| 73 | // Determine amount and size of bins. The bins are made to be always of equal size. |
| 74 | this->nx = max(1, int((this->xmax - this->xmin)/this->cutoff)); |
| 75 | this->ny = max(1, int((this->ymax - this->ymin)/this->cutoff)); |
| 76 | this->nz = max(1, int((this->zmax - this->zmin)/this->cutoff)); |
| 77 | this->dx = max(this->cutoff, (this->xmax - this->xmin)/this->nx); |
| 78 | this->dy = max(this->cutoff, (this->ymax - this->ymin)/this->ny); |
| 79 | this->dz = max(this->cutoff, (this->zmax - this->zmin)/this->nz); |
| 80 | |
| 81 | // Initialize the bin data structure. It is a 4D vector. |
| 82 | this->bins = vector<vector<vector<vector<int>>>>(this->nx, vector<vector<vector<int>>>(this->ny, vector<vector<int>>(this->nz, vector<int>()))); |
| 83 | |
| 84 | // Fill the bins with atom indices |
| 85 | for (int idx = 0; idx < this->positions.shape(0); idx++) { |
| 86 | double x = this->positions(idx, 0); |
| 87 | double y = this->positions(idx, 1); |
| 88 | double z = this->positions(idx, 2); |
| 89 | |
| 90 | // Get bin index |
| 91 | int i = (x - this->xmin)/this->dx; |
no test coverage detected