*DECK DRC3JM SUBROUTINE DRC3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM, + IER) C***BEGIN PROLOGUE DRC3JM C***PURPOSE Evaluate the 3j symbol g(M2) = (L1 L2 L3 ) C (M1 M2 -M1-M2) C for all allowed values of M2, the other parameters C being held fixed. C***LIBRARY SLATEC C***CATEGORY C19 C***TYPE DOUBLE PRECISION (RC3JM-S, DRC3JM-D) C***KEYWORDS 3J COEFFICIENTS, 3J SYMBOLS, CLEBSCH-GORDAN COEFFICIENTS, C RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS, C WIGNER COEFFICIENTS C***AUTHOR Gordon, R. G., Harvard University C Schulten, K., Max Planck Institute C***DESCRIPTION C C *Usage: C C DOUBLE PRECISION L1, L2, L3, M1, M2MIN, M2MAX, THRCOF(NDIM) C INTEGER NDIM, IER C C CALL DRC3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM, IER) C C *Arguments: C C L1 :IN Parameter in 3j symbol. C C L2 :IN Parameter in 3j symbol. C C L3 :IN Parameter in 3j symbol. C C M1 :IN Parameter in 3j symbol. C C M2MIN :OUT Smallest allowable M2 in 3j symbol. C C M2MAX :OUT Largest allowable M2 in 3j symbol. C C THRCOF :OUT Set of 3j coefficients generated by evaluating the C 3j symbol for all allowed values of M2. THRCOF(I) C will contain g(M2MIN+I-1), I=1,2,...,M2MAX-M2MIN+1. C C NDIM :IN Declared length of THRCOF in calling program. C C IER :OUT Error flag. C IER=0 No errors. C IER=1 Either L1.LT.ABS(M1) or L1+ABS(M1) non-integer. C IER=2 ABS(L1-L2).LE.L3.LE.L1+L2 not satisfied. C IER=3 L1+L2+L3 not an integer. C IER=4 M2MAX-M2MIN not an integer. C IER=5 M2MAX less than M2MIN. C IER=6 NDIM less than M2MAX-M2MIN+1. C C *Description: C C Although conventionally the parameters of the vector addition C coefficients satisfy certain restrictions, such as being integers C or integers plus 1/2, the restrictions imposed on input to this C subroutine are somewhat weaker. See, for example, Section 27.9 of C Abramowitz and Stegun or Appendix C of Volume II of A. Messiah. C The restrictions imposed by this subroutine are C 1. L1.GE.ABS(M1) and L1+ABS(M1) must be an integer; C 2. ABS(L1-L2).LE.L3.LE.L1+L2; C 3. L1+L2+L3 must be an integer; C 4. M2MAX-M2MIN must be an integer, where C M2MAX=MIN(L2,L3-M1) and M2MIN=MAX(-L2,-L3-M1). C If the conventional restrictions are satisfied, then these C restrictions are met. C C The user should be cautious in using input parameters that do C not satisfy the conventional restrictions. For example, the C the subroutine produces values of C g(M2) = (0.75 1.50 1.75 ) C (0.25 M2 -0.25-M2) C for M2=-1.5,-0.5,0.5,1.5 but none of the symmetry properties of the C 3j symbol, set forth on page 1056 of Messiah, is satisfied. C C The subroutine generates g(M2MIN), g(M2MIN+1), ..., g(M2MAX) C where M2MIN and M2MAX are defined above. The sequence g(M2) is C generated by a three-term recurrence algorithm with scaling to C control overflow. Both backward and forward recurrence are used to C maintain numerical stability. The two recurrence sequences are C matched at an interior point and are normalized from the unitary C property of 3j coefficients and Wigner's phase convention. C C The algorithm is suited to applications in which large quantum C numbers arise, such as in molecular dynamics. C C***REFERENCES 1. Abramowitz, M., and Stegun, I. A., Eds., Handbook C of Mathematical Functions with Formulas, Graphs C and Mathematical Tables, NBS Applied Mathematics C Series 55, June 1964 and subsequent printings. C 2. Messiah, Albert., Quantum Mechanics, Volume II, C North-Holland Publishing Company, 1963. C 3. Schulten, Klaus and Gordon, Roy G., Exact recursive C evaluation of 3j and 6j coefficients for quantum- C mechanical coupling of angular momenta, J Math C Phys, v 16, no. 10, October 1975, pp. 1961-1970. C 4. Schulten, Klaus and Gordon, Roy G., Semiclassical C approximations to 3j and 6j coefficients for C quantum-mechanical coupling of angular momenta, C J Math Phys, v 16, no. 10, October 1975, C pp. 1971-1988. C 5. Schulten, Klaus and Gordon, Roy G., Recursive C evaluation of 3j and 6j coefficients, Computer C Phys Comm, v 11, 1976, pp. 269-278. C***ROUTINES CALLED D1MACH, XERMSG C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 880515 SLATEC prologue added by G. C. Nielson, NBS; parameters C HUGE and TINY revised to depend on D1MACH. C 891229 Prologue description rewritten; other prologue sections C revised; MMATCH (location of match point for recurrences) C removed from argument list; argument IER changed to serve C only as an error flag (previously, in cases without error, C it returned the number of scalings); number of error codes C increased to provide more precise error information; C program comments revised; SLATEC error handler calls C introduced to enable printing of error messages to meet C SLATEC standards. These changes were done by D. W. Lozier, C M. A. McClain and J. M. Smith of the National Institute C of Standards and Technology, formerly NBS. C 910415 Mixed type expressions eliminated; variable C1 initialized; C description of THRCOF expanded. These changes were done by C D. W. Lozier. C***END PROLOGUE DRC3JM C INTEGER NDIM, IER DOUBLE PRECISION L1, L2, L3, M1, M2MIN, M2MAX, THRCOF(NDIM) C INTEGER I, INDEX, LSTEP, N, NFIN, NFINP1, NFINP2, NFINP3, NLIM, + NSTEP2 DOUBLE PRECISION A1, A1S, C1, C1OLD, C2, CNORM, D1MACH, DV, EPS, + HUGE, M2, M3, NEWFAC, OLDFAC, ONE, RATIO, SIGN1, + SIGN2, SRHUGE, SRTINY, SUM1, SUM2, SUMBAC, + SUMFOR, SUMUNI, THRESH, TINY, TWO, X, X1, X2, X3, + Y, Y1, Y2, Y3, ZERO C DATA ZERO,EPS,ONE,TWO /0.0D0,0.01D0,1.0D0,2.0D0/ C C***FIRST EXECUTABLE STATEMENT DRC3JM IER=0 C HUGE is the square root of one twentieth of the largest floating C point number, approximately. HUGE = SQRT(D1MACH(2)/20.0D0) SRHUGE = SQRT(HUGE) TINY = 1.0D0/HUGE SRTINY = 1.0D0/SRHUGE C C MMATCH = ZERO C C C Check error conditions 1, 2, and 3. IF((L1-ABS(M1)+EPS.LT.ZERO).OR. + (MOD(L1+ABS(M1)+EPS,ONE).GE.EPS+EPS))THEN IER=1 CALL XERMSG('SLATEC','DRC3JM','L1-ABS(M1) less than zero or '// + 'L1+ABS(M1) not integer.',IER,1) RETURN ELSEIF((L1+L2-L3.LT.-EPS).OR.(L1-L2+L3.LT.-EPS).OR. + (-L1+L2+L3.LT.-EPS))THEN IER=2 CALL XERMSG('SLATEC','DRC3JM','L1, L2, L3 do not satisfy '// + 'triangular condition.',IER,1) RETURN ELSEIF(MOD(L1+L2+L3+EPS,ONE).GE.EPS+EPS)THEN IER=3 CALL XERMSG('SLATEC','DRC3JM','L1+L2+L3 not integer.',IER,1) RETURN ENDIF C C C Limits for M2 M2MIN = MAX(-L2,-L3-M1) M2MAX = MIN(L2,L3-M1) C C Check error condition 4. IF(MOD(M2MAX-M2MIN+EPS,ONE).GE.EPS+EPS)THEN IER=4 CALL XERMSG('SLATEC','DRC3JM','M2MAX-M2MIN not integer.',IER,1) RETURN ENDIF IF(M2MIN.LT.M2MAX-EPS) GO TO 20 IF(M2MIN.LT.M2MAX+EPS) GO TO 10 C C Check error condition 5. IER=5 CALL XERMSG('SLATEC','DRC3JM','M2MIN greater than M2MAX.',IER,1) RETURN C C C This is reached in case that M2 and M3 can take only one value. 10 CONTINUE C MSCALE = 0 THRCOF(1) = (-ONE) ** INT(ABS(L2-L3-M1)+EPS) / 1 SQRT(L1+L2+L3+ONE) RETURN C C This is reached in case that M1 and M2 take more than one value. 20 CONTINUE C MSCALE = 0 NFIN = INT(M2MAX-M2MIN+ONE+EPS) IF(NDIM-NFIN) 21, 23, 23 C C Check error condition 6. 21 IER = 6 CALL XERMSG('SLATEC','DRC3JM','Dimension of result array for '// + '3j coefficients too small.',IER,1) RETURN C C C C Start of forward recursion from M2 = M2MIN C 23 M2 = M2MIN THRCOF(1) = SRTINY NEWFAC = 0.0D0 C1 = 0.0D0 SUM1 = TINY C C LSTEP = 1 30 LSTEP = LSTEP + 1 M2 = M2 + ONE M3 = - M1 - M2 C C OLDFAC = NEWFAC A1 = (L2-M2+ONE) * (L2+M2) * (L3+M3+ONE) * (L3-M3) NEWFAC = SQRT(A1) C C DV = (L1+L2+L3+ONE)*(L2+L3-L1) - (L2-M2+ONE)*(L3+M3+ONE) 1 - (L2+M2-ONE)*(L3-M3-ONE) C IF(LSTEP-2) 32, 32, 31 C 31 C1OLD = ABS(C1) 32 C1 = - DV / NEWFAC C IF(LSTEP.GT.2) GO TO 60 C C C If M2 = M2MIN + 1, the third term in the recursion equation vanishes, C hence C X = SRTINY * C1 THRCOF(2) = X SUM1 = SUM1 + TINY * C1*C1 IF(LSTEP.EQ.NFIN) GO TO 220 GO TO 30 C C 60 C2 = - OLDFAC / NEWFAC C C Recursion to the next 3j coefficient X = C1 * THRCOF(LSTEP-1) + C2 * THRCOF(LSTEP-2) THRCOF(LSTEP) = X SUMFOR = SUM1 SUM1 = SUM1 + X*X IF(LSTEP.EQ.NFIN) GO TO 100 C C See if last unnormalized 3j coefficient exceeds SRHUGE C IF(ABS(X).LT.SRHUGE) GO TO 80 C C This is reached if last 3j coefficient larger than SRHUGE, C so that the recursion series THRCOF(1), ... , THRCOF(LSTEP) C has to be rescaled to prevent overflow C C MSCALE = MSCALE + 1 DO 70 I=1,LSTEP IF(ABS(THRCOF(I)).LT.SRTINY) THRCOF(I) = ZERO 70 THRCOF(I) = THRCOF(I) / SRHUGE SUM1 = SUM1 / HUGE SUMFOR = SUMFOR / HUGE X = X / SRHUGE C C C As long as ABS(C1) is decreasing, the recursion proceeds towards C increasing 3j values and, hence, is numerically stable. Once C an increase of ABS(C1) is detected, the recursion direction is C reversed. C 80 IF(C1OLD-ABS(C1)) 100, 100, 30 C C C Keep three 3j coefficients around MMATCH for comparison later C with backward recursion values. C 100 CONTINUE C MMATCH = M2 - 1 NSTEP2 = NFIN - LSTEP + 3 X1 = X X2 = THRCOF(LSTEP-1) X3 = THRCOF(LSTEP-2) C C Starting backward recursion from M2MAX taking NSTEP2 steps, so C that forwards and backwards recursion overlap at the three points C M2 = MMATCH+1, MMATCH, MMATCH-1. C NFINP1 = NFIN + 1 NFINP2 = NFIN + 2 NFINP3 = NFIN + 3 THRCOF(NFIN) = SRTINY SUM2 = TINY C C C M2 = M2MAX + TWO LSTEP = 1 110 LSTEP = LSTEP + 1 M2 = M2 - ONE M3 = - M1 - M2 OLDFAC = NEWFAC A1S = (L2-M2+TWO) * (L2+M2-ONE) * (L3+M3+TWO) * (L3-M3-ONE) NEWFAC = SQRT(A1S) DV = (L1+L2+L3+ONE)*(L2+L3-L1) - (L2-M2+ONE)*(L3+M3+ONE) 1 - (L2+M2-ONE)*(L3-M3-ONE) C1 = - DV / NEWFAC IF(LSTEP.GT.2) GO TO 120 C C If M2 = M2MAX + 1 the third term in the recursion equation vanishes C Y = SRTINY * C1 THRCOF(NFIN-1) = Y IF(LSTEP.EQ.NSTEP2) GO TO 200 SUMBAC = SUM2 SUM2 = SUM2 + Y*Y GO TO 110 C 120 C2 = - OLDFAC / NEWFAC C C Recursion to the next 3j coefficient C Y = C1 * THRCOF(NFINP2-LSTEP) + C2 * THRCOF(NFINP3-LSTEP) C IF(LSTEP.EQ.NSTEP2) GO TO 200 C THRCOF(NFINP1-LSTEP) = Y SUMBAC = SUM2 SUM2 = SUM2 + Y*Y C C C See if last 3j coefficient exceeds SRHUGE C IF(ABS(Y).LT.SRHUGE) GO TO 110 C C This is reached if last 3j coefficient larger than SRHUGE, C so that the recursion series THRCOF(NFIN), ... , THRCOF(NFIN-LSTEP+1) C has to be rescaled to prevent overflow. C C MSCALE = MSCALE + 1 DO 111 I=1,LSTEP INDEX = NFIN - I + 1 IF(ABS(THRCOF(INDEX)).LT.SRTINY) 1 THRCOF(INDEX) = ZERO 111 THRCOF(INDEX) = THRCOF(INDEX) / SRHUGE SUM2 = SUM2 / HUGE SUMBAC = SUMBAC / HUGE C GO TO 110 C C C C The forward recursion 3j coefficients X1, X2, X3 are to be matched C with the corresponding backward recursion values Y1, Y2, Y3. C 200 Y3 = Y Y2 = THRCOF(NFINP2-LSTEP) Y1 = THRCOF(NFINP3-LSTEP) C C C Determine now RATIO such that YI = RATIO * XI (I=1,2,3) holds C with minimal error. C RATIO = ( X1*Y1 + X2*Y2 + X3*Y3 ) / ( X1*X1 + X2*X2 + X3*X3 ) NLIM = NFIN - NSTEP2 + 1 C IF(ABS(RATIO).LT.ONE) GO TO 211 C DO 210 N=1,NLIM 210 THRCOF(N) = RATIO * THRCOF(N) SUMUNI = RATIO * RATIO * SUMFOR + SUMBAC GO TO 230 C 211 NLIM = NLIM + 1 RATIO = ONE / RATIO DO 212 N=NLIM,NFIN 212 THRCOF(N) = RATIO * THRCOF(N) SUMUNI = SUMFOR + RATIO*RATIO*SUMBAC GO TO 230 C 220 SUMUNI = SUM1 C C C Normalize 3j coefficients C 230 CNORM = ONE / SQRT((L1+L1+ONE) * SUMUNI) C C Sign convention for last 3j coefficient determines overall phase C SIGN1 = SIGN(ONE,THRCOF(NFIN)) SIGN2 = (-ONE) ** INT(ABS(L2-L3-M1)+EPS) IF(SIGN1*SIGN2) 235,235,236 235 CNORM = - CNORM C 236 IF(ABS(CNORM).LT.ONE) GO TO 250 C DO 240 N=1,NFIN 240 THRCOF(N) = CNORM * THRCOF(N) RETURN C 250 THRESH = TINY / ABS(CNORM) DO 251 N=1,NFIN IF(ABS(THRCOF(N)).LT.THRESH) THRCOF(N) = ZERO 251 THRCOF(N) = CNORM * THRCOF(N) C C C RETURN END