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: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
8 //                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 //                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
10 //                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
11 //
12 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
13 //////////////////////////////////////////////////////////////////////////////////////
14 
15 #ifndef QMCPLUSPLUS_QUADRATURE_H
16 #define QMCPLUSPLUS_QUADRATURE_H
17 
18 #include <assert.h>
19 #include "Numerics/Ylm.h"
20 #include "type_traits/scalar_traits.h"
21 #include "QMCWaveFunctions/LCAO/SoaSphericalTensor.h"
22 
23 namespace qmcplusplus
24 {
25 template<class T>
26 struct Quadrature3D
27 {
28   typedef T RealType;
29   typedef TinyVector<T, 3> PosType;
30   typedef OHMMS_PRECISION_FULL mRealType;
31 
32   int nk;
33   typedef enum
34   {
35     SINGLE,
36     TETRA,
37     OCTA,
38     ICOSA
39   } SymmType;
40   SymmType symmetry;
41   int lexact;
42   RealType A, B, C, D;
43   std::vector<PosType> xyz_m;
44   std::vector<RealType> weight_m;
45   bool quad_ok;
46   const bool fail_abort;
47 
quad_okQuadrature3D48   Quadrature3D(int rule, bool request_abort = true) : quad_ok(true), fail_abort(request_abort)
49   {
50     A = B = C = D = 0;
51     switch (rule)
52     {
53     case 1:
54       nk       = 1;
55       symmetry = SINGLE;
56       lexact   = 0;
57       A        = 1.0;
58       break;
59     case 2:
60       nk       = 4;
61       symmetry = TETRA;
62       lexact   = 2;
63       A        = 0.25;
64       break;
65     case 3:
66       nk       = 6;
67       symmetry = OCTA;
68       lexact   = 3;
69       A        = 1.0 / 6.0;
70       break;
71     case 4:
72       nk       = 12;
73       symmetry = ICOSA;
74       lexact   = 5;
75       A        = 1.0 / 12.0;
76       B        = 1.0 / 12.0;
77       break;
78     case 5:
79       nk       = 18;
80       symmetry = OCTA;
81       lexact   = 5;
82       A        = 1.0 / 30.0;
83       B        = 1.0 / 15.0;
84       break;
85     case 6:
86       nk       = 26;
87       symmetry = OCTA;
88       lexact   = 7;
89       A        = 1.0 / 21.0;
90       B        = 4.0 / 105.0;
91       C        = 27.0 / 840.0;
92       break;
93     case 7:
94       nk       = 50;
95       symmetry = OCTA;
96       lexact   = 11;
97       A        = 4.0 / 315.0;
98       B        = 64.0 / 2835.0;
99       C        = 27.0 / 1280.0;
100       D        = 14641.0 / 725760.0;
101       break;
102     default:
103       ERRORMSG("Unrecognized spherical quadrature rule " << rule << ".");
104       abort();
105     }
106     // First, build a_i, b_i, and c_i points
107     std::vector<PosType> a, b, c, d;
108     RealType p = 1.0 / std::sqrt(2.0);
109     RealType q = 1.0 / std::sqrt(3.0);
110     RealType r = 1.0 / std::sqrt(11.0);
111     RealType s = 3.0 / std::sqrt(11.0);
112     if (symmetry == SINGLE)
113     {
114       a.push_back(PosType(1.0, 0.0, 0.0));
115     }
116     else if (symmetry == TETRA)
117     {
118       a.push_back(PosType(q, q, q));
119       a.push_back(PosType(q, -q, -q));
120       a.push_back(PosType(-q, q, -q));
121       a.push_back(PosType(-q, -q, q));
122     }
123     else if (symmetry == OCTA)
124     {
125       a.push_back(PosType(1.0, 0.0, 0.0));
126       a.push_back(PosType(-1.0, 0.0, 0.0));
127       a.push_back(PosType(0.0, 1.0, 0.0));
128       a.push_back(PosType(0.0, -1.0, 0.0));
129       a.push_back(PosType(0.0, 0.0, 1.0));
130       a.push_back(PosType(0.0, 0.0, -1.0));
131       b.push_back(PosType(p, p, 0.0));
132       b.push_back(PosType(p, -p, 0.0));
133       b.push_back(PosType(-p, p, 0.0));
134       b.push_back(PosType(-p, -p, 0.0));
135       b.push_back(PosType(p, 0.0, p));
136       b.push_back(PosType(p, 0.0, -p));
137       b.push_back(PosType(-p, 0.0, p));
138       b.push_back(PosType(-p, 0.0, -p));
139       b.push_back(PosType(0.0, p, p));
140       b.push_back(PosType(0.0, p, -p));
141       b.push_back(PosType(0.0, -p, p));
142       b.push_back(PosType(0.0, -p, -p));
143       c.push_back(PosType(q, q, q));
144       c.push_back(PosType(q, q, -q));
145       c.push_back(PosType(q, -q, q));
146       c.push_back(PosType(q, -q, -q));
147       c.push_back(PosType(-q, q, q));
148       c.push_back(PosType(-q, q, -q));
149       c.push_back(PosType(-q, -q, q));
150       c.push_back(PosType(-q, -q, -q));
151       d.push_back(PosType(r, r, s));
152       d.push_back(PosType(r, r, -s));
153       d.push_back(PosType(r, -r, s));
154       d.push_back(PosType(r, -r, -s));
155       d.push_back(PosType(-r, r, s));
156       d.push_back(PosType(-r, r, -s));
157       d.push_back(PosType(-r, -r, s));
158       d.push_back(PosType(-r, -r, -s));
159       d.push_back(PosType(r, s, r));
160       d.push_back(PosType(r, s, -r));
161       d.push_back(PosType(r, -s, r));
162       d.push_back(PosType(r, -s, -r));
163       d.push_back(PosType(-r, s, r));
164       d.push_back(PosType(-r, s, -r));
165       d.push_back(PosType(-r, -s, r));
166       d.push_back(PosType(-r, -s, -r));
167       d.push_back(PosType(s, r, r));
168       d.push_back(PosType(s, r, -r));
169       d.push_back(PosType(s, -r, r));
170       d.push_back(PosType(s, -r, -r));
171       d.push_back(PosType(-s, r, r));
172       d.push_back(PosType(-s, r, -r));
173       d.push_back(PosType(-s, -r, r));
174       d.push_back(PosType(-s, -r, -r));
175     }
176     else if (symmetry == ICOSA)
177     {
178       mRealType t, p; // theta and phi
179       // a points
180       t = 0.0;
181       p = 0.0;
182       a.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
183       t = M_PI;
184       p = 0.0;
185       a.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
186       // b points
187       for (int k = 0; k < 5; k++)
188       {
189         t = std::atan(2.0);
190         p = (mRealType)(2 * k + 0) * M_PI / 5.0;
191         b.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
192         t = M_PI - std::atan(2.0);
193         p = (mRealType)(2 * k + 1) * M_PI / 5.0;
194         b.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
195       }
196       // c points
197       mRealType t1 = std::acos((2.0 + std::sqrt(5.0)) / std::sqrt(15.0 + 6.0 * std::sqrt(5.0)));
198       mRealType t2 = std::acos(1.0 / std::sqrt(15.0 + 6.0 * std::sqrt(5.0)));
199       for (int k = 0; k < 5; k++)
200       {
201         t = t1;
202         p = (mRealType)(2 * k + 1) * M_PI / 5.0;
203         c.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
204         t = t2;
205         p = (mRealType)(2 * k + 1) * M_PI / 5.0;
206         c.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
207         t = M_PI - t1;
208         p = (mRealType)(2 * k + 0) * M_PI / 5.0;
209         c.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
210         t = M_PI - t2;
211         p = (mRealType)(2 * k + 0) * M_PI / 5.0;
212         c.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
213       }
214     }
215     // Now, construct rule
216     if (std::abs(A) > 1.0e-10)
217       for (int i = 0; i < a.size(); i++)
218       {
219         xyz_m.push_back(a[i]);
220         weight_m.push_back(A);
221       }
222     if (std::abs(B) > 1.0e-10)
223       for (int i = 0; i < b.size(); i++)
224       {
225         xyz_m.push_back(b[i]);
226         weight_m.push_back(B);
227       }
228     if (std::abs(C) > 1.0e-10)
229       for (int i = 0; i < c.size(); i++)
230       {
231         xyz_m.push_back(c[i]);
232         weight_m.push_back(C);
233       }
234     if (std::abs(D) > 1.0e-10)
235       for (int i = 0; i < d.size(); i++)
236       {
237         xyz_m.push_back(d[i]);
238         weight_m.push_back(D);
239       }
240     // Finally, check the rule for correctness
241     assert(xyz_m.size() == nk);
242     assert(weight_m.size() == nk);
243     double wSum = 0.0;
244 
245     for (int k = 0; k < nk; k++)
246     {
247       PosType r  = xyz_m[k];
248       double nrm = dot(r, r);
249       assert(std::abs(nrm - 1.0) < 2 * std::numeric_limits<float>::epsilon());
250       wSum += weight_m[k];
251       //cout << pp_nonloc->xyz_m[k] << " " << pp_nonloc->weight_m[k] << std::endl;
252     }
253     assert(std::abs(wSum - 1.0) < 2 * std::numeric_limits<float>::epsilon());
254     // Check the quadrature rule
255     // using complex spherical harmonics
256     CheckQuadratureRule(lexact);
257     // using real spherical harmonics
258     CheckQuadratureRuleReal(lexact);
259   }
260 
CheckQuadratureRuleQuadrature3D261   void CheckQuadratureRule(int lexact)
262   {
263     std::vector<PosType>& grid = xyz_m;
264     std::vector<RealType>& w   = weight_m;
265     for (int l1 = 0; l1 <= lexact; l1++)
266       for (int l2 = 0; l2 <= (lexact - l1); l2++)
267         for (int m1 = -l1; m1 <= l1; m1++)
268           for (int m2 = -l2; m2 <= l2; m2++)
269           {
270             std::complex<mRealType> sum(0.0, 0.0);
271             for (int k = 0; k < grid.size(); k++)
272             {
273               std::complex<mRealType> v1 = Ylm(l1, m1, grid[k]);
274               std::complex<mRealType> v2 = Ylm(l2, m2, grid[k]);
275               sum += 4.0 * M_PI * w[k] * qmcplusplus::conj(v1) * v2;
276             }
277             mRealType re = real(sum);
278             mRealType im = imag(sum);
279             if ((l1 == l2) && (m1 == m2))
280               re -= 1.0;
281             if ((std::abs(im) > 7 * std::numeric_limits<float>::epsilon()) ||
282                 (std::abs(re) > 7 * std::numeric_limits<float>::epsilon()))
283             {
284               app_error() << "Broken spherical quadrature for " << grid.size() << "-point rule.\n" << std::endl;
285               app_error() << "  Should be zero:  Real part = " << re << " Imaginary part = " << im << std::endl;
286               quad_ok = false;
287               if (fail_abort)
288                 APP_ABORT("Give up");
289             }
290             //   	    fprintf (stderr, "(l1,m1,l2m,m2) = (%2d,%2d,%2d,%2d)  sum = (%20.16f %20.16f)\n",
291             //   	     l1, m1, l2, m2, real(sum), imag(sum));
292           }
293   }
294 
CheckQuadratureRuleRealQuadrature3D295   void CheckQuadratureRuleReal(int lexact)
296   {
297     std::vector<PosType>& grid = xyz_m;
298     std::vector<RealType>& w   = weight_m;
299     SoaSphericalTensor<RealType> Ylm(lexact);
300     const RealType* restrict Ylm_v = Ylm[0];
301     for (int l1 = 0; l1 <= lexact; l1++)
302       for (int l2 = 0; l2 <= (lexact - l1); l2++)
303         for (int m1 = -l1; m1 <= l1; m1++)
304           for (int m2 = -l2; m2 <= l2; m2++)
305           {
306             mRealType sum(0.0);
307             for (int k = 0; k < grid.size(); k++)
308             {
309               Ylm.evaluateV(grid[k][0], grid[k][1], grid[k][2]);
310               RealType v1 = Ylm_v[Ylm.index(l1, m1)];
311               RealType v2 = Ylm_v[Ylm.index(l2, m2)];
312               sum += 4.0 * M_PI * w[k] * v1 * v2;
313             }
314             if ((l1 == l2) && (m1 == m2))
315               sum -= 1.0;
316             if (std::abs(sum) > 16 * std::numeric_limits<float>::epsilon())
317             {
318               app_error() << "Broken real spherical quadrature for " << grid.size() << "-point rule.\n" << std::endl;
319               app_error() << "  Should be zero:  " << sum << std::endl;
320               quad_ok = false;
321               if (fail_abort)
322                 APP_ABORT("Give up");
323             }
324           }
325   }
326 };
327 
328 } // namespace qmcplusplus
329 
330 #endif
331