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 "Vector.h"
24 : #include "Pbc.h"
25 : #include "AtomNumber.h"
26 : #include "Communicator.h"
27 : #include "OpenMP.h"
28 : #include "Tools.h"
29 : #include <vector>
30 : #include <algorithm>
31 : #include <numeric>
32 :
33 : namespace PLMD {
34 :
35 255 : NeighborList::NeighborList(const std::vector<AtomNumber>& list0, const std::vector<AtomNumber>& list1,
36 : const bool& serial, const bool& do_pair, const bool& do_pbc, const Pbc& pbc, Communicator& cm,
37 255 : const double& distance, const unsigned& stride): reduced(false),
38 255 : serial_(serial), do_pair_(do_pair), do_pbc_(do_pbc), pbc_(&pbc), comm(cm),
39 255 : distance_(distance), stride_(stride) {
40 : // store full list of atoms needed
41 255 : fullatomlist_=list0;
42 255 : fullatomlist_.insert(fullatomlist_.end(),list1.begin(),list1.end());
43 255 : nlist0_=list0.size();
44 255 : nlist1_=list1.size();
45 255 : twolists_=true;
46 255 : if(!do_pair) {
47 220 : nallpairs_=nlist0_*nlist1_;
48 : } else {
49 35 : plumed_assert(nlist0_==nlist1_) << "when using PAIR option, the two groups should have the same number of elements\n"
50 0 : << "the groups you specified have size "<<nlist0_<<" and "<<nlist1_;
51 35 : nallpairs_=nlist0_;
52 : }
53 255 : initialize();
54 255 : lastupdate_=0;
55 255 : }
56 :
57 39 : NeighborList::NeighborList(const std::vector<AtomNumber>& list0, const bool& serial, const bool& do_pbc,
58 : const Pbc& pbc, Communicator& cm, const double& distance,
59 39 : const unsigned& stride): reduced(false),
60 39 : serial_(serial), do_pbc_(do_pbc), pbc_(&pbc), comm(cm),
61 39 : distance_(distance), stride_(stride) {
62 39 : fullatomlist_=list0;
63 39 : nlist0_=list0.size();
64 39 : twolists_=false;
65 39 : nallpairs_=nlist0_*(nlist0_-1)/2;
66 39 : initialize();
67 39 : lastupdate_=0;
68 39 : }
69 :
70 294 : void NeighborList::initialize() {
71 294 : neighbors_.clear();
72 42954435 : for(unsigned int i=0; i<nallpairs_; ++i) {
73 85908282 : neighbors_.push_back(getIndexPair(i));
74 : }
75 294 : }
76 :
77 69984 : std::vector<AtomNumber>& NeighborList::getFullAtomList() {
78 69984 : return fullatomlist_;
79 : }
80 :
81 49191861 : std::pair<unsigned,unsigned> NeighborList::getIndexPair(unsigned ipair) {
82 : std::pair<unsigned,unsigned> index;
83 49191861 : if(twolists_ && do_pair_) {
84 946 : index=std::pair<unsigned,unsigned>(ipair,ipair+nlist0_);
85 49190915 : } else if (twolists_ && !do_pair_) {
86 8866469 : index=std::pair<unsigned,unsigned>(ipair/nlist1_,ipair%nlist1_+nlist0_);
87 40324446 : } else if (!twolists_) {
88 40324446 : unsigned ii = nallpairs_-1-ipair;
89 40324446 : unsigned K = unsigned(std::floor((std::sqrt(double(8*ii+1))+1)/2));
90 40324446 : unsigned jj = ii-K*(K-1)/2;
91 40324446 : index=std::pair<unsigned,unsigned>(nlist0_-1-K,nlist0_-1-jj);
92 : }
93 49191861 : return index;
94 : }
95 :
96 3138 : void NeighborList::update(const std::vector<Vector>& positions) {
97 3138 : neighbors_.clear();
98 3138 : const double d2=distance_*distance_;
99 : // check if positions array has the correct length
100 3138 : plumed_assert(positions.size()==fullatomlist_.size());
101 :
102 3138 : unsigned stride=comm.Get_size();
103 3138 : unsigned rank=comm.Get_rank();
104 3138 : unsigned nt=OpenMP::getNumThreads();
105 3138 : if(serial_) {
106 : stride=1;
107 : rank=0;
108 : nt=1;
109 : }
110 : std::vector<unsigned> local_flat_nl;
111 :
112 3138 : #pragma omp parallel num_threads(nt)
113 : {
114 : std::vector<unsigned> private_flat_nl;
115 : #pragma omp for nowait
116 : for(unsigned int i=rank; i<nallpairs_; i+=stride) {
117 : std::pair<unsigned,unsigned> index=getIndexPair(i);
118 : unsigned index0=index.first;
119 : unsigned index1=index.second;
120 : Vector distance;
121 : if(do_pbc_) {
122 : distance=pbc_->distance(positions[index0],positions[index1]);
123 : } else {
124 : distance=delta(positions[index0],positions[index1]);
125 : }
126 : double value=modulo2(distance);
127 : if(value<=d2) {
128 : private_flat_nl.push_back(index0);
129 : private_flat_nl.push_back(index1);
130 : }
131 : }
132 : #pragma omp critical
133 : local_flat_nl.insert(local_flat_nl.end(), private_flat_nl.begin(), private_flat_nl.end());
134 : }
135 :
136 : // find total dimension of neighborlist
137 3138 : std::vector <int> local_nl_size(stride, 0);
138 3138 : local_nl_size[rank] = local_flat_nl.size();
139 3138 : if(!serial_) {
140 3126 : comm.Sum(&local_nl_size[0], stride);
141 : }
142 : int tot_size = std::accumulate(local_nl_size.begin(), local_nl_size.end(), 0);
143 3138 : if(tot_size==0) {
144 18 : setRequestList();
145 : return;
146 : }
147 : // merge
148 3120 : std::vector<unsigned> merge_nl(tot_size, 0);
149 : // calculate vector of displacement
150 3120 : std::vector<int> disp(stride);
151 3120 : disp[0] = 0;
152 : int rank_size = 0;
153 9216 : for(unsigned i=0; i<stride-1; ++i) {
154 6096 : rank_size += local_nl_size[i];
155 6096 : disp[i+1] = rank_size;
156 : }
157 : // Allgather neighbor list
158 3120 : if(comm.initialized()&&!serial_) {
159 4181 : comm.Allgatherv((!local_flat_nl.empty()?&local_flat_nl[0]:NULL), local_nl_size[rank], &merge_nl[0], &local_nl_size[0], &disp[0]);
160 : } else {
161 1016 : merge_nl = local_flat_nl;
162 : }
163 : // resize neighbor stuff
164 3120 : neighbors_.resize(tot_size/2);
165 6114036 : for(unsigned i=0; i<tot_size/2; i++) {
166 6110916 : unsigned j=2*i;
167 6110916 : neighbors_[i] = std::make_pair(merge_nl[j],merge_nl[j+1]);
168 : }
169 :
170 3120 : setRequestList();
171 : }
172 :
173 3138 : void NeighborList::setRequestList() {
174 3138 : requestlist_.clear();
175 6114054 : for(unsigned int i=0; i<size(); ++i) {
176 6110916 : requestlist_.push_back(fullatomlist_[neighbors_[i].first]);
177 6110916 : requestlist_.push_back(fullatomlist_[neighbors_[i].second]);
178 : }
179 3138 : Tools::removeDuplicates(requestlist_);
180 3138 : reduced=false;
181 3138 : }
182 :
183 240 : std::vector<AtomNumber>& NeighborList::getReducedAtomList() {
184 240 : if(!reduced)
185 168162 : for(unsigned int i=0; i<size(); ++i) {
186 : unsigned newindex0=0,newindex1=0;
187 168119 : AtomNumber index0=fullatomlist_[neighbors_[i].first];
188 168119 : AtomNumber index1=fullatomlist_[neighbors_[i].second];
189 : // I exploit the fact that requestlist_ is an ordered vector
190 : // And I assume that index0 and index1 actually exists in the requestlist_ (see setRequestList())
191 : // so I can use lower_bond that uses binary seach instead of find
192 : plumed_dbg_assert(std::is_sorted(requestlist_.begin(),requestlist_.end()));
193 168119 : auto p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index0);
194 : plumed_dbg_assert(p!=requestlist_.end());
195 168119 : newindex0=p-requestlist_.begin();
196 168119 : p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index1);
197 : plumed_dbg_assert(p!=requestlist_.end());
198 168119 : newindex1=p-requestlist_.begin();
199 : neighbors_[i]=std::pair<unsigned,unsigned>(newindex0,newindex1);
200 : }
201 240 : reduced=true;
202 240 : return requestlist_;
203 : }
204 :
205 13616 : unsigned NeighborList::getStride() const {
206 13616 : return stride_;
207 : }
208 :
209 0 : unsigned NeighborList::getLastUpdate() const {
210 0 : return lastupdate_;
211 : }
212 :
213 0 : void NeighborList::setLastUpdate(unsigned step) {
214 0 : lastupdate_=step;
215 0 : }
216 :
217 23864320 : unsigned NeighborList::size() const {
218 23864320 : return neighbors_.size();
219 : }
220 :
221 82146772 : std::pair<unsigned,unsigned> NeighborList::getClosePair(unsigned i) const {
222 82146772 : return neighbors_[i];
223 : }
224 :
225 34665936 : std::pair<AtomNumber,AtomNumber> NeighborList::getClosePairAtomNumber(unsigned i) const {
226 : std::pair<AtomNumber,AtomNumber> Aneigh;
227 34665936 : Aneigh=std::pair<AtomNumber,AtomNumber>(fullatomlist_[neighbors_[i].first],fullatomlist_[neighbors_[i].second]);
228 34665936 : return Aneigh;
229 : }
230 :
231 0 : std::vector<unsigned> NeighborList::getNeighbors(unsigned index) {
232 : std::vector<unsigned> neighbors;
233 0 : for(unsigned int i=0; i<size(); ++i) {
234 0 : if(neighbors_[i].first==index) {
235 0 : neighbors.push_back(neighbors_[i].second);
236 : }
237 0 : if(neighbors_[i].second==index) {
238 0 : neighbors.push_back(neighbors_[i].first);
239 : }
240 : }
241 0 : return neighbors;
242 : }
243 :
244 : }
|
