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 "Pbc.h"
23 : #include "Tools.h"
24 : #include "Exception.h"
25 : #include "LatticeReduction.h"
26 : #include <iostream>
27 : #include "Random.h"
28 : #include <cmath>
29 :
30 : namespace PLMD {
31 :
32 26492 : Pbc::Pbc():
33 397380 : type(unset)
34 : {
35 26492 : box.zero();
36 26492 : invBox.zero();
37 26492 : }
38 :
39 20180 : void Pbc::buildShifts(gch::small_vector<Vector,maxshiftsize> shifts[2][2][2])const {
40 : const double small=1e-28;
41 :
42 : // clear all shifts
43 302700 : for(int i=0; i<2; i++) for(int j=0; j<2; j++) for(int k=0; k<2; k++) shifts[i][j][k].clear();
44 :
45 : // enumerate all possible shifts
46 : // since box is reduced, only 27 shifts have to be attempted
47 807200 : for(int l=-1; l<=1; l++) for(int m=-1; m<=1; m++) for(int n=-1; n<=1; n++) {
48 :
49 : // int/double shift vectors
50 544860 : const int ishift[3]= {l,m,n};
51 544860 : Vector dshift(l,m,n);
52 :
53 : // count how many components are != 0
54 : unsigned count=0;
55 2179440 : for(int s=0; s<3; s++) if(ishift[s]!=0) count++;
56 :
57 : // skips trivial (0,0,0) and cases with three shifts
58 : // only 18 shifts survive past this point
59 624164 : if(count==0 || count==3) continue;
60 :
61 : // check if that Wigner-Seitz face is perpendicular to the axis.
62 : // this allows to eliminate shifts in symmetric cells.
63 : // e.g., if one lactice vector is orthogonal to the plane spanned
64 : // by the other two vectors, that shift should never be tried
65 363240 : Vector cosdir=matmul(reduced,transpose(reduced),dshift);
66 363240 : double dp=dotProduct(dshift,cosdir);
67 363240 : double ref=modulo2(dshift)*modulo2(cosdir);
68 363240 : if(std::fabs(ref-dp*dp)<small) continue;
69 :
70 : // here we start pruning depending on the sign of the scaled coordinate
71 4259040 : for(int i=0; i<2; i++) for(int j=0; j<2; j++) for(int k=0; k<2; k++) {
72 :
73 2271488 : const int block[3]= {2*i-1,2*j-1,2*k-1};
74 :
75 : // skip cases where shift would bring too far from origin
76 : bool skip=false;
77 9085952 : for(int s=0; s<3; s++) if(ishift[s]*block[s]>0) skip=true;
78 2566066 : if(skip) continue;
79 : skip=true;
80 3103824 : for(int s=0; s<3; s++) {
81 : // check that the components of cosdir along the non-shifted directions
82 : // have the proper sign
83 2327868 : if(((1-ishift[s]*ishift[s])*block[s])*cosdir[s]<-small) skip=false;
84 : }
85 775956 : if(skip)continue;
86 :
87 : // if we arrive to this point, shift is eligible and is added to the list
88 962756 : shifts[i][j][k].push_back(matmul(transpose(reduced),dshift));
89 : }
90 : }
91 20180 : }
92 :
93 600000 : void Pbc::fullSearch(Vector&d)const {
94 600000 : if(type==unset) return;
95 528100 : Vector s=matmul(invReduced.transpose(),d);
96 2112400 : for(int i=0; i<3; i++) s[i]=Tools::pbc(s[i]);
97 528100 : d=matmul(reduced.transpose(),s);
98 : const int smax=4;
99 528100 : Vector a0(reduced.getRow(0));
100 528100 : Vector a1(reduced.getRow(1));
101 528100 : Vector a2(reduced.getRow(2));
102 528100 : Vector best(d);
103 528100 : double lbest=d.modulo2();
104 433042000 : for(int i=-smax; i<=smax; i++) for(int j=-smax; j<=smax; j++) for(int k=-smax; k<=smax; k++) {
105 384984900 : Vector trial=d+i*a0+j*a1+k*a2;
106 384984900 : double ltrial=trial.modulo2();
107 384984900 : if(ltrial<lbest) {
108 29897 : best=trial;
109 : lbest=ltrial;
110 : }
111 : }
112 528100 : d=best;
113 : }
114 :
115 90932 : void Pbc::setBox(const Tensor&b) {
116 90932 : box=b;
117 : // detect type:
118 : const double epsilon=1e-28;
119 :
120 90932 : type=unset;
121 90932 : double det=box.determinant();
122 90932 : if(det*det<epsilon) return;
123 :
124 : bool cxy=false;
125 : bool cxz=false;
126 : bool cyz=false;
127 88887 : if(box(0,1)*box(0,1)<epsilon && box(1,0)*box(1,0)<epsilon) cxy=true;
128 88887 : if(box(0,2)*box(0,2)<epsilon && box(2,0)*box(2,0)<epsilon) cxz=true;
129 88887 : if(box(1,2)*box(1,2)<epsilon && box(2,1)*box(2,1)<epsilon) cyz=true;
130 :
131 88887 : invBox=box.inverse();
132 :
133 88887 : if(cxy && cxz && cyz) type=orthorombic;
134 20180 : else type=generic;
135 :
136 88887 : if(type==orthorombic) {
137 68707 : reduced=box;
138 68707 : invReduced=inverse(reduced);
139 274828 : for(unsigned i=0; i<3; i++) {
140 206121 : diag[i]=box[i][i];
141 206121 : hdiag[i]=0.5*box[i][i];
142 206121 : mdiag[i]=-0.5*box[i][i];
143 : }
144 : } else {
145 20180 : reduced=box;
146 20180 : LatticeReduction::reduce(reduced);
147 20180 : invReduced=inverse(reduced);
148 20180 : buildShifts(shifts);
149 : }
150 :
151 : }
152 :
153 0 : double Pbc::distance( const bool pbc, const Vector& v1, const Vector& v2 ) const {
154 0 : if(pbc) { return ( distance(v1,v2) ).modulo(); }
155 0 : else { return ( delta(v1,v2) ).modulo(); }
156 : }
157 :
158 405876 : void Pbc::apply(std::vector<Vector>& dlist, unsigned max_index) const {
159 405876 : apply(VectorView(&dlist[0][0],dlist.size()),max_index);
160 405876 : }
161 :
162 405876 : void Pbc::apply(VectorView dlist, unsigned max_index) const {
163 405876 : if (max_index==0) max_index=dlist.size();
164 405876 : if(type==unset) {
165 : // do nothing
166 405872 : } else if(type==orthorombic) {
167 : #ifdef __PLUMED_PBC_WHILE
168 : for(unsigned k=0; k<max_index; ++k) {
169 : while(dlist[k][0]>hdiag[0]) dlist[k][0]-=diag[0];
170 : while(dlist[k][0]<=mdiag[0]) dlist[k][0]+=diag[0];
171 : while(dlist[k][1]>hdiag[1]) dlist[k][1]-=diag[1];
172 : while(dlist[k][1]<=mdiag[1]) dlist[k][1]+=diag[1];
173 : while(dlist[k][2]>hdiag[2]) dlist[k][2]-=diag[2];
174 : while(dlist[k][2]<=mdiag[2]) dlist[k][2]+=diag[2];
175 : }
176 : #else
177 26205115 : for(unsigned k=0; k<max_index; ++k) {
178 103197924 : for(int i=0; i<3; i++) {
179 77398443 : dlist[k][i]=Tools::pbc(dlist[k][i]*invBox(i,i))*box(i,i);
180 : }
181 : }
182 : #endif
183 238 : } else if(type==generic) {
184 21116 : for(unsigned k=0; k<max_index; ++k) {
185 : //Inlining by hand this part of function from distance speeds up by about 20-30%
186 : //against the previos version, and 60-80% agains this version non inlined.
187 : //I do not think is the `if(nshifts) *nshifts+=myshifts.size();`,
188 : //but that the compiler now see how we are juggling with the memory and it
189 : //does its magic
190 :
191 : //I tried writing VectorGeneric<3> matmul(const MemoryView<3UL> a,const TensorGeneric<3,3>&b)
192 : // by copy-pasting the original vector-tensor, but slows down this method by 10%... (on gcc9)
193 20878 : Vector s=matmul(Vector{dlist[k][0],dlist[k][1],dlist[k][2]},invReduced);
194 : // bring to -0.5,+0.5 region in scaled coordinates:
195 83512 : for(int i=0; i<3; ++i) {
196 62634 : s[i]=Tools::pbc(s[i]);
197 : }
198 20878 : Vector best(matmul(s,reduced));
199 : // check if shifts have to be attempted:
200 20878 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
201 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
202 34605 : const auto & myshifts(shifts[(s[0]>0?1:0)][(s[1]>0?1:0)][(s[2]>0?1:0)]);
203 17298 : Vector reference = best;
204 17298 : double lbest(modulo2(best));
205 : // loop over possible shifts:
206 81547 : for(unsigned i=0; i<myshifts.size(); ++i) {
207 64249 : Vector trial=reference+myshifts[i];
208 64249 : double ltrial=modulo2(trial);
209 64249 : if(ltrial<lbest) {
210 : lbest=ltrial;
211 1433 : best=trial;
212 : }
213 : }
214 : }
215 20878 : dlist[k][0] = best[0];
216 20878 : dlist[k][1] = best[1];
217 20878 : dlist[k][2] = best[2];
218 : }
219 0 : } else plumed_merror("unknown pbc type");
220 405876 : }
221 :
222 558983578 : Vector Pbc::distance(const Vector&v1,const Vector&v2,int*nshifts)const {
223 558983578 : Vector d=delta(v1,v2);
224 558983578 : if(type==unset) {
225 : // do nothing
226 532632348 : } else if(type==orthorombic) {
227 : #ifdef __PLUMED_PBC_WHILE
228 : for(unsigned i=0; i<3; i++) {
229 : while(d[i]>hdiag[i]) d[i]-=diag[i];
230 : while(d[i]<=mdiag[i]) d[i]+=diag[i];
231 : }
232 : #else
233 1937607024 : for(int i=0; i<3; i++) d[i]=Tools::pbc(d[i]*invBox(i,i))*box(i,i);
234 : #endif
235 48230592 : } else if(type==generic) {
236 48230592 : Vector s=matmul(d,invReduced);
237 : // check if images have to be computed:
238 : // if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)){
239 : // NOTICE: the check in the previous line, albeit correct, is breaking many regtest
240 : // since it does not apply Tools::pbc in many cases. Moreover, it does not
241 : // introduce a significant gain. I thus leave it out for the moment.
242 : if constexpr (true) {
243 : // bring to -0.5,+0.5 region in scaled coordinates:
244 192922368 : for(int i=0; i<3; i++) {
245 144691776 : s[i]=Tools::pbc(s[i]);
246 : }
247 48230592 : d=matmul(s,reduced);
248 : // check if shifts have to be attempted:
249 48230592 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
250 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
251 78664302 : const auto & myshifts(shifts[(s[0]>0?1:0)][(s[1]>0?1:0)][(s[2]>0?1:0)]);
252 39816653 : Vector best(d);
253 39816653 : double lbest(modulo2(best));
254 : // loop over possible shifts:
255 39816653 : if(nshifts) *nshifts+=myshifts.size();
256 179244081 : for(unsigned i=0; i<myshifts.size(); i++) {
257 139427428 : Vector trial=d+myshifts[i];
258 139427428 : double ltrial=modulo2(trial);
259 139427428 : if(ltrial<lbest) {
260 : lbest=ltrial;
261 3759325 : best=trial;
262 : }
263 : }
264 39816653 : d=best;
265 : }
266 : }
267 0 : } else plumed_merror("unknown pbc type");
268 558983578 : return d;
269 : }
270 :
271 1958610 : Vector Pbc::realToScaled(const Vector&d)const {
272 1958610 : return matmul(invBox.transpose(),d);
273 : }
274 :
275 201929 : Vector Pbc::scaledToReal(const Vector&d)const {
276 201929 : return matmul(box.transpose(),d);
277 : }
278 :
279 8569 : bool Pbc::isOrthorombic()const {
280 8569 : return type==orthorombic;
281 : }
282 :
283 28686 : const Tensor& Pbc::getBox()const {
284 28686 : return box;
285 : }
286 :
287 45429 : const Tensor& Pbc::getInvBox()const {
288 45429 : return invBox;
289 : }
290 :
291 : }
|
