1*DECK RC6J 2 SUBROUTINE RC6J (L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF, NDIM, 3 + IER) 4C***BEGIN PROLOGUE RC6J 5C***PURPOSE Evaluate the 6j symbol h(L1) = {L1 L2 L3} 6C {L4 L5 L6} 7C for all allowed values of L1, the other parameters 8C being held fixed. 9C***LIBRARY SLATEC 10C***CATEGORY C19 11C***TYPE SINGLE PRECISION (RC6J-S, DRC6J-D) 12C***KEYWORDS 6J COEFFICIENTS, 6J SYMBOLS, CLEBSCH-GORDAN COEFFICIENTS, 13C RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS, 14C WIGNER COEFFICIENTS 15C***AUTHOR Gordon, R. G., Harvard University 16C Schulten, K., Max Planck Institute 17C***DESCRIPTION 18C 19C *Usage: 20C 21C REAL L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF(NDIM) 22C INTEGER NDIM, IER 23C 24C CALL RC6J(L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF, NDIM, IER) 25C 26C *Arguments: 27C 28C L2 :IN Parameter in 6j symbol. 29C 30C L3 :IN Parameter in 6j symbol. 31C 32C L4 :IN Parameter in 6j symbol. 33C 34C L5 :IN Parameter in 6j symbol. 35C 36C L6 :IN Parameter in 6j symbol. 37C 38C L1MIN :OUT Smallest allowable L1 in 6j symbol. 39C 40C L1MAX :OUT Largest allowable L1 in 6j symbol. 41C 42C SIXCOF :OUT Set of 6j coefficients generated by evaluating the 43C 6j symbol for all allowed values of L1. SIXCOF(I) 44C will contain h(L1MIN+I-1), I=1,2,...,L1MAX-L1MIN+1. 45C 46C NDIM :IN Declared length of SIXCOF in calling program. 47C 48C IER :OUT Error flag. 49C IER=0 No errors. 50C IER=1 L2+L3+L5+L6 or L4+L2+L6 not an integer. 51C IER=2 L4, L2, L6 triangular condition not satisfied. 52C IER=3 L4, L5, L3 triangular condition not satisfied. 53C IER=4 L1MAX-L1MIN not an integer. 54C IER=5 L1MAX less than L1MIN. 55C IER=6 NDIM less than L1MAX-L1MIN+1. 56C 57C *Description: 58C 59C The definition and properties of 6j symbols can be found, for 60C example, in Appendix C of Volume II of A. Messiah. Although the 61C parameters of the vector addition coefficients satisfy certain 62C conventional restrictions, the restriction that they be non-negative 63C integers or non-negative integers plus 1/2 is not imposed on input 64C to this subroutine. The restrictions imposed are 65C 1. L2+L3+L5+L6 and L2+L4+L6 must be integers; 66C 2. ABS(L2-L4).LE.L6.LE.L2+L4 must be satisfied; 67C 3. ABS(L4-L5).LE.L3.LE.L4+L5 must be satisfied; 68C 4. L1MAX-L1MIN must be a non-negative integer, where 69C L1MAX=MIN(L2+L3,L5+L6) and L1MIN=MAX(ABS(L2-L3),ABS(L5-L6)). 70C If all the conventional restrictions are satisfied, then these 71C restrictions are met. Conversely, if input to this subroutine meets 72C all of these restrictions and the conventional restriction stated 73C above, then all the conventional restrictions are satisfied. 74C 75C The user should be cautious in using input parameters that do 76C not satisfy the conventional restrictions. For example, the 77C the subroutine produces values of 78C h(L1) = { L1 2/3 1 } 79C {2/3 2/3 2/3} 80C for L1=1/3 and 4/3 but none of the symmetry properties of the 6j 81C symbol, set forth on pages 1063 and 1064 of Messiah, is satisfied. 82C 83C The subroutine generates h(L1MIN), h(L1MIN+1), ..., h(L1MAX) 84C where L1MIN and L1MAX are defined above. The sequence h(L1) is 85C generated by a three-term recurrence algorithm with scaling to 86C control overflow. Both backward and forward recurrence are used to 87C maintain numerical stability. The two recurrence sequences are 88C matched at an interior point and are normalized from the unitary 89C property of 6j coefficients and Wigner's phase convention. 90C 91C The algorithm is suited to applications in which large quantum 92C numbers arise, such as in molecular dynamics. 93C 94C***REFERENCES 1. Messiah, Albert., Quantum Mechanics, Volume II, 95C North-Holland Publishing Company, 1963. 96C 2. Schulten, Klaus and Gordon, Roy G., Exact recursive 97C evaluation of 3j and 6j coefficients for quantum- 98C mechanical coupling of angular momenta, J Math 99C Phys, v 16, no. 10, October 1975, pp. 1961-1970. 100C 3. Schulten, Klaus and Gordon, Roy G., Semiclassical 101C approximations to 3j and 6j coefficients for 102C quantum-mechanical coupling of angular momenta, 103C J Math Phys, v 16, no. 10, October 1975, 104C pp. 1971-1988. 105C 4. Schulten, Klaus and Gordon, Roy G., Recursive 106C evaluation of 3j and 6j coefficients, Computer 107C Phys Comm, v 11, 1976, pp. 269-278. 108C***ROUTINES CALLED R1MACH, XERMSG 109C***REVISION HISTORY (YYMMDD) 110C 750101 DATE WRITTEN 111C 880515 SLATEC prologue added by G. C. Nielson, NBS; parameters 112C HUGE and TINY revised to depend on R1MACH. 113C 891229 Prologue description rewritten; other prologue sections 114C revised; LMATCH (location of match point for recurrences) 115C removed from argument list; argument IER changed to serve 116C only as an error flag (previously, in cases without error, 117C it returned the number of scalings); number of error codes 118C increased to provide more precise error information; 119C program comments revised; SLATEC error handler calls 120C introduced to enable printing of error messages to meet 121C SLATEC standards. These changes were done by D. W. Lozier, 122C M. A. McClain and J. M. Smith of the National Institute 123C of Standards and Technology, formerly NBS. 124C 910415 Mixed type expressions eliminated; variable C1 initialized; 125C description of SIXCOF expanded. These changes were done by 126C D. W. Lozier. 127C***END PROLOGUE RC6J 128C 129 INTEGER NDIM, IER 130 REAL L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF(NDIM) 131C 132 INTEGER I, INDEX, LSTEP, N, NFIN, NFINP1, NFINP2, NFINP3, NLIM, 133 + NSTEP2 134 REAL A1, A1S, A2, A2S, C1, C1OLD, C2, CNORM, R1MACH, 135 + DENOM, DV, EPS, HUGE, L1, NEWFAC, OLDFAC, ONE, 136 + RATIO, SIGN1, SIGN2, SRHUGE, SRTINY, SUM1, SUM2, 137 + SUMBAC, SUMFOR, SUMUNI, THREE, THRESH, TINY, TWO, 138 + X, X1, X2, X3, Y, Y1, Y2, Y3, ZERO 139C 140 DATA ZERO,EPS,ONE,TWO,THREE /0.0,0.01,1.0,2.0,3.0/ 141C 142C***FIRST EXECUTABLE STATEMENT RC6J 143 IER=0 144C HUGE is the square root of one twentieth of the largest floating 145C point number, approximately. 146 HUGE = SQRT(R1MACH(2)/20.0) 147 SRHUGE = SQRT(HUGE) 148 TINY = 1.0/HUGE 149 SRTINY = 1.0/SRHUGE 150C 151C LMATCH = ZERO 152C 153C Check error conditions 1, 2, and 3. 154 IF((MOD(L2+L3+L5+L6+EPS,ONE).GE.EPS+EPS).OR. 155 + (MOD(L4+L2+L6+EPS,ONE).GE.EPS+EPS))THEN 156 IER=1 157 CALL XERMSG('SLATEC','RC6J','L2+L3+L5+L6 or L4+L2+L6 not '// 158 + 'integer.',IER,1) 159 RETURN 160 ELSEIF((L4+L2-L6.LT.ZERO).OR.(L4-L2+L6.LT.ZERO).OR. 161 + (-L4+L2+L6.LT.ZERO))THEN 162 IER=2 163 CALL XERMSG('SLATEC','RC6J','L4, L2, L6 triangular '// 164 + 'condition not satisfied.',IER,1) 165 RETURN 166 ELSEIF((L4-L5+L3.LT.ZERO).OR.(L4+L5-L3.LT.ZERO).OR. 167 + (-L4+L5+L3.LT.ZERO))THEN 168 IER=3 169 CALL XERMSG('SLATEC','RC6J','L4, L5, L3 triangular '// 170 + 'condition not satisfied.',IER,1) 171 RETURN 172 ENDIF 173C 174C Limits for L1 175C 176 L1MIN = MAX(ABS(L2-L3),ABS(L5-L6)) 177 L1MAX = MIN(L2+L3,L5+L6) 178C 179C Check error condition 4. 180 IF(MOD(L1MAX-L1MIN+EPS,ONE).GE.EPS+EPS)THEN 181 IER=4 182 CALL XERMSG('SLATEC','RC6J','L1MAX-L1MIN not integer.',IER,1) 183 RETURN 184 ENDIF 185 IF(L1MIN.LT.L1MAX-EPS) GO TO 20 186 IF(L1MIN.LT.L1MAX+EPS) GO TO 10 187C 188C Check error condition 5. 189 IER=5 190 CALL XERMSG('SLATEC','RC6J','L1MIN greater than L1MAX.',IER,1) 191 RETURN 192C 193C 194C This is reached in case that L1 can take only one value 195C 196 10 CONTINUE 197C LSCALE = 0 198 SIXCOF(1) = (-ONE) ** INT(L2+L3+L5+L6+EPS) / 199 1 SQRT((L1MIN+L1MIN+ONE)*(L4+L4+ONE)) 200 RETURN 201C 202C 203C This is reached in case that L1 can take more than one value. 204C 205 20 CONTINUE 206C LSCALE = 0 207 NFIN = INT(L1MAX-L1MIN+ONE+EPS) 208 IF(NDIM-NFIN) 21, 23, 23 209C 210C Check error condition 6. 211 21 IER = 6 212 CALL XERMSG('SLATEC','RC6J','Dimension of result array for 6j '// 213 + 'coefficients too small.',IER,1) 214 RETURN 215C 216C 217C Start of forward recursion 218C 219 23 L1 = L1MIN 220 NEWFAC = 0.0 221 C1 = 0.0 222 SIXCOF(1) = SRTINY 223 SUM1 = (L1+L1+ONE) * TINY 224C 225 LSTEP = 1 226 30 LSTEP = LSTEP + 1 227 L1 = L1 + ONE 228C 229 OLDFAC = NEWFAC 230 A1 = (L1+L2+L3+ONE) * (L1-L2+L3) * (L1+L2-L3) * (-L1+L2+L3+ONE) 231 A2 = (L1+L5+L6+ONE) * (L1-L5+L6) * (L1+L5-L6) * (-L1+L5+L6+ONE) 232 NEWFAC = SQRT(A1*A2) 233C 234 IF(L1.LT.ONE+EPS) GO TO 40 235C 236 DV = TWO * ( L2*(L2+ONE)*L5*(L5+ONE) + L3*(L3+ONE)*L6*(L6+ONE) 237 1 - L1*(L1-ONE)*L4*(L4+ONE) ) 238 2 - (L2*(L2+ONE) + L3*(L3+ONE) - L1*(L1-ONE)) 239 3 * (L5*(L5+ONE) + L6*(L6+ONE) - L1*(L1-ONE)) 240C 241 DENOM = (L1-ONE) * NEWFAC 242C 243 IF(LSTEP-2) 32, 32, 31 244C 245 31 C1OLD = ABS(C1) 246 32 C1 = - (L1+L1-ONE) * DV / DENOM 247 GO TO 50 248C 249C If L1 = 1, (L1 - 1) has to be factored out of DV, hence 250C 251 40 C1 = - TWO * ( L2*(L2+ONE) + L5*(L5+ONE) - L4*(L4+ONE) ) 252 1 / NEWFAC 253C 254 50 IF(LSTEP.GT.2) GO TO 60 255C 256C If L1 = L1MIN + 1, the third term in recursion equation vanishes 257C 258 X = SRTINY * C1 259 SIXCOF(2) = X 260 SUM1 = SUM1 + TINY * (L1+L1+ONE) * C1 * C1 261C 262 IF(LSTEP.EQ.NFIN) GO TO 220 263 GO TO 30 264C 265C 266 60 C2 = - L1 * OLDFAC / DENOM 267C 268C Recursion to the next 6j coefficient X 269C 270 X = C1 * SIXCOF(LSTEP-1) + C2 * SIXCOF(LSTEP-2) 271 SIXCOF(LSTEP) = X 272C 273 SUMFOR = SUM1 274 SUM1 = SUM1 + (L1+L1+ONE) * X * X 275 IF(LSTEP.EQ.NFIN) GO TO 100 276C 277C See if last unnormalized 6j coefficient exceeds SRHUGE 278C 279 IF(ABS(X).LT.SRHUGE) GO TO 80 280C 281C This is reached if last 6j coefficient larger than SRHUGE, 282C so that the recursion series SIXCOF(1), ... ,SIXCOF(LSTEP) 283C has to be rescaled to prevent overflow 284C 285C LSCALE = LSCALE + 1 286 DO 70 I=1,LSTEP 287 IF(ABS(SIXCOF(I)).LT.SRTINY) SIXCOF(I) = ZERO 288 70 SIXCOF(I) = SIXCOF(I) / SRHUGE 289 SUM1 = SUM1 / HUGE 290 SUMFOR = SUMFOR / HUGE 291 X = X / SRHUGE 292C 293C 294C As long as the coefficient ABS(C1) is decreasing, the recursion 295C proceeds towards increasing 6j values and, hence, is numerically 296C stable. Once an increase of ABS(C1) is detected, the recursion 297C direction is reversed. 298C 299 80 IF(C1OLD-ABS(C1)) 100, 100, 30 300C 301C 302C Keep three 6j coefficients around LMATCH for comparison later 303C with backward recursion. 304C 305 100 CONTINUE 306C LMATCH = L1 - 1 307 X1 = X 308 X2 = SIXCOF(LSTEP-1) 309 X3 = SIXCOF(LSTEP-2) 310C 311C 312C 313C Starting backward recursion from L1MAX taking NSTEP2 steps, so 314C that forward and backward recursion overlap at the three points 315C L1 = LMATCH+1, LMATCH, LMATCH-1. 316C 317 NFINP1 = NFIN + 1 318 NFINP2 = NFIN + 2 319 NFINP3 = NFIN + 3 320 NSTEP2 = NFIN - LSTEP + 3 321 L1 = L1MAX 322C 323 SIXCOF(NFIN) = SRTINY 324 SUM2 = (L1+L1+ONE) * TINY 325C 326C 327 L1 = L1 + TWO 328 LSTEP = 1 329 110 LSTEP = LSTEP + 1 330 L1 = L1 - ONE 331C 332 OLDFAC = NEWFAC 333 A1S = (L1+L2+L3)*(L1-L2+L3-ONE)*(L1+L2-L3-ONE)*(-L1+L2+L3+TWO) 334 A2S = (L1+L5+L6)*(L1-L5+L6-ONE)*(L1+L5-L6-ONE)*(-L1+L5+L6+TWO) 335 NEWFAC = SQRT(A1S*A2S) 336C 337 DV = TWO * ( L2*(L2+ONE)*L5*(L5+ONE) + L3*(L3+ONE)*L6*(L6+ONE) 338 1 - L1*(L1-ONE)*L4*(L4+ONE) ) 339 2 - (L2*(L2+ONE) + L3*(L3+ONE) - L1*(L1-ONE)) 340 3 * (L5*(L5+ONE) + L6*(L6+ONE) - L1*(L1-ONE)) 341C 342 DENOM = L1 * NEWFAC 343 C1 = - (L1+L1-ONE) * DV / DENOM 344 IF(LSTEP.GT.2) GO TO 120 345C 346C If L1 = L1MAX + 1 the third term in the recursion equation vanishes 347C 348 Y = SRTINY * C1 349 SIXCOF(NFIN-1) = Y 350 IF(LSTEP.EQ.NSTEP2) GO TO 200 351 SUMBAC = SUM2 352 SUM2 = SUM2 + (L1+L1-THREE) * C1 * C1 * TINY 353 GO TO 110 354C 355C 356 120 C2 = - (L1-ONE) * OLDFAC / DENOM 357C 358C Recursion to the next 6j coefficient Y 359C 360 Y = C1 * SIXCOF(NFINP2-LSTEP) + C2 * SIXCOF(NFINP3-LSTEP) 361 IF(LSTEP.EQ.NSTEP2) GO TO 200 362 SIXCOF(NFINP1-LSTEP) = Y 363 SUMBAC = SUM2 364 SUM2 = SUM2 + (L1+L1-THREE) * Y * Y 365C 366C See if last unnormalized 6j coefficient exceeds SRHUGE 367C 368 IF(ABS(Y).LT.SRHUGE) GO TO 110 369C 370C This is reached if last 6j coefficient larger than SRHUGE, 371C so that the recursion series SIXCOF(NFIN), ... ,SIXCOF(NFIN-LSTEP+1) 372C has to be rescaled to prevent overflow 373C 374C LSCALE = LSCALE + 1 375 DO 130 I=1,LSTEP 376 INDEX = NFIN-I+1 377 IF(ABS(SIXCOF(INDEX)).LT.SRTINY) SIXCOF(INDEX) = ZERO 378 130 SIXCOF(INDEX) = SIXCOF(INDEX) / SRHUGE 379 SUMBAC = SUMBAC / HUGE 380 SUM2 = SUM2 / HUGE 381C 382 GO TO 110 383C 384C 385C The forward recursion 6j coefficients X1, X2, X3 are to be matched 386C with the corresponding backward recursion values Y1, Y2, Y3. 387C 388 200 Y3 = Y 389 Y2 = SIXCOF(NFINP2-LSTEP) 390 Y1 = SIXCOF(NFINP3-LSTEP) 391C 392C 393C Determine now RATIO such that YI = RATIO * XI (I=1,2,3) holds 394C with minimal error. 395C 396 RATIO = ( X1*Y1 + X2*Y2 + X3*Y3 ) / ( X1*X1 + X2*X2 + X3*X3 ) 397 NLIM = NFIN - NSTEP2 + 1 398C 399 IF(ABS(RATIO).LT.ONE) GO TO 211 400C 401 DO 210 N=1,NLIM 402 210 SIXCOF(N) = RATIO * SIXCOF(N) 403 SUMUNI = RATIO * RATIO * SUMFOR + SUMBAC 404 GO TO 230 405C 406 211 NLIM = NLIM + 1 407 RATIO = ONE / RATIO 408 DO 212 N=NLIM,NFIN 409 212 SIXCOF(N) = RATIO * SIXCOF(N) 410 SUMUNI = SUMFOR + RATIO*RATIO*SUMBAC 411 GO TO 230 412C 413 220 SUMUNI = SUM1 414C 415C 416C Normalize 6j coefficients 417C 418 230 CNORM = ONE / SQRT((L4+L4+ONE)*SUMUNI) 419C 420C Sign convention for last 6j coefficient determines overall phase 421C 422 SIGN1 = SIGN(ONE,SIXCOF(NFIN)) 423 SIGN2 = (-ONE) ** INT(L2+L3+L5+L6+EPS) 424 IF(SIGN1*SIGN2) 235,235,236 425 235 CNORM = - CNORM 426C 427 236 IF(ABS(CNORM).LT.ONE) GO TO 250 428C 429 DO 240 N=1,NFIN 430 240 SIXCOF(N) = CNORM * SIXCOF(N) 431 RETURN 432C 433 250 THRESH = TINY / ABS(CNORM) 434 DO 251 N=1,NFIN 435 IF(ABS(SIXCOF(N)).LT.THRESH) SIXCOF(N) = ZERO 436 251 SIXCOF(N) = CNORM * SIXCOF(N) 437C 438 RETURN 439 END 440