1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
8 // Luke Shulenburger, lshulen@sandia.gov, Sandia National Laboratories
9 // Anouar Benali, benali@anl.gov, Argonne National Laboratory
10 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
11 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
12 //
13 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
14 //////////////////////////////////////////////////////////////////////////////////////
15
16
17 #include <iostream>
18 #include <cmath>
19 #include <cstdlib>
20 #include <stdio.h>
21 #include <string>
22 #include <cstring>
23
24
25 template<typename T>
getDet(T * mat)26 T getDet(T* mat)
27 {
28 return mat[0] * (mat[4] * mat[8] - mat[7] * mat[5]) - mat[1] * (mat[3] * mat[8] - mat[5] * mat[6]) +
29 mat[2] * (mat[3] * mat[7] - mat[4] * mat[6]);
30 }
31
getSupercell(double * prim,int * tile,double * super)32 void getSupercell(double* prim, int* tile, double* super)
33 {
34 super[0] = tile[0] * prim[0] + tile[1] * prim[3] + tile[2] * prim[6];
35 super[1] = tile[0] * prim[1] + tile[1] * prim[4] + tile[2] * prim[7];
36 super[2] = tile[0] * prim[2] + tile[1] * prim[5] + tile[2] * prim[8];
37 super[3] = tile[3] * prim[0] + tile[4] * prim[3] + tile[5] * prim[6];
38 super[4] = tile[3] * prim[1] + tile[4] * prim[4] + tile[5] * prim[7];
39 super[5] = tile[3] * prim[2] + tile[4] * prim[5] + tile[5] * prim[8];
40 super[6] = tile[6] * prim[0] + tile[7] * prim[3] + tile[8] * prim[6];
41 super[7] = tile[6] * prim[1] + tile[7] * prim[4] + tile[8] * prim[7];
42 super[8] = tile[6] * prim[2] + tile[7] * prim[5] + tile[8] * prim[8];
43 }
44
45 template<typename T>
dot(T * a,T * b)46 T dot(T* a, T* b)
47 {
48 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
49 }
50
51 template<typename T>
cross(T * a,T * b,T * c)52 void cross(T* a, T* b, T* c)
53 {
54 c[0] = a[1] * b[2] - a[2] * b[1];
55 c[1] = -a[0] * b[2] + a[2] * b[0];
56 c[2] = a[0] * b[1] - a[1] * b[0];
57 }
58
SimCellRad(double * mat)59 double SimCellRad(double* mat)
60 {
61 double A[3];
62 double B[3];
63 double C[3];
64 double BxC[3];
65 double radius = 5000000000000000.0;
66 for (int i = 0; i < 3; i++)
67 {
68 const int astart = i * 3;
69 const int bstart = ((i + 1) % 3) * 3;
70 const int cstart = ((i + 2) % 3) * 3;
71 for (int j = 0; j < 3; j++)
72 {
73 A[j] = mat[astart + j];
74 B[j] = mat[bstart + j];
75 C[j] = mat[cstart + j];
76 }
77 cross(B, C, BxC);
78 double val = std::abs(0.5 * dot(A, BxC) / sqrt(dot(BxC, BxC)));
79 if (val < radius)
80 radius = val;
81 }
82 return radius;
83 }
84
getScore(double * mat)85 double getScore(double* mat)
86 {
87 double score = 0;
88 static const double tol = 0.001;
89 // Highest preference for diagonal matrices
90 const double abssumoffidag =
91 std::abs(mat[1]) + std::abs(mat[2]) + std::abs(mat[3]) + std::abs(mat[5]) + std::abs(mat[6]) + std::abs(mat[7]);
92 if (abssumoffidag < tol)
93 {
94 // std::cout << "Found a diagonal supercell!" << std::endl;
95 score += 50000;
96 }
97 // Now prefer positive elements that come as early as possible
98 for (int i = 0; i < 9; i++)
99 {
100 if (mat[i] > 0.0)
101 {
102 score += 10;
103 double v = (9.0 - static_cast<double>(i)) * 0.1;
104 score += v * v;
105 }
106 }
107 return score;
108 }
109
WigSeitzRad(double * mat)110 double WigSeitzRad(double* mat)
111 {
112 double rmin = 1000000000000000;
113 for (int i = -1; i <= 1; i++)
114 {
115 for (int j = -1; j <= 1; j++)
116 {
117 for (int k = -1; k <= 1; k++)
118 {
119 if ((i != 0) || (j != 0) || (k != 0))
120 {
121 double d[3];
122 d[0] = i * mat[0] + j * mat[3] + k * mat[6];
123 d[1] = i * mat[1] + j * mat[4] + k * mat[7];
124 d[2] = i * mat[2] + j * mat[5] + k * mat[8];
125 double dist = 0.5 * sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
126 if (dist < rmin)
127 rmin = dist;
128 }
129 }
130 }
131 }
132 return rmin;
133 }
134
135
getBestTile(double * primcell,int target,int * tilemat,double & radius,int range=7)136 void getBestTile(double* primcell, int target, int* tilemat, double& radius, int range = 7)
137 {
138 double largest = 0.0;
139 double bestScore = 0.0;
140 static const double tol = 0.0000001;
141 double detPrim = getDet(primcell);
142 if (detPrim < 0)
143 {
144 target *= -1;
145 }
146 #pragma omp parallel
147 {
148 double my_largest = 0.0;
149 int my_besttile[9];
150 double localBestScore = 0.0;
151 #pragma omp for
152 for (int i = -range; i <= range; i++)
153 {
154 int d[9];
155 double super[9];
156 d[0] = i;
157 for (int j = -range; j <= range; j++)
158 {
159 d[1] = j;
160 for (int k = -range; k <= range; k++)
161 {
162 d[2] = k;
163 for (int l = -range; l <= range; l++)
164 {
165 d[3] = l;
166 for (int m = -range; m <= range; m++)
167 {
168 d[4] = m;
169 int denominator = j * l - i * m;
170 for (int n = -range; n <= range; n++)
171 {
172 d[5] = n;
173 int fpp = k * l - i * n;
174 for (int o = -range; o <= range; o++)
175 {
176 d[6] = o;
177 int sp = o * (n * j - k * m);
178 for (int p = -range; p <= range; p++)
179 {
180 d[7] = p;
181 int numpart = p * fpp + sp;
182 if (denominator != 0)
183 {
184 int rem = 5;
185 rem = (numpart - target) % denominator;
186 if (rem == 0)
187 {
188 d[8] = (numpart - target) / denominator;
189 getSupercell(primcell, d, super);
190 double score = getScore(super);
191 double rad = SimCellRad(super);
192 //double rad = WigSeitzRad(super);
193 if (rad > my_largest + tol || (rad > my_largest - tol && score > localBestScore))
194 {
195 my_largest = rad;
196 localBestScore = score;
197 std::memcpy(my_besttile, d, 9 * sizeof(int));
198 }
199 }
200 }
201 else
202 // Handle case where denominator is 0
203 {
204 if (numpart == target)
205 {
206 // cout << "Got here" << std::endl;
207 for (int q = -range; q <= range; q++)
208 {
209 d[8] = q;
210 getSupercell(primcell, d, super);
211 double score = getScore(super);
212 double rad = SimCellRad(super);
213 //double rad = WigSeitzRad(super);
214 if (rad > my_largest + tol || (rad > my_largest - tol && score > localBestScore))
215 {
216 my_largest = rad;
217 localBestScore = score;
218 std::memcpy(my_besttile, d, 9 * sizeof(int));
219 }
220 }
221 }
222 }
223 }
224 }
225 }
226 }
227 }
228 }
229 }
230 }
231 if (my_largest > largest + tol || (my_largest > largest - tol && localBestScore > bestScore))
232 {
233 #pragma omp critical
234 {
235 if (my_largest > largest + tol || (my_largest > largest - tol && localBestScore > bestScore))
236 {
237 largest = my_largest;
238 bestScore = localBestScore;
239 std::memcpy(tilemat, my_besttile, 9 * sizeof(int));
240 }
241 }
242 }
243 }
244 radius = largest;
245 }
246
247
main(int argc,char * argv[])248 int main(int argc, char* argv[])
249 {
250 double prim[9];
251 int target = 0;
252 char* pend;
253 int verbose = 0;
254 int maxentry = 4;
255 for (int i = 1; i < argc; i++)
256 {
257 if (i <= argc)
258 {
259 if (!std::strcmp(argv[i], "--ptvs"))
260 {
261 for (int j = 0; j < 9; j++)
262 {
263 prim[j] = strtod(argv[i + j + 1], &pend);
264 }
265 i += 9;
266 }
267 if (!std::strcmp(argv[i], "--target"))
268 {
269 target = strtol(argv[i + 1], &pend, 10);
270 i++;
271 }
272 if (!std::strcmp(argv[i], "--verbose"))
273 {
274 verbose = 1;
275 }
276 if (!std::strcmp(argv[i], "--maxentry"))
277 {
278 maxentry = strtol(argv[i + 1], &pend, 10);
279 i++;
280 }
281 }
282 }
283
284 int besttile[9];
285 double super[9];
286 double radius;
287 getBestTile(prim, target, besttile, radius, maxentry);
288 getSupercell(prim, besttile, super);
289 if (verbose)
290 {
291 std::cout << "Best Tiling Matrix = " << std::endl;
292 for (int i = 0; i < 3; i++)
293 {
294 for (int j = 0; j < 3; j++)
295 {
296 std::cout << besttile[i * 3 + j] << " ";
297 }
298 std::cout << std::endl;
299 }
300 std::cout << "Best Supercell = " << std::endl;
301 for (int i = 0; i < 3; i++)
302 {
303 for (int j = 0; j < 3; j++)
304 {
305 std::cout << super[i * 3 + j] << " ";
306 }
307 std::cout << std::endl;
308 }
309 std::cout << "radius = " << radius << std::endl;
310 double score = getScore(super);
311 std::cout << "score = " << score << std::endl;
312 int trialtile[9];
313 trialtile[0] = -1;
314 trialtile[1] = 3;
315 trialtile[2] = -1;
316 trialtile[3] = -3;
317 trialtile[4] = 1;
318 trialtile[5] = 1;
319 trialtile[6] = 1;
320 trialtile[7] = 1;
321 trialtile[8] = 1;
322 double ts[9];
323 getSupercell(prim, trialtile, ts);
324 std::cout << "Trial Supercell = " << std::endl;
325 for (int i = 0; i < 3; i++)
326 {
327 for (int j = 0; j < 3; j++)
328 {
329 std::cout << ts[i * 3 + j] << " ";
330 }
331 std::cout << std::endl;
332 }
333 std::cout << "radius = " << SimCellRad(ts) << std::endl;
334 score = getScore(ts);
335 std::cout << "score = " << score << std::endl;
336 }
337 else
338 {
339 std::cout << radius << " ";
340 for (int i = 0; i < 9; i++)
341 {
342 std::cout << besttile[i] << " ";
343 }
344 for (int i = 0; i < 9; i++)
345 {
346 std::cout << super[i] << " ";
347 }
348 std::cout << std::endl;
349 }
350 }
351