1 /*-------------------------------------------------------------------
2 Copyright 2012 Ravishankar Sundararaman
3
4 This file is part of JDFTx.
5
6 JDFTx is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 JDFTx is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with JDFTx. If not, see <http://www.gnu.org/licenses/>.
18 -------------------------------------------------------------------*/
19
20 #ifndef JDFTX_CORE_SPHERICALHARMONICS_H
21 #define JDFTX_CORE_SPHERICALHARMONICS_H
22
23 #include <core/vector3.h>
24
25 //! @addtogroup Utilities
26 //! @{
27 //! @file SphericalHarmonics.h Real spherical Harmonics and spherical bessel functions
28
29 //! Spherical bessel function
bessel_jl(int l,double x)30 inline double bessel_jl(int l, double x)
31 { if(fabs(x) > 1.+0.1*l)
32 { double s, c; sincos(x, &s, &c);
33 double xInv=1./x, xInvSq = xInv * xInv;
34 switch(l)
35 { case 0: return xInv * s;
36 case 1: return xInv * (xInv*s - c);
37 case 2: return xInv * ((3*xInvSq-1)*s - 3*xInv*c);
38 case 3: return xInv * ((15*xInvSq-6)*xInv*s + (1-15*xInvSq)*c);
39 case 4: return xInv * ((1+xInvSq*(-45+xInvSq*105))*s + xInv*(10-105*xInvSq)*c);
40 case 5: return xInv * (xInv*(15+xInvSq*(-420+xInvSq*945))*s + (-1+xInvSq*(105-945*xInvSq))*c);
41 case 6: return xInv * ((-1+xInvSq*(210+xInvSq*(-4725 + 10395*xInvSq)))*s + xInv*(-21+xInvSq*(1260 - 10395*xInvSq))*c);
42 default: return 0.; //unsupported l
43 }
44 }
45 else //Series expansions about 0 to prevent roundoff errors (accurate to 1 in 1e15 at the crossover for each l)
46 { double term = 1.;
47 for(int i=3; i<=2*l+1; i+=2)
48 term *= x/i;
49 double ret = term;
50 double xSq = x*x;
51 for(int i=2; i<=14; i+=2)
52 { term *= -xSq/(i*(i+2*l+1));
53 ret += term;
54 }
55 return ret;
56 }
57 }
58
59
60 //! Auto-generated and hand-tweaked code for computing real spherical harmonics and their Clebsch-Gordon coefficients
61 //! The external interface to these functions follows the end of this namespace block
62 namespace YlmInternal
63 {
64 template<int lm> __hostanddev__ double Ylm(double x, double y, double z); //flat-indexed by lm := l*(l+1) + m
65 #define X2 double x2 = x*x;
66 #define Y2 double y2 = y*y;
67 #define Z2 double z2 = z*z;
68 #define X2Y2 double x2y2 = x*x + y*y;
69 #define X2Y2i X2 Y2 double x2y2 = x2 + y2;
70 #define DECLARE_Ylm(lm,code) \
71 template<> __hostanddev__ double Ylm<lm>(double x, double y, double z) { code; }
72 DECLARE_Ylm(-2, return 0.) //Corner case in YlmPrime; needed for compilation but never used
73 DECLARE_Ylm(-1, return 0.) //Corner case in YlmPrime; needed for compilation but never used
74 DECLARE_Ylm(0, return 0.28209479177387814)
75 DECLARE_Ylm(1, return 0.4886025119029199*y)
76 DECLARE_Ylm(2, return 0.4886025119029199*z)
77 DECLARE_Ylm(3, return 0.4886025119029199*x)
78 DECLARE_Ylm(4, return 1.0925484305920792*x*y)
79 DECLARE_Ylm(5, return 1.0925484305920792*y*z)
80 DECLARE_Ylm(6, return -0.31539156525252005*(x*x + y*y - 2.*z*z))
81 DECLARE_Ylm(7, return 1.0925484305920792*x*z)
82 DECLARE_Ylm(8, return 0.5462742152960396*(x-y)*(x+y))
83 DECLARE_Ylm(9, return -0.5900435899266435*y*(y*y-3.*x*x))
84 DECLARE_Ylm(10, return 2.890611442640554*x*y*z)
85 DECLARE_Ylm(11, return -0.4570457994644658*y*(x*x + y*y - 4.*z*z))
86 DECLARE_Ylm(12, return 0.3731763325901154*z*(2.*z*z -3.*(x*x + y*y)))
87 DECLARE_Ylm(13, return -0.4570457994644658*x*(x*x + y*y - 4.*z*z))
88 DECLARE_Ylm(14, return 1.445305721320277*(x-y)*(x+y)*z)
89 DECLARE_Ylm(15, return 0.5900435899266435*x*(x*x - 3.*y*y))
90 DECLARE_Ylm(16, return 2.5033429417967046*x*y*(x-y)*(x+y))
91 DECLARE_Ylm(17, return -1.7701307697799304*y*z*(y*y - 3.*x*x))
92 DECLARE_Ylm(18, return -0.9461746957575601*x*y*(x*x + y*y - 6.*z*z))
93 DECLARE_Ylm(19, return -0.6690465435572892*y*z*(3.*(x*x + y*y) - 4.*z*z))
94 DECLARE_Ylm(20, X2Y2 Z2 return 0.03526184897173477*(9.*x2y2*(x2y2 - 8*z2) + 24.*z2*z2))
95 DECLARE_Ylm(21, return -0.6690465435572892*x*z*(3.*(x*x + y*y) - 4.*z*z))
96 DECLARE_Ylm(22, X2 Y2 return -0.47308734787878004*(x2 - y2)*(x2 + y2 - 6.*z*z))
97 DECLARE_Ylm(23, return 1.7701307697799304*x*z*(x*x - 3.*y*y))
98 DECLARE_Ylm(24, X2 Y2 return 0.6258357354491761*(x2*(x2 - 6.*y2) + y2*y2))
99 DECLARE_Ylm(25, X2 Y2 return 0.6563820568401701*y*(5.*x2*(x2 - 2.*y2) + y2*y2))
100 DECLARE_Ylm(26, return 8.302649259524166*x*y*z*(x-y)*(x+y))
101 DECLARE_Ylm(27, X2 Y2 return 0.4892382994352504*y*(y2 -3.*x2)*(x2 + y2 - 8.*z*z))
102 DECLARE_Ylm(28, return -4.793536784973324*x*y*z*(x*x + y*y - 2.*z*z))
103 DECLARE_Ylm(29, X2Y2 Z2 return 0.45294665119569694*y*(x2y2*(x2y2 - 12.*z2) + 8.*z2*z2))
104 DECLARE_Ylm(30, X2Y2 Z2 return 0.1169503224534236*z*(15.*x2y2*x2y2 - 8.*z2*(5.*x2y2 - z2)))
105 DECLARE_Ylm(31, X2Y2 Z2 return 0.45294665119569694*x*(x2y2*(x2y2 - 12.*z2) + 8.*z2*z2))
106 DECLARE_Ylm(32, X2 Y2 return -2.396768392486662*(x2-y2)*z*(x2 + y2 - 2.*z*z))
107 DECLARE_Ylm(33, X2 Y2 return -0.4892382994352504*x*(x2 - 3.*y2)*(x2 + y2 - 8.*z*z))
108 DECLARE_Ylm(34, X2 Y2 return 2.0756623148810416*z*(x2*(x2 - 6.*y2) + y2*y2))
109 DECLARE_Ylm(35, X2 Y2 return 0.6563820568401701*x*(x2*(x2 - 10.*y2) + 5.*y2*y2))
110 DECLARE_Ylm(36, X2 Y2 return 1.3663682103838286*x*y*(x2*(3.*x2 - 10.*y2) + 3.*y2*y2))
111 DECLARE_Ylm(37, X2 Y2 return 2.366619162231752*y*z*(5.*x2*(x2 - 2.*y2) + y2*y2))
112 DECLARE_Ylm(38, X2 Y2 return -2.0182596029148967*x*y*(x2 - y2)*(x2 + y2 - 10.*z*z))
113 DECLARE_Ylm(39, X2 Y2 return 0.9212052595149236*y*z*(y2 - 3.*x2)*(3.*(x2 + y2) - 8.*z*z))
114 DECLARE_Ylm(40, X2Y2 Z2 return 0.9212052595149236*x*y*(x2y2*(x2y2 - 16.*z2) + 16.*z2*z2))
115 DECLARE_Ylm(41, X2Y2 Z2 return 0.5826213625187314*y*z*(5.*x2y2*(x2y2 - 4.*z2) + 8.*z2*z2))
116 DECLARE_Ylm(42, X2Y2 Z2 return 0.06356920226762842*(5.*x2y2*x2y2*(18.*z2 - x2y2) + 8.*z2*z2*(2.*z2 - 15.*x2y2)))
117 DECLARE_Ylm(43, X2Y2 Z2 return 0.5826213625187314*x*z*(5.*x2y2*(x2y2 - 4.*z2) + 8.*z2*z2))
118 DECLARE_Ylm(44, X2Y2i Z2 return 0.4606026297574618*(x2-y2)*(x2y2*(x2y2 - 16.*z2) + 16.*z2*z2))
119 DECLARE_Ylm(45, X2 Y2 return -0.9212052595149236*x*z*(x2 - 3.*y2)*(3.*(x2 + y2) - 8.*z*z))
120 DECLARE_Ylm(46, X2 Y2 return -0.5045649007287242*(x2*(x2 - 6.*y2) + y2*y2)*(x2 + y2 - 10.*z*z))
121 DECLARE_Ylm(47, X2 Y2 return 2.366619162231752*x*z*(x2*(x2 - 10.*y2) + 5.*y2*y2))
122 DECLARE_Ylm(48, X2 Y2 return 0.6831841051919143*(x2*x2*(x2 - 15.*y2) + y2*y2*(15.*x2 - y2)))
123 #undef DECLARE_Ylm
124 #undef X2
125 #undef Y2
126 #undef Z2
127 #undef X2Y2
128 #undef X2Y2i
129 }
130
131 //! Index by combined lm := l*(l+1)+m index (useful when static-looping over all l,m)
Ylm(const vector3<> & qhat)132 template<int lm> __hostanddev__ double Ylm(const vector3<>& qhat)
133 { return YlmInternal::Ylm<lm>(qhat[0],qhat[1],qhat[2]);
134 }
135
136 //! Index by l and m separately
Ylm(const vector3<> & qhat)137 template<int l, int m> __hostanddev__ double Ylm(const vector3<>& qhat)
138 { return Ylm<l*(l+1)+m>(qhat);
139 }
140
141 //! Switch a function templated over l,m for all supported l,m with parenthesis enclosed argument list argList
142 #define SwitchTemplate_lm(l,m,fTemplate,argList) \
143 switch(l*(l+1)+m) \
144 { case 0: fTemplate<0,0> argList; break; \
145 case 1: fTemplate<1,-1> argList; break; \
146 case 2: fTemplate<1,0> argList; break; \
147 case 3: fTemplate<1,1> argList; break; \
148 case 4: fTemplate<2,-2> argList; break; \
149 case 5: fTemplate<2,-1> argList; break; \
150 case 6: fTemplate<2,0> argList; break; \
151 case 7: fTemplate<2,1> argList; break; \
152 case 8: fTemplate<2,2> argList; break; \
153 case 9: fTemplate<3,-3> argList; break; \
154 case 10: fTemplate<3,-2> argList; break; \
155 case 11: fTemplate<3,-1> argList; break; \
156 case 12: fTemplate<3,0> argList; break; \
157 case 13: fTemplate<3,1> argList; break; \
158 case 14: fTemplate<3,2> argList; break; \
159 case 15: fTemplate<3,3> argList; break; \
160 case 16: fTemplate<4,-4> argList; break; \
161 case 17: fTemplate<4,-3> argList; break; \
162 case 18: fTemplate<4,-2> argList; break; \
163 case 19: fTemplate<4,-1> argList; break; \
164 case 20: fTemplate<4,0> argList; break; \
165 case 21: fTemplate<4,1> argList; break; \
166 case 22: fTemplate<4,2> argList; break; \
167 case 23: fTemplate<4,3> argList; break; \
168 case 24: fTemplate<4,4> argList; break; \
169 case 25: fTemplate<5,-5> argList; break; \
170 case 26: fTemplate<5,-4> argList; break; \
171 case 27: fTemplate<5,-3> argList; break; \
172 case 28: fTemplate<5,-2> argList; break; \
173 case 29: fTemplate<5,-1> argList; break; \
174 case 30: fTemplate<5,0> argList; break; \
175 case 31: fTemplate<5,1> argList; break; \
176 case 32: fTemplate<5,2> argList; break; \
177 case 33: fTemplate<5,3> argList; break; \
178 case 34: fTemplate<5,4> argList; break; \
179 case 35: fTemplate<5,5> argList; break; \
180 case 36: fTemplate<6,-6> argList; break; \
181 case 37: fTemplate<6,-5> argList; break; \
182 case 38: fTemplate<6,-4> argList; break; \
183 case 39: fTemplate<6,-3> argList; break; \
184 case 40: fTemplate<6,-2> argList; break; \
185 case 41: fTemplate<6,-1> argList; break; \
186 case 42: fTemplate<6,0> argList; break; \
187 case 43: fTemplate<6,1> argList; break; \
188 case 44: fTemplate<6,2> argList; break; \
189 case 45: fTemplate<6,3> argList; break; \
190 case 46: fTemplate<6,4> argList; break; \
191 case 47: fTemplate<6,5> argList; break; \
192 case 48: fTemplate<6,6> argList; break; \
193 }
194
195 //! Use above macro to provide a non-templated version of the function
set_Ylm(const vector3<> qHat,double & result)196 template<int l, int m> void set_Ylm(const vector3<> qHat, double& result) { result = Ylm<l,m>(qHat); }
Ylm(int l,int m,const vector3<> & qHat)197 inline double Ylm(int l, int m, const vector3<>& qHat) { double result=0.; SwitchTemplate_lm(l,m, set_Ylm, (qHat, result)); return result; }
198
199 //! Term in real spherical harmonic expansion of a product of two real spherical harmonics
200 struct YlmProdTerm
201 { int l, m; //!< angular quantum numbers of current term
202 double coeff; //!< coefficient of Ylm with current l and m
YlmProdTermYlmProdTerm203 YlmProdTerm(int l, int m, double coeff) : l(l), m(m), coeff(coeff) {}
204 };
205
206 //! Real spherical harmonic expansion of product of two real spherical harmonics
207 //! (effectively returns list of non-zero Clebsch-Gordon coefficients in the modified real Ylm basis)
expandYlmProd(int lm1,int lm2)208 inline std::vector<YlmProdTerm> expandYlmProd(int lm1, int lm2)
209 { if(lm2 > lm1) std::swap(lm1, lm2);
210 std::vector<YlmProdTerm> result;
211 #define ADD(l,m,coeff) result.push_back(YlmProdTerm(l,m,coeff)) //shorthand for use in Mathematica generated list below:
212 switch(lm2 + (lm1*(lm1+1))/2)
213 { case 0: ADD(0,0,0.28209479177387814); break;
214 case 1: ADD(1,-1,0.28209479177387814); break;
215 case 2: ADD(0,0,0.28209479177387814); ADD(2,0,-0.126156626101008); ADD(2,2,-0.2185096861184158); break;
216 case 3: ADD(1,0,0.28209479177387814); break;
217 case 4: ADD(2,-1,0.2185096861184158); break;
218 case 5: ADD(0,0,0.28209479177387814); ADD(2,0,0.252313252202016); break;
219 case 6: ADD(1,1,0.28209479177387814); break;
220 case 7: ADD(2,-2,0.2185096861184158); break;
221 case 8: ADD(2,1,0.2185096861184158); break;
222 case 9: ADD(0,0,0.28209479177387814); ADD(2,0,-0.126156626101008); ADD(2,2,0.2185096861184158); break;
223 case 10: ADD(2,-2,0.28209479177387814); break;
224 case 11: ADD(1,1,0.2185096861184158); ADD(3,1,-0.058399170081901854); ADD(3,3,-0.2261790131595403); break;
225 case 12: ADD(3,-2,0.18467439092237178); break;
226 case 13: ADD(1,-1,0.2185096861184158); ADD(3,-3,0.2261790131595403); ADD(3,-1,-0.058399170081901854); break;
227 case 14: ADD(0,0,0.28209479177387814); ADD(2,0,-0.18022375157286857); ADD(4,0,0.04029925596769687); ADD(4,4,-0.23841361350444806); break;
228 case 15: ADD(2,-1,0.28209479177387814); break;
229 case 16: ADD(1,0,0.2185096861184158); ADD(3,0,-0.14304816810266882); ADD(3,2,-0.18467439092237178); break;
230 case 17: ADD(1,-1,0.2185096861184158); ADD(3,-1,0.23359668032760741); break;
231 case 18: ADD(3,-2,0.18467439092237178); break;
232 case 19: ADD(2,1,0.15607834722743988); ADD(4,1,-0.06371871843402754); ADD(4,3,-0.16858388283618386); break;
233 case 20: ADD(0,0,0.28209479177387814); ADD(2,0,0.09011187578643429); ADD(2,2,-0.15607834722743988); ADD(4,0,-0.1611970238707875); ADD(4,2,-0.18022375157286857); break;
234 case 21: ADD(2,0,0.28209479177387814); break;
235 case 22: ADD(1,-1,-0.126156626101008); ADD(3,-1,0.20230065940342062); break;
236 case 23: ADD(1,0,0.252313252202016); ADD(3,0,0.2477666950834761); break;
237 case 24: ADD(1,1,-0.126156626101008); ADD(3,1,0.20230065940342062); break;
238 case 25: ADD(2,-2,-0.18022375157286857); ADD(4,-2,0.15607834722743988); break;
239 case 26: ADD(2,-1,0.09011187578643429); ADD(4,-1,0.2207281154418226); break;
240 case 27: ADD(0,0,0.28209479177387814); ADD(2,0,0.18022375157286857); ADD(4,0,0.24179553580618124); break;
241 case 28: ADD(2,1,0.28209479177387814); break;
242 case 29: ADD(3,-2,0.18467439092237178); break;
243 case 30: ADD(1,1,0.2185096861184158); ADD(3,1,0.23359668032760741); break;
244 case 31: ADD(1,0,0.2185096861184158); ADD(3,0,-0.14304816810266882); ADD(3,2,0.18467439092237178); break;
245 case 32: ADD(2,-1,0.15607834722743988); ADD(4,-3,0.16858388283618386); ADD(4,-1,-0.06371871843402754); break;
246 case 33: ADD(2,-2,0.15607834722743988); ADD(4,-2,0.18022375157286857); break;
247 case 34: ADD(2,1,0.09011187578643429); ADD(4,1,0.2207281154418226); break;
248 case 35: ADD(0,0,0.28209479177387814); ADD(2,0,0.09011187578643429); ADD(2,2,0.15607834722743988); ADD(4,0,-0.1611970238707875); ADD(4,2,0.18022375157286857); break;
249 case 36: ADD(2,2,0.28209479177387814); break;
250 case 37: ADD(1,-1,-0.2185096861184158); ADD(3,-3,0.2261790131595403); ADD(3,-1,0.058399170081901854); break;
251 case 38: ADD(3,2,0.18467439092237178); break;
252 case 39: ADD(1,1,0.2185096861184158); ADD(3,1,-0.058399170081901854); ADD(3,3,0.2261790131595403); break;
253 case 40: ADD(4,-4,0.23841361350444806); break;
254 case 41: ADD(2,-1,-0.15607834722743988); ADD(4,-3,0.16858388283618386); ADD(4,-1,0.06371871843402754); break;
255 case 42: ADD(2,2,-0.18022375157286857); ADD(4,2,0.15607834722743988); break;
256 case 43: ADD(2,1,0.15607834722743988); ADD(4,1,-0.06371871843402754); ADD(4,3,0.16858388283618386); break;
257 case 44: ADD(0,0,0.28209479177387814); ADD(2,0,-0.18022375157286857); ADD(4,0,0.04029925596769687); ADD(4,4,0.23841361350444806); break;
258 case 45: ADD(3,-3,0.28209479177387814); break;
259 case 46: ADD(2,2,0.2261790131595403); ADD(4,2,-0.04352817137756816); ADD(4,4,-0.23032943298089034); break;
260 case 47: ADD(4,-3,0.16286750396763996); break;
261 case 48: ADD(2,-2,0.2261790131595403); ADD(4,-4,0.23032943298089034); ADD(4,-2,-0.04352817137756816); break;
262 case 49: ADD(1,1,0.2261790131595403); ADD(3,1,-0.09403159725795937); ADD(5,1,0.01694331772935932); ADD(5,5,-0.2455320005465369); break;
263 case 50: ADD(3,2,0.1486770096793976); ADD(5,2,-0.04482780509623635); ADD(5,4,-0.1552880720369528); break;
264 case 51: ADD(3,-3,-0.21026104350168); ADD(5,-3,0.12679217987703037); break;
265 case 52: ADD(3,-2,0.1486770096793976); ADD(5,-4,0.1552880720369528); ADD(5,-2,-0.04482780509623635); break;
266 case 53: ADD(1,-1,0.2261790131595403); ADD(3,-1,-0.09403159725795937); ADD(5,-5,0.2455320005465369); ADD(5,-1,0.01694331772935932); break;
267 case 54: ADD(0,0,0.28209479177387814); ADD(2,0,-0.21026104350168); ADD(4,0,0.07693494321105766); ADD(6,0,-0.011854396693264043); ADD(6,6,-0.25480059867297505); break;
268 case 55: ADD(3,-2,0.28209479177387814); break;
269 case 56: ADD(2,1,0.18467439092237178); ADD(4,1,-0.07539300438651343); ADD(4,3,-0.19947114020071635); break;
270 case 57: ADD(2,-2,0.18467439092237178); ADD(4,-2,0.21324361862292307); break;
271 case 58: ADD(2,-1,0.18467439092237178); ADD(4,-3,0.19947114020071635); ADD(4,-1,-0.07539300438651343); break;
272 case 59: ADD(1,0,0.18467439092237178); ADD(3,0,-0.18806319451591874); ADD(5,0,0.05357947514468781); ADD(5,4,-0.19018826981554557); break;
273 case 60: ADD(1,1,0.18467439092237178); ADD(3,1,0.11516471649044517); ADD(3,3,-0.1486770096793976); ADD(5,1,-0.08300496597356405); ADD(5,3,-0.1793112203849454); break;
274 case 61: ADD(5,-2,0.19018826981554557); break;
275 case 62: ADD(1,-1,0.18467439092237178); ADD(3,-3,0.1486770096793976); ADD(3,-1,0.11516471649044517); ADD(5,-3,0.1793112203849454); ADD(5,-1,-0.08300496597356405); break;
276 case 63: ADD(5,-4,0.19018826981554557); break;
277 case 64: ADD(2,1,0.1486770096793976); ADD(4,1,-0.09932258459927992); ADD(6,1,0.022177545476549994); ADD(6,5,-0.1801712311720527); break;
278 case 65: ADD(0,0,0.28209479177387814); ADD(4,0,-0.1795148674924679); ADD(4,4,-0.15171775404828514); ADD(6,0,0.07112638015958425); ADD(6,4,-0.18818271355849853); break;
279 case 66: ADD(3,-1,0.28209479177387814); break;
280 case 67: ADD(2,0,0.20230065940342062); ADD(2,2,0.058399170081901854); ADD(4,0,-0.15078600877302686); ADD(4,2,-0.16858388283618386); break;
281 case 68: ADD(2,-1,0.23359668032760741); ADD(4,-1,0.23841361350444806); break;
282 case 69: ADD(2,-2,-0.058399170081901854); ADD(4,-2,0.16858388283618386); break;
283 case 70: ADD(1,1,-0.058399170081901854); ADD(3,1,0.1456731240789439); ADD(3,3,0.09403159725795937); ADD(5,1,-0.0656211873953095); ADD(5,3,-0.14175796661021045); break;
284 case 71: ADD(1,0,0.23359668032760741); ADD(3,0,0.05947080387175903); ADD(3,2,-0.11516471649044517); ADD(5,0,-0.1694331772935932); ADD(5,2,-0.17361734258475534); break;
285 case 72: ADD(1,-1,0.20230065940342062); ADD(3,-1,0.126156626101008); ADD(5,-1,0.22731846124334898); break;
286 case 73: ADD(3,-2,0.11516471649044517); ADD(5,-2,0.17361734258475534); break;
287 case 74: ADD(1,-1,0.058399170081901854); ADD(3,-3,-0.09403159725795937); ADD(3,-1,-0.1456731240789439); ADD(5,-3,0.14175796661021045); ADD(5,-1,0.0656211873953095); break;
288 case 75: ADD(2,2,-0.09403159725795937); ADD(4,2,0.13325523051897814); ADD(4,4,0.11752006695060024); ADD(6,2,-0.04435509095309999); ADD(6,4,-0.1214714192760309); break;
289 case 76: ADD(2,1,0.11516471649044517); ADD(4,1,0.10257992428141023); ADD(4,3,-0.06785024228911189); ADD(6,1,-0.08589326429043577); ADD(6,3,-0.16297101049475005); break;
290 case 77: ADD(0,0,0.28209479177387814); ADD(2,0,0.126156626101008); ADD(2,2,-0.1456731240789439); ADD(4,0,0.025644981070352558); ADD(4,2,-0.11468784191000729); ADD(6,0,-0.17781595039896067); ADD(6,2,-0.17178652858087154); break;
291 case 78: ADD(3,0,0.28209479177387814); break;
292 case 79: ADD(2,-1,-0.14304816810266882); ADD(4,-1,0.19466390027300617); break;
293 case 80: ADD(2,0,0.2477666950834761); ADD(4,0,0.24623252122982908); break;
294 case 81: ADD(2,1,-0.14304816810266882); ADD(4,1,0.19466390027300617); break;
295 case 82: ADD(3,-2,-0.18806319451591874); ADD(5,-2,0.14175796661021045); break;
296 case 83: ADD(1,-1,-0.14304816810266882); ADD(3,-1,0.05947080387175903); ADD(5,-1,0.21431790057875125); break;
297 case 84: ADD(1,0,0.2477666950834761); ADD(3,0,0.168208834801344); ADD(5,0,0.23961469724456466); break;
298 case 85: ADD(1,1,-0.14304816810266882); ADD(3,1,0.05947080387175903); ADD(5,1,0.21431790057875125); break;
299 case 86: ADD(3,2,-0.18806319451591874); ADD(5,2,0.14175796661021045); break;
300 case 87: ADD(4,-3,-0.20355072686733566); ADD(6,-3,0.10864734032983336); break;
301 case 88: ADD(2,-2,-0.18806319451591874); ADD(4,-2,-0.04441841017299272); ADD(6,-2,0.17742036381239995); break;
302 case 89: ADD(2,-1,0.05947080387175903); ADD(4,-1,0.09932258459927992); ADD(6,-1,0.22177545476549995); break;
303 case 90: ADD(0,0,0.28209479177387814); ADD(2,0,0.168208834801344); ADD(4,0,0.15386988642211533); ADD(6,0,0.23708793386528085); break;
304 case 91: ADD(3,1,0.28209479177387814); break;
305 case 92: ADD(2,-2,-0.058399170081901854); ADD(4,-2,0.16858388283618386); break;
306 case 93: ADD(2,1,0.23359668032760741); ADD(4,1,0.23841361350444806); break;
307 case 94: ADD(2,0,0.20230065940342062); ADD(2,2,-0.058399170081901854); ADD(4,0,-0.15078600877302686); ADD(4,2,0.16858388283618386); break;
308 case 95: ADD(1,-1,-0.058399170081901854); ADD(3,-3,-0.09403159725795937); ADD(3,-1,0.1456731240789439); ADD(5,-3,0.14175796661021045); ADD(5,-1,-0.0656211873953095); break;
309 case 96: ADD(3,-2,0.11516471649044517); ADD(5,-2,0.17361734258475534); break;
310 case 97: ADD(1,1,0.20230065940342062); ADD(3,1,0.126156626101008); ADD(5,1,0.22731846124334898); break;
311 case 98: ADD(1,0,0.23359668032760741); ADD(3,0,0.05947080387175903); ADD(3,2,0.11516471649044517); ADD(5,0,-0.1694331772935932); ADD(5,2,0.17361734258475534); break;
312 case 99: ADD(1,1,-0.058399170081901854); ADD(3,1,0.1456731240789439); ADD(3,3,-0.09403159725795937); ADD(5,1,-0.0656211873953095); ADD(5,3,0.14175796661021045); break;
313 case 100: ADD(2,-2,-0.09403159725795937); ADD(4,-4,-0.11752006695060024); ADD(4,-2,0.13325523051897814); ADD(6,-4,0.1214714192760309); ADD(6,-2,-0.04435509095309999); break;
314 case 101: ADD(2,-1,0.11516471649044517); ADD(4,-3,0.06785024228911189); ADD(4,-1,0.10257992428141023); ADD(6,-3,0.16297101049475005); ADD(6,-1,-0.08589326429043577); break;
315 case 102: ADD(2,-2,0.1456731240789439); ADD(4,-2,0.11468784191000729); ADD(6,-2,0.17178652858087154); break;
316 case 103: ADD(2,1,0.05947080387175903); ADD(4,1,0.09932258459927992); ADD(6,1,0.22177545476549995); break;
317 case 104: ADD(0,0,0.28209479177387814); ADD(2,0,0.126156626101008); ADD(2,2,0.1456731240789439); ADD(4,0,0.025644981070352558); ADD(4,2,0.11468784191000729); ADD(6,0,-0.17781595039896067); ADD(6,2,0.17178652858087154); break;
318 case 105: ADD(3,2,0.28209479177387814); break;
319 case 106: ADD(2,-1,-0.18467439092237178); ADD(4,-3,0.19947114020071635); ADD(4,-1,0.07539300438651343); break;
320 case 107: ADD(2,2,0.18467439092237178); ADD(4,2,0.21324361862292307); break;
321 case 108: ADD(2,1,0.18467439092237178); ADD(4,1,-0.07539300438651343); ADD(4,3,0.19947114020071635); break;
322 case 109: ADD(5,-4,0.19018826981554557); break;
323 case 110: ADD(1,-1,-0.18467439092237178); ADD(3,-3,0.1486770096793976); ADD(3,-1,-0.11516471649044517); ADD(5,-3,0.1793112203849454); ADD(5,-1,0.08300496597356405); break;
324 case 111: ADD(5,2,0.19018826981554557); break;
325 case 112: ADD(1,1,0.18467439092237178); ADD(3,1,0.11516471649044517); ADD(3,3,0.1486770096793976); ADD(5,1,-0.08300496597356405); ADD(5,3,0.1793112203849454); break;
326 case 113: ADD(1,0,0.18467439092237178); ADD(3,0,-0.18806319451591874); ADD(5,0,0.05357947514468781); ADD(5,4,0.19018826981554557); break;
327 case 114: ADD(2,-1,0.1486770096793976); ADD(4,-1,-0.09932258459927992); ADD(6,-5,0.1801712311720527); ADD(6,-1,0.022177545476549994); break;
328 case 115: ADD(4,-4,0.15171775404828514); ADD(6,-4,0.18818271355849853); break;
329 case 116: ADD(2,-1,-0.11516471649044517); ADD(4,-3,0.06785024228911189); ADD(4,-1,-0.10257992428141023); ADD(6,-3,0.16297101049475005); ADD(6,-1,0.08589326429043577); break;
330 case 117: ADD(2,2,-0.18806319451591874); ADD(4,2,-0.04441841017299272); ADD(6,2,0.17742036381239995); break;
331 case 118: ADD(2,1,0.11516471649044517); ADD(4,1,0.10257992428141023); ADD(4,3,0.06785024228911189); ADD(6,1,-0.08589326429043577); ADD(6,3,0.16297101049475005); break;
332 case 119: ADD(0,0,0.28209479177387814); ADD(4,0,-0.1795148674924679); ADD(4,4,0.15171775404828514); ADD(6,0,0.07112638015958425); ADD(6,4,0.18818271355849853); break;
333 case 120: ADD(3,3,0.28209479177387814); break;
334 case 121: ADD(2,-2,-0.2261790131595403); ADD(4,-4,0.23032943298089034); ADD(4,-2,0.04352817137756816); break;
335 case 122: ADD(4,3,0.16286750396763996); break;
336 case 123: ADD(2,2,0.2261790131595403); ADD(4,2,-0.04352817137756816); ADD(4,4,0.23032943298089034); break;
337 case 124: ADD(1,-1,-0.2261790131595403); ADD(3,-1,0.09403159725795937); ADD(5,-5,0.2455320005465369); ADD(5,-1,-0.01694331772935932); break;
338 case 125: ADD(3,-2,-0.1486770096793976); ADD(5,-4,0.1552880720369528); ADD(5,-2,0.04482780509623635); break;
339 case 126: ADD(3,3,-0.21026104350168); ADD(5,3,0.12679217987703037); break;
340 case 127: ADD(3,2,0.1486770096793976); ADD(5,2,-0.04482780509623635); ADD(5,4,0.1552880720369528); break;
341 case 128: ADD(1,1,0.2261790131595403); ADD(3,1,-0.09403159725795937); ADD(5,1,0.01694331772935932); ADD(5,5,0.2455320005465369); break;
342 case 129: ADD(6,-6,0.25480059867297505); break;
343 case 130: ADD(2,-1,-0.1486770096793976); ADD(4,-1,0.09932258459927992); ADD(6,-5,0.1801712311720527); ADD(6,-1,-0.022177545476549994); break;
344 case 131: ADD(2,-2,0.09403159725795937); ADD(4,-4,-0.11752006695060024); ADD(4,-2,-0.13325523051897814); ADD(6,-4,0.1214714192760309); ADD(6,-2,0.04435509095309999); break;
345 case 132: ADD(4,3,-0.20355072686733566); ADD(6,3,0.10864734032983336); break;
346 case 133: ADD(2,2,-0.09403159725795937); ADD(4,2,0.13325523051897814); ADD(4,4,-0.11752006695060024); ADD(6,2,-0.04435509095309999); ADD(6,4,0.1214714192760309); break;
347 case 134: ADD(2,1,0.1486770096793976); ADD(4,1,-0.09932258459927992); ADD(6,1,0.022177545476549994); ADD(6,5,0.1801712311720527); break;
348 case 135: ADD(0,0,0.28209479177387814); ADD(2,0,-0.21026104350168); ADD(4,0,0.07693494321105766); ADD(6,0,-0.011854396693264043); ADD(6,6,0.25480059867297505); break;
349 }
350 #undef ADD
351 return result;
352 }
353 //! Wrapper function expandYlmProd with individual indices
expandYlmProd(int l1,int m1,int l2,int m2)354 inline std::vector<YlmProdTerm> expandYlmProd(int l1, int m1, int l2, int m2)
355 { int lm1 = l1*(l1+1) + m1;
356 int lm2 = l2*(l2+1) + m2;
357 return expandYlmProd(lm1, lm2);
358 }
359
360 //! Derivative of spherical Harmonic with repect to iDir'th component of qHat with separate l,m indices
YlmPrime(const vector3<> & qHat)361 template<int l, int m> __hostanddev__ vector3<> YlmPrime(const vector3<>& qHat)
362 { vector3<> result;
363 //z-component:
364 if(l > 0)
365 result[2] = sqrt((l*l-m*m)*(2*l+1.)/(2*l-1.)) * Ylm<l-1,m>(qHat);
366 //m-dependent expressions for x and y components:
367 if(m == 0)
368 { if(l>1)
369 { double alpha = sqrt(0.5*l*(l-1)*(2*l+1.)/(2*l-1.));
370 result[0] = -alpha * Ylm<l-1,+1>(qHat);
371 result[1] = -alpha * Ylm<l-1,-1>(qHat);
372 }
373 }
374 else
375 { double alphaM = 0.5*sqrt((l-m)*(l-m-1)*(2*l+1.)/(2*l-1.));
376 double alphaP = 0.5*sqrt((l+m)*(l+m-1)*(2*l+1.)/(2*l-1.));
377 if(m > 0)
378 { if(l > m+1)
379 { result[0] -= Ylm<l-1, m+1 >(qHat) * alphaM;
380 result[1] -= Ylm<l-1,-(m+1)>(qHat) * alphaM;
381 }
382 if(l >= m)
383 { result[0] += Ylm<l-1, m-1 >(qHat) * (alphaP * (m==1 ? sqrt(2.) : 1.));
384 result[1] -= Ylm<l-1,-(m-1)>(qHat) * (alphaP * (m==1 ? 0. : 1.));
385 }
386 }
387 else //m < 0
388 { if(-l <= m)
389 { result[0] += Ylm<l-1, m+1 >(qHat) * (alphaM * (m==-1 ? 0. : 1.));
390 result[1] += Ylm<l-1,-(m+1)>(qHat) * (alphaM * (m==-1 ? sqrt(2.) : 1.));
391 }
392 if(-l < m-1)
393 { result[0] -= Ylm<l-1, m-1 >(qHat) * alphaP;
394 result[1] += Ylm<l-1,-(m-1)>(qHat) * alphaP;
395 }
396 }
397 }
398 return result;
399 }
400
401 //! Derivative of spherical Harmonic with repect to iDir'th component of qHat with combined lm index
402 template<int lm> __hostanddev__ vector3<> YlmPrime(const vector3<>& qHat);
403 #define DECLARE_YlmPrime(l,m) template<> __hostanddev__ vector3<> YlmPrime<l*(l+1)+m>(const vector3<>& qHat) { return YlmPrime<l,m>(qHat); }
404 DECLARE_YlmPrime(0,0)
405 DECLARE_YlmPrime(1,-1)
406 DECLARE_YlmPrime(1,0)
407 DECLARE_YlmPrime(1,1)
408 DECLARE_YlmPrime(2,-2)
409 DECLARE_YlmPrime(2,-1)
410 DECLARE_YlmPrime(2,0)
411 DECLARE_YlmPrime(2,1)
412 DECLARE_YlmPrime(2,2)
413 DECLARE_YlmPrime(3,-3)
414 DECLARE_YlmPrime(3,-2)
415 DECLARE_YlmPrime(3,-1)
416 DECLARE_YlmPrime(3,0)
417 DECLARE_YlmPrime(3,1)
418 DECLARE_YlmPrime(3,2)
419 DECLARE_YlmPrime(3,3)
420 DECLARE_YlmPrime(4,-4)
421 DECLARE_YlmPrime(4,-3)
422 DECLARE_YlmPrime(4,-2)
423 DECLARE_YlmPrime(4,-1)
424 DECLARE_YlmPrime(4,0)
425 DECLARE_YlmPrime(4,1)
426 DECLARE_YlmPrime(4,2)
427 DECLARE_YlmPrime(4,3)
428 DECLARE_YlmPrime(4,4)
429 DECLARE_YlmPrime(5,-5)
430 DECLARE_YlmPrime(5,-4)
431 DECLARE_YlmPrime(5,-3)
432 DECLARE_YlmPrime(5,-2)
433 DECLARE_YlmPrime(5,-1)
434 DECLARE_YlmPrime(5,0)
435 DECLARE_YlmPrime(5,1)
436 DECLARE_YlmPrime(5,2)
437 DECLARE_YlmPrime(5,3)
438 DECLARE_YlmPrime(5,4)
439 DECLARE_YlmPrime(5,5)
440 DECLARE_YlmPrime(6,-6)
441 DECLARE_YlmPrime(6,-5)
442 DECLARE_YlmPrime(6,-4)
443 DECLARE_YlmPrime(6,-3)
444 DECLARE_YlmPrime(6,-2)
445 DECLARE_YlmPrime(6,-1)
446 DECLARE_YlmPrime(6,0)
447 DECLARE_YlmPrime(6,1)
448 DECLARE_YlmPrime(6,2)
449 DECLARE_YlmPrime(6,3)
450 DECLARE_YlmPrime(6,4)
451 DECLARE_YlmPrime(6,5)
452 DECLARE_YlmPrime(6,6)
453 #undef DECLARE_YlmPrime
454
455 //! Non-templated version of YlmPrime (for debugging)
set_YlmPrime(const vector3<> qHat,vector3<> & result)456 template<int l, int m> void set_YlmPrime(const vector3<> qHat, vector3<>& result) { result = YlmPrime<l,m>(qHat); }
YlmPrime(int l,int m,const vector3<> & qHat)457 inline vector3<> YlmPrime(int l, int m, const vector3<>& qHat)
458 { vector3<> result; SwitchTemplate_lm(l,m, set_YlmPrime, (qHat, result)); return result;
459 }
460
461 //! @}
462 #endif // JDFTX_CORE_SPHERICALHARMONICS_H
463
464