Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "LinkCells.h"
23 : #include "OpenMP.h"
24 : #include "Exception.h"
25 : #include "Tensor.h"
26 : #include "Tools.h"
27 : #include "View.h"
28 : #include <algorithm>
29 : #include <cmath>
30 : #include <functional>
31 : #include <numeric>
32 : #include <vector>
33 :
34 : namespace PLMD {
35 :
36 751 : LinkCells::LinkCells( Communicator& ) {
37 751 : }
38 :
39 791 : void LinkCells::setCutoff( const double lcut ) {
40 791 : cutoffwasset=true;
41 791 : link_cutoff=lcut;
42 791 : }
43 :
44 504 : double LinkCells::getCutoff() const {
45 504 : plumed_assert( cutoffwasset );
46 504 : return link_cutoff;
47 : }
48 :
49 139 : void LinkCells::setupCells(View<const Vector> pos) {
50 :
51 139 : nopbc=true;
52 139 : Tensor box;
53 : box.zero();
54 : Vector minp, maxp;
55 139 : minp = maxp = pos[0];
56 556 : for(unsigned k=0; k<3; ++k) {
57 26262 : for(unsigned i=1; i<pos.size(); ++i) {
58 25845 : if( pos[i][k]>maxp[k] ) {
59 2031 : maxp[k] = pos[i][k];
60 : }
61 25845 : if( pos[i][k]<minp[k] ) {
62 134 : minp[k] = pos[i][k];
63 : }
64 : }
65 417 : if( link_cutoff<std::sqrt(std::numeric_limits<double>::max()) ) {
66 411 : box[k][k] = link_cutoff*( 1 + std::ceil( (maxp[k] - minp[k])/link_cutoff ) );
67 : } else {
68 6 : box[k][k] = maxp[k] - minp[k] + 1;
69 : }
70 417 : origin[k] = ( minp[k] + maxp[k] ) / 2;
71 : }
72 139 : createCells(box);
73 139 : }
74 :
75 12101 : void LinkCells::setupCells(const Pbc& pbc ) {
76 :
77 12101 : nopbc=false;
78 :
79 12101 : auto box = pbc.getBox();
80 12101 : auto determinant = box.determinant();
81 12101 : plumed_assert(determinant > epsilon) <<"Cell lists cannot be built when passing a box with null volume. Volume is "<<determinant;
82 12101 : createCells(box);
83 12101 : }
84 :
85 12230 : void LinkCells::setupCells(View<const Vector> pos,
86 : const Pbc& pbc ) {
87 12230 : auto box = pbc.getBox();
88 129 : if(box(0,0)==0.0 && box(0,1)==0.0 && box(0,2)==0.0 &&
89 129 : box(1,0)==0.0 && box(1,1)==0.0 && box(1,2)==0.0 &&
90 12359 : box(2,0)==0.0 && box(2,1)==0.0 && box(2,2)==0.0) {
91 : // Create an orthorhombic box around the atomic positions that encompasses every atomic position if there are no pbc
92 129 : setupCells(pos);
93 : } else {
94 : // Setup the pbc object by copying it from action if the Pbc were set
95 12101 : setupCells(pbc);
96 : }
97 12230 : }
98 :
99 12240 : void LinkCells::createCells(const PLMD::Tensor& box) {
100 12240 : mypbc.setBox( box );
101 : // This is the reciprocal lattice
102 : // notice that reciprocal.getRow(0) is a vector that is orthogonal to b and c
103 : // This allows to use linked cells in non orthorhomic boxes
104 12240 : const auto reciprocal = transpose(mypbc.getInvBox());
105 : // the width of the cells along each direction must be at least the cutoff
106 : std::array<double,3> nc;
107 48960 : for(unsigned k=0; k<3; ++k) {
108 36720 : nc[k]=std::floor( 1.0/ reciprocal.getRow(k).modulo() / link_cutoff );
109 : // also catches NaNs
110 36720 : if( !(nc[k]>=1.0) ) {
111 33937 : nc[k]=1.0;
112 : }
113 : }
114 : // Limit the total number of cells: with a small cutoff in a large box the number of
115 : // cells (and the memory and the time spent looping over empty cells) would grow
116 : // with the volume. Larger cells are still correct, they only contain more atoms.
117 12240 : while(nc[0]*nc[1]*nc[2]>maxNumberOfCells) {
118 : const unsigned k=std::distance(nc.begin(),std::max_element(nc.begin(),nc.end()));
119 0 : const double factor=std::cbrt(nc[0]*nc[1]*nc[2]/maxNumberOfCells);
120 0 : nc[k]=std::max(1.0,std::min(nc[k]-1.0,std::floor(nc[k]/factor)));
121 : }
122 48960 : for(unsigned k=0; k<3; ++k) {
123 36720 : ncells[k]=static_cast<unsigned>(nc[k]);
124 : }
125 : // Setup the strides
126 12240 : nstride[0]=1;
127 12240 : nstride[1]=ncells[0];
128 12240 : nstride[2]=ncells[0]*ncells[1];
129 12240 : }
130 :
131 12136 : void LinkCells::buildCellLists( View<const Vector> pos,
132 : View<const unsigned> indices,
133 : const Pbc& pbc ) {
134 12136 : setupCells(pos,pbc);
135 12136 : resetCollection(innerCollection,pos, indices);
136 12136 : }
137 :
138 70 : LinkCells::CellCollection LinkCells::getCollection( View<const Vector> pos,
139 : View<const unsigned> indices) {
140 70 : CellCollection collection;
141 70 : resetCollection(collection,pos, indices);
142 70 : return collection;
143 0 : }
144 :
145 12206 : void LinkCells::resetCollection(LinkCells::CellCollection &collection,
146 : View<const Vector> pos,
147 : View<const unsigned> indices) {
148 12206 : plumed_assert( cutoffwasset && pos.size()==indices.size() );
149 12206 : const auto nat=pos.size();
150 : // Resize and resets the lists
151 12206 : allcells.assign( nat, 0 );
152 12206 : collection.lcell_lists.resize( nat );
153 : // Setup the storage for link cells
154 12206 : const unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
155 12206 : collection.lcell_starts.resize( ncellstot );
156 : // Clear nlcells
157 12206 : collection.lcell_tots.assign( ncellstot, 0 );
158 :
159 : // Find out what cell everyone is in.
160 : // This is done locally on every process: finding the cell of an atom costs a few
161 : // nanoseconds, much less than the two MPI_Allreduce (of the cell of each atom and of
162 : // the number of atoms in each cell) that were needed when the atoms were split among processes.
163 : #ifdef _OPENMP
164 : // without OpenMP nt would be unused (an error with -Werror)
165 12453 : const unsigned nt=(nat>=minAtomsPerThread)?
166 247 : std::min<unsigned>(OpenMP::getNumThreads(),nat/minAtomsPerThread) : 1;
167 : #endif
168 12206 : #pragma omp parallel for num_threads(nt) schedule(static)
169 : for(unsigned i=0; i<nat; ++i) {
170 : // no exceptions inside the parallel region: invalid cells are checked below
171 : allcells[i]=findCellUnchecked( pos[i] );
172 : }
173 1392359 : for(unsigned i=0; i<nat; ++i) {
174 1380153 : if(allcells[i]>=ncellstot) {
175 : // this raises an error with an explanation
176 0 : findMyCell( pos[i] );
177 0 : plumed_error()<<"atom "<<i<<" is not inside a link cell";
178 : }
179 1380153 : collection.lcell_tots[allcells[i]]++;
180 : }
181 : // Now prepare the link cell lists
182 12206 : unsigned tot=0;
183 1639362 : for(unsigned i=0; i<collection.lcell_tots.size(); ++i) {
184 1627156 : collection.lcell_starts[i]=tot;
185 1627156 : tot+=collection.lcell_tots[i];
186 1627156 : collection.lcell_tots[i]=0;
187 : }
188 12206 : plumed_assert( tot==nat ) <<"Total number of atoms found in link cells is "<<tot<<" number of atoms is "<<nat;
189 :
190 : // And setup the link cells properly
191 1392359 : for(unsigned j=0; j<nat; ++j) {
192 1380153 : unsigned myind = collection.lcell_starts[ allcells[j] ] + collection.lcell_tots[ allcells[j] ];
193 1380153 : collection.lcell_lists[ myind ] = indices[j];
194 1380153 : collection.lcell_tots[allcells[j]]++;
195 : }
196 12206 : }
197 :
198 : constexpr inline int min_cell(int celn, const int ncells, const bool usePBC) {
199 4648938 : celn+=(ncells<2)?0:-1;
200 2324983 : return (usePBC) ? celn : std::max(celn,0);
201 : }
202 :
203 : constexpr inline int max_cell(int celn, const int ncells, const bool usePBC) {
204 : //if we are not using pbcs the minimum is 0 and not -1 (if there are only 2 cells)
205 : //so the maximum should be 2, so that cell 0 can check cell 1
206 575061 : celn+=(ncells<3 && usePBC)?1:2;
207 2326593 : return (usePBC) ? celn : std::min(celn,ncells);
208 : }
209 :
210 774823 : void LinkCells::addRequiredCells( const std::array<unsigned,3>& celn,
211 : unsigned& ncells_required,
212 : std::vector<unsigned>& cells_required,
213 : const bool usePbc) const {
214 : #define LINKC_MIN(n) min_cell(celn[n],ncells[n],usePbc)
215 : #define LINKC_MAX(n) max_cell(celn[n],ncells[n],usePbc)
216 : #define LINKC_PBC(n,num) ((n<0)? num-1 : n%num )
217 : unsigned nnew_cells=0;
218 : plumed_dbg_assert((celn[0] < ncells[0]) && (celn[1] < ncells[1]) && (celn[2] < ncells[2]));
219 :
220 1549646 : const int minX=LINKC_MIN(0);
221 774823 : const int maxX=LINKC_MAX(0);
222 1549646 : const int minY=LINKC_MIN(1);
223 774823 : const int maxY=LINKC_MAX(1);
224 1549646 : const int minZ=LINKC_MIN(2);
225 774823 : const int maxZ=LINKC_MAX(2);
226 :
227 2722541 : for(int nx=minX; nx<maxX; ++nx) {
228 1947718 : int xval=LINKC_PBC(nx,ncells[0])*nstride[0];
229 7405855 : for(int ny=minY; ny<maxY; ++ny) {
230 5458137 : int yval=LINKC_PBC(ny,ncells[1])*nstride[1];
231 21404469 : for(int nz=minZ; nz<maxZ; ++nz) {
232 15946332 : int zval=LINKC_PBC(nz,ncells[2])*nstride[2];
233 :
234 15946332 : unsigned mybox=xval+yval+zval;
235 15946332 : plumed_assert(mybox < getNumberOfCells());
236 : bool added=false;
237 15946332 : for(unsigned k=0; k<ncells_required; ++k) {
238 0 : if( mybox==cells_required[k] ) {
239 : added=true;
240 : break;
241 : }
242 : }
243 15946332 : if( !added ) {
244 15946332 : cells_required[ncells_required+nnew_cells]=mybox;
245 15946332 : nnew_cells++;
246 : }
247 : }
248 : }
249 : }
250 774823 : ncells_required += nnew_cells;
251 : #undef LINKC_MIN
252 : #undef LINKC_MAX
253 : #undef LINKC_PBC
254 774823 : }
255 :
256 6750 : void LinkCells::retrieveNeighboringAtoms( const Vector& pos,
257 : std::vector<unsigned>& cell_list,
258 : unsigned& natomsper,
259 : std::vector<unsigned>& atoms ) const {
260 6750 : if( cell_list.size()!=getNumberOfCells() ) {
261 243 : cell_list.resize( getNumberOfCells() );
262 : }
263 6750 : unsigned ncellt=0;
264 6750 : addRequiredCells( findMyCell( pos ), ncellt, cell_list );
265 6750 : retrieveAtomsInCells( ncellt, make_const_view(cell_list), natomsper, atoms );
266 6750 : }
267 :
268 773606 : void LinkCells::retrieveAtomsInCells( const unsigned ncells_required,
269 : View<const unsigned> cells_required,
270 : unsigned& natomsper,
271 : std::vector<unsigned>& atoms,
272 : const unsigned avoidIndex) const {
273 16700757 : for(unsigned i=0; i<ncells_required; ++i) {
274 15927151 : unsigned mybox=cells_required[i];
275 15927151 : auto boxList = innerCollection.getCellIndexes(mybox);
276 15927151 : if (avoidIndex!=std::numeric_limits<unsigned>::max()) {
277 81241016 : for(const unsigned myatom : boxList) {
278 65460176 : if( myatom!=avoidIndex ) { // Ideally would provide an option to not do this
279 64716911 : atoms[natomsper]=myatom;
280 64716911 : natomsper++;
281 : }
282 : }
283 : } else {
284 378097 : for(const unsigned myatom : boxList) {
285 231786 : atoms[natomsper]=myatom;
286 231786 : natomsper++;
287 : }
288 : }
289 : }
290 773606 : }
291 :
292 1380153 : unsigned LinkCells::findCellUnchecked( Vector mypos ) const noexcept {
293 1380153 : if( nopbc ) {
294 8754 : mypos = mypos - origin;
295 : }
296 1380153 : const Vector fpos=mypbc.realToScaled( mypos );
297 1380153 : const unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
298 : unsigned index=0;
299 5520612 : for(unsigned j=0; j<3; ++j) {
300 4140459 : const double c=std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
301 4140459 : if(!(c>=0.0 && c<ncells[j])) {
302 : return ncellstot;
303 : }
304 4140459 : index+=static_cast<unsigned>(c)*nstride[j];
305 : }
306 : return index;
307 : }
308 :
309 774106 : std::array<unsigned,3> LinkCells::findMyCell( Vector mypos ) const {
310 : std::array<unsigned,3> celn;
311 774106 : if( nopbc ) {
312 3381 : mypos = mypos - origin;
313 : }
314 774106 : Vector fpos=mypbc.realToScaled( mypos );
315 3096424 : for(unsigned j=0; j<3; ++j) {
316 2322318 : celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
317 2322318 : plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ) <<"in link cell "<<celn[j]
318 : <<" but should be between 0 and "<<ncells[j]
319 0 : <<" link cell cutoff is "<<link_cutoff
320 : <<" position is "<<fpos[0]<<" "<<fpos[1]<<" "<<fpos[2]
321 0 : <<" box is "<<mypbc.getBox()(0,0)<<" "<<mypbc.getBox()(1,1)<<" "<<mypbc.getBox()(2,2);
322 : }
323 774106 : return celn;
324 : }
325 :
326 2601 : std::array<unsigned,3> LinkCells::findMyCell( unsigned cellIndex ) const {
327 : std::array<unsigned,3> cell;
328 : //std::div returns quotient and remainder of an integer division
329 2601 : auto r=std::div(static_cast<int>(cellIndex),nstride[2]);
330 2601 : cell[2]=r.quot;
331 2601 : r=std::div(r.rem,nstride[1]);
332 2601 : cell[1]=r.quot;
333 : //nstride[0] is hardcoded to be 1;
334 2601 : cell[0]=r.rem;
335 2601 : return cell;
336 : }
337 :
338 4020 : unsigned LinkCells::convertIndicesToIndex( std::array<unsigned,3> cellCoord) const {
339 : //nstride[0] is hardcoded to be 1;
340 : return cellCoord[0]//*nstride[0]
341 4020 : + cellCoord[1]*nstride[1]
342 4020 : + cellCoord[2]*nstride[2];
343 : }
344 :
345 500 : unsigned LinkCells::findCell( const Vector& pos ) const {
346 500 : return convertIndicesToIndex(findMyCell(pos));
347 : }
348 :
349 0 : unsigned LinkCells::getMaxInCell() const {
350 0 : unsigned maxn = innerCollection.lcell_tots[0];
351 0 : for(unsigned i=1; i<innerCollection.lcell_tots.size(); ++i) {
352 0 : if( innerCollection.lcell_tots[i]>maxn ) {
353 : maxn=innerCollection.lcell_tots[i];
354 : }
355 : }
356 0 : return maxn;
357 : }
358 :
359 11886 : void LinkCells::createNeighborList( unsigned nat,
360 : View<const Vector> pos,
361 : View<const unsigned> ind,
362 : View<const unsigned> tind,
363 : View<const Vector> neigh_pos,
364 : View<const unsigned> neigh_ind,
365 : const Pbc& pbc,
366 : unsigned& natoms_per_list,
367 : std::vector<std::size_t>& nlist,
368 : const bool applyCutoff ) {
369 11886 : buildCellLists( neigh_pos, neigh_ind, pbc );
370 : // natoms_per_list = 27*getMaxInCell();
371 : //this should save a little memory
372 11886 : natoms_per_list=innerCollection.getMaximimumCombination(27);
373 11886 : if( natoms_per_list>innerCollection.lcell_lists.size() ) {
374 0 : natoms_per_list = innerCollection.lcell_lists.size();
375 : }
376 11886 : const std::size_t nlist_sz = std::size_t(nat)*( 2 + std::size_t(natoms_per_list) );
377 11886 : nlist.resize( nlist_sz );
378 : // The distance test needs the position of a candidate from its index: the
379 : // cell lists hold the indices of neigh_ind, the positions are in neigh_pos
380 11886 : const bool useCutoff = applyCutoff && link_cutoff<std::sqrt(std::numeric_limits<double>::max());
381 11886 : const double cutoff2 = link_cutoff*link_cutoff;
382 11886 : std::vector<unsigned> posOfIndex;
383 11886 : if( useCutoff ) {
384 497 : unsigned maxind=0;
385 1101961 : for(unsigned j=0; j<neigh_ind.size(); ++j) {
386 1101464 : maxind = std::max( maxind, neigh_ind[j] );
387 : }
388 497 : posOfIndex.resize( maxind+1 );
389 1101961 : for(unsigned j=0; j<neigh_ind.size(); ++j) {
390 1101464 : posOfIndex[neigh_ind[j]] = j;
391 : }
392 : }
393 : // Each central atom writes its own slots of nlist, so the loop is parallel
394 : #ifdef _OPENMP
395 : // without OpenMP nt would be unused (an error with -Werror)
396 12120 : const unsigned nt=(pos.size()>=minAtomsPerThread/16)?
397 234 : std::min<unsigned>(OpenMP::getNumThreads(),pos.size()/(minAtomsPerThread/16)) : 1;
398 : #endif
399 11886 : #pragma omp parallel num_threads(nt)
400 : {
401 : std::vector<unsigned> indices( 1+natoms_per_list );
402 : std::vector<unsigned> cells_required( getNumberOfCells() );
403 : #pragma omp for schedule(static)
404 : for(unsigned i=0; i<pos.size(); ++i) {
405 : unsigned ncells_required=0;
406 : addRequiredCells( findMyCell( pos[i] ), ncells_required, cells_required );
407 : unsigned natoms=1;
408 : indices[0] = ind[i];
409 : retrieveAtomsInCells( ncells_required,
410 : make_const_view(cells_required),
411 : natoms,
412 : indices,
413 : ind[i]);
414 : const std::size_t lstart = nat + std::size_t(tind[i])*(1+std::size_t(natoms_per_list));
415 : if( natoms>neigh_ind.size() ) {
416 : // the cells wrap around the box: the same atom was found more than once,
417 : // use the whole list of candidates instead
418 : natoms = 1 + neigh_ind.size();
419 : indices[0] = ind[i];
420 : std::copy( neigh_ind.begin(), neigh_ind.end(), indices.begin()+1 );
421 : }
422 : if( useCutoff ) {
423 : unsigned nkept=1;
424 : for(unsigned j=1; j<natoms; ++j) {
425 : const Vector d = pbc.distance( pos[i], neigh_pos[posOfIndex[indices[j]]] );
426 : if( d.modulo2()<=cutoff2 ) {
427 : indices[nkept] = indices[j];
428 : nkept++;
429 : }
430 : }
431 : natoms = nkept;
432 : }
433 : nlist[tind[i]] = natoms;
434 : std::copy(indices.begin(),
435 : indices.begin()+natoms,
436 : nlist.begin()+lstart);
437 : }
438 : }
439 11886 : }
440 :
441 11886 : unsigned LinkCells::CellCollection::getMaximimumCombination(const unsigned numCells) const {
442 : //this is not efficient (there is a copy), but in principle it should be called not much times:
443 : //nth_element order by partition until the array is partially sorted
444 : //(it stops when the pivot is at the asked point, so all the elements to the right of the nth element are higher/lower than it)
445 11886 : if (numCells >= lcell_tots.size()) {
446 11496 : return std::accumulate(lcell_tots.begin(),lcell_tots.end(),0);
447 : }
448 390 : auto tmpCopy=lcell_tots;
449 : auto nth_place=tmpCopy.begin() + numCells-1;
450 390 : std::nth_element(tmpCopy.begin(),nth_place,tmpCopy.end(),std::greater<>());
451 :
452 : //the upper limit is not included, hence the +1
453 390 : return std::accumulate(tmpCopy.begin(),nth_place+1,0);
454 390 : }
455 :
456 : } //namespace PLMD
|
