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