*DECK DQC36J SUBROUTINE DQC36J (LUN, KPRINT, IPASS) C***BEGIN PROLOGUE DQC36J C***SUBSIDIARY C***PURPOSE THIS IS A QUICK CHECK PROGRAM FOR THE SUBROUTINES DRC3JJ, C DRC3JM, AND DRC6J, WHICH CALCULATE THE WIGNER COEFFICIENTS, C 3J AND 6J. C***LIBRARY SLATEC C***CATEGORY C19 C***TYPE DOUBLE PRECISION (QC36J-S, DQC36J-D) C***KEYWORDS 3J COEFFICIENTS, 3J SYMBOLS, 6J COEFFICIENTS, 6J SYMBOLS, C CLEBSCH-GORDAN COEFFICIENTS, QUICK CHECK, C RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS, C WIGNER COEFFICIENTS C***AUTHOR LOZIER, DANIEL W., (NIST) C MCCLAIN, MARJORIE A., (NIST) C SMITH, JOHN M., (NIST AND GEORGE MASON UNIVERSITY) C***REFERENCES MESSIAH, ALBERT., QUANTUM MECHANICS, VOLUME II, C NORTH-HOLLAND PUBLISHING COMPANY, 1963. C***ROUTINES CALLED D1MACH, DRC3JJ, DRC3JM, DRC6J, NUMXER, XERCLR, C XSETF C***REVISION HISTORY (YYMMDD) C 891129 DATE WRITTEN C 910415 Mixed type expressions eliminated; precision of output C formats made uniform for all tests; detail added to output C when KPRINT=2 and a test fails; name of quick check added C to output when KPRINT=3 or KPRINT=2 and a test fails; some C output formats modified for clarity or adherence to SLATEC C guidelines. These changes were done by D. W. Lozier. C 930115 Replaced direct calculation of 3j-6j symbols in tests 1, 2, C and 4 with values stored in data statements. This involved C removing all calls to subroutine DRACAH. These changes were C made by M. McClain. C***END PROLOGUE DQC36J C INTEGER LUN, KPRINT, IPASS C CHARACTER STRING*36, FMT*30, FMT2*13 INTEGER IPASS1, IPASS2, IPASS3, IPASS4, IPASS5, NDIM, IER, INDEX, + I, FIRST, LAST, NSIG, NUMXER, NERR, IERJJ, IERJM PARAMETER(NDIM=15) DOUBLE PRECISION TOL, L1, L2, L3, M1, M2, M3, L1MIN, L1MAX, M2MIN, + M2MAX, DIFF(NDIM), D1MACH, X, JJVAL, JMVAL, THRCOF(NDIM), + SIXCOF(NDIM), R3JJ(8), R3JM(14), R6J(15) C DATA R3JJ / 2.7888667551135851599272400859506249646427D-1, + -9.5346258924559231544677592152721599861388D-2, + -6.7419986246324208624649067643642846008908D-2, + 1.5331103516796664129653641995122585402552D-1, + -1.5644655469368596972508355755184909201031D-1, + 1.0994504121565505107947893271429777505797D-1, + -5.5362356931317194333395729256559987745156D-2, + 1.7998354511377858329814092962590761537262D-2/ C DATA R3JM / 2.0915897328861524261384476677886072045904D-2, + 8.5375655532152472212727551895879778672762D-2, + 9.0829537086869251694343772675676175068677D-2, + -3.8905437784649939169989036459327065796765D-2, + -6.6373497016568063569146501153397525003444D-2, + 6.4952404052838939503061387831391216401903D-2, + 2.1589431059540375939250708046202926313913D-2, + -7.7891271178523921999229618972588887261359D-2, + 3.5976437105954340188005810512211794384411D-2, + 5.4730150002126342307937096038252488407360D-2, + -7.5967866595676151462927617736745078548338D-2, + -2.1922444553989211377558215380002910257762D-2, + 1.0116774428077220242411199686231560525497D-1, + 7.3482572624471970469595137204530687381176D-2/ C DATA R6J / 3.4909051383732997774596981092927782159095D-2, + -3.7430250396597916085929064401358002747549D-2, + 1.8908663909595601841537964135129184202064D-2, + 7.3424482549286434570947151839589351100581D-3, + -2.3589351850817944585847816357296508528608D-2, + 1.9134769552154365200026782557432864615918D-2, + 1.2880173977241722084434864685591278730958D-3, + -1.9300183662905265397749119277519305417805D-2, + 1.6773059493828887697413611251392749162229D-2, + 5.5011472748509487167380502058890639729979D-3, + -2.1354397908968309742136976853078409839580D-2, + 3.4603644514353873082775312319159137043869D-3, + 2.5209500547955845860442730268272167527589D-2, + 1.4839905612217133028540464232557124565509D-2, + 2.7085776806331855972407001825016114677027D-3/ C C***FIRST EXECUTABLE STATEMENT DQC36J C C --- INITIALIZATION OF TESTS TOL=100.0D0*D1MACH(3) IF(KPRINT.GE.2)THEN WRITE(LUN,*)' THIS IS DQC36J, A TEST PROGRAM FOR THE ' // + 'DOUBLE PRECISION 3J6J PACKAGE.' WRITE(LUN,*)' AN EXPLANATION OF THE VARIOUS ' // + 'TESTS CAN BE FOUND IN THE PROGRAM COMMENTS.' WRITE(LUN,*) ENDIF C C --- FIND NUMBER OF SIGNIFICANT FIGURES FOR FORMATTING X=1.D0/3.D0 WRITE(STRING,100)X 100 FORMAT(F35.25) DO 200 I=1,35 IF(STRING(I:I).EQ.'3')THEN FIRST=I GOTO 300 ENDIF 200 CONTINUE 300 CONTINUE DO 400 I=FIRST,35 IF(STRING(I:I).NE.'3')THEN LAST=I-1 GOTO 500 ENDIF 400 CONTINUE LAST=36 500 CONTINUE NSIG=LAST-FIRST+1 FMT(1:16)='(1X,F5.1,T8,G35.' WRITE(FMT(17:18),'(I2)')NSIG FMT(19:27)=',T45,G35.' WRITE(FMT(28:29),'(I2)')NSIG FMT(30:30)=')' FMT2(1:10)='(1X,A,G35.' WRITE(FMT2(11:12),'(I2)')NSIG FMT2(13:13)=')' C C --- TEST 1: COMPARE DRC3JJ VALUES WITH FORMULA IPASS1=1 L2=4.5D0 L3=3.5D0 M2=-3.5D0 M3=2.5D0 CALL DRC3JJ(L2,L3,M2,M3,L1MIN,L1MAX,THRCOF,NDIM,IER) IF(IER.NE.0)THEN IPASS1=0 ELSE DO 550 L1=L1MIN,L1MAX INDEX=INT(L1-L1MIN)+1 M1=1.0D0 DIFF(INDEX)=ABS(THRCOF(INDEX)-R3JJ(INDEX)) IF(DIFF(INDEX).GT.ABS(R3JJ(INDEX))*TOL)IPASS1=0 550 CONTINUE ENDIF IF(KPRINT.GE.3 .OR. (KPRINT.EQ.2.AND.IPASS1.EQ.0))THEN WRITE(LUN,*)' TEST 1, RECURRENCE IN L1, COMPARE VALUES OF 3J ', + 'CALCULATED BY DRC3JJ TO' WRITE(LUN,*)' VALUES CALCULATED BY EXPLICIT FORMULA FROM ', + 'MESSIAH''S QUANTUM MECHANICS' WRITE(LUN,600)L2,L3 600 FORMAT(' L2 = ',F5.1,' L3 = ',F5.1) WRITE(LUN,700)M1,M2,M3 700 FORMAT(' M1 = ',F5.1,' M2 = ',F5.1,' M3 = ',F5.1) IF(IER.NE.0)THEN WRITE(LUN,*)' ERROR RETURNED FROM SUBROUTINE ', + 'DRC3JJ: IER =',IER ELSE WRITE(LUN,800) 800 FORMAT(' L1',T31,'DRC3JJ VALUE',T67,'FORMULA VALUE') DO 900 L1=L1MIN,L1MAX INDEX=INT(L1-L1MIN)+1 WRITE(LUN,FMT)L1,THRCOF(INDEX),R3JJ(INDEX) IF(DIFF(INDEX).GT.ABS(R3JJ(INDEX))*TOL)THEN WRITE(LUN,'(1X,A,F5.1)')'DIFFERENCE EXCEEDS ERROR '// + 'TOLERANCE FOR L1 =',L1 ENDIF 900 CONTINUE ENDIF ENDIF IF(IPASS1.EQ.0)THEN IF(KPRINT.GE.1)THEN WRITE(LUN,*)' ***** ***** TEST 1 FAILED ***** *****' WRITE(LUN,*) ENDIF ELSE IF(KPRINT.GE.2)THEN WRITE(LUN,*)' ***** ***** TEST 1 PASSED ***** *****' WRITE(LUN,*) ENDIF ENDIF C C --- TEST 2: COMPARE DRC3JM VALUES WITH FORMULA IPASS2=1 L1=8.0D0 L2=7.5D0 L3=6.5D0 M1=1.0D0 CALL DRC3JM(L1,L2,L3,M1,M2MIN,M2MAX,THRCOF,NDIM,IER) IF(IER.NE.0)THEN IPASS2=0 ELSE DO 950 M2=M2MIN,M2MAX INDEX=INT(M2-M2MIN)+1 M3=-M1-M2 DIFF(INDEX)=ABS(THRCOF(INDEX)-R3JM(INDEX)) IF(DIFF(INDEX).GT.ABS(R3JM(INDEX))*TOL)IPASS2=0 950 CONTINUE ENDIF IF(KPRINT.GE.3 .OR. (KPRINT.EQ.2.AND.IPASS2.EQ.0))THEN WRITE(LUN,*)' TEST 2, RECURRENCE IN M2, COMPARE VALUES OF 3J ', + 'CALCULATED BY DRC3JM TO' WRITE(LUN,*)' VALUES CALCULATED BY EXPLICIT FORMULA FROM ', + 'MESSIAH''S QUANTUM MECHANICS' WRITE(LUN,1000)L1,L2,L3 1000 FORMAT(' L1 = ',F5.1,' L2 = ',F5.1,' L3 = ',F5.1) WRITE(LUN,1100)M1 1100 FORMAT(' M1 = ',F5.1,' M3 = -(M1+M2)') IF(IER.NE.0)THEN WRITE(LUN,*)' ERROR RETURNED FROM SUBROUTINE ', + 'DRC3JM: IER =',IER ELSE WRITE(LUN,1200) 1200 FORMAT(' M2',T31,'DRC3JM VALUE',T67,'FORMULA VALUE') DO 1300 M2=M2MIN,M2MAX INDEX=INT(M2-M2MIN)+1 WRITE(LUN,FMT)M2,THRCOF(INDEX),R3JM(INDEX) IF(DIFF(INDEX).GT.ABS(R3JM(INDEX))*TOL)THEN WRITE(LUN,'(1X,A,F5.1)')'DIFFERENCE EXCEEDS ERROR '// + 'TOLERANCE FOR M2 =',M2 ENDIF 1300 CONTINUE ENDIF ENDIF IF(IPASS2.EQ.0)THEN IF(KPRINT.GE.1)THEN WRITE(LUN,*)' ***** ***** TEST 2 FAILED ***** *****' WRITE(LUN,*) ENDIF ELSE IF(KPRINT.GE.2)THEN WRITE(LUN,*)' ***** ***** TEST 2 PASSED ***** *****' WRITE(LUN,*) ENDIF ENDIF C C --- TEST3: COMPARE COMMON VALUE OF DRC3JJ AND DRC3JM IPASS3=1 L1=100.0D0 L2=2.0D0 L3=100.0D0 M1=-10.0D0 M2=0.0D0 M3=10.0D0 CALL DRC3JJ(L2,L3,M2,M3,L1MIN,L1MAX,THRCOF,NDIM,IERJJ) JJVAL=THRCOF(3) CALL DRC3JM(L1,L2,L3,M1,M2MIN,M2MAX,THRCOF,NDIM,IERJM) JMVAL=THRCOF(3) IF(IERJJ.NE.0 .OR. IERJM.NE.0)THEN IPASS3=0 ELSE DIFF(1)=ABS(JJVAL-JMVAL) IF(DIFF(1).GT.0.5*ABS(JJVAL+JMVAL)*TOL)IPASS3=0 ENDIF IF(KPRINT.GE.3 .OR. (KPRINT.EQ.2.AND.IPASS3.EQ.0))THEN WRITE(LUN,*)' TEST 3, COMPARE A COMMON VALUE CALCULATED BY ', + 'BOTH DRC3JJ AND DRC3JM' WRITE(LUN,*)' L1 = 100.0 L2 = 2.0 L3 = 100.0' WRITE(LUN,*)' M1 = -10.0 M2 = 0.0 M3 = 10.0' IF(IERJJ.NE.0)THEN WRITE(LUN,*)' ERROR RETURNED FROM SUBROUTINE ', + 'DRC3JJ: IER =',IERJJ ELSEIF(IERJM.NE.0)THEN WRITE(LUN,*)' ERROR RETURNED FROM SUBROUTINE ', + 'DRC3JM: IER =',IERJM ELSE WRITE(LUN,FMT2)'DRC3JJ VALUE =',JJVAL WRITE(LUN,FMT2)'DRC3JM VALUE =',JMVAL IF(DIFF(1).GT.0.5*ABS(JJVAL+JMVAL)*TOL)THEN WRITE(LUN,'(1X,A)')'DIFFERENCE EXCEEDS ERROR TOLERANCE' ENDIF ENDIF ENDIF IF(IPASS3.EQ.0)THEN IF(KPRINT.GE.1)THEN WRITE(LUN,*)' ***** ***** TEST 3 FAILED ***** *****' WRITE(LUN,*) ENDIF ELSE IF(KPRINT.GE.2)THEN WRITE(LUN,*)' ***** ***** TEST 3 PASSED ***** *****' WRITE(LUN,*) ENDIF ENDIF C C --- TEST 4: COMPARE DRC6J VALUES WITH FORMULA IPASS4=1 L2=8.0D0 L3=7.0D0 M1=6.5D0 M2=7.5D0 M3=7.5D0 CALL DRC6J(L2,L3,M1,M2,M3,L1MIN,L1MAX,SIXCOF,NDIM,IER) IF(IER.NE.0)THEN IPASS4=0 ELSE DO 1310 L1=L1MIN,L1MAX INDEX=INT(L1-L1MIN)+1 DIFF(INDEX)=ABS(SIXCOF(INDEX)-R6J(INDEX)) IF(DIFF(INDEX).GT.ABS(R6J(INDEX))*TOL)IPASS4=0 1310 CONTINUE ENDIF IF(KPRINT.GE.3 .OR. (KPRINT.EQ.2.AND.IPASS4.EQ.0))THEN WRITE(LUN,*)' TEST 4, RECURRENCE IN L1, COMPARE VALUES OF 6J ', + 'CALCULATED BY DRC6J TO' WRITE(LUN,*)' VALUES CALCULATED BY EXPLICIT FORMULA FROM ', + 'MESSIAH''S QUANTUM MECHANICS' WRITE(LUN,600)L2,L3 WRITE(LUN,700)M1,M2,M3 IF(IER.NE.0)THEN WRITE(LUN,*)' ERROR RETURNED FROM SUBROUTINE ', + 'DRC6J: IER =',IER ELSE WRITE(LUN,1320) 1320 FORMAT(' L1',T32,'DRC6J VALUE',T67,'FORMULA VALUE') DO 1350 L1=L1MIN,L1MAX INDEX=INT(L1-L1MIN)+1 WRITE(LUN,FMT)L1,SIXCOF(INDEX),R6J(INDEX) IF(DIFF(INDEX).GT.ABS(R6J(INDEX))*TOL)THEN WRITE(LUN,'(1X,A,F5.1)')'DIFFERENCE EXCEEDS ERROR '// + 'TOLERANCE FOR L1 =',L1 ENDIF 1350 CONTINUE ENDIF ENDIF IF(IPASS4.EQ.0)THEN IF(KPRINT.GE.1)THEN WRITE(LUN,*)' ***** ***** TEST 4 FAILED ***** *****' WRITE(LUN,*) ENDIF ELSE IF(KPRINT.GE.2)THEN WRITE(LUN,*)' ***** ***** TEST 4 PASSED ***** *****' WRITE(LUN,*) ENDIF ENDIF C C --- TEST 5: CHECK INVALID INPUT IPASS5=1 IF(KPRINT.LE.2)THEN CALL XSETF(0) ELSE CALL XSETF(-1) ENDIF IF(KPRINT.GE.3)WRITE(LUN,*)' TEST 5, CHECK FOR PROPER HANDLING ', + 'OF INVALID INPUT' C --- DRC3JJ: L2-ABS(M2) OR L3-ABS(M3) LESS THAN ZERO (IER=1) L2=2.0D0 L3=100.0D0 M1=-6.0D0 M2=-4.0D0 M3=10.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JJ(L2,L3,M2,M3,L1MIN,L1MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC3JJ: L2+ABS(M2) OR L3+ABS(M3) NOT INTEGER (IER=2) L2=2.0D0 L3=99.5D0 M1=-10.0D0 M2=0.0D0 M3=10.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JJ(L2,L3,M2,M3,L1MIN,L1MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC3JJ: L1MAX-L1MIN NOT INTEGER (IER=3) L2=3.2D0 L3=4.5D0 M1=-1.3D0 M2=0.8D0 M3=0.5D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JJ(L2,L3,M2,M3,L1MIN,L1MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC3JJ: L1MIN GREATER THAN L1MAX (IER=4) C (NO TEST -- THIS ERROR SHOULD NEVER OCCUR) C --- DRC3JJ: DIMENSION OF THRCOF TOO SMALL (IER=5) L2=10.0D0 L3=150.0D0 M1=-10.0D0 M2=0.0D0 M3=10.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JJ(L2,L3,M2,M3,L1MIN,L1MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC3JM: L1-ABS(M1) .LT. ZERO OR L1+ABS(M1) NOT INTEGER (IER=1) L1=100.0D0 L2=2.0D0 L3=100.0D0 M1=150.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JM(L1,L2,L3,M1,M2MIN,M2MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC3JM: L1, L2, L3 DO NOT SATISFY TRIANGULAR CONDITION (IER=2) L1=20.0D0 L2=5.0D0 L3=10.0D0 M1=-10.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JM(L1,L2,L3,M1,M2MIN,M2MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC3JM: L1+L2+L3 NOT INTEGER (IER=3) L1=1.0D0 L2=1.3D0 L3=1.5D0 M1=0.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JM(L1,L2,L3,M1,M2MIN,M2MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC3JM: M2MAX-M2MIN NOT INTEGER (IER=4) L1=1.0D0 L2=1.3D0 L3=1.7D0 M1=0.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JM(L1,L2,L3,M1,M2MIN,M2MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC3JM: M2MIN GREATER THAN M2MAX (IER=5) C (NO TEST -- THIS ERROR SHOULD NEVER OCCUR) C --- DRC3JM: DIMENSION OF THRCOF TOO SMALL (IER=6) L1=100.0D0 L2=10.0D0 L3=110.0D0 M1=-10.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC3JM(L1,L2,L3,M1,M2MIN,M2MAX,THRCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC6J: L2+L3+L5+L6 OR L4+L2+L6 NOT INTEGER (IER=1) L2=0.5D0 L3=1.0D0 M1=0.5D0 M2=2.0D0 M3=3.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC6J(L2,L3,M1,M2,M3,L1MIN,L1MAX,SIXCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC6J: L4, L2, L6 TRIANGULAR CONDITION NOT SATISFIED (IER=2) L2=1.0D0 L3=3.0D0 M1=5.0D0 M2=6.0D0 M3=2.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC6J(L2,L3,M1,M2,M3,L1MIN,L1MAX,SIXCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC6J: L4, L5, L3 TRIANGULAR CONDITION NOT SATISFIED (IER=3) L2=4.0D0 L3=1.0D0 M1=5.0D0 M2=3.0D0 M3=2.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC6J(L2,L3,M1,M2,M3,L1MIN,L1MAX,SIXCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC6J: L1MAX-L1MIN NOT INTEGER (IER=4) L2=0.9D0 L3=0.5D0 M1=0.9D0 M2=0.4D0 M3=0.2D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC6J(L2,L3,M1,M2,M3,L1MIN,L1MAX,SIXCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 C --- DRC6J: L1MIN GREATER THAN L1MAX (IER=5) C (NO TEST -- THIS ERROR SHOULD NEVER OCCUR) C --- DRC6J: DIMENSION OF SIXCOF TOO SMALL (IER=6) L2=50.0D0 L3=25.0D0 M1=15.0D0 M2=30.0D0 M3=40.0D0 IF(KPRINT.GE.3)WRITE(LUN,*) CALL XERCLR CALL DRC6J(L2,L3,M1,M2,M3,L1MIN,L1MAX,SIXCOF,NDIM,IER) IF(NUMXER(NERR).NE.IER)IPASS5=0 IF(IPASS5.EQ.0)THEN IF(KPRINT.GE.1)THEN WRITE(LUN,*)' ***** ***** TEST 5 FAILED ***** *****' WRITE(LUN,*) ENDIF ELSE IF(KPRINT.GE.2)THEN WRITE(LUN,*)' ***** ***** TEST 5 PASSED ***** *****' WRITE(LUN,*) ENDIF ENDIF C C --- END OF TESTS IF((IPASS1.EQ.0).OR.(IPASS2.EQ.0).OR.(IPASS3.EQ.0).OR. + (IPASS4.EQ.0).OR.(IPASS5.EQ.0))THEN IPASS=0 IF(KPRINT.GE.1)WRITE(LUN,1500) ELSE IPASS=1 IF(KPRINT.GE.2)WRITE(LUN,1600) ENDIF 1500 FORMAT(' ***** DQC36J FAILED SOME TESTS *****') 1600 FORMAT(' ***** DQC36J PASSED ALL TESTS *****') C RETURN END