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 <string>
31 : #include <vector>
32 : #include <algorithm>
33 : #include <numeric>
34 : #include <cstdlib>
35 :
36 : namespace PLMD {
37 :
38 310 : NeighborList::NeighborList(const std::vector<AtomNumber>& list0,
39 : const std::vector<AtomNumber>& list1,
40 : const bool& serial,
41 : const bool& do_pair,
42 : const bool& do_pbc,
43 : const Pbc& pbc,
44 : Communicator& cm,
45 : const double& distance,
46 310 : const unsigned& stride)
47 310 : : reduced(false),
48 310 : serial_(serial),
49 310 : do_pair_(do_pair),
50 310 : do_pbc_(do_pbc),
51 310 : twolists_(true),
52 310 : pbc_(&pbc),
53 310 : comm(cm),
54 : //copy-initialize fullatomlist_
55 310 : fullatomlist_(list0),
56 310 : distance_(distance),
57 310 : nlist0_(list0.size()),
58 310 : nlist1_(list1.size()),
59 310 : stride_(stride) {
60 : // store the rest of the atoms into fullatomlist_
61 311 : fullatomlist_.insert(fullatomlist_.end(),list1.begin(),list1.end());
62 310 : if(!do_pair) {
63 271 : nallpairs_=nlist0_*nlist1_;
64 : } else {
65 39 : plumed_assert(nlist0_==nlist1_)
66 : << "when using PAIR option, the two groups should have the same number"
67 : " of elements\n" << "the groups you specified have size "
68 0 : <<nlist0_<<" and "<<nlist1_;
69 39 : nallpairs_=nlist0_;
70 : }
71 310 : initialize();
72 309 : }
73 :
74 44 : NeighborList::NeighborList(const std::vector<AtomNumber>& list0,
75 : const bool& serial, const bool& do_pbc,
76 : const Pbc& pbc,
77 : Communicator& cm,
78 : const double& distance,
79 44 : const unsigned& stride)
80 44 : : reduced(false),
81 44 : serial_(serial),
82 44 : do_pbc_(do_pbc),
83 44 : twolists_(false),
84 44 : pbc_(&pbc),
85 44 : comm(cm),
86 : //copy-initialize fullatomlist_
87 44 : fullatomlist_(list0),
88 44 : distance_(distance),
89 44 : nlist0_(list0.size()),
90 44 : nallpairs_(nlist0_*(nlist0_-1)/2),
91 44 : stride_(stride) {
92 44 : initialize();
93 43 : }
94 :
95 352 : NeighborList::~NeighborList()=default;
96 :
97 354 : void NeighborList::initialize() {
98 : constexpr const char* envKey="PLUMED_IGNORE_NL_MEMORY_ERROR";
99 354 : if(!std::getenv(envKey)) {
100 : //blocking memory allocation on more than 10 GB of memory
101 : //A single list of more than 50000 atoms
102 : //two different lists of more than 35355 atoms each (lista*listb < max, see below)
103 : //that is more than 1250000000 pairs
104 : //each pairIDs occupies 64 bit (where unsigned are 32bit integers)
105 : //4294967296 is max(uint32)+1 and is more than 34 GB (correspond to a system of 65536 atoms)
106 354 : if(nallpairs_ > 1250000000 ) {
107 2 : const unsigned GB = sizeof(decltype(neighbors_)::value_type) * nallpairs_ / 1000000000;
108 2 : plumed_error() << "A NeighborList is trying to allocate "
109 4 : + std::to_string( GB ) +" GB of data for the list of neighbors\n"
110 8 : "You can skip this error by exporting \""+envKey+"\"";
111 : }
112 : }
113 : try {
114 352 : neighbors_.resize(nallpairs_);
115 0 : } catch (...) {
116 0 : plumed_error_nested() << "An error happened while allocating the neighbor "
117 : "list, please decrease the number of atoms used";
118 0 : }
119 : //TODO: test if this is feasible for accelerating the loop
120 : //#pragma omp parallel for default(shared)
121 228150421 : for(unsigned int i=0; i<nallpairs_; ++i) {
122 228150069 : neighbors_[i]=getIndexPair(i);
123 : }
124 352 : }
125 :
126 70018 : std::vector<AtomNumber>& NeighborList::getFullAtomList() {
127 70018 : return fullatomlist_;
128 : }
129 :
130 234387789 : NeighborList::pairIDs NeighborList::getIndexPair(const unsigned ipair) {
131 : pairIDs index;
132 234387789 : if(twolists_ && do_pair_) {
133 12546 : index=pairIDs(ipair,ipair+nlist0_);
134 234375243 : } else if (twolists_ && !do_pair_) {
135 143426597 : index=pairIDs(ipair/nlist1_,ipair%nlist1_+nlist0_);
136 90948646 : } else if (!twolists_) {
137 90948646 : unsigned ii = nallpairs_-1-ipair;
138 90948646 : unsigned K = unsigned(std::floor((std::sqrt(double(8*ii+1))+1)/2));
139 90948646 : unsigned jj = ii-K*(K-1)/2;
140 90948646 : index=pairIDs(nlist0_-1-K,nlist0_-1-jj);
141 : }
142 234387789 : return index;
143 : }
144 :
145 3138 : void NeighborList::update(const std::vector<Vector>& positions) {
146 3138 : neighbors_.clear();
147 3138 : const double d2=distance_*distance_;
148 : // check if positions array has the correct length
149 3138 : plumed_assert(positions.size()==fullatomlist_.size());
150 :
151 3138 : unsigned stride=comm.Get_size();
152 3138 : unsigned rank=comm.Get_rank();
153 3138 : unsigned nt=OpenMP::getNumThreads();
154 3138 : if(serial_) {
155 : stride=1;
156 : rank=0;
157 : nt=1;
158 : }
159 : std::vector<unsigned> local_flat_nl;
160 :
161 3138 : #pragma omp parallel num_threads(nt)
162 : {
163 : std::vector<unsigned> private_flat_nl;
164 : #pragma omp for nowait
165 : for(unsigned int i=rank; i<nallpairs_; i+=stride) {
166 : pairIDs index=getIndexPair(i);
167 : unsigned index0=index.first;
168 : unsigned index1=index.second;
169 : Vector distance;
170 : if(do_pbc_) {
171 : distance=pbc_->distance(positions[index0],positions[index1]);
172 : } else {
173 : distance=delta(positions[index0],positions[index1]);
174 : }
175 : double value=modulo2(distance);
176 : if(value<=d2) {
177 : private_flat_nl.push_back(index0);
178 : private_flat_nl.push_back(index1);
179 : }
180 : }
181 : #pragma omp critical
182 : local_flat_nl.insert(local_flat_nl.end(),
183 : private_flat_nl.begin(),
184 : private_flat_nl.end());
185 : }
186 :
187 : // find total dimension of neighborlist
188 3138 : std::vector <int> local_nl_size(stride, 0);
189 3138 : local_nl_size[rank] = local_flat_nl.size();
190 3138 : if(!serial_) {
191 3126 : comm.Sum(&local_nl_size[0], stride);
192 : }
193 : int tot_size = std::accumulate(local_nl_size.begin(), local_nl_size.end(), 0);
194 3138 : if(tot_size==0) {
195 18 : setRequestList();
196 : return;
197 : }
198 : // merge
199 3120 : std::vector<unsigned> merge_nl(tot_size, 0);
200 : // calculate vector of displacement
201 3120 : std::vector<int> disp(stride);
202 3120 : disp[0] = 0;
203 : int rank_size = 0;
204 9216 : for(unsigned i=0; i<stride-1; ++i) {
205 6096 : rank_size += local_nl_size[i];
206 6096 : disp[i+1] = rank_size;
207 : }
208 : // Allgather neighbor list
209 3120 : if(comm.initialized()&&!serial_) {
210 4181 : comm.Allgatherv((!local_flat_nl.empty()?&local_flat_nl[0]:NULL),
211 : local_nl_size[rank],
212 : &merge_nl[0],
213 : &local_nl_size[0],
214 : &disp[0]);
215 : } else {
216 1016 : merge_nl = local_flat_nl;
217 : }
218 : // resize neighbor stuff
219 3120 : neighbors_.resize(tot_size/2);
220 6114036 : for(unsigned i=0; i<tot_size/2; i++) {
221 6110916 : unsigned j=2*i;
222 6110916 : neighbors_[i] = std::make_pair(merge_nl[j],merge_nl[j+1]);
223 : }
224 :
225 3120 : setRequestList();
226 : }
227 :
228 3138 : void NeighborList::setRequestList() {
229 3138 : requestlist_.clear();
230 6114054 : for(unsigned int i=0; i<size(); ++i) {
231 6110916 : requestlist_.push_back(fullatomlist_[neighbors_[i].first]);
232 6110916 : requestlist_.push_back(fullatomlist_[neighbors_[i].second]);
233 : }
234 3138 : Tools::removeDuplicates(requestlist_);
235 3138 : reduced=false;
236 3138 : }
237 :
238 240 : std::vector<AtomNumber>& NeighborList::getReducedAtomList() {
239 240 : if(!reduced) {
240 168162 : for(unsigned int i=0; i<size(); ++i) {
241 168119 : AtomNumber index0=fullatomlist_[neighbors_[i].first];
242 168119 : AtomNumber index1=fullatomlist_[neighbors_[i].second];
243 : // I exploit the fact that requestlist_ is an ordered vector
244 : // And I assume that index0 and index1 actually exists in the requestlist_ (see setRequestList())
245 : // so I can use lower_bond that uses binary seach instead of find
246 : plumed_dbg_assert(std::is_sorted(requestlist_.begin(),requestlist_.end()));
247 168119 : auto p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index0);
248 : plumed_dbg_assert(p!=requestlist_.end());
249 168119 : unsigned newindex0=p-requestlist_.begin();
250 168119 : p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index1);
251 : plumed_dbg_assert(p!=requestlist_.end());
252 168119 : unsigned newindex1=p-requestlist_.begin();
253 : neighbors_[i]=pairIDs(newindex0,newindex1);
254 : }
255 : }
256 240 : reduced=true;
257 240 : return requestlist_;
258 : }
259 :
260 14950 : unsigned NeighborList::getStride() const {
261 14950 : return stride_;
262 : }
263 :
264 24 : unsigned NeighborList::getLastUpdate() const {
265 24 : return lastupdate_;
266 : }
267 :
268 0 : void NeighborList::setLastUpdate(const unsigned step) {
269 0 : lastupdate_=step;
270 0 : }
271 :
272 23863719 : unsigned NeighborList::size() const {
273 23863719 : return neighbors_.size();
274 : }
275 :
276 262285472 : NeighborList::pairIDs NeighborList::getClosePair(const unsigned i) const {
277 262285472 : return neighbors_[i];
278 : }
279 :
280 : NeighborList::pairAtomNumbers
281 34665936 : NeighborList::getClosePairAtomNumber(const unsigned i) const {
282 : pairAtomNumbers Aneigh=pairAtomNumbers(
283 34665936 : fullatomlist_[neighbors_[i].first],
284 34665936 : fullatomlist_[neighbors_[i].second]);
285 34665936 : return Aneigh;
286 : }
287 :
288 0 : std::vector<unsigned> NeighborList::getNeighbors(const unsigned index) {
289 : std::vector<unsigned> neighbors;
290 0 : for(unsigned int i=0; i<size(); ++i) {
291 0 : if(neighbors_[i].first==index) {
292 0 : neighbors.push_back(neighbors_[i].second);
293 : }
294 0 : if(neighbors_[i].second==index) {
295 0 : neighbors.push_back(neighbors_[i].first);
296 : }
297 : }
298 0 : return neighbors;
299 : }
300 :
301 : } // namespace PLMD
|
