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 |