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 26321 : Pbc::Pbc():
33 394815 : type(unset) {
34 26321 : box.zero();
35 26321 : invBox.zero();
36 26321 : }
37 :
38 20240 : void Pbc::buildShifts(gch::small_vector<Vector,maxshiftsize> shifts[2][2][2])const {
39 : const double small=1e-28;
40 :
41 : // clear all shifts
42 60720 : for(int i=0; i<2; i++)
43 121440 : for(int j=0; j<2; j++)
44 242880 : for(int k=0; k<2; k++) {
45 161920 : shifts[i][j][k].clear();
46 : }
47 :
48 : // enumerate all possible shifts
49 : // since box is reduced, only 27 shifts have to be attempted
50 80960 : for(int l=-1; l<=1; l++)
51 242880 : for(int m=-1; m<=1; m++)
52 728640 : for(int n=-1; n<=1; n++) {
53 :
54 : // int/double shift vectors
55 546480 : const int ishift[3]= {l,m,n};
56 546480 : Vector dshift(l,m,n);
57 :
58 : // count how many components are != 0
59 : unsigned count=0;
60 2185920 : for(int s=0; s<3; s++)
61 1639440 : if(ishift[s]!=0) {
62 1092960 : count++;
63 : }
64 :
65 : // skips trivial (0,0,0) and cases with three shifts
66 : // only 18 shifts survive past this point
67 546480 : if(count==0 || count==3) {
68 261680 : continue;
69 : }
70 :
71 : // check if that Wigner-Seitz face is perpendicular to the axis.
72 : // this allows to eliminate shifts in symmetric cells.
73 : // e.g., if one lactice vector is orthogonal to the plane spanned
74 : // by the other two vectors, that shift should never be tried
75 364320 : Vector cosdir=matmul(reduced,transpose(reduced),dshift);
76 364320 : double dp=dotProduct(dshift,cosdir);
77 364320 : double ref=modulo2(dshift)*modulo2(cosdir);
78 364320 : if(std::fabs(ref-dp*dp)<small) {
79 79520 : continue;
80 : }
81 :
82 : // here we start pruning depending on the sign of the scaled coordinate
83 854400 : for(int i=0; i<2; i++)
84 1708800 : for(int j=0; j<2; j++)
85 3417600 : for(int k=0; k<2; k++) {
86 :
87 2278400 : const int block[3]= {2*i-1,2*j-1,2*k-1};
88 :
89 : // skip cases where shift would bring too far from origin
90 : bool skip=false;
91 9113600 : for(int s=0; s<3; s++)
92 6835200 : if(ishift[s]*block[s]>0) {
93 : skip=true;
94 : }
95 2278400 : if(skip) {
96 1795510 : continue;
97 : }
98 : skip=true;
99 3113328 : for(int s=0; s<3; s++) {
100 : // check that the components of cosdir along the non-shifted directions
101 : // have the proper sign
102 2334996 : if(((1-ishift[s]*ishift[s])*block[s])*cosdir[s]<-small) {
103 : skip=false;
104 : }
105 : }
106 778332 : if(skip) {
107 295442 : continue;
108 : }
109 :
110 : // if we arrive to this point, shift is eligible and is added to the list
111 965780 : shifts[i][j][k].push_back(matmul(transpose(reduced),dshift));
112 : }
113 : }
114 20240 : }
115 :
116 600000 : void Pbc::fullSearch(Vector&d)const {
117 600000 : if(type==unset) {
118 71900 : return;
119 : }
120 528100 : Vector s=matmul(invReduced.transpose(),d);
121 2112400 : for(int i=0; i<3; i++) {
122 1584300 : s[i]=Tools::pbc(s[i]);
123 : }
124 528100 : d=matmul(reduced.transpose(),s);
125 : const int smax=4;
126 528100 : Vector a0(reduced.getRow(0));
127 528100 : Vector a1(reduced.getRow(1));
128 528100 : Vector a2(reduced.getRow(2));
129 528100 : Vector best(d);
130 528100 : double lbest=d.modulo2();
131 5281000 : for(int i=-smax; i<=smax; i++)
132 47529000 : for(int j=-smax; j<=smax; j++)
133 427761000 : for(int k=-smax; k<=smax; k++) {
134 384984900 : Vector trial=d+i*a0+j*a1+k*a2;
135 384984900 : double ltrial=trial.modulo2();
136 384984900 : if(ltrial<lbest) {
137 29897 : best=trial;
138 : lbest=ltrial;
139 : }
140 : }
141 528100 : d=best;
142 : }
143 :
144 93182 : void Pbc::setBox(const Tensor&b) {
145 93182 : box=b;
146 : // detect type:
147 : const double epsilon=1e-28;
148 :
149 93182 : type=unset;
150 93182 : double det=box.determinant();
151 93182 : if(det*det<epsilon) {
152 : return;
153 : }
154 :
155 : bool cxy=false;
156 : bool cxz=false;
157 : bool cyz=false;
158 91141 : if(box(0,1)*box(0,1)<epsilon && box(1,0)*box(1,0)<epsilon) {
159 : cxy=true;
160 : }
161 91141 : if(box(0,2)*box(0,2)<epsilon && box(2,0)*box(2,0)<epsilon) {
162 : cxz=true;
163 : }
164 91141 : if(box(1,2)*box(1,2)<epsilon && box(2,1)*box(2,1)<epsilon) {
165 : cyz=true;
166 : }
167 :
168 91141 : invBox=box.inverse();
169 :
170 91141 : if(cxy && cxz && cyz) {
171 70901 : type=orthorombic;
172 : } else {
173 20240 : type=generic;
174 : }
175 :
176 91141 : if(type==orthorombic) {
177 70901 : reduced=box;
178 70901 : invReduced=inverse(reduced);
179 283604 : for(unsigned i=0; i<3; i++) {
180 212703 : diag[i]=box[i][i];
181 212703 : hdiag[i]=0.5*box[i][i];
182 212703 : mdiag[i]=-0.5*box[i][i];
183 : }
184 : } else {
185 20240 : reduced=box;
186 20240 : LatticeReduction::reduce(reduced);
187 20240 : invReduced=inverse(reduced);
188 20240 : buildShifts(shifts);
189 : }
190 :
191 : }
192 :
193 0 : double Pbc::distance( const bool pbc, const Vector& v1, const Vector& v2 ) const {
194 0 : if(pbc) {
195 0 : return ( distance(v1,v2) ).modulo();
196 : } else {
197 0 : return ( delta(v1,v2) ).modulo();
198 : }
199 : }
200 :
201 405876 : void Pbc::apply(std::vector<Vector>& dlist, unsigned max_index) const {
202 405876 : apply(VectorView(&dlist[0][0],dlist.size()),max_index);
203 405876 : }
204 :
205 405876 : void Pbc::apply(VectorView dlist, unsigned max_index) const {
206 405876 : if (max_index==0) {
207 0 : max_index=dlist.size();
208 : }
209 405876 : if(type==unset) {
210 : // do nothing
211 405872 : } else if(type==orthorombic) {
212 : #ifdef __PLUMED_PBC_WHILE
213 : for(unsigned k=0; k<max_index; ++k) {
214 : while(dlist[k][0]>hdiag[0]) {
215 : dlist[k][0]-=diag[0];
216 : }
217 : while(dlist[k][0]<=mdiag[0]) {
218 : dlist[k][0]+=diag[0];
219 : }
220 : while(dlist[k][1]>hdiag[1]) {
221 : dlist[k][1]-=diag[1];
222 : }
223 : while(dlist[k][1]<=mdiag[1]) {
224 : dlist[k][1]+=diag[1];
225 : }
226 : while(dlist[k][2]>hdiag[2]) {
227 : dlist[k][2]-=diag[2];
228 : }
229 : while(dlist[k][2]<=mdiag[2]) {
230 : dlist[k][2]+=diag[2];
231 : }
232 : }
233 : #else
234 26205115 : for(unsigned k=0; k<max_index; ++k) {
235 103197924 : for(int i=0; i<3; i++) {
236 77398443 : dlist[k][i]=Tools::pbc(dlist[k][i]*invBox(i,i))*box(i,i);
237 : }
238 : }
239 : #endif
240 238 : } else if(type==generic) {
241 21116 : for(unsigned k=0; k<max_index; ++k) {
242 : //Inlining by hand this part of function from distance speeds up by about 20-30%
243 : //against the previos version, and 60-80% agains this version non inlined.
244 : //I do not think is the `if(nshifts) *nshifts+=myshifts.size();`,
245 : //but that the compiler now see how we are juggling with the memory and it
246 : //does its magic
247 :
248 : //I tried writing VectorGeneric<3> matmul(const MemoryView<3UL> a,const TensorGeneric<3,3>&b)
249 : // by copy-pasting the original vector-tensor, but slows down this method by 10%... (on gcc9)
250 20878 : Vector s=matmul(Vector{dlist[k][0],dlist[k][1],dlist[k][2]},invReduced);
251 : // bring to -0.5,+0.5 region in scaled coordinates:
252 83512 : for(int i=0; i<3; ++i) {
253 62634 : s[i]=Tools::pbc(s[i]);
254 : }
255 20878 : Vector best(matmul(s,reduced));
256 : // check if shifts have to be attempted:
257 20878 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
258 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
259 34605 : const auto & myshifts(shifts[(s[0]>0?1:0)][(s[1]>0?1:0)][(s[2]>0?1:0)]);
260 17298 : Vector reference = best;
261 17298 : double lbest(modulo2(best));
262 : // loop over possible shifts:
263 81547 : for(unsigned i=0; i<myshifts.size(); ++i) {
264 64249 : Vector trial=reference+myshifts[i];
265 64249 : double ltrial=modulo2(trial);
266 64249 : if(ltrial<lbest) {
267 : lbest=ltrial;
268 1433 : best=trial;
269 : }
270 : }
271 : }
272 20878 : dlist[k][0] = best[0];
273 20878 : dlist[k][1] = best[1];
274 20878 : dlist[k][2] = best[2];
275 : }
276 : } else {
277 0 : plumed_merror("unknown pbc type");
278 : }
279 405876 : }
280 :
281 559421392 : Vector Pbc::distance(const Vector&v1,const Vector&v2,int*nshifts)const {
282 559421392 : Vector d=delta(v1,v2);
283 559421392 : if(type==unset) {
284 : // do nothing
285 533228202 : } else if(type==orthorombic) {
286 : #ifdef __PLUMED_PBC_WHILE
287 : for(unsigned i=0; i<3; i++) {
288 : while(d[i]>hdiag[i]) {
289 : d[i]-=diag[i];
290 : }
291 : while(d[i]<=mdiag[i]) {
292 : d[i]+=diag[i];
293 : }
294 : }
295 : #else
296 1939966680 : for(int i=0; i<3; i++) {
297 1454975010 : d[i]=Tools::pbc(d[i]*invBox(i,i))*box(i,i);
298 : }
299 : #endif
300 48236532 : } else if(type==generic) {
301 48236532 : Vector s=matmul(d,invReduced);
302 : // check if images have to be computed:
303 : // if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)){
304 : // NOTICE: the check in the previous line, albeit correct, is breaking many regtest
305 : // since it does not apply Tools::pbc in many cases. Moreover, it does not
306 : // introduce a significant gain. I thus leave it out for the moment.
307 : if constexpr (true) {
308 : // bring to -0.5,+0.5 region in scaled coordinates:
309 192946128 : for(int i=0; i<3; i++) {
310 144709596 : s[i]=Tools::pbc(s[i]);
311 : }
312 48236532 : d=matmul(s,reduced);
313 : // check if shifts have to be attempted:
314 48236532 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
315 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
316 78664302 : const auto & myshifts(shifts[(s[0]>0?1:0)][(s[1]>0?1:0)][(s[2]>0?1:0)]);
317 39816653 : Vector best(d);
318 39816653 : double lbest(modulo2(best));
319 : // loop over possible shifts:
320 39816653 : if(nshifts) {
321 109939 : *nshifts+=myshifts.size();
322 : }
323 179244081 : for(unsigned i=0; i<myshifts.size(); i++) {
324 139427428 : Vector trial=d+myshifts[i];
325 139427428 : double ltrial=modulo2(trial);
326 139427428 : if(ltrial<lbest) {
327 : lbest=ltrial;
328 3759325 : best=trial;
329 : }
330 : }
331 39816653 : d=best;
332 : }
333 : }
334 : } else {
335 0 : plumed_merror("unknown pbc type");
336 : }
337 559421392 : return d;
338 : }
339 :
340 1974153 : Vector Pbc::realToScaled(const Vector&d)const {
341 1974153 : return matmul(invBox.transpose(),d);
342 : }
343 :
344 207347 : Vector Pbc::scaledToReal(const Vector&d)const {
345 207347 : return matmul(box.transpose(),d);
346 : }
347 :
348 8569 : bool Pbc::isOrthorombic()const {
349 8569 : return type==orthorombic;
350 : }
351 :
352 28812 : const Tensor& Pbc::getBox()const {
353 28812 : return box;
354 : }
355 :
356 45554 : const Tensor& Pbc::getInvBox()const {
357 45554 : return invBox;
358 : }
359 :
360 : }
|
