Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 127 : 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 127 : distance_(distance), stride_(stride)
38 : {
39 : // store full list of atoms needed
40 127 : fullatomlist_=list0;
41 127 : fullatomlist_.insert(fullatomlist_.end(),list1.begin(),list1.end());
42 127 : nlist0_=list0.size();
43 127 : nlist1_=list1.size();
44 127 : twolists_=true;
45 127 : if(!do_pair) {
46 115 : nallpairs_=nlist0_*nlist1_;
47 : } else {
48 12 : plumed_massert(nlist0_==nlist1_,"when using PAIR option, the two groups should have the same number of elements");
49 12 : nallpairs_=nlist0_;
50 : }
51 127 : initialize();
52 127 : lastupdate_=0;
53 127 : }
54 :
55 9 : NeighborList::NeighborList(const vector<AtomNumber>& list0, const bool& do_pbc,
56 : const Pbc& pbc, const double& distance,
57 : const unsigned& stride): reduced(false),
58 : do_pbc_(do_pbc), pbc_(&pbc),
59 9 : 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 136 : void NeighborList::initialize() {
69 136 : neighbors_.clear();
70 258592 : for(unsigned int i=0; i<nallpairs_; ++i) {
71 258456 : neighbors_.push_back(getIndexPair(i));
72 : }
73 136 : }
74 :
75 201 : vector<AtomNumber>& NeighborList::getFullAtomList() {
76 201 : return fullatomlist_;
77 : }
78 :
79 260175 : pair<unsigned,unsigned> NeighborList::getIndexPair(unsigned ipair) {
80 260175 : pair<unsigned,unsigned> index;
81 260175 : if(twolists_ && do_pair_) {
82 319 : index=pair<unsigned,unsigned>(ipair,ipair+nlist0_);
83 259856 : } else if (twolists_ && !do_pair_) {
84 254030 : 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 260175 : return index;
92 : }
93 :
94 53 : void NeighborList::update(const vector<Vector>& positions) {
95 53 : neighbors_.clear();
96 53 : const double d2=distance_*distance_;
97 : // check if positions array has the correct length
98 53 : plumed_assert(positions.size()==fullatomlist_.size());
99 1772 : for(unsigned int i=0; i<nallpairs_; ++i) {
100 1719 : pair<unsigned,unsigned> index=getIndexPair(i);
101 1719 : unsigned index0=index.first;
102 1719 : unsigned index1=index.second;
103 1719 : Vector distance;
104 1719 : if(do_pbc_) {
105 1719 : distance=pbc_->distance(positions[index0],positions[index1]);
106 : } else {
107 0 : distance=delta(positions[index0],positions[index1]);
108 : }
109 1719 : double value=modulo2(distance);
110 1719 : if(value<=d2) {neighbors_.push_back(index);}
111 : }
112 53 : setRequestList();
113 53 : }
114 :
115 53 : void NeighborList::setRequestList() {
116 53 : requestlist_.clear();
117 949 : for(unsigned int i=0; i<size(); ++i) {
118 896 : requestlist_.push_back(fullatomlist_[neighbors_[i].first]);
119 896 : requestlist_.push_back(fullatomlist_[neighbors_[i].second]);
120 : }
121 53 : Tools::removeDuplicates(requestlist_);
122 53 : reduced=false;
123 53 : }
124 :
125 13 : vector<AtomNumber>& NeighborList::getReducedAtomList() {
126 626 : if(!reduced)for(unsigned int i=0; i<size(); ++i) {
127 613 : unsigned newindex0=0,newindex1=0;
128 613 : AtomNumber index0=fullatomlist_[neighbors_[i].first];
129 613 : AtomNumber index1=fullatomlist_[neighbors_[i].second];
130 : // I exploit the fact that requestlist_ is an ordered vector
131 613 : vector<AtomNumber>::iterator p;
132 613 : p = std::find(requestlist_.begin(), requestlist_.end(), index0); plumed_assert(p!=requestlist_.end()); newindex0=p-requestlist_.begin();
133 613 : p = std::find(requestlist_.begin(), requestlist_.end(), index1); plumed_assert(p!=requestlist_.end()); newindex1=p-requestlist_.begin();
134 613 : neighbors_[i]=pair<unsigned,unsigned>(newindex0,newindex1);
135 : }
136 13 : reduced=true;
137 13 : return requestlist_;
138 : }
139 :
140 2982 : unsigned NeighborList::getStride() const {
141 2982 : 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 4505 : unsigned NeighborList::size() const {
153 4505 : return neighbors_.size();
154 : }
155 :
156 11509693 : pair<unsigned,unsigned> NeighborList::getClosePair(unsigned i) const {
157 11509693 : return neighbors_[i];
158 : }
159 :
160 0 : vector<unsigned> NeighborList::getNeighbors(unsigned index) {
161 0 : vector<unsigned> neighbors;
162 0 : for(unsigned int i=0; i<size(); ++i) {
163 0 : if(neighbors_[i].first==index) neighbors.push_back(neighbors_[i].second);
164 0 : if(neighbors_[i].second==index) neighbors.push_back(neighbors_[i].first);
165 : }
166 0 : return neighbors;
167 : }
168 :
169 2523 : }
|
