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 214 : NeighborList::NeighborList(const vector<AtomNumber>& list0, const vector<AtomNumber>& list1,
34 : const bool& do_pair, const bool& do_pbc, const Pbc& pbc,
35 : const double& distance, const unsigned& stride): reduced(false),
36 : do_pair_(do_pair), do_pbc_(do_pbc), pbc_(&pbc),
37 428 : distance_(distance), stride_(stride)
38 : {
39 : // store full list of atoms needed
40 214 : fullatomlist_=list0;
41 214 : fullatomlist_.insert(fullatomlist_.end(),list1.begin(),list1.end());
42 214 : nlist0_=list0.size();
43 214 : nlist1_=list1.size();
44 214 : twolists_=true;
45 214 : if(!do_pair) {
46 191 : nallpairs_=nlist0_*nlist1_;
47 : } else {
48 23 : plumed_assert(nlist0_==nlist1_) << "when using PAIR option, the two groups should have the same number of elements\n"
49 0 : << "the groups you specified have size "<<nlist0_<<" and "<<nlist1_;
50 23 : nallpairs_=nlist0_;
51 : }
52 214 : initialize();
53 214 : lastupdate_=0;
54 214 : }
55 :
56 39 : NeighborList::NeighborList(const vector<AtomNumber>& list0, const bool& do_pbc,
57 : const Pbc& pbc, const double& distance,
58 : const unsigned& stride): reduced(false),
59 : do_pbc_(do_pbc), pbc_(&pbc),
60 78 : distance_(distance), stride_(stride) {
61 39 : fullatomlist_=list0;
62 39 : nlist0_=list0.size();
63 39 : twolists_=false;
64 39 : nallpairs_=nlist0_*(nlist0_-1)/2;
65 39 : initialize();
66 39 : lastupdate_=0;
67 39 : }
68 :
69 253 : void NeighborList::initialize() {
70 : neighbors_.clear();
71 42950381 : for(unsigned int i=0; i<nallpairs_; ++i) {
72 85900256 : neighbors_.push_back(getIndexPair(i));
73 : }
74 253 : }
75 :
76 69901 : vector<AtomNumber>& NeighborList::getFullAtomList() {
77 69901 : return fullatomlist_;
78 : }
79 :
80 48976615 : pair<unsigned,unsigned> NeighborList::getIndexPair(unsigned ipair) {
81 : pair<unsigned,unsigned> index;
82 48976615 : if(twolists_ && do_pair_) {
83 634 : index=pair<unsigned,unsigned>(ipair,ipair+nlist0_);
84 48975981 : } else if (twolists_ && !do_pair_) {
85 8651535 : index=pair<unsigned,unsigned>(ipair/nlist1_,ipair%nlist1_+nlist0_);
86 40324446 : } else if (!twolists_) {
87 40324446 : unsigned ii = nallpairs_-1-ipair;
88 40324446 : unsigned K = unsigned(floor((sqrt(double(8*ii+1))+1)/2));
89 40324446 : unsigned jj = ii-K*(K-1)/2;
90 40324446 : index=pair<unsigned,unsigned>(nlist0_-1-K,nlist0_-1-jj);
91 : }
92 48976615 : return index;
93 : }
94 :
95 144 : void NeighborList::update(const vector<Vector>& positions) {
96 : neighbors_.clear();
97 144 : const double d2=distance_*distance_;
98 : // check if positions array has the correct length
99 144 : plumed_assert(positions.size()==fullatomlist_.size());
100 6026487 : for(unsigned int i=0; i<nallpairs_; ++i) {
101 6026487 : pair<unsigned,unsigned> index=getIndexPair(i);
102 6026487 : unsigned index0=index.first;
103 6026487 : unsigned index1=index.second;
104 6026487 : Vector distance;
105 6026487 : if(do_pbc_) {
106 24105948 : distance=pbc_->distance(positions[index0],positions[index1]);
107 : } else {
108 0 : distance=delta(positions[index0],positions[index1]);
109 : }
110 6026487 : double value=modulo2(distance);
111 6026487 : if(value<=d2) {neighbors_.push_back(index);}
112 : }
113 144 : setRequestList();
114 144 : }
115 :
116 144 : void NeighborList::setRequestList() {
117 : requestlist_.clear();
118 5920914 : for(unsigned int i=0; i<size(); ++i) {
119 17762310 : requestlist_.push_back(fullatomlist_[neighbors_[i].first]);
120 11841540 : requestlist_.push_back(fullatomlist_[neighbors_[i].second]);
121 : }
122 144 : Tools::removeDuplicates(requestlist_);
123 144 : reduced=false;
124 144 : }
125 :
126 222 : vector<AtomNumber>& NeighborList::getReducedAtomList() {
127 162491 : if(!reduced)for(unsigned int i=0; i<size(); ++i) {
128 : unsigned newindex0=0,newindex1=0;
129 486807 : AtomNumber index0=fullatomlist_[neighbors_[i].first];
130 324538 : AtomNumber index1=fullatomlist_[neighbors_[i].second];
131 : // I exploit the fact that requestlist_ is an ordered vector
132 324538 : auto p = std::find(requestlist_.begin(), requestlist_.end(), index0); plumed_assert(p!=requestlist_.end()); newindex0=p-requestlist_.begin();
133 324538 : p = std::find(requestlist_.begin(), requestlist_.end(), index1); plumed_assert(p!=requestlist_.end()); newindex1=p-requestlist_.begin();
134 : neighbors_[i]=pair<unsigned,unsigned>(newindex0,newindex1);
135 : }
136 222 : reduced=true;
137 222 : return requestlist_;
138 : }
139 :
140 6776 : unsigned NeighborList::getStride() const {
141 6776 : return stride_;
142 : }
143 :
144 0 : unsigned NeighborList::getLastUpdate() const {
145 0 : return lastupdate_;
146 : }
147 :
148 0 : void NeighborList::setLastUpdate(unsigned step) {
149 0 : lastupdate_=step;
150 0 : }
151 :
152 23510634 : unsigned NeighborList::size() const {
153 23510634 : return neighbors_.size();
154 : }
155 :
156 80546775 : pair<unsigned,unsigned> NeighborList::getClosePair(unsigned i) const {
157 161093550 : return neighbors_[i];
158 : }
159 :
160 34665936 : pair<AtomNumber,AtomNumber> NeighborList::getClosePairAtomNumber(unsigned i) const {
161 : pair<AtomNumber,AtomNumber> Aneigh;
162 103997808 : Aneigh=pair<AtomNumber,AtomNumber>(fullatomlist_[neighbors_[i].first],fullatomlist_[neighbors_[i].second]);
163 34665936 : return Aneigh;
164 : }
165 :
166 0 : vector<unsigned> NeighborList::getNeighbors(unsigned index) {
167 : vector<unsigned> neighbors;
168 0 : for(unsigned int i=0; i<size(); ++i) {
169 0 : if(neighbors_[i].first==index) neighbors.push_back(neighbors_[i].second);
170 0 : if(neighbors_[i].second==index) neighbors.push_back(neighbors_[i].first);
171 : }
172 0 : return neighbors;
173 : }
174 :
175 : }
|
