1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13
14 Contributing author: David Nicholson (MIT)
15 -------------------------------------------------------------------------
16
17 This class contains functions to calculate the evolution of the periodic
18 simulation box under elongational flow as described by Matthew Dobson
19 in the arXiv preprint at https://arxiv.org/abs/1408.7078
20
21 Additionally, there are methods to do a lattice reduction to further
22 reduce the simulation box using the method of Igor Semaev at
23 https://link.springer.com/chapter/10.1007%2F3-540-44670-2_13
24 */
25
26 #include <cmath>
27 #include "uef_utils.h"
28
29 namespace LAMMPS_NS {
30 namespace UEF_utils{
31
UEFBox()32 UEFBox::UEFBox()
33 {
34
35 // initial box (also an inverse eigenvector matrix of automorphisms)
36
37 double x = 0.327985277605681;
38 double y = 0.591009048506103;
39 double z = 0.736976229099578;
40 l0[0][0]= z; l0[0][1]= y; l0[0][2]= x;
41 l0[1][0]=-x; l0[1][1]= z; l0[1][2]=-y;
42 l0[2][0]=-y; l0[2][1]= x; l0[2][2]= z;
43
44 // spectra of the two automorpisms (log of eigenvalues)
45
46 w1[0]=-1.177725211523360;
47 w1[1]=-0.441448620566067;
48 w1[2]= 1.619173832089425;
49 w2[0]= w1[1];
50 w2[1]= w1[2];
51 w2[2]= w1[0];
52
53 // initialize theta
54 // strain = w1 * theta1 + w2 * theta2
55
56 theta[0]=theta[1]=0;
57
58 //set up the initial box l and change of basis matrix r
59
60 for (int k=0;k<3;k++)
61 for (int j=0;j<3;j++) {
62 l[k][j] = l0[k][j];
63 r[j][k]=(j==k);
64 ri[j][k]=(j==k);
65 }
66
67 // get the initial rotation and upper triangular matrix
68
69 rotation_matrix(rot, lrot ,l);
70
71 // this is just a way to calculate the automorphisms
72 // themselves, which play a minor role in the calculations
73 // it's overkill, but only called once
74
75 double t1[3][3];
76 double t1i[3][3];
77 double t2[3][3];
78 double t2i[3][3];
79 double l0t[3][3];
80 for (int k=0; k<3; ++k)
81 for (int j=0; j<3; ++j) {
82 t1[k][j] = exp(w1[k])*l0[k][j];
83 t1i[k][j] = exp(-w1[k])*l0[k][j];
84 t2[k][j] = exp(w2[k])*l0[k][j];
85 t2i[k][j] = exp(-w2[k])*l0[k][j];
86 l0t[k][j] = l0[j][k];
87 }
88 mul_m2(l0t,t1);
89 mul_m2(l0t,t1i);
90 mul_m2(l0t,t2);
91 mul_m2(l0t,t2i);
92 for (int k=0; k<3; ++k)
93 for (int j=0; j<3; ++j) {
94 a1[k][j] = round(t1[k][j]);
95 a1i[k][j] = round(t1i[k][j]);
96 a2[k][j] = round(t2[k][j]);
97 a2i[k][j] = round(t2i[k][j]);
98 }
99
100 // winv used to transform between
101 // strain increments and theta increments
102
103 winv[0][0] = w2[1];
104 winv[0][1] = -w2[0];
105 winv[1][0] = -w1[1];
106 winv[1][1] = w1[0];
107 double d = w1[0]*w2[1] - w2[0]*w1[1];
108 for (int k=0;k<2;k++)
109 for (int j=0;j<2;j++)
110 winv[k][j] /= d;
111 }
112
113 /* ----------------------------------------------------------------------
114 get volume-correct r basis in: basis*cbrt(vol) = q*r
115 ------------------------------------------------------------------------- */
get_box(double x[3][3],double v)116 void UEFBox::get_box(double x[3][3], double v)
117 {
118 v = cbrtf(v);
119 for (int k=0;k<3;k++)
120 for (int j=0;j<3;j++)
121 x[k][j] = lrot[k][j]*v;
122 }
123
124 /* ----------------------------------------------------------------------
125 get rotation matrix q in: basis = q*r
126 ------------------------------------------------------------------------- */
get_rot(double x[3][3])127 void UEFBox::get_rot(double x[3][3])
128 {
129 for (int k=0;k<3;k++)
130 for (int j=0;j<3;j++)
131 x[k][j]=rot[k][j];
132 }
133
134 /* ----------------------------------------------------------------------
135 get inverse change of basis matrix
136 ------------------------------------------------------------------------- */
get_inverse_cob(int x[3][3])137 void UEFBox::get_inverse_cob(int x[3][3])
138 {
139 for (int k=0;k<3;k++)
140 for (int j=0;j<3;j++)
141 x[k][j]=ri[k][j];
142 }
143
144 /* ----------------------------------------------------------------------
145 apply diagonal, incompressible deformation
146 ------------------------------------------------------------------------- */
step_deform(const double ex,const double ey)147 void UEFBox::step_deform(const double ex, const double ey)
148 {
149 // increment theta values used in the reduction
150
151 theta[0] +=winv[0][0]*ex + winv[0][1]*ey;
152 theta[1] +=winv[1][0]*ex + winv[1][1]*ey;
153
154 // deformation of the box. reduce() needs to be called regularly or
155 // calculation will become unstable
156
157 double eps[3];
158 eps[0]=ex; eps[1] = ey; eps[2] = -ex-ey;
159 for (int k=0;k<3;k++) {
160 eps[k] = exp(eps[k]);
161 l[k][0] = eps[k]*l[k][0];
162 l[k][1] = eps[k]*l[k][1];
163 l[k][2] = eps[k]*l[k][2];
164 }
165 rotation_matrix(rot,lrot, l);
166 }
167
168 /* ----------------------------------------------------------------------
169 reduce the current basis
170 ------------------------------------------------------------------------- */
reduce()171 bool UEFBox::reduce()
172 {
173 // determine how many times to apply the automorphisms and find new theta
174 // values
175
176 int f1 = round(theta[0]);
177 int f2 = round(theta[1]);
178 theta[0] -= f1;
179 theta[1] -= f2;
180
181 // store old change or basis matrix to determine if it changes
182
183 int r0[3][3];
184 for (int k=0;k<3;k++)
185 for (int j=0;j<3;j++)
186 r0[k][j]=r[k][j];
187
188 // this modifies the old change basis matrix to handle the case where the
189 // automorphism transforms the box but the reduced basis doesn't change
190 // (r0 should still equal r at the end)
191
192 if (f1 > 0) for (int k=0;k<f1;k++) mul_m2 (a1,r0);
193 if (f1 < 0) for (int k=0;k<-f1;k++) mul_m2 (a1i,r0);
194 if (f2 > 0) for (int k=0;k<f2;k++) mul_m2 (a2,r0);
195 if (f2 < 0) for (int k=0;k<-f2;k++) mul_m2 (a2i,r0);
196
197 // robust reduction to the box defined by Dobson
198
199 for (int k=0;k<3;k++) {
200 double eps = exp(theta[0]*w1[k]+theta[1]*w2[k]);
201 l[k][0] = eps*l0[k][0];
202 l[k][1] = eps*l0[k][1];
203 l[k][2] = eps*l0[k][2];
204 }
205
206 // further reduce the box using greedy reduction and check
207 // if it changed from the last step using the change of basis
208 // matrices r and r0
209
210 greedy(l,r,ri);
211
212 // multiplying the inverse by the old change of basis matrix gives
213 // the inverse of the transformation itself (should be identity if
214 // no reduction takes place). This is used for image flags only.
215
216 mul_m1(ri,r0);
217 rotation_matrix(rot,lrot, l);
218 return !mat_same(r,r0);
219 }
220
221 /* ----------------------------------------------------------------------
222 set the strain to a specific value
223 ------------------------------------------------------------------------- */
set_strain(const double ex,const double ey)224 void UEFBox::set_strain(const double ex, const double ey)
225 {
226 theta[0] = winv[0][0]*ex + winv[0][1]*ey;
227 theta[1] = winv[1][0]*ex + winv[1][1]*ey;
228 theta[0] -= round(theta[0]);
229 theta[1] -= round(theta[1]);
230
231 for (int k=0;k<3;k++) {
232 double eps = exp(theta[0]*w1[k]+theta[1]*w2[k]);
233 l[k][0] = eps*l0[k][0];
234 l[k][1] = eps*l0[k][1];
235 l[k][2] = eps*l0[k][2];
236 }
237 greedy(l,r,ri);
238 rotation_matrix(rot,lrot, l);
239 }
240
241 /* ----------------------------------------------------------------------
242 qr reduction using householder reflections
243 q*m = r. q is orthogonal. m is input matrix. r is upper triangular
244 ------------------------------------------------------------------------- */
rotation_matrix(double q[3][3],double r[3][3],const double m[3][3])245 void rotation_matrix(double q[3][3], double r[3][3], const double m[3][3])
246 {
247 for (int k=0;k<3;k++)
248 for (int j=0;j<3;j++)
249 r[k][j] = m[k][j];
250
251 double a = -sqrt(col_prod(r,0,0))*r[0][0]/fabs(r[0][0]);
252 double v[3];
253 v[0] = r[0][0]-a;
254 v[1] = r[1][0];
255 v[2] = r[2][0];
256 a = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
257 v[0] /= a; v[1] /= a; v[2] /= a;
258 double qt[3][3];
259 for (int k=0;k<3;k++)
260 for (int j=0;j<3;j++) {
261 qt[k][j] = (k==j) - 2*v[k]*v[j];
262 q[k][j]= qt[k][j];
263 }
264 mul_m2(qt,r);
265 a = -sqrt(r[1][1]*r[1][1] + r[2][1]*r[2][1])*r[1][1]/fabs(r[1][1]);
266 v[0] = 0;
267 v[1] = r[1][1] - a;
268 v[2] = r[2][1];
269 a = sqrt(v[1]*v[1]+v[2]*v[2]);
270 v[1] /= a;
271 v[2] /= a;
272 for (int k=0;k<3;k++)
273 for (int j=0;j<3;j++)
274 qt[k][j] = (k==j) - 2*v[k]*v[j];
275 mul_m2(qt,r);
276 mul_m2(qt,q);
277
278 // this makes r have positive diagonals
279 // q*m = r <==> (-q)*m = (-r) will hold row-wise
280
281 if (r[0][0] < 0) { neg_row(q,0); neg_row(r,0); }
282 if (r[1][1] < 0) { neg_row(q,1); neg_row(r,1); }
283 if (r[2][2] < 0) { neg_row(q,2); neg_row(r,2); }
284 }
285
286 /* ----------------------------------------------------------------------
287 sort columns of b in order of increasing length
288 mimic column operations on ri and r
289 ------------------------------------------------------------------------- */
col_sort(double b[3][3],int r[3][3],int ri[3][3])290 void col_sort(double b[3][3],int r[3][3],int ri[3][3])
291 {
292 if (col_prod(b,0,0)>col_prod(b,1,1)) {
293 col_swap(b,0,1);
294 col_swap(r,0,1);
295 col_swap(ri,0,1);
296 }
297 if (col_prod(b,0,0)>col_prod(b,2,2)) {
298 col_swap(b,0,2);
299 col_swap(r,0,2);
300 col_swap(ri,0,2);
301 }
302 if (col_prod(b,1,1)>col_prod(b,2,2)) {
303 col_swap(b,1,2);
304 col_swap(r,1,2);
305 col_swap(ri,1,2);
306 }
307 }
308
309 /* ----------------------------------------------------------------------
310 1-2 reduction (Graham-Schmidt)
311 ------------------------------------------------------------------------- */
red12(double b[3][3],int r[3][3],int ri[3][3])312 void red12(double b[3][3],int r[3][3],int ri[3][3])
313 {
314 int y = round(col_prod(b,0,1)/col_prod(b,0,0));
315 b[0][1] -= y*b[0][0];
316 b[1][1] -= y*b[1][0];
317 b[2][1] -= y*b[2][0];
318
319 r[0][1] -= y*r[0][0];
320 r[1][1] -= y*r[1][0];
321 r[2][1] -= y*r[2][0];
322
323 ri[0][0] += y*ri[0][1];
324 ri[1][0] += y*ri[1][1];
325 ri[2][0] += y*ri[2][1];
326
327 if (col_prod(b,1,1) < col_prod(b,0,0)) {
328 col_swap(b,0,1);
329 col_swap(r,0,1);
330 col_swap(ri,0,1);
331 red12(b,r,ri);
332 }
333 }
334
335 /* ----------------------------------------------------------------------
336 Apply the Semaev condition for a 3-reduced basis
337 ------------------------------------------------------------------------- */
red3(double b[3][3],int r[3][3],int ri[3][3])338 void red3(double b[3][3], int r[3][3], int ri[3][3])
339 {
340 double b11 = col_prod(b,0,0);
341 double b22 = col_prod(b,1,1);
342 double b12 = col_prod(b,0,1);
343 double b13 = col_prod(b,0,2);
344 double b23 = col_prod(b,1,2);
345
346 double y2 =-(b23/b22-b12/b22*b13/b11)/(1-b12/b11*b12/b22);
347 double y1 =-(b13/b11-b12/b11*b23/b22)/(1-b12/b11*b12/b22);
348
349 int x1=0;
350 int x2=0;
351 double min = col_prod(b,2,2);
352 int x1v[2];
353 int x2v[2];
354 x1v[0] = floor(y1); x1v[1] = x1v[0]+1;
355 x2v[0] = floor(y2); x2v[1] = x2v[0]+1;
356 for (int k=0;k<2;k++)
357 for (int j=0;j<2;j++) {
358 double a[3];
359 a[0] = b[0][2] + x1v[k]*b[0][0] + x2v[j]*b[0][1];
360 a[1] = b[1][2] + x1v[k]*b[1][0] + x2v[j]*b[1][1];
361 a[2] = b[2][2] + x1v[k]*b[2][0] + x2v[j]*b[2][1];
362 double val=a[0]*a[0]+a[1]*a[1]+a[2]*a[2];
363 if (val<min) {
364 min = val;
365 x1 = x1v[k];
366 x2 = x2v[j];
367 }
368 }
369 if (x1 || x2) {
370 b[0][2] += x1*b[0][0] + x2*b[0][1];
371 b[1][2] += x1*b[1][0] + x2*b[1][1];
372 b[2][2] += x1*b[2][0] + x2*b[2][1];
373 r[0][2] += x1*r[0][0] + x2*r[0][1];
374 r[1][2] += x1*r[1][0] + x2*r[1][1];
375 r[2][2] += x1*r[2][0] + x2*r[2][1];
376 ri[0][0] += -x1*ri[0][2];
377 ri[1][0] += -x1*ri[1][2];
378 ri[2][0] += -x1*ri[2][2];
379 ri[0][1] += -x2*ri[0][2];
380 ri[1][1] += -x2*ri[1][2];
381 ri[2][1] += -x2*ri[2][2];
382 greedy_recurse(b,r,ri); // note the recursion step is here
383 }
384 }
385
386 /* ----------------------------------------------------------------------
387 the meat of the greedy reduction algorithm
388 ------------------------------------------------------------------------- */
greedy_recurse(double b[3][3],int r[3][3],int ri[3][3])389 void greedy_recurse(double b[3][3], int r[3][3], int ri[3][3])
390 {
391 col_sort(b,r,ri);
392 red12(b,r,ri);
393 red3(b,r,ri); // recursive caller
394 }
395
396 /* ----------------------------------------------------------------------
397 reduce the basis b. also output the change of basis matrix r and its
398 inverse ri
399 ------------------------------------------------------------------------- */
greedy(double b[3][3],int r[3][3],int ri[3][3])400 void greedy(double b[3][3],int r[3][3],int ri[3][3])
401 {
402 r[0][1]=r[0][2]=r[1][0]=r[1][2]=r[2][0]=r[2][1]=0;
403 r[0][0]=r[1][1]=r[2][2]=1;
404 ri[0][1]=ri[0][2]=ri[1][0]=ri[1][2]=ri[2][0]=ri[2][1]=0;
405 ri[0][0]=ri[1][1]=ri[2][2]=1;
406 greedy_recurse(b,r,ri);
407 make_unique(b,r,ri);
408 transpose(ri);
409 }
410
411 /* ----------------------------------------------------------------------
412 A reduced basis isn't unique. This procedure will make it
413 "more" unique. Degenerate cases are possible, but unlikely
414 with floating point math.
415 ------------------------------------------------------------------------- */
make_unique(double b[3][3],int r[3][3],int ri[3][3])416 void make_unique(double b[3][3], int r[3][3], int ri[3][3])
417 {
418 if (fabs(b[0][0]) < fabs(b[0][1])) {
419 col_swap(b,0,1);
420 col_swap(r,0,1);
421 col_swap(ri,0,1);
422 }
423 if (fabs(b[0][0]) < fabs(b[0][2])) {
424 col_swap(b,0,2);
425 col_swap(r,0,2);
426 col_swap(ri,0,2);
427 }
428 if (fabs(b[1][1]) < fabs(b[1][2])) {
429 col_swap(b,1,2);
430 col_swap(r,1,2);
431 col_swap(ri,1,2);
432 }
433
434 if (b[0][0] < 0) {
435 neg_col(b,0);
436 neg_col(r,0);
437 neg_col(ri,0);
438 }
439 if (b[1][1] < 0) {
440 neg_col(b,1);
441 neg_col(r,1);
442 neg_col(ri,1);
443 }
444 if (det(b) < 0) {
445 neg_col(b,2);
446 neg_col(r,2);
447 neg_col(ri,2);
448 }
449 }
450 }}
451