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 "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 9280 : Pbc::Pbc():
33 9280 : type(unset)
34 : {
35 9280 : box.zero();
36 9280 : invBox.zero();
37 9280 : }
38 :
39 7540 : void Pbc::buildShifts(std::vector<Vector> shifts[2][2][2])const {
40 7540 : const double small=1e-28;
41 :
42 : // clear all shifts
43 7540 : 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 211120 : 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 203580 : int ishift[3]= {l,m,n};
51 203580 : Vector dshift(l,m,n);
52 :
53 : // count how many components are != 0
54 203580 : unsigned count=0;
55 203580 : 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 327988 : 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 135720 : Vector cosdir=matmul(reduced,transpose(reduced),dshift);
66 135720 : double dp=dotProduct(dshift,cosdir);
67 135720 : double ref=modulo2(dshift)*modulo2(cosdir);
68 135720 : if(std::fabs(ref-dp*dp)<small) continue;
69 :
70 : // here we start pruning depending on the sign of the scaled coordinate
71 712548 : for(int i=0; i<2; i++) for(int j=0; j<2; j++) for(int k=0; k<2; k++) {
72 :
73 633376 : 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 633376 : bool skip=false;
77 633376 : for(int s=0; s<3; s++) if(ishift[s]*block[s]>0) skip=true;
78 1140826 : if(skip) continue;
79 214916 : skip=true;
80 859664 : 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 644748 : if(((1-ishift[s]*ishift[s])*block[s])*cosdir[s]<-small) skip=false;
84 : }
85 214916 : if(skip)continue;
86 :
87 : // if we arrive to this point, shift is eligible and is added to the list
88 125926 : shifts[i][j][k].push_back(matmul(transpose(reduced),dshift));
89 : }
90 : }
91 7540 : }
92 :
93 600000 : void Pbc::fullSearch(Vector&d)const {
94 1200000 : if(type==unset) return;
95 528100 : Vector s=matmul(invReduced.transpose(),d);
96 528100 : for(int i=0; i<3; i++) s[i]=Tools::pbc(s[i]);
97 528100 : d=matmul(reduced.transpose(),s);
98 528100 : 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 385513000 : 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 29897 : lbest=ltrial;
110 : }
111 : }
112 528100 : d=best;
113 : }
114 :
115 34808 : void Pbc::setBox(const Tensor&b) {
116 34808 : box=b;
117 : // detect type:
118 34808 : const double epsilon=1e-28;
119 :
120 34808 : type=unset;
121 34808 : double det=box.determinant();
122 69616 : if(det*det<epsilon) return;
123 :
124 33714 : bool cxy=false;
125 33714 : bool cxz=false;
126 33714 : bool cyz=false;
127 33714 : if(box(0,1)*box(0,1)<epsilon && box(1,0)*box(1,0)<epsilon) cxy=true;
128 33714 : if(box(0,2)*box(0,2)<epsilon && box(2,0)*box(2,0)<epsilon) cxz=true;
129 33714 : if(box(1,2)*box(1,2)<epsilon && box(2,1)*box(2,1)<epsilon) cyz=true;
130 :
131 33714 : invBox=box.inverse();
132 :
133 33714 : if(cxy && cxz && cyz) type=orthorombic;
134 7540 : else type=generic;
135 :
136 33714 : if(type==orthorombic) {
137 26174 : reduced=box;
138 26174 : invReduced=inverse(reduced);
139 104696 : for(unsigned i=0; i<3; i++) {
140 78522 : diag[i]=box[i][i];
141 78522 : hdiag[i]=0.5*box[i][i];
142 78522 : mdiag[i]=-0.5*box[i][i];
143 : }
144 : } else {
145 7540 : reduced=box;
146 7540 : LatticeReduction::reduce(reduced);
147 7540 : invReduced=inverse(reduced);
148 7540 : 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 191388 : void Pbc::apply(std::vector<Vector>& dlist, unsigned max_index) const {
159 191388 : if (max_index==0) max_index=dlist.size();
160 191388 : if(type==unset) {
161 192690 : } else if(type==orthorombic) {
162 : #ifdef __PLUMED_PBC_WHILE
163 : for(unsigned k=0; k<max_index; ++k) {
164 : while(dlist[k][0]>hdiag[0]) dlist[k][0]-=diag[0];
165 : while(dlist[k][0]<=mdiag[0]) dlist[k][0]+=diag[0];
166 : while(dlist[k][1]>hdiag[1]) dlist[k][1]-=diag[1];
167 : while(dlist[k][1]<=mdiag[1]) dlist[k][1]+=diag[1];
168 : while(dlist[k][2]>hdiag[2]) dlist[k][2]-=diag[2];
169 : while(dlist[k][2]<=mdiag[2]) dlist[k][2]+=diag[2];
170 : }
171 : #else
172 188762 : for(unsigned k=0; k<max_index; ++k) for(int i=0; i<3; i++) dlist[k][i]=Tools::pbc(dlist[k][i]*invBox(i,i))*box(i,i);
173 : #endif
174 3928 : } else if(type==generic) {
175 3928 : for(unsigned k=0; k<max_index; ++k) dlist[k]=distance(Vector(0.0,0.0,0.0),dlist[k]);
176 0 : } else plumed_merror("unknown pbc type");
177 190109 : }
178 :
179 17286620 : Vector Pbc::distance(const Vector&v1,const Vector&v2,int*nshifts)const {
180 17286620 : Vector d=delta(v1,v2);
181 17322401 : if(type==unset) {
182 17237281 : } else if(type==orthorombic) {
183 : #ifdef __PLUMED_PBC_WHILE
184 : for(unsigned i=0; i<3; i++) {
185 : while(d[i]>hdiag[i]) d[i]-=diag[i];
186 : while(d[i]<=mdiag[i]) d[i]+=diag[i];
187 : }
188 : #else
189 12272103 : for(int i=0; i<3; i++) d[i]=Tools::pbc(d[i]*invBox(i,i))*box(i,i);
190 : #endif
191 4965178 : } else if(type==generic) {
192 4965178 : Vector s=matmul(d,invReduced);
193 : // check if images have to be computed:
194 : // if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)){
195 : // NOTICE: the check in the previous line, albeit correct, is breaking many regtest
196 : // since it does not apply Tools::pbc in many cases. Moreover, it does not
197 : // introduce a significant gain. I thus leave it out for the moment.
198 : if(true) {
199 : // bring to -0.5,+0.5 region in scaled coordinates:
200 4953366 : for(int i=0; i<3; i++) s[i]=Tools::pbc(s[i]);
201 4964950 : d=matmul(s,reduced);
202 : // check if shifts have to be attempted:
203 4964343 : if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)) {
204 : // list of shifts is specific for that "octant" (depends on signs of s[i]):
205 3783366 : const std::vector<Vector> & myshifts(shifts[(s[0]>0?1:0)][(s[1]>0?1:0)][(s[2]>0?1:0)]);
206 3783430 : Vector best(d);
207 3783430 : double lbest(modulo2(best));
208 : // loop over possible shifts:
209 3783293 : if(nshifts) *nshifts+=myshifts.size();
210 14322612 : for(unsigned i=0; i<myshifts.size(); i++) {
211 10554490 : Vector trial=d+myshifts[i];
212 10553211 : double ltrial=modulo2(trial);
213 10539319 : if(ltrial<lbest) {
214 738446 : lbest=ltrial;
215 738446 : best=trial;
216 : }
217 : }
218 3782482 : d=best;
219 : }
220 : }
221 0 : } else plumed_merror("unknown pbc type");
222 17300198 : return d;
223 : }
224 :
225 636232 : Vector Pbc::realToScaled(const Vector&d)const {
226 636232 : return matmul(invBox.transpose(),d);
227 : }
228 :
229 152079 : Vector Pbc::scaledToReal(const Vector&d)const {
230 152079 : return matmul(box.transpose(),d);
231 : }
232 :
233 338 : bool Pbc::isOrthorombic()const {
234 338 : return type==orthorombic;
235 : }
236 :
237 21407 : const Tensor& Pbc::getBox()const {
238 21407 : return box;
239 : }
240 :
241 1709 : const Tensor& Pbc::getInvBox()const {
242 1709 : return invBox;
243 : }
244 :
245 0 : void Pbc::test() {
246 0 : Random r;
247 0 : r.setSeed(-20);
248 0 : for(int i=0; i<1000; i++) {
249 : // random matrix with some zero element
250 0 : Tensor box;
251 0 : for(int j=0; j<3; j++) for(int k=0; k<3; k++) if(r.U01()>0.2) {
252 0 : box[j][k]=2.0*r.U01()-1.0;
253 : }
254 0 : int boxtype=i%10;
255 0 : switch(boxtype) {
256 : case 0:
257 : // cubic
258 0 : for(int j=0; j<3; j++) for(int k=0; k<3; k++) if(j!=k) box[j][k]=0.0;
259 0 : for(int j=1; j<3; j++) box[j][j]=box[0][0];
260 0 : break;
261 : case 1:
262 : // orthorombic
263 0 : for(int j=0; j<3; j++) for(int k=0; k<3; k++) if(j!=k) box[j][k]=0.0;
264 0 : break;
265 : case 2:
266 : // hexagonal
267 : {
268 0 : int perm=r.U01()*100;
269 0 : Vector a;
270 0 : a(0)=r.U01()*2-2; a(1)=0.0; a(2)=0.0;
271 0 : double d=r.U01()*2-2;
272 0 : Vector b(0.0,d,0.0);
273 0 : Vector c(0.0,0.5*d,sqrt(3.0)*d*0.5);
274 0 : box.setRow((perm+0)%3,a);
275 0 : box.setRow((perm+1)%3,b);
276 0 : box.setRow((perm+2)%3,c);
277 : }
278 0 : break;
279 : case 3:
280 : // bcc
281 : {
282 0 : int perm=r.U01()*100;
283 0 : double d=r.U01()*2-2;
284 0 : Vector a(d,d,d);
285 0 : Vector b(d,-d,d);
286 0 : Vector c(d,d,-d);
287 0 : box.setRow((perm+0)%3,a);
288 0 : box.setRow((perm+1)%3,b);
289 0 : box.setRow((perm+2)%3,c);
290 : }
291 0 : break;
292 : case 4:
293 : // fcc
294 : {
295 0 : int perm=r.U01()*100;
296 0 : double d=r.U01()*2-2;
297 0 : Vector a(d,d,0);
298 0 : Vector b(d,0,d);
299 0 : Vector c(0,d,d);
300 0 : box.setRow((perm+0)%3,a);
301 0 : box.setRow((perm+1)%3,b);
302 0 : box.setRow((perm+2)%3,c);
303 : }
304 0 : break;
305 : default:
306 : // triclinic
307 0 : break;
308 : }
309 :
310 0 : Pbc pbc;
311 0 : pbc.setBox(box);
312 0 : std::cerr<<"( "<<boxtype<<" )\n";
313 0 : std::cerr<<"Box:";
314 0 : for(int j=0; j<3; j++) for(int k=0; k<3; k++) std::cerr<<" "<<box[j][k];
315 0 : std::cerr<<"\n";
316 0 : std::cerr<<"Determinant: "<<determinant(box)<<"\n";
317 0 : std::cerr<<"Shifts:";
318 0 : for(int j=0; j<2; j++) for(int k=0; k<2; k++) for(int l=0; l<2; l++) std::cerr<<" "<<pbc.shifts[j][k][l].size();
319 0 : std::cerr<<"\n";
320 0 : int nshifts=0;
321 0 : int ntot=10000;
322 0 : for(int j=0; j<ntot; j++) {
323 0 : Vector v(r.U01()-0.5,r.U01()-0.5,r.U01()-0.5);
324 0 : v*=5;
325 0 : for(int j=0; j<3; j++) if(r.U01()>0.2) v(j)=0.0;
326 0 : Vector full(v);
327 0 : Vector fast=pbc.distance(Vector(0,0,0),v,&nshifts);
328 0 : full=fast;
329 :
330 0 : pbc.fullSearch(full);
331 :
332 0 : if(modulo2(fast-full)>1e-10) {
333 0 : std::cerr<<"orig "<<v[0]<<" "<<v[1]<<" "<<v[2]<<"\n";
334 0 : std::cerr<<"fast "<<fast[0]<<" "<<fast[1]<<" "<<fast[2]<<"\n";
335 0 : std::cerr<<"full "<<full[0]<<" "<<full[1]<<" "<<full[2]<<"\n";
336 0 : std::cerr<<"diff "<<modulo2(fast)-modulo2(full)<<std::endl;
337 0 : if(std::fabs(modulo2(fast)-modulo2(full))>1e-15) plumed_error();
338 : }
339 : }
340 0 : std::cerr<<"Average number of shifts: "<<double(nshifts)/double(ntot)<<"\n";
341 0 : }
342 0 : }
343 :
344 :
345 :
346 2523 : }
|
