JDFTx  1.2.0
SphericalHarmonics.h
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_ELECTRONIC_SPHERICALHARMONICS_H
21 #define JDFTX_ELECTRONIC_SPHERICALHARMONICS_H
22 
23 #include <core/vector3.h>
24 
26 inline double bessel_jl(int l, double x)
27 { double xInv=1./x, xInvSq = xInv * xInv;
28  if(fabs(x) > 1.+0.1*l)
29  { double s, c; sincos(x, &s, &c);
30  switch(l)
31  { case 0: return xInv * s;
32  case 1: return xInv * (xInv*s - c);
33  case 2: return xInv * ((3*xInvSq-1)*s - 3*xInv*c);
34  case 3: return xInv * ((15*xInvSq-6)*xInv*s + (1-15*xInvSq)*c);
35  case 4: return xInv * ((1+xInvSq*(-45+xInvSq*105))*s + xInv*(10-105*xInvSq)*c);
36  case 5: return xInv * (xInv*(15+xInvSq*(-420+xInvSq*945))*s + (-1+xInvSq*(105-945*xInvSq))*c);
37  case 6: return xInv * ((-1+xInvSq*(210+xInvSq*(-4725 + 10395*xInvSq)))*s + xInv*(-21+xInvSq*(1260 - 10395*xInvSq))*c);
38  default: return 0.; //unsupported l
39  }
40  }
41  else //Series expansions about 0 to prevent roundoff errors (accurate to 1 in 1e15 at the crossover for each l)
42  { double term = 1.;
43  for(int i=3; i<=2*l+1; i+=2)
44  term *= x/i;
45  double ret = term;
46  double xSq = x*x;
47  for(int i=2; i<=14; i+=2)
48  { term *= -xSq/(i*(i+2*l+1));
49  ret += term;
50  }
51  return ret;
52  }
53 }
54 
55 
58 namespace YlmInternal
59 {
60  template<int lm> __hostanddev__ double Ylm(double x, double y, double z); //flat-indexed by lm := l*(l+1) + m
61  #define Power pow //For code auto-generated by Mathematica
62  #define DECLARE_Ylm(lm,code) \
63  template<> __hostanddev__ double Ylm<lm>(double x, double y, double z) { return code; }
64  DECLARE_Ylm(0, 0.28209479177387814)
65  DECLARE_Ylm(1, 0.4886025119029199*y)
66  DECLARE_Ylm(2, 0.4886025119029199*z)
67  DECLARE_Ylm(3, 0.4886025119029199*x)
68  DECLARE_Ylm(4, 1.0925484305920792*x*y)
69  DECLARE_Ylm(5, 1.0925484305920792*y*z)
70  DECLARE_Ylm(6, -0.31539156525252005*(Power(x,2) + Power(y,2) - 2.*Power(z,2)))
71  DECLARE_Ylm(7, 1.0925484305920792*x*z)
72  DECLARE_Ylm(8, 0.5462742152960396*(x - y)*(x + y))
73  DECLARE_Ylm(9, -0.5900435899266435*y*(-3.*Power(x,2) + Power(y,2)))
74  DECLARE_Ylm(10, 2.890611442640554*x*y*z)
75  DECLARE_Ylm(11, -0.4570457994644658*y*(Power(x,2) + Power(y,2) - 4.*Power(z,2)))
76  DECLARE_Ylm(12, 0.3731763325901154*z*(-3.*(Power(x,2) + Power(y,2)) + 2.*Power(z,2)))
77  DECLARE_Ylm(13, -0.4570457994644658*x*(Power(x,2) + Power(y,2) - 4.*Power(z,2)))
78  DECLARE_Ylm(14, 1.445305721320277*(x - y)*(x + y)*z)
79  DECLARE_Ylm(15, 0.5900435899266435*x*(Power(x,2) - 3.*Power(y,2)))
80  DECLARE_Ylm(16, 2.5033429417967046*x*(x - y)*y*(x + y))
81  DECLARE_Ylm(17, -1.7701307697799304*y*(-3.*Power(x,2) + Power(y,2))*z)
82  DECLARE_Ylm(18, -0.9461746957575601*x*y*(Power(x,2) + Power(y,2) - 6.*Power(z,2)))
83  DECLARE_Ylm(19, -0.6690465435572892*y*z*(3.*(Power(x,2) + Power(y,2)) - 4.*Power(z,2)))
84  DECLARE_Ylm(20, 0.03526184897173477*(9.*Power(Power(x,2) + Power(y,2),2) - 72.*(Power(x,2) + Power(y,2))*Power(z,2) + 24.*Power(z,4)))
85  DECLARE_Ylm(21, -0.6690465435572892*x*z*(3.*(Power(x,2) + Power(y,2)) - 4.*Power(z,2)))
86  DECLARE_Ylm(22, -0.47308734787878004*(x - y)*(x + y)*(Power(x,2) + Power(y,2) - 6.*Power(z,2)))
87  DECLARE_Ylm(23, 1.7701307697799304*x*(Power(x,2) - 3.*Power(y,2))*z)
88  DECLARE_Ylm(24, 0.6258357354491761*(Power(x,4) - 6.*Power(x,2)*Power(y,2) + Power(y,4)))
89  DECLARE_Ylm(25, 0.6563820568401701*y*(5.*Power(x,4) - 10.*Power(x,2)*Power(y,2) + Power(y,4)))
90  DECLARE_Ylm(26, 8.302649259524166*x*(x - y)*y*(x + y)*z)
91  DECLARE_Ylm(27, 0.4892382994352504*y*(-3.*Power(x,2) + Power(y,2))*(Power(x,2) + Power(y,2) - 8.*Power(z,2)))
92  DECLARE_Ylm(28, -4.793536784973324*x*y*z*(Power(x,2) + Power(y,2) - 2.*Power(z,2)))
93  DECLARE_Ylm(29, 0.45294665119569694*y*(Power(Power(x,2) + Power(y,2),2) - 12.*(Power(x,2) + Power(y,2))*Power(z,2) + 8.*Power(z,4)))
94  DECLARE_Ylm(30, 0.1169503224534236*z*(15.*Power(Power(x,2) + Power(y,2),2) - 40.*(Power(x,2) + Power(y,2))*Power(z,2) + 8.*Power(z,4)))
95  DECLARE_Ylm(31, 0.45294665119569694*x*(Power(Power(x,2) + Power(y,2),2) - 12.*(Power(x,2) + Power(y,2))*Power(z,2) + 8.*Power(z,4)))
96  DECLARE_Ylm(32, -2.396768392486662*(x - y)*(x + y)*z*(Power(x,2) + Power(y,2) - 2.*Power(z,2)))
97  DECLARE_Ylm(33, -0.4892382994352504*x*(Power(x,2) - 3.*Power(y,2))*(Power(x,2) + Power(y,2) - 8.*Power(z,2)))
98  DECLARE_Ylm(34, 2.0756623148810416*(Power(x,4) - 6.*Power(x,2)*Power(y,2) + Power(y,4))*z)
99  DECLARE_Ylm(35, 0.6563820568401701*x*(Power(x,4) - 10.*Power(x,2)*Power(y,2) + 5.*Power(y,4)))
100  DECLARE_Ylm(36, 1.3663682103838286*x*y*(3.*Power(x,4) - 10.*Power(x,2)*Power(y,2) + 3.*Power(y,4)))
101  DECLARE_Ylm(37, 2.366619162231752*y*(5.*Power(x,4) - 10.*Power(x,2)*Power(y,2) + Power(y,4))*z)
102  DECLARE_Ylm(38, -2.0182596029148967*x*(x - y)*y*(x + y)*(Power(x,2) + Power(y,2) - 10.*Power(z,2)))
103  DECLARE_Ylm(39, 0.9212052595149236*y*(-3.*Power(x,2) + Power(y,2))*z*(3.*(Power(x,2) + Power(y,2)) - 8.*Power(z,2)))
104  DECLARE_Ylm(40, 0.9212052595149236*x*y*(Power(Power(x,2) + Power(y,2),2) - 16.*(Power(x,2) + Power(y,2))*Power(z,2) + 16.*Power(z,4)))
105  DECLARE_Ylm(41, 0.5826213625187314*y*z*(5.*Power(Power(x,2) + Power(y,2),2) - 20.*(Power(x,2) + Power(y,2))*Power(z,2) + 8.*Power(z,4)))
106  DECLARE_Ylm(42, 0.06356920226762842*(-5.*Power(Power(x,2) + Power(y,2),3) + 90.*Power(Power(x,2) + Power(y,2),2)*Power(z,2) - 120.*(Power(x,2) + Power(y,2))*Power(z,4) + 16.*Power(z,6)))
107  DECLARE_Ylm(43, 0.5826213625187314*x*z*(5.*Power(Power(x,2) + Power(y,2),2) - 20.*(Power(x,2) + Power(y,2))*Power(z,2) + 8.*Power(z,4)))
108  DECLARE_Ylm(44, 0.4606026297574618*(x - y)*(x + y)*(Power(Power(x,2) + Power(y,2),2) - 16.*(Power(x,2) + Power(y,2))*Power(z,2) + 16.*Power(z,4)))
109  DECLARE_Ylm(45, -0.9212052595149236*x*(Power(x,2) - 3.*Power(y,2))*z*(3.*(Power(x,2) + Power(y,2)) - 8.*Power(z,2)))
110  DECLARE_Ylm(46, -0.5045649007287242*(Power(x,4) - 6.*Power(x,2)*Power(y,2) + Power(y,4))*(Power(x,2) + Power(y,2) - 10.*Power(z,2)))
111  DECLARE_Ylm(47, 2.366619162231752*x*(Power(x,4) - 10.*Power(x,2)*Power(y,2) + 5.*Power(y,4))*z)
112  DECLARE_Ylm(48, 0.6831841051919143*(Power(x,6) - 15.*Power(x,4)*Power(y,2) + 15.*Power(x,2)*Power(y,4) - Power(y,6)))
113  #undef DECLARE_Ylm
114  #undef Power
115 }
116 
118 template<int lm> __hostanddev__ double Ylm(const vector3<>& qhat)
119 { return YlmInternal::Ylm<lm>(qhat[0],qhat[1],qhat[2]);
120 }
121 
123 template<int l, int m> __hostanddev__ double Ylm(const vector3<>& qhat)
124 { return Ylm<l*(l+1)+m>(qhat);
125 }
126 
128 #define SwitchTemplate_lm(l,m,fTemplate,argList) \
129  switch(l*(l+1)+m) \
130  { case 0: fTemplate<0,0> argList; break; \
131  case 1: fTemplate<1,-1> argList; break; \
132  case 2: fTemplate<1,0> argList; break; \
133  case 3: fTemplate<1,1> argList; break; \
134  case 4: fTemplate<2,-2> argList; break; \
135  case 5: fTemplate<2,-1> argList; break; \
136  case 6: fTemplate<2,0> argList; break; \
137  case 7: fTemplate<2,1> argList; break; \
138  case 8: fTemplate<2,2> argList; break; \
139  case 9: fTemplate<3,-3> argList; break; \
140  case 10: fTemplate<3,-2> argList; break; \
141  case 11: fTemplate<3,-1> argList; break; \
142  case 12: fTemplate<3,0> argList; break; \
143  case 13: fTemplate<3,1> argList; break; \
144  case 14: fTemplate<3,2> argList; break; \
145  case 15: fTemplate<3,3> argList; break; \
146  case 16: fTemplate<4,-4> argList; break; \
147  case 17: fTemplate<4,-3> argList; break; \
148  case 18: fTemplate<4,-2> argList; break; \
149  case 19: fTemplate<4,-1> argList; break; \
150  case 20: fTemplate<4,0> argList; break; \
151  case 21: fTemplate<4,1> argList; break; \
152  case 22: fTemplate<4,2> argList; break; \
153  case 23: fTemplate<4,3> argList; break; \
154  case 24: fTemplate<4,4> argList; break; \
155  case 25: fTemplate<5,-5> argList; break; \
156  case 26: fTemplate<5,-4> argList; break; \
157  case 27: fTemplate<5,-3> argList; break; \
158  case 28: fTemplate<5,-2> argList; break; \
159  case 29: fTemplate<5,-1> argList; break; \
160  case 30: fTemplate<5,0> argList; break; \
161  case 31: fTemplate<5,1> argList; break; \
162  case 32: fTemplate<5,2> argList; break; \
163  case 33: fTemplate<5,3> argList; break; \
164  case 34: fTemplate<5,4> argList; break; \
165  case 35: fTemplate<5,5> argList; break; \
166  case 36: fTemplate<6,-6> argList; break; \
167  case 37: fTemplate<6,-5> argList; break; \
168  case 38: fTemplate<6,-4> argList; break; \
169  case 39: fTemplate<6,-3> argList; break; \
170  case 40: fTemplate<6,-2> argList; break; \
171  case 41: fTemplate<6,-1> argList; break; \
172  case 42: fTemplate<6,0> argList; break; \
173  case 43: fTemplate<6,1> argList; break; \
174  case 44: fTemplate<6,2> argList; break; \
175  case 45: fTemplate<6,3> argList; break; \
176  case 46: fTemplate<6,4> argList; break; \
177  case 47: fTemplate<6,5> argList; break; \
178  case 48: fTemplate<6,6> argList; break; \
179  }
180 
182 template<int l, int m> void set_Ylm(const vector3<> qHat, double& result) { result = Ylm<l,m>(qHat); }
183 inline double Ylm(int l, int m, const vector3<>& qHat) { double result=0.; SwitchTemplate_lm(l,m, set_Ylm, (qHat, result)); return result; }
184 
186 struct YlmProdTerm
187 { int l, m;
188  double coeff;
189  YlmProdTerm(int l, int m, double coeff) : l(l), m(m), coeff(coeff) {}
190 };
191 
194 inline std::vector<YlmProdTerm> expandYlmProd(int lm1, int lm2)
195 { if(lm2 > lm1) std::swap(lm1, lm2);
196  std::vector<YlmProdTerm> result;
197  #define ADD(l,m,coeff) result.push_back(YlmProdTerm(l,m,coeff)) //shorthand for use in Mathematica generated list below:
198  switch(lm2 + (lm1*(lm1+1))/2)
199  { case 0: ADD(0,0,0.28209479177387814); break;
200  case 1: ADD(1,-1,0.28209479177387814); break;
201  case 2: ADD(0,0,0.28209479177387814); ADD(2,0,-0.126156626101008); ADD(2,2,-0.2185096861184158); break;
202  case 3: ADD(1,0,0.28209479177387814); break;
203  case 4: ADD(2,-1,0.2185096861184158); break;
204  case 5: ADD(0,0,0.28209479177387814); ADD(2,0,0.252313252202016); break;
205  case 6: ADD(1,1,0.28209479177387814); break;
206  case 7: ADD(2,-2,0.2185096861184158); break;
207  case 8: ADD(2,1,0.2185096861184158); break;
208  case 9: ADD(0,0,0.28209479177387814); ADD(2,0,-0.126156626101008); ADD(2,2,0.2185096861184158); break;
209  case 10: ADD(2,-2,0.28209479177387814); break;
210  case 11: ADD(1,1,0.2185096861184158); ADD(3,1,-0.058399170081901854); ADD(3,3,-0.2261790131595403); break;
211  case 12: ADD(3,-2,0.18467439092237178); break;
212  case 13: ADD(1,-1,0.2185096861184158); ADD(3,-3,0.2261790131595403); ADD(3,-1,-0.058399170081901854); break;
213  case 14: ADD(0,0,0.28209479177387814); ADD(2,0,-0.18022375157286857); ADD(4,0,0.04029925596769687); ADD(4,4,-0.23841361350444806); break;
214  case 15: ADD(2,-1,0.28209479177387814); break;
215  case 16: ADD(1,0,0.2185096861184158); ADD(3,0,-0.14304816810266882); ADD(3,2,-0.18467439092237178); break;
216  case 17: ADD(1,-1,0.2185096861184158); ADD(3,-1,0.23359668032760741); break;
217  case 18: ADD(3,-2,0.18467439092237178); break;
218  case 19: ADD(2,1,0.15607834722743988); ADD(4,1,-0.06371871843402754); ADD(4,3,-0.16858388283618386); break;
219  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;
220  case 21: ADD(2,0,0.28209479177387814); break;
221  case 22: ADD(1,-1,-0.126156626101008); ADD(3,-1,0.20230065940342062); break;
222  case 23: ADD(1,0,0.252313252202016); ADD(3,0,0.2477666950834761); break;
223  case 24: ADD(1,1,-0.126156626101008); ADD(3,1,0.20230065940342062); break;
224  case 25: ADD(2,-2,-0.18022375157286857); ADD(4,-2,0.15607834722743988); break;
225  case 26: ADD(2,-1,0.09011187578643429); ADD(4,-1,0.2207281154418226); break;
226  case 27: ADD(0,0,0.28209479177387814); ADD(2,0,0.18022375157286857); ADD(4,0,0.24179553580618124); break;
227  case 28: ADD(2,1,0.28209479177387814); break;
228  case 29: ADD(3,-2,0.18467439092237178); break;
229  case 30: ADD(1,1,0.2185096861184158); ADD(3,1,0.23359668032760741); break;
230  case 31: ADD(1,0,0.2185096861184158); ADD(3,0,-0.14304816810266882); ADD(3,2,0.18467439092237178); break;
231  case 32: ADD(2,-1,0.15607834722743988); ADD(4,-3,0.16858388283618386); ADD(4,-1,-0.06371871843402754); break;
232  case 33: ADD(2,-2,0.15607834722743988); ADD(4,-2,0.18022375157286857); break;
233  case 34: ADD(2,1,0.09011187578643429); ADD(4,1,0.2207281154418226); break;
234  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;
235  case 36: ADD(2,2,0.28209479177387814); break;
236  case 37: ADD(1,-1,-0.2185096861184158); ADD(3,-3,0.2261790131595403); ADD(3,-1,0.058399170081901854); break;
237  case 38: ADD(3,2,0.18467439092237178); break;
238  case 39: ADD(1,1,0.2185096861184158); ADD(3,1,-0.058399170081901854); ADD(3,3,0.2261790131595403); break;
239  case 40: ADD(4,-4,0.23841361350444806); break;
240  case 41: ADD(2,-1,-0.15607834722743988); ADD(4,-3,0.16858388283618386); ADD(4,-1,0.06371871843402754); break;
241  case 42: ADD(2,2,-0.18022375157286857); ADD(4,2,0.15607834722743988); break;
242  case 43: ADD(2,1,0.15607834722743988); ADD(4,1,-0.06371871843402754); ADD(4,3,0.16858388283618386); break;
243  case 44: ADD(0,0,0.28209479177387814); ADD(2,0,-0.18022375157286857); ADD(4,0,0.04029925596769687); ADD(4,4,0.23841361350444806); break;
244  case 45: ADD(3,-3,0.28209479177387814); break;
245  case 46: ADD(2,2,0.2261790131595403); ADD(4,2,-0.04352817137756816); ADD(4,4,-0.23032943298089034); break;
246  case 47: ADD(4,-3,0.16286750396763996); break;
247  case 48: ADD(2,-2,0.2261790131595403); ADD(4,-4,0.23032943298089034); ADD(4,-2,-0.04352817137756816); break;
248  case 49: ADD(1,1,0.2261790131595403); ADD(3,1,-0.09403159725795937); ADD(5,1,0.01694331772935932); ADD(5,5,-0.2455320005465369); break;
249  case 50: ADD(3,2,0.1486770096793976); ADD(5,2,-0.04482780509623635); ADD(5,4,-0.1552880720369528); break;
250  case 51: ADD(3,-3,-0.21026104350168); ADD(5,-3,0.12679217987703037); break;
251  case 52: ADD(3,-2,0.1486770096793976); ADD(5,-4,0.1552880720369528); ADD(5,-2,-0.04482780509623635); break;
252  case 53: ADD(1,-1,0.2261790131595403); ADD(3,-1,-0.09403159725795937); ADD(5,-5,0.2455320005465369); ADD(5,-1,0.01694331772935932); break;
253  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;
254  case 55: ADD(3,-2,0.28209479177387814); break;
255  case 56: ADD(2,1,0.18467439092237178); ADD(4,1,-0.07539300438651343); ADD(4,3,-0.19947114020071635); break;
256  case 57: ADD(2,-2,0.18467439092237178); ADD(4,-2,0.21324361862292307); break;
257  case 58: ADD(2,-1,0.18467439092237178); ADD(4,-3,0.19947114020071635); ADD(4,-1,-0.07539300438651343); break;
258  case 59: ADD(1,0,0.18467439092237178); ADD(3,0,-0.18806319451591874); ADD(5,0,0.05357947514468781); ADD(5,4,-0.19018826981554557); break;
259  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;
260  case 61: ADD(5,-2,0.19018826981554557); break;
261  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;
262  case 63: ADD(5,-4,0.19018826981554557); break;
263  case 64: ADD(2,1,0.1486770096793976); ADD(4,1,-0.09932258459927992); ADD(6,1,0.022177545476549994); ADD(6,5,-0.1801712311720527); break;
264  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;
265  case 66: ADD(3,-1,0.28209479177387814); break;
266  case 67: ADD(2,0,0.20230065940342062); ADD(2,2,0.058399170081901854); ADD(4,0,-0.15078600877302686); ADD(4,2,-0.16858388283618386); break;
267  case 68: ADD(2,-1,0.23359668032760741); ADD(4,-1,0.23841361350444806); break;
268  case 69: ADD(2,-2,-0.058399170081901854); ADD(4,-2,0.16858388283618386); break;
269  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;
270  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;
271  case 72: ADD(1,-1,0.20230065940342062); ADD(3,-1,0.126156626101008); ADD(5,-1,0.22731846124334898); break;
272  case 73: ADD(3,-2,0.11516471649044517); ADD(5,-2,0.17361734258475534); break;
273  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;
274  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;
275  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;
276  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;
277  case 78: ADD(3,0,0.28209479177387814); break;
278  case 79: ADD(2,-1,-0.14304816810266882); ADD(4,-1,0.19466390027300617); break;
279  case 80: ADD(2,0,0.2477666950834761); ADD(4,0,0.24623252122982908); break;
280  case 81: ADD(2,1,-0.14304816810266882); ADD(4,1,0.19466390027300617); break;
281  case 82: ADD(3,-2,-0.18806319451591874); ADD(5,-2,0.14175796661021045); break;
282  case 83: ADD(1,-1,-0.14304816810266882); ADD(3,-1,0.05947080387175903); ADD(5,-1,0.21431790057875125); break;
283  case 84: ADD(1,0,0.2477666950834761); ADD(3,0,0.168208834801344); ADD(5,0,0.23961469724456466); break;
284  case 85: ADD(1,1,-0.14304816810266882); ADD(3,1,0.05947080387175903); ADD(5,1,0.21431790057875125); break;
285  case 86: ADD(3,2,-0.18806319451591874); ADD(5,2,0.14175796661021045); break;
286  case 87: ADD(4,-3,-0.20355072686733566); ADD(6,-3,0.10864734032983336); break;
287  case 88: ADD(2,-2,-0.18806319451591874); ADD(4,-2,-0.04441841017299272); ADD(6,-2,0.17742036381239995); break;
288  case 89: ADD(2,-1,0.05947080387175903); ADD(4,-1,0.09932258459927992); ADD(6,-1,0.22177545476549995); break;
289  case 90: ADD(0,0,0.28209479177387814); ADD(2,0,0.168208834801344); ADD(4,0,0.15386988642211533); ADD(6,0,0.23708793386528085); break;
290  case 91: ADD(3,1,0.28209479177387814); break;
291  case 92: ADD(2,-2,-0.058399170081901854); ADD(4,-2,0.16858388283618386); break;
292  case 93: ADD(2,1,0.23359668032760741); ADD(4,1,0.23841361350444806); break;
293  case 94: ADD(2,0,0.20230065940342062); ADD(2,2,-0.058399170081901854); ADD(4,0,-0.15078600877302686); ADD(4,2,0.16858388283618386); break;
294  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;
295  case 96: ADD(3,-2,0.11516471649044517); ADD(5,-2,0.17361734258475534); break;
296  case 97: ADD(1,1,0.20230065940342062); ADD(3,1,0.126156626101008); ADD(5,1,0.22731846124334898); break;
297  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;
298  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;
299  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;
300  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;
301  case 102: ADD(2,-2,0.1456731240789439); ADD(4,-2,0.11468784191000729); ADD(6,-2,0.17178652858087154); break;
302  case 103: ADD(2,1,0.05947080387175903); ADD(4,1,0.09932258459927992); ADD(6,1,0.22177545476549995); break;
303  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;
304  case 105: ADD(3,2,0.28209479177387814); break;
305  case 106: ADD(2,-1,-0.18467439092237178); ADD(4,-3,0.19947114020071635); ADD(4,-1,0.07539300438651343); break;
306  case 107: ADD(2,2,0.18467439092237178); ADD(4,2,0.21324361862292307); break;
307  case 108: ADD(2,1,0.18467439092237178); ADD(4,1,-0.07539300438651343); ADD(4,3,0.19947114020071635); break;
308  case 109: ADD(5,-4,0.19018826981554557); break;
309  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;
310  case 111: ADD(5,2,0.19018826981554557); break;
311  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;
312  case 113: ADD(1,0,0.18467439092237178); ADD(3,0,-0.18806319451591874); ADD(5,0,0.05357947514468781); ADD(5,4,0.19018826981554557); break;
313  case 114: ADD(2,-1,0.1486770096793976); ADD(4,-1,-0.09932258459927992); ADD(6,-5,0.1801712311720527); ADD(6,-1,0.022177545476549994); break;
314  case 115: ADD(4,-4,0.15171775404828514); ADD(6,-4,0.18818271355849853); break;
315  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;
316  case 117: ADD(2,2,-0.18806319451591874); ADD(4,2,-0.04441841017299272); ADD(6,2,0.17742036381239995); break;
317  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;
318  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;
319  case 120: ADD(3,3,0.28209479177387814); break;
320  case 121: ADD(2,-2,-0.2261790131595403); ADD(4,-4,0.23032943298089034); ADD(4,-2,0.04352817137756816); break;
321  case 122: ADD(4,3,0.16286750396763996); break;
322  case 123: ADD(2,2,0.2261790131595403); ADD(4,2,-0.04352817137756816); ADD(4,4,0.23032943298089034); break;
323  case 124: ADD(1,-1,-0.2261790131595403); ADD(3,-1,0.09403159725795937); ADD(5,-5,0.2455320005465369); ADD(5,-1,-0.01694331772935932); break;
324  case 125: ADD(3,-2,-0.1486770096793976); ADD(5,-4,0.1552880720369528); ADD(5,-2,0.04482780509623635); break;
325  case 126: ADD(3,3,-0.21026104350168); ADD(5,3,0.12679217987703037); break;
326  case 127: ADD(3,2,0.1486770096793976); ADD(5,2,-0.04482780509623635); ADD(5,4,0.1552880720369528); break;
327  case 128: ADD(1,1,0.2261790131595403); ADD(3,1,-0.09403159725795937); ADD(5,1,0.01694331772935932); ADD(5,5,0.2455320005465369); break;
328  case 129: ADD(6,-6,0.25480059867297505); break;
329  case 130: ADD(2,-1,-0.1486770096793976); ADD(4,-1,0.09932258459927992); ADD(6,-5,0.1801712311720527); ADD(6,-1,-0.022177545476549994); break;
330  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;
331  case 132: ADD(4,3,-0.20355072686733566); ADD(6,3,0.10864734032983336); break;
332  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;
333  case 134: ADD(2,1,0.1486770096793976); ADD(4,1,-0.09932258459927992); ADD(6,1,0.022177545476549994); ADD(6,5,0.1801712311720527); break;
334  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;
335  }
336  #undef ADD
337  return result;
338 }
339 inline std::vector<YlmProdTerm> expandYlmProd(int l1, int m1, int l2, int m2)
340 { int lm1 = l1*(l1+1) + m1;
341  int lm2 = l2*(l2+1) + m2;
342  return expandYlmProd(lm1, lm2);
343 }
344 
345 #endif // JDFTX_ELECTRONIC_SPHERICALHARMONICS_H
346 
YlmProdTerm::m
int m
angular quantum numbers of current term
Definition: SphericalHarmonics.h:187
YlmInternal
Definition: SphericalHarmonics.h:58
YlmProdTerm
Term in real spherical harmonic expansion of a product of two real spherical harmonics.
Definition: SphericalHarmonics.h:186
YlmProdTerm::coeff
double coeff
coefficient of Ylm with current l and m
Definition: SphericalHarmonics.h:188