Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 "Tools.h"
27 : #include <vector>
28 : #include <algorithm>
29 :
30 : namespace PLMD {
31 : using namespace std;
32 :
33 169 : NeighborList::NeighborList(const vector<AtomNumber>& list0, const vector<AtomNumber>& list1,
34 : const bool& do_pair, const bool& do_pbc, const Pbc& pbc,
35 169 : const double& distance, const unsigned& stride): reduced(false),
36 338 : do_pair_(do_pair), do_pbc_(do_pbc), pbc_(&pbc),
37 676 : distance_(distance), stride_(stride)
38 : {
39 : // store full list of atoms needed
40 169 : fullatomlist_=list0;
41 169 : fullatomlist_.insert(fullatomlist_.end(),list1.begin(),list1.end());
42 169 : nlist0_=list0.size();
43 169 : nlist1_=list1.size();
44 169 : twolists_=true;
45 169 : if(!do_pair) {
46 152 : nallpairs_=nlist0_*nlist1_;
47 : } else {
48 17 : plumed_massert(nlist0_==nlist1_,"when using PAIR option, the two groups should have the same number of elements");
49 17 : nallpairs_=nlist0_;
50 : }
51 169 : initialize();
52 169 : lastupdate_=0;
53 169 : }
54 :
55 9 : NeighborList::NeighborList(const vector<AtomNumber>& list0, const bool& do_pbc,
56 : const Pbc& pbc, const double& distance,
57 9 : const unsigned& stride): reduced(false),
58 9 : do_pbc_(do_pbc), pbc_(&pbc),
59 27 : distance_(distance), stride_(stride) {
60 9 : fullatomlist_=list0;
61 9 : nlist0_=list0.size();
62 9 : twolists_=false;
63 9 : nallpairs_=nlist0_*(nlist0_-1)/2;
64 9 : initialize();
65 9 : lastupdate_=0;
66 9 : }
67 :
68 178 : void NeighborList::initialize() {
69 178 : neighbors_.clear();
70 679798 : for(unsigned int i=0; i<nallpairs_; ++i) {
71 679620 : neighbors_.push_back(getIndexPair(i));
72 : }
73 178 : }
74 :
75 339 : vector<AtomNumber>& NeighborList::getFullAtomList() {
76 339 : return fullatomlist_;
77 : }
78 :
79 386889 : pair<unsigned,unsigned> NeighborList::getIndexPair(unsigned ipair) {
80 : pair<unsigned,unsigned> index;
81 386889 : if(twolists_ && do_pair_) {
82 336 : index=pair<unsigned,unsigned>(ipair,ipair+nlist0_);
83 386553 : } else if (twolists_ && !do_pair_) {
84 380727 : index=pair<unsigned,unsigned>(ipair/nlist1_,ipair%nlist1_+nlist0_);
85 5826 : } else if (!twolists_) {
86 5826 : unsigned ii = nallpairs_-1-ipair;
87 5826 : unsigned K = unsigned(floor((sqrt(double(8*ii+1))+1)/2));
88 5826 : unsigned jj = ii-K*(K-1)/2;
89 5826 : index=pair<unsigned,unsigned>(nlist0_-1-K,nlist0_-1-jj);
90 : }
91 386889 : return index;
92 : }
93 :
94 125 : void NeighborList::update(const vector<Vector>& positions) {
95 125 : neighbors_.clear();
96 125 : const double d2=distance_*distance_;
97 : // check if positions array has the correct length
98 125 : plumed_assert(positions.size()==fullatomlist_.size());
99 94283 : for(unsigned int i=0; i<nallpairs_; ++i) {
100 47079 : pair<unsigned,unsigned> index=getIndexPair(i);
101 47079 : unsigned index0=index.first;
102 47079 : unsigned index1=index.second;
103 47079 : Vector distance;
104 47079 : if(do_pbc_) {
105 141237 : distance=pbc_->distance(positions[index0],positions[index1]);
106 : } else {
107 0 : distance=delta(positions[index0],positions[index1]);
108 : }
109 47079 : double value=modulo2(distance);
110 47079 : if(value<=d2) {neighbors_.push_back(index);}
111 : }
112 125 : setRequestList();
113 125 : }
114 :
115 125 : void NeighborList::setRequestList() {
116 125 : requestlist_.clear();
117 2061 : for(unsigned int i=0; i<size(); ++i) {
118 2904 : requestlist_.push_back(fullatomlist_[neighbors_[i].first]);
119 1936 : requestlist_.push_back(fullatomlist_[neighbors_[i].second]);
120 : }
121 125 : Tools::removeDuplicates(requestlist_);
122 125 : reduced=false;
123 125 : }
124 :
125 19 : vector<AtomNumber>& NeighborList::getReducedAtomList() {
126 638 : if(!reduced)for(unsigned int i=0; i<size(); ++i) {
127 : unsigned newindex0=0,newindex1=0;
128 1857 : AtomNumber index0=fullatomlist_[neighbors_[i].first];
129 1238 : AtomNumber index1=fullatomlist_[neighbors_[i].second];
130 : // I exploit the fact that requestlist_ is an ordered vector
131 1857 : auto p = std::find(requestlist_.begin(), requestlist_.end(), index0); plumed_assert(p!=requestlist_.end()); newindex0=p-requestlist_.begin();
132 1857 : p = std::find(requestlist_.begin(), requestlist_.end(), index1); plumed_assert(p!=requestlist_.end()); newindex1=p-requestlist_.begin();
133 : neighbors_[i]=pair<unsigned,unsigned>(newindex0,newindex1);
134 : }
135 19 : reduced=true;
136 19 : return requestlist_;
137 : }
138 :
139 3636 : unsigned NeighborList::getStride() const {
140 3636 : return stride_;
141 : }
142 :
143 0 : unsigned NeighborList::getLastUpdate() const {
144 0 : return lastupdate_;
145 : }
146 :
147 0 : void NeighborList::setLastUpdate(unsigned step) {
148 0 : lastupdate_=step;
149 0 : }
150 :
151 4931 : unsigned NeighborList::size() const {
152 4931 : return neighbors_.size();
153 : }
154 :
155 12261460 : pair<unsigned,unsigned> NeighborList::getClosePair(unsigned i) const {
156 24522920 : return neighbors_[i];
157 : }
158 :
159 0 : vector<unsigned> NeighborList::getNeighbors(unsigned index) {
160 : vector<unsigned> neighbors;
161 0 : for(unsigned int i=0; i<size(); ++i) {
162 0 : if(neighbors_[i].first==index) neighbors.push_back(neighbors_[i].second);
163 0 : if(neighbors_[i].second==index) neighbors.push_back(neighbors_[i].first);
164 : }
165 0 : return neighbors;
166 : }
167 :
168 : }
|
