1 /*----------------------------------------------------------------------
2 eri_sf.c
3
4 Special functions used in LIBERI.
5
6 This is based on Spherical_Bessel.c and openmx_common.c in OpenMX
7 3.4 package.
8
9 22 Nov. 2001, opnemx_common.c by T. Ozaki
10 08 Nov. 2005 Spherical_Bessel.c by T.Ozaki
11 16 Feb. 2009, M. Toyoda
12 ----------------------------------------------------------------------*/
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <float.h>
16 #include <math.h>
17 #include <assert.h>
18 #include "eri_sf.h"
19 #include "eri_def.h"
20
21
22 #define LMAX_MAX 512
23 #define MAX_ITR 50
24 #define ERR_MAX 1e-12
25
26 #define CG_JMAX 30
27
28
29 static void sph_bessel_asym(double x, int l, double *sb, double *dsb);
30 static void sph_bessel_drec(double x, int l, double *sb, double *dsb);
31 static void sph_bessel_arec(double x, int l, double *sb, double *dsb);
32
33 static double CG(int j1, int m1, int j2, int m2, int j, int m);
34
35
36
37 /*----------------------------------------------------------------------
38 rotation matrix for Complex -> Real conversion of SH
39
40 Z(l, m) = a(m)*Y(l, m) + b(m)*Y(l, -m)
41
42 where
43
44 a(m) = 1/sqrt(2) (m>0), 1 (m=0), and i(-1)^m/sqrt(2) (m<0)
45 b(m) = (-1)^m/sqrt(2) (m>0), 0 (m=0), and -i/sqrt(2) (m<0)
46
47 a(m) = (-1)^m/sqrt(2) (m>0), 1 (m=0), and i/sqrt(2) (m<0)
48 b(m) = 1/sqrt(2) (m>0), 0 (m=0), and -i(-1)^m/sqrt(2) (m<0)
49 ----------------------------------------------------------------------*/
ERI_RCSH_Coeff(int m,double a[2],double b[2])50 static void ERI_RCSH_Coeff(int m, double a[2], double b[2])
51 {
52 double oost;
53
54 oost = 1.0/sqrt(2.0);
55
56 if (m>0) {
57 if (m%2==0) {
58 a[0] = oost; a[1] = 0.0;
59 b[0] = oost; b[1] = 0.0;
60 } else {
61 a[0] = -oost; a[1] = 0.0;
62 b[0] = oost; b[1] = 0.0;
63 }
64 } else if (m<0) {
65 if (abs(m)%2==0) {
66 a[0] = 0.0; a[1] = oost;
67 b[0] = 0.0; b[1] = -oost;
68 } else {
69 a[0] = 0.0; a[1] = oost;
70 b[0] = 0.0; b[1] = oost;
71 }
72 } else {
73 /* m == 0 */
74 a[0] = 1.0; a[1] = 0.0;
75 b[0] = 0.0; b[1] = 0.0;
76 }
77
78 #if 0
79 if (m>0) {
80 if (m%2==0) {
81 a[0] = oost; a[1] = 0.0;
82 b[0] = oost; b[1] = 0.0;
83 } else {
84 a[0] = oost; a[1] = 0.0;
85 b[0] = -oost; b[1] = 0.0;
86 }
87 } else if (m<0) {
88 if (m%2==0) {
89 a[0] = 0.0; a[1] = oost;
90 b[0] = 0.0; b[1] = -oost;
91 } else {
92 a[0] = 0.0; a[1] = -oost;
93 b[0] = 0.0; b[1] = -oost;
94 }
95 } else {
96 /* m == 0 */
97 a[0] = 1.0; a[1] = 0.0;
98 b[0] = 0.0; b[1] = 0.0;
99 }
100 #endif
101 }
102
103
104 /*----------------------------------------------------------------------
105 rotation matrix for Real -> Complex conversion of SH
106
107 Y(l, m) = c(m)*Z(l, m) + d(m)*Z(l, -m)
108
109 where
110
111 c(m) = 1/sqrt(2) (m>0), 1 (m=0), and -i(-1)^m/sqrt(2) (m<0)
112 f(m) = i/sqrt(2) (m>0), 0 (m=0), and (-1)^m/sqrt(2) (m<0)
113 ----------------------------------------------------------------------*/
ERI_RCSH_Coeff_Inverse(int m,double c[2],double d[2])114 static void ERI_RCSH_Coeff_Inverse(int m, double c[2], double d[2])
115 {
116 double na, nb, a[2], b[2];
117
118 if (0==m) {
119 c[0] = 1.0; c[1] = 0.0;
120 d[0] = 0.0; d[1] = 0.0;
121 } else {
122 ERI_RCSH_Coeff(m, a, b);
123 na = a[0]*a[0] + a[1]*a[1];
124 nb = b[0]*b[0] + b[1]*b[1];
125 c[0] = 0.5*a[0]/na;
126 c[1] = -0.5*a[1]/na;
127 d[0] = 0.5*b[0]/nb;
128 d[1] = -0.5*b[1]/nb;
129 }
130 }
131
132
133
134
135 /*----------------------------------------------------------------------
136 ERI_Spherical_Bessel
137
138 Routine to calculate the spherical Bessel function of the first kind
139 j_l(x) of the order from 0 to l and their derivatives.
140
141 Descending recurrence algoritm is used for small x in order to avoid
142 accumulation of rounding errors, whereas ascending recurrence is used
143 for large x because rounding error can be negligible and the
144 descending recurrence becomes too slow for large x.
145
146 Reference:
147 A. jablonski, "Numerical Evaluation of Spherical Bessel Functions of
148 the First Kind", Journal of Computational Physics Vol. 111, pp. 256
149 (1994).
150 ----------------------------------------------------------------------*/
ERI_Spherical_Bessel(double x,int l,double * sb,double * dsb)151 void ERI_Spherical_Bessel(
152 double x, /* (IN) argument x >= 0*/
153 int l, /* (IN) positive integer */
154 double *sb, /* (OUT) function values of order from 0 to l */
155 double *dsb /* (OUT) derivatives of order from 0 to l */
156 )
157 {
158 const double xmin = 1e-10;
159 const double xthresh = 100.0;
160
161 /* negative x value is invalid */
162 assert( x >= 0.0 );
163 assert( l >= 0 );
164
165 if (x<xmin) {
166 /* x is near zero */
167 sph_bessel_asym(x, l, sb, dsb);
168 } else if (x<xthresh) {
169 /* descending recurrence */
170 sph_bessel_drec(x, l, sb, dsb);
171 } else {
172 /* ascending recurrence */
173 sph_bessel_arec(x, l, sb, dsb);
174 }
175 }
176
177
178
179
sph_bessel_asym(double x,int l,double * sb,double * dsb)180 static void sph_bessel_asym(
181 double x,
182 int l,
183 double *sb,
184 double *dsb
185 )
186 {
187 int i;
188 double f = 1.0;
189
190 sb[0] = 1.0;
191 for (i=1; i<=l+1; i++) {
192 f += 2.0;
193 sb[i] = sb[i-1]*x/f;
194 }
195
196 dsb[0] = 0.0;
197 for (i=1; i<=l; i++) {
198 dsb[i] = ((double)i*sb[i-1]-(double)(i+1)*sb[i+1])
199 /(double)(2*i+1);
200 }
201 }
202
203
204
205
sph_bessel_drec(double x,int l,double * sb,double * dsb)206 static void sph_bessel_drec(
207 double x,
208 int l,
209 double *sb,
210 double *dsb
211 )
212 {
213 int i, j, lmax;
214 double j0, j1, sf, tmp, si, co, ix, ix2, tsb[1024], huge;
215
216 huge = sqrt(DBL_MAX);
217
218 /* find an appropriate nmax */
219 lmax = l + 3*(int)x + 20;
220 assert( lmax+1 < LMAX_MAX );
221 if (lmax<LMAX_MAX) lmax = LMAX_MAX;
222
223 /* initial values*/
224 tsb[lmax] = 0.0;
225 tsb[lmax-1] = 1.0e-14;
226
227 /* downward recurrence from nmax-2 to lmax+2 */
228 for (i=lmax-1; (l+2)<i; i--) {
229 tsb[i-1] = (2.0*(double)i+1.0)/x*tsb[i]-tsb[i+1];
230 if (huge<tsb[i-1]){
231 tsb[i] /= tsb[i-1];
232 tsb[i-1] = 1.0;
233 }
234 }
235
236 /* downward recurrence from lmax+1 to 0 */
237 tsb[l+3] /= tsb[l+2];
238 tsb[l+2] = 1.0;
239 for (i=l+2; 0<i; i--) {
240 tsb[i-1] = (2.0*(double)i+1.0)/x*tsb[i] - tsb[i+1];
241 if (huge<tsb[i-1]){
242 tmp = tsb[i-1];
243 for (j=i-1; j<=l+1; j++) { tsb[j] /= tmp; }
244 }
245 }
246
247 /* normalization */
248 si = sin(x);
249 co = cos(x);
250 ix = 1.0/x;
251 ix2 = ix*ix;
252 j0 = si*ix;
253 j1 = si*ix*ix - co*ix;
254
255 if (fabs(tsb[1])<fabs(tsb[0])) {
256 sf = j0/tsb[0];
257 } else {
258 sf = j1/tsb[1];
259 }
260
261 /* tsb to sb */
262 for (i=0; i<=l+1; i++) {
263 sb[i] = tsb[i]*sf;
264 }
265
266 /* derivative of sb */
267 dsb[0] = co*ix - si*ix*ix;
268 for (i=1; i<=l; i++) {
269 dsb[i] = ((double)i*sb[i-1] - (double)(i+1.0)*sb[i+1])
270 /(2.0*(double)i + 1.0);
271 }
272 }
273
274
275
276
sph_bessel_arec(double x,int l,double * sb,double * dsb)277 static void sph_bessel_arec(
278 double x,
279 int l,
280 double *sb,
281 double *dsb
282 )
283 {
284 int i;
285 double si = sin(x), co = cos(x);
286
287 sb[0] = si/x;
288 dsb[0] = co/x-si/x/x;
289 if (l==0) { return; }
290
291 sb[1] = -dsb[0];
292 dsb[1] = sb[0] - 2.0*sb[1]/x;
293 if (l==1) { return; }
294
295 for (i=2; i<=l; i++) {
296 sb[i] = (double)(2*i-1)*sb[i-1]/x-sb[i-2];
297 dsb[i] = sb[i-1]-(double)(i+1)*sb[i]/x;
298 }
299 } /* end if */
300
301
302
303
304 /*----------------------------------------------------------------------
305 ERI_Spherical_Harmonics
306 ----------------------------------------------------------------------*/
ERI_Spherical_Harmonics(int l,int m,double theta,double phi,double Y[2],double dYdt[2],double dYdp[2])307 void ERI_Spherical_Harmonics(
308 int l, /* (IN) azimuthal quantum number */
309 int m, /* (IN) magnetuc quantum number */
310 double theta, /* (IN) argument */
311 double phi, /* (IN) argument */
312 double Y[2], /* (OUT) function value (complex) */
313 double dYdt[2], /* (OUT) derivative w.r.t. theta (complex) */
314 double dYdp[2] /* (OUT) derivative w.r.t. phi (complex) */
315 )
316 {
317 unsigned int i, am;
318 long double f0, f1;
319 double norm, co, si, P, dP;
320
321 am = abs(m);
322
323 /* (l-|m|)! */
324 f0 = 1.0;
325 for (i=2; i<=(l-am); i++) { f0 *= (long double)i; }
326
327 /* (l+|m|)! */
328 f1 = 1.0;
329 for (i=2; i<=(l+am); i++) { f1 *= (long double)i; }
330
331 /* sqrt((2*l+1)/(4*PI)*(l-|m|)!/(l+|m|)!) */
332 norm = sqrt((2.0*(double)l+1.0)/(4.0*PI)*f0/f1);
333
334 /* i^(m-|m|) */
335 if (m<0 && am%2==1) { norm = -norm; }
336
337 /* P_l^|m| */
338 ERI_Associated_Legendre_Polynomial(l, am, cos(theta), &P, &dP);
339
340 /* exp(i*m*phi) */
341 co = cos((double)m*phi);
342 si = sin((double)m*phi);
343
344 Y[0] = norm*P*co;
345 Y[1] = norm*P*si;
346 dYdt[0] = norm*dP*co;
347 dYdt[1] = norm*dP*si;
348 dYdp[0] = -(double)m*Y[1];
349 dYdp[1] = (double)m*Y[0];
350 }
351
352
353
354
355 /*----------------------------------------------------------------------
356 ERI_Real_Spherical_Harmonics
357
358 Real-valued spherical harmonics which is defined as:
359
360 Z(l, m; theta, phi) = a(m)*Y(l, m; theta, phi)
361 + b(m)*Y(l, -m; theta, phi)
362
363 where a(m) and b(m) are defined by real_coeff.
364 ----------------------------------------------------------------------*/
ERI_Real_Spherical_Harmonics(int l,int m,double theta,double phi,double * Z,double * dZdt,double * dZdp)365 void ERI_Real_Spherical_Harmonics(
366 int l, /* (IN) azimuthal quantum number */
367 int m, /* (IN) magnetuc quantum number */
368 double theta, /* (IN) argument */
369 double phi, /* (IN) argument */
370 double *Z, /* (OUT) function value (optional) */
371 double *dZdt, /* (OUT) derivative w.r.t. theta (optional)*/
372 double *dZdp /* (OUT) derivative w.r.t. phi (optional) */
373 )
374 {
375 double Y1[2], dY1dt[2], dY1dp[2];
376 double Y2[2], dY2dt[2], dY2dp[2];
377 double a[2], b[2];
378
379 ERI_Spherical_Harmonics(l, m, theta, phi, Y1, dY1dt, dY1dp);
380 ERI_Spherical_Harmonics(l, -m, theta, phi, Y2, dY2dt, dY2dp);
381
382 ERI_RCSH_Coeff(m, a, b);
383
384 if (Z) {
385 *Z = a[0]*Y1[0] - a[1]*Y1[1] + b[0]*Y2[0] - b[1]*Y2[1];
386 }
387 if (dZdt) {
388 *dZdt = a[0]*dY1dt[0] - a[1]*dY1dt[1] + b[0]*dY2dt[0] - b[1]*dY2dt[1];
389 }
390 if (dZdp) {
391 *dZdp = a[0]*dY1dp[0] - a[1]*dY1dp[1] + b[0]*dY2dp[0] - b[1]*dY2dp[1];
392 }
393 }
394
395
396
397
398 /*----------------------------------------------------------------------
399 ERI_Associated_Legendre_Polynomial
400
401 Associated Legendre Polynomial and its derivative.
402 ----------------------------------------------------------------------*/
ERI_Associated_Legendre_Polynomial(int l,int m,double x,double * P,double * dP)403 void ERI_Associated_Legendre_Polynomial(
404 int l, /* (IN) positive integer */
405 int m, /* (IN) integer 0 < m < l */
406 double x, /* (IN) argument -1 <= x <= 1 */
407 double *P, /* (OUT) function value */
408 double *dP /* (OUT) derivative */
409 )
410 {
411 const double cut0 = 1.0e-24;
412 const double cut1 = 1.0e-12;
413 double Pm, Pm1, f, p0, p1, tmp;
414 int i, ll;
415
416 assert( fabs(x) <= 1.0 );
417 assert( l >= 0 && m >= 0 && m <= l );
418
419 if ((1.0-cut0)<fabs(x)){
420 /* x = sgn(x)*(1.0-cut0); */
421 if (x < 0.0) {
422 x = -(1.0-cut0);
423 } else {
424 x = (1.0-cut0);
425 }
426 }
427
428 /* calculate Pm */
429
430 /* (-1)^m (2m-1)!! (1-x^2)^(m/2) */
431 Pm = 1.0;
432 if (m>0){
433 f = 1.0;
434 tmp = -sqrt((1.0-x)*(1.0+x));
435 for (i=1; i<=m; i++){
436 Pm *= f*tmp;
437 f += 2.0;
438 }
439 }
440
441 if (l==m){
442 /* m = l */
443 p0 = Pm;
444 p1 = 0.0;
445 } else{
446 /* m < l */
447 Pm1 = x*(2.0*(double)m + 1.0)*Pm;
448 if (m<l-1) {
449 /* m < l-1 */
450 for (ll=m+2; ll<=l; ll++) {
451 tmp = (x*(2.0*(double)ll-1.0)*Pm1 -
452 ((double)ll+(double)m-1.0)*Pm)/(double)(ll-m);
453 Pm = Pm1;
454 Pm1 = tmp;
455 }
456 } /* end if */
457 p0 = Pm1;
458 p1 = Pm;
459 } /* end if */
460
461 *P = p0;
462
463 *dP = 0.0;
464 tmp = sqrt(1.0-x*x);
465 if (cut1<tmp) {
466 *dP = ((double)l*x*p0-(double)(l+m)*p1)/tmp;
467 }
468 }
469
470
471
472
473 /*----------------------------------------------------------------------
474 ERI_Associated_Laguerre_Polynomial
475
476 The associated Laguerre polynomials in the Rodridues representation.
477 ----------------------------------------------------------------------*/
ERI_Associated_Laguerre_Polynomial(double x,int n,double alpha,double * L,double * dL)478 void ERI_Associated_Laguerre_Polynomial(
479 double x, /* (IN) argument */
480 int n, /* (IN) zero or positive integer */
481 double alpha, /* (IN) */
482 double *L, /* (OUT) function value */
483 double *dL /* (OUT) derivative */
484 )
485 {
486 int i;
487 double L0, L1, L2;
488
489 assert( n >= 0 );
490
491 if (n==0) {
492 *L = 1.0;
493 *dL = 0.0;
494 return;
495 } else if (n==1) {
496 *L = (1.0 + alpha)-x;
497 *dL = -1.0;
498 return;
499 }
500
501 /* n > 1 */
502 L1 = 1.0; /* L_0 */
503 L2 = (1.0+alpha)-x; /* L_1 */
504
505 /* recurrence */
506 for (i=1; i<n; i++) {
507 L0 = L1; /* L_{i-1} */
508 L1 = L2; /* L_i */
509 L2 = (((double)(2*i+1)+(alpha-x))*L1
510 - ((double)i+alpha)*L0)/(double)(i+1); /* L_{i+1} */
511 }
512 *L = L2;
513 *dL = ((double)n*L2-((double)n+alpha)*L1)/x;
514 }
515
516
517
518
519 /*----------------------------------------------------------------------
520 ERI_Gaunt
521
522 Routine to calculate the Gaunt coefficients.
523
524 See Eq. (3.7.73) in Modern Quantum Mechanics by J. J. Sakurai.
525 ----------------------------------------------------------------------*/
ERI_Gaunt(int l,int m,int l1,int m1,int l2,int m2)526 double ERI_Gaunt(int l, int m, int l1, int m1, int l2, int m2)
527 {
528 double tmp0, tmp1, tmp2, tmp3;
529 double result, cg1, cg2;
530
531 cg1 = CG(l1, 0, l2, 0, l, 0);
532 cg2 = CG(l1, m1, l2, m2, l, m);
533
534 return cg1*cg2*sqrt( (double)((2*l1+1)*(2*l2+1))/(double)(8*l+4)/PI );
535 }
536
537
538 /*----------------------------------------------------------------------
539 ERI_Gaunt_R
540
541 Gaunt coefficients for real spherical harmonics.
542
543 GC_R(L,L1,L2) = int Z(L) Z(L1) Z(L2)
544 ----------------------------------------------------------------------*/
ERI_Gaunt_R(int l,int m,int l1,int m1,int l2,int m2)545 double ERI_Gaunt_R(int l, int m, int l1, int m1, int l2, int m2)
546 {
547 double a[2], b[2], a1[2], b1[2], a2[2], b2[2];
548 double aaa, aab, aba, abb, baa, bab, bba, bbb, s;
549
550 ERI_RCSH_Coeff( m, a, b);
551 ERI_RCSH_Coeff(m1, a1, b1);
552 ERI_RCSH_Coeff(m2, a2, b2);
553
554 s = (0==m%2) ? 1.0 : -1.0;
555
556 aaa = (a[0]*a1[0]-a[1]*a1[1])*a2[0] - (a[0]*a1[1]+a[1]*a1[0])*a2[1];
557 aab = (a[0]*a1[0]-a[1]*a1[1])*b2[0] - (a[0]*a1[1]+a[1]*a1[0])*b2[1];
558 aba = (a[0]*b1[0]-a[1]*b1[1])*a2[0] - (a[0]*b1[1]+a[1]*b1[0])*a2[1];
559 abb = (a[0]*b1[0]-a[1]*b1[1])*b2[0] - (a[0]*b1[1]+a[1]*b1[0])*b2[1];
560 baa = (b[0]*a1[0]-b[1]*a1[1])*a2[0] - (b[0]*a1[1]+b[1]*a1[0])*a2[1];
561 bab = (b[0]*a1[0]-b[1]*a1[1])*b2[0] - (b[0]*a1[1]+b[1]*a1[0])*b2[1];
562 bba = (b[0]*b1[0]-b[1]*b1[1])*a2[0] - (b[0]*b1[1]+b[1]*b1[0])*a2[1];
563 bbb = (b[0]*b1[0]-b[1]*b1[1])*b2[0] - (b[0]*b1[1]+b[1]*b1[0])*b2[1];
564
565 return s*( aaa*ERI_Gaunt(l, -m, l1, m1, l2, m2)
566 + aab*ERI_Gaunt(l, -m, l1, m1, l2, -m2)
567 + aba*ERI_Gaunt(l, -m, l1, -m1, l2, m2)
568 + abb*ERI_Gaunt(l, -m, l1, -m1, l2, -m2)
569 + baa*ERI_Gaunt(l, m, l1, m1, l2, m2)
570 + bab*ERI_Gaunt(l, m, l1, m1, l2, -m2)
571 + bba*ERI_Gaunt(l, m, l1, -m1, l2, m2)
572 + bbb*ERI_Gaunt(l, m, l1, -m1, l2, -m2) );
573 }
574
575
576
577
578
579 /*----------------------------------------------------------------------
580 Clebsch-Gordan coefficients
581
582 <j1 j2; m1 m2 | j1 j2; j m>
583 =
584 delta(m,m1+m2)*sqrt(2j+1)
585 *sqrt{ (j1+j2-j)! (j+j1-j2)! (j+j2-j1)! / (j+j1+j2+1)! }
586 *sqrt{ (j1+m1)! (j1-m1)! (j2+m2)! (j2-m2)! (j+m)! (j-m!) }
587 *sum_k (-1)^k / ( k! (j1+j2-j-k)! (j1-m1-k)! (j2+m2-k)!
588 (j-j2+m1+k)! (j-j1-m2+k)! )
589
590 See M. Avbramowitz and I.A. Stegun (eds.), "Handbook of Mathematical
591 Functions with Formulas, Graphs, and Tables,"
592 Dover Publications, Inc. (Mineola, NY), 1972; an electronic copy of
593 this book is available at http://www.math.sfu.ca/~cbm/aands/.
594 ----------------------------------------------------------------------*/
CG(int j1,int m1,int j2,int m2,int j,int m)595 static double CG(int j1, int m1, int j2, int m2, int j, int m)
596 {
597 int kmin, kmax, k;
598 double sgn, cg;
599 const int fmax = 3*CG_JMAX;
600 double f[3*CG_JMAX];
601
602 /* this routine safely calculates when j1, j2, j < CG_JMAX */
603 if (j1>CG_JMAX || j2>CG_JMAX || j>CG_JMAX) {
604 fprintf(stderr, "*** out of bound (%s, %d)\n", __FILE__, __LINE__);
605 abort();
606 }
607
608 /* factorials */
609 f[0] = f[1] = 1.0;
610 for (k=2; k<fmax; k++) { f[k] = f[k-1]*(double)k; }
611
612 /* delta(m,m1+m2) */
613 if (m != m1+m2) return 0.0;
614
615 /* conditions for j1, j2, and j */
616 if (j1+j2 < j || j2+j < j1 || j+j1 < j2) return 0.0;
617
618 /* conditions for m1, m2 and m */
619 if (abs(m1)>j1 || abs(m2)>j2 || abs(m)>j) return 0.0;
620
621 /* determin the range of k
622 max(0, -j+j2-m1, -j+j1+m2) <= k <= min(j1+j2-j, j1-m1, j2+m2) */
623 kmin = 0;
624 if (kmin < -j+j2-m1) kmin = -j+j2-m1;
625 if (kmin < -j+j1+m2) kmin = -j+j1+m2;
626
627 kmax = j1+j2-j;
628 if (kmax > j1-m1) kmax = j1-m1;
629 if (kmax > j2+m2) kmax = j2+m2;
630
631 if (kmin>kmax) return 0.0;
632
633 cg = 0.0;
634 sgn = kmin%2==0 ? 1.0 : -1.0;
635 for (k=kmin; k<=kmax; k++) {
636 cg += sgn/(f[k]*f[j1+j2-j-k]*f[j1-m1-k]*f[j2+m2-k]
637 *f[j-j2+m1+k]*f[j-j1-m2+k]);
638 sgn = -sgn;
639 }
640
641 cg *= sqrt( (double)(2*j+1)*f[j1+j2-j]*f[j+j1-j2]*f[j+j2-j1]
642 *f[j1+m1]*f[j1-m1]*f[j2+m2]*f[j2-m2]*f[j+m]*f[j-m]
643 /f[j+j1+j2+1] );
644
645 return cg;
646 }
647
648
649
ERI_GLQ(double * x,double * w,int n)650 void ERI_GLQ(double *x, double *w, int n)
651 {
652 int i, itr;
653 double err, xi, i1, L, dL;
654
655 for (i=0; i<n; i++) {
656
657 /* Initial guess for the roots of Laguerre polynomials
658 See A. Stroud and D. Secrest, "Gaussian Quadrature
659 Formulas," Prentice-Hall, Englewood, NJ, 1966. */
660 if (i==0) {
661 xi = 3.0/(1.0+2.4*(double)n);
662 } else if (i==1) {
663 xi = x[i-1] + 15.0/(1.0+2.5*(double)n);
664 } else {
665 i1 = (double)(i-1);
666 xi = x[i-1] + (x[i-1]-x[i-2])*(1.0+2.55*i1)/(1.9*i1);
667 }
668
669 /* Find the roots by Newton's method */
670 err = DBL_MAX;
671 for (itr=0; (itr<MAX_ITR) && (fabs(err)>ERR_MAX); itr++) {
672 ERI_Associated_Laguerre_Polynomial(xi, n, 0.0, &L, &dL);
673 err = L/dL;
674 xi -= err;
675 }
676
677 #ifndef NDEBUG
678 if (fabs(err)>ERR_MAX) {
679 fprintf(stderr, "*** error in ERI_GLQ: no convergence\n");
680 fprintf(stderr, " err = %12.4e\n", err);
681 }
682 #endif
683
684 x[i] = xi;
685 w[i] = 1.0/dL/dL/xi;
686 } /* end loop of i */
687 }
688
689
690