Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "NeighborList.h"
23 : #include "Exception.h"
24 : #include "Vector.h"
25 : #include "Pbc.h"
26 : #include "AtomNumber.h"
27 : #include "Communicator.h"
28 : #include "OpenMP.h"
29 : #include "Tools.h"
30 : #include "core/Colvar.h"
31 : #include <string>
32 : #include "LinkCells.h"
33 : #include "View.h"
34 : #include <vector>
35 : #include <array>
36 : #include <algorithm>
37 : #include <numeric>
38 : #include <cstdlib>
39 : #include <limits>
40 :
41 : #pragma GCC diagnostic error "-Wswitch"
42 :
43 : namespace PLMD {
44 :
45 2473 : NeighborList::NeighborList(const std::vector<AtomNumber>& list0,
46 : const std::vector<AtomNumber>& list1,
47 : const bool serial,
48 : const bool do_pair,
49 : const bool do_pbc,
50 : const Pbc& pbc,
51 : Communicator& cm,
52 : const double distance,
53 : const unsigned stride,
54 2473 : const bool doCells)
55 2473 : : serial_(serial),
56 2473 : do_pbc_(do_pbc),
57 2473 : useCellList_(doCells),
58 4946 : style_(do_pair ? NNStyle::Pair : NNStyle::TwoList),
59 2473 : pbc_(&pbc),
60 2473 : comm(cm),
61 : //copy-initialize fullatomlist_
62 2473 : fullatomlist_(list0),
63 2473 : distance_(distance),
64 2473 : nlist0_(list0.size()),
65 2473 : nlist1_(list1.size()),
66 4946 : stride_(stride) {
67 : // store the rest of the atoms into fullatomlist_
68 2473 : fullatomlist_.insert(fullatomlist_.end(),list1.begin(),list1.end());
69 2473 : if(style_ != NNStyle::Pair) {
70 1366 : nallpairs_=nlist0_*nlist1_;
71 : } else {
72 1107 : plumed_assert(!useCellList_ ) << "You can't use the cell implementation in Pair mode";
73 1107 : plumed_assert(nlist0_==nlist1_)
74 : << "when using PAIR option, the two groups should have the same number"
75 : " of elements\n" << "the groups you specified have size "
76 0 : <<nlist0_<<" and "<<nlist1_;
77 1107 : nallpairs_=nlist0_;
78 : }
79 2473 : initialize();
80 2475 : }
81 :
82 1112 : NeighborList::NeighborList(const std::vector<AtomNumber>& list0,
83 : const bool serial,
84 : const bool do_pbc,
85 : const Pbc& pbc,
86 : Communicator& cm,
87 : const double distance,
88 : const unsigned stride,
89 1112 : const bool doCells)
90 1112 : : serial_(serial),
91 1112 : do_pbc_(do_pbc),
92 1112 : useCellList_(doCells),
93 1112 : style_(NNStyle::SingleList),
94 1112 : pbc_(&pbc),
95 1112 : comm(cm),
96 : //copy-initialize fullatomlist_
97 1112 : fullatomlist_(list0),
98 1112 : distance_(distance),
99 1112 : nlist0_(list0.size()),
100 1112 : nallpairs_(nlist0_*(nlist0_-1)/2),
101 1112 : stride_(stride) {
102 1112 : initialize();
103 1114 : }
104 :
105 3583 : NeighborList::~NeighborList()=default;
106 :
107 3585 : void NeighborList::initialize() {
108 : constexpr const char* envKey="PLUMED_IGNORE_NL_MEMORY_ERROR";
109 : //this checks the upper limit of the memory
110 3585 : if(!std::getenv(envKey)) {
111 : //blocking memory allocation on more than 10 GB of memory
112 : //A single list of more than 50000 atoms
113 : //two different lists of more than 35355 atoms each (lista*listb < max, see below)
114 : //that is more than 1250000000 pairs
115 : //each pairIDs occupies 64 bit (where unsigned are 32bit integers)
116 : //4294967296 is max(uint32)+1 and is more than 34 GB (correspond to a system of 65536 atoms)
117 3585 : if(nallpairs_ > 1250000000 ) {
118 2 : const unsigned GB = sizeof(decltype(neighbors_)::value_type) * nallpairs_ / 1000000000;
119 2 : plumed_error() << "A NeighborList is trying to allocate "
120 2 : + std::to_string( GB ) +" GB of data for the list of neighbors\n"
121 4 : "You can skip this error by exporting \""+envKey+"\"";
122 : }
123 : }
124 : // pairs are indexed with unsigned integers (getClosePair, size): this limit cannot be skipped
125 3583 : plumed_assert(nallpairs_ <= std::numeric_limits<unsigned>::max())
126 0 : << "A NeighborList cannot handle " << nallpairs_ << " pairs (the maximum is "
127 0 : << std::numeric_limits<unsigned>::max() << ")";
128 : try {
129 : //preallocating some memory now to not spend too much time in the first update
130 : // If the memory is underestimated it should be less than doubled if needed, the factor should depend on the implementation
131 : //Asserting 20 maximun neighbors to occupy less memory
132 5773 : neighbors_.resize((stride_==0 || style_ == NNStyle::Pair) ?
133 : nallpairs_
134 4374 : : (std::min((nlist0_*((useCellList_)?100:30)),nallpairs_)));
135 0 : } catch (...) {
136 0 : plumed_error_nested() << "An error happened while allocating the neighbor "
137 : "list, please decrease the number of atoms used";
138 0 : }
139 3583 : if (nallpairs_ == neighbors_.size()) {
140 : //TODO: test if this is feasible for accelerating the loop
141 : //#pragma omp parallel for default(shared)
142 188180038 : for(unsigned int i=0; i<nallpairs_; ++i) {
143 188177455 : neighbors_[i]=getIndexPair(i);
144 : }
145 2583 : listBuilded=true;
146 : }
147 3583 : }
148 :
149 70155 : std::vector<AtomNumber>& NeighborList::getFullAtomList() {
150 70155 : return fullatomlist_;
151 : }
152 :
153 231031963 : NeighborList::pairIDs NeighborList::getIndexPair(const unsigned ipair) const {
154 : pairIDs index;
155 231031963 : switch (style_) {
156 5591521 : case NNStyle::Pair : {
157 5591521 : index=pairIDs(ipair,ipair+nlist0_);
158 : }
159 5591521 : break;
160 155119056 : case NNStyle::TwoList : {
161 155119056 : index=pairIDs(ipair/nlist1_,ipair%nlist1_+nlist0_);
162 : }
163 155119056 : break;
164 70321386 : case NNStyle::SingleList: {
165 : // 64-bit arithmetic: with more than 2^15 atoms 8*ii does not fit in 32 bits
166 70321386 : const std::size_t ii = nallpairs_-1-std::size_t(ipair);
167 : // K is the largest integer with K*(K-1)/2 <= ii
168 70321386 : std::size_t K = std::size_t(std::floor((std::sqrt(8.0*double(ii)+1.0)+1.0)/2.0));
169 : // guard against rounding in the square root
170 70321386 : while(K*(K-1)/2>ii) {
171 : --K;
172 : }
173 70321386 : while((K+1)*K/2<=ii) {
174 : ++K;
175 : }
176 70321386 : const std::size_t jj = ii-K*(K-1)/2;
177 70321386 : index=pairIDs(nlist0_-1-K,nlist0_-1-jj);
178 : }
179 : }
180 231031963 : return index;
181 : }
182 :
183 9454 : void NeighborList::update(const std::vector<Vector>& positions) {
184 9454 : neighbors_.clear();
185 : // check if positions array has the correct length
186 9454 : plumed_assert(positions.size()==fullatomlist_.size());
187 :
188 9454 : const unsigned nt=(serial_)? 1 : OpenMP::getNumThreads();
189 9454 : const unsigned stride=(serial_)? 1 : comm.Get_size();
190 9454 : const unsigned rank =(serial_)? 0 : comm.Get_rank();
191 9454 : const double d2=distance_*distance_;
192 : // Each thread stores its pairs (as a flat list of indices) in its own buffer.
193 : // Threads work on contiguous blocks (static schedule) and the buffers are
194 : // concatenated in thread order, and the lists of the processes in rank order:
195 : // the order of the pairs does not depend on the number of threads and processes,
196 : // and the results are reproducible.
197 9454 : std::vector<std::vector<unsigned>> thread_flat_nl(nt);
198 :
199 9454 : if(useCellList_) {
200 36 : std::vector<unsigned> indexesForCells(fullatomlist_.size());
201 : std::iota(indexesForCells.begin(),indexesForCells.end(),0);
202 36 : LinkCells cells(comm);
203 36 : cells.setCutoff(distance_);
204 36 : if(do_pbc_) {
205 26 : cells.setupCells(make_const_view(positions),*pbc_);
206 : } else {
207 : // without pbc the cells must enclose the atoms: cells built on the periodic
208 : // box would split close pairs across the planes where scaled coordinates wrap
209 10 : cells.setupCells(make_const_view(positions));
210 : }
211 36 : plumed_assert(style_!=NNStyle::Pair) << "Cell list should not be active with a Pair NL";
212 36 : const bool twoLists = style_==NNStyle::TwoList;
213 36 : const std::size_t nA=nlist0_;
214 36 : LinkCells::CellCollection listA;
215 36 : LinkCells::CellCollection listB;
216 36 : if(twoLists) {
217 52 : listA = cells.getCollection(View{positions.data(),nlist0_},
218 26 : View<const unsigned> {indexesForCells.data(),nlist0_});
219 26 : listB = cells.getCollection(View{positions.data()+nlist0_,nlist1_},
220 26 : View<const unsigned> {indexesForCells.data()+nlist0_,nlist1_});
221 26 : plumed_assert((listA.lcell_lists.size()+listB.lcell_lists.size()) == positions.size())
222 0 : << listA.lcell_lists.size()<<"+"<<listB.lcell_lists.size() <<"==" <<positions.size();
223 : } else {
224 20 : listA = cells.getCollection(View{positions.data(),nA},
225 10 : View<const unsigned> {indexesForCells.data(),nA});
226 : }
227 36 : const LinkCells::CellCollection& candidates = twoLists ? listB : listA;
228 : // cells are split among processes
229 36 : const unsigned ncells=cells.getNumberOfCells();
230 36 : const unsigned cellsPerRank=(ncells+stride-1)/stride;
231 36 : const unsigned cstart=std::min(rank*cellsPerRank,ncells);
232 36 : const unsigned cend=std::min(cstart+cellsPerRank,ncells);
233 36 : #pragma omp parallel num_threads(nt)
234 : {
235 : auto & private_flat_nl=thread_flat_nl[OpenMP::getThreadNum()];
236 : std::vector<unsigned> cells_required(27);
237 : #pragma omp for schedule(static) nowait
238 : for(unsigned c=cstart; c<cend; ++c) {
239 : auto atomsInC= listA.getCellIndexes(c);
240 : if(atomsInC.size()==0) {
241 : continue;
242 : }
243 : auto cell = cells.findMyCell(c);
244 : unsigned ncells_required=0;
245 : cells.addRequiredCells(cell,ncells_required, cells_required,do_pbc_);
246 : for (auto A : atomsInC) {
247 : for (unsigned cb=0; cb <ncells_required ; ++cb) {
248 : for (auto B : candidates.getCellIndexes(cells_required[cb])) {
249 : if(!twoLists && B<=A) {
250 : continue;
251 : }
252 : // the cells contain all the pairs within the cutoff, and many more:
253 : // only the pairs within the cutoff are kept
254 : const Vector distance=do_pbc_ ? pbc_->distance(positions[A],positions[B]) :
255 : delta(positions[A],positions[B]);
256 : if(modulo2(distance)<=d2) {
257 : private_flat_nl.push_back(A);
258 : private_flat_nl.push_back(B);
259 : }
260 : }
261 : }
262 : }
263 : }
264 : }
265 36 : } else {
266 9418 : const std::size_t elementsPerRank = (nallpairs_+stride-1)/stride;
267 9418 : const unsigned start= std::min<std::size_t>(rank*elementsPerRank,nallpairs_);
268 9418 : const unsigned end = std::min<std::size_t>(start + elementsPerRank, nallpairs_);
269 9418 : constexpr unsigned batchSize=1024;
270 9418 : const unsigned work=end>start ? end-start : 0;
271 : // Small lists do not amortize initializing the per-thread batch buffers.
272 : // PAIR lists are linear searches and also perform better without batching.
273 9418 : const bool batchPbc=do_pbc_ && style_!=NNStyle::Pair && work/nt>=batchSize;
274 :
275 9418 : #pragma omp parallel num_threads(nt)
276 : {
277 : auto & private_flat_nl=thread_flat_nl[OpenMP::getThreadNum()];
278 : if(batchPbc) {
279 : // Bound the scratch space and keep enough batches for all threads.
280 : const unsigned batchStep=nt>1 ? std::min(batchSize,std::max(1u,work/(4*nt))) : batchSize;
281 : std::array<Vector,batchSize> distances;
282 : std::array<pairIDs,batchSize> pairs;
283 : #pragma omp for schedule(static) nowait
284 : for(unsigned begin=start; begin<end; begin+=batchStep) {
285 : const unsigned count=std::min(batchStep,end-begin);
286 : for(unsigned j=0; j<count; ++j) {
287 : const auto [index0,index1]=getIndexPair(begin+j);
288 : pairs[j]= {index0,index1};
289 : distances[j]=delta(positions[index0],positions[index1]);
290 : }
291 : pbc_->apply(VectorView(distances[0].data(),count));
292 : for(unsigned j=0; j<count; ++j) {
293 : if(modulo2(distances[j])<=d2) {
294 : private_flat_nl.push_back(pairs[j].first);
295 : private_flat_nl.push_back(pairs[j].second);
296 : }
297 : }
298 : }
299 : } else {
300 : #pragma omp for schedule(static) nowait
301 : for(unsigned int i=start; i<end; ++i) {
302 : auto [index0, index1 ] = getIndexPair(i);
303 : const Vector distance=do_pbc_ ? pbc_->distance(positions[index0],positions[index1]) :
304 : delta(positions[index0],positions[index1]);
305 : double value=modulo2(distance);
306 : if(value<=d2) {
307 : private_flat_nl.push_back(index0);
308 : private_flat_nl.push_back(index1);
309 : }
310 : }
311 : }
312 : }
313 : }
314 :
315 : // concatenate the lists of the threads, in thread order
316 9454 : std::vector<unsigned> local_flat_nl;
317 9454 : if(nt==1) {
318 : local_flat_nl=std::move(thread_flat_nl[0]);
319 : } else {
320 : std::size_t local_size=0;
321 25290 : for(const auto & t : thread_flat_nl) {
322 18920 : local_size+=t.size();
323 : }
324 6370 : local_flat_nl.reserve(local_size);
325 25290 : for(const auto & t : thread_flat_nl) {
326 18920 : local_flat_nl.insert(local_flat_nl.end(),t.begin(),t.end());
327 : }
328 : }
329 :
330 : // find total dimension of neighborlist
331 9454 : std::vector <int> local_nl_size(stride, 0);
332 9454 : local_nl_size[rank] = local_flat_nl.size();
333 9454 : if(!serial_) {
334 6370 : comm.Sum(&local_nl_size[0], stride);
335 : }
336 : int tot_size = std::accumulate(local_nl_size.begin(), local_nl_size.end(), 0);
337 9454 : if(tot_size!=0) {
338 : // merge
339 6692 : std::vector<unsigned> merge_nl;
340 6692 : if(comm.initialized()&&!serial_ && stride>1) {
341 3415 : merge_nl.assign(tot_size, 0);
342 : // calculate vector of displacement
343 3415 : std::vector<int> disp(stride);
344 3415 : disp[0] = 0;
345 : int rank_size = 0;
346 12133 : for(unsigned i=0; i<stride-1; ++i) {
347 8718 : rank_size += local_nl_size[i];
348 8718 : disp[i+1] = rank_size;
349 : }
350 : // Allgather neighbor list
351 6466 : comm.Allgatherv((!local_flat_nl.empty()?&local_flat_nl[0]:NULL),
352 : local_nl_size[rank],
353 : &merge_nl[0],
354 : &local_nl_size[0],
355 : &disp[0]);
356 : } else {
357 : merge_nl = std::move(local_flat_nl);
358 : }
359 : // resize neighbor stuff
360 6692 : neighbors_.resize(tot_size/2);
361 12581732 : for(int i=0; i<tot_size/2; i++) {
362 12575040 : unsigned j=2*i;
363 12575040 : neighbors_[i] = std::make_pair(merge_nl[j],merge_nl[j+1]);
364 : }
365 : }
366 9454 : listBuilded=true;
367 : // the reduced list changed, so it must be requested again
368 9454 : reducedRequested=false;
369 9454 : if (stride_ >1) {
370 230 : setRequestList();
371 : } else {
372 9224 : reduced=true;
373 : }
374 9454 : }
375 :
376 230 : void NeighborList::setRequestList() {
377 : // at time of adding the `if (stride_>1)` in `update()`
378 : // this function is called only from `update()` and it is private
379 : // so as now it is not necessary to add extra logic in this function
380 230 : requestlist_.clear();
381 267940 : for(unsigned int i=0; i<size(); ++i) {
382 267710 : requestlist_.push_back(fullatomlist_[neighbors_[i].first]);
383 267710 : requestlist_.push_back(fullatomlist_[neighbors_[i].second]);
384 : }
385 230 : Tools::removeDuplicates(requestlist_);
386 230 : reduced=false;
387 230 : }
388 :
389 121 : std::vector<AtomNumber>& NeighborList::getReducedAtomList() {
390 121 : if(stride_>1) {
391 121 : if(!reduced) {
392 255792 : for(unsigned int i=0; i<size(); ++i) {
393 255671 : AtomNumber index0=fullatomlist_[neighbors_[i].first];
394 255671 : AtomNumber index1=fullatomlist_[neighbors_[i].second];
395 : // I exploit the fact that requestlist_ is an ordered vector
396 : // And I assume that index0 and index1 actually exists in the requestlist_ (see setRequestList())
397 : // so I can use lower_bond that uses binary seach instead of find
398 : plumed_dbg_assert(std::is_sorted(requestlist_.begin(),requestlist_.end()));
399 255671 : auto p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index0);
400 : plumed_dbg_assert(p!=requestlist_.end());
401 255671 : unsigned newindex0=p-requestlist_.begin();
402 255671 : p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index1);
403 : plumed_dbg_assert(p!=requestlist_.end());
404 255671 : unsigned newindex1=p-requestlist_.begin();
405 : neighbors_[i]=pairIDs(newindex0,newindex1);
406 : }
407 : }
408 121 : reduced=true;
409 121 : return requestlist_;
410 : } else {
411 : //safeguard in case the stride is set to 0 or 1
412 0 : return getFullAtomList();
413 : }
414 : }
415 :
416 11056 : unsigned NeighborList::getStride() const {
417 11056 : return stride_;
418 : }
419 :
420 24 : unsigned NeighborList::getLastUpdate() const {
421 24 : return lastupdate_;
422 : }
423 :
424 0 : void NeighborList::setLastUpdate(const unsigned step) {
425 0 : lastupdate_=step;
426 0 : }
427 :
428 29911459 : unsigned NeighborList::size() const {
429 29911459 : if(!listBuilded) {
430 11515245 : return nallpairs_;
431 : }
432 18396214 : return neighbors_.size();
433 :
434 : }
435 :
436 239402863 : NeighborList::pairIDs NeighborList::getClosePair(const unsigned i) const {
437 239402863 : if(listBuilded) {
438 239339647 : return neighbors_[i];
439 : } else {
440 63216 : return getIndexPair(i);
441 : }
442 : }
443 :
444 15827028 : NeighborList::pairIDs NeighborList::getUpdatedPair(const unsigned i) const {
445 : plumed_dbg_assert(listBuilded) << "NeighborList::getUpdatedPair should be called after update()";
446 15827028 : return neighbors_[i];
447 : }
448 :
449 140 : bool NeighborList::ready() const {
450 140 : return listBuilded || neighbors_.size()==nallpairs_;
451 : }
452 :
453 : NeighborList::pairAtomNumbers
454 34665936 : NeighborList::getClosePairAtomNumber(const unsigned i) const {
455 34665936 : if(listBuilded)
456 : return pairAtomNumbers{
457 11635536 : fullatomlist_[neighbors_[i].first],
458 11635536 : fullatomlist_[neighbors_[i].second]};
459 : else {
460 23030400 : auto p = getIndexPair(i);
461 : return pairAtomNumbers{
462 23030400 : fullatomlist_[p.first],
463 23030400 : fullatomlist_[p.second]};
464 : }
465 : }
466 :
467 8000 : std::vector<unsigned> NeighborList::getNeighbors(const unsigned index) const {
468 8000 : std::vector<unsigned> neighbors;
469 8000 : if (listBuilded) {
470 5377590 : for(unsigned int i=0; i<size(); ++i) {
471 5369590 : if(neighbors_[i].first==index) {
472 48580 : neighbors.push_back(neighbors_[i].second);
473 : }
474 5369590 : if(neighbors_[i].second==index) {
475 48580 : neighbors.push_back(neighbors_[i].first);
476 : }
477 : }
478 : } else {
479 0 : for(unsigned int i=0; i<size(); ++i) {
480 0 : auto neigh = getClosePair(i);
481 0 : if(neigh.first==index) {
482 0 : neighbors.push_back(neigh.second);
483 : }
484 0 : if(neigh.second==index) {
485 0 : neighbors.push_back(neigh.first);
486 : }
487 : }
488 : }
489 8000 : return neighbors;
490 : }
491 :
492 4874 : NeighborList::preparestatus NeighborList::prepare(Colvar* const aa,
493 : bool firsttime,
494 : bool invalidateList) {
495 4874 : if(stride_>0) {
496 714 : if(stride_==1) {
497 : invalidateList=true;
498 : firsttime=false;
499 660 : } else if(firsttime || (aa->getStep()%stride_==0 )) {
500 264 : aa->requestAtoms(getFullAtomList());
501 264 : reducedRequested=false;
502 : invalidateList=true;
503 : firsttime=false;
504 : } else {
505 : // The reduced list only changes when the list is updated, so it is requested once after each update.
506 : // Requesting the same atoms at every step would make the action rebuild its lists of atoms,
507 : // and PLUMED rebuild the lists of the atoms that are shared with the MD code, at every step.
508 396 : if( !reducedRequested ) {
509 121 : aa->requestAtoms(getReducedAtomList());
510 121 : reducedRequested=true;
511 : }
512 : invalidateList=false;
513 396 : if(aa->getExchangeStep()) {
514 0 : aa->error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
515 : }
516 : }
517 714 : if(aa->getExchangeStep()) {
518 : firsttime=true;
519 : }
520 : }
521 4874 : return {firsttime,invalidateList};
522 : }
523 : } // namespace PLMD
|
