C LIBRARY ROUTINES FOR COMBINED STG2 OR RECUPD OR STGH - RL'97May14 C '05May31: ORDTRI corrected as by SMWB'98Apr29 '05JUN06 C*********************************************************************** C C SUBROUTINES INCLUDED IN LIBRARY: C C ORTHOG C FANO C GENSUM, CHVAR C RME REAL FUNCTION C NTAB1 INTEGER FUNCTION C MUMDAD C CFP C CFPF C SETJ1 C J23SPN C MODJ23 C J23ANG C BLOCK DATA - CONTAINS /BPSIZE/ AND NOW CFPD BLOCK WITH /FRPAR2/ C REDUCE C MEKEEP C MEREST C H0WTS C SETM, SETDM C CFPP C CFPD C TENSOR C SETUPE C NJSYM rather NJGRAF C BUBBLE, LOLPOP, ORDTRI, SQUARE, SPRATE,VAR C CUT1L, CUT2L, CUTNL, SEARCH C SETTAB C DIAGRM, INTAB, WAY C ZERO C CHANGE C DELTA C FACTT C DRACAH C CG REAL FUNCTION C INTACT C CHOP C BAKSUB C EIGEN C EIGVEC C HOUSE C HSLDR C NORM C VECTOR C C COMMON BLOCKS USED IN LIBRARY: C C COMMON/BPSIZE/KFLN,KFL2,KFLM,KDUMMY(6) C " " " DELIMITING /MDEFN/ AND /MSTATE/! C COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), C + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP C COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS C COMMON/COUPLE/NJ1S,NJ23S,J1(&L76),J2(&L74,3),J3(&L74,3),FREE(&L76) C COMMON/CUTDIG/CUT C COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 C COMMON/DIAGNL/IDIAG,JA,JB C COMMON/ENAV/COEFCT(5),NINTS,KVALUE(5) C COMMON/FACT/ GAMMA( 250) C COMMON/FACTS/ GAM(199) C COMMON/FRPAR2/K(5),IV(5,16),IL(5,16),IS(5,16),ITAB1(5,1),ITAB2(8,5 C 1 ),ITAB3(16,8),ITAB4(16,16),NORM1(5),NORM2(8),NORM3(16),NORM4(16) C COMMON/HOLD/J2SPIN(&L76),J3SPIN(&L76),J2ANG(&L76),J3ANG(&L76) C COMMON/INFORM/IREAD,IWRITE,IPUNCH C COMMON/INTERM/J1BAR1(&L75,3),J1BAR2(&L75,3),J1TLD1(3),J1TLD2(3) C COMMON/KRON/ IDEL(&L75,&L75) C COMMON/MEDEFN/IHSH,NJ(&L75),LJ(&L75),NOSH1(&L75),NOSH2(&L75), C 1 J1QN1(&L43,3),J1QN2(&L43,3),IJFUL(&L75) C COMMON/MSTATE/MCFG,MOCCSH(2500),MOCORB(&L61,2500), C 1 MELCSH(&L61,2500),M1QNRD(&L42,3,2500),KCFG, C 2 KOCCSH(2500),KOCORB(&L61,2500),KELCSH(&L61,2500), C 3 K1QNRD(&L42,3,2500),MAXOR C COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15, C 1 M16,M17,M18,M19,M20 C COMMON/NAM/NAMSUB C COMMON/NJLJ/NRHO,LRHO,NSIG,LSIG,NRHOP,LRHOP,NSIGP,LSIGP C COMMON/STORE/I(&L99),I1,I2,I3,I4 C COMMON/TERMS/ NROWS,ITAB(18),JTAB(18),NTAB(189) C COMMON/XATION/ AMULT(21),BMULT(21),KD1,KD2,KE1,KE2,MULTD,MULTE C (99 99) --- TEMPFIX RUB'95JUL9! C COMMON/ZER/NZERO,JZERO(KFLZ) C*********************************************************************** C C INSERT THE FOLLOWING SEGMENTS OF THE GENERAL PROGRAM TO CALCULATE C ANGULAR MOMENTUM INTEGRALS IN ATOMIC STRUCTURE BY A HIBBERT. C CATALOGUE NUMBER ACQV (CPC 2 (1971) 180) MODIFIED BY THE C CORRECTION DECK CATALOGUE NUMBER ACQV000A (CPC 6 (1973) 59) AND C THE ADAPTATION DECK CATALOGUE NUMBER ACQV0001 (CPC 7 (1974) 318) C C ORTHOG,RKWTS,FANO,RME,NTAB1,MUMDAD,CFP,CFPF,SETJ1,J23SPN,MODJ23, C J23ANG,BLOCK DATA,INTACT,CHOP,REDUCE,MEKEEP,MEREST,H0WTS,SETM C C CARDS ACQV0573-ACQV2499,ACQV001A-ACQV098A,ACQVA006-ACQVA856 C C*********************************************************************** C C C*********************************************************************** C C +++ REPLACE CARD ACQV0861 BY THE NEXT CARD WITHOUT THE C C DIMENSION K6(100),K7(200),K8(100),KW(6,20) C +++ REPLACE CARDS ACQV873,ACQV1499,ACQV1553,ACQV1097,ACQV1751,AND C +++ ACQVA143 BY A COPY OF THE NEXT CARD WITHOUT THE C C COMMON/COUPLE/NJ1S,NJ23S,J1(100),J2(&L74,3),J3(&L74,3),FREE(100) C +++ REPLACE CARD ACQVA141 BY THE NEXT CARD WITHOUT THE C C DIMENSION L6(100),L7(200),L8(100),LW(6,20) C +++ A DIMENSION CHANGE WAS NECESSARY C C +++ REPLACE CARD ACQV070A BY A COPY OF THE NEXT CARD WITHOUT THE C C 2 7.0E 00,1.1E 01,1.0E-05/ C +++ A DIFFERENT ACCURACY PARAMETER WAS USED C C*********************************************************************** C SUBROUTINE ORTHOG(LET) C C CHECKS FOR POSSIBLE ORTHOGONALITY DUE TO C COUPLING DIFFERENCES OR UNEVEN PARITY C PARAMETER (LL43= 21*2+3, LL75= 21+2) COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75), CC REPLACE CARDS 548,577,673,867,1494,2262,2309,2355 BY A COPY OF C 1 ... J1QN2(45,3),IJFUL(23) 1 J1QN1(LL43,3),J1QN2(LL43,3),IJFUL(LL75) C 101 FORMAT(37H DIFFERING RESULTANT ANGULAR MOMENTUM) 102 FORMAT(52H ORTHOGONALITY IN COUPLING SCHEMES OF CONFIGURATIONS) 103 FORMAT(59H THE TWO CONFIGURATIONS HAVE DIFFERING NUMBERS OF ELECTR *ONS) 104 FORMAT(51H THE TWO CONFIGURATIONS HAVE DIFFERING TOTAL PARITY) C C --- DO PSI AND PSIP CONTAIN THE SAME NUMBERS OF ELECTRONS C DO PSI AND PSIP HAVE THE SAME TOTAL PARITY C N5=0 N6=0 N7=0 DO 20 I=1,IHSH L1=LJ(I) L2=NOSH1(I) L3=NOSH2(I) N5=N5+L2 N6=N6+L3 20 N7=N7+L1*(L2-L3) C C CHECK ON NUMBER OF ELECTRONS C IF (N5.EQ.N6) GO TO 22 IF(IBUG1.GE.1) WRITE(IWRITE,103) GO TO 11 C C CHECK ON PARITY C 22 IF(MOD(N7,2).EQ.0) GO TO 24 IF(IBUG1.GE.1) WRITE(IWRITE,104) GO TO 11 24 N1=2*IHSH-1 N2=IHSH+1 N3=IHSH-1 N4=IHSH-2 C C --- IS THE FINAL STATE THE SAME FOR PSI AND PSIP C DO 1 K=2,3 IF(J1QN1(N1,K).NE.J1QN2(N1,K)) GO TO 2 1 CONTINUE GO TO 12 2 IF(IBUG1.GE.1) WRITE(IWRITE,101) C C --- THE TWO CONFIGURATIONS WILL HAVE ZERO HAMILTONIAN MATRIX ELEMENT C 11 LET=0 RETURN C REPLACE CARDS 0634 THROUGH 0662 BY THE FOLLOWING CARD C C --- NO OBVIOUS ANGULAR MOMENTUM ORTHOGONALITY C 12 LET=1 RETURN END C*********************************************************************** REAL FUNCTION RME(L,LP,K) C C --- EVALUATES THE REDUCED MATRIX ELEMENT (L//C(K)//LP) - SEE FANO C AND RACAH, IRREDUCIBLE TENSORIAL SETS, CHAP. 14, P. 81 C C IMPLICIT REAL*8(A-H,O-Z) C INSERT A COPY OF THE NEXT CARD AFTER CARDS 314,861,1285,1389, C 1481,1837,1883 AND 2257 COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 COMMON/FACT/ GAMMA( 250) C REPLACE CARDS 1293 AND 1294 BY THE FOLLOWING THREE CARDS 200 FORMAT(//4H L =,I3,5H LP =,I3,4H K =,I3,44H: RME SET ZERO SINCE AN * ANGLE DOES NOT MATCH) IF(K.GT.L+LP.OR.K.LT.ABS(L-LP)) GO TO 100 I2G=L+LP+K IG=I2G/2 IF(I2G.EQ.2*IG) GO TO 2 C REPLACE CARD NUMBER 1298 BY THE FOLLOWING CARD 1 RME=0.0 RETURN 100 IF(IBUG1.GT.1) WRITE(IWRITE,200) L,LP,K GO TO 1 2 IF(IG) 100,13,12 C REPLACE CARD NUMBER 1304 BY THE FOLLOWING CARD 13 RME=1.0 RETURN 12 I1=IG-L I2=IG-LP I3=IG-K RME=SQRT((2*L+1)*(2*LP+1)*GAMMA(2*I1+1)*GAMMA(2*I2+1)* * GAMMA(2*I3+1)/GAMMA(I2G+2)) * * GAMMA(IG+1)/(GAMMA(I1+1)*GAMMA(I2+1)*GAMMA(I3+1)) RETURN END INTEGER FUNCTION NTAB1(NELCTS,K) C C - CALCULATES THE ROW OF NTAB CORRESPONDING TO THE PARENTS C WHICH MAY GIVE RISE TO THE TERM ASSOCIATED WITH SHELL LAMBDA. C E.G. IF WE SEEK THE ROW OF NTAB CONTAINING THE PARENTS C OF ONE OF THE P**3 TERMS, THE ROW = VALUE OF NTAB1 IS THAT C CONTAINING THE P**2 TERMS C C USE IS MADE OF THE FACT THAT THE LIST OF POSSIBLE PARENTS (SEE C WHITE - ATOMIC SPECTRA - APPENDIX) IS SYMMETRICAL ABOUT THE C CONFIGURATION L**(2L+1) C C C --- FOR ONE ELECTRON IN A TERM, THE PARENT IS ALWAYS A SINGLET S TERM C COMMON/INFORM/IREAD,IWRITE,IPUNCH IF(NELCTS-1) 1,2,1 2 NTAB1=2 RETURN C C OTHERWISE THE VALUE OF NTAB1 DEPENDS ON THE L VALUE OF THE C ELECTRONS C 1 GO TO (3,4,5,6,14,25,26,27,28,29,30),K C C --- FOR S ELECTRONS, THE ONLY OTHER POSSIBILITY IS THAT NELCTS=2 C 3 NTAB1=1 RETURN C C --- P ELECTRONS - ARE WE BEYOND P**4 C 4 IF(NELCTS.LE.4) GO TO 7 8 NELCTS=8-NELCTS 7 NTAB1=1+NELCTS RETURN C C --- D ELECTRONS - ARE WE BEYOND D**6 C 5 IF(NELCTS.LE.6) GO TO 9 NELCTS=12-NELCTS 9 NTAB1=4+NELCTS RETURN C C --- F ELECTRONS - ARE THERE MORE THAN TWO. IF SO, THE PROGRAMME NEEDS C AN F-SHELL COEFFICIENT-OF-FRACTIONAL-PARENTAGE ROUTINE, AND THE C ARRAYS IN /TERMS/ NEED EXTENDING C 6 IF(NELCTS-2) 2,11,12 11 NTAB1 = 11 RETURN C C C L=3 ELECTRONS SHOULD NOT BE MORE THAN TWO C C 25 IF(NELCTS-2)2,35,45 35 NTAB1=13 RETURN C 26 IF(NELCTS-2)2,36,45 36 NTAB1=14 RETURN C 27 IF(NELCTS-2)2,37,45 37 NTAB1=15 RETURN C 28 IF(NELCTS-2)2,38,45 38 NTAB1=16 RETURN C 29 IF(NELCTS-2)2,39,45 39 NTAB1=17 RETURN C 30 IF(NELCTS-2)2,40,45 40 NTAB1=18 45 RETURN C C --- G ELECTRONS - ARE THERE MORE THAN TWO. IF SO, THE PROGRAMME C NEEDS A G-SHELL COEFFICIENT-OF-FRACTIONAL-PARENTAGE ROUTINE, AND C THE ARRAYS IN /TERMS/ NEED EXTENDING C 14 IF(NELCTS-2) 2,15,16 15 NTAB1=12 RETURN 12 WRITE(IWRITE,13) 13 FORMAT(////67H STOP AND EXTEND THE NTAB AND ITAB ARRAYS TO ALLOW C *ORE F-ELECTRONS/78H YOU WILL ALSO REQUIRE A COMPLETE FRACTIONAL PA 2RENTAGE ROUTINE FOR F-ELECTRONS//) GO TO 17 16 WRITE(IWRITE,18) 18 FORMAT(////67H STOP AND EXTEND THE NTAB AND ITAB ARRAYS TO ALLOW C *ORE G-ELECTRONS/78H YOU WILL ALSO REQUIRE A COMPLETE FRACTIONAL PA 2RENTAGE ROUTINE FOR G-ELECTRONS//) 17 STOP END SUBROUTINE MUMDAD(II,IJ,IK,M,X) C C INSERT A COPY OF THE NEXT CARD AFTER CARDS 314,861,1285,1389, C 1481,1837,1883 AND 2257 C IMPLICIT REAL*8(A-H,O-Z) PARAMETER (LL43= 21*2+3, LL75= 21+2) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS C REPLACE CARDS 384 AND 1390 BY A COPY OF THE NEXT CARD COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH(LL75,2), * J1QN(LL43,3,2),IJFUL(LL75) COMMON/INTERM/J1B(LL75,3,2),J1T(3,2) C C NOTICE THE NAMES IN THE COMMON BLOCKS. SEE SETUP FOR DESCRIPTION C C --- CALLS AND EVALUATES FRACTIONAL PARENTAGE COEFFICIENTS C C REPLACE CARD NUMBER 1397 BY THE FOLLOWING CARD X=1.0 LIJ=LJ(IJ) IF(LIJ.GT.0) GO TO 11 IF(M)4,5,4 11 N=NOSH(IJ,II) IVI=J1QN(IJ,1,II) ILI=(J1QN(IJ,2,II)-1)/2 ISI=J1QN(IJ,3,II) C C IF M=0 THERE ARE QUANTUM NUMBERS WITH TILDES TO CONSIDER C IF(M.EQ.0) GO TO 2 1 IVJ=J1B(IJ,1,II) ILJ=(J1B(IJ,2,II)-1)/2 ISJ= J1B(IJ,3,II) GO TO 3 2 IVJ=J1T(1,II) ILJ=(J1T(2,II)-1)/2 ISJ=J1T(3,II) 3 CALL CFP(LIJ,N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP) X=X*COEFP C REPLACE CARD NUMBER 1418 BY THE FOLLOWING CARD IF(ABS(X).LT.EPS) GO TO 5 4 LIJ=LJ(IK) IF(LIJ.LE.0) GO TO 5 IF(M.EQ.0) GO TO 7 N=NOSH(IK,II) IVI=J1QN(IK,1,II) ILI=(J1QN(IK,2,II)-1)/2 ISI=J1QN(IK,3,II) IVJ = J1B(IK,1,II) ILJ =(J1B(IK,2,II)-1)/2 ISJ = J1B(IK,3,II) GO TO 8 7 N=NOSH(IJ,II)-1 IVI=IVJ ILI=ILJ ISI=ISJ IVJ=J1B(IJ,1,II) ILJ=(J1B(IJ,2,II)-1)/2 ISJ = J1B(IJ,3,II) 8 CALL CFP(LIJ,N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP) X=X*COEFP 5 RETURN END SUBROUTINE CFP(LIJ,N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP) C C === CHOOSES APPROPRIATE FRACTIONAL PARENTAGE SUBROUTINE C C IMPLICIT REAL*8(A-H,O-Z) COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 10 FORMAT(8H CFP(LL=,I2,3H) =,F13.9,'; (N: SLI, SLJ) =',I4,2(I3,I2)) C C IF F-SHELL OR G-SHELL COEFFICIENT-OF-FRACTIONAL-PARENTAGE ROUTINES C ARE INCLUDED, THIS COMPUTED GO TO NEEDS MODIFYING TO ACCOUNT FOR C THIS C GO TO (1,2,3,4), LIJ+1 C C --- FALSE CALL FOR S-SHELLS C 1 COEFP=1.0 C FORMAT(69H UNNECESSARY ATTEMPT TO FORM CFP OF AN S-ELECTRON - ???? GO TO 5 C C --- P-SHELLS C 2 CALL CFPP(N,ILI,ISI,ILJ,ISJ,COEFP) GO TO 5 C C --- D-SHELLS C 3 CALL CFPD(N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP) GO TO 5 C C --- F-SHELLS, G-SHELLS ETC. WITH UP TO TWO ELECTRONS C 4 CALL CFPF(N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP) 5 IF(IBUG1.GT.1 .OR. COEFP.EQ.9.9) WRITE(IWRITE,10) LIJ,COEFP, * N, ISI,ILI, ISJ,ILJ 6 RETURN END SUBROUTINE CFPF(N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP) C C THIS IS A DUMMY SUBROUTINE TO CALCULATE CFP OF F-ELECTRONS. IT IS C VALID ONLY FOR ONE OR TWO ELECTRONS IN THE F-SHELL UNDER C CONSIDERATION. C C INSERT A COPY OF THE NEXT CARD AFTER CARDS 314,861,1285,1389, C 1481,1837,1883 AND 2257 C IMPLICIT REAL*8(A-H,O-Z) C REPLACE CARD NUMBER 1487 BY THE FOLLOWING CARD COEFP=1.0 RETURN END ***** SUBROUTINE J23SPN(IRHO,ISIG,IRHOP,ISIGP,JSNDIR) C C === SET UP THE J2 AND J3 ARRAYS FOR THE DIRECT SPIN INTEGRAL CALL C OF NJSYM C LOGICAL FREE PARAMETER (LL74= 21+4, LL76= 21*3+19) COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15, 1 M16,M17,M18,M19,M20 C DELETE CARD NUMBER 1552 - COMMON BLOCK XATION NOT NEEDED COMMON/COUPLE/NJ1S,NJ23S,J1(LL76),J2(LL74,3),J3(LL74,3),FREE(LL76) COMMON/HOLD/J2SPIN(LL76),J3SPIN(LL76),J2ANG(LL76),J3ANG(LL76) C 304 FORMAT(3H J2,18X,2HJ3,16X,'from J23SPN'/(3I5,5X,3I5)) C C IF(JSNDIR.GT.0) GO TO 1 ! ARRAYS J2 AND J3 ALREADY SET -- C C --- SET THIRD ROW OF J2 AND J3 C M = M9 IF(M1.EQ.0) M = M10 ! case ISIG=IRHO J2(3,1)=1 IF(IRHO.NE.1) GO TO 273 J2(3,1)=M 273 J2(3,2)=2 IF(ISIG.NE.2.AND.IRHO.NE.2) GO TO 284 J2(3,2)=M IF(IRHO.NE.1) GO TO 284 J2(3,2)=M10 !M+1 284 J2(3,3)=M4 ! IHSH+1 C M = M11 ! IHSH*3+1 IF(M2.EQ.0) M = M12 ! ISIGP=IRHOP J3(3,1)=1 IF(IRHOP.NE.1) GO TO 274 J3(3,1)=M 274 J3(3,2)=2 IF(ISIGP.NE.2.AND.IRHOP.NE.2) GO TO 285 J3(3,2)=M IF(IRHOP.NE.1) GO TO 285 J3(3,2)=M12 285 J3(3,3)=M7 ! IHSH*2 C C --- SET ROWS 4,5,.. UP TO IHSH+1 C IF(M4.LT.4) GO TO 203 DO 470 J=4,M4 J2(J,1) = M4+J-4 J2(J,2)=M10 IF(ISIG.EQ.J-1) GO TO 472 J2(J,2)=J-1 IF(M1.EQ.0 .OR. IRHO.NE.J-1) GO TO 472 J2(J,2)=M9 472 J2(J,3) = M4+J-3 J3(J,1)=M7+J-4 J3(J,2)=M12 IF(ISIGP.EQ.J-1) GO TO 470 J3(J,2)=J-1 IF(M2.EQ.0 .OR. IRHOP.NE.J-1) GO TO 470 J3(J,2)=M11 470 J3(J,3)=M7+J-3 J3(M4,3)=J2(M4,3) ! 2*M4-3 C C --- SET FIRST TWO ROWS, CORRESPONDING TO COUPLING OF INTERACTING C ELECTRONS WITHIN THEIR SHELLS C 203 J2(2,3) = M10 J2(1,2) = M13 J2(2,2) = M14 J2(1,3) = M9 IF(M1.EQ.0) GO TO 83 J2(1,1) = IRHO J2(2,1) = ISIG GO TO 84 83 J2(1,1) = ISIG J2(2,1) = M9 84 J3(2,3) = M12 J3(1,2) = M13 J3(2,2) = M14 J3(1,3) = M11 IF(M2.EQ.0) GO TO 86 J3(1,1) = IRHOP J3(2,1) = ISIGP GO TO 187 86 J3(1,1) = ISIGP J3(2,1) = M11 C C --- STORE J2,J3 ARRAYS FOR USE IN CALCULATING EXCHANGE INTEGRAL C 187 I1=0 DO 451 J=1,M4 DO 452 K=1,3 J2SPIN(K+I1)=J2(J,K) 452 J3SPIN(K+I1)=J3(J,K) 451 I1=I1+3 JSNDIR=1 3 IF(IBUG1.LE.1) GO TO 570 C C PRINT VALUES IN NJSYM IF IBUG3=1 C IF(IBUG3.NE.1) GO TO 570 WRITE(IWRITE,304) ((J2(J,K),K=1,3),(J3(J,K),K=1,3),J=1,M4) 570 RETURN C C --- SET J2 AND J3 ARRAYS FROM STORE OF PREVIOUS CALCULATIONS C 1 I1=0 DO 4 J=1,M4 DO 5 K=1,3 J2(J,K)=J2SPIN(K+I1) 5 J3(J,K)=J3SPIN(K+I1) 4 I1=I1+3 GO TO 3 END SUBROUTINE MODJ23(K) C C === MODIFIES THE DIRECT J2 AND J3 ARRAYS FOR EXCHANGE CALL OF NJSYM C LOGICAL FREE PARAMETER (LL74= 21+4, LL76= 21*3+19) COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15, 1 M16,M17,M18,M19,M20 COMMON/COUPLE/NJ1S,NJ23S,J1(LL76),J2(LL74,3),J3(LL74,3),FREE(LL76) COMMON/HOLD/J2SPIN(LL76),J3SPIN(LL76),J2ANG(LL76),J3ANG(LL76) C 8 FORMAT(3H J2,18X,2HJ3,16X,'from MODJ23'/(3I5,5X,3I5)) IF (K.GT.1) GO TO 2 C C --- K=1 - SPIN INTEGRALS C MK=M4 I1=0 DO 11 J=1,MK DO 12 K=1,3 J2(J,K)=J2SPIN(K+I1) 12 J3(J,K)=J3SPIN(K+I1) 11 I1=I1+3 J3(1,2)=M14 J3(2,2)=M13 GO TO 3 C C --- K=2 - ANGULAR INTEGRALS C 2 MK=M5 I1=0 DO 21 J=1,MK DO 22 K=1,3 J2(J,K)=J2ANG(K+I1) 22 J3(J,K)=J3ANG(K+I1) 21 I1=I1+3 J2(1,1)=M15 J3(1,3)=M16 3 IF(IBUG1.LE.1) GO TO 4 C C PRINT-OUT OF VALUES IN NJSYM IF IBUG3=1 C IF(IBUG3.EQ.1 ) GO TO 4 WRITE(IWRITE,8) ((J2(J,K),K=1,3),(J3(J,K),K=1,3),J=1,MK) 4 RETURN END SUBROUTINE J23ANG(IRHO,ISIG,IRHOP,ISIGP,JANGDI) C C === SETS UP J2 AND J3 ARRAYS FOR DIRECT ANGULAR INTEGRAL CALL OF NJSYM C LOGICAL FREE PARAMETER (LL74= 21+4, LL76= 21*3+19) COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15, 1 M16,M17,M18,M19,M20 C DELETE CARD NUMBER 1750 - COMMON BLOCK XATION NOT NEEDED COMMON/COUPLE/NJ1S,NJ23S,J1(LL76),J2(LL74,3),J3(LL74,3),FREE(LL76) COMMON/HOLD/J2SPIN(LL76),J3SPIN(LL76),J2ANG(LL76),J3ANG(LL76) C 304 FORMAT(3H J2,18X,2HJ3/(3I5,5X,3I5)) C C HAVE THE J2 AND J3 ARRAYS ALREADY BEEN SET. IF NOT, THEN GO TO 2 C IF(JANGDI) 2,2,1 C C --- ROWS 3 TO M4 OF SPIN J2 AND J3 ARE SAME AS ROWS 4 TO (M4+1) OF C ANGULAR J2 AND J3 C 2 I1=6 DO 103 J=3,M4 DO 104 K=1,3 J2(J+1,K)=J2SPIN(K+I1) 104 J3(J+1,K)=J3SPIN(K+I1) 103 I1=I1+3 C C --- SET ROWS 1, 2 AND 3 C J2(3,1)=ISIG IF(M1.EQ.0) J2(3,1)=M9 J3(3,1)=ISIGP IF(M2.EQ.0) J3(3,1)=M11 J2(2,3)=M9 J2(2,1)=IRHO J2(2,2)=M13 J2(1,3)=M14 J2(3,2)=M14 J2(3,3)=M10 J2(1,1)=M16 J2(1,2)=M17 J3(3,2)=M16 J3(3,3)=M12 J3(1,2)=M13 J3(1,1)=M17 J3(1,3)=M15 J3(2,3)= M11 J3(2,1)=IRHOP J3(2,2)=M15 C C --- STORE J2 AND J3 FOR USE IN CALCULATING THE EXCHANGE TERM C I1=0 DO 535 J=1,M5 DO 536 K=1,3 J2ANG(K+I1)=J2(J,K) 536 J3ANG(K+I1)=J3(J,K) 535 I1=I1+3 JANGDI=1 3 IF(IBUG1.LE.1) GO TO 209 C C PRINT-OUT OF VALUES IN NJSYM IF IBUG3=1 C IF(IBUG3.EQ.1) GO TO 209 WRITE(IWRITE,304) ((J2(J,K),K=1,3),(J3(J,K),K=1,3),J=1,M5) 209 RETURN C C --- SET J2 AND J3 ARRAYS FROM STORE OF PREVIOUS CALCULATIONS C 1 I1=0 DO 4 J=1,M5 DO 5 K=1,3 J2(J,K)=J2ANG(K+I1) 5 J3(J,K)=J3ANG(K+I1) 4 I1=I1+3 GO TO 3 END BLOCK DATA C IMPLICIT REAL*8(A-H,O-Z) COMMON/TERMS/ NROWS,ITAB(18),J(18),N(189) C INSERT A COPY OF THE NEXT CARD AFTER CARDS 314,861,1285,1389, C 1481,1837,1883 AND 2257 PARAMETER (LL74= 21+4) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/BPSIZE/KFLN,KFL2,KFLM,KDUMMY(6) C C --- SET QUANTUM NUMBERS OF TERMS WHICH CAN BE FORMED FROM C CONFIGURATIONS L**Q. ONLY THE FIRST HALF OF THAT PART OF THE C TABLE, CORRESPONDING TO A GIVEN L, IS INCLUDED, BECAUSE OF THE C SYMMETRY OF THE TABLE. E.G. D**7 FORMS THE SAME TERMS AS D**3 C C THE ARRAYS ITAB,J,N CORRESPOND TO THE ARRAYS ITAB,JTAB,NTAB C DATA NROWS/18/ DATA ITAB/ 1, 1, 1, 3, 3, 1, 5, 8,16,16, 1, 1, 1, 1, 1, 1, 1, 1/ DATA J( 1),J( 2),J( 3),J( 4),J( 5),J( 6)/ 0, 3, 6, 9, 18, 27/ DATA J( 7),J( 8),J( 9),J(10),J(11),J(12)/ 30, 45, 69,117,165,168/ DATA J(13),J(14),J(15),J(16),J(17),J(18)/171,174,177,180,183,186/ DATA N( 1),N( 2),N( 3),N( 4),N( 5),N( 6)/ 1, 1, 2, 0, 1, 1/ DATA N( 7),N( 8),N( 9),N( 10),N( 11),N( 12)/ 1, 3, 2, 0, 1, 1/ DATA N( 13),N( 14),N( 15),N( 16),N( 17),N( 18)/ 2, 5, 1, 2, 3, 3/ DATA N( 19),N( 20),N( 21),N( 22),N( 23),N( 24)/ 1, 3, 2, 3, 5, 2/ DATA N( 25),N( 26),N( 27),N( 28),N( 29),N( 30)/ 3, 1, 4, 1, 5, 2/ DATA N( 31),N( 32),N( 33),N( 34),N( 35),N( 36)/ 0, 1, 1, 2, 5, 1/ DATA N( 37),N( 38),N( 39),N( 40),N( 41),N( 42)/ 2, 9, 1, 2, 3, 3/ DATA N( 43),N( 44),N( 45),N( 46),N( 47),N( 48)/ 2, 7, 3, 1, 5, 2/ DATA N( 49),N( 50),N( 51),N( 52),N( 53),N( 54)/ 3, 3, 2, 3, 5, 2/ DATA N( 55),N( 56),N( 57),N( 58),N( 59),N( 60)/ 3, 7, 2, 3, 9, 2/ DATA N( 61),N( 62),N( 63),N( 64),N( 65),N( 66)/ 3,11, 2, 3, 3, 4/ DATA N( 67),N( 68),N( 69),N( 70),N( 71),N( 72)/ 3, 7, 4, 0, 1, 1/ DATA N( 73),N( 74),N( 75),N( 76),N( 77),N( 78)/ 2, 5, 1, 2, 9, 1/ DATA N( 79),N( 80),N( 81),N( 82),N( 83),N( 84)/ 2, 3, 3, 2, 7, 3/ DATA N( 85),N( 86),N( 87),N( 88),N( 89),N( 90)/ 4, 1, 1, 4, 5, 1/ DATA N( 91),N( 92),N( 93),N( 94),N( 95),N( 96)/ 4, 7, 1, 4, 9, 1/ DATA N( 97),N( 98),N( 99),N(100),N(101),N(102)/ 4,13, 1, 4, 3, 3/ DATA N(103),N(104),N(105),N(106),N(107),N(108)/ 4, 5, 3, 4, 7, 3/ DATA N(109),N(110),N(111),N(112),N(113),N(114)/ 4, 9, 3, 4,11, 3/ DATA N(115),N(116),N(117),N(118),N(119),N(120)/ 4, 5, 5, 1, 5, 2/ DATA N(121),N(122),N(123),N(124),N(125),N(126)/ 3, 3, 2, 3, 5, 2/ DATA N(127),N(128),N(129),N(130),N(131),N(132)/ 3, 7, 2, 3, 9, 2/ DATA N(133),N(134),N(135),N(136),N(137),N(138)/ 3,11, 2, 3, 3, 4/ DATA N(139),N(140),N(141),N(142),N(143),N(144)/ 3, 7, 4, 5, 1, 2/ DATA N(145),N(146),N(147),N(148),N(149),N(150)/ 5, 5, 2, 5, 7, 2/ DATA N(151),N(152),N(153),N(154),N(155),N(156)/ 5, 9, 2, 5,13, 2/ DATA N(157),N(158),N(159),N(160),N(161),N(162)/ 5, 5, 4, 5, 9, 4/ DATA N(163),N(164),N(165),N(166),N(167),N(168)/ 5, 1, 6, 1, 7, 2/ DATA N(169),N(170),N(171) / 1, 9, 2/ DATA N(172),N(173),N(174),N(175),N(176),N(177)/ 1,11, 2, 1,13, 2/ DATA N(178),N(179),N(180),N(181),N(182),N(183)/ 1,15, 2, 1,17, 2/ DATA N(184),N(185),N(186),N(187),N(188),N(189)/ 1,19, 2, 1,21, 2/ C INSERT THE NEXT SEVEN CARDS AFTER CARD NUMBER 1880 C C SET GLOBAL REAL CONSTANTS C DATA ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS/ 1 0.0,0.1,0.5,1.0,2.0,3.0,4.0, 7.0,11.0,1.0E-05/ C C SET COMPILE PARAMETERS USED FOR RECUPD, IN /MDEFN/ AND /MSTATE/: DATA KFLN,KFL2,KFLM/ 15,LL74,2500/ C C C SET BLOCK DATA FOR CFPD SUBROUTINE C COMMON/FRPAR2/I(719) C DATA I( 1),I( 2),I( 3),I( 4),I( 5),I( 6),I( 7),I( 8), 1 I( 9),I( 10),I( 11),I( 12),I( 13),I( 14),I( 15),I( 16), 1 I( 17),I( 18),I( 19),I( 20),I( 21),I( 22),I( 23),I( 24), 2 I( 25),I( 26),I( 27),I( 28),I( 29),I( 30),I( 31),I( 32), 3 I( 33),I( 34),I( 35),I( 36),I( 37),I( 38),I( 39),I( 40), 4 I( 41),I( 42),I( 43),I( 44),I( 45),I( 46),I( 47),I( 48), 5 I( 49),I( 50),I( 51),I( 52),I( 53),I( 54),I( 55),I( 56), 6 I( 57),I( 58),I( 59),I( 60),I( 61),I( 62),I( 63),I( 64), 7 I( 65),I( 66),I( 67),I( 68),I( 69),I( 70),I( 71),I( 72), 8 I( 73),I( 74),I( 75),I( 76),I( 77),I( 78),I( 79),I( 80)/ 1 1, 5, 8, 16, 16, 1, 2, 3, 1 4, 5, 0, 2, 3, 4, 5, 0, 1 2, 3, 4, 3, 0, 2, 3, 2, 2 5, 0, 0, 3, 4, 3, 0, 0, 3 1, 4, 5, 0, 0, 3, 2, 3, 4 0, 0, 3, 4, 3, 0, 0, 0, 5 4, 5, 0, 0, 0, 2, 3, 0, 6 0, 0, 4, 5, 0, 0, 0, 4, 7 1, 0, 0, 0, 2, 3, 0, 0, 8 0, 4, 5, 0, 0, 0, 0, 3/ DATA I( 81),I( 82),I( 83),I( 84),I( 85),I( 86),I( 87),I( 88), 1 I( 89),I( 90),I( 91),I( 92),I( 93),I( 94),I( 95),I( 96), 1 I( 97),I( 98),I( 99),I(100),I(101),I(102),I(103),I(104), 2 I(105),I(106),I(107),I(108),I(109),I(110),I(111),I(112), 3 I(113),I(114),I(115),I(116),I(117),I(118),I(119),I(120), 4 I(121),I(122),I(123),I(124),I(125),I(126),I(127),I(128), 5 I(129),I(130),I(131),I(132),I(133),I(134),I(135),I(136), 6 I(137),I(138),I(139),I(140),I(141),I(142),I(143),I(144), 7 I(145)/ 1 0, 0, 0, 4, 5, 2, 3, 3, 1 2, 0, 0, 1, 1, 5, 4, 0, 1 4, 5, 4, 3, 0, 2, 4, 3, 2 2, 0, 0, 3, 3, 1, 0, 0, 3 2, 2, 6, 0, 0, 2, 1, 5, 4 0, 0, 1, 1, 4, 0, 0, 0, 5 6, 4, 0, 0, 0, 4, 3, 0, 6 0, 0, 4, 3, 0, 0, 0, 3, 7 2/ DATA I(146),I(147),I(148),I(149),I(150),I(151),I(152),I(153), 1 I(154),I(155),I(156),I(157),I(158),I(159),I(160),I(161), 1 I(162),I(163),I(164),I(165),I(166),I(167),I(168),I(169), 2 I(170),I(171),I(172),I(173),I(174),I(175),I(176),I(177), 3 I(178),I(179),I(180),I(181),I(182),I(183),I(184),I(185), 4 I(186),I(187),I(188),I(189),I(190),I(191),I(192),I(193), 5 I(194),I(195),I(196),I(197),I(198),I(199),I(200),I(201), 6 I(202),I(203),I(204),I(205),I(206),I(207),I(208),I(209), 7 I(210),I(211),I(212),I(213),I(214),I(215),I(216),I(217), 8 I(218),I(219),I(220),I(221),I(222),I(223),I(224),I(225)/ 1 0, 0, 0, 2, 2, 0, 0, 0, 1 2, 2, 0, 0, 0, 0, 1, 0, 1 0, 0, 0, 0, 2, 3, 4, 5, 2 6, 0, 3, 4, 3, 4, 0, 1, 3 2, 3, 4, 0, 1, 2, 3, 4, 4 0, 1, 2, 3, 4, 0, 0, 2, 5 3, 2, 0, 0, 2, 3, 2, 0, 6 0, 2, 3, 2, 0, 0, 0, 1, 7 2, 0, 0, 0, 1, 2, 0, 0, 8 0, 1, 2, 0, 0, 0, 1, 2/ DATA I(226),I(227),I(228),I(229),I(230),I(231),I(232),I(233), 1 I(234),I(235),I(236),I(237),I(238),I(239),I(240),I(241), 1 I(242),I(243),I(244),I(245),I(246),I(247),I(248),I(249), 2 I(250),I(251),I(252),I(253),I(254),I(255),I(256),I(257), 3 I(258),I(259),I(260),I(261),I(262),I(263),I(264),I(265), 4 I(266),I(267),I(268),I(269),I(270),I(271),I(272),I(273), 5 I(274),I(275),I(276),I(277),I(278),I(279),I(280),I(281), 6 I(282),I(283),I(284),I(285),I(286),I(287),I(288),I(289), 7 I(290)/ 1 0, 0, 0, 1, 2, 0, 0, 0, 1 1, 2, 0, 0, 0, 1, 2, 0, 1 0, 0, 1, 2, 1, 1, 1, 1, 2 1, 4, -7, -1, 21, 7, -21, 21, 3 -8, -1, -8, 0, 0, 28, -9, -49, 4 7, 0, 0, 1, 11, -25, -9, -25, 5 0, 0, 0, 0, -10, -10, -5, 45, 6 15, 0, 0, 0, 0, 0, 16, 0, 7 0/ DATA I(291),I(292),I(293),I(294),I(295),I(296),I(297),I(298), 1 I(299),I(300),I(301),I(302),I(303),I(304),I(305),I(306), 1 I(307),I(308),I(309),I(310),I(311),I(312),I(313),I(314), 2 I(315),I(316),I(317),I(318),I(319),I(320),I(321),I(322), 3 I(323),I(324),I(325),I(326),I(327),I(328),I(329),I(330), 4 I(331),I(332),I(333),I(334),I(335),I(336),I(337),I(338), 5 I(339),I(340),I(341),I(342),I(343),I(344),I(345),I(346), 6 I(347),I(348),I(349),I(350),I(351),I(352),I(353),I(354), 7 I(355),I(356),I(357),I(358),I(359),I(360),I(361),I(362), 8 I(363),I(364),I(365),I(366),I(367),I(368),I(369),I(370)/ 1 7, 20, -560, 224, -112, -21, -56, 16, 1 0, 0, 0, 0, 0, 0, 0, 0, 1 3, 0, 0, -56, -448, 49, -64, -14, 2 0, 0, 0, 0, 0, 0, 0, 0, 3 0, 26, 308, 110, 220, 0, 0, 0, 4 7, -154, -28, -132, 0, 0, 0, 0, 5 0, -9, 297, 90, -405, 45, 0, 0, 6 3, 66, -507, -3, -60, 15, 0, 0, 7 0, 5, 315, -14, -175, -21, -56, -25, 8 0, 70, 385, -105, 28, 63, 0, 0/ DATA I(371),I(372),I(373),I(374),I(375),I(376),I(377),I(378), 1 I(379),I(380),I(381),I(382),I(383),I(384),I(385),I(386), 1 I(387),I(388),I(389),I(390),I(391),I(392),I(393),I(394), 2 I(395),I(396),I(397),I(398),I(399),I(400),I(401),I(402), 3 I(403),I(404),I(405),I(406),I(407),I(408),I(409),I(410), 4 I(411),I(412),I(413),I(414),I(415),I(416),I(417),I(418), 5 I(419),I(420),I(421),I(422),I(423),I(424),I(425),I(426), 6 I(427),I(428),I(429),I(430),I(431),I(432),I(433),I(434), 7 I(435)/ 1 0, 0, 0, 315, 0, 0, 135, 0, 1 0, 189, 0, 0, 105, 0, 1, 0, 1 0, 0, 200, 15, 120, 60, -35, 10, 2 0, -25, 88, 200, 45, 20, 0, 1, 3 0, 0, 0, 16, -200, -14, -14, 25, 4 0, 0, 0, 120, -42, 42, 0, 0, 5 1, -105, -175, -175, -75, 0, 0, 0, 6 0, 0, 0, 0, 0, 0, 0, 0, 7 0/ DATA I(436),I(437),I(438),I(439),I(440),I(441),I(442),I(443), 1 I(444),I(445),I(446),I(447),I(448),I(449),I(450),I(451), 1 I(452),I(453),I(454),I(455),I(456),I(457),I(458),I(459), 2 I(460),I(461),I(462),I(463),I(464),I(465),I(466),I(467), 3 I(468),I(469),I(470),I(471),I(472),I(473),I(474),I(475), 4 I(476),I(477),I(478),I(479),I(480),I(481),I(482),I(483), 5 I(484),I(485),I(486),I(487),I(488),I(489),I(490),I(491), 6 I(492),I(493),I(494),I(495),I(496),I(497),I(498),I(499), 7 I(500),I(501),I(502),I(503),I(504),I(505),I(506),I(507), 8 I(508),I(509),I(510),I(511),I(512),I(513),I(514),I(515)/ 1 154, -110, 0, 0, 231, 286, 924, -308, 1 220, -396, 0, 0, 0, 0, 0, 0, 1 -66, -90, 180, 0, 99, -99, 891,-5577, 2 -405, -9, 0, 45, 45, 0, 0, 0, 3 0, 224, 0, -56, 0, -220, 1680, 0, 4 112, 0, -21, 21, 0, -16, 0, 0, 5 -70, 14, -84, 56, 0, 55, 945, 4235, 6 -175, -315, 0, -21, 189, -25, 0, 0, 7 25, -15, -135, 35, 0, 0, 600, 968, 8 120, 600, 0, 60, 60, 10, 3, 0/ DATA I(516),I(517),I(518),I(519),I(520),I(521),I(522),I(523), 1 I(524),I(525),I(526),I(527),I(528),I(529),I(530),I(531), 1 I(532),I(533),I(534),I(535),I(536),I(537),I(538),I(539), 2 I(540),I(541),I(542),I(543),I(544),I(545),I(546),I(547), 3 I(548),I(549),I(550),I(551),I(552),I(553),I(554),I(555), 4 I(556),I(557),I(558),I(559),I(560),I(561),I(562),I(563), 5 I(564),I(565),I(566),I(567),I(568),I(569),I(570),I(571), 6 I(572),I(573),I(574),I(575),I(576),I(577),I(578),I(579), 7 I(580)/ 1 0, -56, 0, -64, 0, 0, 0, 0, 1 448, 0, -9, -49, 0, 14, 0, 0, 1 0, -16, 126, 14, 0, 0, 0, 0, 2 -200, 360, 0, -14, 126, 25, 0, 0, 3 0, 0, 0, 0, -175, 182, -728,-2184, 4 0, 0, 0, 0, 0, 0, 0, 0, 5 0, 0, 0, 0, 0, 220, 880, 0, 6 -400, 0, -9, -25, 0, 0, 0, 0, 7 0/ DATA I(581),I(582),I(583),I(584),I(585),I(586),I(587),I(588), 1 I(589),I(590),I(591),I(592),I(593),I(594),I(595),I(596), 1 I(597),I(598),I(599),I(600),I(601),I(602),I(603),I(604), 2 I(605),I(606),I(607),I(608),I(609),I(610),I(611),I(612), 3 I(613),I(614),I(615),I(616),I(617),I(618),I(619),I(620), 4 I(621),I(622),I(623),I(624),I(625),I(626),I(627),I(628), 5 I(629),I(630),I(631),I(632),I(633),I(634),I(635),I(636), 6 I(637),I(638),I(639),I(640),I(641),I(642),I(643),I(644), 7 I(645),I(646),I(647),I(648),I(649),I(650),I(651),I(652), 8 I(653),I(654),I(655),I(656),I(657),I(658),I(659),I(660)/ 1 0, 0, 0, -45, -5, 845,-1215, 275, 1 495, 0, -11, 99, 0, 0, 0, 0, 1 0, 0, 0, 0, 33, -7,-2541, 105, 2 -525, 0, 35, 35, -15, 0, 0, 0, 3 0, 0, 0, 0, 0, -800, 0, -160, 4 0, -5, 45, 0, 30, 0, 0, 0, 5 0, 0, 0, 0, 0, -100, 1452, 180, 6 -100, 0, -10, 90, 15, -2, 0, 0, 7 0, 0, 0, 0, 0, 0, 0, 0, 8 0, 6, 0, 0, 0, 0, 0, 0/ DATA I(661),I(662),I(663),I(664),I(665),I(666),I(667),I(668), 1 I(669),I(670),I(671),I(672),I(673),I(674),I(675),I(676), 1 I(677),I(678),I(679),I(680),I(681),I(682),I(683),I(684), 2 I(685),I(686),I(687),I(688),I(689),I(690),I(691),I(692), 3 I(693),I(694),I(695),I(696),I(697),I(698),I(699),I(700), 4 I(701),I(702),I(703),I(704),I(705),I(706),I(707),I(708), 5 I(709),I(710),I(711),I(712),I(713),I(714),I(715),I(716), 6 I(717),I(718),I(719)/ 1 0, 0, 0, 0, 0, 0, 0, 0, 1 0, 0, -14, -56, 0, 0, 1, 1, 1 1, 1, 1, 5, 15, 2, 42, 70, 2 60, 140, 30, 10, 60, 1680, 840, 1680, 3 210, 360, 90, 10, 504, 1008, 560, 280, 4 140, 1, 1, 1, 420, 700, 700, 300, 5 550, 1100, 8400,18480, 2800, 2800, 50, 350, 6 700, 150, 5/ C END SUBROUTINE INTACT(L,LP,IEQUIV) C IMPLICIT REAL*8(A-H,O-Z) COMMON/INFORM/IREAD,IWRITE,IPUNCH C REPLACE CARDS 1886,2039 AND 2260 BY A COPY OF THE FOLLOWING CARD COMMON/ENAV/COEFCT(5),NINTS,KVALUE(5) C C INSERT THE FOLLOWING SIX CARDS AFTER CARD NUMBER 2040 C ********************************************************* C * * C * THIS SUBROUTINE CONTAINS REAL CONSTANTS IN THE TEXT * C * * C ********************************************************* C C THIS SUBROUTINE GIVES THE INTERACTION ENERGY BETWEEN TWO SHELLS, C ONE WITH ORBITAL ANGULAR MOMENTUM L , THE OTHER WITH ORBITAL C ANGULAR MOMENTUM LP. NOTICE THAT THE FIRST TERM OF THIS C INTERACTION ENERGY IS ALWAYS F0(L,LP) AND THIS IS NOT GIVEN C IN THIS SUBROUTINE. THUS ONLY THE EXTRA TERMS ARE HERE PRODUCED. C FOR EQUIVALENT ELECTRONS (IEQUIV = 1), THERE WILL BE FK C INTEGRALS ONLY. FOR NON-EQUIVALENT ELECTRONS (IEQUIV = 2) C THERE WILL BE GK INTEGRALS ONLY. C C THE EXPRESSIONS FOR THE INTERACTION ENERGIES INVOLVING SHELLS WITH C L.LE.3 ARE GIVEN BY J.C. SLATER, QUANTUM THEORY OF ATOMIC C STRUCTURE, VOL. I, EQUATIONS (14.20) AND (14.22). IN THE LAST C OF HIS EQUATIONS (14.22) A TERM -1/14 G0(F,FP) IS OMITTED. C THIS IS INCLUDED BELOW. THE INTERACTION ENERGIES FOR G-ELECTRON C SHELLS MAY BE EVALUATED USING HIS EQUATIONS (13.12), (13.17), C (14.19), AND (14,21) C IF(L.GT.4.OR.LP.GT.4) GO TO 3 IF(IEQUIV.GT.1) GO TO 2 C C === EQUIVALENT ELECTRONS C IF(L.NE.LP) GO TO 4 GO TO 5 4 WRITE(IWRITE,11) 11 FORMAT(85H ERROR IN INTACT - EQUIVALENT ELECTRONS CALLED FOR DIFFE *RING ORBITAL ANGULAR MOMENTUM) STOP 5 IF(L.GT.0) GO TO 12 C C S ELECTRONS C NINTS=0 RETURN 12 GO TO (13,14,15,16),L C C P ELECTRONS C 13 NINTS=1 KVALUE(1)=2 COEFCT(1)=-0.08 RETURN C C D ELECTRONS C 14 NINTS=2 KVALUE(1)=2 KVALUE(2)=4 COEFCT(1)=-0.031746032 COEFCT(2)=-0.031746032 RETURN C C F ELECTRONS C 15 NINTS=3 KVALUE(1)=2 KVALUE(2)=4 KVALUE(3)=6 COEFCT(1)=-0.020512821 COEFCT(2)=-0.013986014 COEFCT(3)=-0.017930787 RETURN C C G ELECTRONS C 16 NINTS=4 KVALUE(1)=2 KVALUE(2)=4 KVALUE(3)=6 KVALUE(4)=8 COEFCT(1)=-0.015278839 COEFCT(2)=-0.009519892 COEFCT(3)=-0.008227067 COEFCT(4)=-0.011856655 RETURN C C --- NON-EQUIVALENT ELECTRONS C 2 IF(L.GT.LP) GO TO 21 L1=L L2=LP GO TO 22 21 L1=LP L2=L 22 L1D=L1+1 L2D=L2+1 GO TO (30,40,50,60,70),L1D 30 NINTS=1 KVALUE(1)=L2 GO TO (31,32,33,34,35),L2D C C S - S INTERACTION C 31 COEFCT(1)=-0.5 RETURN C C S - P INTERACTION C 32 COEFCT(1)=-0.166666667 RETURN C C S - D INTERACTION C 33 COEFCT(1)=-0.1 RETURN C C S - F INTERACTION C 34 COEFCT(1)=-0.071428571 RETURN C C S - G INTERACTION C 35 COEFCT(1)=-0.055555556 RETURN 40 NINTS=2 KVALUE(1)=L2-1 KVALUE(2)=L2D GO TO (41,42,43,44),L2 C C P - P INTERACTION C 41 COEFCT(1)=-0.166666667 COEFCT(2)=-0.066666667 RETURN C C P - D INTERACTION C 42 COEFCT(1)=-0.066666667 COEFCT(2)=-0.042857143 RETURN C C P - F INTERACTION C 43 COEFCT(1)=-0.042857143 COEFCT(2)=-0.031746032 RETURN C C P - G INTERACTION C 44 COEFCT(1)=-0.031746032 COEFCT(2)=-0.025252525 RETURN 50 NINTS=3 KVALUE(1)=L2-2 KVALUE(2)=L2 KVALUE(3)=L2+2 IF(L2-3) 51,52,53 C C D - D INTERACTION C 51 COEFCT(1)=-0.1 COEFCT(2)=-0.028571429 COEFCT(3)=-0.028571429 RETURN C C D - F INTERACTION C 52 COEFCT(1)=-0.042857143 COEFCT(2)=-0.019047619 COEFCT(3)=-0.021645022 RETURN C C D - G INTERACTION C 53 COEFCT(1)=-0.028571429 COEFCT(2)=-0.014430014 COEFCT(3)=-0.017482517 RETURN 60 NINTS=4 KVALUE(1)=L2-3 KVALUE(2)=L2-1 KVALUE(3)=L2+1 KVALUE(4)=L2+3 IF(L2.GT.3) GO TO 62 C C F - F INTERACTION C COEFCT(1)=-0.071428571 COEFCT(2)=-0.019047619 COEFCT(3)=-0.012987013 COEFCT(4)=-0.016650017 RETURN C C F - G INTERACTION C 62 COEFCT(1)=-0.031746032 COEFCT(2)=-0.012987013 COEFCT(3)=-0.009990010 COEFCT(4)=-0.013597514 RETURN C C G - G INTERACTION C 70 NINTS=5 KVALUE(1)=0 KVALUE(2)=2 KVALUE(3)=4 KVALUE(4)=6 KVALUE(5)=8 COEFCT(1)=-0.055555556 COEFCT(2)=-0.014430014 COEFCT(3)=-0.008991009 COEFCT(4)=-0.007770008 COEFCT(5)=-0.011197952 RETURN C C --- IF ANGULAR MOMENTUM VALUES ARE TOO LARGE ----- C 3 WRITE(IWRITE,6) L,LP 6 FORMAT(//47H THE ORBITAL ANGULAR MOMENTUM VALUES, WHICH ARE,2I5/ * 9X,38HARE TOO LARGE FOR THE CODING OF INTACT//) STOP END SUBROUTINE CHOP C IMPLICIT REAL*8(A-H,O-Z) PARAMETER (LL43= 21*2+3, LL75= 21+2) COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75), C REPLACE CARDS 548,577,673,867,1494,2262,2309,2355 BY A COPY OF 1 J1QN1(LL43,3),J1QN2(LL43,3),IJFUL(LL75) C REPLACE CARD NUMBERS 179,1888,2311 BY A COPY OF THE NEXT CARD COMMON/DIAGNL/IDIAG,JA,JB DIMENSION ICHOP(LL75) SAVE ICHOP ! /REMOVE/ C C --- ZEROIZE THE OUTPUT ARRAY C DO 1 I=1,IHSH 1 ICHOP(I)=0 C C NO AVERAGE ENERGY TERMS FOR OFF-DIAGONAL MATRIX ELEMENTS C IF(IDIAG.EQ.0) RETURN JSTO=0 ICOUNT=0 DO 3 J=1,IHSH NFULL=4*LJ(J)+2 I2=NOSH1(J) C C IS THE SHELL FULL OR EMPTY C IF(I2.EQ.NFULL.OR.I2.EQ.0) GO TO 4 C C IF NOT, DOES IT CONTAIN ONLY ONE ELECTRON, OR ONLY ONE =HOLE= C IF(I2.EQ.1.OR.I2.EQ.NFULL-1) JSTO=J GO TO 3 4 ICOUNT=ICOUNT+1 C C ICHOP SET UNITY FOR CLOSED SHELLS C ICHOP(J)=1 3 CONTINUE C C IF ALL BUT ONE SHELL IS CLOSED, AND THIS CONTAINS ONE ELECTRON OR C =HOLE= , THEN IT CAN BE TREATED PURELY BY AVERAGE ENERGY C IF (ICOUNT.NE.IHSH-1 .OR. JSTO.EQ.0) RETURN ICHOP(JSTO)=1 RETURN END SUBROUTINE REDUCE(IRHO,ISIG,IRHOP,ISIGP,LESSEN) C C THIS SUBROUTINE REMOVES SPECTATOR SINGLET S SHELLS, WHICH HAVE C NO EFFECT IN ANGULAR OR SPIN INTEGRALS C PARAMETER (LL43= 21*2+3, LL75= 21+2) DIMENSION LEAVE(LL75) ! '05Feb8: not (23) -- VKL's Li+/ 21.ge.24! COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75), C REPLACE CARDS 548,577,673,867,1494,2262,2309,2355 BY A COPY OF 1 J1QN1(LL43,3),J1QN2(LL43,3),IJFUL(LL75) C C LMIN INITIALLY SET LARGE C LMIN=99 ICOUNT=0 DO 1 I=1,IHSH C C NO INTERACTING SHELL MAY BE REMOVED C IF(I.EQ.IRHO.OR.I.EQ.ISIG.OR.I.EQ.IRHOP.OR.I.EQ.ISIGP) GO TO 2 C C IF A SPECTATOR SHELL HAS SINGLET S COUPLING ON BOTH SIDES OF C THE MATRIX ELEMENT, IT MAY IN GENERAL BE REMOVED, AS IT HAS NO C EFFECT IN FANO C IF(J1QN1(I,1).EQ.0.AND.J1QN2(I,1).EQ.0) GO TO 7 2 ICOUNT=ICOUNT+1 LEAVE(ICOUNT)=I GO TO 1 7 IF(LJ(I).GE.LMIN) GO TO 1 LMIN=LJ(I) ILMIN=I 1 CONTINUE IF(ICOUNT.EQ.IHSH) GO TO 8 C C IF A CHANGE IN THE COMMON BLOCK MEDEFN IS TO BE MADE, C ITS PRESENT SITUATION MUST BE PRESERVED BY A CALL OF MEKEEP C CALL MEKEEP(IRHO,ISIG,IRHOP,ISIGP) C C IF ONLY ONE SHELL WOULD BE LEFT IN THIS WAY, THE ONE DESTINED FOR C REMOVAL, WITH THE LOWEST L-VALUE, MUST BE RETAINED TO DEFINE A C COUPLING C IF(ICOUNT.EQ.1) GO TO 10 C C --- MODIFY THE COMMON BLOCK MEDEFN C C INSERT THE FOLLOWING CARD AFTER CARD 2396 13 DO 3 I=1,ICOUNT J=LEAVE(I) IF(J.EQ.IRHO) IRHO=I IF(J.EQ.ISIG) ISIG=I IF(J.EQ.IRHOP) IRHOP=I IF(J.EQ.ISIGP) ISIGP=I NJ(I)=NJ(J) LJ(I)=LJ(J) NOSH1(I)=NOSH1(J) NOSH2(I)=NOSH2(J) DO 4 K=1,3 J1QN1(I,K)=J1QN1(J,K) 4 J1QN2(I,K)=J1QN2(J,K) 3 CONTINUE ISUBH=IHSH-1 DO 5 I=2,ICOUNT J=LEAVE(I) II=ICOUNT+I-1 IJ=ISUBH+J DO 6 K=1,3 J1QN1(II,K)=J1QN1(IJ,K) 6 J1QN2(II,K)=J1QN2(IJ,K) 5 CONTINUE IHSH=ICOUNT GO TO 20 C C THIS SITUATION ONLY OCCURS IF IRHO=ISIG=IRHOP=ISIGP C 10 J=LEAVE(1) C REPLACE CARDS 2428 THROUGH 2449 BY THE FOLLOWING FIVE CARDS II=MIN(J,ILMIN) LEAVE(1)=II LEAVE(2)=J+ILMIN-II ICOUNT=2 GO TO 13 20 IF(IBUG1.LE.1) GO TO 9 I2HSH=IHSH+IHSH-1 WRITE(IWRITE,35) ((J1QN1(J,K),K=1,3),J=1,I2HSH) WRITE(IWRITE,36) ((J1QN2(J,K),K=1,3),J=1,I2HSH) 35 FORMAT(/43H REDUCE: NEW DEFINITION OF COUPLING SCHEMES/ * 9X,37HFOR THIS SET OF RHO, SIG, RHOP, SIGP// * 11X,47HL.H.S. OF HAMILTONIAN MATRIX ELEMENT DEFINED BY/ * 11X,4HJ1QN,(T18,3I3)) 36 FORMAT(11X,47HR.H.S. OF HAMILTONIAN MATRIX ELEMENT DEFINED BY/ * 11X,4HJ1QN,(T18,3I3)) C C LESSEN = 0 IF NO CHANGE IN MEDEFN C = 1 OTHERWISE C 9 LESSEN=1 RETURN 8 LESSEN=0 RETURN END SUBROUTINE MEKEEP(IRHO,ISIG,IRHOP,ISIGP) C C STORES THE COMMON BLOCK MEDEFN , AND IRHO,ISIG,IRHOP,ISIGP C PARAMETER (LL43= 21*2+3, LL75= 21+2) PARAMETER(MED=5*LL75+6*LL43+1) C REPLACE CARDS 2471 AND 2486 BY A COPY OF THE NEXT CARD COMMON/MEDEFN/J(MED) C REPLACE CARDS 2472 AND 2487 BY A COPY OF THE NEXT CARD COMMON/STORE/I(MED),I1,I2,I3,I4 C DO 1 K=1,MED 1 I(K)=J(K) I1=IRHO I2=ISIG I3=IRHOP I4=ISIGP RETURN END SUBROUTINE MEREST(IRHO,ISIG,IRHOP,ISIGP) C C RESTORES THE COMMON BLOCK MEDEFN, AND IRHO,ISIG,IRHOP,ISIGP PARAMETER (LL43= 21*2+3, LL75= 21+2) PARAMETER(MED=5*LL75+6*LL43+1) C C REPLACE CARDS 2471 AND 2486 BY A COPY OF THE NEXT CARD COMMON/MEDEFN/J(MED) C REPLACE CARDS 2472 AND 2487 BY A COPY OF THE NEXT CARD COMMON/STORE/I(MED),I1,I2,I3,I4 C DO 1 K=1,MED 1 J(K)=I(K) IRHO=I1 ISIG=I2 IRHOP=I3 ISIGP=I4 RETURN END C INSERT THE FOLLOWING SUBROUTINES - H0WTS,SETM,DH0,ODH0 - C AFTER CARD 2499 C*********************************************************************** SUBROUTINE SETM C C --- SET CONSTANTS USEFUL IN INNER SUBROUTINES C C IMPLICIT REAL*8(A-H,O-Z) PARAMETER (LL43= 21*2+3, LL75= 21+2) COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75), 1 J1QN1(LL43,3),J1QN2(LL43,3),IJFUL(LL75) COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15, 1 M16,M17,M18,M19,M20 M3=IHSH-1 M4=IHSH+1 M5=IHSH+2 M6=2*IHSH-1 M7=M6+1 M8=M3+M6 ! IHSH*3-2 M9=M8+1 M10=M8+2 M11=M8+3 M12=M8+4 M13=M8+5 M14=M8+6 M15=M8+7 M16=M8+8 M17=M8+9 M18=IHSH+3 RETURN END C C*********************************************************************** C C INSERT THE PROGRAM PSHELL C F P BY D C S ALLISON CATALOGUE C NUMBER ACQB (CPC 1 (1969) 15) AND THE PROGRAM A NEW D SHELL C F P C BY A T CHIVERS CATALOGUE NUMBER ACRN (CPC 6 (1973) 88) C C CFPP,BLOCK DATA,CFPD C C CARDS ACQB0050-ACQB0171,ACRN0065-ACRN0410 C C*********************************************************************** C SUBROUTINE CFPP(N,LI,ISI,LJ,ISJ,COEFP) C IMPLICIT REAL*8(A-H,O-Z) C C THIS SUBROUTINE EVALUATES THE COEFFICIENTS OF FRACTIONAL PARENTAGE C FOR EQUIVALENT P SHELL ELECTRONS FROM TABLES GIVEN IN J.C.SLATER C QUANTUM THEORY OF ATOMIC STRUCTURE,VOLUME2,P350(1960) C IN THE SUBROUTINE LIST N,THE NO. OF ELECTRONS,L THE ANGULAR C MOMENTUM QUANTUM NO.,(2S+1) THE SPIN QUANTUM NO. OF BOTH THE STATE C IN QUESTION AND ITS PARENT STATE ARE INPUT PARAMETERS. THE RESULT C IS OUTPUT AS COEFP C DIMENSION IL(3,3),IS(3,3),ITAB1(3,1),ITAB2(3,3),NORM1(3),NORM2(3) C C C SET UP P SHELL PARAMETERS AND TABLES C DATA IL(1,1),IL(2,1),IL(2,2),IL(2,3),IL(3,1),IL(3,2),IL(3,3)/1,1,2 1 ,0,0,2,1/ DATA IS(1,1),IS(2,1),IS(2,2),IS(2,3),IS(3,1),IS(3,2),IS(3,3)/2,3,1 1 ,1,4,2,2/ DATA ITAB1(1,1),ITAB1(2,1),ITAB1(3,1)/1,1,1/ DATA ITAB2(1,1),ITAB2(1,2),ITAB2(1,3),ITAB2(2,1),ITAB2(2,2),ITAB2( 1 2,3),ITAB2(3,1),ITAB2(3,2),ITAB2(3,3)/1,0,0,1,-1,0,-9,-5,4/ DATA NORM1(1),NORM1(2),NORM1(3)/1,1,1/ DATA NORM2(1),NORM2(2),NORM2(3)/1,2,18/ C C TEST IF N IS IN THE FIRST HALF OF SHELL C 99 IF(N.GE.4) GO TO 103 C C TEST IF STATE IN QUESTION IS ALLOWED C IF IT IS, IDENTIFY THE ROW OF THE TABLE BY J1 C J = 0 101 J = J+1 IF(J.GE.4) GO TO 8 IF(IL(N,J).NE.LI) GO TO 101 IF(IS(N,J).NE.ISI) GO TO 101 J1 = J C C TEST IF PARENT STATE IS ALLOWED C IF IT IS, IDENTIFY THE COLUMN OF THE TABLE BY J2 C IF(N.NE.1) GO TO 44 IF(LJ.NE.0) GO TO 8 IF(ISJ-1) 8,1,8 44 J = 0 102 J = J+1 IF(J.GE.4) GO TO 8 IF(IL(N-1,J).NE.LJ) GO TO 102 IF(IS(N-1,J).NE.ISJ) GO TO 102 J2=J GO TO 100 C C SIMILAR SETTING OF J1 AND J2 IF N IS IN SECOND HALF OF SHELL C 103 M =6-N IF(M) 72,73,72 73 IF(LI) 8,74,8 74 IF(ISI-1) 8,75,8 72 J = 0 104 J = J+1 IF(J.GE.4) GO TO 8 IF(IL(M,J).NE.LI) GO TO 104 IF(IS(M,J).NE.ISI) GO TO 104 J1 = J 75 J = 0 105 J = J+1 IF(J.GE.4) GO TO 8 IF(IL(M+1,J).NE.LJ) GO TO 105 IF(IS(M+1,J).NE.ISJ) GO TO 105 J2 = J C C C IDENTIFY THE F.P.C AS A UNIQUE ELEMENT OF ITABN(J1,J2) C 100 GO TO (1,2,3,4,4,1),N 1 COEFP = 1.0 GO TO 10 2 COEFP = ITAB1(J1,J2) IF(COEFP) 54,10,31 54 COEFP = -SQRT(-COEFP/NORM1(J1)) GO TO 10 31 COEFP = SQRT(COEFP/NORM1(J1)) GO TO 10 3 COEFP = ITAB2(J1,J2) IF(COEFP) 55,10,32 55 COEFP = -SQRT(-COEFP/NORM2(J1)) GO TO 10 32 COEFP =SQRT(COEFP/NORM2(J1)) GO TO 10 C C USE RECURRENCE RELATION EQUATION (19) OF RACAH FOR SECOND HALF OF C SHELL C 4 ISIGN = (-1)**((ISI+ISJ-5)/2+LI+LJ) FACTOR = ((7-N)*ISJ*(2*LJ+1))/REAL(N*ISI*(2*LI+1)) IF(N-5) 56,5,8 56 COEFP = ITAB2(J2,J1) IF(COEFP) 57,10,33 57 COEFP = -SQRT(-COEFP/NORM2(J2)) GO TO 34 33 COEFP = SQRT(COEFP/NORM2(J2)) 34 COEFP = COEFP * ISIGN * SQRT(FACTOR) IF(LJ-1) 35,10,35 35 COEFP = -COEFP GO TO 10 5 COEFP = ITAB1(J2,J1) IF(COEFP) 58,10,36 58 COEFP = -SQRT(-COEFP/NORM1(J2)) GO TO 37 36 COEFP = SQRT(COEFP/NORM1(J2)) 37 COEFP = COEFP * ISIGN * SQRT(FACTOR) GO TO 10 C C FOR AN UNALLOWED STATE THE F.P.COEFFICIENT IS SET EQUAL TO AN C ERRONEOUS VALUE. THIS STATEMENT COULD BE REPLACED BY AN ERROR C STATEMENT C 8 COEFP = 9.9 10 CONTINUE RETURN END SUBROUTINE CFPD(N,IVI,LI,ISI,IVJ,LJ,ISJ,COEFP) C PRECEEDING BLOCK DATA MERGED WITH FIRST BLOCK DATA - RUB'95JUN13 C C THIS SUBROUTINE EVALUATES THE COEFFICIENTS OF FRACTIONAL PARENTAGE C FOR EQUIVALENT D SHELL ELECTRONS FROM TABLES GIVEN IN J.C.SLATER C QUANTUM THEORY OF ATOMIC STRUCTURE,VOLUME2,P350(1960) C IN THE SUBROUTINE LIST N THE NO OF ELECTRONS V THE SENIORITY QUAN- C TUM NO, L THE ANGULAR MOMENTUM QUANTUM NO, (2S+1) THE SPIN QUANTUM C NO. OF BOTH THE STATE IN QUESTION AND ITS PARENT STATE ARE INPUT C PARAMETERS THE RESULT IS OUTPUT AS COEFP C C IMPLICIT REAL*8(A-H,O-Z) C COMMON/FRPAR2/K(5),IV(5,16),IL(5,16),IS(5,16),ITAB1(5,1),ITAB2(8,5 1 ),ITAB3(16,8),ITAB4(16,16),NORM1(5),NORM2(8),NORM3(16),NORM4(16) COMMON/INFORM/IREAD,IWRITE,IPUNCH C C C TEST IF N IS IN THE FIRST HALF OF SHELL C 99 IF(N.GE.6) GO TO 103 C C TEST IF STATE IN QUESTION IS ALLOWED C IF IT IS, IDENTIFY THE ROW OF THE TABLE BY J1 C J = 0 101 J = J+1 IF(J.GE.17) GO TO 11 IF(IV(N,J).NE.IVI) GO TO 101 IF(IL(N,J).NE.LI) GO TO 101 IF(IS(N,J).NE.ISI) GO TO 101 J1=J C C TEST IF PARENT STATE IS ALLOWED C IF IT IS, IDENTIFY THE COLUMN OF THE TABLE BY J2 C IF(N-1) 45,30,45 30 IF(IVJ) 11,31,11 31 IF(LJ) 11,32,11 32 IF(ISJ-1) 11,1,11 45 J = 0 102 J = J+1 IF(J.GE.17) GO TO 11 IF(IV(N-1,J).NE.IVJ) GO TO 102 IF(IL(N-1,J).NE.LJ) GO TO 102 IF(IS(N-1,J).NE.ISJ) GO TO 102 J2=J GO TO 100 C C SIMILAR SETTING OF J1 AND J2 IF N IS IN SECOND HALF OF SHELL C 103 M = 10-N IF(M) 36,33,36 33 IF(IVI) 11,34,11 34 IF(LI) 11,35,11 35 IF(ISI-1) 11,37,11 36 J = 0 104 J = J+1 IF(J.GE.17) GO TO 11 IF(IV(M,J).NE.IVI) GO TO 104 IF(IL(M,J).NE.LI) GO TO 104 IF(IS(M,J).NE.ISI) GO TO 104 J1=J 37 J = 0 105 J = J+1 IF(J.GE.17) GO TO 11 IF(IV(M+1,J).NE.IVJ) GO TO 105 IF(IL(M+1,J).NE.LJ) GO TO 105 IF(IS(M+1,J).NE.ISJ) GO TO 105 J2=J C C IDENTIFY THE F.P.C AS A UNIQUE ELEMENT OF ITABN(J1,J2) C 100 GO TO (1,2,3,4,5,12,12,12,12,1),N 1 COEFP = 1.0 GO TO 10 2 COEFP = ITAB1(J1,J2) IF(COEFP) 60,10,81 60 COEFP = - SQRT(-COEFP/NORM1(J1)) GO TO 10 81 COEFP = SQRT(COEFP/NORM1(J1)) GO TO 10 3 COEFP = ITAB2(J1,J2) IF(COEFP) 61,10,82 61 COEFP = -SQRT(-COEFP/NORM2(J1)) GO TO 10 82 COEFP = SQRT(COEFP/NORM2(J1)) GO TO 10 4 COEFP = ITAB3(J1,J2) IF(COEFP) 62,10,83 62 COEFP = -SQRT(-COEFP/NORM3(J1)) GO TO 10 83 COEFP = SQRT(COEFP/NORM3(J1)) GO TO 10 5 COEFP = ITAB4(J1,J2) IF(COEFP) 63,10,84 63 COEFP = -SQRT(-COEFP/NORM4(J1)) GO TO 10 84 COEFP = SQRT(COEFP/NORM4(J1)) GO TO 10 C C USE RECURRENCE RELATION EQUATION (19) OF RACAH FOR SECOND HALF OF C SHELL C 12 ISIGN = (-1)**((ISI+ISJ-7)/2 +LI +LJ) FACTOR = SQRT(((2*LJ+1)*(11-N)*ISJ)/REAL((2*LI+1)*ISI*N)) M1 =N-5 GO TO(6,7,8,9),M1 6 COEFP = ITAB4(J2,J1) IF(COEFP) 64,10,85 64 COEFP = -SQRT(-COEFP/NORM4(J2)) GO TO 86 85 COEFP = SQRT(COEFP/NORM4(J2)) 86 COEFP = COEFP*ISIGN*FACTOR IF(MOD((IVJ-1)/2,2)) 87,10,87 87 COEFP = -COEFP GO TO 10 7 COEFP = ITAB3(J2,J1) IF(COEFP) 65,10,88 65 COEFP = -SQRT(-COEFP/NORM3(J2)) GO TO 89 88 COEFP = SQRT(COEFP/NORM3(J2)) 89 COEFP = COEFP * ISIGN * FACTOR GO TO 10 8 COEFP = ITAB2(J2,J1) IF(COEFP) 66,10,90 66 COEFP = -SQRT(-COEFP/NORM2(J2)) GO TO 91 90 COEFP = SQRT(COEFP/NORM2(J2)) 91 COEFP = COEFP * ISIGN * FACTOR GO TO 10 9 COEFP = ITAB1(J2,J1) IF(COEFP) 67,10,92 67 COEFP = -SQRT(-COEFP/NORM1(J2)) GO TO 93 92 COEFP = SQRT(COEFP/NORM1(J2)) 93 COEFP = COEFP * ISIGN * FACTOR GO TO 10 C C AN UNALLOWED STATE C FOR AN UNALLOWED STATE THE F.P. COEFFICIENT IS SET EQUAL TO AN C ERRONEOUS VALUE.BY REPLACING THE 11 COEFP=9.9 STATEMENT BY THE C FOLLWING 3 CARDS THE USER CAN TERMINATE THE RUN WHEN AN C UNALLOWED STATE OCCURS C 106 FORMAT(37H FAIL IN CFPD AT 11 UNALLOWED STATE) C 11 WRITE(IWRITE,106) C PAUSE C 11 COEFP=9.9 10 RETURN END C C*********************************************************************** C C INSERT THE FOLLOWING SEGMENTS OF THE PROGRAM TO EVALUATE THE C REDUCED MATRIX ELEMENT OF SUMMATIONS OF ONE PARTICLE TENSOR C OPERATORS BY W D ROBB CATALOGUE NUMBER AAKF (CPC 6 (1973) 132) C AND CORRECTION DECK CATALOGUE NUMBER AAKF000A TO BE PUBLISHED C C TENSOR,SETUPE C C CARDS AAKF0113-AAKF0550,AAKF0645-AAKF0832 C C*********************************************************************** C C*********************************************************************** C C +++ REPLACE CARDS 132,133 BY THE NEXT TWO CARDS WITHOUT C C DIMENSION L6(100),L7(200),L8(100),LW(6,20),J2STO(10,3), C 1 J3STO(10,3),JMEM(5),VSHELL(20) C +++ REPLACE CARD 137 BY A COPY OF THE NEXT CARD WITHOUT THE C C COMMON/COUPLE/MN1,M0,J1(100),J2(&L74,3),J3(&L74,3),FREE(100) C +++ A DIMENSION INCREASE WAS NECESSARY C C +++ REPLACE CARD 141 BY A COPY OF THE NEXT CARD WITH THE C C COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 C +++ A DIFFERENT DEBUG PRINT BLOCK IS USED C C +++ INSERT A COPY OF THE NEXT CARD WITHOUT THE C BEFORE CARD 163 C NBUG6=IBUG6 C +++ TO SET UP DEBUG PRINT PARAMETER C C +++ REPLACE CARD 480 BY A COPY OF THE NEXT CARD WITHOUT THE C C 172 IF(ABS(RECUPS).LT.1.0E-5)GO TO 183 C C +++ REPLACE CARD 647 BY A COPY OF THE NEXT CARD WITHOUT C C DIMENSION NJCOMP(23),LJCOMP(23) C +++ A DIMENSION DECREASE WAS NECESSARY C C +++ REPLACE CARDS 649-650 BY THE NEXT THREE CARDS WOTHOUT C C COMMON/STATE/ C 1 MCFG,MOCCSH(40),MOCORB(15,40),MELCSH(15,40),M1QNRD(29,3,40), C 2 KCFG,KOCCSH(40),KOCORB(15,40),KELCSH(15,40),K1QNRD(29,3,40),MAXOR C +++ A DIFFERENT VERSION OF MSTATE IS USED C C*********************************************************************** C REAL FUNCTION CG(J1,J2,J3,M1,M2,M3) C C CALCULATES A CLEBSCH-GORDAN COEFFICIENT. C C J1,J2,J3 ARE THE J VALUES OF THE COEFFICIENT C M1,M2,M3 ARE THE M VALUES OF THE COEFFICIENT C C FORMULA FOR VECTOR COUPLING COEFFICIENT C (CONDON AND SHORTLEY PAGE 75 FORMULA (3.14.5)) C C IMPLICIT REAL*8(A-H,O-Z) COMMON /FACT/ GAMM( 250) C C = REAL(0) IF (ABS(M1).GT.J1) GO TO 9 IF (ABS(M2).GT.J2) GO TO 9 IF (ABS(M3).GT.J3) GO TO 9 IF (M1+M2.NE.M3) GO TO 9 J=MAX(J1,J2,J3) IF (2*J.GT.J1+J2+J3) GO TO 9 C IA=-J1+M1 IB=J2-J1+M3 IC=J2+J3+M1 ID=-J1+J2+J3 IE=J3+M3 NUMIN=MAX(0,IB)+1 NUMAX=MIN(IC,ID,IE)+1 C IF(NUMAX.LT.NUMIN) GO TO 9 DO 8 N=NUMIN,NUMAX 8 C=(GAMM(N-IA)*GAMM(IC-N+2))/(GAMM(ID-N+2)*GAMM(IE-N+2) * * GAMM(N-IB)*GAMM(N)) * (MOD(J2+M2+N,2)*2-1) + C C C = SQRT((GAMM(-J1+J2+J3+1)/GAMM(J2+M2+1))* 1 (GAMM(J1-J2+J3+1)/GAMM(J1+M1+1))* 2 (GAMM(J1+J2-J3+1)/GAMM(J1-M1+1))* 3 (GAMM(J3-M3+1)/GAMM(J2-M2+1))* 4 ((2*J3+1)*GAMM(J3+M3+1)/GAMM(J1+J2+J3+2))) * C 9 CG=C RETURN END C*********************************************************************** SUBROUTINE BAKSUB(N,X,B,U,V,W) C C SOLVES FOR X THE VECTOR EQUATION UPP X=B C WHERE UPP IS AN UPPER TRIANGULAR MATRIX WITH ONLY THREE NON-ZERO C DIAGONALS, X AND B ARE COLUMN VECTORS. THE ARRAYS U, V AND W ARE C AS DESCRIBED IN SUBROUTINE VECTOR. C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N), B(N),U(N),V(N),W(N) C X(N)=B(N)/U(N) X(N-1)=(B(N-1)-V(N-1)*X(N))/U(N-1) IF(N.EQ.2) GO TO 2 N1=N-1 DO 1 I=2,N1 1 X(N-I)=(B(N-I)-V(N-I)*X(N-I+1)-W(N-I)*X(N-I+2))/U(N-I) C 2 RETURN END C*********************************************************************** SUBROUTINE EIGEN(N,EIG,EPSI,P,R,POLY,BETA) C C ACCEPTS THE ARRAYS R AND P OF MAIN AND SUPER DIAGONAL ELEMENTS C RESPECTIVELY. USING THE STURM SEQUENCE PROPERTY C A BISECTION METHOD IS APPLIED TO DETERMINE THE EIGENVALUES C (STORED IN THE ARRAY EIG ON RETURN) TO AN ACCURACY C SPECIFIED BY EPSI. N IS AS DEFINED IN THE SUBROUTINE HSLDR. C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION EIG(N), P(N),R(N),POLY(N),BETA(N) DATA HALF/0.5/, TINY/1.0E-5/ C C CALCULATE THE AVERAGE OF THE GREATEST AND SMALLEST MAIN C DIAGONAL ELEMENTS STORING THE RESULT IN AMID. C ASMALL=R(1) ALARG=R(1) DO 1 I=2,N IF(R(I).GT.ALARG) ALARG=R(I) IF(R(I).LT.ASMALL) ASMALL=R(I) 1 CONTINUE AMID=(ALARG+ASMALL)*HALF C C REDUCE EACH MAIN DIAGONAL ELEMENT BY AMID AND CALCUATE, USING C THE STURM SEQUENCE PROPERTY, THE EIGENVALUES OF THE CORRESPONDING C REDUCED TRI-DIAGONAL MATRIX. C DO 2 I=1,N 2 R(I)=R(I)-AMID C C CALCULATE THE MAXIMUM INFINITY NORM G OF THE MATRIX. THE C EIGENVALUES LIE IN THE RANGE -G TO +G. C G=0.0 DO 3 I=1,N G1=ABS(R(I)) IF(I.GT.1)G1=G1+ABS(P(I-1)) IF(I.LT.N)G1=G1+ABS(P(I)) IF(G.LT.G1)G=G1 3 CONTINUE C C CALCULATE THE SQUARES OF THE SUPER DIAGONAL ELEMENTS AND STORE C THESE IN THE ARRAY BETA. C N1=N-1 DO 4 I=1,N1 4 BETA(I)=P(I)*P(I) C C THIS LOOP DETERMINES THE EIGENVALUES ONE AT A TIME IN ORDER OF C ALGEBRAIC SIZE DOWNWARDS. C DO 11 K=1,N AL=-G BL=G C C ONCE THROUGH THIS LOOP IS ONE BISECTION OF THE RANGE. CL1 IS THE C CURRENT ESTIMATE, CL THE IMMEDIATELY PREVIOUS ESTIMATE OF THE C EIGENVALUE. C DO 9 J=1,100 CL1=(AL+BL)*HALF IF(J.EQ.1) GO TO 5 C C IF THE EIGENVALUE HAS BEEN DETERMINED TO A SPECIFIED ACCURACY C EPSI, THE CALCULATION IS COMPLETE. C IF(ABS(CL1-CL).LT.EPSI) GO TO 10 C C LSUM STORES THE NUMBER OF AGREEMENTS IN SIGN IN THE STURM C SEQUENCE. C 5 LSUM=0 DO 51 I=1,N 51 POLY(I)=R(I)-CL1 X=POLY(1) IF(POLY(1).GT.0.0) LSUM=1 C C THIS LOOP CALCULATES ALL THE REMAINING MEMBERS OF THE STURM C SEQUENCE. THE NUMBER OF AGREEMENTS IN SIGN IS ALSO DETERMINED. C DO 6 I=2,N IF(X.EQ.0.0) THEN X=POLY(I)-ABS(P(I-1))/(1.0+TINY) ELSE X=POLY(I)-BETA(I-1)/X ENDIF IF(X.GT.0.0) LSUM=LSUM+1 6 CONTINUE C C THE NEW RANGE FOR THE EIGENVALUE (DEPENDENT ON THE VALUE OF LSUM) C IS DETERMINED. C CL=CL1 IF(LSUM.GE.K) GO TO 8 BL=CL1 GO TO 9 8 AL=CL1 9 CONTINUE C C THE EIGENVALUE IS STORED IN THE ARRAY EIG C 10 EIG(K)=CL1 C C RETURN TO CALCULATE THE NEXT EIGENVALUE. C 11 CONTINUE C C THE ELEMENTS OF THE ORIGINAL TRI-DIAGONAL MATRIX ARE REGAINED C AND ITS EIGENVALUES FOUND. C DO 12 I=1,N R(I)=R(I)+AMID 12 EIG(I)=EIG(I)+AMID C RETURN END C*********************************************************************** SUBROUTINE EIGVEC(N,A,LENGTH,X,P) C C TAKES THE EIGENVECTOR OF THE TRI-DIAGONAL MATRIX C STORED IN X AND DETAILS OF THE MATRICES USED IN TRANSFORMING C THE ORIGINAL MATRIX TO TRI-DIAGONAL FORM, STORED IN A, C AND OBTAINS THE CORRESPONDING EIGENVECTOR OF THE ORIGINAL C MATRIX. N AND LENGTH ARE AS DEFINED IN THE SUBROUTINE HSLDR. C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(LENGTH),X(N), P(N) C N2=N-2 N22=(N*(N+1))/2 C C N2=N-2 TRANSFORMATIONS TO OBTAIN EACH EIGENVECTOR. C DO 3 K=1,N2 C C K1 IS THE NUMBER OF ELEMENTS IN THE FIRST (K-1) ROWS OF THE C UPPER TRIANGLE STORED IN A. C K1=N22-((K+2)*(K+3))/2+1 SOP=0.0 KP1=K+1 NK1=N-K-1 DO 1 I=1,KP1 1 SOP=SOP+A(K1+I)*X(NK1+I) IF(SOP.EQ.0.0) GO TO 3 C C FROM INFORMATION STORED IN THE ARRAY A BKR IS DETERMINED AS IN C SUBROUTINE HOUSE. C BKR=-(P(NK1)*A(K1+1)) SOP=SOP/BKR DO 2 J=1,KP1 2 X(NK1+J)=X(NK1+J)-A(K1+J)*SOP 3 CONTINUE C C THE EIGENVECTOR OF THE ORIGINAL MATRIX IS NORMALISED. C CALL NORM(N,X) RETURN END C*********************************************************************** SUBROUTINE HOUSE(N,A,LENGTH,P,R,ARRAY1) C C ACCEPTS THE UPPER TRIANGLE OF ELEMENTS OF A SYMMETRIC MATRIX, C STORED IN THE LINEAR ARRAY A, AND USING THE HOUSEHOLDER METHOD C REDUCES THIS TO TRI-DIAGONAL FORM, STORING THE NEW MAIN C DIAGONAL ELEMENTS IN POSITION IN A AND ALSO IN THE ARRAY R, C AND THE SUPER-DIAGONAL ELEMENTS IN THE ARRAY P. DETAILS C OF THE TRANSFORMING MATRICES ARE OVERWRITTEN IN THE REDUNDANT C SPACE OF A. N AND LENGTH ARE AS DEFINED IN SUBROUTINE HSLDR. C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(LENGTH), P(N),R(N),ARRAY1(N) DATA EPSI/1.0E-9/, HALF/0.5/ C N1=N-1 N2=N-2 C C EACH TIME THROUGH THIS LOOP ONE MORE ROW OF THE MATRIX IS C TRANSFORMED TO TRI-DIAGONAL FORM C DO 11 K=1,N2 DO 1 J=K,N 1 P(J)=0.0 C C KB=KC-1 IS THE NUMBER OF ELEMENTS IN THE FIRST (K-1) ROWS OF THE C UPPER TRIANGLE C NJ=N-K KC = ((K-1)*(NJ+N+2))/2 + 1 C C THE SQUARE ROOT OF THE SUM OF THE SQUARES OF THE REMAINING C OFF DIAGONAL ELEMENTS IN ROW K IS FOUND AND STORED IN SUM C SUM=0.0 DO 2 J=1,NJ 2 SUM=A(KC+J)*A(KC+J)+SUM SUM=SQRT(SUM) IF(SUM.LT.EPSI) GO TO 100 C C SUM IS GIVEN THE SAME SIGN AS THE SUPER DIAGONAL ELEMENT IN ROW K C IF(A(KC+1).LT.0.0) SUM=-SUM BKR=SUM*(SUM+A(KC+1)) C C THE FIRST ELEMENT OF THE VECTOR FROM WHICH THE TRANSFORMING C MATRIX IS DERIVED IS OVERWRITTEN ON THE OLD SUPER DIAGONAL C ELEMENT IN ROW K. THE REMAINING ELEMENTS ARE ALREADY IN POSITION C IN ROW K. C A(KC+1)=A(KC+1)+SUM C C THE SUPER DIAGONAL ELEMENT IN ROW K OF THE NEW TRI-DIAGONAL C MATRIX IS STORED IN ARRAY P C P(K)=-SUM C C KD IS THE NUMBER OF ELEMENTS IN THE FIRST K ROWS OF THE UPPER C TRIANGLE C KD=(K*(NJ+N+1))/2 C C THE TRANSFORMATION DERIVED FROM THE NJ=(N-K) VECTOR ELEMENTS C STORED IN A IS NOW APPLIED TO THE LAST NJ ROWS OF THE MATRIX C (THE LAST NJ SPACES OF THE ARRAY P ARE SUCCESSIVELY C OVERWRITTEN IN THE PROCESS). C LOL=KC-K+1 DO 6 M=K,N1 NM=N-M MO=((NM+N+1)*M)/2 M1=M-1 LO=LOL+M1 SUM=P(1+M) DO 3 L=1,NM 3 SUM=SUM+A(MO+L)*A(LO+L) IF(M.EQ.K) GO TO 6 JM1=KD+M-K+1 DO 4 L=K,M1 ARRAY1(L)=A(JM1) 4 JM1=JM1+N1-L DO 5 L=K,M1 5 SUM=SUM+ARRAY1(L)*A(LOL+L) 6 P(1+M)=SUM/BKR C SUM=0.0 DO 7 L=1,NJ 7 SUM=SUM+A(KC+L)*P(K+L) SUM=SUM/BKR DO 8 L=1,NJ 8 P(K+L)=P(K+L)-A(KC+L)*SUM*HALF DO 10 I=1,NJ AIM=A(KC+I) PIN=P(K+I) DO 9 J=I,NJ 9 ARRAY1(J)=-AIM*P(K+J)-PIN*A(KC+J) KE=KD+(I-1)*NJ-(I*(I-1))/2 DO 10 J=I,NJ 10 A(KE+J)=A(KE+J)+ARRAY1(J) GO TO 11 C 100 P(K)=-EPSI DO 101 I=1,NJ 101 A(I+KC)=0.0 11 CONTINUE C C THE LAST SUPER DIAGONAL ELEMENT IS ENTERED INTO THE ARRAY P. C ILK=((N+1)*N)/2-1 P(N1)=A(ILK) C C THE MAIN DIAGONAL ELEMENTS PICKED OUT FROM THE ARRAY A ARE C STORED IN THE ARRAY R. C DO 12 I=1,N ILK=((I-1)*(2*N-I+2))/2+1 12 R(I)=A(ILK) C RETURN END C*********************************************************************** SUBROUTINE HSLDR(N,A,LENGTH,EPSI,EIG,X,NO,P,MSV) C C THIS SUBROUTINE ACCEPTS THE UPPER TRIANGLE OF AN N*N SYMMETRIC C MATRIX AND ON THE FIRST CALL DETERMINES ALL THE EIGENVALUES AND C THE FIRST EIGENVECTOR. ON EACH FURTHER CALL ONE MORE OF THE C REMAINING EIGENVECTORS IS CALCULATED. C C DEFINITION OF THE ARGUMENTS. C N........... THE DEGREE OF THE SYMMETRIC MATRIX TO BE C DIAGONALISED. C A........... THE LINEAR ARRAY CONTAINING THE UPPER TRIANGLE OF C THE ORIGINAL MATRIX, OVERWRITTEN ON RETURN BY THE C MAIN DIAGONAL ELEMENTS OF THE TRI-DIAGONAL MATRIX C AND DETAILS OF THE TRANSFORMING MATRICES. C LENGTH...... =(N*(N+1))/2 , THE SIZE OF THE ARRAY A. C EPSI........ THE ACCURACY TO WHICH THE EIGENVALUES ARE TO BE C DETERMINED. C EIG......... THIS ARRAY CONTAINS THE EIGENVALUES ON RETURN. C X........... THIS ARRAY CONTAINS ONE EIGENVECTOR ON RETURN. C NO.......... THIS RUNS FROM 1 TO N AND SPECIFIES WHICH C EIGENVECTOR IS STORED IN X ON RETURN C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(LENGTH),EIG(N),X(N), P(MSV,9) CC COMMON/SAVE1/ P(900),R(900),ARRAY1(900) C C THE SIZE OF THE MATRIX WAS LIMITED BY THE SIZE OF THE INTERNAL C ARRAYS IN THE SAVE1, SAVE2 AND SAVE3 COMMON BLOCKS. C IF(N-2) 1,2,3 1 X(1)=1.0 EIG(1)=A(1) P(1,2)=A(1) GO TO 9 2 P(1,1)=A(2) P(1,2)=A(1) P(2,2)=A(3) GO TO 4 C C IF THE FIRST EIGENVECTOR HAS ALREADY BEEN FOUND THE C TRI-DIAGONALISING AND EIGENVALUE SUBROUTINES ARE SKIPPED ROUND. C 3 IF(NO.NE.1) GO TO 5 C C THE TRI-DIAGONALISING SUBROUTINE IS ENTERED. C CALL HOUSE(N,A,LENGTH,P(1,1),P(1,2),P(1,9)) C C THE ELEMENTS OF THE TRI-DIAGONAL MATRIX ARE USED TO DETERMINE C THE EIGENVALUES. C 4 CALL EIGEN(N,EIG,EPSI,P(1,1),P(1,2),P(1,3),P(1,4)) C C THE EIGENVECTOR OF THE TRI-DIAGONAL MATRIX CORRESPONDING TO A C PARTICULAR EIGENVALUE IS DETERMINED C 5 CALL VECTOR(N,EIG,X,NO,P(1,1),P(1,2),P(1,3), * P(1,4),P(1,5),P(1,6),P(1,7),P(1,8),P(1,9)) IF(N.EQ.2) GO TO 6 C C THE CORRESPONDING EIGENVECTOR OF THE ORIGINAL MATRIX IS FOUND C CALL EIGVEC(N,A,LENGTH,X,P) C C NORMALIZING THE EIGENVECTOR C 6 APP=0.0 DO 7 I=1,N 7 APP=APP+X(I)*X(I) APP=SQRT(APP) DO 8 I=1,N 8 X(I)=X(I)/APP C 9 RETURN END C*********************************************************************** SUBROUTINE NORM(N,X) C C NORMALISES THE VECTOR X OF DIMENSION N C SUCH THAT THE LARGEST COMPONENT IS UNITY C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N) C G=ABS(X(1)) DO 1 I=2,N GA=ABS(X(I)) IF(GA.GT.G)G=GA 1 CONTINUE DO 2 I=1,N 2 X(I)=X(I)/G C RETURN END C*********************************************************************** SUBROUTINE VECTOR(N,EIG,X,NO,P,R,RN,XM,LXCH,B,U,V,W) C C TAKES ARRAYS R OF MAIN DIAGONAL ELEMENTS, P OF SUPER DIAGONAL C ELEMENTS, EIG OF EIGENVALUES, OF THE TRI-DIAGONAL MATRIX C AND BY MEANS OF INVERSE ITERATIONS DETERMINES C AN EIGENVECTOR OF THE TRI-DIAGONAL MATRIX. C C IMPLICIT REAL*8(A-H,O-Z) LOGICAL ITER DIMENSION EIG(N),X(N),P(N),R(N), * RN(N),XM(N),LXCH(N),B(N),U(N),V(N),W(N) DATA EPSI/1.0E-9/ C C THE ARRAY RN HOLDS THE MAIN DIAGONAL ELEMENTS OF A NEW C TRI-DIAGONAL MATRIX. C DO 1 K=1,N 1 RN(K)=R(K)-EIG(NO) C C BY MEANS OF GAUSSIAN ELIMINATION THE NEW TRI-DIAGONAL MATRIX C IS TRANSFORMED INTO UPPER TRIANGULAR FORM. THE ROW MULTIPLIERS C ARE STORED IN ARRAY XM. IF ROW I IS INTERCHANGED WITH ROW I+1 C LXCH(I)=1 THE MAIN DIAGONAL ELEMENTS OF THE UPPER TRIANGULAR C MATRIX ARE STORED IN THE ARRAY U, THE NEXT DIAGONAL IN THE C ARRAY V AND THE LAST DIAGONAL IN THE ARRAY W. C PA=RN(1) QA=P(1) N1=N-1 DO 4 I=1,N1 C C DETERMINE IF A ROW INTERCHANGE IS NECESSARY C GA=ABS(P(I)) PPA=ABS(PA) IF(GA.LE.EPSI) GO TO 20 IF(GA.GT.PPA) GO TO 2 C C NO INTERCHANGE C 20 U(I)=PA V(I)=QA W(I)=0.0 HA=P(I) PA=RN(I+1) QA=P(I+1) LXCH(I)=0 GO TO 3 C C INTERCHANGE C 2 U(I)=P(I) V(I)=RN(I+1) W(I)=P(I+1) HA=PA PA=QA QA=0.0 LXCH(I)=1 C C THE ROW MULTIPLIER IS DETERMINED C 3 XM(I)=HA/U(I) C C ROW I IS MULTIPLIED BY XM(I) AND SUBTRACTED FROM ROW I+1 C HA=0.0 PA=PA-XM(I)*V(I) 4 QA=QA-XM(I)*W(I) C C THE SINGLE ELEMENT IN THE LAST ROW IS PLACED IN U(N) C U(N)=PA IF(ABS(U(N)).LT.EPSI) U(N)=EPSI C C THE ARRAY X IS RESERVED FOR THE CURRENT ESTIMATE OF THE C EIGENVECTOR. THE SUBROUTINE BAKSUB IS CALLED AND AN ESTIMATE OF C THE EIGENVECTOR OBTAINED. C 5 DO 6 I=1,N 6 B(I)=1.0 CALL BAKSUB(N,X,B,U,V,W) C C THE EIGENVECTOR STORED IN X IS NORMALISED C 7 CALL NORM(N,X) C C TO SAVE COMPUTING TIME, A FURTHER ITERATION TO YIELD A BETTER C ESTIMATE OF THE VECTOR CAN BE SUPPRESSED BY SETTING ITER=.TRUE. C 8 ITER=.FALSE. IF(ITER) GO TO 13 C C A NEW COLUMN VECTOR RELATED TO X BY EXACTLY THE SAME ROW C INTERCHANGES AND MULTIPLICATIONS BY WHICH THE UPPER TRIANGULAR C MATRIX WAS OBTAINED FROM THE NEW TRI-DIAGONAL MATRIX, IS C CALCULATED AND STORED IN B. C DO 9 I=1,N 9 B(I)=X(I) N1=N-1 DO 10 I=1,N1 IF(LXCH(I).EQ.0) GO TO 10 AC=B(I) BC=B(I+1) B(I)=BC B(I+1)=AC 10 B(I+1)=B(I+1)-XM(I)*B(I) C C THE NEW VECTOR STORED IN B IS NORMALISED C CALL NORM(N,B) C C THE SUBROUTINE BAKSUB IS CALLED AND A NEW ESTIMATE OF THE C EIGENVECTOR IS OBTAINED C CALL BAKSUB(N,X,B,U,V,W) C C THE NEW ESTIMATE OF THE EIGENVECTOR IS NORMALISED C CALL NORM(N,X) 13 RETURN END C*********************************************************************** SUBROUTINE FANO(IRHO,ISIG,IRHOP,ISIGP) C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C PARAMETER(LL43= 21*2+3, LL74= 21+4, LL75= 21+2, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20, + KFLS=12,KFLN=10,KFLV=40) C DIMENSION RMEDIR(99),RMEEX(99),NBAR(LL75) ! '05Feb20 from NBAR(21) C LOGICAL FAILSD,FAILSE,FAILAD,FAILAE,FREE C LOGICAL SMVRSD,SMVRSE,SMVRAD,SMVRAE DIMENSION K6SD(KFL6),K7SD(KFL7),K8SD(KFL8),K9SD(KFL9),KWSD(6,KFLW) + ,LDELSD(KFLW,2),SMVRSD(KFL1) DIMENSION K6SE(KFL6),K7SE(KFL7),K8SE(KFL8),K9SE(KFL9),KWSE(6,KFLW) + ,LDELSE(KFLW,2),SMVRSE(KFL1) DIMENSION K6AD(KFL6),K7AD(KFL7),K8AD(KFL8),K9AD(KFL9),KWAD(6,KFLW) + ,LDELAD(KFLW,2),SMVRAD(KFL1) DIMENSION K6AE(KFL6),K7AE(KFL7),K8AE(KFL8),K9AE(KFL9),KWAE(6,KFLW) + ,LDELAE(KFLW,2),SMVRAE(KFL1) C DIMENSION J6PSD(KFLV),J7PSD(KFLV),J8PSD(KFLV),J9PSD(KFLV), + JWRDSD(6,KFLW), + NBJSD(KFLN),NB6JSD(KFLN),K6CPSD(KFLN),K7CPSD(KFLN),K8CPSD(KFLN), + K9CPSD(KFLN),JSM6SD(KFLS),JSM4SD(KFLS,KFLW),JSM5SD(KFLS,KFLW), + IN6JSD(KFLW) DIMENSION J6PSE(KFLV),J7PSE(KFLV),J8PSE(KFLV),J9PSE(KFLV), + JWRDSE(6,KFLW), + NBJSE(KFLN),NB6JSE(KFLN),K6CPSE(KFLN),K7CPSE(KFLN),K8CPSE(KFLN), + K9CPSE(KFLN),JSM6SE(KFLS),JSM4SE(KFLS,KFLW),JSM5SE(KFLS,KFLW), + IN6JSE(KFLW) DIMENSION J6PAD(KFLV),J7PAD(KFLV),J8PAD(KFLV),J9PAD(KFLV), + JWRDAD(6,KFLW), + NBJAD(KFLN),NB6JAD(KFLN),K6CPAD(KFLN),K7CPAD(KFLN),K8CPAD(KFLN), + K9CPAD(KFLN),JSM6AD(KFLS),JSM4AD(KFLS,KFLW),JSM5AD(KFLS,KFLW), + IN6JAD(KFLW) DIMENSION J6PAE(KFLV),J7PAE(KFLV),J8PAE(KFLV),J9PAE(KFLV), + JWRDAE(6,KFLW), + NBJAE(KFLN),NB6JAE(KFLN),K6CPAE(KFLN),K7CPAE(KFLN),K8CPAE(KFLN), + K9CPAE(KFLN),JSM6AE(KFLS),JSM4AE(KFLS,KFLW),JSM5AE(KFLS,KFLW), + IN6JAE(KFLW) C COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 COMMON/KRON/IDEL(LL75,LL75) COMMON/TERMS/ NROWS,ITAB(18),JTAB(18),NTAB(189) COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75), 1 J1QN1(LL43,3),J1QN2(LL43,3),IJFUL(LL75) COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15, 1 M16,M17,M18,M19,M20 COMMON/XATION/ AMULT(99),BMULT(99),KD1,KD2,KE1,KE2,MULTD,MULTE COMMON/NJLJ/NRHO,LRHO,NSIG,LSIG,NRHOP,LRHOP,NSIGP,LSIGP COMMON/INTERM/J1BAR1(LL75,3),J1BAR2(LL75,3),J1TLD1(3),J1TLD2(3) COMMON/COUPLE/NJ1S,NJ23S,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) C LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP C COMMON/SUMARG/J6P(KFLV),J7P(KFLV),J8P(KFLV),J9P(KFLV), + JWORD(6,KFLW),NLSUM, + NBJ(KFLN),NB6J(KFLN),K6CP(KFLN),K7CP(KFLN),K8CP(KFLN),K9CP(KFLN), + JSUM6(KFLS),JSUM4(KFLS,KFLW),JSUM5(KFLS,KFLW),INV6J(KFLW) C C C === DETERMINES POTENTIAL MATRIX ELEMENT FOR GIVEN RHO,SIG,RHOP,SIGP C C 301 FORMAT(21H NO POSSIBLE K-VALUES) 302 FORMAT(66H SPECTATOR QUANTUM NUMBERS NOT DIAGONAL FOR NON-INTERACT 1ING SHELLS) 305 FORMAT(23H DIRECT SPIN INTEGRAL =,F10.6) 306 FORMAT(25H EXCHANGE SPIN INTEGRAL =,F10.6) 307 FORMAT(6H NJ,LJ,4(I8,I4)) 308 FORMAT(6H KD1 =,I4,6H KD2 =,I4,6H KE1 =,I4,6H KE2 =,I4,6X,'IHSH =' *,I2) 309 FORMAT(56H ROWS OF TERM TABLE CONTAINING PARENTS ARE RESPECTIVELY 1,2(I6,I3)) 310 FORMAT(26H DIRECT ANGULAR INTEGRAL =,F10.6) 312 FORMAT(28H EXCHANGE ANGULAR INTEGRAL =,F10.6) C C --- SET USEFUL CONSTANTS C C TST IBUG1=ABS(JBUG1) M1=ISIG-IRHO M2=ISIGP-IRHOP M19=IRHO-IRHOP M20=ISIG-ISIGP MULTD=0 MULTE=0 C C JSNDIR,JANGDI=0 IMPLIES APPROPRIATE J2,J3 ARRAYS FOR CALL OF C NJGRAF HAVE NOT BEEN SET C JSNDIR=0 JANGDI=0 ISPIND=0 ISPINE=0 IANGD=0 IANGE=0 C C SET THE FAIL PARAMETERS .FALSE. INITIALLY SO THAT NJGRAF CAN C BE CALLED. C FAILSD=.FALSE. FAILSE=.FALSE. FAILAD=.FALSE. FAILAE=.FALSE. C C --- SET N,L VALUES OF INTERACTING SHELLS C NRHO=NJ(IRHO) LRHO=LJ(IRHO) NSIG=NJ(ISIG) LSIG=LJ(ISIG) NRHOP=NJ(IRHOP) LRHOP=LJ(IRHOP) NSIGP=NJ(ISIGP) LSIGP=LJ(ISIGP) IF(IBUG1.LE.1) GO TO 160 WRITE(IWRITE,307) NRHO,LRHO,NSIG,LSIG,NRHOP,LRHOP,NSIGP,LSIGP C C --- EVALUATE MULTIPLICATIVE FACTORS C 160 IL=IDEL(IRHO,ISIG) IR=IDEL(IRHOP,ISIGP) Q=NOSH1(IRHO)*(NOSH1(ISIG)-IL)*NOSH2(IRHOP)*(NOSH2(ISIGP)-IR) XMULT=SQRT(Q)*0.5 ADIRCT=((IL-1)*(IR-1)+1)/SQRT(REAL((LSIG+LSIG+1)*(LRHOP+LRHOP+1))) IEXCHG=2-IL-IR AEXCHG=IEXCHG/SQRT(REAL((LSIG+LSIG+1)*(LSIGP+LSIGP+1))) DO 13 J=1,IHSH 13 NBAR(J)=NOSH1(J)-IDEL(J,IRHO)-IDEL(J,ISIG) IDELP=0 IF(M1.EQ.0) GO TO 15 K1=IRHO+1 DO 16 J=K1,ISIG 16 IDELP=IDELP+NBAR(J) 15 IF(M2.EQ.0) GO TO 18 K1=IRHOP+1 DO 19 J=K1,ISIGP 19 IDELP=IDELP+NBAR(J) 18 XMULT = (1-MOD(IDELP,2)*2)*XMULT C C --- DETERMINES RANGES OF K FOR DIRECT AND EXCHANGE INTEGRALS C TRIANGULAR RELATIONS LIMIT (K+1) VALUES TO LIE BETWEEN KD1 AND KD2 C FOR =DIRECT= INTEGRALS, KE1 AND KE2 FOR =EXCHANGE= INTEGRALS C K1=ABS(LSIG-LSIGP) K2=LSIG+LSIGP K3=ABS(LRHO-LRHOP) K4=LRHO+LRHOP KD1=MAX(K1,K3)+1 KD2=MIN(K2,K4)+1 K1=ABS(LRHOP-LSIG) K2=LRHOP+LSIG K3=ABS(LRHO-LSIGP) K4=LRHO+LSIGP KE1=MAX(K1,K3)+1 KE2=MIN(K2,K4)+1 IF(IBUG1.GT.1) WRITE(IWRITE,308) KD1,KD2,KE1,KE2, IHSH IF(KD1.LE.KD2) GO TO 213 IF(KE1.LE.KE2) GO TO 213 IF(IBUG1.GT.1) WRITE(IWRITE,301) C C --- ZEROIZE MULTIPLYING FACTORS FOR ALLOWED K-VALUES. THE LOWEST C VALUES KD1 AND KD2 ARE ALWAYS ALLOWED (UNLESS THEY ARE C GREATER THAN KD2 AND KE2 RESPECTIVELY). ALLOWED K-VALUES THEN C STEP BY 2 TO SATISFY THE EVEN CONDITION OF THE REDUCED MATRIX C ELEMENTS, WHICH ARE THEN CALCULATED AND STORED C 213 IF(KD1.GT.KD2) GO TO 232 IF(KD2.GE.99) GO TO 999 DO 230 K=KD1,KD2,2 AMULT(K)=0.0 230 RMEDIR(K)=RME(LRHO,LRHOP,K-1)*RME(LSIG,LSIGP,K-1) 232 IF(KE1.GT.KE2) GO TO 27 IF(KE2.GE.99) GO TO 999 DO 234 K=KE1,KE2,2 BMULT(K)=0.0 234 RMEEX(K)=RME(LRHO,LSIGP,K-1)*RME(LSIG,LRHOP,K-1) C C --- SET SENIORITY, S AND L VALUES OF SPECTATOR SHELLS C 27 DO 26 J=1,IHSH IF(IRHO.EQ.J) GO TO 29 IF(ISIG.EQ.J) GO TO 29 DO 28 K=1,3 28 J1BAR1(J,K)=J1QN1(J,K) 29 IF(IRHOP.EQ.J) GO TO 26 IF(ISIGP.EQ.J) GO TO 26 DO 30 K=1,3 30 J1BAR2(J,K)=J1QN2(J,K) IF(IRHO.EQ.J) GO TO 26 IF(ISIG.EQ.J) GO TO 26 C CHECK THAT SPECTATOR SHELLS HAVE SAME RESPECTIVE QUANTUM NUMBERS DO 31 K=1,3 IF(J1BAR1(J,K).EQ.J1BAR2(J,K)) GO TO 31 IF(IBUG1.GT.1) WRITE(IWRITE,302) GO TO 527 31 CONTINUE 26 CONTINUE C C --- FROM WHICH ROWS OF NTAB DO WE FIND THE QUANTUM NUMBERS WITH BARS C OR TILDES C N=NOSH1(ISIG) K2=NTAB1(N,LSIG+1) N=NOSH1(IRHO) IF(M1.EQ.0) N=NOSH1(ISIG)-1 K1=NTAB1(N,LRHO+1) N=NOSH2(ISIGP) K4=NTAB1(N,LSIGP+1) N=NOSH2(IRHOP) IF(M2.EQ.0) N=NOSH2(ISIGP)-1 K3=NTAB1(N,LRHOP+1) IF(IBUG1.GT.1) WRITE(IWRITE,309) K1,K2,K3,K4 KK1=ITAB(K1) KK2=ITAB(K2) KK3=ITAB(K3) KK4=ITAB(K4) C C C === SUM OVER QUANTUM NUMBERS WITH BARS OR TILDES C DO 151 JJ2=1,KK2 C C --- TEST TO SEE WHICH PARENT TERMS ARE ALLOWABLE. ONLY TEST THIS ON C L AND S VALUES AT THIS STAGE, BY MEANS OF TRIANGULAR CONDITIONS C FOR TWICE THE QUANTUM NUMBERS (THUS ALWAYS FOR INTEGER QUANTITIES) C IN3=2*LSIG IJK2=3*(JJ2-1)+JTAB(K2) DO 140 K=2,3 IN1=NTAB(IJK2+K)-1 IN2=J1QN1(ISIG,K)-1 IF(IN1.GT.IN2+IN3) GO TO 151 IF(IN1.LT.ABS(IN2-IN3)) GO TO 151 140 IN3=1 DO 152 JJ1=1,KK1 IN3=2*LRHO IJK1=3*(JJ1-1)+JTAB(K1) 162 DO 145 K=2,3 IN1=NTAB(IJK1+K)-1 IF(M1.EQ.0) GO TO 142 IN2=J1QN1(IRHO,K)-1 GO TO 143 142 IN2=NTAB(IJK2+K)-1 143 IF(IN1.GT.IN2+IN3) GO TO 152 IF(IN1.LT.ABS(IN2-IN3)) GO TO 152 145 IN3=1 DO 153 JJ4=1,KK4 IN3=2*LSIGP IJK4=3*(JJ4-1)+JTAB(K4) DO 133 K=2,3 IN1=NTAB(IJK4+K)-1 IN2=J1QN2(ISIGP,K)-1 IF(IN1.GT.IN2+IN3) GO TO 153 IF(IN1.LT.ABS(IN2-IN3)) GO TO 153 133 IN3=1 DO 154 JJ3=1,KK3 IN3=2*LRHOP IJK3=3*(JJ3-1)+JTAB(K3) 137 DO 150 K=2,3 IN1=NTAB(IJK3+K)-1 IF(M2.EQ.0) GO TO 139 IN2=J1QN2(IRHOP,K)-1 GO TO 148 139 IN2=NTAB(IJK4+K)-1 148 IF(IN1.GT.IN2+IN3) GO TO 154 IF(IN1.LT.ABS(IN2-IN3)) GO TO 154 150 IN3=1 C C SUMMATIONS NOW PERFORMED OVER ALLOWED QUANTUM NUMBERS C THE TILDES CORRESPOND TO IRHO=ISIG AND/OR IRHOP=ISIGP C C --- SET THE REMAINING QUANTUM NUMBERS WITH BARS OR TILDES C DO 35 K=1,3 J1BAR1(IRHO,K)=NTAB(IJK1+K) J1BAR2(IRHOP,K)=NTAB(IJK3+K) J1TLD1(K)=NTAB(IJK2+K) IF(M1.NE.0) J1BAR1(ISIG,K)=NTAB(IJK2+K) J1TLD2(K)=NTAB(IJK4+K) IF(M2.NE.0) J1BAR2(ISIGP,K)=NTAB(IJK4+K) 35 CONTINUE C C --- IS POTENTIAL DIAG. IN BARRED QU. NOS. FOR INTERACTING SHELLS C I5=0 I=ISIG GO TO 50 42 I=IRHO IF(M1.NE.0) GO TO 50 I5=I5+1 45 I=ISIGP GO TO 50 46 I=IRHOP IF(M2.EQ.0) GO TO 48 50 I5=I5+1 DO 51 K=1,3 IF(J1BAR1(I,K).NE.J1BAR2(I,K)) GO TO 154 51 CONTINUE GO TO (42,45,46,48),I5 48 PICFP=1.0 C C --- EVALUATE FRACTIONAL PARENTAGE COEFFICIENTS C I=1 CALL MUMDAD (I,ISIG,IRHO,M1,CFPLHS) PICFP=PICFP*CFPLHS IF(ABS(PICFP).LT.EPS) GO TO 154 53 I=2 CALL MUMDAD(I,ISIGP,IRHOP,M2,CFPRHS) PICFP=PICFP*CFPRHS IF(ABS(PICFP).LT.EPS) GO TO 154 C C === SET UP J1,J2,J3 AND EVALUATE RECOUPLING COEFFICIENTS C C --- FIRST OF ALL, THE SPIN COEFFICIENTS C 55 I=3 CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,ISPIND,ISPINE,KK1,KK2,KK3,KK4) IF (ISPIND.NE.0) GO TO 90 CALL J23SPN(IRHO,ISIG,IRHOP,ISIGP,JSNDIR) ISPIND=1 C C --- DIRECT SPIN INTEGRAL C 90 SPINDT=0.0 IF(KD1.GT.KD2) GO TO 78 IF (FAILSD) GO TO 78 IF (ISPIND.NE.2) THEN CALL NJGRAF(SPINDT,FAILSD) ISPIND=2 IF(FAILSD) GO TO 78 CALL KNJ(J6CSD,J7CSD,J8CSD,J9CSD,JWCSD,K6SD,K7SD,K8SD,K9SD, + KWSD,JDELSD,LDELSD,SMVRSD,MPSD, + J6PSD,J7PSD,J8PSD,J9PSD,JWRDSD,NLSMSD,NBJSD,NB6JSD, + K6CPSD,K7CPSD,K8CPSD,K9CPSD,JSM4SD,JSM5SD,JSM6SD, + IN6JSD) ELSE CALL GENSUM(J6CSD,J7CSD,J8CSD,J9CSD,JWCSD,K6SD,K7SD,K8SD,K9SD, + KWSD,JDELSD,LDELSD,SMVRSD,MPSD, + J6PSD,J7PSD,J8PSD,J9PSD,JWRDSD,NLSMSD,NBJSD,NB6JSD, + K6CPSD,K7CPSD,K8CPSD,K9CPSD,JSM4SD,JSM5SD,JSM6SD, + IN6JSD,SPINDT) C ENDIF 78 IF(IBUG1.GT.1) WRITE(IWRITE,305) SPINDT C C IEXCHG IS ZERO WHENEVER M1=0=M2, IN WHICH CASE THE EXCHANGE C INTEGRAL HAS ZERO COEFFICIENT. THERE IS THEN NO POINT IN C CALCULATING THIS INTEGRAL, AND SPINEX IS SET ZERO (AT LABEL 93) C AS A MARKER OF THIS SITUATION C SPINEX=0.0 IF(IEXCHG.EQ.0) GO TO 94 C C --- MODIFY J2 AND J3 TO CALCULATE THE EXCHANGE SPIN INTEGRAL C IF (ISPINE.EQ.0) THEN I=1 CALL MODJ23(I) ISPINE=1 ENDIF C C --- EXCHANGE SPIN INTEGRAL C IF(KE1.GT.KE2) GO TO 94 IF (FAILSE) GO TO 94 IF (ISPINE.NE.2) THEN CALL NJGRAF(SPINEX,FAILSE) ISPINE=2 IF (FAILSE) GO TO 94 CALL KNJ(J6CSE,J7CSE,J8CSE,J9CSE,JWCSE,K6SE,K7SE,K8SE,K9SE, + KWSE,JDELSE,LDELSE,SMVRSE,MPSE, + J6PSE,J7PSE,J8PSE,J9PSE,JWRDSE,NLSMSE,NBJSE,NB6JSE, + K6CPSE,K7CPSE,K8CPSE,K9CPSE,JSM4SE,JSM5SE,JSM6SE, + IN6JSE) ELSE CALL GENSUM(J6CSE,J7CSE,J8CSE,J9CSE,JWCSE,K6SE,K7SE,K8SE,K9SE, + KWSE,JDELSE,LDELSE,SMVRSE,MPSE, + J6PSE,J7PSE,J8PSE,J9PSE,JWRDSE,NLSMSE,NBJSE,NB6JSE, + K6CPSE,K7CPSE,K8CPSE,K9CPSE,JSM4SE,JSM5SE,JSM6SE, + IN6JSE,SPINEX) ENDIF 94 IF(IBUG1.GT.1) WRITE(IWRITE,306) SPINEX C C --- MULTIPLY SPIN INTEGRALS BY PRODUCT OF FRACTIONAL PARENTAGE COEFFS C BDIRCT=SPINDT*PICFP BEXCHG=SPINEX*PICFP C C --- THE ANGULAR RECOUPLING COEFFICIENTS C SET J1,J2,J3 (COMPARE SPIN INTEGRALS) C C IF BOTH SPIN INTEGRALS ARE ZERO, THERE IS NO PURPOSE IN C CALCULATING THE ANGULAR INTEGRALS C IF(ABS(SPINDT).LT.EPS.AND.ABS(SPINEX).LT.EPS) GO TO 155 C I=2 CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,IANGD,IANGE,KK1,KK2,KK3,KK4) IF (IANGD.EQ.0) CALL J23ANG(IRHO,ISIG,IRHOP,ISIGP,JANGDI) C C IF THE DIRECT SPIN RECOUPLING COEFFICIENT IS ZERO, WE NEED NOT C CALCULATE THE CORRESPONDING ORBITAL RECOUPLING COEFFICIENT C IF(ABS(SPINDT).LT.EPS) GO TO 121 IF (IANGD.EQ.0) IANGD=1 C C --- DIRECT ANGULAR INTEGRAL C C CONSIDER ALL ALLOWED K-VALUES C IF (.NOT.FAILAD) THEN DO 114 JK1=KD1,KD2,2 J1(NJ1S)=2*JK1-1 IF (IANGD.NE.2) THEN FREE(NJ1S) = KD2.NE.KD1 CALL NJGRAF(ANGDIR,FAILAD) IANGD=2 IF (FAILAD) GO TO 121 CALL KNJ(J6CAD,J7CAD,J8CAD,J9CAD,JWCAD,K6AD,K7AD,K8AD,K9AD, + KWAD,JDELAD,LDELAD,SMVRAD,MPAD, + J6PAD,J7PAD,J8PAD,J9PAD,JWRDAD,NLSMAD,NBJAD,NB6JAD, + K6CPAD,K7CPAD,K8CPAD,K9CPAD,JSM4AD,JSM5AD,JSM6AD, + IN6JAD) ELSE CALL GENSUM(J6CAD,J7CAD,J8CAD,J9CAD,JWCAD,K6AD,K7AD,K8AD,K9AD, + KWAD,JDELAD,LDELAD,SMVRAD,MPAD, + J6PAD,J7PAD,J8PAD,J9PAD,JWRDAD,NLSMAD,NBJAD,NB6JAD, + K6CPAD,K7CPAD,K8CPAD,K9CPAD,JSM4AD,JSM5AD,JSM6AD, + IN6JAD,ANGDIR) ENDIF C C ADD INTO THE COEFFICIENT OF THE SLATER INTEGRAL C MULTD=1 WHEN A DIRECT INTEGRAL COEFFICIENT HAS BEEN CALCULATED - C FOR USE, SEE PRNTWT C MULTD=1 IF(IBUG1.GT.1) WRITE(IWRITE,310) ANGDIR 114 AMULT(JK1)=AMULT(JK1)+ANGDIR*BDIRCT ENDIF C C IF THE EXCHANGE SPIN RECOUPLING COEFFICIENT IS ZERO, WE NEED NOT C CALCULATE THE CORRESPONDING ORBITAL RECOUPLING COEFFICIENT C 121 IF(ABS(SPINEX).LT.EPS) GO TO 155 C C --- EXCHANGE ANGULAR INTEGRAL C C --- MODIFY J2 AND J3 ARRAYS TO CALCULATE EXCHANGE TERMS C IF (IANGE.NE.0) GO TO 279 I=2 CALL MODJ23(I) IANGE=1 C C CONSIDER ALL ALLOWED K-VALUES C 279 IF (.NOT.FAILAE) THEN DO 115 JK1=KE1,KE2,2 J1(NJ1S) = 2*JK1-1 IF(IANGE.NE.2) THEN FREE(NJ1S) = KE2.NE.KE1 CALL NJGRAF(ANGEX,FAILAE) IANGE=2 IF (FAILAE) GO TO 155 CALL KNJ(J6CAE,J7CAE,J8CAE,J9CAE,JWCAE,K6AE,K7AE,K8AE,K9AE, + KWAE,JDELAE,LDELAE,SMVRAE,MPAE, + J6PAE,J7PAE,J8PAE,J9PAE,JWRDAE,NLSMAE,NBJAE,NB6JAE, + K6CPAE,K7CPAE,K8CPAE,K9CPAE,JSM4AE,JSM5AE,JSM6AE, + IN6JAE) ELSE CALL GENSUM(J6CAE,J7CAE,J8CAE,J9CAE,JWCAE,K6AE,K7AE,K8AE,K9AE, + KWAE,JDELAE,LDELAE,SMVRAE,MPAE, + J6PAE,J7PAE,J8PAE,J9PAE,JWRDAE,NLSMAE,NBJAE,NB6JAE, + K6CPAE,K7CPAE,K8CPAE,K9CPAE,JSM4AE,JSM5AE,JSM6AE, + IN6JAE,ANGEX) ENDIF C C MULTE=1 WHEN AN EXCHANGE INTEGRAL COEFFICIENT HAS BEEN CALCULATED C MULTE=1 IF(IBUG1.GT.1) WRITE(IWRITE,312) ANGEX 115 BMULT(JK1)=BMULT(JK1)-ANGEX*BEXCHG ENDIF C C IF FAILSD OR FAILAD IS TRUE, AND FAILSE OR FAILAE IS TRUE, C THEN BOTH THE DIRECT AND EXCHANGE COEFFICIENTS WILL BE ZERO. C 155 IF ((FAILSD.OR.FAILAD).AND.(FAILSE.OR.FAILAE)) GO TO 527 C 154 CONTINUE 153 CONTINUE 152 CONTINUE 151 CONTINUE C C === INCLUDE MULTIPLICATIVE FACTORS COMMON TO ALL TERMS WITHIN THIS C FOUR-FOLD SUMMATION C IF(MULTD.EQ.0) GO TO 525 DO 518 JK1=KD1,KD2,2 518 AMULT(JK1)=AMULT(JK1)*XMULT*RMEDIR(JK1)*ADIRCT 525 IF(MULTE.EQ.0) GO TO 527 DO 519 JK1=KE1,KE2,2 519 BMULT(JK1)=BMULT(JK1)*XMULT*RMEEX(JK1)*AEXCHG C C --- PRINT OUT THE VALUES OF THE COEFFICIENTS OF THE SLATER INTEGRALS C C THE SUBROUTINE PRNTWT IS CALLED FROM RKWTS C 527 RETURN C C *** DEFINITION OF DIMENSION LIST C C RMEDIR(K),K=KD1,KD2,2 - DIRECT REDUCED MATRIX ELEMENT PRODUCT C RMEEX(K),K=KE1,KE2,2 - EXCHANGE REDUCED MATRIX ELEMENT PRODUCT C KD1,KE1 ARE ALWAYS .GE. 1 C KD2,KE2 ARE .LE. 1+2*MAX(L-VALUE) C NBAR(I), I=1,IHSH - NUMBER OF SPECTATOR ELECTRONS IN EACH SHELL C THE OTHER ARRAYS ARE DEFINED IN NJGRAF. C 999 PRINT *,' FANO: FIXED LENGTH 99 TOO SHORT IN /XATION/.' STOP C END C*********************************************************************** SUBROUTINE KNJ(JD6C,JD7C,JD8C,JD9C,JDWC,JD6,JD7,JD8,JD9,KDW,JDDEL, + LDDEL,DSUMVR,MDP, + JD6P,JD7P,JD8P,JD9P,JDWORD,NDLSUM,NDBJ,NDB6J, + KD6CP,KD7CP,KD8CP,KD9CP,JDSUM4,JDSUM5,JDSUM6,INVD6J) C PARAMETER (LL76= 21*3+19) PARAMETER(KFL1=LL76, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20, + KFLS=12,KFLN=10,KFLV=40) C LOGICAL DSUMVR DIMENSION JD6(KFL6),JD7(KFL7),JD8(KFL8), + JD9(KFL9),KDW(6,KFLW),LDDEL(KFLW,2),DSUMVR(KFL1) C DIMENSION JD6P(KFLV),JD7P(KFLV),JD8P(KFLV),JD9P(KFLV), + JDWORD(6,KFLW), + NDBJ(KFLN),NDB6J(KFLN),KD6CP(KFLN),KD7CP(KFLN),KD8CP(KFLN), + KD9CP(KFLN),JDSUM6(KFLS),JDSUM4(KFLS,KFLW),JDSUM5(KFLS,KFLW), + INVD6J(KFLW) C LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP C COMMON/SUMARG/J6P(KFLV),J7P(KFLV),J8P(KFLV),J9P(KFLV), + JWORD(6,KFLW),NLSUM, + NBJ(KFLN),NB6J(KFLN),K6CP(KFLN),K7CP(KFLN),K8CP(KFLN),K9CP(KFLN), + JSUM6(KFLS),JSUM4(KFLS,KFLW),JSUM5(KFLS,KFLW),INV6J(KFLW) C C C THIS ROUTINE STORES DATA FOR FUTURE CALLS TO GENSUM C JD6C=J6C JD7C=J7C JD8C=J8C JD9C=J9C JDWC=JWC JDDEL=JDEL MDP=MP NDLSUM=NLSUM IF(J6C.NE.0) THEN DO 800 I=1,J6C 800 JD6(I)=J6(I) ENDIF IF(J7C.NE.0) THEN DO 801 I=1,J7C 801 JD7(I)=J7(I) ENDIF IF(J8C.NE.0) THEN DO 802 I=1,J8C 802 JD8(I)=J8(I) ENDIF IF(J9C.NE.0) THEN DO 803 I=1,J9C 803 JD9(I)=J9(I) ENDIF IF(JWC.NE.0) THEN DO 804 J=1,JWC DO 804 I=1,6 804 KDW(I,J)=KW(I,J) DO 805 I=1,JWC 805 INVD6J(I)=INV6J(I) ENDIF IF(JDEL.NE.0) THEN DO 806 J=1,JDEL LDDEL(J,1)=LDEL(J,1) 806 LDDEL(J,2)=LDEL(J,2) ENDIF IF(MP.NE.0) THEN DO 807 I=1,MP 807 DSUMVR(I)=SUMVAR(I) ENDIF IF(NLSUM.NE.0) THEN DO 808 I=1,NLSUM NDBJ(I)=NBJ(I) NDB6J(I)=NB6J(I) KD6CP(I)=K6CP(I) KD7CP(I)=K7CP(I) KD8CP(I)=K8CP(I) 808 KD9CP(I)=K9CP(I) ENDIF DO 809 I=1,KFLV JD6P(I)=J6P(I) JD7P(I)=J7P(I) JD8P(I)=J8P(I) 809 JD9P(I)=J9P(I) DO 810 I=1,KFLS JDSUM6(I)=JSUM6(I) DO 810 J=1,KFLW JDSUM4(I,J)=JSUM4(I,J) 810 JDSUM5(I,J)=JSUM5(I,J) DO 811 J=1,KFLW DO 811 I=1,6 811 JDWORD(I,J)=JWORD(I,J) C C RETURN END C*********************************************************************** SUBROUTINE H0WTS(ISIG,ISIGP,TIMES,ICAL) C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C PARAMETER(LL43= 21*2+3, LL74= 21+4, LL75= 21+2, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2L=3*KFL2,KFL2S=3*KFL2-3) C LOGICAL FAIL,FREE C COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/COUPLE/NJ1S,NJ23S,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 COMMON/HOLD/J2STO(LL76),J3STO(LL76),J2ANG(LL76),J3ANG(LL76) COMMON/INFORM/IREAD,IWRITE,IPUNCH COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75), 1 J1QN1(LL43,3),J1QN2(LL43,3),IJFUL(LL75) COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15, 1 M16,M17,M18,M19,M20 303 FORMAT(7H ISIG =,I3,8H ISIGP =,I3) 306 FORMAT(5H J1 =,I8,36I3) 308 FORMAT(3H J2,18X,2HJ3/(3I5,I10,2I5)) 309 FORMAT(8H TIMES =,F10.6,8H RECUP =,F10.6) 310 FORMAT(8H TIMES =,F10.6,8H IDELP =,I3) C C --- OFF-DIAGONAL ONE-ELECTRON HAMILTONIAN MATRIX ELEMENT C ICAL=0 CALL SETM ICOUNT=0 C C TEST THAT FINAL ANGULAR MOMENTA ARE EQUAL C DO 8 K=2,3 IF(J1QN1(M6,K).NE.J1QN2(M6,K)) GO TO 7 8 CONTINUE C C --- DETERMINE INTERACTING SHELLS, ISIG ON L.H.S., ISIGP ON R.H.S., C FOR NON-ZERO ONE-ELECTRON MATRIX ELEMENT, N-1 ELECTRONS MUST BE C COMMON TO BOTH SIDES. THUS THE SUM OF N(I) = NOSH1(I)-NOSH2(I), C I=1,IHSH MUST BE EQUAL TO 0 OR 2 . THUS AT NO STAGE CAN C N(I) BE GREATER THAN 1 . IF THIS SUM IS ZERO, THE TWO C CONFIGURATIONS ARE MADE UP FROM THE SAME ELECTRONS, WITH TWO C DIFFERENT COUPLING SCHEMES. SINCE THE SPHERICAL HARMONICS ARE C EIGENFUNCTIONS OF DEL**2 , THE ORTHOGONALITY OF THE TWO COUPLING C SCHEMES WILL BE MAINTAINED AND ORTHOGONALITY GIVES A ZERO RESULT. C DO 9 I=1,IHSH N=NOSH1(I)-NOSH2(I) IF(ABS(N).GT.1) GO TO 7 IF(N) 11,9,12 11 ISIGP=I GO TO 13 12 ISIG=I 13 ICOUNT=ICOUNT+1 9 CONTINUE IF(ICOUNT.NE.2) GO TO 7 IF(IBUG4.GE.2) WRITE(IWRITE,303) ISIG,ISIGP LSIG=LJ(ISIG) LSIGP=LJ(ISIGP) C C THE ANGULAR MOMENTA OF THE INTERACTING ELECTRONS MUST BE EQUAL C IF(LSIG.NE.LSIGP) GO TO 7 C C THE SPECTATOR SHELLS MUST HAVE MATCHING QUANTUM NUMBERS C DO 16 J=1,IHSH IF(J.EQ.ISIG.OR.J.EQ.ISIGP) GO TO 16 DO 19 K=1,3 IF(J1QN1(J,K).NE.J1QN2(J,K)) GO TO 7 19 CONTINUE 16 CONTINUE C C --- TEST ON TRIANGULAR CONDITIONS C IN3=2*LSIG DO 20 K=2,3 IN1=J1QN1(ISIG,K) IN2=J1QN2(ISIG,K) IN4=J1QN1(ISIGP,K) IN5=J1QN2(ISIGP,K) IF(IN1.GT.IN2+IN3 .OR.IN1.LT.ABS(IN2-IN3)) GO TO 7 IF(IN4.GT.IN5+IN3 .OR.IN4.LT.ABS(IN5-IN3)) GO TO 7 20 IN3=1 C C --- CALCULATE FRACTIONAL PARENTAGE COEFFICIENTS C TIMES=1.0 IF(LSIG.EQ.0) GO TO 26 N=NOSH1(ISIG) IVI=J1QN1(ISIG,1) ILI=(J1QN1(ISIG,2)-1)/2 ISI=J1QN1(ISIG,3) IVJ=J1QN2(ISIG,1) ILJ=(J1QN2(ISIG,2)-1)/2 ISJ=J1QN2(ISIG,3) CALL CFP(LSIG,N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP) TIMES=TIMES*COEFP IF(ABS(TIMES).LT.EPS) GO TO 7 N=NOSH2(ISIGP) IVI=J1QN2(ISIGP,1) ILI=(J1QN2(ISIGP,2)-1)/2 ISI=J1QN2(ISIGP,3) IVJ=J1QN1(ISIGP,1) ILJ=(J1QN1(ISIGP,2)-1)/2 ISJ=J1QN1(ISIGP,3) CALL CFP(LSIG,N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP) TIMES=TIMES*COEFP IF(ABS(TIMES).LT.EPS) GO TO 7 C C --- SET UP J2 AND J3 ARRAYS C 26 M1=IHSH-2 M2=M6-2 J2(1,1)=ISIG J2(1,2)=M11 J2(1,3)=M9 J3(1,1)=ISIGP J3(1,2)=M11 J3(1,3)=M10 IF(ISIG.EQ.1) GO TO 29 J2(2,1)=1 GO TO 30 29 J2(2,1)=M9 30 IF(ISIG.EQ.2) GO TO 32 J2(2,2)=2 GO TO 33 32 J2(2,2)=M9 33 J2(2,3)=M4 IF(ISIGP.EQ.1) GO TO 35 J3(2,1)=1 GO TO 36 35 J3(2,1)=M10 36 IF(ISIGP.EQ.2) GO TO 38 J3(2,2)=2 GO TO 39 38 J3(2,2)=M10 39 J3(2,3)=M7 IF(IHSH.LT.3) GO TO 40 DO 42 J=3,IHSH J2(J,1)=M1+J J2(J,3)=M1+J+1 J3(J,1)=M2+J IF(J.EQ.IHSH) GO TO 44 J3(J,3)=M2+J+1 GO TO 45 44 J3(J,3)=M1+J+1 45 IF(J.EQ.ISIG) GO TO 47 J2(J,2)=J GO TO 48 47 J2(J,2)=M9 48 IF(J.EQ.ISIGP) GO TO 50 J3(J,2)=J GO TO 42 50 J3(J,2)=M10 42 CONTINUE C C --- STORE J2 AND J3 ARRAYS FOR USE IN SPIN RECOUPLING COEFFICIENT C 40 I1=0 DO 51 J=1,IHSH DO 52 K=1,3 J2STO(K+I1)=J2(J,K) 52 J3STO(K+I1)=J3(J,K) 51 I1=I1+3 C C --- ORBITAL RECOUPLING COEFFICIENT C J1(M11)=LSIG+LSIG+1 K=2 C C --- SET J1 ARRAY C 64 DO 53 J=1,IHSH IF(ISIG.EQ.J) GO TO 55 J1(J)=J1QN1(J,K) GO TO 53 55 J1(J)=J1QN2(ISIG,K) 53 CONTINUE DO 56 J=M4,M6 56 J1(J)=J1QN1(J,K) DO 57 J=M7,M8 57 J1(J)=J1QN2(J-M3,K) J1(M9)=J1QN1(ISIG,K) J1(M10)=J1QN2(ISIGP,K) NJ1S=M11 NJ23S=M4 DO 100 J=1,M11 100 FREE(J)=.FALSE. IF(IBUG4.LT.2.OR.IBUG3.EQ.1) GO TO 77 WRITE(IWRITE,306) (J1(J),J=1,M11) WRITE(IWRITE,308) ((J2(J,KL),KL=1,3),(J3(J,KL),KL=1,3),J=1,IHSH) C C --- CALCULATE RECOUPLING COEFFICIENT C 77 CALL NJGRAF(RECUP,FAIL) TIMES=TIMES*RECUP IF(IBUG4.LT.2) GO TO 58 WRITE(IWRITE,309) TIMES,RECUP 58 IF(ABS(TIMES).LT.EPS) GO TO 7 C C --- SPIN RECOUPLING COEFFICIENT C IF(K.EQ.3) GO TO 60 J1(M11 )=2 K=3 I1=0 DO 62 J=1,IHSH DO 63 KK=1,3 J2(J,KK)=J2STO(KK+I1) 63 J3(J,KK)=J3STO(KK+I1) 62 I1=I1+3 GO TO 64 C C --- INCLUDE MULTIPLICATIVE FACTORS C 60 K=0 IF(ISIG-ISIGP) 65,70,66 65 JSIG=ISIG+1 DO 67 J=JSIG,ISIGP 67 K=NOSH1(J)+K GO TO 70 66 JSIGP=ISIGP+1 DO 68 J=JSIGP,ISIG 68 K=NOSH2(J)+K 70 TIMES = (1-MOD(K,2)*2)*SQRT(REAL(NOSH1(ISIG)*NOSH2(ISIGP)))*TIMES ICAL=1 IF(IBUG4.LT.2) GO TO 7 WRITE(IWRITE,310) TIMES,K 7 RETURN END C*********************************************************************** SUBROUTINE NJGRAF(RECUP,FAIL) C*********************************************************************** C C WRITTEN BY C A. BAR-SHALOM AND M. KLAPISCH C RACAH INSTITUTE OF PHYSICS C HEBREW UNIVERSITY C 91904 JERUSALEM C ISRAEL. C C*********************************************************************** C C-----THIS IS THE MAIN PROGRAM. IT HANDLES ALL THE ANALYSIS OF THE C RECOUPLING COEFFICIENT WITHOUT REFERRING EXPLICITLY TO THE VALUES C OF ANGULAR MOMENTA WHICH ARE IN J1(J) EXCEPT FOR ZERO IN CASE FREE C =.FALSE. LIKE NJSYM IT PREPARES ARRAYS OF ARGUMENTS FOR PHASE C FACTORS, (2*J+1) FACTORS AND 6J COEFFICIENTS TO BE COMPUTED IN C GENSUM,WHICH CAN ALSO BE CALLED SEPARATELY WHEN ONLY THE NUMERICAL C VALUES OF ANGULAR MOMENTA CHANGE.THESE VARIABLE ANGULAR MOMENTA C SHOULD BE DECLARED FREE(J)=.TRUE. SO THAT THE FORMULA PREPARED FOR C GENSUM SHOULD BE CORRECT WHEN J1 IS NOT ZERO. C FAIL WILL BE TRUE WHEN THE RECOUPLING COEFFICIENT IS ZERO BECAUSE C OF UNSATISFIED DELTA OR OTHER SIMILAR CAUSES. C C*********************************************************************** C C THIS VERSION HAS BEEN MODIFIED TO HOLD THE ARRAY DIMENSIONS IN C PARAMETER STATEMENTS. THE DIMENSIONS ARE LABELLED: C C KFL1 - DIMENSION OF THE J1 AND FREE ARRAYS IN /COUPLE/, AND THE C FIRST DIMENSION OF THE LINE AND LCOL ARRAYS IN /TREE/. C ALSO THE DIMENSION OF THE SUMVAR ARRAY IN /ARGU/, AND C OF THE INVER ARRAY IN ROUTINE SPRATE. IT IS TESTED FOR M C ON ENTRY TO NJGRAF, AND FOR MP IN ROUTINE SPRATE. C C KFL2 - DIMENSION OF THE J2 AND J3 ARRAYS IN /COUPLE/. C THE LENGTH OF THESE ARRAYS IS CHECKED ON ENTRY TO NJGRAF. C C KFL2A = (2*KFL2) - DIMENSION OF THE J23, ARROW AND TABS ARRAYS IN C /TREE/. C KFL2B = (4*KFL2) - DIMENSION OF THE JDIAG,ARR,IL AND IH ARRAYS C IN /GRAPH/, AND OF THE IAL ARRAY IN /BUILD/ C KFL2C = (2*KFL2+2) - DIMENSION OF THE NPOINT ARRAY IN /GRAPH/. C C KFL6 - DIMENSION OF THE J6 ARRAY IN /ARGU/. TESTED IN SPRATE. C C KFL7 - DIMENSION OF THE J7 ARRAY IN /ARGU/. TESTED IN SPRATE. C C KFL8 - DIMENSION OF THE J8 ARRAY IN /ARGU/. TESTED IN SPRATE. C C KFL9 - DIMENSION OF THE J9 ARRAY IN /ARGU/. TESTED IN SPRATE. C C KFLW - DIMENSION OF THE JW(OR KW) AND LDEL ARRAYS IN /ARGU/, AND C OF THE JWORD AND INV6J ARRAYS IN /SUMARG/. ALSO THE SECOND C DIMENSION OF THE JSUM4 AND JSUM5 ARRAYS IN /SUMARG/. C IN ADDITION IT GIVES THE DIMENSIONS OF A NUMBER OF C TEMPORARY WORKING ARRAYS IN ROUTINES SPRATE AND GENSUM. C KFLW IS TESTED IN SPRATE. C C KFLS - DIMENSION OF THE JSUM6 ARRAY AND FIRST DIMENSION OF THE C JSUM4 AND JSUM5 ARRAYS IN /SUMARG/. ALSO GIVES THE C DIMENSIONS OF SOME TEMPORARY WORKING ARRAYS IN SPRATE C AND GENSUM. KFLS IS THE MAXIMUM NUMBER OF SUMMATION C VARIABLES IN A PARTICULAR SUM, AND IS TESTED IN SPRATE. C C KFLN - DIMENSION OF THE NBJ, NB6J, K6CP, K7CP, K8CP AND K9CP C ARRAYS IN /SUMARG/. KFLN IS THE MAXIMUM NUMBER OF C SUMS ALLOWED, AND IS TESTED IN ROUTINE SPRATE. C C KFLV - DIMENSION OF THE J6P, J7P, J8P AND J9P ARRAYS IN C /SUMARG/, AND OF THE JNS ARRAY IN ROUTINE VAR. C KFLV IS TESTED IN VAR. C C KFLZ - DIMENSION OF THE JZERO ARRAY IN /ZER/. KFLZ IS TESTED IN C ROUTINE ZERO. C C****************************************************************** C C OTHER CHANGES: C C (1) THIS VERSION HAS BEEN ALTERED TO USE THE RACAH C COEFFICIENT ROUTINE WRITTEN BY STAN SCOTT AND ALAN HIBBERT C (CPC 28 189-200 1982). SUBROUTINE FACTT MUST BE CALLED TO SET UP C THE ARRAY GAM, WHICH HOLDS LN(N!) BEFORE NJGRAF IS CALLED. C (2) THE SUBROUTINES ORDER AND SETDIM HAVE BEEN RENAMED AS C ORDTRI AND SETDM TO AVOID CONFLICT IN THE QUB CODES. C (3) COMMON BLOCKS /ARGU/ AND /SUMARG/ HAVE BEEN REMOVED FROM C SUBROUTINE GENSUM, AND THEIR CONTENTS ARE TRANSFERRED VIA C THE ARGUMENT LIST INSTEAD. /SUMARG/ HAS THUS BEEN INCLUDED C IN SUBROUTINE NJGRAF TO ACHIEVE THIS. C (4) THE TIMING ROUTINES HAVE BEEN REMOVED FROM THE NJGRAF C PACKAGE AS THEY WERE CDC SPECIFIC AND SEEMED UNNECESSARY. C C**************************************************************** C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20, + KFLS=12,KFLN=10,KFLV=40) C LOGICAL FAIL,FIND, TABS,CUT, FREE, SUMVAR C INTEGER ARROW, ARR,TAB1 COMMON/TREE/J23(KFL2A,3),ARROW(KFL2A,3),LINE(KFL1,2),LCOL(KFL1,2), + TABS(KFL2A),NBTR COMMON/CUTDIG/CUT COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C COMMON/COUPLE/M,N,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) C COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP C COMMON/SUMARG/J6P(KFLV),J7P(KFLV),J8P(KFLV),J9P(KFLV), + JWORD(6,KFLW),NLSUM, + NBJ(KFLN),NB6J(KFLN),K6CP(KFLN),K7CP(KFLN),K8CP(KFLN),K9CP(KFLN), + JSUM6(KFLS),JSUM4(KFLS,KFLW),JSUM5(KFLS,KFLW),INV6J(KFLW) C CHARACTER*6 NAMSUB COMMON /NAM/NAMSUB C C TEST THE LENGTH OF THE J1 ARRAY C IF(M+1.GT.KFL1) THEN PRINT 100,M+1,KFL1 100 FORMAT(' DIMENSION ERROR IN NJGRAF. M+1=',I5,' KFL1=',I5) STOP ENDIF C C TEST THE LENGTH OF THE J2 AND J3 ARRAYS C IF(N-1.GT.KFL2) THEN PRINT 101,N-1,KFL2 101 FORMAT(' DIMENSION ERROR IN NJGRAF: N-1=',I5,' KFL2=',I5) STOP ENDIF C C 1.0 C BUILDING UP OF THE UNSTRUCTURED GRAPH C FAIL=.FALSE. J6C=0 J7C=0 J8C=0 J9C=0 JWC=0 JDEL=0 CALL SETDM NFIN=0 CUT=.FALSE. CALL SETTAB(FAIL) M=M+1 IF(FAIL)GO TO 7 M=M-1 JF=0 JF1=0 C C 1.1 C LOCATING AND HANDLING OF ZEROS C CALL ZERO(JF1,JF,FAIL) IF(FAIL)GO TO 7 MP=M IF(NBTR.EQ.0)GO TO 6 JUMP=1 C C 2. C BUILDING OF A FLAT DIAGRAM OUT OF THE UNSTRUCTURED GRAPH. C THERE MAY BE SEVERAL FLAT DIAGRAMS OUT OF THE ORIGINAL C GRAPH IN CASE OF POSSIBLE CUTS. THEN THE FLAT DIAGRAMS C WILL HAVE FREE ENDS. C 1 CALL DIAGRM(JUMP) NFIN=MAX(0,NFREE-2) IF(NFIN.EQ.0) GO TO 2 JUMP=3 C C 2.1 C HANDLING OF FREE ENDS IF A CUT WAS FOUND C CALL CUTNL(FAIL) IF(FAIL) GO TO 7 GO TO 5 2 JUMP=2 IF(NFREE-1)5,3,4 3 CALL CUT1L(FAIL) IF(FAIL) GO TO 7 GO TO 5 4 CALL CUT2L(FAIL) IF(FAIL) GO TO 7 5 NBTR=NBTR+NFIN IF(NBTR.NE.0) CUT=.TRUE. C C 3. C ANALYSIS OF THE FLAT DIAGRAM C CLOSED CIRCUITS OF INCREASING ORDER NC ARE SEARCHED, ANALYSED, AND C TAKEN OUT OF THE FLAT DIAGRAM, THUS REDUCING THE NUMBER NBNODE OF C NODES: C NC=0 10 NC=NC+1 CALL SEARCH(FIND) IF(.NOT.FIND) GO TO 10 NCP=NC-2 JPOL=0 IF(M.EQ.MP.AND.NC.GT.3)CALL SETDM IF(IPARTL.GT.2)CALL POLYGN(JPOL) GO TO (11,12,13,14),NC 11 CALL LOLPOP(FAIL) IF(FAIL) GO TO 7 GO TO 15 12 CALL BUBBLE(JPOL,FAIL) IF(FAIL)GO TO 7 GO TO 15 13 CALL TRIANG(FAIL) IF(FAIL)GO TO 7 GO TO 15 14 CALL SQUARE 15 NBNODE=NBNODE-2 IF(NBNODE.EQ.0) GO TO 9 IFIRST=IH(1) ILAST=IH(NBNODE) C PRINTJ IS AN ALL PURPOSE PRINTING SUBROUTINE CALLED FROM MANY PLA CALL PRINTJ(NAMSUB,8) IF(NBNODE.EQ.NFIN) GO TO 9 NC=NCP C C PROCEED TO OTHER CIRCUITS OF ORDER NC-1 C GO TO 10 9 IF(NBTR.EQ.0)GO TO 6 IF(JUMP.EQ.3) CALL ORDTRI C C AT THIS STAGE THE FLAT DIAGRAM HAS BEEN REDUCED TO NODES C INVOLVING FREE ENDS. PROCEED TO BUILD OTHER FLAT DIAGRAMS C IF NECESSARY. C GO TO 1 ! '05May3-4 tmpmod (before correcting ORDTRI 'o5May31): C if(IFIRST.le.ILAST) GO TO 1 ! not IF(NBNODE.GT.0) GO TO 1 ! C print"(' NJGRAF: tempmod, setting RECUP=0 because of FIRST>LAST')" C C ALL PARTS OF THE ORIGINAL GRAPH HAVE BEEN REDUCED C 7 RECUP=0. M=M-1 RETURN 6 CALL PRINTJ('NJGRAF',0) C C 4. PREPARATION OF THE RESULTS, AND SEPARATION IN SEVERAL SUMS C IF CUTS HAVE BEEN DETECTED, ALSO IN THE FLAT DIAGRAM ITSELF C CALL SPRATE(M) M=M-1 C C 5. GENSUM COMPUTES THE NUMERICAL VALUE OF THE RECOUPLING COEFFICIEN C CALL GENSUM(J6C,J7C,J8C,J9C,JWC,J6,J7,J8,J9,KW,JDEL,LDEL, * SUMVAR,MP, J6P,J7P,J8P,J9P,JWORD,NLSUM,NBJ,NB6J, + K6CP,K7CP,K8CP,K9CP,JSUM4,JSUM5,JSUM6,INV6J, RECUP) RETURN END C*********************************************************************** SUBROUTINE BUBBLE(JPOL,FAIL) C C -- REDUCES A CIRCUIT OF ORDER 2, GIVING DELTA FUNCTION AND PHASE C FACTORS C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C LOGICAL FAIL C INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP C CHARACTER*6 NAMSUB COMMON /NAM/NAMSUB C C NAMSUB = 'BUBBLE' K2=2 K23=3 I1=1 I2=1 IT1=NPOINT(1) IT2=NPOINT(2) IF(IT2.NE.ILAST) GO TO 2 IF(IT1.EQ.IFIRST) GO TO 1 IT2=IT1 IT1=ILAST 1 I1=-1 K23=2 I2=2 C C 2 CALL PHASE(IT1,JDIAG,KFL2B) C ARISING FROM NON-CYCLIC PERMUTATION OF ARGUMENTS IN TRIAD L: 2 J7(J7C+1) = JDIAG(IT1,1) J7(J7C+2) = JDIAG(IT1,2) J7C = J7C + 3 J7(J7C) = JDIAG(IT1,3) K = ABS((3*ARR(IT2,1)+2*ARR(IT2,2)+ARR(IT2,3))/2) + 1 IF (K.NE.4) THEN C CALL PHASE2(JDIAG(IT2,K)) ! ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = JDIAG(IT2,K) ENDIF C IF(NBNODE.EQ.2) GO TO 7 IL1=IL(IT2)+I1 IT=IH(IL1) ARR(IT,K23)=ARR(IT1,K23) L=JDIAG(IT1,K23) L1=JDIAG(IT,K23) JDIAG(IT,K23)=L IF(JPOL.NE.1)GO TO 3 MP=MP-1 KW(2,JWC)=L J6(J6C-1)=L J6(J6C)=L IF(K.EQ.2)J8(J8C)=L GO TO 4 3 CALL DELTA(L,L1,FAIL) IF(FAIL)GO TO 7 4 TAB1(L,I2)=IT IF(IT1.EQ.ILAST)GO TO 6 IF(IT2.NE.ILAST)GO TO 9 TAB1(L,1)=IH(2) IL1=2 K2=1 9 DO 5 I=IL1,NBNODE IT=IH(I) IL(IT)=I-K2 5 IH(I-K2)=IT 6 J9(J9C+1)=L J9C=J9C+2 J9(J9C)=L 7 RETURN END C*********************************************************************** SUBROUTINE CHANGE(L,K) C C EXCHANGES THE FREE ENDS IN EITHER FIRST OR LAST TRIAD OF JDIAG. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2) PARAMETER (KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C LOGICAL SUMVAR INTEGER ARR,TAB1 COMMON /ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), A J9(KFL9),JW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C C CALL PHASE(L,JDIAG,KFL2B) -- PHASE FACTOR C ARISING FROM NON-CYCLIC PERMUTATION OF ARGUMENTS IN TRIAD L: J7(J7C+1) = JDIAG(L,1) J7(J7C+2) = JDIAG(L,2) J7C = J7C + 3 J7(J7C) = JDIAG(L,3) C JP=JDIAG(L,K) JDIAG(L,K)=JDIAG(L,1) JDIAG(L,1)=JP JAR=ARR(L,K) ARR(L,K)=ARR(L,1) ARR(L,1)=JAR RETURN END C*********************************************************************** SUBROUTINE CHVAR(JP,NBC,KBC,JT,JINV,NSUM) C C - CHANGES THE ORDER OF SUMMATION VARIABLE TO BE ABLE TO PERFORM C SEPARATELY THE SUMMATIONS IN GENSUM. C LOGICAL JT(NSUM) DIMENSION JP(NBC),JINV(NSUM) KB=KBC+1 IF(KB.GT.NBC)GO TO 2 DO 1 I=KB,NBC JK=JP(I) IF(.NOT.JT(JK))GO TO 1 KBC=KBC+1 JP(I)=JP(KBC) JP(KBC)=JINV(JK) 1 CONTINUE 2 RETURN END C*********************************************************************** SUBROUTINE CUT1L(FAIL) C C CUT ON ONE LINE THAT WAS LEFT AS A FREE END IN JDIAG. PUTS C CORRESPONDING DELTA IN J23. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C LOGICAL FAIL C INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP LOGICAL FREE COMMON/COUPLE/M,N,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) C C IT=ITFREE(1) J0=JDIAG(IT,1) CALL DELTA (J0,M,FAIL) IF(FAIL) GO TO 2 CALL DELTA(JDIAG(IT,3),JDIAG(IT,2),FAIL) IF(FAIL)GO TO 2 JDIAG(IT+1,3)=JDIAG(IT,3) IF(ARR(IT,2)-ARR(IT,3)) 4,3,5 3 ARR(IT+1,3)=1 ARR(IT-1,2)=-1 GO TO 5 4 ARR(IT+1,3)=-1 ARR(IT-1,2)=1 5 J9C=J9C+1 J9(J9C)=JDIAG(IT,3) J=2 CALL ZERO(J,J0,FAIL) IF(FAIL) GO TO 2 IL1=IL(IT+1) DO 1 I=IL1,NBNODE IT=IH(I) ILP=I-1 IL(IT)=ILP 1 IH(ILP)=IT NBNODE=NBNODE-1 2 CALL PRINTJ('CUT1L ',12) RETURN END C*********************************************************************** SUBROUTINE CUT2L(FAIL) C C CUT ON TWO LINES THAT WERE LEFT AS FREE ENDS IN JDIAG. PUTS C CORRESPONDING DELTA IN J23. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C LOGICAL FAIL C INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C LOGICAL TABS INTEGER ARROW COMMON/TREE/J23(KFL2A,3),ARROW(KFL2A,3),LINE(KFL1,2),LCOL(KFL1,2), + TABS(KFL2A),NBTR C LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP C C IT1=ITFREE(1) IT2=ITFREE(2) JT1=JDIAG(IT1,1) JT2=JDIAG(IT2,1) CALL DELTA(JT1,JT2,FAIL) IF(FAIL) GO TO 1 IF(ARR(IT1,1).EQ.ARR(IT2,1)) THEN C CALL PHASE2(JT1) ! ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = JT1 ENDIF ARR(IT2,1)=-ARR(IT1,1) JDIAG(IT2,1)=JT1 TAB1(JT1,2)=IT2 J9(J9C+1)=JT1 J9C=J9C+2 J9(J9C)=JT1 C CALL OTHERJ(0,JT1,L1,LC1,K1) C GIVES THE OTHER TRIAD (WHERE A GIVEN J OCCURS) AND ITS POSITION: C L1 = LINE(JT1,1) IF (L1.EQ.0 .OR. TABS(L1)) THEN K1 = 1 L1 = LINE(JT1,2) LC1 = LCOL(JT1,2) C ELSE K1 = 2 LC1 = LCOL(JT1,1) ENDIF C C CALL OTHERJ(0,JT2,L2,LC2,K2) C GIVES THE OTHER TRIAD WHERE A GIVEN J OCCURS AND ITS POSITION: C L2 = LINE(JT2,1) IF (L2.EQ.0 .OR. TABS(L2)) THEN L2 = LINE(JT2,2) LC2 = LCOL(JT2,2) C ELSE LC2 = LCOL(JT2,1) ENDIF C J23(L2,LC2)=JT1 LINE(JT1,K1)=L2 LCOL(JT1,K1)=LC2 ARROW(L2,LC2)=-ARROW(L1,LC1) 1 CALL PRINTJ('CUT2L ',12) RETURN END C*********************************************************************** SUBROUTINE CUTNL(FAIL) C C - EXAMINES THE CASE WHERE THERE ARE MORE THAN C TWO FREE ENDS, BUT THEY ARE CONTIGUOUS, SO THAT THE GRAPH CAN C BE CUT WITHOUT DESTROYING THE FLAT STRUCTURE. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C INTEGER ARROW LOGICAL TABS COMMON/TREE/J23(KFL2A,3),ARROW(KFL2A,3),LINE(KFL1,2),LCOL(KFL1,2), + TABS(KFL2A),NBTR C INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C COMMON/KEEP/JKP(2,3),JARR(2,3),IT2,IT3,IT5 ! for ORDTRI C LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP C LOGICAL FAIL C C NTF=ITFREE(NFREE)-ITFREE(1) IF(NTF.GT.NFREE) GO TO 8 ! 9 ??? IT2=ITFREE(1) IT3=ITFREE(NFREE) IT1=IT2-1 IT4=IT3+1 IF(NTF.EQ.NFREE)GO TO 2 JT=JDIAG(IT2,3) CALL DELTA(JT,JDIAG(IT3,2),FAIL) IF(FAIL) GO TO 8 ! will do instead of 9 IF(ARR(IT2,3).NE.ARR(IT3,2))GO TO 1 C CALL PHASE2(JT) ! ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = JT ARR(IT2,3)=-ARR(IT2,3) ARR(IT1,2)=-ARR(IT1,2) 1 JDIAG(IT3,2)=JT JDIAG(IT4,3)=JT J9(J9C+1)=JT J9C=J9C+2 J9(J9C)=JT NBTR=NBTR+NFREE IT5=0 GO TO 6 2 NFR=0 DO 3 IT5=IT2,IT3 NFR=NFR+1 IF(ITFREE(NFR).GT.IT5)GO TO 4 3 CONTINUE 4 JKP(1,1)=JDIAG(IT5,1) JARR(1,1)=-ARR(IT5,1) JKP(1,2)=JDIAG(IT2,3) JARR(1,2)=-ARR(IT2,3) JKP(1,3)=JDIAG(IT3,2) JARR(1,3)=-ARR(IT3,2) DO 5 J=1,3 JKP(2,J)=JDIAG(IT5,J) 5 JARR(2,J)=ARR(IT5,J) JDIAG(IT5,2)=JDIAG(IT3,2) ARR(IT5,2)=ARR(IT3,2) JDIAG(IT5,3)=JDIAG(IT2,3) ARR(IT5,3)=ARR(IT2,3) ILP=IL(IT2) IL(IT5)=ILP IH(ILP)=IT5 NBTR=NBTR+NFREE+2 C C CALL PHASE(IT5,JDIAG,KFL2B) -- PHASE FACTOR C ARISING FROM NON-CYCLIC PERMUTATION OF ARGUMENTS IN TRIAD L: J7(J7C+1) = JDIAG(IT5,1) J7(J7C+2) = JDIAG(IT5,2) J7C = J7C + 3 J7(J7C) = JDIAG(IT5,3) K = ABS((3*ARR(IT5,1)+2*ARR(IT5,2)+ARR(IT5,3))/2) + 1 IF(K.NE.4) THEN C CALL PHASE2(JDIAG(IT5,K)) ! ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = JDIAG(IT5,K) ENDIF C 6 IL1=IL(IT4) DO 7 I=IL1,NBNODE IT=IH(I) ILP=I-NFREE IL(IT)=ILP 7 IH(ILP)=IT NBNODE=NBNODE-NFREE NFIN=0 GO TO 8 9 FAIL=.TRUE. 8 CALL PRINTJ('CUTNL ',8) RETURN END C*********************************************************************** SUBROUTINE DELTA(JA,JB,FAIL) C C TEST FOR DELTA(JA,JB). IF THEY ARE SUMMATION VARIABLES THE SECOND C IS CHANGED INTO THE FIRST EVERYWHERE. IF THEY ARE FIXED THEIR C VALUE IS CHECKED AND FAIL PUT TO .TRUE. IF THEY DIFFER. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C LOGICAL FAIL C LOGICAL CUT COMMON/CUTDIG/CUT LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 C COMMON/DIM/J6CC,J7CC,J8CC,J9CC,JWCC,JDELC C LOGICAL FREE COMMON/COUPLE/M,N,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) C C IF(IBUG3.EQ.1) PRINT 1000, JA,SUMVAR(JA),JB,SUMVAR(JB) 1000 FORMAT(/' from DELTA JA=',I2,L2,5X,'JB=',I2,L2) IF(SUMVAR(JA).AND.SUMVAR(JB)) GO TO 2 IF(FREE(JA).OR.FREE(JB)) GO TO 1 IF(J1(JA).NE.J1(JB)) FAIL=.TRUE. CUT=.TRUE. GO TO 18 1 JDEL=JDEL+1 LDEL(JDEL,1)=JA LDEL(JDEL,2)=JB SUMVAR(JA)=.FALSE. SUMVAR(JB)=.FALSE. GO TO 18 2 IF(J6C.EQ.J6CC) GO TO 4 DO 3 I=J6CC+1,J6C 3 IF(J6(I).EQ.JB) J6(I)=JA 4 IF(J7C.EQ.J7CC) GO TO 6 DO 5 I=J7CC+1,J7C IF(J7(I).EQ.JB) J7(I)=JA 5 CONTINUE 6 IF(J8C.EQ.J8CC) GO TO 8 DO 7 I=J8CC+1,J8C IF(J8(I).EQ.JB) J8(I)=JA 7 CONTINUE 8 IF(J9C.EQ.J9CC) GO TO 10 DO 9 I=J9CC+1,J9C IF(J9(I).EQ.JB) J9(I)=JA 9 CONTINUE 10 IF(JWC.EQ.JWCC) GO TO 15 DO 14 I=JWCC+1,JWC DO 13 J=1,6 IF(KW(J,I).EQ.JB) KW(J,I)=JA 13 CONTINUE 14 CONTINUE 15 IF(JDEL.EQ.JDELC) GO TO 18 DO 16 I=JDELC+1,JDEL IF(LDEL(I,1).EQ.JB) LDEL(I,1)=JA IF(LDEL(I,2).EQ.JB) LDEL(I,2)=JA 16 CONTINUE SUMVAR(JB)=.FALSE. 18 RETURN END C*********************************************************************** SUBROUTINE DIAGRM(JUMP) C C - BUILDS UP A FLAT DIAGRAM FROM THE TRIADS J23 AND PLACES THEM C IN JDIAG. ARROWS ARE IN ARR (INTEGER). THE DIAGRAM IS BUILT SO C AS TO MAXIMIZE THE NUMBER OF TRIADS INVOLVED, WITHIN A ONE- C -STEP-FORWARD-CHECK PROCESS. IF THE DIAGRAM DOES NOT INCLUDE C ALL THE NBTR TRIADS IT WILL HAVE 'FREE ENDS'. JDIAG HAS SIZE C TWICE THAT OF J23, BECAUSE THE PATH MAY PROCEED EITHER WAY. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP LOGICAL TABS INTEGER ARROW COMMON/TREE/J23(KFL2A,3),ARROW(KFL2A,3),LINE(KFL1,2),LCOL(KFL1,2), + TABS(KFL2A),NBTR INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C COMMON/BUILD/IAL(KFL2B),IF1,IF2,NODE C C LOGICAL FREE COMMON/COUPLE/M,N,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) C SAVE NB C C INITIALIZATION IF(JUMP-2) 16,1,17 16 NB=0 1 NB=NB+1 IF(TABS(NB)) GO TO 1 NODE=NBTR ILAST=NBTR DO 2 J=1,3 JDIAG(NODE,J)=J23(NB,J) 2 ARR(NODE,J)=ARROW(NB,J) TABS(NB)=.TRUE. DO 15 I=1,MP 15 IAL(I)=0 IF1=JDIAG(NODE,1) IF2=JDIAG(NODE,3) IAL(IF1)=1 IAL(IF2)=1 17 NTIME=0 I1=1 K1=1 K2=2 K3=3 3 JB=JDIAG(NODE,K2) C C CALL OTHERJ(0,JB,L,LC,KP) -- IP version (johnp 1993Aug18): C GIVES THE OTHER TRIAD, WHERE A GIVEN J OCCURS, AND ITS POSITION L = LINE(JB,1) IF (L.EQ.0 .OR. TABS(L)) THEN L = LINE(JB,2) LC = LCOL(JB,2) ELSE LC = LCOL(JB,1) ENDIF C CALL NEIBOR(LC,L1,L2) -- same C GIVES THE POSITIONS OF THE OTHER TWO ARGUMENTS IN THE TRIAD IF (LC.LT.2) THEN L1 = 2 L2 = 3 ELSE IF (LC.EQ.2) THEN L1 = 3 L2 = 1 ELSE L1 = 1 L2 = 2 ENDIF CALL WAY(L,L1,L2,ICH,ND) NODE=NODE+I1 TABS(L)=.TRUE. JDIAG(NODE,K3)=J23(L,LC) ARR(NODE,K3)=ARROW(L,LC) ICT=ICH*I1 IF(ICH.GT.0) GO TO 18 LP=L1 L1=L2 L2=LP 18 IF(ICT.GT.0) GO TO 5 C C CALL PHASE(L,J23,KFL2A) -- PHASE FACTOR C ARISING FROM NON-CYCLIC PERMUTATION OF ARGUMENTS IN TRIAD L: J7(J7C+1) = J23(L,1) J7(J7C+2) = J23(L,2) J7C = J7C + 3 J7(J7C) = J23(L,3) C 5 JDIAG(NODE,K1)=J23(L,L1) ARR(NODE,K1)=ARROW(L,L1) JDIAG(NODE,K2)=J23(L,L2) ARR(NODE,K2)=ARROW(L,L2) J=J23(L,L1) IAL(J)=IAL(J)+1 J=J23(L,L2) IAL(J)=IAL(J)+1 IF(ND.LT.1)GO TO 3 NTIME=NTIME+1 ILAST=MAX(NODE,ILAST) IFIRST=MIN(NODE,NBTR) NBP=IAL(IF1)+IAL(IF2) IF(NBP.GT.3) GO TO 12 IF(NTIME.GT.1) GO TO 12 IF(NBP.LE.2) GO TO 11 IF(IAL(IF1).GT.IAL(IF2)) GO TO 11 JT=JDIAG(NBTR,1) JAR=ARR(NBTR,1) JDIAG(NBTR,1)=JDIAG(NBTR,3) ARR(NBTR,1)=ARR(NBTR,3) JDIAG(NBTR,3)=JT ARR(NBTR,3)=JAR C C CALL PHASE(NBTR,JDIAG,KFL2B) -- PHASE FACTOR C ARISING FROM NON-CYCLIC PERMUTATION OF ARGUMENTS IN TRIAD L: J7(J7C+1) = JDIAG(NBTR,1) J7(J7C+2) = JDIAG(NBTR,2) J7C = J7C + 3 J7(J7C) = JDIAG(NBTR,3) C 11 NODE=NBTR I1=-1 K2=3 K3=2 GO TO 3 12 NBNODE=ILAST-IFIRST+1 NBTR=NBTR-NBNODE C C DEFINITION OF FREE ENDS AND OTHER QUANTITIES C CALL INTAB CALL PRINTJ('DIAGRM',12) RETURN END C*********************************************************************** SUBROUTINE FACTT C C - CALCULATES THE LOGS OF FACTORIALS REQUIRED BY THE RACAH C COEFFICIENT ROUTINE DRACAH C WRITTEN BY N.S. SCOTT; MXFCT CODE ADDED WE'89JUN19TH. C C IMPLICIT DOUBLE PRECISION(G,L,O,X) COMMON /FACTS/ GAM( 199) DATA MXFCT,ONE/ 199,1.0/ GAM(1)=ONE X=ONE M=MIN(MXFCT,34) DO 10 I=2,M GAM(I)=GAM(I-1)*X 10 X=X+ONE DO 20 I=1,M 20 GAM(I)=LOG(GAM(I)) M=M+1 DO 30 I=M,MXFCT GAM(I)=GAM(I-1)+LOG(X) 30 X=X+ONE RETURN END C*********************************************************************** SUBROUTINE DRACAH(I,J,K,L,M,N,RAC) C C SUBROUTINE TO CALCULATE RACAH COEFFICIENTS C THE ARGUMENTS I,J,K,L,M,N SHOULD BE TWICE THEIR ACTUAL VALUE C WORKS FOR INTEGER AND HALF-INTEGER VALUES OF ANGULAR MOMENTA. C THE ROUTINE MAKES USE OF THE GAM ARRAY, THUS SUBROUTINE FACTT C MUST BE CALLED BEFORE THIS ROUTINE IS USED. C WRITTEN BY N S SCOTT; CHECK IF...PRINT6 ADDED WE'89JUN19TH C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /FACTS/ GAM( 199) C RAC=0.0 J1=I+J+M J2=K+L+M J3=I+K+N J4=J+L+N IF(2*MAX(I,J,M).GT.J1 .OR. MOD(J1,2).NE.0) GO TO 2 IF(2*MAX(K,L,M).GT.J2 .OR. MOD(J2,2).NE.0) GO TO 2 IF(2*MAX(I,K,N).GT.J3 .OR. MOD(J3,2).NE.0) GO TO 2 IF(2*MAX(J,L,N).GT.J4 .OR. MOD(J4,2).NE.0) GO TO 2 J5=(I+J+K+L)/2 J6=(I+L+M+N)/2 J7=(J+K+M+N)/2 NUMAX=MIN(J5,J6,J7)+1 IF(NUMAX.GE. 199) GO TO 5 RAC=1.0 J1=J1/2 J2=J2/2 J3=J3/2 J4=J4/2 NUMIN=MAX(J1,J2,J3,J4)+1 IF(NUMIN.EQ.NUMAX) GO TO 4 KF=NUMIN+1 DO 3 KI=NUMAX,KF,-1 3 RAC = - (KI*(J5-KI+2)*(J6-KI+2)*(J7-KI+2))*RAC/ * ((KI-1-J1)*(KI-1-J2)*(KI-1-J3)*(KI-1-J4)) + 1.0 4 RAC = (2*MOD(J5+NUMIN,2)-1) * EXP(GAM(NUMIN+1)-GAM(NUMIN-J1) * -GAM(NUMIN-J2)-GAM(NUMIN-J3)-GAM(NUMIN-J4)-GAM(J5+2-NUMIN) * -GAM(J6+2-NUMIN)-GAM(J7+2-NUMIN) +(GAM(J1+1-I)+GAM(J1+1-J) * +GAM(J1+1-M)-GAM(J1+2)+GAM(J2+1-K)+GAM(J2+1-L)+GAM(J2+1-M) * -GAM(J2+2)+GAM(J3+1-I)+GAM(J3+1-K)+GAM(J3+1-N)-GAM(J3+2) * +GAM(J4+1-J)+GAM(J4+1-L)+GAM(J4+1-N)-GAM(J4+2))/2) * RAC 2 RETURN 5 PRINT 6, NUMAX STOP 6 FORMAT(/' STOP IN DRACAH, L.8 >',I4,' NEEDED FOR FACTORIAL ARRAY') END C*********************************************************************** SUBROUTINE GENSUM(J6C,J7C,J8C,J9C,JWC,J6,J7,J8,J9,JW,JDEL, + LDEL,SUMVAR,MP, + J6P,J7P,J8P,J9P,JWORD,NLSUM,NBJ,NB6J, + K6CP,K7CP,K8CP,K9CP,JSUM4,JSUM5,JSUM6,INV6J, + RECUP) C C -- CARRIES OUT THE SUMMATION OVER COEFFICIENTS DEFINED BY THE ARRAYS C J6,J7,J8,LDEL AND JW TO GIVE RECUP C THE ENTRY IS EITHER MADE FROM NJGRAF OR DIRECTLY, ASSUMING THAT THE C ARRAYS J6,...JW HAVE ALREADY BEEN DETERMINED BY A PREVIOUS C ENTRY TO NJGRAF AND THAT THE SUMMATION IS REQUIRED FOR ANOTHER SET C OF J VALUES DEFINED BY THE ARRAY J1. TIDIED UP MEUDON'89AOUT - WE. C C RECUP IS THE RECOUPLING COEFFICIENT C C SUBROUTINE CALLED: DRACAH C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER (LL74= 21+4, LL76= 21*3+19, LL01=50) ! from 5 PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20, + KFLS=12,KFLN=10,KFLV=40) C LOGICAL SUMVAR DIMENSION J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),JW(6,KFLW),LDEL(KFLW,2),SUMVAR(KFL1) C DIMENSION J6P(KFLV),J7P(KFLV),J8P(KFLV),J9P(KFLV),JWORD(6,KFLW), + NBJ(KFLN),NB6J(KFLN),K6CP(KFLN),K7CP(KFLN),K8CP(KFLN),K9CP(KFLN), + JSUM6(KFLS),JSUM4(KFLS,KFLW),JSUM5(KFLS,KFLW),INV6J(KFLW) C LOGICAL LDIAG,NOEL DIMENSION MAT(KFLS,KFLS),JMNP(LL01),JMXP(LL01),NOEL(KFLS), * MAXLP(KFLS),IST(6),JSUM2(KFLS),JSUM3(KFLS),JSUM8(KFLS), * JSUM(2,KFLW),JWTEST(KFLW),WSTOR(KFLW),IPAIR(2,2), * LDIAG(KFLS),J12(LL01,KFLS,KFLS), XJ1(KFL1) C LOGICAL FREE COMMON/COUPLE/M,N,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9 C DATA EPSIL/1.E-10/, MXCSVR/4/ C C FORMAT STATEMENTS USED IN GENSUM C 302 FORMAT(' SUM NR.',I3) 303 FORMAT(' NO SUMMATION. RECOUPLING COEFFICIENT =',G16.8) 304 FORMAT(' RECOUPLING COEFFICIENT=',G15.8) 305 FORMAT(6F5.1,10X,G15.8) 306 FORMAT(' NUMBER OF INDEPENDENT SUMS:',I3//' RACAH W FUNCTIONS(6JS) *'/' ARGUMENTS IN *REAL* FORMAT',18X,'VALUE') 307 FORMAT(' SUM NR.',I2,' SUM VALUE=',G15.8,' RECUP=',G15.8) 309 FORMAT(' NOT INVOLVING SUMMATION VARIABLE') 400 FORMAT(//' PRINTOUT FROM SUBROUTINE GENSUM'/ +' VALUES OF ANGULAR MOMENTA IN *REAL* FORMAT'/(14F5.1)) C C C 1. C EVALUATES ALL TERMS IN J6,J7,J8,J9,LDEL, AND JW WHICH DO NOT C INVOLVE SUMMATION. THE RESULT IS STORED IN RECUP AND IASTOR C IF(IBUG3.NE.1) GO TO 140 DO 139 I=1,M 139 XJ1(I)=(J1(I)-1.0)/2.0 PRINT 400, (XJ1(I),I=1,M) PRINT 306, NLSUM 140 MM=M+1 J1(MM)=1 C C TEST DELTA FUNCTIONS C J1(MM)=1 IF(JDEL.LE.0) GO TO 180 DO 141 I=1,JDEL I1=LDEL(I,1) I2=LDEL(I,2) IF(I1.LE.MM.AND.I2.LE.MM)GO TO 143 IF(I1.GT.MM)J1(I1)=J1(I2) IF(I2.GT.MM)J1(I2)=J1(I1) GO TO 141 143 IF(J1(I1).EQ.J1(I2)) GO TO 141 RECUP=0.0 GO TO 110 141 CONTINUE 180 RECUP=1.0 IF(JWC.EQ.0) GO TO 8 C C MULTIPLY RECUP BY ALL RACAH COEFFICIENTS WHICH DO NOT INVOLVE A C SUMMATION C IF(IBUG3.EQ.1) PRINT 309 1 DO 7 I=1,JWC IF(INV6J(I).GT.0)GO TO 7 DO 3 J=1,6 I1=JW(J,I) 3 IST(J) = J1(I1) - 1 CALL DRACAH(IST(1),IST(2),IST(3),IST(4),IST(5),IST(6),X1) IF(IBUG3.EQ.1) PRINT 305,(XJ1(K),K=1,6),X1 RECUP = RECUP*X1 7 CONTINUE 8 SQR=1.0 IF(J6C.EQ.0) GO TO 10 C 9 DO 12 I=1,J6C I1=J6(I) 12 SQR=SQR*J1(I1) 10 SPR=1.0 IF(J9C.EQ.0)GO TO 13 DO 144 I=1,J9C I1=J9(I) 144 SPR=SPR*J1(I1) 13 RECUP=RECUP*SQRT(SQR/SPR) IF(ABS(RECUP).LT.EPSIL)GO TO 145 IASTOR = 0 IF(J7C.EQ.0) GO TO 18 C 14 DO 17 I=1,J7C I1=J7(I) 17 IASTOR = IASTOR + J1(I1) -1 18 IF(J8C.EQ.0) GO TO 23 C 19 DO 22 I=1,J8C I1=J8(I) 22 IASTOR = IASTOR +2*(J1(I1)-1) 23 IF(NLSUM.GT.0)GO TO 24 C IASTOR=IASTOR/2 C NO SUMMATION INVOLVED. END OF COMPUTATION C STOR1=1.0 STOR=1.0 IF(MOD(IASTOR,2).NE.0) RECUP=-RECUP IF(IBUG3.EQ.1) PRINT 303,RECUP GO TO 110 C C C C 2. C EVALUATION OF THE PART INVOLVING SUMMATIONS C C 24 NFS=0 JWR=0 J6F=0 J7F=0 J8F=0 J9F=0 NPS=0 25 NPS=NPS+1 IF(IBUG3.EQ.1) PRINT 302,NPS C C 2.0 LOOP ON THE DISCONNECTED SUMMATIONS C C IAS=0 NSUM=NBJ(NPS)-NFS JWRD=NB6J(NPS)-JWR J6CP=K6CP(NPS) J7CP=K7CP(NPS) J8CP=K8CP(NPS) J9CP=K9CP(NPS) C C 2.1 THE RANGE OF VALUES OF EACH SUMMATION VARIABLE IS C DEFINED BY ESTABLISHING A MATRIX OF THE LINKS BETWEEN C VARIABLES. MAT(I,J) CONTAINS C I=J ,NUMBER OF POSSIBLE VALUES OF I DUE TO TRIANGULAR C RELATIONS WITH NON-VARIABLES, I.E. CONSTANTS C I.GT.J NUMBER OF LINKS BETWEEN I AND J THROUGH CONSTANTS C I.LT.J VALUE OF THE CONSTANT, IF THE ABOVE IS 1. IF NOT, C THESE VALUES ARE STORED IN J12(L,I,J) WHERE THERE C IS ROOM FOR MXCSVR SUCH VALUES (L.LE.4) C C DO 52 J=1,NSUM DO 52 I=1,NSUM 52 MAT(I,J)=0 DO 66 I1=1,NSUM I1T=I1+NFS I2=JSUM6(I1T) DO 65 I3=1,I2 I=JSUM5(I1T,I3) J=JSUM4(I1T,I3) GO TO (54,55,56,57,58,59),J C C THE ROWS OF THE IPAIR ARRAYS GIVE LIMITS OF SUMMATION IMPOSED C BY THE TRIANGULAR CONDITION C 54 IPAIR(1,1) = JWORD(2,I) IPAIR(1,2) = JWORD(5,I) IPAIR(2,1) = JWORD(3,I) IPAIR(2,2) = JWORD(6,I) GO TO 60 55 IPAIR(1,1) = JWORD(1,I) IPAIR(1,2) = JWORD(5,I) IPAIR(2,1) = JWORD(4,I) IPAIR(2,2) = JWORD(6,I) GO TO 60 56 IPAIR(1,1) = JWORD(1,I) IPAIR(1,2) = JWORD(6,I) IPAIR(2,1) = JWORD(4,I) IPAIR(2,2) = JWORD(5,I) GO TO 60 57 IPAIR(1,1) = JWORD(2,I) IPAIR(1,2) = JWORD(6,I) IPAIR(2,1) = JWORD(3,I) IPAIR(2,2) = JWORD(5,I) GO TO 60 58 IPAIR(1,1) = JWORD(1,I) IPAIR(1,2) = JWORD(2,I) IPAIR(2,1) = JWORD(3,I) IPAIR(2,2) = JWORD(4,I) GO TO 60 59 IPAIR(1,1) = JWORD(1,I) IPAIR(1,2) = JWORD(3,I) IPAIR(2,1) = JWORD(2,I) IPAIR(2,2) = JWORD(4,I) 60 DO 63 I4=1,2 KM=0 DO 62 I5=1,2 IF(IPAIR(I4,I5).GT.MP)KM=KM+1 62 CONTINUE JJ1=IPAIR(I4,1) JJ2=IPAIR(I4,2) IF(KM-1)64,67,63 C C ONE VARIABLE LINKED TO TWO CONSTANTS. FIX THE DIAGONAL MAT(I,I) C 64 JT1=J1(JJ1)-1 JT2=J1(JJ2)-1 JMIN=ABS(JT1-JT2) JMAX=JT1+JT2 IF(MAT(I1,I1)-1) 68,63,69 C C FIRST TIME C 68 MAT(I1,I1)=(JMAX-JMIN)/2+1 JSUM(1,I1)=JMIN JSUM(2,I1)=JMAX GO TO 63 C C IF THERE ARE SEVERAL COUPLES OF CONSTANTS, TAKE THE MORE C STRINGENT COMBINATION C 69 JMIN=MAX(JMIN,JSUM(1,I1)) JMAX=MIN(JMAX,JSUM(2,I1)) IF(JMAX.GE.JMIN) GO TO 70 RECUP=0.0 GO TO 110 70 JSUM(1,I1)=JMIN JSUM(2,I1)=JMAX MAT(I1,I1)=(JMAX-JMIN)/2+1 GO TO 63 C C ONE VARIABLE LINKED TO ONE CONSTANT AND ONE VARIABLE NON DIAGONAL C ELEMENT C 67 JT1=MIN(JJ1,JJ2) JT2=MAX(JJ1,JJ2)-MP IF(JT2.GT.I1) GO TO 63 JT4=J1(JT1)-1 K=MAT(I1,JT2) IF(K.EQ.0) GO TO 107 DO 71 LL=1,K IF(JT4.EQ.J12(LL,JT2,I1)) GO TO 63 71 CONTINUE 107 K=K+1 IF(K.GT.MXCSVR) GO TO 63 MAT(I1,JT2)=K J12(K,JT2,I1)=JT4 63 CONTINUE 65 CONTINUE 66 CONTINUE C C REDUCE THE DIAGONAL ELEMENTS BY TAKING INTO ACCOUNT THE NON C DIAGONAL ELEMENTS, AND KEEP THE LATTER ONLY IF NEEDED C LOOP STRUCTURE MODIFIED AVOIDING JUMPS BACK IN WE'89AOUT24 C 150 ICHAN=0 DO 74 I=1,NSUM NOEL(I)=.TRUE. I2=1 I1=I-1 IF(I1.NE.0) GO TO 72 160 IF(I.EQ.NSUM) GO TO 74 I1=NSUM I2=I+1 72 DO 73 J=I2,I1 IF(MAT(J,J).EQ.0) GO TO 73 IF(I.EQ.1) GO TO 170 IF(MAT(I,J).EQ.0) GO TO 73 IK1=I IK2=J NOEL(I)=.FALSE. GO TO 201 170 IF(MAT(J,I).EQ.0) GO TO 73 IK1=J IK2=I C 201 JMIN1=0 JMAX1=1000 K=MAT(IK1,IK2) if(K.gt.LL01) then ! Lyngby'05Jul07 print"(/' stglib stopping in GENSUM: LL01 .ge.',I3,' needed'/)",K stop endif DO 203 L1=1,K L3=MAT(J,J) JJ1=JSUM(1,J) JND=J12(L1,IK2,IK1) JMIN=1000 JMAX=0 JMNP(L1)=0 JMXP(L1)=1000 DO 204 L2=1,L3 JMN=ABS(JND-JJ1) JMX=JND+JJ1 JMIN=MIN(JMN,JMIN) JMAX=MAX(JMX,JMAX) JMNP(L1)=MAX(JMN,JMNP(L1)) JMXP(L1)=MIN(JMX,JMXP(L1)) 204 JJ1=JJ1+2 JMIN1=MAX(JMIN1,JMIN) 203 JMAX1=MIN(JMAX1,JMAX) IF(MAT(I,I).NE.0) GO TO 205 JSUM(1,I)=JMIN1 JSUM(2,I)=JMAX1 MAT(I,I)=(JMAX1-JMIN1)/2+1 ICHAN=ICHAN+1 GO TO 206 205 IF(JSUM(1,I).GE.JMIN1) GO TO 211 JSUM(1,I)=JMIN1 ICHAN=ICHAN+1 211 IF(JSUM(2,I).LE.JMAX1) GO TO 206 JSUM(2,I)=JMAX1 ICHAN=ICHAN+1 206 K1=0 DO 207 L1=1,K IF(JMNP(L1).LE.JSUM(1,I).AND.JMXP(L1).GE.JSUM(2,I))GO TO 207 K1=K1+1 J12(K1,IK2,IK1)=J12(L1,IK2,IK1) 207 CONTINUE IF(K1.EQ.K)GO TO 208 MAT(IK1,IK2)=K1 ICHAN=ICHAN+1 208 MAT(IK2,IK1)=J12(1,IK2,IK1) C OUT GO TO JKM(171,172) -- SECTION REWRITTEN AT CECAM MEUDON 1989. W.E. C 73 CONTINUE IF(I1.NE.NSUM) GO TO 160 74 CONTINUE IF(ICHAN.NE.0) GO TO 150 C C 2.3 CARRY OUT THE SUMMATIONS C DO 230 I=1,NSUM LDIAG(I)=MAT(I,I).EQ.1 230 JSUM3(I)=1 DO 231 I=1,JWRD 231 JWTEST(I)=1 STOR=0.0 STOR1=1.0 NOLP=0 IP=1 240 NOLP=NOLP+1 C 2.3.1 C FIND THE RANGE OF JSUM2(NOLP) C NOLP IS THE INDEX OF THE SUMMATION VARIABLE C JMIN=JSUM(1,NOLP) JMAX=JSUM(2,NOLP) IF(NOEL(NOLP))GO TO 241 NO1=NOLP-1 DO 242 NJ=1,NO1 IF(MAT(NOLP,NJ)-1)242,243,244 243 JJ1=MAT(NJ,NOLP) JJ2=JSUM2(NJ) JMIN=MAX(JMIN,ABS(JJ2-JJ1)) JMAX=MIN(JMAX,JJ1+JJ2) GO TO 242 244 K=MAT(NOLP,NJ) JJ2=JSUM2(NJ) DO 245 I=1,K JJ1=J12(I,NJ,NOLP) JMIN=MAX(JMIN,ABS(JJ2-JJ1)) 245 JMAX=MIN(JMAX,JJ1+JJ2) 242 CONTINUE 241 JSUM2(NOLP)=JMIN MAXLP(NOLP)=JMAX IF(LDIAG(NOLP))JSUM3(NOLP)=0 IF(NOLP.LT.NSUM)GO TO 240 DO 260 JJ=JMIN,JMAX,2 JSUM2(NSUM)=JJ C C 2.3.2 DETERMINE WHICH RACAH COEFFICIENTS NEED RE-EVALUATING AND C SET JWTEST APPROPRIATELY C 111 DO 114 J=IP,NSUM IF(JSUM3(J).LE.0) GO TO 114 I2=JSUM6(J) DO 113 I1=1,I2 I3=JSUM5(J,I1) 113 JWTEST(I3)=1 114 CONTINUE DO 98 J=1,JWRD IF(JWTEST(J).EQ.0)GO TO 98 JWJ=J+JWR DO 90 I=1,6 I1=JWORD(I,JWJ) IF(I1.GT.MP) GO TO 89 IST(I) = J1(I1) - 1 GO TO 90 89 I1=I1-MP-NFS IST(I) = JSUM2(I1) 90 CONTINUE CALL DRACAH(IST(1),IST(2),IST(3),IST(4),IST(5),IST(6),X1) WSTOR(J)=X1 IF(IBUG3.NE.1) GO TO 98 DO 99 I=1,6 99 XJ1(I)=IST(I)/2.0 PRINT 305,(XJ1(I),I=1,6),X1 98 CONTINUE C C 2.3.3 C FORM PRODUCT OF RACAH COEFFICIENTS,(2J+1) FACTORS AND (-1) C FACTORS IN STOR1 C DO 126 I=1,JWRD 126 STOR1 = STOR1*WSTOR(I) C C IASTOR CONTAINS THE POWER OF (-1) WHICH IS COMMON TO ALL TERMS C 123 IX2 = 0 IJ6CP=1 IF(J6CP.EQ.J6F) GO TO 146 JB=J6F+1 DO 128 I=JB,J6CP I1=J6P(I)-NFS 128 IJ6CP=IJ6CP*(JSUM2(I1)+1) 146 IF(J9CP.EQ.J9F) GO TO 129 JB=J9F+1 DO 147 I=JB,J9CP I1=J9P(I)-NFS 147 IJ6CP=IJ6CP/(JSUM2(I1)+1) 129 STOR1 = STOR1*SQRT(REAL(IJ6CP)) IF(J7CP.EQ.J7F) GO TO 132 130 JB=J7F+1 DO 131 I=JB,J7CP I1=J7P(I)-NFS 131 IX2 = IX2 + JSUM2(I1) 132 IF(J8CP.EQ.J8F) GO TO 135 JB=J8F+1 DO 134 I=JB,J8CP I1=J8P(I)-NFS 134 IX2 = IX2 + 2*(JSUM2(I1)) 135 IF(MOD(IX2,2).NE.1)GO TO 137 IAS=-1 IX2=IX2+1 137 IX2 = IX2/2 C C 2.3.4 C ADD TERM INTO STOR AND RESET STOR1 TO 1 READY FOR NEXT TERM C IF (MOD(IX2,2) .EQ. 1) STOR1 = -STOR1 STOR = STOR + STOR1 STOR1=1.0 NSUM1 =NSUM-1 IF(NSUM1.EQ.0)GO TO 260 DO 261 IK=1,NSUM1 261 JSUM3(IK)=0 DO 262 IK=1,JWRD 262 JWTEST(IK)=0 260 CONTINUE 250 NOLP=NOLP-1 IF(NOLP.EQ.0)GO TO 108 IF(LDIAG(NOLP))GO TO 250 JSUM3(NOLP)=1 JSUM2(NOLP)=JSUM2(NOLP)+2 IF(JSUM2(NOLP).GT.MAXLP(NOLP))GO TO 250 IP=NOLP C CC 2.3.5 PROCEED TO NEXT VARIABLE C GO TO 240 108 RECUP=RECUP*STOR IF(IBUG3.EQ.1) PRINT 307,NPS,STOR,RECUP IF(ABS(RECUP).LT.EPSIL)GO TO 145 JWR=JWRD + JWR NFS=NSUM + NFS J6F=J6CP J7F=J7CP J8F=J8CP J9F=J9CP IASTOR=IASTOR+IAS C C 2.4 PROCEED TO NEXT SUM C IF(NPS.LT.NLSUM)GO TO 25 IASTOR=IASTOR/2 IF(MOD(IASTOR,2).NE.0)RECUP=-RECUP 136 IF(IBUG3.EQ.1) PRINT 304,RECUP 110 RETURN C NO SUMMATIONS. CHECK THAT THERE ARE NO INCONSISTENCIES. THEN C MULTIPLY BY (-1) FACTOR AND EXIT C 145 RECUP=0.0 GO TO 110 END C*********************************************************************** SUBROUTINE INTAB C C - IS CALLED AT THE END OF DIAGRM, FIXES THE ARRAYS IH AND IL -- C SO TO SPEAK HARDWARE AND LOGICAL ADDRESSES OF TRIADS IN JDIAG; C ALSO DETERMINES THE NUMBER OF FREE ENDS NFREE AND THEIR C LOCATION ITFREE. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2) C INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C COMMON/BUILD/IAL(KFL2B),IF1,IF2,NODE C C LOGICAL FREE COMMON/COUPLE/M,N,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) C DO 1 I=1,M 1 IAL(I)=1 DO 3 I=IFIRST,ILAST J=JDIAG(I,1) K=IAL(J) C!! if (K.gt.2) then ! C(-) trouble zhl 2005 Feb-Mar C!! print "(/' M,first,last,I,J,K =',6I5)", m,ifirst,ilast,i,j,k C!! stop ! cured in ORDTRI '05May31 C!! endif TAB1(J,K)=I 3 IAL(J)=K+1 IFR=IFIRST-1 DO 4 I=IFIRST,ILAST IT=I-IFR IL(I)=IT 4 IH(IT)=I J=JDIAG(IFIRST,3) K=IAL(J) IF(K.LE.1) GO TO 6 TAB1(J,2)=TAB1(J,1) 6 TAB1(J,1)=IFIRST IAL(J)=3 J=JDIAG(ILAST,2) TAB1(J,2)=ILAST IAL(J)=3 NFREE=0 DO 7 I=IFIRST,ILAST J=JDIAG(I,1) IF(IAL(J).EQ.3)GO TO 7 NFREE=NFREE+1 ITT=ILAST+NFREE TAB1(J,2)=ITT IL(ITT)=NFREE*1000 ITFREE(NFREE)=I 7 CONTINUE RETURN END C*********************************************************************** SUBROUTINE LOLPOP(FAIL) C C - REDUCES A LOOP WITH ONE LINE AND ONE NODE IN THE FLAT GRAPH. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C LOGICAL FAIL C INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP LOGICAL FREE COMMON/COUPLE/M,N,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) C DIMENSION KP(3),KS(3) C CHARACTER*6 NAMSUB COMMON /NAM/NAMSUB C DATA KP/2,3,1/, KS/0,1,-1/ C NAMSUB = 'LOLPOP' I1=NPOINT(1) K3=2 IF(I1.EQ.ILAST)K3=3 L=JDIAG(I1,K3) CALL DELTA(L,M,FAIL) IF(FAIL)GO TO 2 K=KP(K3) IF(ARR(I1,K).LT.0) THEN C CALL PHASE2(JDIAG(I1,K)) ! ADDS FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = JDIAG(I1,K) ENDIF K1=KS(K3) IL1=IL(I1)+K1 I2=IH(IL1) L1=JDIAG(I2,1) CALL DELTA(L1,JDIAG(I2,K3),FAIL) IF(FAIL)GO TO 2 IF(ARR(I2,K3).EQ.K1) THEN C CALL PHASE2(L1) ! ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = L1 ENDIF IL2=IL(I2)+K1 I3=IH(IL2) K2=K3+K1 JDIAG(I3,K2)=L1 ARR(I3,K2)=ARR(I2,1) J9C=J9C+1 J9(J9C)=L1 J6C=J6C+1 J6(J6C)=JDIAG(I1,1) IF(K3.EQ.3)GO TO 2 DO 1 I=3,NBNODE IT=IH(I) ILP=I-2 IL(IT)=ILP 1 IH(ILP)=IT 2 RETURN END C*********************************************************************** C SUBROUTINE NEIBOR(LC,L1,L2) C GIVES THE POSITIONS OF THE OTHER TWO ARGUMENTS IN THE TRIAD C*********************************************************************** SUBROUTINE ORDTRI C C -- ORDERS THE TRIADS WHICH WERE LEFT WITH FREE ENDS AS C CONSEQUENCE OF CUTTING, SO THAT THE NEW GRAPH WILL START THERE C C '05May31 -- Corrections by SMW Bailey, DAMTP QUB, 29 Apr 1998: C the value of NBT1 has been changed at two places to put a triad C in the correct position, and DO loop 21 has been moved to avoid C overwriting a triad at position NM -- cures zhl's C(-) trouble! C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C INTEGER ARROW LOGICAL TABS COMMON/TREE/J23(KFL2A,3),ARROW(KFL2A,3),LINE(KFL1,2),LCOL(KFL1,2), + TABS(KFL2A),NBTR C INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC C COMMON/BUILD/IAL(KFL2B),IF1,IF2,NODE COMMON/KEEP/JKP(2,3),JARR(2,3),IT2,IT3,IT5 ! from CUTNL C * print"(/' ORDTRI: NC,NFIN,NFREE,NBTR =',4I5/29X,5I5,' =BNODE -> 2' * *)", NC,NFIN,NFREE,NBTR, IT2,IT3,IT5,ITFREE(1),NBNODE !'05Apr-May DO 10 I=1,KFL2B ! MPWE95JAN8ZHL: IF(J1QN1(J,1).NE.J1QN2(J,1)) GO TO 182 IF(J1QN1(J,2).NE.J1QN2(J,2)) GO TO 182 IF(J1QN1(J,3).NE.J1QN2(J,3)) GO TO 182 156 CONTINUE IF(K.EQ.3) GO TO 157 J1(M5)=2*LJ(IRHO)+1 J1(M10)=2*LJ(ISIG)+1 J1(MN1)=2*KA+1 IF(ISPIN.EQ.1) J1(MN1)=1 J1(M3)=JMEM2 J1(M4)=JMEM5 IF(IRHO.EQ.ISIG) GO TO 158 J1(M3)=J1QN1(IRHO,K) J1(M4)=J1QN2(ISIG,K) GO TO 158 157 J1(M5)=2 J1(M10)=2 J1(MN1)=1 IF(ISPIN.NE.0) J1(MN1)=2*KA+1 J1(M3)=JMEM3 J1(M4)=JMEM6 IF(IRHO.EQ.ISIG) GO TO 158 J1(M3)=J1QN1(IRHO,K) J1(M4)=J1QN2(ISIG,K) 158 DO 161 J=1,IHSH IF(IRHO.NE.J) GO TO 160 J1(J)=J1QN2(IRHO,K ) GO TO 161 160 J1(J)=J1QN1(J,K) 161 CONTINUE IF(IHSH.EQ.1) GO TO 197 DO 162 J=M6,M7 162 J1(J)=J1QN1(J,K) DO 163 J=M9,M11 163 J1(J)=J1QN2(J-M12,K) C C PRINT OUT THE J1,J2 AND J3 ARRAYS C C 10 STATEMENTS REPLACED BY THE FOLLOWING 4 -- WE'90MAR15: 197 IF(K.GE.3) GO TO 164 IF(.NOT.NBUG6) GO TO 164 WRITE(IWRITE,210) (J1(J),J=1,MN1) WRITE(IWRITE,212) ((J2(I,J),J=1,3),(J3(I,J),J=1,3),I=1,IHSHP1) C C EVALUATE ORBITAL AND SPIN RECOUPLING COEFFICIENTS C NJGRAF IS CALLED EACH TIME SO SET THE ELEMENTS OF THE FREE ARRAY C TO .FALSE. C 164 DO 500 I=1,MN1 500 FREE(I)=.FALSE. C CALL NJGRAF(RECUP,FAIL) C RECUPS=RECUPS*RECUP IF(K.GE.3) GO TO 170 IF(NBUG6) WRITE(IWRITE,213) RECUP 170 K=K+1 DO 169 KK=1,3 DO 169 J=1,IHSHP1 J2(J,KK)=J2STO(J,KK) 169 J3(J,KK)=J3STO(J,KK) IF(K.EQ.3) GO TO 155 IF(NBUG6) WRITE(IWRITE,214) RECUP C C FIRST FRACTIONAL PARENTAGE COEFFICIENT C LIJ=LJ(IRHO) COEFP=1.0 IF(LIJ.EQ.0) GO TO 272 N=NOSH1(IRHO) IV1=JMEM1 IL1=(JMEM2-1)/2 IS1= JMEM3 IV2=J1QN2(IRHO,1) IL2=(J1QN2(IRHO,2)-1 )/2 IS2=J1QN2(IRHO,3) CALL CFP(LIJ,N,IV1,IL1,IS1,IV2,IL2,IS2,COEFP) RECUPS=RECUPS*COEFP 272 IF(IRHO.EQ.ISIG) GO TO 173 IF(ABS(RECUPS).LT.1.0E-5) GO TO 183 C C SECOND FRACTIONAL PARENTAGE COEFFICIENT C 173 LIJ=LJ(ISIG) COEFP=1.0 IF(LIJ.LE.0) GO TO 176 N=NOSH2(ISIG) IV1=JMEM4 IL1=(JMEM5-1)/2 IS1=JMEM6 IV2=J1QN1(ISIG,1) IL2=(J1QN1(ISIG,2)-1)/2 IS2=J1QN1(ISIG,3) CALL CFP(LIJ,N,IV1,IL1,IS1,IV2,IL2,IS2,COEFP) 176 RECUPS=RECUPS*COEFP IF(ABS(RECUPS).LT.1.0E-8.AND.IRHO.NE.ISIG) GO TO 183 C C PERMUTATION FACTOR C 175 N=0 IF(IRHO-ISIG) 177,185,179 177 JRHO = IRHO+1 DO 178 J=JRHO,ISIG 178 N = NOSH1(J)+N GO TO 181 179 JSIG = ISIG+1 DO 180 J=JSIG,IRHO 180 N = NOSH2(J)+N C C MULTIPLICATIVE FACTOR C 181 VALML = SQRT(REAL(NOSH1(IRHO)*NOSH2(ISIG)))*(1-MOD(N,2)*2)*RECUPS GO TO 184 182 RECUPS=0.0 IF(IRHO.EQ.ISIG) GO TO 185 183 VALML = 0.0 184 RML = RML+VALML C RESULT STORED IN VSHELL IF(NTOT.EQ.0) NTOT=1 VSHELL(NTOT)=RML*SQRT(AJF) GO TO 190 185 VALUML=RECUPS IF(IBUG6.NE.0) WRITE(IWRITE,219) JJ1,VALUML RPL = RPL+VALUML 186 JJ1=JJ1+1 IF(JJ1.LE.KK1) GO TO 152 C END OF IRHO.EQ.ISIG LOOP J1QN1(IRHO,1)=JMEM1 J1QN1(IRHO,2)=JMEM2 J1QN1(IRHO,3)=JMEM3 J1QN2(ISIG,1)=JMEM4 J1QN2(ISIG,2)=JMEM5 J1QN2(ISIG,3)=JMEM6 ANL=NOSH1(IRHO)*RPL C C RESULTS STORED IN VSHELL C IF(NTOT.EQ.0) NTOT=1 VSHELL(NTOT)=ANL*SQRT(AJF) IF(NBUG6) WRITE(IWRITE,215) IRHO,ANL IRHO=IRHO+1 ISIG=ISIG+1 RPL=0.0 IF(IRHO.LE.IHSH) GO TO 109 190 IF(NBUG6) WRITE(IWRITE,217) (VSHELL(N),N=1,NTOT) 192 RETURN END C*********************************************************************** SUBROUTINE SETUPE(JA,JB,NJCOMP,LJCOMP) C PARAMETER (LL42= 15*2+1, LL43= 21*2+3, LL61= 15+1, LL75= 21+2) DIMENSION NJCOMP(*),LJCOMP(*) ! '05Feb8: not (23), perhaps (LL75) COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH(LL75,2), 1 J1QN(LL43,3,2),IJFUL(LL75) COMMON/MSTATE/MCFG,MOCCSH(2500),MOCORB(LL61,2500), 1 MELCSH(LL61,2500),M1QNRD(LL42,3,2500),KCFG, 2 KOCCSH(2500),KOCORB(LL61,2500),KELCSH(LL61,2500), 3 K1QNRD(LL42,3,2500),MAXOR C C NOTICE THE DIFFERENT NAMES IN THE COMMON BLOCK MEDEFN - WE C STORE NOSH1(I=1,10) AS NOSH((I=1,10),1) AND NOSH2(I=1,10) AS C NOSH((I=1,10),2) AND USE THE FACT THAT NOSH1 AND NOSH2 WILL THEN C BE EQUIVALENT TO THE SINGLE 2-DIMENSIONAL ARRAY NOSH. SIMILARLY C FOR J1QN C C === GENERATES THE ARRAYS NJ,LJ - DEFINING THE QUANTUM NUMBERS OF THE C SHELLS, NOSH - DEFINING THE OCCUPATION OF THE SHELLS, J1QN - C DEFINING THE COUPLING OF THE SHELLS, FOR EACH OF THE TWO C CONFIGURATIONS CONSIDERED. ONLY THOSE SHELLS OCCURRING IN AT C LEAST ONE CONFIGURATION ARE INCLUDED. C AT LEAST TWO SHELLS MUST BE CONSIDERED OCCUPIED. C THUS (1S)**2 HELIUM MUST BE TREATED AS E.G. (1S)**2(2S)**0 C THE SIZE OF THE ARRAYS HERE CALCULATED IS ARRANGED TO BE NO C GREATER THAN IS NECESSARY TO INCLUDE ALL ORBITALS WHICH ARE C DEEMED TO BE OCCUPIED IN EITHER OR BOTH OF THE CONFIGURATIONS C JA,JB C C --- INITIALIZE BASIC QUANTITIES - (I1+1) RUNS OVER 1,MAXORB, IHSH IS C THE CURRENT VALUE OF THE HIGHEST OCCUPIED SHELL YET CONSIDERED, C WHILE I2HSH=2*IHSH-1 C I1=0 IHSH=0 I2HSH=-1 IA=MOCCSH(JA) IB=KOCCSH(JB) C C --- TEST ON WHETHER LIMIT OF I1 HAS BEEN REACHED C 1 IF(I1.GE.MAXOR ) GO TO 100 C C --- INCREASE BASIC QUANTITIES C I1=I1+1 I3=IHSH+1 I5=I2HSH+I3 C C --- IS THE SHELL I1 OCCUPIED IN JA C DO 2 J=1,IA IF(I1.EQ.MOCORB(J,JA)) GO TO 3 2 CONTINUE NA=1 GO TO 4 3 NA=2 J1=J C C --- IS THE SHELL I1 OCCUPIED IN JB C 4 DO 5 J=1,IB IF(I1.EQ.KOCORB(J,JB)) GO TO 6 5 CONTINUE NB=1 GO TO 7 6 NB=2 J2=J C C IF THE SHELL I1 IS NOT OCCUPIED IN EITHER JA OR JB, IGNORE THE C SHELL, DO NOT INCREASE IHSH, AND CONSIDER NEXT SHELL BY INCREASING C I1 C 7 IF(NA.GT.1) GO TO 9 IF(NB.LE.1) GO TO 1 C C --- IF THE SHELL I1 IS OCCUPIED IN EITHER JA OR JB - C (1) IF IHSH.GT.1, THEN ALREADY AT LEAST TWO SHELLS AND THE C RESULTING COUPLINGS HAVE BEEN STORED. WE MUST THUS MAKE ROOM FOR C THE QUANTUM NUMBERS OF THIS NEW SHELL BETWEEN THE QUANTUM NUMBERS C OF THE PREVIOUS SHELLS AND THE QUANTUM NUMBERS OF THE INTERMEDIATE C COUPLINGS OF THE CONFIGURATIONS. THUS THE LATTER SET ARE =MOVED C ALONG= TO MAKE ROOM FOR THE NEW SHELL C (2) IF IHSH.LE.1, THERE ARE NO INTERMEDIATE COUPLING QUANTUM C NUMBERS, AND SO THERE IS NOTHING TO MOVE C 9 IF(IHSH.LE.1) GO TO 11 DO 12 I=1,2 DO 13 J=I3,I2HSH I4=I5-J DO 14 K=1,3 14 J1QN(I4+1,K,I)=J1QN(I4,K,I) 13 CONTINUE 12 CONTINUE 11 IHSH=I3 I2HSH=I2HSH+2 NC=NA I=1 IC=J1 JC=JA C C --- FIRST CONSIDER THE L.H.S. (I=1) OF THE MATRIX ELEMENT. NC=1 MEANS C UNOCCUPIED, REPRESENTED BY A DUMMY SINGLET S SHELL, AND THE C ADDITIONAL SET OF COUPLING QUANTUM NUMBERS WILL BE THE SAME AS THE C LAST SET OF COUPLING QUANTUM NUMBERS ALREADY OBTAINED. C NC=2 MEANS OCCUPIED. THEN ALL THE NEW QUANTUM NUMBERS (BOTH FOR C THE SHELL AND FOR THE COUPLING OF THIS SHELL TO THE RESULTANT OF C THE PREVIOUS ONES) ARE DEFINED IN THE CORRESPONDING J1QNRD ARRAY. C NOSH - THE NUMBER OF ELECTRONS IN THIS SHELL, IS DEFINED BY THE C APPROPRIATE ENTRY IN NELCSH. THE R.H.S. IS THEN CONSIDERED C SIMILARLY (I=2) C 25 IF(NC.GT.1) GO TO 16 NOSH(IHSH,I)=0 J1QN(IHSH,1,I)=0 J1QN(IHSH,2,I)=1 J1QN(IHSH,3,I)=1 IF(IHSH-2) 22,18,19 18 J1QN(3,1,I)=0 J1QN(3,2,I)=J1QN(1,2,I) J1QN(3,3,I)=J1QN(1,3,I) GO TO 22 19 DO 27 K=1,3 27 J1QN(I2HSH,K,I)=J1QN(I2HSH-1,K,I) GO TO 22 16 IF(I.GE.2) GO TO 38 NOSH(IHSH,I)=MELCSH(IC,JC) DO 20 K=1,3 J1QN(IHSH,K,I)=M1QNRD(IC,K,JC) C C IS THIS THE FIRST OCCUPIED SHELL OF EITHER CONFIGURATION. IF SO, C THEN THERE ARE NO INTERMEDIATE COUPLINGS TO CONSIDER AT THIS STAGE C IF(IHSH.LE.1) GO TO 20 C C IS THIS THE FIRST OCCUPIED SHELL OF THIS CONFIGURATION, THOUGH NOT C THE FIRST OF THE OTHER CONFIGURATION. IF SO, THE INTERMEDIATE C COUPLING FORMED HAS THE SAME L,S VALUES AS THIS OCCUPIED SHELL, C SINCE WE COUPLE THE SHELL TO A DUMMY SINGLET S. C IF(IC.GT.1) GO TO 29 I2=1 IF(K.NE.1) GO TO 28 C C SENIORITY SET (ARBITRARILY) ZERO FOR INTERMEDIATE COUPLING C J1QN(I2HSH,1,I)=0 GO TO 20 29 I2=MOCCSH(JC)+IC-1 28 J1QN(I2HSH,K,I)=M1QNRD(I2,K,JC) 20 CONTINUE GO TO 22 38 NOSH(IHSH,I)=KELCSH(ICE,JCE) DO 30 K=1,3 J1QN(IHSH,K,I)=K1QNRD(ICE,K,JCE) IF(IHSH.LE.1) GO TO 30 C REPLACE CARDS 799,800 BY THE FOLLOWING CARDS.CARDS WERE OMITTED C C IS THIS THE FIRST OCCUPIED SHELL OF THIS CONFIGURATION, THOUGH NOT C THE FIRST OF THE OTHER CONFIGURATION. IF SO, THE INTERMEDIATE C COUPLING FORMED HAS THE SAME L,S VALUES AS THIS OCCUPIED SHELL, C SINCE WE COUPLE THE SHELL TO A DUMMY SINGLET S. C IF(ICE.GT.1) GO TO 31 I2=1 IF(K.NE.1) GO TO 32 J1QN(I2HSH,1,I)=0 GO TO 30 31 I2=KOCCSH(JCE)+ICE-1 32 J1QN(I2HSH,K,I)=K1QNRD(I2,K,JCE) 30 CONTINUE 22 IF(I.GE.2) GO TO 24 NC=NB I=2 ICE=J2 JCE=JB GO TO 25 C C --- SET THE NJ AND LJ VALUES OF THE OCCUPIED SHELLS C 24 NJ(IHSH)=NJCOMP(I1) LJ(IHSH)=LJCOMP(I1) C C --- RETURN TO 1 TO SEE IF MAXORB HAS BEEN REACHED C GO TO 1 100 RETURN END C*********************************************************************** SUBROUTINE SEARCH(FIND) C C - LOCATES CIRCUITS OR LOOPS OF ORDER NC. NPOINT(NC) ARE THE C INDICES OF THE POINTS(TRIADS) PERTAINING TO THE FIRST C SUCH LOOP FOUND. C NPART IS THE NUMBER OF SEPARATE PARTS(GROUPS OF CONTIGUOUS POINTS) C ON THE AXIS OF THE FLAT GRAPH. IPARTS IS THE NUMBER OF POINTS IN C THE SMALLEST PART, IPARTL IS THE NUMBER OF POINTS IN THE LARGEST. C C THIS SUBROUTINE FINDS ALL THE POSSIBLE LOOPS OF ORDER 3 AND 4. FOR C NC.GE.5 IT LOOKS FOR ONLY THOSE WHO ARE PARTITIONED IN NPART.LE.2, C WHICH CAN EVENTUALLY REDUCE TO A LOOP OF ORDER 4 WITHOUT BREAKING C THE BASIC STRUCTURE OF THE FLAT GRAPH. ICROSS=-1 IF LINES CROSS C INTEGER ARR,TAB1 LOGICAL FIND, SUMVAR C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2) PARAMETER (KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C COMMON /ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), A J9(KFL9),JW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC COMMON/INFORM/IRD,IPD,IXX C C C INITIALIZATION C FIND = .FALSE. NCM1 = NC-1 NCM = NC-2 ICROSS = 0 C C FIRST TREAT TWO CASES THAT DO NOT INVOLVE DO LOOPS: C C 1. ONE ISOLATED POINT, EITHER THE FIRST OR THE LAST C NPART = 1 IPARTL = NC-1 IPARTS = 1 C C A. FIRST C I1 = IFIRST K3 = 3 K2 = 2 200 JA = JDIAG(I1,1) JC = JDIAG(I1,K3) C IF (JA .EQ. JC) THEN IF (NC .GT. 1) THEN WRITE(IPD,300) I1,K3,JA,JC,NC STOP ENDIF NPOINT(1) = I1 GO TO 900 ENDIF C I2 = TAB1(JA,K2) I3 = TAB1(JC,K2) C IF (ABS(IL(I3)-IL(I2))-NCM) 2,4,3 2 WRITE(IPD,301) I2,I3,JA,JC,K2,NC STOP C C C B. LAST C 3 IF (I1 .NE. IFIRST) GO TO 250 I1 = ILAST K3 = 2 K2 = 1 GO TO 200 C 4 NPOINT(1) = I1 I20 = MIN (I2,I3) I21 = IL(I20) I31 = I21+NCM1 IC = I21-2 C DO 203 II = I21,I31 203 NPOINT(II-IC) = IH(II) C IF (NC .LE. 2) THEN C : CALL PHASE(I1,JDIAG,KFL2B) -- PHASE FACTOR C ARISING FROM NON-CYCLIC PERMUTATION OF ARGUMENTS IN TRIAD L: IF (JDIAG(IFIRST,1).NE.JDIAG(ILAST,1)) THEN J7(J7C+1) = JDIAG(I1,1) J7(J7C+2) = JDIAG(I1,2) J7C = J7C + 3 J7(J7C) = JDIAG(I1,3) GO TO 900 ENDIF ENDIF C IF (I1 .NE. ILAST) THEN IT = I2 JT = JDIAG(ILAST,2) K4 = 2 I4 = ILAST ELSE IT = I3 JT = JDIAG(IFIRST,3) K4 = 3 I4 = IFIRST ENDIF C C IF (IT .EQ. I20) CALL PHASE (I1,JDIAG,KFL2B) IF (IT.EQ.I20) THEN C PHASE(I1,JDIAG,KFL2B) -- PHASE FACTOR C ARISING FROM NON-CYCLIC PERMUTATION OF ARGUMENTS IN TRIAD L: J7(J7C+1) = JDIAG(I1,1) J7(J7C+2) = JDIAG(I1,2) J7C = J7C + 3 J7(J7C) = JDIAG(I1,3) ENDIF C IF (JT.EQ.JA .OR. JT.EQ.JC) CALL CHANGE(I4,K4) GO TO 900 C C 2. TWO ISOLATED POINTS, FIRST AND LAST C 250 IF (NC .EQ. 1) RETURN IF (NC .LE. 3) GO TO 100 IPARTL = NC-2 IPARTS = 1 I1 = IFIRST I2 = ILAST JA = JDIAG(I1,1) JB = JDIAG(I1,3) C IF (TAB1(JA,2) .NE. I2) THEN JA = JDIAG(I1,3) JB = JDIAG(I1,1) IF (TAB1(JA,2) .NE. I2) GO TO 100 ENDIF C IF (JA .EQ. JDIAG(I2,1)) THEN JC = JDIAG(I2,2) ELSE JC = JDIAG(ILAST,1) ENDIF C I3 = TAB1(JB,2) I4 = TAB1(JC,1) IDIST = IL(I4)-IL(I3) C IF (ABS(IDIST)-NCM+1) 263,262,100 263 WRITE(IPD,302) I3,I4,JB,JC,IDIST,NC STOP C 262 NPOINT(1) = ILAST NPOINT(2) = IFIRST ICROSS = SIGN(1,IDIST) I20 = MIN(I3,I4) I21 = IL(I20) I31 = I21+NCM IC = I21-3 DO 261 II = I21,I31 261 NPOINT(II-IC) = IH(II) IF (JA .EQ. JDIAG(IFIRST,1)) CALL CHANGE(IFIRST,3) IF (JA .EQ. JDIAG(ILAST,1)) CALL CHANGE(ILAST,2) GO TO 900 C C FIRST GENERAL CASE: ALL POINTS IN ONE GROUP C 100 NPART = 1 IPARTS = 0 IPARTL = NC K3 = 1 C DO 101 IN = 1,NBNODE I = IH(IN) 108 JA = JDIAG(I,K3) I2 = TAB1(JA,2) IF (I .EQ. I2) GO TO 102 IF (IL(I2)-IN-NCM1) 105,104,102 C 105 WRITE(IPD,303) IN,I,I2,IL(I2),JA,NC STOP C 104 I21 = IL(I2) IC = IN-1 DO 103 II = IN,I21 103 NPOINT(II-IC) = IH(II) IF (JA .EQ. JDIAG(IFIRST,3)) CALL CHANGE(IFIRST,3) IF (JA .EQ. JDIAG(ILAST,2)) CALL CHANGE(ILAST,2) GO TO 900 C 102 IF (IN .NE. 1) GO TO 101 IF (K3 .NE. 3) THEN K3 = 3 GO TO 108 ELSE K3 = 1 ENDIF C 101 CONTINUE C C SEARCH DID NOT FIND LOOP NC .LE. 3 C IF (NC .LE. 3) RETURN C C GENERAL CASE OF LOOP PARTITIONNED IN 2 GROUPS, DO LOOP ON IPARTS C NPART = 2 NC2 = NC/2 K3 = 1 K2 = 1 C DO 400 IPS = 2,NC2 JPS = IPS-1 NBN = NBNODE-JPS C DO 401 I1 = 1,NBN I = IH(I1) I2 = IH(I1+JPS) 402 JA = JDIAG(I,K3) JD = JDIAG(I2,K2) C IF (I .EQ. TAB1(JA,1)) THEN II2 = TAB1(JD,2) II1 = TAB1(JA,2) ELSE II1 = TAB1(JA,1) II2 = TAB1(JD,1) ENDIF C IDIST = IL(II1)-IL(II2) C IF (ABS(IDIST)-NCM+JPS) 404,406,420 404 WRITE(IPD,304) JPS,I1,I,I2,JA,JD,II1,II2,IDIST,NC STOP 406 ICROSS = SIGN(1,IDIST) IC = 0 I21 = IL(I2) C DO 410 II = I1,I21 IC = IC+1 410 NPOINT(IC) = IH(II) C I20 = MIN(II1,II2) I30 = MAX(II1,II2) I21 = IL(I20) I31 = IL(I30) C DO 411 II = I21,I31 IC = IC+1 411 NPOINT(IC) = IH(II) C IPARTS = IPS IPARTL = NC-IPS IF (JDIAG(IFIRST,3).EQ.JA .OR. : JDIAG(IFIRST,3).EQ.JD) CALL CHANGE(IFIRST,3) IF (JDIAG(ILAST,2).EQ.JA .OR. : JDIAG(ILAST,2).EQ.JD) CALL CHANGE(ILAST,2) GO TO 900 C 420 IF (I1.EQ.1) THEN IF (K3 .EQ. 3) THEN K3 = 1 GO TO 401 ELSE K3 = 3 GO TO 402 ENDIF ENDIF C IF (I2.NE.ILAST) GO TO 401 IF (K2.NE.2) THEN K2 = 2 GO TO 402 ENDIF C 401 CONTINUE 400 CONTINUE C C SEARCH DID NOT FIND CIRCUIT OF ORDER NC C RETURN C C LOOP FOUND C 900 FIND = .TRUE. CALL PRINTJ('SEARCH',10) C RETURN C C ERROR PRINTOUT C C 300 FORMAT (' ERROR IN SEARCH. I1,K3,JA,JC,NC = ',5I5) 301 FORMAT (' ERROR IN SEARCH. I2,I3,JA,JC,K2,NC = ',6I5) 302 FORMAT (' ERROR IN SEARCH. I3,I4,JB,JC,IDIST,NC = ',6I5) 303 FORMAT (' ERROR IN SEARCH. IN,I,I2,IL(I2),JA,NC = ',6I5) 304 FORMAT (' ERROR IN SEARCH. JPS,I1,I,I2,JA,JD,II1,II2,IDIST,NC = ' : ,10I5) C END C*********************************************************************** SUBROUTINE SQUARE C C REDUCES A CIRCUIT OF ORDER 4 IN THE TWO CASES WHICH ARE LEFT C OVER BY POLYGN, NAMELY TWO DISCONNECTED GROUPS OF TWO POINTS C AND ONE GROUP OF TWO POINTS PLUS THE TWO ENDS OF THE AXIS. IN C THE LATTER THE END OF THE AXIS IS TRANSFERRED TO THE BEGINNING. C IN THIS PROCESS ONE SUMMATION VARIABLE AND TWO 6J SYMBOLS ARE C INTRODUCED. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFL2C=2*KFL2+2, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C INTEGER ARR,TAB1 COMMON/GRAPH/JDIAG(KFL2B,3),ARR(KFL2B,3),TAB1(KFL1,2),IL(KFL2B), + IH(KFL2B),NPOINT(KFL2C),NBNODE,IFIRST,ILAST,IPARTS,IPARTL,NPART, + ICROSS,NFREE,ITFREE(KFL2A),NFIN,NC LOGICAL SUMVAR COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP C CHARACTER*6 NAMSUB COMMON /NAM/NAMSUB C NAMSUB = 'SQUARE' MP = MP+1 SUMVAR(MP) = .TRUE. IT1 = NPOINT(1) IT2 = NPOINT(2) C IF (ICROSS .EQ. 1) THEN IT3 = NPOINT(3) IT4 = NPOINT(4) K23 = 3 K32 = 2 ELSE IT3 = NPOINT(4) IT4 = NPOINT(3) K23 = 2 K32 = 3 ENDIF C L4 = JDIAG(IT2,1) C IF (ARR(IT2,1) .LE. 0) THEN C CALL PHASE2(L4) ! ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = L4 ARR(IT2,1) = 1 ARR(IT3,1) = -1 ENDIF C L2 = JDIAG(IT1,1) IF (ARR(IT1,1) .GT. 0) THEN C CALL PHASE2(L2) ! ADDS A FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = L2 ENDIF C JWC = JWC+1 KW(1,JWC) = L4 KW(2,JWC) = L2 KW(3,JWC) = JDIAG(IT2,2) JJ1 = JDIAG(IT1,3) KW(4,JWC) = JJ1 KW(5,JWC) = MP KW(6,JWC) = JDIAG(IT1,2) IF (ARR(IT1,2).LT.0) THEN C CALL PHASE2 (JDIAG(IT1,2)) ! ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = JDIAG(IT1,2) ENDIF JWC = JWC+1 KW(1,JWC) = L4 KW(2,JWC) = L2 JJ3 = JDIAG(IT3,K23) JJ2 = JDIAG(IT4,K32) KW(3,JWC) = JJ3 KW(4,JWC) = JJ2 KW(5,JWC) = MP KW(6,JWC) = JDIAG(IT3,K32) IF (ARR(IT3,K32).LT.0) THEN C CALL PHASE2 (JDIAG(IT3,K32)) ! ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = JDIAG(IT3,K32) ENDIF J6(J6C+1) = MP J6C = J6C+2 J6(J6C) = MP C IF (NPART .EQ. 1) THEN ITMIN = IT2 ITMAX = IT3 ELSE ITMIN = MIN (IT2,IT3) ITMAX = MAX (IT2,IT3) ENDIF ITMN = MIN (IT1,IT4) ITMX = MAX (IT1,IT4) C TAB1(MP,1) = ITMIN TAB1(MP,2) = ITMAX JDIAG(IT2,1) = MP JDIAG(IT3,1) = MP JDIAG(IT2,3) = JJ1 ARR(IT2,3) = ARR(IT1,3) JDIAG(IT3,K32) = JJ2 ARR(IT3,K32) = ARR(IT4,K32) C J8C = J8C + 1 IF (ICROSS .EQ. 1) THEN J7(J7C+1) = L2 J7(J7C+2) = L4 C CALL PHASE2(L4): J8C=J8C+1 & J8(J8C) = L4 J7C = J7C+3 J7(J7C) = MP ELSE C CALL PHASE2(JJ2): J8C=J8C+1 & J8(J8C) = JJ2 ENDIF C ITLL = IL(ITMN) ITHL = IL(ITMX) C DO 5 I = ITLL+1,ITHL-1 IT = IH(I) ILP = I-1 IL(IT) = ILP IH(ILP) = IT 5 CONTINUE IF (ITHL .NE. NBNODE) THEN DO 6 I = ITHL+1,NBNODE IT = IH(I) ILP = I-2 IL(IT) = ILP IH(ILP) = IT 6 CONTINUE ENDIF C IF (NPART .NE. 2) THEN TAB1(JJ1,1) = IH(1) TAB1(JJ1,2) = IH(NBNODE-2) ENDIF C RETURN END C*********************************************************************** SUBROUTINE ZERO(J,JZ,FAIL) C C SUPPRESSES ONE LINE AND TWO NODES OF THE UNSTRUCTURED GRAPH, C INTRODUCES ZEROS IN THE TRIADS J23. AS A CONSEQUENCE THE OTHER C TWO ARGUMENTS OF THE TRIAD ARE PUT EQUAL. IF THERE WAS ALREADY C A ZERO IN THE TRIAD WHICH IS CHANGED IT IS A SPECIAL CASE. C PARAMETER (LL74= 21+4, LL76= 21*3+19) PARAMETER(KFL1=LL76,KFL2=LL74,KFL2A=2*KFL2,KFL2B=4*KFL2, + KFLZ=50, + KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20) C LOGICAL FAIL COMMON/ZER/NZERO,JZERO(KFLZ) COMMON/BUILD/IAL(KFL2B),IF1,IF2,NODE C LOGICAL TABS,FREE,SUMVAR,CUT,NOCUT INTEGER ARROW COMMON/TREE/J23(KFL2A,3),ARROW(KFL2A,3),LINE(KFL1,2),LCOL(KFL1,2), + TABS(KFL2A),NBTR COMMON/COUPLE/M,N,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1) COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8), + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP COMMON/CUTDIG/CUT C C NOCUT=.FALSE. NZERO=0 IF(J.LT.1) GO TO 12 C CALL OTHERJ(0,JZ,LIN,LC,K1) C GIVES THE OTHER TRIAD (WHERE A GIVEN J OCCURS) AND ITS POSITION: LIN = LINE(JZ,1) IF (LIN.EQ.0 .OR. TABS(LIN)) THEN K1 = 1 LIN = LINE(JZ,2) LC = LCOL(JZ,2) ELSE K1 = 2 LC = LCOL(JZ,1) ENDIF I=NZERO GO TO 8 12 DO 11 I=1,M IF(J1(I).NE.1.OR.FREE(I).OR.IAL(I).LE.1) GO TO 11 NZERO=NZERO+1 IF(NZERO.GT.KFLZ) THEN PRINT 100, NZERO,KFLZ 100 FORMAT(' DIMENSION ERROR IN ZERO: NZERO=',I3,' > KFLZ=',I3) STOP ENDIF JZERO(NZERO)=I 11 CONTINUE NOCUT=.TRUE. M=M+1 J1(M)=1 SUMVAR(M)=.FALSE. FREE(M)=.FALSE. IF(NZERO.EQ.0) GO TO 7 CALL PRINTJ('ZERO_1',1) C I=0 1 I=I+1 JZ=JZERO(I) J=0 13 J=J+1 LIN=LINE(JZ,J) IF(TABS(LIN)) GO TO 2 LC=LCOL(JZ,J) C 8 CALL NEIBOR(LC,L1,L2) C GIVES THE POSITIONS OF THE OTHER TWO ARGUMENTS IN THE TRIAD: 8 IF (LC.LT.2) THEN L1 = 2 L2 = 3 ELSE IF (LC.EQ.2) THEN L1 = 3 L2 = 1 ELSE L1 = 1 L2 = 2 ENDIF JJ1=J23(LIN,L1) JJ2=J23(LIN,L2) IF(JJ1.NE.JJ2) GO TO 9 J6C=J6C+1 J6(J6C)=JJ1 GO TO 10 9 CALL DELTA(JJ1,JJ2,FAIL) IF(FAIL) GO TO 7 IF(J1(JJ1).NE.1.AND.J1(JJ2).NE.1) GO TO 15 IF(J1(JJ1)-J1(JJ2)) 15,16,19 16 IF(NZERO.EQ.0) GO TO 15 DO 17 K=I,NZERO IF(JJ1.EQ.JZERO(K)) GO TO 15 IF(JJ2.EQ.JZERO(K)) GO TO 19 17 CONTINUE GO TO 15 19 L=JJ2 JJ2=JJ1 JJ1=L C 15 CALL OTHERJ(LIN,JJ1,LO1,LCO1,K1) C GIVES THE OTHER TRIAD... AND ITS POSITION: 15 LO1 = LINE(JJ1,1) IF (LO1.EQ.LIN .OR. TABS(LO1)) THEN K1 = 1 LO1 = LINE(JJ1,2) LCO1 = LCOL(JJ1,2) ELSE K1 = 2 LCO1 = LCOL(JJ1,1) ENDIF C CALL OTHERJ(LIN,JJ2,LO2,LCO2,K2) C GIVES THE OTHER TRIAD... AND ITS POSITION: LO2 = LINE(JJ2,1) IF (LO2.EQ.LIN .OR. TABS(LO2)) THEN LO2 = LINE(JJ2,2) LCO2 = LCOL(JJ2,2) ELSE LCO2 = LCOL(JJ2,1) ENDIF J9C=J9C+1 J9(J9C)=JJ1 J23(LO2,LCO2)=JJ1 LINE(JJ1,K1)=LO2 LCOL(JJ1,K1)=LCO2 10 IF(ARROW(LIN,L1)-ARROW(LIN,L2)) 3,4,5 C 3 CALL PHASE2(JJ1) -- ADDS A PHASE FACTOR (-1)**2J: 3 J8C = J8C + 1 J8(J8C) = JJ1 GO TO 5 4 ARROW(LO1,LCO1)=1 ARROW(LO2,LCO2)=-1 5 TABS(LIN)=.TRUE. NBTR=NBTR-1 IF(NBTR.EQ.0) GO TO 7 IF(LO1.NE.LO2) GO TO 2 L=6-LCO1-LCO2 JT=J23(LO1,L) IF(J1(JT).EQ.1.AND..NOT.FREE(JT)) GO TO 2 CALL DELTA(JT,M,FAIL) IF(FAIL) GO TO 7 C CALL NEIBOR(L,L1,L2) C GIVES THE POSITIONS OF THE OTHER TWO ARGUMENTS IN THE TRIAD: IF (L.LT.2) THEN L1 = 2 L2 = 3 ELSE IF (L.EQ.2) THEN L1 = 3 L2 = 1 ELSE L1 = 1 L2 = 2 ENDIF JTF=J23(LO1,L1) IF(ARROW(LO1,L1).LT.ARROW(LO1,L2)) THEN C CALL PHASE2(JTF) -- ADDS A PHASE FACTOR (-1)**2J: J8C = J8C + 1 J8(J8C) = JTF ENDIF J6C=J6C+1 J6(J6C)=JTF NBTR=NBTR-1 TABS(LO1)=.TRUE. C CALL OTHERJ(LO1,JT,LIN,LC,K) C GIVES THE OTHER TRIAD... AND ITS POSITION: LIN = LINE(JT,1) IF (LIN.EQ.LO1 .OR. TABS(LIN)) THEN LIN = LINE(JT,2) LC = LCOL(JT,2) ELSE LC = LCOL(JT,1) ENDIF GO TO 8 2 IF(J.EQ.1) GO TO 13 IF(NBTR.EQ.0) GO TO 7 IF(I.LT.NZERO) GO TO 1 C 7 CALL PRINTJ('ZERO_2',4) if(NBTR.lt.0.and..not.FAIL) then ! FAIL-tmpmod '05Aug13-14 print "(' ZERO_2 stopgap: FAIL=.t. since NBTR =',I2)",NBTR FAIL=.true. endif IF(NOCUT) CUT=.FALSE. RETURN END