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 <cmath>
27 :
28 : namespace PLMD {
29 :
30 28418 : Pbc::Pbc():
31 28418 : type(unset),
32 28418 : nShiftsPerOctant{} {
33 : box.zero();
34 : invBox.zero();
35 28418 : }
36 :
37 0 : void Pbc::toACCDevice()const {
38 : // all the data (including the reduced box and the shifts used for generic
39 : // boxes) is stored in fixed-size members, so copying the object is enough
40 : #pragma acc enter data copyin(this[0:1], type)
41 : box.toACCDevice();
42 : invBox.toACCDevice();
43 0 : }
44 0 : void Pbc::removeFromACCDevice() const {
45 : invBox.removeFromACCDevice();
46 : box.removeFromACCDevice();
47 : // just delete
48 : #pragma acc exit data delete(type, this[0:1])
49 0 : }
50 :
51 : // Pbc::Pbc(const Pbc& other):
52 : // type(other.type),
53 : // box(other.box),
54 : // invBox(other.invBox),
55 : // reduced(other.reduced),
56 : // invReduced(other.invReduced)
57 : // {}
58 :
59 22626 : void Pbc::buildShifts() {
60 : const double small=1e-28;
61 :
62 : // clear all shifts
63 67878 : for(int i=0; i<2; i++)
64 135756 : for(int j=0; j<2; j++)
65 271512 : for(int k=0; k<2; k++) {
66 181008 : nShiftsPerOctant[i][j][k]=0;
67 : }
68 :
69 : // enumerate all possible shifts
70 : // since box is reduced, only 27 shifts have to be attempted
71 90504 : for(int l=-1; l<=1; l++)
72 271512 : for(int m=-1; m<=1; m++)
73 814536 : for(int n=-1; n<=1; n++) {
74 :
75 : // int/double shift vectors
76 610902 : const int ishift[3]= {l,m,n};
77 610902 : Vector dshift(l,m,n);
78 :
79 : // count how many components are != 0
80 : unsigned count=0;
81 2443608 : for(int s=0; s<3; s++)
82 1832706 : if(ishift[s]!=0) {
83 1221804 : count++;
84 : }
85 :
86 : // skips trivial (0,0,0) and cases with three shifts
87 : // only 18 shifts survive past this point
88 610902 : if(count==0 || count==3) {
89 301744 : continue;
90 : }
91 :
92 : // check if that Wigner-Seitz face is perpendicular to the axis.
93 : // this allows to eliminate shifts in symmetric cells.
94 : // e.g., if one lactice vector is orthogonal to the plane spanned
95 : // by the other two vectors, that shift should never be tried
96 407268 : Vector cosdir=matmul(reduced,transpose(reduced),dshift);
97 : double dp=dotProduct(dshift,cosdir);
98 407268 : double ref=modulo2(dshift)*modulo2(cosdir);
99 407268 : if(std::fabs(ref-dp*dp)<small) {
100 98110 : continue;
101 : }
102 :
103 : // here we start pruning depending on the sign of the scaled coordinate
104 927474 : for(int i=0; i<2; i++)
105 1854948 : for(int j=0; j<2; j++)
106 3709896 : for(int k=0; k<2; k++) {
107 :
108 2473264 : const int block[3]= {2*i-1,2*j-1,2*k-1};
109 :
110 : // skip cases where shift would bring too far from origin
111 : bool skip=false;
112 9893056 : for(int s=0; s<3; s++)
113 7419792 : if(ishift[s]*block[s]>0) {
114 : skip=true;
115 : }
116 2473264 : if(skip) {
117 1949770 : continue;
118 : }
119 : skip=true;
120 3373744 : for(int s=0; s<3; s++) {
121 : // check that the components of cosdir along the non-shifted directions
122 : // have the proper sign
123 2530308 : if(((1-ishift[s]*ishift[s])*block[s])*cosdir[s]<-small) {
124 : skip=false;
125 : }
126 : }
127 843436 : if(skip) {
128 319942 : continue;
129 : }
130 :
131 : // if we arrive to this point, shift is eligible and is added to the list
132 523494 : plumed_assert(nShiftsPerOctant[i][j][k]<maxshiftsize) << "too many shifts for a generic box";
133 523494 : shifts[i][j][k][nShiftsPerOctant[i][j][k]++]=matmul(transpose(reduced),dshift);
134 : }
135 : }
136 22626 : }
137 :
138 600000 : void Pbc::fullSearch(Vector&d)const {
139 600000 : if(type==unset) {
140 71900 : return;
141 : }
142 528100 : Vector s=matmul(invReduced.transpose(),d);
143 2112400 : for(int i=0; i<3; i++) {
144 1584300 : s[i]=Tools::pbc(s[i]);
145 : }
146 528100 : d=matmul(reduced.transpose(),s);
147 : const int smax=4;
148 : Vector a0(reduced.getRow(0));
149 : Vector a1(reduced.getRow(1));
150 : Vector a2(reduced.getRow(2));
151 528100 : Vector best(d);
152 : double lbest=d.modulo2();
153 5281000 : for(int i=-smax; i<=smax; i++)
154 47529000 : for(int j=-smax; j<=smax; j++)
155 427761000 : for(int k=-smax; k<=smax; k++) {
156 : Vector trial=d+i*a0+j*a1+k*a2;
157 : double ltrial=trial.modulo2();
158 384984900 : if(ltrial<lbest) {
159 : best=trial;
160 : lbest=ltrial;
161 : }
162 : }
163 528100 : d=best;
164 : }
165 :
166 186602 : void Pbc::setBox(const Tensor&b) {
167 186602 : box=b;
168 : // detect type:
169 : constexpr double boxEpsilon=1e-28;
170 :
171 186602 : type=unset;
172 186602 : double det=box.determinant();
173 186602 : if(det*det<boxEpsilon) {
174 : return;
175 : }
176 :
177 : bool cxy=false;
178 : bool cxz=false;
179 : bool cyz=false;
180 182105 : if(box(0,1)*box(0,1)<boxEpsilon && box(1,0)*box(1,0)<boxEpsilon) {
181 : cxy=true;
182 : }
183 182105 : if(box(0,2)*box(0,2)<boxEpsilon && box(2,0)*box(2,0)<boxEpsilon) {
184 : cxz=true;
185 : }
186 182105 : if(box(1,2)*box(1,2)<boxEpsilon && box(2,1)*box(2,1)<boxEpsilon) {
187 : cyz=true;
188 : }
189 :
190 182105 : invBox=box.inverse();
191 :
192 182105 : if(cxy && cxz && cyz) {
193 159479 : type=orthorombic;
194 : } else {
195 22626 : type=generic;
196 : }
197 :
198 182105 : if(type==orthorombic) {
199 159479 : reduced=box;
200 159479 : invReduced=inverse(reduced);
201 637916 : for(unsigned i=0; i<3; i++) {
202 478437 : diag[i]=box[i][i];
203 478437 : hdiag[i]=0.5*box[i][i];
204 478437 : mdiag[i]=-0.5*box[i][i];
205 : }
206 : } else {
207 22626 : reduced=box;
208 22626 : LatticeReduction::reduce(reduced);
209 22626 : invReduced=inverse(reduced);
210 22626 : buildShifts();
211 : }
212 :
213 : }
214 :
215 0 : double Pbc::distance( const bool pbc, const Vector& v1, const Vector& v2 ) const {
216 0 : if(pbc) {
217 0 : return ( distance(v1,v2) ).modulo();
218 : } else {
219 0 : return ( delta(v1,v2) ).modulo();
220 : }
221 : }
222 :
223 5571 : void Pbc::apply(std::vector<Vector>& dlist, unsigned max_index) const {
224 5571 : apply(VectorView(&dlist[0][0],dlist.size()),max_index);
225 5571 : }
226 :
227 1123886 : void Pbc::apply(VectorView dlist, unsigned max_index) const {
228 1123886 : if (max_index==0) {
229 46356 : max_index=dlist.size();
230 : }
231 1123886 : if(type==unset) {
232 : // do nothing
233 1120740 : } else if(type==orthorombic) {
234 : #ifdef __PLUMED_PBC_WHILE
235 : for(unsigned k=0; k<max_index; ++k) {
236 : while(dlist[k][0]>hdiag[0]) {
237 : dlist[k][0]-=diag[0];
238 : }
239 : while(dlist[k][0]<=mdiag[0]) {
240 : dlist[k][0]+=diag[0];
241 : }
242 : while(dlist[k][1]>hdiag[1]) {
243 : dlist[k][1]-=diag[1];
244 : }
245 : while(dlist[k][1]<=mdiag[1]) {
246 : dlist[k][1]+=diag[1];
247 : }
248 : while(dlist[k][2]>hdiag[2]) {
249 : dlist[k][2]-=diag[2];
250 : }
251 : while(dlist[k][2]<=mdiag[2]) {
252 : dlist[k][2]+=diag[2];
253 : }
254 : }
255 : #else
256 : // Keep invariant box values independent of writes through the distance view.
257 1090851 : const Vector inverseDiagonal(invBox(0,0),invBox(1,1),invBox(2,2));
258 1090851 : const Vector diagonal(box(0,0),box(1,1),box(2,2));
259 : // VectorView has exactly three components per row. Express the fixed
260 : // stride directly so the compiler can prove that rows do not overlap.
261 : double* const values=dlist.data();
262 32500173 : for(unsigned k=0; k<max_index; ++k) {
263 125637288 : for(int i=0; i<3; i++) {
264 94227966 : double& value=values[std::size_t(k)*3+i];
265 94227966 : value=Tools::pbc(value*inverseDiagonal[i])*diagonal[i];
266 : }
267 : }
268 : #endif
269 29889 : } else if(type==generic) {
270 21632147 : for(unsigned k=0; k<max_index; ++k) {
271 : //Inlining by hand this part of function from distance speeds up by about 20-30%
272 : //against the previos version, and 60-80% agains this version non inlined.
273 : //I do not think is the `if(nshifts) *nshifts+=nmyshifts;`,
274 : //but that the compiler now see how we are juggling with the memory and it
275 : //does its magic
276 :
277 : //I tried writing VectorGeneric<3> matmul(const MemoryView<3UL> a,const TensorGeneric<3,3>&b)
278 : // by copy-pasting the original vector-tensor, but slows down this method by 10%... (on gcc9)
279 21602258 : Vector s=matmul(Vector{dlist[k][0],dlist[k][1],dlist[k][2]},invReduced);
280 : // bring to -0.5,+0.5 region in scaled coordinates:
281 86409032 : for(int i=0; i<3; ++i) {
282 64806774 : s[i]=Tools::pbc(s[i]);
283 : }
284 21602258 : Vector best(matmul(s,reduced));
285 : // check if shifts have to be attempted:
286 21602258 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
287 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
288 35705576 : const unsigned o0=(s[0]>0?1:0), o1=(s[1]>0?1:0), o2=(s[2]>0?1:0);
289 17619700 : const Vector* myshifts=shifts[o0][o1][o2];
290 17619700 : const unsigned nmyshifts=nShiftsPerOctant[o0][o1][o2];
291 : Vector reference = best;
292 : double lbest(modulo2(best));
293 : // loop over possible shifts:
294 72656308 : for(unsigned i=0; i<nmyshifts; ++i) {
295 55036608 : Vector trial=reference+myshifts[i];
296 : double ltrial=modulo2(trial);
297 55036608 : if(ltrial<lbest) {
298 : lbest=ltrial;
299 : best=trial;
300 : }
301 : }
302 : }
303 21602258 : dlist[k][0] = best[0];
304 21602258 : dlist[k][1] = best[1];
305 21602258 : dlist[k][2] = best[2];
306 : }
307 : }
308 : // throws are not compatible with GPUs
309 : // so the `plumed_merror("unknown pbc type");` must be anticipated by
310 : // a check before running the whole calculation
311 1123886 : }
312 :
313 600000 : Vector Pbc::distance(const Vector&v1,const Vector&v2,int*nshifts)const {
314 : Vector d=delta(v1,v2);
315 600000 : if(type==unset) {
316 : // do nothing
317 528100 : } else if(type==orthorombic) {
318 : #ifdef __PLUMED_PBC_WHILE
319 : for(unsigned i=0; i<3; i++) {
320 : while(d[i]>hdiag[i]) {
321 : d[i]-=diag[i];
322 : }
323 : while(d[i]<=mdiag[i]) {
324 : d[i]+=diag[i];
325 : }
326 : }
327 : #else
328 535200 : for(int i=0; i<3; i++) {
329 401400 : d[i]=Tools::pbc(d[i]*invBox(i,i))*box(i,i);
330 : }
331 : #endif
332 394300 : } else if(type==generic) {
333 394300 : Vector s=matmul(d,invReduced);
334 : // check if images have to be computed:
335 : // if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)){
336 : // NOTICE: the check in the previous line, albeit correct, is breaking many regtest
337 : // since it does not apply Tools::pbc in many cases. Moreover, it does not
338 : // introduce a significant gain. I thus leave it out for the moment.
339 : if constexpr (true) {
340 : // bring to -0.5,+0.5 region in scaled coordinates:
341 1577200 : for(int i=0; i<3; i++) {
342 1182900 : s[i]=Tools::pbc(s[i]);
343 : }
344 394300 : d=matmul(s,reduced);
345 : // check if shifts have to be attempted:
346 394300 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
347 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
348 238157 : const unsigned o0=(s[0]>0?1:0), o1=(s[1]>0?1:0), o2=(s[2]>0?1:0);
349 109939 : const Vector* myshifts=shifts[o0][o1][o2];
350 109939 : const unsigned nmyshifts=nShiftsPerOctant[o0][o1][o2];
351 109939 : Vector best(d);
352 : double lbest(modulo2(best));
353 : // loop over possible shifts:
354 109939 : if(nshifts) {
355 109939 : *nshifts+=nmyshifts;
356 : }
357 443770 : for(unsigned i=0; i<nmyshifts; i++) {
358 333831 : Vector trial=d+myshifts[i];
359 : double ltrial=modulo2(trial);
360 333831 : if(ltrial<lbest) {
361 : lbest=ltrial;
362 : best=trial;
363 : }
364 : }
365 109939 : d=best;
366 : }
367 : }
368 : } else {
369 0 : plumed_merror("unknown pbc type");
370 : }
371 600000 : return d;
372 : }
373 :
374 593321624 : Vector Pbc::distance(const Vector&v1,const Vector&v2)const {
375 : Vector d=delta(v1,v2);
376 593321624 : if(type==unset) {
377 : // do nothing
378 563901882 : } else if(type==orthorombic) {
379 : #ifdef __PLUMED_PBC_WHILE
380 : for(unsigned i=0; i<3; i++) {
381 : while(d[i]>hdiag[i]) {
382 : d[i]-=diag[i];
383 : }
384 : while(d[i]<=mdiag[i]) {
385 : d[i]+=diag[i];
386 : }
387 : }
388 : #else
389 2111850072 : for(int i=0; i<3; i++) {
390 1583887554 : d[i]=Tools::pbc(d[i]*invBox(i,i))*box(i,i);
391 : }
392 : #endif
393 35939364 : } else if(type==generic) {
394 35939364 : Vector s=matmul(d,invReduced);
395 : // check if images have to be computed:
396 : // if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)){
397 : // NOTICE: the check in the previous line, albeit correct, is breaking many regtest
398 : // since it does not apply Tools::pbc in many cases. Moreover, it does not
399 : // introduce a significant gain. I thus leave it out for the moment.
400 : if constexpr (true) {
401 : // bring to -0.5,+0.5 region in scaled coordinates:
402 143757456 : for(int i=0; i<3; i++) {
403 107818092 : s[i]=Tools::pbc(s[i]);
404 : }
405 35939364 : d=matmul(s,reduced);
406 : // check if shifts have to be attempted:
407 35939364 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
408 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
409 58758253 : const unsigned o0=(s[0]>0?1:0), o1=(s[1]>0?1:0), o2=(s[2]>0?1:0);
410 29507878 : const Vector* myshifts=shifts[o0][o1][o2];
411 29507878 : const unsigned nmyshifts=nShiftsPerOctant[o0][o1][o2];
412 29507878 : Vector best(d);
413 : double lbest(modulo2(best));
414 : // loop over possible shifts:
415 136630219 : for(unsigned i=0; i<nmyshifts; i++) {
416 107122341 : Vector trial=d+myshifts[i];
417 : double ltrial=modulo2(trial);
418 107122341 : if(ltrial<lbest) {
419 : lbest=ltrial;
420 : best=trial;
421 : }
422 : }
423 29507878 : d=best;
424 : }
425 : }
426 : // no throws on the GPU
427 : } //else plumed_merror("unknown pbc type");
428 593321624 : return d;
429 : }
430 :
431 100000 : Vector3f Pbc::distance(const Vector3f&v1,const Vector3f&v2)const {
432 100000 : Vector3f d=delta(v1,v2);
433 100000 : if(type==unset) {
434 : // do nothing
435 100000 : } else if(type==orthorombic) {
436 : #ifdef __PLUMED_PBC_WHILE
437 : for(unsigned i=0; i<3; i++) {
438 : while(d[i]>hdiag[i]) {
439 : d[i]-=diag[i];
440 : }
441 : while(d[i]<=mdiag[i]) {
442 : d[i]+=diag[i];
443 : }
444 : }
445 : #else
446 0 : for(int i=0; i<3; i++) {
447 0 : d[i]=Tools::pbc(d[i]*invBox(i,i))*box(i,i);
448 : }
449 : #endif
450 100000 : } else if(type==generic) {
451 100000 : Vector3f s=matmul(d,invReduced);
452 : // check if images have to be computed:
453 : // if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)){
454 : // NOTICE: the check in the previous line, albeit correct, is breaking many regtest
455 : // since it does not apply Tools::pbc in many cases. Moreover, it does not
456 : // introduce a significant gain. I thus leave it out for the moment.
457 : if constexpr (true) {
458 : // bring to -0.5,+0.5 region in scaled coordinates:
459 400000 : for(int i=0; i<3; i++) {
460 300000 : s[i]=Tools::pbc(s[i]);
461 : }
462 100000 : d=matmul(s,reduced);
463 : // check if shifts have to be attempted:
464 100000 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
465 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
466 166576 : const unsigned o0=(s[0]>0?1:0), o1=(s[1]>0?1:0), o2=(s[2]>0?1:0);
467 83177 : const Vector* myshifts=shifts[o0][o1][o2];
468 83177 : const unsigned nmyshifts=nShiftsPerOctant[o0][o1][o2];
469 83177 : Vector3f best(d);
470 83177 : double lbest(modulo2(best));
471 : // loop over possible shifts:
472 394535 : for(unsigned i=0; i<nmyshifts; i++) {
473 311358 : Vector3f trial=sumT(d,myshifts[i]);
474 311358 : double ltrial=modulo2(trial);
475 311358 : if(ltrial<lbest) {
476 : lbest=ltrial;
477 : best=trial;
478 : }
479 : }
480 83177 : d=best;
481 : }
482 : }
483 : // no throws on the GPU
484 : } //else plumed_merror("unknown pbc type");
485 100000 : return d;
486 : }
487 :
488 2711285 : Vector Pbc::realToScaled(const Vector&d)const {
489 2711285 : return matmul(invBox.transpose(),d);
490 : }
491 :
492 0 : Vector3f Pbc::realToScaled(const Vector3f&d)const {
493 0 : return matmul(invBox.transpose(),d).convert<float>();
494 : }
495 :
496 165488 : Vector Pbc::scaledToReal(const Vector&d)const {
497 165488 : return matmul(box.transpose(),d);
498 : }
499 0 : Vector3f Pbc::scaledToReal(const Vector3f&d)const {
500 0 : return matmul(box.transpose(),d).convert<float>();
501 : }
502 :
503 9361 : bool Pbc::isOrthorombic()const {
504 9361 : return type==orthorombic;
505 : }
506 :
507 43429 : const Tensor& Pbc::getBox()const {
508 43429 : return box;
509 : }
510 :
511 50581 : const Tensor& Pbc::getInvBox()const {
512 50581 : return invBox;
513 : }
514 :
515 : }
|
