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