drc3jm.f 13 KB

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