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