REAL(KIND=DP) FUNCTION cleb (ja, jb, jc, ma, mb, mc) USE Numeric_Kinds_Module IMPLICIT NONE INTEGER ja, jb, jc, ma, mb, mc, j1, j2, j3, m1, m2, m3, jx, mx INTEGER j1max, j1min, j2old, j3old, m1old, m2old, ncoef, mid INTEGER n2, mid1, i, ip1, im1, j, j3m3 !---------------------------------------------------------------------- ! this program was written by g. a. parker 1/3/77. for a ! description fo the algorithm used see klaus schulten and r. g. ! gordon jmp 16,1961,(1975). ! this program is most efficient when a series of clebsch- ! gordan coefficients are desired !---------------------------------------------------------------------- REAL(KIND=DP) f(500), sj2mj3, sj2j31, xj2m, xj3m, xm1, y REAL(KIND=DP) sj1, xj1, xm, x, xj1p1, sj1p1, g, coef, sum, sig, sm1 DATA j2old, j3old, m1old, m2old / -1000000, -1000000, -10000000, -1000000 / !---------------------------------------------------------------------- ! statement functions used in recurrence relations. !---------------------------------------------------------------------- !a(sj1) = sqrt ( (sj1 - sj2mj3) * (sj2j31 - sj1) * (sj1 - sm1) ) !b(xj1) = - (2.0d0 * xj1 + 1) * (xj2m - xj3m - xj1 * (xj1 + 1) * xm1) !---------------------------------------------------------------------- ! check triangle inequalities. !---------------------------------------------------------------------- IF(iabs (mc) .gt.jc) goto 130 IF(iabs (mb) .gt.jb) goto 130 IF(iabs (ma) .gt.ja) goto 130 IF(mc.ne.ma + mb) goto 130 !---------------------------------------------------------------------- ! calculate maximum and minimum allowed values of j1 !---------------------------------------------------------------------- j1 = ja j2 = jb j3 = jc m1 = ma m2 = mb m3 = - mc jx = iabs (j2 - j3) mx = iabs (m1) j1max = j2 + j3 j1min = max0 (jx, mx) IF(j1max - j1min + 1.ge.500) goto 140 IF(j1.gt.j1max.or.j1.lt.j1min) goto 130 !---------------------------------------------------------------------- ! check to see if the clebsch-gordan coefficient has already ! determined. !---------------------------------------------------------------------- IF(j2old.ne.j2) goto 10 IF(j3old.ne.j3) goto 10 IF(m1old.ne.m1) goto 10 IF(m2old.ne.m2) goto 10 goto 100 !---------------------------------------------------------------------- ! coefficients have not been previously determined. !---------------------------------------------------------------------- 10 j2old = j2 j3old = j3 m1old = m1 m2old = m2 !---------------------------------------------------------------------- ! determine constants for statement functions. !---------------------------------------------------------------------- sj2mj3 = jx * jx sj2j31 = (j1max + 1) * (j1max + 1) sm1 = mx * mx xj2m = j2 * (j2 + 1) * m1 xj3m = j3 * (j3 + 1) * m1 xm1 = - m1 - m2 - m2 !---------------------------------------------------------------------- ! determine where the middle coefficient is located and the ! number of coefficients. !---------------------------------------------------------------------- ncoef = j1max - j1min + 1 mid = ncoef / 2 + 1 IF(ncoef.le.2) goto 110 !---------------------------------------------------------------------- ! set initial values for the coefficietts f(1) and f(ncoef). !---------------------------------------------------------------------- f (1) = 1.d0 f (ncoef) = 1.d0 !---------------------------------------------------------------------- ! use equations 16 and 17 from the reference to determine ! f(2) and f(ncoef-1) !---------------------------------------------------------------------- IF(j1min.ne.0) goto 20 xm = 2 * m2 x = (2 * j2 + 1) **2 - 1 f (2) = - f (1) * xm / sqrt (x) goto 30 20 xj1 = j1min xj1p1 = j1min + 1 sj1p1 = xj1p1 * xj1p1 f (2) = f (1) * b (xj1) / (xj1 * a (sj1p1) ) 30 g = f (2) n2 = ncoef - 1 xj1 = j1max sj1 = xj1 * xj1 f (n2) = f (ncoef) * b (xj1) / ( (xj1 + 1.0d0) * a (sj1) ) mid1 = mid-1 IF(ncoef.eq.3) goto 50 !---------------------------------------------------------------------- ! start recursive evaluation of the coefficients using ! equation 6'. !---------------------------------------------------------------------- do i = 2, mid1 ip1 = i + 1 im1 = i - 1 xj1 = j1min + i - 1 xj1p1 = xj1 + 1.0d0 sj1 = xj1 * xj1 sj1p1 = xj1p1 * xj1p1 g = (b (xj1) * f (i) - xj1p1 * a (sj1) * f (im1) ) / (xj1 * a ( & sj1p1) ) n2 = ncoef - i + 1 IF(n2.eq.mid) goto 50 f (ip1) = g xj1 = j1max - i + 1 xj1p1 = xj1 + 1.0d0 ip1 = n2 + 1 im1 = n2 - 1 sj1 = xj1 * xj1 sj1p1 = xj1p1 * xj1p1 f (im1) = ( - xj1 * a (sj1p1) * f (ip1) + b (xj1) * f (n2) ) / (xj1p1 * a (sj1) ) ENDDO !---------------------------------------------------------------------- ! match at the midpoint. !---------------------------------------------------------------------- 50 IF(dabs (g) .gt.1.d-8) goto 60 i = mid im1 = mid1 xj1 = j1min + mid1 xj1p1 = xj1 + 1.0d0 sj1 = xj1 * xj1 sj1p1 = xj1p1 * xj1p1 g = (b (xj1) * f (i) - xj1p1 * a (sj1) * f (im1) ) / (xj1 * a ( & sj1p1) ) mid1 = mid mid = mid+1 60 coef = f (mid) / g do i = 1, mid1 f (i) = coef * f (i) ENDDO !---------------------------------------------------------------------- ! find magnitude of the normalization coefficient. !---------------------------------------------------------------------- sum = 0.0d0 x = j1min * 2 - 1 do i = 1, ncoef x = x + 2.0d0 sum = sum + x * f (i) * f (i) ENDDO coef = sqrt (1.0d0 / sum) !---------------------------------------------------------------------- ! determine sign of the normalization coefficient. !---------------------------------------------------------------------- sig = 1.0d0 j = j2 - j3 - m1 IF( (j / 2) * 2.ne.j) sig = - 1.0d0 coef = dsign (coef, sig) !---------------------------------------------------------------------- ! normalize to obtain the clebsch-gordan coefficients. !---------------------------------------------------------------------- y = 2 * j3 + 1 y = sqrt (y) j3m3 = j3 + m3 IF( (j3m3 / 2) * 2.ne.j3m3) y = - y do i = 1, ncoef f (i) = coef * f (i) * y ENDDO 100 i = j1 - j1min + 1 cleb = f (i) RETURN 110 IF(ncoef.eq.2) goto 120 !---------------------------------------------------------------------- ! special case for only 1 coefficient. !---------------------------------------------------------------------- y = 2 * j3 + 1 x = 2 * j1 + 1 j = j2 - m2 sig = 1.0d0 IF( (j / 2) * 2.ne.j) sig = - 1.0d0 f (1) = sqrt (y / x) * sig goto 100 !---------------------------------------------------------------------- ! special case for only 2 coefficients. !---------------------------------------------------------------------- 120 f (1) = - 1.0d0 xj1 = j1min xj1p1 = j1min + 1 sj1p1 = xj1p1 * xj1p1 f (2) = f (1) * b (xj1) / (xj1 * a (sj1p1) ) x = j1min * 2 + 1 sum = x * f (1) * f (1) + (x + 2.0d0) * f (2) * f (2) coef = sqrt (1.0d0 / sum) sig = 1.0d0 j = j2 - j3 - m1 IF( (j / 2) * 2.ne.j) sig = - 1.0d0 sig = sig * f (2) coef = dsign(coef, sig) y = 2 * j3 + 1 y = sqrt (y) j3m3 = j3 + m3 IF( (j3m3 / 2) * 2.ne.j3m3) y = - y f (1) = coef * y * f (1) f (2) = coef * y * f (2) goto 100 !---------------------------------------------------------------------- ! triangle inequalities were not satisfied. !---------------------------------------------------------------------- 130 cleb = 0.0d0 RETURN 140 WRITE(6, 150) j1min, j1max 150 FORMAT(/,'error in cleb. the number of j1 values exceeds', & & ' dimensions.',5x,'j1min=',5x,i5,'j1max=',5x,i5) !----------------***end-cleb***--------------------------------------- CONTAINS REAL(KIND=DP) FUNCTION A(sj1) USE Numeric_Kinds_Module REAL(KIND=DP) sj1 a = sqrt ( (sj1 - sj2mj3) * (sj2j31 - sj1) * (sj1 - sm1) ) END FUNCTION A REAL(KIND=DP) FUNCTION B(xj1) USE Numeric_Kinds_Module REAL(KIND=DP) xj1 b = - (2.0d0 * xj1 + 1) * (xj2m - xj3m - xj1 * (xj1 + 1) * xm1) END FUNCTION B END FUNCTION cleb