C N. R. BADNELL UoS v3.7 29/06/21 C PROGRAM MNSTG3 IMPLICIT REAL*8 (A-H,O-Z) C C C C*********************************************************************** C C Belfast Originated Atomic R-matrix Codes C C*********************************************************************** C C THE THIRD PART OF C C A GENERAL PROGRAM TO CALCULATE ATOMIC CONTINUUM C C PROCESSES USING THE R-MATRIX METHOD C C S T G 3 C C ORIGINALLY DISTRIBUTED BY C C QUEEN'S UNIVERSITY BELFAST C C*********************************************************************** C C DIAGONALIZES THE HAMILTONIANS AND PROCESSES DIPOLE MATRIX C ELEMENTS USING AS INPUT THE H (AND PERHAPS THE D) FILES C CREATED IN STG2 OR RECUPD. C C CAN ALSO HANDLE H AND D FILES PRODUCED BY DARC, & DTO3 INTERFACE. C C WILL USE A PARTIONED R-MATRIX IF IPRCENT SPECIFIED (40-100). C C REQUIRES TO BE LINKED WITH STGLIB LIBRARY ROUTINES C OR OTHER SUITABLE PACKAGE FOR DIAGONALIZING. C C C*********************************************************************** C C ROUTINES USED IN STG3 C C*********************************************************************** C C MNSTG3 C STG3 DRIVER C DA2 C DIAG C DMAT C DMAT1 C DMAT2 C DMAT2P C HAMNEW C HAMNEW2 C ORDER C RECOV2 C RSCT C STG3RD C TAPERD C VMUL C WRITOP C C*********************************************************************** C C INPUT/OUTPUT CHANNELS USED IN STG3 C C*********************************************************************** C C IREAD USER INPUT FILE (FILE dstg3) C IWRITE OUTPUT TO LINE PRINTER (FILE rout3r) C C IDISC1 OPTIONAL SCRATCH FILE FOR THE H-MATRIX C IDISC2 SCRATCH DIRECT ACCESS FILE OF H EIGENVECTORS C IDISC3 SCRATCH DIRECT ACCESS FILE OF PARTIALLY C PROCESSED DIPOLE MATRIX ELEMENTS C C ITAPE1 INPUT FILE OF STG2 DIPOLE MATRICES C ITAPE2 INPUT FILE OF STG2 H-MATRICES C ITAPE3 OUTPUT FILE OF DIAGONALIZED H-MATRICES C ITAPE4 OUTPUT FILE OF PROCESSED DIPOLE MATRICES. C C IREAD (5) .. input data .. dstg3 C IWRITE (6) .. printed output .. rout3r C C IPUNCH .. NOT USED C C IDISC1 (11) .. scratch C IDISC2 (12) .. scratch (DA2) C IDISC3 (13) .. scratch C IDISC4 .. NOT USED C C ITAPE1 (1) .. RECUPD/STG2 dump (dipole) .. RECUPD/STG2D.DAT C if IPOLPH=2 C ITAPE2 (2) .. RECUPD/STG2 dump (hamiltonians) .. RECUPH/STG2H.DAT C always used C ITAPE3 (3) .. STG3 dump (hamiltonians) .. H.DAT C always used C ITAPE4 (4) .. STG3 dump (dipole matrices) .. D00 etc C if IPOLPH=2 C C*********************************************************************** C C DIMENSIONING PARAMETERS USED IN STG3 C C*********************************************************************** C C INCLUDE PARAMETERS: C C CHF (75) HIGHEST NUMBER OF CHANNELS (NCHAN) C LMX (8) MULTIPOLES IN POTENTIAL (LAMAX) C LR2 (20) HIGHEST L+1 FOR CONTINUUM ORBITALS (LRANG2) C MEG (1) MEGA-WORDS OF MEMORY TO REDUCE I/O (CAN BE 0) CFIXED KIL (100) KILO-WORDS OF MEMORY TO REDUCE I/O (CAN BE 0) C NC2 (600) N+1 ELECTRON CONFIGURATIONS FOR GIVEN S L PI (NCFGP) C NR2 (40) NUMBER OF CONTINUUM ORBITALS FOR GIVEN L (NRANG2) C SLP (80) NUMBER OF DIFFERENT N+1 ELECTRON SYMMETRIES (INAST) C TAR (63) TARGET STATES TOTAL (NAST) C INCLUDE 'PARAM' C parameter (nsizex=3) !min 1 PARAMETER (MXMEM=MZMEG*1000000+MZKIL*1000+1) PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) PARAMETER (MXTRI=(MXMAT*(MXMAT+1))/2) PARAMETER (MXN21=MZNR2+1) C C MXDM=MAX(MZNC2,MZNR2): C PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) C C MXMATX=MAX(MXDMAT,MXRSCT)+1 = SIZE OF /MATRIX/ IN DMAT1, RSCT: C PARAMETER (MXDMAT=4*MXMAT*MZCHF+2*MXMAT*MXDM+2*MZNR2*MXDM, A MXRSCT=MXTRI+MZCHF*MXMAT+MZNR2*MXMAT) PARAMETER (MXM1=MXDMAT/MXRSCT,MXM2=MXRSCT/MXDMAT,MXM3=MXM1+MXM2, A MXMATX=MXDMAT*MXM1/MXM3+MXRSCT*MXM2/MXM3+1) C C*********************************************************************** C C COMMON BLOCKS USED IN STG3 C C*********************************************************************** C common/adip/ad(mznc2,mznc2,nsizex),nlad(mzslp),nrad(mzslp), c ixsym(mzslp),nxsym COMMON /ADDE/EST(MZTAR),EPOLE,IPOLE,NEST,AST(MZTAR),MERGE(MZTAR) COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /BASIC2/CF(MZCHF,MZCHF,MZLMX),ET(MZCHF),KJ(MZCHF), A L2P(MZCHF),LSTARG(MZCHF),NCONAT(MZTAR),LRGL,NSPN,NPTY COMMON /BASIN/EIGENS(MZNR2,MZLR2),ENDS(MXN21,MZLR2),DELTA,ETA COMMON /BUTT/COEFF(3,MZLR2),EK2MAX,EK2MIN,MAXNCB(MZLR2),NELCOR COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST COMMON /COPY/ICOPY1,ICOPY2,ITOTAL,ICOUNT COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /HAMNUW/A(MZNC2,MZNC2),XO(MZNC2,MZNC2),XL(MZNC2) X ,WORK(MZNC2) COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /MATRIX/DUME(MXMATX) COMMON /MEMORY/ARRAY(MXMEM),MEM1,MREC1 COMMON /MORDER/LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECORD/ILSPI(4,MZSLP),IORDER(MZCHF,MZSLP) COMMON /RECOV/IPLACE COMMON /REDTOL/TOLER,iparinit COMMON /SPACES/XX(MZNR2,MXMAT) COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST C C SUN REAL*4 TARRY(2),TIME C C C*********************************************************************** C C C STG3 MAIN PROGRAM C C*********************************************************************** C C MEM1 AND MREC1 ARE THE MEMORY AND DA FILE POINTERS C CMICD GETCPUS MEM1 = 0 MREC1 = -1 CALL STG3 C C SUN DUM=DTIME(TARRY) TIME=TARRY(1) C CRAY CRAY CALL SECOND(TIME) C TIME=TIME/60.0 WRITE(IWRITE,999) TIME 999 FORMAT(//1X,'CPU TIME=',F9.3,' MIN') C STOP END SUBROUTINE DA2(KEY,IREC,JDISC,LENGTH,ARRAY) IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C TO STORE A LARGE ARRAY IN A DA FILE OF RECORD LENGTH 8*LREC BYTES. C C KEY = 1 FOR READ, C = 2 FOR WRITE, C = 0 FOR FINDING NUMBER OF DA RECORDS A GIVEN ARRAY TAKES. C C IREC= (ON CALL) POINTER TO FIRST DA RECORD FOR ARRAY, C = 0 FOR OPENING DA FILE (BY NAME), C =-1 FOR OPENING DA FILE (SCRATCH), C = (ON RETURN) POINTER TO NEXT AVAILABLE DA RECORD. C C JDISC = DA FILE UNIT NUMBER. C C ARRAY(LENGTH) = ARRAY TO READ OR WRITE. C C----------------------------------------------------------------------- PARAMETER (LREC=512) C DIMENSION ARRAY(LENGTH) C----------------------------------------------------------------------- IF (IREC.GT.0) GOTO 20 C IRECL = 8*LREC C IF (IREC.LT.0) THEN OPEN (JDISC,STATUS='SCRATCH',ACCESS='DIRECT',FORM='UNFORMATTED', X RECL=IRECL) GOTO 10 ENDIF C IF (KEY.EQ.2) THEN OPEN (JDISC,STATUS='UNKNOWN',FILE='NAME.DAT',ACCESS='DIRECT', A RECL=IRECL) ELSE OPEN (JDISC,STATUS='OLD',FILE='NAME.DAT',ACCESS='DIRECT', A RECL=IRECL) ENDIF C 10 CONTINUE IREC = 1 C 20 CONTINUE IF (LENGTH.EQ.0) RETURN I2 = 0 30 CONTINUE I1 = I2 + 1 I2 = MIN(I2+LREC,LENGTH) IF (KEY.EQ.0) GOTO 40 C IF (KEY.EQ.2) THEN WRITE (JDISC,REC=IREC) (ARRAY(I),I=I1,I2) ELSE READ (JDISC,REC=IREC) (ARRAY(I),I=I1,I2) ENDIF C 40 CONTINUE IREC = IREC + 1 IF (I2.LT.LENGTH) GOTO 30 C END SUBROUTINE DIAG(N,IUP,Z,D,E,MXMAT) C IMPLICIT REAL*8 (A-H,O-Z) C C BADNELL & BURGESS D.A.M.T.P. CAMBRIDGE C C DIAGONALIZATION OF REAL SYMMETRIC N-BY-N MATRIX Z. C C METHOD: HOUSEHOLDER REDUCTION TO TRI-DIAGONAL FORM AND SHIFTED C QL ALGORITHM TO DETERMINE THE E-VALUES AND E-VECTORS. C C BASED ON MARTIN, REINSCH & WILKINSON: NUM. MATH. 11, 181-95 (1968). C C INPUT REQUIRED. N, IUP AND Z. ONLY LOWER TRIANGLE OF Z NEED BE SUPPLIED. C MATRIX Z OVERWRITTEN BY EIGENVECTORS OF Z. C IUP=1/-1 ASC/DESCENDING SORT, 0 NO SORT. C MXMAT, IS THE ROW DIMENSION OF Z IN THE CALLING ROUTINE. C C OUTPUT. Z AND D, WHERE Z CONSISTS OF COLUMN EIGENVECTORS C AND D CONSISTS OF CORRESPONDING EIGENVALUES. C C NOTE: E IS A WORKING ARRAY. C PARAMETER (TOL = 1.0D-75) PARAMETER (EPS = 1.0D-15) PARAMETER (ZERO = 0.0D0) PARAMETER (ONE = 1.0D0) PARAMETER (JMAX = 30) C DIMENSION D(N),E(N),Z(MXMAT,N) C C DO 1 I = 1,N D(I) = Z(N,I) 1 CONTINUE IF (N.LE.1) GO TO 20 C C HOUSEHOLDER REDUCTION TO TRI-DIAGONAL FORM C DO 19 I = N,2,-1 L = I - 1 F = D(I-1) G = ZERO DO 5 K = 1,I-2 G = G + D(K)*D(K) 5 CONTINUE H = G + F*F IF (G.GT.TOL) GO TO 8 E(I) = F H = ZERO DO 7 J = 1,L D(J) = Z(L,J) Z(I,J) = ZERO Z(J,I) = ZERO 7 CONTINUE GO TO 18 8 G = SQRT(H) IF (F.GE.ZERO) G = -G E(I) = G H = H - F*G D(L) = F - G DO 14 J = 1,L E(J) = ZERO 14 CONTINUE DO 15 J = 1,L Z(J,I) = D(J) G = E(J) + Z(J,J)*D(J) DO 13 K = J+1,L G = G + Z(K,J)*D(K) E(K) = E(K) + Z(K,J)*D(J) 13 CONTINUE E(J) = G 15 CONTINUE F = ZERO DO 12 J = 1,L E(J) = E(J)/H F = F + E(J)*D(J) 12 CONTINUE HH = F/(H+H) DO 11 J = 1,L E(J) = E(J) - HH*D(J) 11 CONTINUE DO 17 J = 1,L F = D(J) G = E(J) DO 16 K = J,L Z(K,J) = Z(K,J) - F*E(K) - G*D(K) 16 CONTINUE D(J) = Z(L,J) Z(I,J) = ZERO 17 CONTINUE 18 D(I) = H 19 CONTINUE C C C ACCUMULATE TRANSFORMATION MATRICES C DO 28 I = 2,N L = I - 1 Z(N,L) = Z(L,L) Z(L,L) = ONE H = D(I) IF (H.EQ.ZERO) GO TO 25 DO 21 K = 1,L D(K) = Z(K,I)/H 21 CONTINUE DO 24 J = 1,L G = ZERO DO 22 K = 1,L G = G + Z(K,I)*Z(K,J) 22 CONTINUE DO 23 K = 1,L Z(K,J) = Z(K,J) - G*D(K) 23 CONTINUE 24 CONTINUE 25 DO 27 J = 1,L Z(J,I) = ZERO 27 CONTINUE 28 CONTINUE DO 29 I = 1,N D(I) = Z(N,I) Z(N,I) = ZERO 29 CONTINUE 20 E(1) = ZERO Z(N,N) = ONE C C C SHIFTED QL ALGORITHM TO DETERMINE E-VALUES & E-VECTORS C DO 32 I = 2,N E(I-1) = E(I) 32 CONTINUE E(N) = ZERO B = ZERO F = ZERO DO 54 L = 1,N J = 0 H = EPS*(ABS(D(L))+ABS(E(L))) IF (B.LT.H) B = H DO 36 M = L,N IF (ABS(E(M)).LE.B) GO TO 37 36 CONTINUE 37 IF (M.EQ.L) GO TO 53 38 IF (J.EQ.JMAX) GO TO 62 J = J+ 1 P = E(L) + E(L) G = D(L) H = D(L+1) - G IF (ABS(H).GE.ABS(E(L))) GO TO 43 P = H/P R = SQRT(P*P+ONE) H = P + R IF (P.LT.ZERO) H = P - R D(L) = E(L)/H GO TO 44 43 P = P/H R = SQRT(P*P+ONE) D(L) = E(L)*P/(R+ONE) 44 H = G - D(L) DO 46 I = L+1,N D(I) = D(I) - H 46 CONTINUE F = F + H P = D(M) C = ONE S = ZERO DO 52 I = M-1,L,-1 G = C*E(I) H = C*P IF (ABS(P).LT.ABS(E(I))) GO TO 49 C = E(I)/P R = SQRT(C*C+ONE) E(I+1) = S*P*R S = C/R C = ONE/R GO TO 50 49 C = P/E(I) R = SQRT(C*C+ONE) E(I+1) = S*E(I)*R S = ONE/R C = C/R 50 P = C*D(I) - S*G D(I+1) = H + S*(C*G+S*D(I)) DO 51 K = 1,N H = Z(K,I+1) Z(K,I+1) = S*Z(K,I) + C*H Z(K,I) = C*Z(K,I) - S*H 51 CONTINUE 52 CONTINUE E(L) = S*P D(L) = C*P IF (ABS(E(L)).GT.B) GO TO 38 53 D(L) = D(L) + F 54 CONTINUE C IF(IUP.EQ.0)RETURN C C BEGIN SORTING INTO ASCENDING E-VALUES C DO 61 I = 1,N K = I P = D(I) DO 57 J = I+1,N IF (IUP.GT.0.AND.D(J).GT.P) GO TO 57 IF (IUP.LT.0.AND.D(J).LT.P) GO TO 57 K = J P = D(J) 57 CONTINUE IF (K.EQ.I) GO TO 61 D(K) = D(I) D(I) = P DO 60 J = 1,N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 60 CONTINUE 61 CONTINUE C RETURN C C 62 WRITE(6,100) 100 FORMAT(' FAILED IN DIAG, TOO MANY ITERATIONS') RETURN C END SUBROUTINE DMAT IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C WRITES 'D FILE' FOR INTERMEDIATE COUPLING ONTO TAPE4, C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C parameter (nsizex=3) !min 1 PARAMETER (MXMEM=MZMEG*1000000+MZKIL*1000+1) PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) PARAMETER (MXTRI=(MXMAT*(MXMAT+1))/2) PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) PARAMETER (MXDMAT=4*MXMAT*MZCHF+2*MXMAT*MXDM+2*MZNR2*MXDM, A MXRSCT=MXTRI+MZCHF*MXMAT+MZNR2*MXMAT) PARAMETER (MXM1=MXDMAT/MXRSCT,MXM2=MXRSCT/MXDMAT,MXM3=MXM1+MXM2, A MXMATX=MXDMAT*MXM1/MXM3+MXRSCT*MXM2/MXM3+1) PARAMETER (MXDUM=MXMATX-4*MXMAT*MZCHF) C CHARACTER FILNAM*7,NUM(0:9)*1 C DIMENSION IMEM(2,99) DIMENSION CGC(MZLR2),AC(MZCHF,MZCHF), A BLC(MZCHF,MZCHF),BVC(MZCHF,MZCHF) DIMENSION DUM1(MXMAT),DUM2(MXMAT) C common/adip/ad(mznc2,mznc2,nsizex),nlad(mzslp),nrad(mzslp), c ixsym(mzslp),nxsym COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /MATRIX/ABUTL(MXMAT,MZCHF),ABUTV(MXMAT,MZCHF), A BBUTL(MXMAT,MZCHF),BBUTV(MXMAT,MZCHF),DUM(MXDUM) COMMON /MEMORY/ARRAY(MXMEM),MEM1,MREC1 COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECORD/ILSPI(4,MZSLP),IORDER(MZCHF,MZSLP) COMMON /REDTOL/TOLER,iparinit COMMON /NRBDIP/MAXLD C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C----------------------------------------------------------------------- C C INITIALIZE C NOTERM=MIN(NOT1,NRANG2) -- NOT1 RESET IN RSCT, HENCE C NOTERM = NOT1 WRITE (IWRITE,3000) KOUNT = 0 LTEST = 1 LRECW3 = NOTERM*NOTERM*2 + 1 C C-----DETERMINE THE NUMBER OF DIPOLE MATRIX ELEMENT SETS C DO 20 JN = 2,INAST C C-----FINAL STATE DATA IG = ILSPI(1,JN) JPTY = ISIGN(1,IG) ITEMP = ABS(IG) JLRGL = MOD(ITEMP,100) - 1 IF(JLRGL.GT.MAXLD)GO TO 20 !NRB JSPN = ITEMP/100 IF (JSPN.EQ.0) LTEST = 2 DO 10 IN = 1,JN - 1 C C-----INITIAL STATE DATA IG = ILSPI(1,IN) IPTY = ISIGN(1,IG) ITEMP = ABS(IG) ILRGL = MOD(ITEMP,100) - 1 IF(ILRGL.GT.MAXLD)GO TO 10 !NRB ISPN = ITEMP/100 C C-----CHECK ELECTRIC DIPOLE SELECTION RULES ARE SATISFIED C IF (ABS(JLRGL-ILRGL).GT.LTEST .OR. A (JLRGL+ILRGL).LT.LTEST) GOTO 10 IF (ISPN.NE.JSPN .OR. JPTY.EQ.IPTY) GOTO 10 IF (KOUNT.GE.99) THEN WRITE (IWRITE,3010) GOTO 30 C ENDIF C KOUNT = KOUNT + 1 IMEM(1,KOUNT) = IN IMEM(2,KOUNT) = JN NCFGP = ILSPI(2,JN) - NOTERM*ILSPI(3,JN) LRECW3 = MAX (LRECW3,NCFGP*NOTERM*2+1) 10 CONTINUE 20 CONTINUE C C-----WRITE THIS BASIC INFORMATION TO OUTPUT FILE C 30 CONTINUE OPEN (ITAPE4,FILE='D00',STATUS='UNKNOWN',FORM='UNFORMATTED') WRITE (ITAPE4) KOUNT IF (KOUNT.EQ.0) WRITE (IWRITE,3020) DO 40 IH = 1,KOUNT ITEMP = ILSPI(1,IMEM(1,IH)) IPTY = 0 IF (ITEMP.LT.0) IPTY = 1 ITEMP = ABS(ITEMP) ILRGL = MOD(ITEMP,100) - 1 ISPN = ITEMP/100 ITEMP = ILSPI(1,IMEM(2,IH)) JPTY = 0 IF (ITEMP.LT.0) JPTY = 1 ITEMP = ABS(ITEMP) JLRGL = MOD(ITEMP,100) - 1 JSPN = ITEMP/100 WRITE (ITAPE4) ISPN,ILRGL,IPTY,JSPN,JLRGL,JPTY 40 CONTINUE C OUT ENDFILE(ITAPE4) CLOSE (ITAPE4) C --- *=*=*=*=*=*=*=* ITAPE = ITAPE1 REWIND ITAPE C C THE DIPOLE MATRIX ELEMENTS ARE READ FROM THE STG2 OUTPUT FILE C DO 120 IH = 1,KOUNT C C-----DEFINE INITIAL STATE DATA IN = IMEM(1,IH) ILRGL = MOD(ABS(ILSPI(1,IN)),100) - 1 MMNP2 = ILSPI(2,IN) MCHAN = ILSPI(3,IN) C C-----DEFINE FINAL STATE DATA JN = IMEM(2,IH) JPTY = 0 ITEMP = ILSPI(1,JN) IF (ITEMP.LT.0) JPTY = 1 ITEMP = ABS(ITEMP) JLRGL = MOD(ITEMP,100) - 1 JSPN = ITEMP/100 MNP2 = ILSPI(2,JN) NCHAN = ILSPI(3,JN) C C WRITE TO OUTPUT FILE DATA DEFINING THE S L PI OR J PI C COMBINATION SATISFYING THE DIPOLE SELECTION RULES C WRITE (6,3040) IN,JN FILNAM = 'D'//NUM(IH/10)//NUM(MOD(IH,10)) C //'.DAT' OPEN (ITAPE4,FILE=FILNAM,STATUS='UNKNOWN',FORM='UNFORMATTED') WRITE (ITAPE4) NOTERM,MNP2,NCHAN,JLRGL,JPTY,JSPN,MMNP2,MCHAN, A ILRGL C C ---- CALCULATE AND WRITE OUT DIPOLE MATRIX ELEMENTS C C CALL DMAT2, WITH SCRATCH SPACE IN /MEMORY/, IF SUFFICIENT. C IF EIGENVECTORS ALREADY IN MEMORY, SET IV AND JV TO POINT TO THEM C C********************************************************************** C NRB: **** DMAT2 REVERSED ROLES OF ABUT AND BBUT WRT DMAT & DMAT1 **** C **** COMPENSATE BY REVERSING THEM IN THE CALL TO DMAT2 **** C********************************************************************** C mnp22=mnp2 jnx=ixsym(jn) if(jnx.gt.0)mnp22=mnp22+nrad(jnx)-nlad(jnx) c IV = ABS(ILSPI(4,IN)) + 1 JV = ABS(ILSPI(4,JN)) + 1 KV = MEM1 + 2*MMNP2*mnp22 IF (ILSPI(4,IN).GT.0) THEN IV = KV + 1 KV = KV + MMNP2**2 ENDIF C IF (ILSPI(4,JN).GT.0) THEN JV = KV + 1 KV = KV + MNP2**2 ENDIF C IF (KV.LE.MXMEM) THEN ML = MEM1 + 1 MV = ML + MMNP2*mnp22 if(mnp22.eq.mnp2)then ! ABUT <-> BBUT !NRB CALL DMAT2(IN,JN,NRANG2,NOTERM,MMNP2,MNP2,bBUTL,bBUTV,aBUTL, A aBUTV,ARRAY(ML),ARRAY(MV),ARRAY(IV),ARRAY(JV),DUM1,DUM2) else !dipole pseudo reduction CALL DMAT2P(IN,JN,NRANG2,NOTERM,MMNP2,MNP2,bBUTL,bBUTV,aBUTL A ,aBUTV,ARRAY(ML),ARRAY(MV),ARRAY(IV),ARRAY(JV),DUM1,DUM2 B ,mnp22) endif C ELSE C C /MEMORY/ INSUFFICIENT, CALL DMAT1, WITH SCRATCH SPACE IN /MATRIX C C NRB: ABUT, BBUT IN DMAT1 ARE UNCHANGED AND CONSISTENT WITH THE WRITES C M = 1 + (KV-1)/1000000 IF(TOLER.LT.1.D-19)THEN WRITE (IWRITE,*) ' *** TO IMPROVE EFFICIENCY, INCREASE ', A 'MZMEG TO ',M CALL DMAT1(IN,JN,NRANG2,NOTERM,MMNP2,MNP2,LRECW3) ELSE WRITE (IWRITE,*)' *** INCREASE MZMEG FOR DIPOLE PSEUDO ', A 'RESONANCE ELIMINATION' CALL RECOV2('DMAT ','MZMEG ',MZMEG,M) ENDIF ENDIF C C-----WRITE OUT BLOCKS INVOLVING BUTTLE CORRECTIONS, REORDERING CHANNELS C C NRB: WRITES UNCHANGED HERE SINCE NEED TO HANDLE DMAT1 AND DMAT2 C BUT NOTE THAT PSTGD_DIP HANDLES THE PROBLEM DIFFERENTLY. C WRITE (ITAPE4) ((BBUTL(I,IORDER(J,IN)),J=1,MCHAN),I=1,MNP2), A ((BBUTV(I,IORDER(J,IN)),J=1,MCHAN),I=1,MNP2) WRITE (ITAPE4) ((ABUTL(I,IORDER(J,JN)),J=1,NCHAN),I=1,MMNP2), A ((ABUTV(I,IORDER(J,JN)),J=1,NCHAN),I=1,MMNP2) C C READ AND WRITE OUT BUTTLE-BUTTLE ELEMENTS: BBUT(NCHAN,MCHAN) C READ (ITAPE) ((BBUTL(I,J),J=1,NCHAN),I=1,MCHAN) READ (ITAPE) ((BBUTV(I,J),J=1,NCHAN),I=1,MCHAN) WRITE (ITAPE4) ((BBUTL(IORDER(I,IN),IORDER(J,JN)),J=1,NCHAN), A I=1,MCHAN), ((BBUTV(IORDER(I,IN),IORDER(J,JN)),J=1,NCHAN),I=1, B MCHAN) C C READ IN THE CLEBSCH GORDAN COEFFICIENTS FOR THE POLARIZABILITY C CALCULATION, WRITE THESE TO THE OUTPUT FILE: C READ (ITAPE) MAXM1, (CGC(M),M=1,MAXM1) WRITE (ITAPE4) MAXM1, (CGC(M),M=1,MAXM1) C C-----READ IN AND WRITE OUT THE ANGULAR INTEGRALS NEEDED FOR OUTER C-----REGION INTEGRATION, REORDERING CHANNELS. C READ (ITAPE) ((AC(I,J),J=1,MCHAN),I=1,NCHAN) READ (ITAPE) ((BLC(I,J),J=1,MCHAN),I=1,NCHAN) READ (ITAPE) ((BVC(I,J),J=1,MCHAN),I=1,NCHAN) WRITE (ITAPE4) ((AC(IORDER(I,JN),IORDER(J,IN)),J=1,MCHAN),I=1, A NCHAN), ((BLC(IORDER(I,JN),IORDER(J,IN)),J=1,MCHAN),I=1, B NCHAN), ((BVC(IORDER(I,JN),IORDER(J,IN)),J=1,MCHAN),I=1,NCHAN) C IF (NBUG7.EQ.0) GOTO 80 WRITE (6,3050) DO 50 I = 1,NCHAN II = IORDER(I,JN) WRITE (6,3030) II, (AC(II,IORDER(J,IN)),J=1,MCHAN) 50 CONTINUE WRITE (6,3060) DO 60 I = 1,NCHAN II = IORDER(I,JN) WRITE (6,3030) II, (BLC(II,IORDER(J,IN)),J=1,MCHAN) 60 CONTINUE WRITE (6,3070) DO 70 I = 1,NCHAN II = IORDER(I,JN) WRITE (6,3030) II, (BVC(II,IORDER(J,IN)),J=1,MCHAN) 70 CONTINUE C 80 IF(ISPN.EQ.0)THEN READ (ITAPE,END=100,ERR=90)((AC(I,J),J=1,MCHAN),I=1,NCHAN) WRITE (ITAPE4)((AC(IORDER(I,JN),IORDER(J,IN)),J=1,MCHAN), A I=1,NCHAN) GO TO 110 90 BACKSPACE (ITAPE) 100 WRITE (6,3080) ENDIF C C OUT ENDFILE(ITAPE4) 110 CLOSE (ITAPE4) C --- *=*=*=*=*=*=*=* C 120 CONTINUE C 3000 FORMAT (/35X,'SUBROUTINE DMAT'/35X,15 ('-')) 3010 FORMAT (' IT IS IMPOSSIBLE TO CREATE D-FILES BEYOND D99.') 3020 FORMAT (/' DMAT HAS NOT ENCOUNTERED ANY PERMITTED TRANSITION!') 3030 FORMAT (I5, (T6,5E15.7)) 3040 FORMAT (' INITIAL STATE NUMBER',I3,6X,'FINAL STATE NUMBER',I3) 3050 FORMAT (/' THE A COEFFICIENTS ARE') 3060 FORMAT (/' THE BL COEFFICIENTS ARE') 3070 FORMAT (/' THE BV COEFFICIENTS ARE') 3080 FORMAT (/' NO ASYMPTOTIC BP-VELOCITY ARRAY FOUND - SKIPPING') END SUBROUTINE DMAT1(IN,JN,NRANG2,NOTERM,MMNP2,MNP2,LRECW3) IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C WRITES 'D FILE' FOR INTERMEDIATE COUPLING ONTO ITAPE4, C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) PARAMETER (MXTRI=(MXMAT*(MXMAT+1))/2) PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) PARAMETER (MXDMAT=4*MXMAT*MZCHF+2*MXMAT*MXDM+2*MZNR2*MXDM, A MXRSCT=MXTRI+MZCHF*MXMAT+MZNR2*MXMAT) PARAMETER (MXM1=MXDMAT/MXRSCT,MXM2=MXRSCT/MXDMAT,MXM3=MXM1+MXM2, A MXMATX=MXDMAT*MXM1/MXM3+MXRSCT*MXM2/MXM3+1) PARAMETER (MXDUM1=MXMATX-MXDMAT,MXPRE=1+MXMAT/MZNR2) C LOGICAL NEWFIL C DIMENSION IPRE(MXPRE) C COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /MATRIX/ABUTL(MXMAT,MZCHF),ABUTV(MXMAT,MZCHF), A BBUTL(MXMAT,MZCHF),BBUTV(MXMAT,MZCHF),DML(MXMAT,MXDM), B DMV(MXMAT,MXDM),DKL(MZNR2,MXDM),DKV(MZNR2,MXDM), C DUM1(MXDUM1) COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RECORD/ILSPI(4,MZSLP),IORDER(MZCHF,MZSLP) COMMON /SPACES/XX(MZNR2,MXMAT) C SAVE NEWFIL C DATA NEWFIL/.TRUE./ C----------------------------------------------------------------------- C C INITIALIZE. OPEN SCRATCH FILE IF NOT ALREADY DONE (NEWFIL) C MCFF = MOD(MMNP2-1,NOTERM) + 1 MCHAN = ILSPI(3,IN) MCONHP = NOTERM*MCHAN MCFGP = MMNP2 - MCONHP NCFF = MOD(MNP2-1,NOTERM) + 1 NCHAN = ILSPI(3,JN) NCONHP = NOTERM*NCHAN NCFGP = MNP2 - NCONHP NDIF = NRANG2 - NOTERM M = MAX(NCHAN,MCHAN) IF (M.GT.MXDM) CALL RECOV2('DMAT1 ','MXDM ',MXDM,M) IF (NEWFIL) THEN OPEN (UNIT=IDISC3,ACCESS='DIRECT',STATUS='SCRATCH', A FORM='UNFORMATTED',RECL=8*LRECW3) NEWFIL = .FALSE. ENDIF C IREC3 = 1 ITAPE = ITAPE1 C C------MULTIPLY BY EIGENVECTORS OF LHS HAMILTONIAN C LAST = NOTERM ML = NDIF 10 CONTINUE NL = 0 MK = ML + 1 - NDIF ML = ML + NOTERM C C READ IN C-C AND C-B MATRIX ELEMENTS: DM(MMNP2,NRANG2) C LOOP = MCHAN IF (MCFGP.GT.0) LOOP = MCHAN + 1 DO 20 KR = 1,LOOP NK = NL + 1 NL = NL + NOTERM IF (KR.EQ.LOOP) NL = MMNP2 READ (ITAPE) ((DML(I,J),J=1,NRANG2),I=NK,NL) READ (ITAPE) ((DMV(I,J),J=1,NRANG2),I=NK,NL) 20 CONTINUE C C----- GIVE JREC THE LOCATION OF FIRST EIGENVECTOR RECORD C READ IN B-BUTTLE AND C-BUTTLE ELEMENTS: ABUT(MNP2,MCHAN) C 30 CONTINUE JREC = ILSPI(4,IN) READ (ITAPE) ((ABUTL(J,I),J=MK,ML),I=1,MCHAN) READ (ITAPE) ((ABUTV(J,I),J=MK,ML),I=1,MCHAN) C C----- READ IN A BLOCK OF EIGENVECTORS C FORM THE BLOCK (SQUARE AS LONG AS LAST=KNT) C KNT = NOTERM LOOP = 1 + (MMNP2-1)/NOTERM DO 40 KR = 1,LOOP IF (KR.EQ.LOOP) KNT = MCFF CALL DA2(1,JREC,IDISC2,MZNR2*MMNP2,XX) CALL VMUL(XX,DML,DKL,KNT,MMNP2,LAST) CALL VMUL(XX,DMV,DKV,KNT,MMNP2,LAST) IF (MK.EQ.1) IPRE(KR) = 0 WRITE (IDISC3,REC=IREC3) ((DKL(I,J),I=1,KNT),J=1,LAST), A ((DKV(I,J),I=1,KNT),J=1,LAST),IPRE(KR) IPRE(KR) = IREC3 IREC3 = (KNT*LAST*2)/LRECW3 + IREC3 + 1 40 CONTINUE C C-----GO BACK UP TO READ MORE DIPOLE MATRIX ELEMENTS (CYCLING NCHAN X) C IF (MK-1+NOTERM.LT.NCONHP) GOTO 10 IF (MK-1+NOTERM.GT.NCONHP) GOTO 80 C C READ IN B-C AND B-B MATRIX ELEMENTS: DM(MMNP2,NCFGP) C IRE = NCHAN IF (NCFGP.EQ.0) GOTO 80 IRE = IRE + 1 LOW = 1 LUP = NCFGP NL = 0 50 CONTINUE NK = NL + 1 NL = NL + NOTERM 60 CONTINUE READ (ITAPE) ((DML(J,I),I=LOW,LUP),J=NK,NL) READ (ITAPE) ((DMV(J,I),I=LOW,LUP),J=NK,NL) IF (NL.LT.MCONHP) GOTO 50 IF (NL.EQ.MCONHP) THEN IF (MCFGP.EQ.0) GOTO 70 NK = MCONHP + 1 NL = MMNP2 LUP = 0 ENDIF C LOW = LUP + 1 LUP = MIN(LUP+NRANG2,NCFGP) IF (LOW.LE.NCFGP) GOTO 60 C C------READ IN SOME MORE MATRIX ELEMENTS INVOLVING BUTTLE CORRECTION IN C INITIAL STATE C 70 CONTINUE MK = NCONHP + 1 ML = MNP2 LAST = NCFGP GOTO 30 C C------NOW MULTIPLY ELEMENTS INVOLVING INITIAL STATE BUTTLE CORRECTION C BY EIGENVECTORS OF THE RHS FINAL STATE HAMILTONIAN C 80 CONTINUE IREC = ILSPI(4,JN) C C------READ IN A BLOCK OF EIGENVECTORS C FORM THE PROCESSED BLOCK C IG = 0 KNT = NOTERM LOOP = 1 + (MNP2-1)/NOTERM DO 110 KR = 1,LOOP IF (KR.EQ.LOOP) KNT = NCFF CALL DA2(1,IREC,IDISC2,MZNR2*MNP2,XX) CALL VMUL(XX,ABUTL,DKL,KNT,MNP2,MCHAN) CALL VMUL(XX,ABUTV,DKV,KNT,MNP2,MCHAN) DO 100 J = 1,MCHAN DO 90 I = 1,KNT BBUTL(I+IG,J) = DKL(I,J) BBUTV(I+IG,J) = DKV(I,J) 90 CONTINUE 100 CONTINUE IG = IG + NOTERM 110 CONTINUE C C------NOW READ IN AND PROCESS WITH EIGENVECTORS OF THE LHS INITIAL C STATE HAMILTONIAN THE ELEMENTS INVOLVING BUTTLE CORRECTION C TO THE FINAL STATE C NL = 0 JREC = ILSPI(4,IN) C C READ BUTTLE-C AND BUTTLE-B ELEMENTS: ABUT(MMNP2,NCHAN) C C LOOP = 1 + (MMNP2-1)/NOTERM LOOP = MCHAN IF(MCFGP.GT.0) LOOP = LOOP + 1 DO 120 KR = 1,LOOP NK = NL + 1 NL = NL + NOTERM IF (KR.EQ.LOOP) NL = MMNP2 READ (ITAPE) ((DML(I,J),J=1,NCHAN),I=NK,NL) READ (ITAPE) ((DMV(I,J),J=1,NCHAN),I=NK,NL) 120 CONTINUE C C------READ IN A BLOCK OF EIGENVECTORS C FORM THE BLOCK C IG = 0 KNT = NOTERM DO 150 KR = 1,LOOP IF (KR.EQ.LOOP) KNT = MCFF CALL DA2(1,JREC,IDISC2,MZNR2*MMNP2,XX) CALL VMUL(XX,DML,DKL,KNT,MMNP2,NCHAN) CALL VMUL(XX,DMV,DKV,KNT,MMNP2,NCHAN) DO 140 J = 1,NCHAN DO 130 I = 1,KNT ABUTL(I+IG,J) = DKL(I,J) ABUTV(I+IG,J) = DKV(I,J) 130 CONTINUE 140 CONTINUE IG = IG + NOTERM 150 CONTINUE C C------NOW MULTIPLY BY EIGENVECTORS OF RHS HAMILTONIAN C LAST = NOTERM KR = (MMNP2-1)/NOTERM + 1 DO 200 IK = 1,KR IREC = ILSPI(4,JN) IPR = IPRE(IK) LUP = MNP2 LOW = NCONHP + 1 IF (NCFGP.EQ.0) LOW = LUP + 1 - NOTERM C C SPECIAL TREATMENT FOR LAST PART C TREAT THE LAST BLOCK OF DIPOLE MATRIX ELEMENTS: C IF (IK.EQ.KR) LAST = MCFF DO 160 JK = 1,IRE READ (IDISC3,REC=IPR) ((DML(I,J),J=1,LAST),I=LOW,LUP), A ((DMV(I,J),J=1,LAST),I=LOW,LUP),IPR LUP = LOW - 1 LOW = LOW - NOTERM 160 CONTINUE C C----- READ IN A BLOCK OF EIGENVECTORS C FORM A SQUARE PROCESSED BLOCK C KNT = NOTERM LOOP = 1 + (MNP2-1)/NOTERM DO 190 JK = 1,LOOP IF (JK.EQ.LOOP) KNT = NCFF CALL DA2(1,IREC,IDISC2,MZNR2*MNP2,XX) CALL VMUL(XX,DML,DKL,KNT,MNP2,LAST) CALL VMUL(XX,DMV,DKV,KNT,MNP2,LAST) WRITE (ITAPE4) ((DKL(J,I),J=1,KNT),I=1,LAST), A ((DKV(J,I),J=1,KNT),I=1,LAST) C IF (NBUG7.LT.1) GOTO 190 WRITE (IWRITE,*) ' LENGTH ',IK,JK DO 170 I = 1,LAST WRITE (IWRITE,3000) (DKL(J,I),J=1,KNT) 170 CONTINUE IF (NBUG7.LT.2) GOTO 190 WRITE (IWRITE,*) ' VELOCITY ',IK,JK DO 180 I = 1,LAST WRITE (IWRITE,3000) (DKV(J,I),J=1,KNT) 180 CONTINUE C 190 CONTINUE 200 CONTINUE C 3000 FORMAT (8F10.6) END SUBROUTINE DMAT2(IN,JN,NRANG2,NOTERM,MMNP2,MNP2,ABUTL,ABUTV,BBUTL, A BBUTV,DML,DMV,VI,VJ,DUM1,DUM2) IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) C DIMENSION ABUTL(MXMAT,MZCHF),ABUTV(MXMAT,MZCHF), A BBUTL(MXMAT,MZCHF),BBUTV(MXMAT,MZCHF) DIMENSION DML(MNP2,MMNP2),DMV(MNP2,MMNP2),VI(MMNP2,MMNP2), A VJ(MNP2,MNP2) DIMENSION DUM1(MXMAT),DUM2(MXMAT) C COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RECORD/ILSPI(4,MZSLP),IORDER(MZCHF,MZSLP) COMMON /SPACES/XX(MZNR2,MXMAT) C----------------------------------------------------------------------- MCHAN = ILSPI(3,IN) NCHAN = ILSPI(3,JN) ITAPE = ITAPE1 C C read in dipole matrix DML,DMV from stg2 C read in RHS Buttle correction into ABUTL,ABUTV C N2 = 0 DO 20 NCH = 1,NCHAN N1 = N2 + 1 N22 = N2 + NRANG2 M2 = 0 DO 10 MCH = 1,MCHAN READ (ITAPE) ((DML(N,M),N=N1,N22),M=M2+1,M2+NRANG2) READ (ITAPE) ((DMV(N,M),N=N1,N22),M=M2+1,M2+NRANG2) M2 = M2 + NOTERM 10 CONTINUE IF (M2.LT.MMNP2) THEN READ (ITAPE) ((DML(N,M),N=N1,N22),M=M2+1,MMNP2) READ (ITAPE) ((DMV(N,M),N=N1,N22),M=M2+1,MMNP2) ENDIF C READ (ITAPE) ((ABUTL(N,M),N=N1,N22),M=1,MCHAN) READ (ITAPE) ((ABUTV(N,M),N=N1,N22),M=1,MCHAN) N2 = N2 + NOTERM 20 CONTINUE IF (N2.LT.MNP2) THEN N1 = N2 + 1 M2 = 0 DO 30 MCH = 1,MCHAN READ (ITAPE) ((DML(N,M),N=N1,MNP2),M=M2+1,M2+NRANG2) READ (ITAPE) ((DMV(N,M),N=N1,MNP2),M=M2+1,M2+NRANG2) M2 = M2 + NOTERM 30 CONTINUE IF (M2.LT.MMNP2) THEN NDIMEN = NRANG2 NTIMES = (MNP2-N2-1)/NDIMEN + 1 DO 40 NT = 1,NTIMES N22 = MIN(N2+NDIMEN,MNP2) READ (ITAPE) ((DML(N,M),N=N2+1,N22),M=M2+1,MMNP2) READ (ITAPE) ((DMV(N,M),N=N2+1,N22),M=M2+1,MMNP2) N2 = N22 40 CONTINUE ENDIF C READ (ITAPE) ((ABUTL(N,M),N=N1,MNP2),M=1,MCHAN) READ (ITAPE) ((ABUTV(N,M),N=N1,MNP2),M=1,MCHAN) ENDIF C C read in RHS Hamiltonian eigenvectors VJ, IF NOT ALREADY THERE C JREC = ILSPI(4,JN) IF (JREC.GT.0) THEN DO 60 NP = 1,MNP2 KOUNT = MOD(NP-1,NOTERM) + 1 IF (KOUNT.EQ.1) THEN CALL DA2(1,JREC,IDISC2,MZNR2*MNP2,XX) ENDIF C DO 50 N = 1,MNP2 VJ(NP,N) = XX(KOUNT,N) 50 CONTINUE 60 CONTINUE ENDIF C C multiply RHS V by dipole matrix, overwriting DML and DMV C DO 110 M = 1,MMNP2 DO 70 NP = 1,MNP2 DUM1(NP) = 0.0D0 DUM2(NP) = 0.0D0 70 CONTINUE DO 90 N = 1,MNP2 DO 80 NP = 1,MNP2 DUM1(NP) = DUM1(NP) + DML(N,M)*VJ(NP,N) DUM2(NP) = DUM2(NP) + DMV(N,M)*VJ(NP,N) 80 CONTINUE 90 CONTINUE DO 100 NP = 1,MNP2 DML(NP,M) = DUM1(NP) DMV(NP,M) = DUM2(NP) 100 CONTINUE 110 CONTINUE C C multiply RHS V by Buttle correction, overwriting ABUTL and ABUTV C DO 160 M = 1,MCHAN DO 120 NP = 1,MNP2 DUM1(NP) = 0.0D0 DUM2(NP) = 0.0D0 120 CONTINUE DO 140 N = 1,MNP2 DO 130 NP = 1,MNP2 DUM1(NP) = DUM1(NP) + ABUTL(N,M)*VJ(NP,N) DUM2(NP) = DUM2(NP) + ABUTV(N,M)*VJ(NP,N) 130 CONTINUE 140 CONTINUE DO 150 NP = 1,MNP2 ABUTL(NP,M) = DUM1(NP) ABUTV(NP,M) = DUM2(NP) 150 CONTINUE 160 CONTINUE C C read in LHS Hamiltonian eigenvectors VI, IF NOT ALREADY THERE C IREC = ILSPI(4,IN) IF (IREC.GT.0) THEN DO 180 MP = 1,MMNP2 KOUNT = MOD(MP-1,NOTERM) + 1 IF (KOUNT.EQ.1) THEN CALL DA2(1,IREC,IDISC2,MZNR2*MMNP2,XX) ENDIF C DO 170 M = 1,MMNP2 VI(MP,M) = XX(KOUNT,M) 170 CONTINUE 180 CONTINUE ENDIF C C multiply LHS V by dipole matrix, overwriting DML and DMV C DO 230 NP = 1,MNP2 DO 190 MP = 1,MMNP2 DUM1(MP) = 0.0D0 DUM2(MP) = 0.0D0 190 CONTINUE DO 210 M = 1,MMNP2 DO 200 MP = 1,MMNP2 DUM1(MP) = DUM1(MP) + VI(MP,M)*DML(NP,M) DUM2(MP) = DUM2(MP) + VI(MP,M)*DMV(NP,M) 200 CONTINUE 210 CONTINUE DO 220 MP = 1,MMNP2 DML(NP,MP) = DUM1(MP) DMV(NP,MP) = DUM2(MP) 220 CONTINUE 230 CONTINUE C C read in LHS Buttle correction into BBUTL,BBUTV C M2 = 0 DO 240 MCH = 1,MCHAN READ (ITAPE) ((BBUTL(M,N),N=1,NCHAN),M=M2+1,M2+NRANG2) READ (ITAPE) ((BBUTV(M,N),N=1,NCHAN),M=M2+1,M2+NRANG2) M2 = M2 + NOTERM 240 CONTINUE IF (M2.LT.MMNP2) THEN READ (ITAPE) ((BBUTL(M,N),N=1,NCHAN),M=M2+1,MMNP2) READ (ITAPE) ((BBUTV(M,N),N=1,NCHAN),M=M2+1,MMNP2) ENDIF C C multiply LHS V by Buttle correction, overwriting BBUTL and BBUTV C DO 290 N = 1,NCHAN DO 250 MP = 1,MMNP2 DUM1(MP) = 0.0D0 DUM2(MP) = 0.0D0 250 CONTINUE DO 270 M = 1,MMNP2 DO 260 MP = 1,MMNP2 DUM1(MP) = DUM1(MP) + VI(MP,M)*BBUTL(M,N) DUM2(MP) = DUM2(MP) + VI(MP,M)*BBUTV(M,N) 260 CONTINUE 270 CONTINUE DO 280 MP = 1,MMNP2 BBUTL(MP,N) = DUM1(MP) BBUTV(MP,N) = DUM2(MP) 280 CONTINUE 290 CONTINUE C C WRITE OUT DIPOLE MATRICES C N2 = 0 DO 330 IK = 1,1 + (MMNP2-1)/NOTERM N1 = N2 + 1 N2 = MIN(N2+NOTERM,MMNP2) M2 = 0 DO 320 JK = 1,1 + (MNP2-1)/NOTERM M1 = M2 + 1 M2 = MIN(M2+NOTERM,MNP2) WRITE (ITAPE4) ((DML(M,N),M=M1,M2),N=N1,N2), A ((DMV(M,N),M=M1,M2),N=N1,N2) C IF (NBUG7.LT.1) GOTO 320 WRITE (IWRITE,*) ' LENGTH ',IK,JK DO 300 N = N1,N2 WRITE (IWRITE,3000) (DML(M,N),M=M1,M2) 300 CONTINUE IF (NBUG7.LT.2) GOTO 320 WRITE (IWRITE,*) ' VELOCITY ',IK,JK DO 310 N = N1,N2 WRITE (IWRITE,3000) (DMV(M,N),M=M1,M2) 310 CONTINUE C 320 CONTINUE 330 CONTINUE C 3000 FORMAT (8F10.6) END SUBROUTINE DMAT2P(IN,JN,NRANG2,NOTERM,MMNP2,MNP2,ABUTL,ABUTV,BBUTL A ,BBUTV,DML,DMV,VI,VJ,DUM1,DUM2,mnp22) IMPLICIT REAL*8 (A-H,O-Z) C C----------------------------------------------------------------------- C C Read dipole matrices and associated Buttle corrections and multiply C by Hamiltonian e-vectors. C C Optionally, transform and reduce so as to eliminate pseudoresonances. C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) parameter (nsizex=3) !min 1 C CHARACTER*1 IFNM(10) CHARACTER*12 IFLNM C DIMENSION ABUTL(MXMAT,MZCHF),ABUTV(MXMAT,MZCHF), A BBUTL(MXMAT,MZCHF),BBUTV(MXMAT,MZCHF) DIMENSION DML(mnp22,MMNP2),DMV(mnp22,MMNP2),VI(MMNP2,MMNP2), A VJ(MNP2,MNP2) DIMENSION DUM1(MXMAT),DUM2(MXMAT) C common/adip/ad(mznc2,mznc2,nsizex),nladq(mzslp),nradq(mzslp), c ixsym(mzslp),nxsym COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /HAMNUW/A(MZNC2,MZNC2),XO(MZNC2,MZNC2),XL(MZNC2) X ,WORK(MZNC2) COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RECORD/ILSPI(4,MZSLP),IORDER(MZCHF,MZSLP) COMMON /SPACES/XX(MZNR2,MXMAT) C DATA IFNM/'0','1','2','3','4','5','6','7','8','9'/ C----------------------------------------------------------------------- MCHAN = ILSPI(3,IN) NCHAN = ILSPI(3,JN) ITAPE = ITAPE1 C C read in dipole matrix DML,DMV from stg2 C read in RHS Buttle correction into ABUTL,ABUTV C N2 = 0 DO 20 NCH = 1,NCHAN N1 = N2 + 1 N22 = N2 + NRANG2 M2 = 0 DO 10 MCH = 1,MCHAN READ (ITAPE) ((DML(N,M),N=N1,N22),M=M2+1,M2+NRANG2) READ (ITAPE) ((DMV(N,M),N=N1,N22),M=M2+1,M2+NRANG2) M2 = M2 + NOTERM 10 CONTINUE IF (M2.LT.MMNP2) THEN READ (ITAPE) ((DML(N,M),N=N1,N22),M=M2+1,MMNP2) READ (ITAPE) ((DMV(N,M),N=N1,N22),M=M2+1,MMNP2) ENDIF C READ (ITAPE) ((ABUTL(N,M),N=N1,N22),M=1,MCHAN) READ (ITAPE) ((ABUTV(N,M),N=N1,N22),M=1,MCHAN) N2 = N2 + NOTERM 20 CONTINUE IF (N2.LT.mnp22) THEN N1 = N2 + 1 M2 = 0 DO 30 MCH = 1,MCHAN READ (ITAPE) ((DML(N,M),N=N1,mnp22),M=M2+1,M2+NRANG2) READ (ITAPE) ((DMV(N,M),N=N1,mnp22),M=M2+1,M2+NRANG2) M2 = M2 + NOTERM 30 CONTINUE IF (M2.LT.MMNP2) THEN NDIMEN = NRANG2 NTIMES = (mnp22-N2-1)/NDIMEN + 1 DO 40 NT = 1,NTIMES N22 = MIN(N2+NDIMEN,mnp22) READ (ITAPE) ((DML(N,M),N=N2+1,N22),M=M2+1,MMNP2) READ (ITAPE) ((DMV(N,M),N=N2+1,N22),M=M2+1,MMNP2) N2 = N22 40 CONTINUE ENDIF C READ (ITAPE) ((ABUTL(N,M),N=N1,mnp22),M=1,MCHAN) READ (ITAPE) ((ABUTV(N,M),N=N1,mnp22),M=1,MCHAN) ENDIF C c c Now enter pseudoresonance elimination section (bail-out if none) c jnx=ixsym(jn) if(jnx.eq.0)go to 45 !No reduction necessary c nrad=nradq(jnx) nlad=nladq(jnx) c if(nrad.lt.nlad)then !should never happen... write(6,*)'Error: nrad.lt.nlad:',nrad,nlad write(0,*)'Error: nrad.lt.nlad:',nrad,nlad stop endif c c write(0,*) 'in,jn,nrad(jn),nlad(jn) = ',in,jn,nrad,nlad c c Load appropriate transformation symmetry, from disc if necessary c if(jnx.le.nsizex)then do j=1,nrad do i=1,nlad a(i,j)=ad(i,j,jnx) enddo enddo else iu=jnx-nsizex iflnm='TEMPU00'//ifnm(iu)//'.TMP' open (unit=76,file=iflnm,status='old',form='unformatted') read(76)((a(i,j),i=1,nlad),j=1,nrad) close(76) endif c c N.R.B. mnp22=mnp2+nrad-nlad c i0=mnp22+1-nrad !=mnp2+1-nlad ii0=nrad-mnp22 !=nlad-mnp2 c c transform and reduce dipole length c do j=1,mmnp2 do i=i0,mnp22 ii=ii0+i !so ii=1,nrad work(ii)=dml(i,j) dml(i,j)=0.0d0 enddo do k=1,nrad t=work(k) do i=i0,mnp2 !=mnp22+nlad-nrad ii=ii0+i !so ii=1,nlad dml(i,j)=dml(i,j)+a(ii,k)*t enddo enddo enddo c c transform and reduce buttle length c do j=1,mchan do i=i0,mnp22 !mnp22 - not minp=mnp2 ii=ii0+i work(ii)=abutl(i,j) abutl(i,j)=0.0d0 enddo do k=1,nrad t=work(k) do i=i0,mnp2 !=mnp22+nlad-nrad ii=ii0+i abutl(i,j)=abutl(i,j)+a(ii,k)*t enddo enddo enddo c c transform and reduce dipole velocity c do j=1,mmnp2 do i=i0,mnp22 ii=ii0+i !so ii=1,nrad work(ii)=dmv(i,j) dmv(i,j)=0.0d0 enddo do k=1,nrad t=work(k) do i=i0,mnp2 !=mnp22+nlad-nrad ii=ii0+i dmv(i,j)=dmv(i,j)+a(ii,k)*t enddo enddo enddo c c transform and reduce buttle velocity c do j=1,mchan do i=i0,mnp22 !mnp22 - not minp=mnp2 ii=ii0+i work(ii)=abutv(i,j) abutv(i,j)=0.0d0 enddo do k=1,nrad t=work(k) do i=i0,mnp2 !=mnp22+nlad-nrad ii=ii0+i abutv(i,j)=abutv(i,j)+a(ii,k)*t enddo enddo enddo c c re-entry point if no reduction. C C read in RHS Hamiltonian eigenvectors VJ, IF NOT ALREADY THERE C 45 JREC = ILSPI(4,JN) IF (JREC.GT.0) THEN DO 60 NP = 1,MNP2 KOUNT = MOD(NP-1,NOTERM) + 1 IF (KOUNT.EQ.1) THEN CALL DA2(1,JREC,IDISC2,MZNR2*MNP2,XX) ENDIF C DO 50 N = 1,MNP2 VJ(NP,N) = XX(KOUNT,N) 50 CONTINUE 60 CONTINUE ENDIF C C multiply RHS V by dipole matrix, overwriting DML and DMV C DO 110 M = 1,MMNP2 DO 70 NP = 1,MNP2 DUM1(NP) = 0.0D0 DUM2(NP) = 0.0D0 70 CONTINUE DO 90 N = 1,MNP2 DO 80 NP = 1,MNP2 DUM1(NP) = DUM1(NP) + DML(N,M)*VJ(NP,N) DUM2(NP) = DUM2(NP) + DMV(N,M)*VJ(NP,N) 80 CONTINUE 90 CONTINUE DO 100 NP = 1,MNP2 DML(NP,M) = DUM1(NP) DMV(NP,M) = DUM2(NP) 100 CONTINUE 110 CONTINUE C C multiply RHS V by Buttle correction, overwriting ABUTL and ABUTV C DO 160 M = 1,MCHAN DO 120 NP = 1,MNP2 DUM1(NP) = 0.0D0 DUM2(NP) = 0.0D0 120 CONTINUE DO 140 N = 1,MNP2 DO 130 NP = 1,MNP2 DUM1(NP) = DUM1(NP) + ABUTL(N,M)*VJ(NP,N) DUM2(NP) = DUM2(NP) + ABUTV(N,M)*VJ(NP,N) 130 CONTINUE 140 CONTINUE DO 150 NP = 1,MNP2 ABUTL(NP,M) = DUM1(NP) ABUTV(NP,M) = DUM2(NP) 150 CONTINUE 160 CONTINUE C C read in LHS Hamiltonian eigenvectors VI, IF NOT ALREADY THERE C IREC = ILSPI(4,IN) IF (IREC.GT.0) THEN DO 180 MP = 1,MMNP2 KOUNT = MOD(MP-1,NOTERM) + 1 IF (KOUNT.EQ.1) THEN CALL DA2(1,IREC,IDISC2,MZNR2*MMNP2,XX) ENDIF C DO 170 M = 1,MMNP2 VI(MP,M) = XX(KOUNT,M) 170 CONTINUE 180 CONTINUE ENDIF C C multiply LHS V by dipole matrix, overwriting DML and DMV C DO 230 NP = 1,MNP2 DO 190 MP = 1,MMNP2 DUM1(MP) = 0.0D0 DUM2(MP) = 0.0D0 190 CONTINUE DO 210 M = 1,MMNP2 DO 200 MP = 1,MMNP2 DUM1(MP) = DUM1(MP) + VI(MP,M)*DML(NP,M) DUM2(MP) = DUM2(MP) + VI(MP,M)*DMV(NP,M) 200 CONTINUE 210 CONTINUE DO 220 MP = 1,MMNP2 DML(NP,MP) = DUM1(MP) DMV(NP,MP) = DUM2(MP) 220 CONTINUE 230 CONTINUE C C read in LHS Buttle correction into BBUTL,BBUTV C M2 = 0 DO 240 MCH = 1,MCHAN READ (ITAPE) ((BBUTL(M,N),N=1,NCHAN),M=M2+1,M2+NRANG2) READ (ITAPE) ((BBUTV(M,N),N=1,NCHAN),M=M2+1,M2+NRANG2) M2 = M2 + NOTERM 240 CONTINUE IF (M2.LT.MMNP2) THEN READ (ITAPE) ((BBUTL(M,N),N=1,NCHAN),M=M2+1,MMNP2) READ (ITAPE) ((BBUTV(M,N),N=1,NCHAN),M=M2+1,MMNP2) ENDIF C C multiply LHS V by Buttle correction, overwriting BBUTL and BBUTV C DO 290 N = 1,NCHAN DO 250 MP = 1,MMNP2 DUM1(MP) = 0.0D0 DUM2(MP) = 0.0D0 250 CONTINUE DO 270 M = 1,MMNP2 DO 260 MP = 1,MMNP2 DUM1(MP) = DUM1(MP) + VI(MP,M)*BBUTL(M,N) DUM2(MP) = DUM2(MP) + VI(MP,M)*BBUTV(M,N) 260 CONTINUE 270 CONTINUE DO 280 MP = 1,MMNP2 BBUTL(MP,N) = DUM1(MP) BBUTV(MP,N) = DUM2(MP) 280 CONTINUE 290 CONTINUE C C WRITE OUT DIPOLE MATRICES C N2 = 0 DO 330 IK = 1,1 + (MMNP2-1)/NOTERM N1 = N2 + 1 N2 = MIN(N2+NOTERM,MMNP2) M2 = 0 DO 320 JK = 1,1 + (MNP2-1)/NOTERM M1 = M2 + 1 M2 = MIN(M2+NOTERM,MNP2) WRITE (ITAPE4) ((DML(M,N),M=M1,M2),N=N1,N2), A ((DMV(M,N),M=M1,M2),N=N1,N2) C IF (NBUG7.LT.1) GOTO 320 WRITE (IWRITE,*) ' LENGTH ',IK,JK DO 300 N = N1,N2 WRITE (IWRITE,3000) (DML(M,N),M=M1,M2) 300 CONTINUE IF (NBUG7.LT.2) GOTO 320 WRITE (IWRITE,*) ' VELOCITY ',IK,JK DO 310 N = N1,N2 WRITE (IWRITE,3000) (DMV(M,N),M=M1,M2) 310 CONTINUE C 320 CONTINUE 330 CONTINUE C 3000 FORMAT (8F10.6) END C********************************************************************** SUBROUTINE HAMNEW(HNP1,MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) C IMPLICIT REAL*8 (A-H,O-Z) C C----------------------------------------------------------------------- C C THIS ROUTINE READS IN THE OVERLAP MATRIX FROM STG2 C AND TRANSFORMS THE HAMILTONIAN, PERHAPS REDUCING THE NUMBER OF C (N+1)-ELECTRON CONFIGURATIONS - TWG & NRB C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C PARAMETER (MXCC=MZCHF*MZNR2) PARAMETER (MXBB=(MZNC2*(MZNC2+1))/2) parameter (nsizex=3) !min 1 C LOGICAL EXA CHARACTER*1 IFNM(10) CHARACTER*12 IFLNM C common/adip/ad(mznc2,mznc2,nsizex),nlad(mzslp),nrad(mzslp), c ixsym(mzslp),nxsym COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST COMMON /REDTOL/TOLER,iparinit COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /HAMNUW/A(MZNC2,MZNC2),XO(MZNC2,MZNC2),XL(MZNC2) X ,WORK(MZNC2) COMMON /INFORM/IREAD,IWRITE,IPUNCH C DIMENSION HNP1(*),HBC(MZNC2,MXCC),HBB(MXBB),B(MZNC2,MZNC2) C SAVE IFL DATA IFL/1/,nxsym/0/ DATA IFNM/'0','1','2','3','4','5','6','7','8','9'/ C----------------------------------------------------------------------- C NPOS(MNP,I,J)=(MNP*(MNP+1))/2-((MNP-I+1)*(MNP-I+2))/2+J-I+1 !I.LE.J C C RETURN (DO NOTHING TO HAMILTONIAN) IF TOLERANCE IS REAL SMALL C IF(TOLER.LT.1.D-19)RETURN C C c 1 READ(77,*,END=3)NSPIN,NLAM,NPAR,NL,NR,NONZER 1 READ(77,END=3)NSPIN,NLAM,NPAR,NL,NR,NONZER C ireturn=0 write(6,*)' ' if(npar.eq.iparinit)then write(6,*)' >>>>>>>>>>>>>>>>>>>>>NOT REDUCING THIS PARTIAL WAVE' c write(0,*)' not reducing this partial wave' ireturn=1 if(mskip.gt.1)then write(6,*)' Error, must have initial parity first in stg2' stop ' Error, must have initial parity first in stg2' endif else write(6,*)' >>>>>>>>>>>>>>>>>>>>>>>>>REDUCING THIS PARTIAL WAVE' c write(0,*)' reducing this partial wave' endif C IF(NSPIN.NE.NSPN2.OR.NLAM.NE.LRGL2.OR.NPAR.NE.NPTY2. C OR.NR.NE.NCFGP2)THEN DO N=1,NONZER c READ(77,*)I,J,DUM READ(77)I,J,DUM ENDDO GO TO 1 ENDIF C IF(NONZER.EQ.0)RETURN CHK IF(MNP1-NCFGP2.GT.MXCC)THEN WRITE(IWRITE,*)'SR.HAMNEW: INCREASE MXCC TO ',MNP1-NCFGP2 ENDIF IF(NCFGP2.GT.MZNC2)THEN WRITE(IWRITE,*)'SR.HAMNEW: INCREASE MZNC2 TO ',NCFGP2 STOP 'SR.HAMNEW: INCREASE MZNC2' ENDIF IF(NL.GT.MZNC2.OR.NR.GT.MZNC2)THEN WRITE(IWRITE,*) X 'SR.HAMNEW: INCREASE MZNC2 TO MAX OF ',NL,NR STOP 'SR.HAMNEW: INCREASE MZNC2' ENDIF CHK C C INITIALIZE A C DO J=1,NR DO I=1,NL A(I,J)=0.0D0 ENDDO ENDDO C C READ IN NON-ZERO ELEMENTS OF A-MATRIX FROM STG2 C DO NZ=1,NONZER c READ(77,*)I,J,ATEMP READ(77)I,J,ATEMP A(I,J)=A(I,J)+ATEMP ENDDO c if(ireturn.ne.0)return C C MAKE XO=A*AT C DO I=1,NL DO J=1,NR B(J,I)=A(I,J) ENDDO ENDDO DO J=1,NL DO I=1,NL XO(I,J)=0.0D0 ENDDO DO K=1,NR DO I=1,NL XO(I,J)=XO(I,J)+A(I,K)*B(K,J) ENDDO ENDDO ENDDO C C MAKE O AND L=OT*A*AT*O (DIAGONAL EIGENVALUES) C CALL DIAG(NL,-1,XO,XL,WORK,MZNC2) C IF(NBUG9.GT.0)THEN WRITE(78,*)' EIGENVALUES FOR PARTIAL WAVE ',NSPN2,LRGL2,NPTY2 DO I=1,NL WRITE(78,100)XL(I) 100 FORMAT(6E13.6) ENDDO ENDIF C DO I=1,NL XL(I)=ABS(XL(I)) IF(XL(I).LT.1.D-18)XL(I)=1.D-18 ENDDO C C MAKE B=OT*A C DO J=1,NR DO I=1,NL SUM=0.0D0 DO K=1,NL SUM=SUM+XO(K,I)*A(K,J) ENDDO B(I,J)=SUM/SQRT(XL(I)) ENDDO ENDDO C C MAKE ANEW=(1/SQRT(L))*OT*A C DO J=1,NR DO I=1,NL A(I,J)=B(I,J) ENDDO ENDDO C C DETERMINE HOW MANY BASIS FUNCTIONS WILL BE KEPT C IF (NL.GT.NR)NL=NR NLNEW=NL DO I=1,NL IF(ABS(XL(NL-I+1)).LT.TOLER)THEN NLNEW=NLNEW-1 ELSE GO TO 51 ENDIF ENDDO 51 NL=NLNEW C C WRITE NEW MATRIX C IF(NBUG9.GT.10)THEN WRITE(78,*)' NEW A-MATRIX OF DIMENSION ',NL,'X',NR DO J=1,NR WRITE(78,100)(A(I,J),I=1,NL) ENDDO ENDIF c c save for dipole transformation c if(ipolph.eq.2)then nxsym=nxsym+1 ixsym(mskip)=nxsym c write(0,*)' nxsym = ',nxsym nlad(nxsym)=nl nrad(nxsym)=nr if(nxsym.le.nsizex)then ! store in memory do j=1,nr do i=1,nl ad(i,j,nxsym)=a(i,j) enddo enddo else ! write to temp file write(6,*)'SR.HAMNEW: WRITING THIS SYMMETRY TO DISC FOR DMAT2' iu=nxsym-nsizex iflnm='TEMPU00'//ifnm(iu)//'.TMP' open (unit=76,file=iflnm,status='unknown',form='unformatted') write(76)((a(i,j),i=1,nl),j=1,nr) close(76) endif endif C C TRANSFORM (AND REDUCE) HAMILTONIAN C C HCC -> HCC C HBC -> A*HBC C HCB -> HCB*AT C HBB -> A*HBB*AT C NCONT=MNP1-NCFGP2 MNPN=NCONT+NL C DO I=1,NCONT DO J=1,NCFGP2 JJ=J+NCONT HBC(J,I)=HNP1(NPOS(MNP1,I,JJ)) ENDDO ENDDO DO I=1,NCFGP2 II=I+NCONT DO J=I,NCFGP2 JJ=J+NCONT HBB(NPOS(NCFGP2,I,J))=HNP1(NPOS(MNP1,II,JJ)) ENDDO ENDDO C C C-C C DO I=1,NCONT DO J=I,NCONT HNP1(NPOS(MNPN,I,J))=HNP1(NPOS(MNP1,I,J)) ENDDO ENDDO C C B-C C DO I=1,NCONT DO J=1,NL WORK(J)=0.0D0 ENDDO DO K=1,NCFGP2 DO J=1,NL WORK(J)=WORK(J)+A(J,K)*HBC(K,I) ENDDO ENDDO DO J=1,NL JJ=J+NCONT HNP1(NPOS(MNPN,I,JJ))=WORK(J) ENDDO ENDDO C C B-B C DO J=1,NL DO I=1,NCFGP2 XO(I,J)=A(J,I) ENDDO ENDDO DO J=1,NL DO I=1,NCFGP2 B(I,J)=0.0D0 DO K=1,NCFGP2 N=NPOS(NCFGP2,MIN(I,K),MAX(I,K)) B(I,J)=B(I,J)+HBB(N)*XO(K,J) ENDDO ENDDO ENDDO C DO J=1,NL DO I=J,NL WORK(I)=0.0D0 ENDDO DO K=1,NCFGP2 DO I=J,NL WORK(I)=WORK(I)+A(I,K)*B(K,J) ENDDO ENDDO DO I=J,NL HNP1(NPOS(MNPN,J+NCONT,I+NCONT))=WORK(I) ENDDO ENDDO C C RE-DEFINE MNP1, NCFGP2 C IF(NL.LT.NCFGP2)THEN WRITE(IWRITE,*) X ' REDUCING # OF (N+1)-CONFIGS FROM ',NCFGP2,' TO ',NL NCFGP2=NL MNP1=NCONT+NL ENDIF WRITE(6,*)' ' C RETURN C C OPEN NEW FILE AMATUXXX.DAT 3 IFL=IFL+1 IFLNM='AMATU00'//IFNM(IFL)//'.DAT' INQUIRE (FILE=IFLNM,EXIST=EXA) IF(EXA) THEN CLOSE(77) OPEN (UNIT=77,FILE=IFLNM,STATUS='OLD',FORM='UNFORMATTED') WRITE(IWRITE,3210) IFLNM 3210 FORMAT(' OPENING FILE ',A12/) GO TO 1 ELSE WRITE(IWRITE,*)' SR.HAMNEW: ERROR, NS,NL,NP =', C NSPN2,LRGL2,NPTY2,' CASE NOT FOUND' STOP 'SR.HAMNEW: ERROR, SLP CASE NOT FOUND' END IF C END C********************************************************************** SUBROUTINE HAMNEW2(HNP1,MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) C IMPLICIT REAL*8 (A-H,O-Z) C C----------------------------------------------------------------------- C C THIS ROUTINE SHIFTS THE EIGENVALUES OF THE BOUND-BOUND C HAMILTONIAN - TWG C TBD: MAKE THE MATRIX MULTIPLICATIONS EFFICIENT! C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST COMMON /SHIFT/ESHIFT,NSHIFT COMMON /HAMNUW/A(MZNC2,MZNC2),XO(MZNC2,MZNC2),XL(MZNC2) X ,WORK(MZNC2) C DIMENSION HNP1(*) C IF(NSHIFT.LE.0)RETURN C READ(81,*)NSPIN,NLAM,NPAR C IF(NSPIN.NE.NSPN2.OR.NLAM.NE.LRGL2.OR.NPAR.NE.NPTY2)THEN WRITE(6,*)' ERROR, DIDNT FIND CORRECT PARTIAL WAVE INFO' WRITE(6,*)NSPIN,NLAM,NPAR,' ON FILE ',NSPN2,LRGL2,NPTY2,' PROGRAM' STOP ENDIF C IF(NCFGP2.GT.MZNC2)THEN WRITE(6,*)'SR.HAMNEW2: INCREASE MZNC2 IN TO ',NCFGP2 STOP 'SR.HAMNEW2: INCREASE MZNC2' ENDIF C NCONT=MNP1-NCFGP2 K=0 DO 10 I=1,MNP1 DO 15 J=I,MNP1 K=K+1 IF(I.GT.NCONT.AND.J.GT.NCONT)THEN XO(I-NCONT,J-NCONT)=HNP1(K) XO(J-NCONT,I-NCONT)=HNP1(K) ENDIF 15 CONTINUE 10 CONTINUE C C C DIAGONALIZE THIS C CALL DIAG(NCFGP2,-1,XO,XL,WORK,MZNC2) C WRITE(79,*)' EIGENVALUES FOR PARTIAL WAVE ',NSPN2,LRGL2,NPTY2, C 'RELATIVE TO GROUND CONTINUUM (IN RYD)' NBELOW=0 DO 33 I=1,NCFGP2 EREL=2.0D0*(XL(I)-ENAT(1)) IF(EREL.LT.0.)NBELOW=NBELOW+1 WRITE(79,111)I,EREL 111 FORMAT(I5,E13.6) 33 CONTINUE C C CORRECT EIGENVALUES C READ(81,*)NUM IF(NUM.LE.0) GO TO 81 DO 40 I=1,NUM READ(81,*)NCORRECT,DELTA J=NCFGP2-(NCORRECT+NBELOW)+1 XL(J)=XL(J)+0.5D0*DELTA 40 CONTINUE C C MAKE NEW BOUND-BOUND HAMILTONIAN C DO 30 I=1,NCFGP2 DO 32 J=1,NCFGP2 SUM=0.0D0 DO 31 K=1,NCFGP2 SUM=SUM+XO(I,K)*XL(K)*XO(J,K) 31 CONTINUE A(I,J)=SUM 32 CONTINUE 30 CONTINUE C C PUT HAMILTONIAN BACK ONTO HNP1 VECTOR C K=0 DO 80 I=1,MNP1 DO 82 J=I,MNP1 K=K+1 IF(I.GT.NCONT.AND.J.GT.NCONT)THEN HNP1(K)=A(I-NCONT,J-NCONT) ENDIF 82 CONTINUE 80 CONTINUE C 81 RETURN END SUBROUTINE ORDER(EN,NORDER,NDIM,IUP) IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C RETURNS NORDER(I)=POINTER TO I-TH ENERGY IN EN ARRAY, C IUP=1 FOR ASCENDING ENERGIES, IUP=-1 FOR DESCENDING ENERGIES C C----------------------------------------------------------------------- DIMENSION EN(NDIM),NORDER(NDIM) C----------------------------------------------------------------------- DO 40 K = 1,NDIM J = K J1 = J - 1 IF (J1.EQ.0) GOTO 30 E = EN(J) + IUP*1.0D-7 DO 20 I = 1,J1 IF (J.LT.K) GOTO 10 N = NORDER(I) IF (IUP.GT.0 .AND. E.GT.EN(N)) GOTO 20 IF (IUP.LT.0 .AND. E.LT.EN(N)) GOTO 20 10 CONTINUE NORDER(J) = NORDER(J-1) J = J - 1 20 CONTINUE 30 NORDER(J) = K 40 CONTINUE C END SUBROUTINE RECOV2(SUBNAM,PARNAM,IDIM,NEED) C C C C----------------------------------------------------------------------- C C THIS ROUTINE IS CALLED ONLY FOR ARRAY OVERFLOW C C SUBNAM = SUBROUTINE NAME C PARNAM = PARAMETER NAME C IDIM = CURRENT DIMENSION C NEED = REQUIRED DIMENSION, RETURN NEED=IDIM C IPLACE = 0 TO STOP PROGRAM, C OTHERWISE RETURN IPLACE=NEED TO THE CALLING ROUTINE. C C----------------------------------------------------------------------- CHARACTER*6 SUBNAM,PARNAM,PREPRO C COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /RECOV/IPLACE C----------------------------------------------------------------------- PREPRO = PARNAM IF (PREPRO(1:3).EQ.' MZ') PREPRO(1:3) = ' MZ' WRITE (IWRITE,3000) SUBNAM,PREPRO,IDIM,PREPRO,NEED C IF (IPLACE.EQ.0) THEN WRITE (IWRITE,3010) STOP ENDIF C IF (IPLACE.LT.0) WRITE (IWRITE,3020) IPLACE = NEED NEED = IDIM C 3000 FORMAT (/' * ARRAY OVERFLOW IN ', A A6/' MUST INCREASE DIMENSION GIVEN BY ',A6,' =',I7, B ' TO AT LEAST ',A6,' =',I9) 3010 FORMAT (/' PROGRAM TERMINATES IN RECOV2'/) 3020 FORMAT (/' CHECK TO SEE IF OTHER ARRAYS ARE GOING TO BE EXCEEDED') END SUBROUTINE RSCT IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C SETS UP AND DIAGONALIZES THE HAMILTONIAN MATRIX C AND DETERMINES THE SURFACE AMPLITUDES FROM THE EIGENVECTORS. C C OPTIONALLY USES LAPACK 3.0 ROUTINES DSYEVR (IEEE O.K.) OR C DSYEVD (IEEE NOT O.K.) WHICH REQUIRE AN ADDITIONAL 2N^2 OR C 3N^2 OF MEMORY. STATEMENTS COMMON TO BOTH ARE MARKED !LAPACK C WHILE THOSE REQUIRED ONLY BY DSYEVR/D ARE MARKED !DSYEVR/D C NOTE, PSYEVD (PACKED STORAGE) IS SIGNIFICANTLY SLOWER THAN DSYEVD C ---- NRB/CPB C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C PARAMETER (MXMEM=MZMEG*1000000+MZKIL*1000+1) PARAMETER (MXN21=MZNR2+1) PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) PARAMETER (MXTRI=(MXMAT*(MXMAT+1))/2) PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) PARAMETER (MXDMAT=4*MXMAT*MZCHF+2*MXMAT*MXDM+2*MZNR2*MXDM, A MXRSCT=MXTRI+MZCHF*MXMAT+MZNR2*MXMAT) PARAMETER (MXM1=MXDMAT/MXRSCT,MXM2=MXRSCT/MXDMAT,MXM3=MXM1+MXM2, A MXMATX=MXDMAT*MXM1/MXM3+MXRSCT*MXM2/MXM3+1) PARAMETER (MXH=MXMATX-MXRSCT+MZNR2*MXMAT) C CL PARAMETER (LIWORK=10*MXMAT) !LAPACK C PARAMETER (LWORK=2*MXMAT*MXMAT+6*MXMAT+1) !DSYEVD C CHARACTER*6 FL8(8),FO8(2) CHARACTER*1 RNG C DIMENSION NORDER(MZCHF),LARGE(5) DIMENSION VALUE(MXMAT),X(MXMAT) DIMENSION AUX(MXMAT,9) !ORIG C CL COMMON /LAP1/A(MXMAT,MXMAT),XA(MXMAT,MXMAT),AUX(MXMAT,6),DUM !LAPACK CL DIMENSION ISUPP(2*MXMAT),IWORK(LIWORK) !LAPACK C DIMENSION WORK(LWORK),B(MXMAT,MXMAT) !DSYEVD C EQUIVALENCE (WORK(1),A(1,1)) !DSYEVD C COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /BASIC2/CF(MZCHF,MZCHF,MZLMX),ET(MZCHF),KJ(MZCHF), A L2P(MZCHF),LSTARG(MZCHF),NCONAT(MZTAR),LRGL,NSPN,NPTY COMMON /BASIN/EIGENS(MZNR2,MZLR2),ENDS(MXN21,MZLR2),DELTA,ETA COMMON /BUTT/COEFF(3,MZLR2),EK2MAX,EK2MIN,MAXNCB(MZLR2),NELCOR COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /MATRIX/HNP1(MXTRI),WMAT(MZCHF,MXMAT),H(MXH) COMMON /MEMORY/ARRAY(MXMEM),MEM1,MREC1 COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECORD/ILSPI(4,MZSLP),IORDER(MZCHF,MZSLP) COMMON /SPACES/XX(MZNR2,MXMAT) COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST COMMON /NRBMNP/MNP2MX,IPRCENT C DATA EPSI/1.0D-9/,ZERO/0.0D0/ DATA FL8/' (8','(14X,7','(28X,6','(42X,5','(56X,4','(70X,3', A '(84X,2',' (98X,'/,FO8(2)/'F14.7)'/ C----------------------------------------------------------------------- WRITE (IWRITE,3000) IF (NOT1.GT.NRANG2) NOT1 = NRANG2 NOTERM = NOT1 C C USE THE NUMBER OF CONTINUUM TERMS TO BE RETAINED IN C THE EXPANSION. THIS WILL REDUCE THE SIZE OF THE C HAMILTONIAN MATRIX FROM MNP1 TO MNP2. C MNP2 = NOTERM*NCHAN + NCFGP MN4 = 2*MNP2 IF (NOT1.EQ.NRANG2) GOTO 60 WRITE (IWRITE,3060) MNP2,NCHAN MN2 = 2*MNP1 C C READ THE FULL HAMILTONIAN MATRIX FROM SCRATCH DISC, C AND PUT THE UPPER TRIANGLE OF THE HAMILTONIAN MATRIX REQUIRED C INTO A ONE DIMENSIONAL ARRAY, HNP1; LOOP OVER EACH CHANNEL: C REWIND IDISC1 C DO 50 M = 1,NCHAN MO = (M-1)*NOTERM MH = MN2 - 2* (M-1)*NRANG2 ML = ((MH-NRANG2+1)*NRANG2)/2 READ (IDISC1) (H(I),I=1,ML) C C LOOP OVER THE NUMBER OF TERMS REQUIRED IN EACH CHANNEL C DO 40 I = 1,NOTERM IA = NOTERM - I + 1 IK = ((MO+I-1)* (MN4-MO-I+2))/2 IL = IK - I + 1 IG = ((I-1)* (MH-I+2))/2 IP = IG - I + 1 C C THE CONTINUUM-CONTINUUM TERMS C LO = IK LP = IG I1 = IA DO 20 L = 1,NCHAN - M + 1 DO 10 K = 1,I1 HNP1(LO+K) = H(LP+K) 10 CONTINUE I1 = NOTERM LO = L*NOTERM + IL LP = L*NRANG2 + IP 20 CONTINUE IF (NCFGP.EQ.0) GOTO 40 C C THE BOUND-CONTINUUM TERMS C DO 30 K = 1,NCFGP HNP1(LO+K) = H(LP+K) 30 CONTINUE 40 CONTINUE 50 CONTINUE C C THE BOUND-BOUND TERMS C IF (NCFGP.EQ.0) GOTO 60 N = NOTERM*NCHAN M = (N* (MN4-N+1))/2 J = M + 1 K = M + (MNP1* (MNP1+1)-NCONHP* (MN2-NCONHP+1))/2 READ (IDISC1) (HNP1(I),I=J,K) C C IF NBUG6=2 WRITE OUT THE HAMILTONIAN MATRIX WHICH IS TO BE C DIAGONALIZED IN THIS LOOP OF RSCT C 60 CONTINUE IF (NBUG6.NE.2) GOTO 100 WRITE (IWRITE,3020) MNP2 M = 0 70 CONTINUE K = M L = M + 1 M = M + 8 N = MIN(M,MNP2) FO8(1) = FL8(1) DO 90 I = 1,N IF (I.LE.K) GOTO 80 FO8(1) = FL8(I-K) L = I 80 CONTINUE MO = ((I-1)* (MN4-I))/2 + L MH = MO + N - L WRITE (IWRITE,FO8) (HNP1(J),J=MO,MH) 90 CONTINUE WRITE (IWRITE,3070) IF (M.LT.MNP2) GOTO 70 100 CONTINUE C C DIAGONALIZE THE HAMILTONIAN MATRIX AND STORE ITS EIGENVALUES IN C THE VALUE ARRAY AND ONE EIGENVECTOR IN THE X ARRAY. FROM WHICH WE C CALCULATE THE SURFACE AMPLITUDES AND THE DIPOLE MATRIX ELEMENTS. C MNP2N IS THE REDUCED NUMBER OF POLES OF ANY PARTITIONED R-MATRIX C MNP2N = (NOTERM*IPRCENT)/100 MNP2N = MNP2N*NCHAN + NCFGP IF(MNP2.NE.MNP2N)WRITE(IWRITE,3025)MNP2,MNP2N C C NLOW = POSITION OF LOWEST EIGENVALUE PRODUCED BY DIAGONALIZATION C NLOW = MNP2N WRITE (IWRITE,3010) IF (NBUG8.EQ.1) NBUG5 = MAX(NBUG5,1) IF (NBUG8.EQ.2) NBUG5 = MNP2N IF (NBUG8.EQ.3) NBUG5 = 0 IF (NBUG5.GT.0) WRITE (IWRITE,*) A ' LOWEST EIGENVALUES (AU) WITH LARGEST EIGENVECTOR COMPONENTS' C C N.B. DIPOLES WORK WITH FULL MNP2 STILL, FOR NOW, XX BEING ZEROED-OUT. C IF (IPOLPH.LT.2) GOTO 130 ILSPI(1,MSKIP) = (100*NSPN+LRGL+1)* (1-2*MOD(NPTY,2)) ILSPI(2,MSKIP) = MNP2 ILSPI(3,MSKIP) = NCHAN ILSPI(4,MSKIP) = -MEM1 MINMAT=MIN(MXMAT,MNP2MX) IF(MNP2.GT.MINMAT)THEN WRITE(IWRITE,*)'TO AID MEMORY ALLOCATION YOU HAVE STATED THAT', X' THE LARGEST VALUE OF MNP2 ON FILE IS',MNP2MX WRITE(IWRITE,*)'HOWEVER A VALUE', X' OF',MNP2,' HAS BEEN ENCOUNTERED. PLEASE RESET MNP2MX.' STOP 'RESET MNP2MX' ENDIF IF(MEM1+MNP2**2+4*MINMAT**2.GT.MXMEM)ILSPI(4,MSKIP) = ABS(MREC1) DO 110 I = 1,MNP2 DO 120 K = 1,MZNR2 XX(K,I) = ZERO 120 CONTINUE 110 CONTINUE 130 CONTINUE C C DIAGONALIZE AND LOOP OVER IPRCENT EIGENVECTORS (NO) C LENGTH = (MNP2* (MNP2+1))/2 !ORIG START NOP=0 IF(MNP2N.NE.MNP2)THEN K=0 TRACE=ZERO DO I=1,MNP2 TRACE=TRACE+HNP1(K+1) K=K+MNP2-I+1 ENDDO ENDIF !ORIG C C DO 210 NO = 1,MNP2N !***ALL*** C C C HSLDR DIAGONALIZATION ROUTINE RETURNS EIGENVECTORS ONE-BY-ONE. C NOP=NOP+1 !ORIG 211 CALL HSLDR(MNP2,HNP1,LENGTH,EPSI,VALUE,X,NOP,AUX,MXMAT) !ORIG C IF(MNP2N.NE.MNP2)THEN IF(NOP.EQ.1)THEN NOP=MNP2-MNP2N+1 GO TO 211 ENDIF IF(NO.EQ.MNP2N)THEN DO N=1,MNP2N VALUE(N)=VALUE(N+MNP2-MNP2N) TRACE=TRACE-VALUE(N) ENDDO TRACE=TRACE/(MNP2-MNP2N) !EZERO WRITE(IWRITE,3045)(TRACE-VALUE(MNP2N))*2 NOP=MNP2N ENDIF ENDIF !ORIG END C C LAPACK 3.0 DRIVERS C CL NOP=NO !LAPACK START CL IF(NO.EQ.1)THEN CL K=0 CL TRACE=ZERO CL DO I=1,MNP2 CL TRACE=TRACE+HNP1(K+1) CL DO J=I,MNP2 CL K=K+1 CL A(J,I)=HNP1(K) !DSYEVR/X CLC B(J,I)=HNP1(K) !DSYEVD CL ENDDO CL ENDDO CLC CL IF(MNP2.EQ.MNP2N)THEN CL RNG='A' CL ELSE !PARTITIONED CL RNG='I' CL IL=1 CL IU=MNP2N CL ENDIF CLC CL IEEEOK=ILAENV(10,'DSYEVR','VL',1,2,3,4) CL IF(IEEEOK.NE.1)THEN CL WRITE(6,*)'WARNING: DSYEVR USING DSTEB/DSTEIN, CALL DSYEVD!' CL STOP 'RSCT: BEST TO CALL DSYEVD!' CL ENDIF CLC CLC DSYEVR (REQUIRES MACHINES TO HANDLE NaNs GRACEFULLY) CLC CL CALL DSYEVR('V',RNG,'L',MNP2,A,MXMAT,VL,VU,IL,IU,-EPSI,ME, CL X X,XA,MXMAT,ISUPP,HNP1,MXTRI,IWORK,LIWORK,INFO) CLC CL IF(ME.NE.MNP2N)THEN CL WRITE(6,*)'RSCT: ERROR IN LAPACK DSYEVR, NOT ALL', CL X ' E-VALUES FOUND:',ME,MNP2N CL STOP 'RSCT: FAILURE IN LAPACK ROUTINE DSYEVR' CL ENDIF CL IF(HNP1(1).GT.MXTRI)THEN CL LWORK=NINT(HNP1(1)) CL WRITE(6,*)'***OPTIMAL USE OF DSYEVR REQUIRES LWORK=',LWORK CL ENDIF CLC CLC DSYEVX, SAFE BUT SLOW. CLC CLC CALL DSYEVX('V',RNG,'L',MNP2,A,MXMAT,VL,VU,IL,IU,-EPSI,ME, CLC X X,XA,MXMAT,HNP1,MXTRI,IWORK,ISUPP,INFO) CLC CLCL IF(ME.NE.MNP2N)THEN CLCL WRITE(6,*)'RSCT: ERROR IN LAPACK DSYEVR, NOT ALL', CLCL X ' E-VALUES FOUND:',ME,MNP2N CLCL STOP 'RSCT: FAILURE IN LAPACK ROUTINE DSYEVR' CLCL ENDIF CLCL IF(HNP1(1).GT.MXTRI)THEN CLCL LWORK=NINT(HNP1(1)) CLCL WRITE(6,*)'***OPTIMAL USE OF DSYEVR REQUIRES LWORK=',LWORK CLCL ENDIF CLC CLC DSYEVD (TO BE USED IF NaNs CAUSE CALCULATION TO ABORT) CLC CLC IF(MNP2N.NE.MNP2)THEN CLC WRITE(6,*)'***ERROR: CANNOT USE PARTITIONED R-MX & DSYEVD' CLC STOP '***ERROR: CANNOT USE PARTITIONED R-MX & DSYEVD' CLC ENDIF CLC CLC CALL DSYEVD('V','L',MNP2,B,MXMAT,X,WORK,LWORK,IWORK,LIWORK CLC X ,INFO) CLC CL IF(INFO.NE.0)THEN CL WRITE(6,*)'RSCT: ERROR IN LAPACK DSYEVR/D: INFO=',INFO CL STOP 'RSCT: FAILURE IN LAPACK ROUTINE DSYEVR/D' CL ENDIF CLC CL DO I=1,MNP2N CL VALUE(I)=X(MNP2N+1-I) CL ENDDO CL IF(MNP2N.NE.MNP2)THEN CL DO I=1,MNP2N CL TRACE=TRACE-VALUE(I) CL ENDDO CL TRACE=TRACE/(MNP2-MNP2N) !EZERO CL WRITE(IWRITE,3045)(TRACE-VALUE(MNP2N))*2 CL ENDIF CL ENDIF CLC CL J=MNP2N+1-NO CL DO I=1,MNP2 CL X(I)=XA(I,J) !DSYEVR/X CLC X(I)=B(I,J) !DSYEVD CL ENDDO !LAPACK END C C DEBUG PRINTOUT OF E-SOLUTIONS C IF (NBUG8.EQ.3) THEN WRITE (IWRITE,3040) VALUE(NOP) WRITE (IWRITE,*) ' EIGENVECTOR' WRITE (IWRITE,3070) (X(I),I=1,MNP2) ENDIF C C PRINT OUT LARGEST EIGENVECTOR COMPONENTS FOR THE NBUG5 C LOWEST EIGENVALUES C IF (ABS(NLOW-NO).LT.NBUG5) THEN NBIG = MIN(4,MNP2) LARGE(1) = 1 TEST = ZERO DO 170 I = 2,MNP2 IF (ABS(X(I)).LT.TEST) GOTO 170 DO 150 M = MIN(I,NBIG) - 1,1,-1 M1 = M + 1 IF (ABS(X(I)).LT.ABS(X(LARGE(M)))) GOTO 160 LARGE(M1) = LARGE(M) 150 CONTINUE M1 = 1 160 CONTINUE LARGE(M1) = I IF (I.GE.NBIG) TEST = ABS(X(LARGE(NBIG))) 170 CONTINUE WRITE (IWRITE,3080) VALUE(NOP),(LARGE(M),X(LARGE(M)),M=1,NBIG) ENDIF C C DETERMINE THE SURFACE AMPLITUDES, WMAT(I,J) C NRB: C TO ALLOW FOR MAXNCB, NEED TO EXTEND WMAT TO COVER NCFGP (NON-TRIVIAL) C M = 0 DO 190 J = 1,NCHAN L = L2P(J) + 1 IE=0 IF(L.LE.LRANG1)IE=MAXNCB(L) SUM = ZERO DO 180 I = 1,NOTERM SUM = SUM + X(I+M)*ENDS(IE+I,L) 180 CONTINUE M = M + NOTERM WMAT(J,NO) = SUM 190 CONTINUE IF (NBUG8.EQ.2) WRITE (IWRITE,3030) (WMAT(K,NO),K=1,NCHAN) IF (IPOLPH.LT.2) GOTO 210 C C IF PHOTOIONIZATION, STORE EIGENVECTORS IN /MEMORY/ OR ON DISK C IF (ILSPI(4,MSKIP).LE.0) THEN DO 200 I = 1,MNP2 ARRAY(MEM1+NO+ (I-1)*MNP2) = X(I) 200 CONTINUE C ELSE KOUNT = MOD(NO-1,NOTERM) + 1 DO 22 II=1,MNP2 XX(KOUNT,II) = X(II) 22 CONTINUE IF (KOUNT.EQ.NOTERM .OR. NO.EQ.MNP2) CALL DA2(2,MREC1,IDISC2, A MZNR2*MNP2,XX) ENDIF C 210 CONTINUE C C WRITE EIGENVALUES AND SURFACE AMPLITUDES TO ITAPE3 C IF (NBUG5.NE.0 .OR. NBUG6.NE.0 .OR. NBUG8.NE.0) THEN ELOW = ENAT(1)*2.0D0 DO 220 I = 1,MNP2N X(I) = VALUE(I)*2.0D0 - ELOW 220 CONTINUE WRITE (IWRITE,3050) ELOW, (X(I),I=1,MNP2N) ENDIF C IF(MNP2N.NE.MNP2)WRITE(ITAPE3)TRACE WRITE (ITAPE3) (VALUE(I),I=1,MNP2N) CALL ORDER(ET,NORDER,NCHAN,-1) WRITE (ITAPE3) ((WMAT(NORDER(K),I),K=1,NCHAN),I=1,MNP2N) IF (IPOLPH.LE.1) RETURN C C FOR PHOTOIONIZATION RUNS, STORE ORDERING ARRAY IN IORDER. C IF STORING EIGENVECTORS IN MEMORY, INCREMENT MEM1 C DO 230 I = 1,NCHAN IORDER(I,MSKIP) = NORDER(I) 230 CONTINUE IF (ILSPI(4,MSKIP).LE.0) MEM1 = MEM1 + MNP2**2 IF (ILSPI(4,MSKIP).GT.0) WRITE (IWRITE, A *) ' VECTORS STORED ON DISK' C 3000 FORMAT (//35X,'SUBROUTINE RSCT'/35X,15 ('-')) 3010 FORMAT (/' DIAGONALIZING THE HAMILTONIAN MATRIX') 3020 FORMAT (/' HAMILTONIAN MATRIX - SIZE =',I9/) 3025 FORMAT (/' HAMILTONIAN MATRIX SIZE REDUCED FROM',I7,' TO',I7/) 3030 FORMAT (' SURFACE AMPLITUDES FOR EACH CHANNEL'/ (8F10.6)) 3040 FORMAT (/' EIGENVALUE/2RY =',F14.6) 3045 FORMAT (/' PARTITIONED R-MATRIX, EZERO=',F10.3,' RYD'/) 3050 FORMAT (/' ENERGY LEVELS WITH RESPECT TO THE GROUND STATE (AT', A F12.5,' RYDBERGS)'/ (1X,10F12.5)) 3060 FORMAT (/' ** WARNING ** NOT1.LT.NRANG2',5X,6 (' -!-')/16X, A 43 ('=')/'HAMILTONIAN MATRIX - SIZE =',I5,I10,' CHANNELS'/) 3070 FORMAT (10F8.4) 3080 FORMAT (F14.6,5 (I5,':',F7.4)) END SUBROUTINE STG3 IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- LOGICAL NBUG10 C COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RECOV/IPLACE COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT C LOGICAL EX C----------------------------------------------------------------------- IREAD = 5 IWRITE = 6 C OPEN (UNIT=IWRITE,FILE='rout3r',STATUS='UNKNOWN', A FORM='FORMATTED') C EX=.TRUE. C INQUIRE (FILE='dstg3',EXIST=EX) IF (.NOT.EX) THEN WRITE (IWRITE,7777) STOP ENDIF 7777 FORMAT (/' dstg3 does not exist....stopping.') C OPEN (UNIT=IREAD,FILE='dstg3',STATUS='OLD',FORM='FORMATTED') C----------------------------------------------------------------------- C C INITIALIZE COUNTER MSKIP: C MSKIP = 1 C C READ THE DATA SPECIFYING THE CASE C 10 CONTINUE C CALL STG3RD C NBUG10 = NBUG5 .NE. 0 .OR. NBUG6 .NE. 0 .OR. NBUG7 .NE. 0 .OR. A NBUG8 .NE. 0 IF (NBUG10 .AND. INAST.NE.0) WRITE (IWRITE,3000) C C READ THE BASIC DATA AND THE HAMILTONIAN MATRIX FROM THE C FILE PRODUCED BY STG2 (NB: ON RETURN NOT2=MIN(NOT2,NRANG2)) C IF (IPLACE.NE.0) GOTO 70 CALL TAPERD IF (IPLACE.NE.0.OR.MSKIP.LT.IABS(IPWINIT)) GOTO 40 CALL WRITOP IF (NBUG7.EQ.3) WRITE (IWRITE,3000) C C---- DIAGONALISE THE HAMILTONIAN (NB: RSCT RETURNS NOT1=MIN(=,NRANG2)) C CALL RSCT IF (NBUG10) WRITE (IWRITE,3000) WRITE (IWRITE,3030) MSKIP C C IF MORE.GT.0 THERE IS MORE DATA, INCREASE MSKIP BY 1, C IF MORE=0 THERE IS NO MORE DATA. C 40 CONTINUE IF (MORE.EQ.0) GOTO 50 IF (NBUG10) WRITE (IWRITE,3010) MSKIP = MSKIP + 1 IF (MSKIP.LE.IPWFINAL) GOTO 10 C C---- PROCESS THE DIPOLE MATRIX ELEMENTS IF IPOLPH EQ 2 C 50 CONTINUE IF (IPOLPH.NE.2 .OR. IPLACE.GT.0 .OR. ITAPE1.EQ.0 X .OR. (IPLACE.LT.0.AND.INAST.NE.0) ) GO TO 60 INAST = MSKIP CALL DMAT 60 CONTINUE WRITE (IWRITE,3020) 70 CONTINUE CLOSE (IDISC1) RETURN C 3000 FORMAT (1X,71 ('*')) 3010 FORMAT (/' REPEAT STG3 WITH NEW DATA') 3020 FORMAT (/35X,'END OF STG3'/35X,11 ('-')) 3030 FORMAT (/' SYMMETRY NUMBER',I4,' COMPLETE.') END SUBROUTINE STG3RD IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C READS THE DATA SPECIFYING THE CASE C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) C CHARACTER TITLE(18)*4,PARITY(0:1)*4,SPIN(0:9)*8,RAD*3 LOGICAL EX,EXA,EXH C DIMENSION NTAPE(5) C COMMON /ADDE/EST(MZTAR),EPOLE,IPOLE,NEST,AST(MZTAR),MERGE(MZTAR) COMMON /BASIC2/CF(MZCHF,MZCHF,MZLMX),ET(MZCHF),KJ(MZCHF), A L2P(MZCHF),LSTARG(MZCHF),NCONAT(MZTAR),LRGL,NSPN,NPTY COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECOV/IPLACE COMMON /REDTOL/TOLER,iparinit COMMON /SHIFT/ESHIFT,NSHIFT COMMON /NRBMNP/MNP2MX,IPRCENT COMMON /NRBSRT/ECORR(MZTAR),NCORR,ISORT COMMON /NRBHDX/EXH C NAMELIST/STG3A/NBUG5,NBUG6,NBUG7,NBUG8,NBUG9,IPOLPH,RAD,MNP2MX X ,ISORT,N2HDAT,IPRCENT NAMELIST/STG3B/NBUT,NCONT,IDIAG,NAST,INAST,MAXC,TOLER,NSHIFT X ,ESHIFT,IPWINIT,IPWFINAL,NCORR,iparinit C DATA NTAPE/1,2,3,4,5/,PARITY/'EVEN',' ODD'/,ZERO/0.0D0/, A SPIN/'(= 2J) ',' SINGLET',' DOUBLET',' TRIPLET',' QUARTET', B ' QUINTET',' SEXTET ',' SEPTET ',' OCTET ',' 2S+1>8 '/ C----------------------------------------------------------------------- C C IF ANY DIMENSION IS EXCEEDED WHEN READING FROM DATA, CALL RECOV2 C WITH IPLACE=0 TO TERMINATE THE PROGRAM C IPLACE = 0 C C IF MSKIP.GT.1, STG3 IS BEING REPEATED FOR NEW (N+1)-ELECTRON DATA C SO IT IS NOT NECESSARY TO READ IN THIS FIRST SET OF DATA AGAIN. C IF (MSKIP.GT.1) GOTO 10 C----------------------------------------------------------------------- MORE = 1 READ (IREAD,3120) TITLE C WRITE (IWRITE,3020) TITLE WRITE (IWRITE,3110) A MZCHF,MZLMX,MZLR2,MZMEG,MZKIL,MZNC2,MZNR2,MZSLP,MZTAR,MXMAT C C ICOPY = 0 ITOTAL = 0 IPUNCH=0 NBUG1=0 NBUG2=0 NBUG3=0 NBUG4=0 C IPOLPH=1 RAD='NO' NBUG5=0 NBUG6=0 NBUG7=0 NBUG8=0 NBUG9=0 C MNP2MX=999999 !MAX MNP2, TO AID MEMORY SET ASIDE ISORT=0 N2HDAT=-1 !DEFAULT READ FROM FILE IPRCENT=100 !NON-PARTITIONED C READ(IREAD,STG3A) C IF(IPRCENT.GT.100)IPRCENT=100 IF(IPRCENT.LE.0)IPRCENT=60 !PARTITIONED IPOLPH = ABS(IPOLPH) IF(RAD.EQ.'YES')IPOLPH=2 C C----------------------------------------------------------------------- C C Set the I/O numbers and open files. C C IREAD (5) .. input data .. dstg3 C IWRITE (6) .. printed output .. rout3r C C IPUNCH .. NOT USED C C IDISC1 (11) .. scratch C IDISC2 (12) .. scratch (DA2) C IDISC3 (13) .. scratch C IDISC4 .. NOT USED C C ITAPE1 (1) .. RECUPD/STG2 dump (dipole) .. RECUPD/STG2D.DAT C if IPOLPH=2 C ITAPE2 (2) .. RECUPD/STG2 dump (hamiltonians) .. C RECUPH/STG2H.DAT/STG2HXXX.DAT C always used C ITAPE3 (3) .. STG3 dump (hamiltonians) .. H.DAT C always used C ITAPE4 (4) .. STG3 dump (dipole matrices) .. D00.DAT etc C if IPOLPH=2 C C----------------------------------------------------------------------- IWRITE = 6 IPUNCH = 0 C IDISC1 = 11 IDISC2 = 12 IDISC3 = 13 IDISC4 = 0 C ITAPE2 = 2 ITAPE3 = 3 IF (IPOLPH.EQ.2) THEN ITAPE1 = 1 ITAPE4 = 4 C FIX FOR DMAT1 TO AVOID DIMENSION ERROR AFTER LARGE CPU TIME USED. C SHOULD REALLY CODE MXDM=MAX(MZNC2,MZNR2,MZCHF) - NRB MZCHF1=MZCHF+1 IF(MZCHF.GT.MXDM)CALL RECOV2('STG3RD','MZNC2 ',MZNC2,MZCHF1) ELSE ITAPE1 = 0 ITAPE4 = 0 ENDIF C OPEN (UNIT=IDISC1,STATUS='SCRATCH',FORM='UNFORMATTED') C IF (ITAPE1.GT.0) THEN INQUIRE (FILE='RECUPD.DAT',EXIST=EX) IF (.NOT.EX) THEN INQUIRE (FILE='STG2D.DAT',EXIST=EX) IF (.NOT.EX) THEN WRITE (IWRITE,7771) STOP 'Neither RECUPD.DAT nor STG2D.DAT exist' ELSE OPEN (UNIT=ITAPE1,FILE='STG2D.DAT',STATUS='OLD', A FORM='UNFORMATTED') ENDIF ELSE OPEN (UNIT=ITAPE1,FILE='RECUPD.DAT',STATUS='OLD', A FORM='UNFORMATTED') ENDIF ENDIF 7771 FORMAT (/' Neither RECUPD.DAT nor STG2D.DAT exist....stopping.') C INQUIRE (FILE='RECUPH.DAT',EXIST=EX) IF (.NOT.EX) THEN INQUIRE (FILE='STG2H.DAT',EXIST=EX) IF (.NOT.EX) THEN INQUIRE (FILE='STG2H000.DAT',EXIST=EX) IF(.NOT.EX) THEN WRITE (IWRITE,7772) STOP 'Neither RECUPH.DAT nor STG2H.DAT/STG2H000.DAT' ELSE OPEN (UNIT=ITAPE2,FILE='STG2H000.DAT',STATUS='OLD', A FORM='UNFORMATTED') WRITE(IWRITE,3210) 3210 FORMAT(' OPENING FILE STG2H000.DAT'/) C INQUIRE (FILE='HDCORR000',EXIST=EXH) IF(EXH)OPEN(31,FILE='HDCORR000',STATUS='OLD', X FORM='FORMATTED') IF(EXH)WRITE(IWRITE,3212) 3212 FORMAT(' OPENING FILE HDCORR000'/) C INQUIRE (FILE='AMATU000.DAT',EXIST=EXA) IF(EXA)OPEN (UNIT=77,FILE='AMATU000.DAT', A STATUS='OLD',FORM='UNFORMATTED') IF(EXA)WRITE(IWRITE,3211) 3211 FORMAT(' OPENING FILE AMATU000.DAT'/) END IF ELSE OPEN (UNIT=ITAPE2,FILE='STG2H.DAT',STATUS='OLD', A FORM='UNFORMATTED') C INQUIRE (FILE='HDCORR.DAT',EXIST=EXH) IF(EXH)OPEN(31,FILE='HDCORR.DAT',STATUS='OLD', X FORM='FORMATTED') C INQUIRE (FILE='AMATU.DAT',EXIST=EXA) IF(EXA)OPEN (UNIT=77,FILE='AMATU.DAT',STATUS='OLD', A FORM='UNFORMATTED') ENDIF ELSE OPEN (UNIT=ITAPE2,FILE='RECUPH.DAT',STATUS='OLD', A FORM='UNFORMATTED') ENDIF 7772 FORMAT (/'Neither RECUPH.DAT nor STG2H.DAT/STG2H000.DAT exist.' A ,' ...stopping.') C C----------------------------------------------------------------------- WRITE (IWRITE,3000) WRITE (IWRITE,3010) A IREAD,IWRITE,IPUNCH, B IDISC1,IDISC2,IDISC3,IDISC4, C ITAPE1,ITAPE2,ITAPE3,ITAPE4 WRITE (IWRITE,3030) NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7, A NBUG8,NBUG9 WRITE (IWRITE,3050) ICOPY,ITOTAL,IPOLPH C IF (ITAPE1.GT.0) WRITE (IWRITE,3070) NTAPE(1),ITAPE1 WRITE (IWRITE,3070) NTAPE(2),ITAPE2 WRITE (IWRITE,3080) NTAPE(3),ITAPE3 IF (ITAPE4.GT.0) WRITE (IWRITE,3080) NTAPE(4),ITAPE4 C IF(IPRCENT.NE.100)WRITE(IWRITE,3085)IPRCENT C WRITE (IWRITE,3040) GOTO 20 C----------------------------------------------------------------------- 10 CONTINUE IF (INAST.EQ.0) GOTO 50 WRITE (IWRITE,3000) IF (INAST.GT.0) GOTO 40 C----------------------------------------------------------------------- C C READ AND WRITE BASIC DATA. C C NOT1 THE NUMBER OF CONTINUUM TERMS TO BE USED (DEFAULT=ALL) C IDIAG .EQ. 1 IF ITAPE2 CONTAINS STG2 HAMILTONIAN MATRICES, C AND THE DIAGONALIZATION DATA IS TO BE STORED ON ITAPE3 C NAST .EQ. 0 IF THE ENERGIES OF THE ATOMIC OR IONIC STATES ARE NOT C TO BE CHANGED *** C .NE. 0 IF THE ENERGIES ARE TO BE CHANGED. C THE VALUE OF NAST IS THEN THE TOTAL NUMBER OF C N-ELECTRON STATES (SEE BELOW FOR INPUT DETAILS). C .GT. 0 THEN IF EST(1).EQ.0, C EST(I)=OBSERVED ENERGIES RELATIVE TO THE GROUND, C IN RYDBERGS (EST(I).GT.0) OR /CM (EST(I).LT.0). C IF EST(1).NE.0, C EST(I)=OBSERVED ABSOLUTE ENERGIES IN A.U. C .LT. 0 MERGE GROUPS OF ENERGIES. C NCORR .GT. 0 READ NCORR TARGET INDEXES AND ENERGY CORRECTIONS C ECORR THAT WILL BE ADDED TO H_DIAG. IF NAST.GT.0 THEY FOLLOW C THOSE OBS. ENERGIES, AND ARE IN THE SAME UNITS, ELSE RYD. C ESHIFT .EQ. SHIFT OF ALL TARGET STATES IN RYDBERGS, CASE OF NAST .GT. 0 C EST(1)=0, TO MOVE THEM RELATIVE TO N+1 CFGS. C NSHIFT .GT.0 THEN READ FROM FILE 'eshift' THE SLPI, THE NUMBER OF C N+1-CFG EIGENVALUES TO BE SHIFTED, AND THEN THE INDEX C REALTIVE TO THE GROUND CONTINUUM OF THE N+1 EIGENVALUE C AND THE SHIFT (IN RY) OF THE EIGENVALUE - TWG C INAST .GE.0 IF STG3 IS TO BE RECYCLED WITH THE SAME BASIC DATA C FOR DIFFERENT (N+1)-ELECTRON SYMMETRIES (FOR ALL IF C .EQ.0, FOR INAST SPECIFIED SYMMETRIES IF .GT.0). C TOLER IS THE TOLERANCE PARAMETER FOR KEEPING/DISCARDING CERTAIN C LINEAR COMBINATIONS OF CONFIGURATIONS. C THE DEFAULT IS -1.0 (OFF) C IF YOU WANT TO KEEP ALL CFGS. SET THIS TO BE LARGER (E.G. 0.15) C IF YOU WANT NO (N+1)-CFGS., SET THIS TO BE 2. OR MORE - TWG & NRB C IPWINIT and C IPWFINAL USED TO SPECIFY AN INITIAL AND FINAL PARTIAL WAVE NUMBER C FROM THE LIST OF PARTIAL WAVES READ FROM STG2. THE CODE C WILL PROCESS ALL PARTIAL WAVES FROM IPWINIT TO IPWFINAL. C DEFAULT: ALL PARTIAL WAVES READ FROM STG2 ARE PROCESSED. C IF IPWINIT IS .GT. 1, THEN THE CODE WILL APPEND THE C NEWLY PROCESSED H.DAT FILE TO THE END OF THE EXISTING C H.DAT FILE. C IPARINIT Specifies the parity of the initial state in c photoionization that does not get reduced in c the pseudoresonance elimination procedure. c This is only necessary for TOLER>0 and IPOLPH=2. c default = 0 C----------------------------------------------------------------------- C 20 CONTINUE C NBUT=1 IDIAG=1 NAST=0 NCONT=999 MAXC=999 INAST=0 C TOLER=0.15 TOLER=-1.0D0 iparinit=0 ESHIFT=0.0D0 NSHIFT=0 IPWINIT=1 IPWFINAL=MZSLP NCORR=0 C READ(IREAD,STG3B) C IF (IPWINIT.EQ.1) THEN OPEN(UNIT=ITAPE3,FILE='H.DAT',STATUS='UNKNOWN', + FORM='UNFORMATTED') ELSE INQUIRE (FILE='H.DAT',EXIST=EX) IF(EX)THEN OPEN(UNIT=ITAPE3,FILE='H.DAT',STATUS='OLD', C + FORM='UNFORMATTED',ACCESS='APPEND') !f77 only + FORM='UNFORMATTED',POSITION='APPEND') !f90 only IPWINIT=-IPWINIT !FLAG SKIP OF HEADER ELSE OPEN(UNIT=ITAPE3,FILE='H.DAT',STATUS='NEW', + FORM='UNFORMATTED') ENDIF INAST=0 ENDIF C ESHIFT=ESHIFT/2.0D0 IF(NSHIFT.GT.0) THEN OPEN(UNIT=81,FILE='ESHIFT.DAT',STATUS='OLD',FORM='FORMATTED') OPEN(UNIT=79,FILE='BBEIGEN.DAT',STATUS='OLD',FORM='FORMATTED') ENDIF C IF(TOLER.GE.1.D-19.AND..NOT.EXA)THEN WRITE(IWRITE,*) X 'FILE AMATU.DAT DOES NOT EXIST BUT TOLER.GT.0' STOP 'FILE AMATU.DAT DOES NOT EXIST BUT TOLER.GT.0' ENDIF C NOT1=MIN(NCONT,MAXC) NOT2=NOT1 C WRITE (IWRITE,3100) NBUT,NOT1,NOT2,IDIAG,NAST,INAST C IF(EXH)write(IWRITE,*)'Using HDCORR to alter H diagonal' C C READ THE OBSERVED ENERGY OF THE ATOMIC OR IONIC TARGET STATES C NEST = NAST IF (NAST.EQ.0) GOTO 30 C C MERGE A GROUP OF ENERGIES TOGETHER C SEE SUBROUTINE TAPERD (GROUPING ENERGIES AROUND MEAN) C IF (NEST.LT.0) THEN READ(IREAD,*) (MERGE(I),I=1,-NEST) WRITE(IWRITE,*) -NEST, * ' target state groups, number of states in each group' WRITE(IWRITE,*) (MERGE(I),I=1,-NEST) GOTO 30 ENDIF C DO N=NAST+1,MZTAR EST(N)=ZERO ENDDO C READ (IREAD,*) (EST(N),N=1,NAST) C IF (EST(1).NE.ZERO) THEN WRITE (IWRITE,*) ' OBSERVED ENERGIES FOR EACH TARGET STATE (AU)' ELSE WRITE (IWRITE,*) ' OBSERVED ENERGIES RELATIVE TO GROUND STATE' IF (EST(NAST).GT.ZERO) WRITE (IWRITE,*) ' (RYD UNITS)' IF (EST(NAST).LT.ZERO) WRITE (IWRITE,*) ' (CM-1 UNITS)' ENDIF C WRITE (IWRITE,3090) (EST(N),N=1,NAST) C 30 CONTINUE C IF(NCORR.GT.0)THEN IF(NAST.EQ.0)THEN DO N=1,MZTAR EST(N)=ZERO ENDDO ENDIF DO N=1,MZTAR ECORR(N)=ZERO ENDDO DO N=1,NCORR READ (IREAD,*)I,ECORR(I) ENDDO ENDIF C----------------------------------------------------------------------- C C READ THE TOTAL ANGULAR MOMENTUM, SPIN, PARITY; C N.B. DO SET INAST=0 WHEN PROCESSING ALL SYMMETRIES! C IF (INAST.EQ.0) GOTO 50 C 40 CONTINUE C READ (IREAD,*) NSPN,LRGL,NPTY C IF (NPTY.LT.0 .OR. NPTY.GT.1) STOP !NSPN.LT.0 .OR. FOR NX LS = MIN(iabs(NSPN),9) WRITE (IWRITE,3060) LRGL,SPIN(LS),PARITY(NPTY) IF (MSKIP.EQ.INAST) THEN MORE = 0 END IF C 50 CONTINUE C 3000 FORMAT (//35X,'SUBROUTINE STG3RD'/35X,17 ('-')) 3010 FORMAT ( A /' INPUT-OUTPUT CHANNEL NUMBERS' B /' ----------------------------' C//' IREAD (',I2,') .. input data .. dstg3' C /' IWRITE (',I2,') .. printed output .. rout3r' C /' IPUNCH (',I2,') .. NOT USED' C /' IDISC1 (',I2,') .. scratch' C /' IDISC2 (',I2,') .. scratch (DA2) ' C /' IDISC3 (',I2,') .. scratch' C /' IDISC4 (',I2,') .. NOT USED' C//' ITAPE1 (',I2,') .. RECUPD/STG2 dump (dipole) ', C '.. RECUPD/STG2D.DAT' C /' if IPOLPH=2' C//' ITAPE2 (',I2,') .. RECUPD/STG2 dump (hamiltonians) ', C '.. RECUPH/STG2H.DAT' C /' always used' C//' ITAPE3 (',I2,') .. STG3 dump (hamiltonians) ', C '.. H.DAT' C /' always used' C//' ITAPE4 (',I2,') .. STG3 dump (dipole matrices) ', C '.. D00 etc' C /' if IPOLPH=2' C ) 3020 FORMAT (//1X,72 ('-')//1X,18A4//1X,72 ('-')////15X, A ' SSSSSSSS TTTTTTTTTT GGGGGGGG 33333333 '/ B 15X, C 'SSSSSSSSSS TTTTTTTTTT GGGGGGGGGG 3333333333'/ D 15X, E 'SS TT GG GG 33 33'/ F 15X, G 'SS TT GG 33'/ H 15X, I 'SS TT GG 33'/ J 15X, K 'SSSSSSSSS TT GG 3333333 '/ L 15X, M ' SSSSSSSSS TT GG GGGG 3333333 '/ N 15X, O ' SS TT GG GGGG 33'/ P 15X, Q ' SS TT GG GG 33'/ R 15X, S ' SS TT GG GG 33 33'/ A 15X, B 'SSSSSSSSSS TT GGGGGGGGGG 33 33'/ C 15X, D ' SSSSSSSS TT GGGGGGGG 33333333 ') 3030 FORMAT (' NBUG PARAMETERS'/9I5) 3040 FORMAT ( A/' BASIC DATA' B/' ----------') 3050 FORMAT ( A/' ICOPY = ',I3 B/' ITOTAL = ',I3 C/' IPOLPH = ',I3) 3060 FORMAT (/' L =',I3,3X,A8,3X,A4/1X,24 ('-')) 3070 FORMAT (38X,' INPUT CHANNEL ITAPE',I1,' =',I5) 3080 FORMAT (38X,'OUTPUT CHANNEL ITAPE',I1,' =',I5) 3085 FORMAT (/'***** PARTITIONED R-MATRIX IN USE: ',I2, X'% OF E-SOLUTIONS IN USE.'/) 3090 FORMAT (5F16.6) 3100 FORMAT (/' NBUT =',I3,' NOT1 =',I3,' NOT2 =',I3,' IDIAG =',I3, A ' NAST =',I3,' INAST =',I3) 3110 FORMAT (//10X,'COMPILED FOR DIMENSIONS -'//15X, A 'MAX. NUMBER OF CHANNELS (NCHAN) MZCHF = ',I6/15X, C 'NO. OF MULTIPOLES IN POTENTIALS (LAMAX) MZLMX = ',I6/15X, E 'NO OF CONTINUUM ANGULAR MOM. (LRANG2) MZLR2 = ',I6/15X, G 'MEGAWORDS OF MEMORY AVAILABLE MZMEG = ',I6/15X, G 'KILOWORDS OF MEMORY AVAILABLE MZKIL = ',I6/15X, I 'N+1 ELECTRON CONFIGURATIONS MZNC2 = ',I6/15X, K 'CONTINUUM ORBITALS (NRANG2) MZNR2 = ',I6/15X, O 'NO OF (N+1) ELECTRON SYMMETRIES MZSLP = ',I6/15X, M 'TARGET STATES (NAST) MZTAR = ',I6/15X, P 'MAX. ORDER OF HAMILTONIAN MATRIX (MNP1) MXMAT = ',I6//) 3120 FORMAT (18A4) END SUBROUTINE TAPERD IMPLICIT REAL*8 (A-H,O-Z) C----------------------------------------------------------------------- C C C C----------------------------------------------------------------------- C C READS THE INPUT DATA FROM STG2 CORRESPONDING TO THE REQUIRED C VALUES OF LRGL, NSPN, NPTY, NCFGP. C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C parameter (nsizex=3) !min 1 PARAMETER (MXN21=MZNR2+1) PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) PARAMETER (MXTRI=(MXMAT*(MXMAT+1))/2) PARAMETER (MXD1=MZNC2/MZNR2,MXD2=MZNR2/MZNC2,MXD3=MXD1+MXD2, A MXDM=MZNC2*MXD1/MXD3+MZNR2*MXD2/MXD3) PARAMETER (MXDMAT=4*MXMAT*MZCHF+2*MXMAT*MXDM+2*MZNR2*MXDM, A MXRSCT=MXTRI+MZCHF*MXMAT+MZNR2*MXMAT) PARAMETER (MXM1=MXDMAT/MXRSCT,MXM2=MXRSCT/MXDMAT,MXM3=MXM1+MXM2, A MXMATX=MXDMAT*MXM1/MXM3+MXRSCT*MXM2/MXM3+1) PARAMETER (MXDUM2=MXMATX-MXTRI-MZNC2*MZNR2) C LOGICAL EX,EXH CHARACTER PARITY(0:1)*4,SPIN(0:9)*8 CHARACTER*6 FL8(8),FO8(2) CHARACTER*1 IFNM(10) CHARACTER*12 IFLNM C common/adip/ad(mznc2,mznc2,nsizex),nlad(mzslp),nrad(mzslp), c ixsym(mzslp),nxsym COMMON /ADDE/EST(MZTAR),EPOLE,IPOLE,NEST,AST(MZTAR),MERGE(MZTAR) COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /BASIC2/CF(MZCHF,MZCHF,MZLMX),ET(MZCHF),KJ(MZCHF), A L2P(MZCHF),LSTARG(MZCHF),NCONAT(MZTAR),LRGL,NSPN,NPTY COMMON /BASIN/EIGENS(MZNR2,MZLR2),ENDS(MXN21,MZLR2),DELTA,ETA COMMON /BUTT/COEFF(3,MZLR2),EK2MAX,EK2MIN,MAXNCB(MZLR2),NELCOR COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /INFORM/IREAD,IWRITE,IPUNCH COMMON /MATRIX/HNP1(MXTRI),STORE(MZNC2,MZNR2),DUM2(MXDUM2) COMMON /MORDER/LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /NBUG/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /RECOV/IPLACE COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST COMMON /SHIFT/ESHIFT,NSHIFT COMMON /NRBSRT/ECORR(MZTAR),NCORR,ISORT COMMON /NRBDIP/MAXLD COMMON /NRBCOD/ICODE COMMON /NRBHDX/EXH C DIMENSION JRELOP(3),NORDER(MZTAR) C SAVE IFL C DATA IFL/1/ DATA PARITY/'EVEN',' ODD'/,TWO/2.0D0/, A ZERO/0.0D0/,SPIN/'(= 2J) ',' SINGLET', A ' DOUBLET',' TRIPLET',' QUARTET',' QUINTET',' SEXTET ', B ' SEPTET ',' OCTET ',' 2S+1>8 '/ DATA FL8/' (8','(14X,7','(28X,6','(42X,5','(56X,4','(70X,3', A '(84X,2',' (98X,'/,FO8(2)/'F14.7)'/ DATA IFNM/'0','1','2','3','4','5','6','7','8','9'/ C----------------------------------------------------------------------- C C IF ANY DIMENSION IS EXCEEDED WHEN READING FROM DATA, CALL RECOV2 C WITH IPLACE=0 TO TERMINATE THE PROGRAM C IPLACE = 0 WRITE (IWRITE,3000) ITAPE = ITAPE2 IF (MSKIP.GT.1 .AND. (INAST.EQ.0.OR. A ABS(IPOLPH).EQ.2)) GOTO 20 REWIND (77) REWIND (ITAPE) C C READ THE BASIC DATA FROM FILE C IF(N2HDAT.GT.0) THEN READ (ITAPE) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,ICODE,LAM,IZESP, A JRELOP,MAXLD ELSE MAXLD=999 READ (ITAPE,ERR=9) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,ICODE,LAM, A IZESP,JRELOP,MAXLD0,N2HDT0 MAXLD=MAXLD0 N2HDAT=N2HDT0 END IF 9 WRITE (IWRITE,3070) NELC,NZ,LRANG1,LRANG2,NRANG2,LAMAX,LAM,JRELOP, A IZESP IF (LRANG1.GT.MZLR2) CALL RECOV2('TAPERD',' MZLR2',MZLR2,LRANG1) IF (LRANG2.GT.MZLR2) CALL RECOV2('TAPERD',' MZLR2',MZLR2,LRANG2) IF (NRANG2.GT.MZNR2) CALL RECOV2('TAPERD',' MZNR2',MZNR2,NRANG2) IF (LAMAX.GT.MZLMX) CALL RECOV2('TAPERD',' MZLMX',MZLMX,LAMAX) NOT2 = MIN(NRANG2,NOT2) READ (ITAPE,ERR=1) (MAXNHF(L),L=1,LRANG1), (MAXNLG(L),L=1,LRANG1), A (MAXNC(L),L=1,LRANG1), (MAXNCB(L),L=1,LRANG1) GO TO 5 CNRB ADD MAXNCB 1 DO L=1,LRANG1 MAXNCB(L)=0 ENDDO 5 WRITE (IWRITE,3050) (MAXNHF(L),L=1,LRANG1) WRITE (IWRITE,3060) (MAXNLG(L),L=1,LRANG1) DO 10 L = 1,LRANG2 NRNG2=NRANG2 IF(L.LE.LRANG1)THEN NRNG2=NRNG2+MAXNCB(L) IF (NRNG2.GT.MZNR2) CALL RECOV2('TAPERD',' MZNR2',MZNR2,NRNG2) ENDIF READ (ITAPE) (EIGENS(N,L),N=1,NRNG2) READ (ITAPE) (ENDS(N,L),N=1,NRNG2) 10 CONTINUE READ (ITAPE) RA,BSTO,HINT,DELTA,ETA,NIX WRITE (IWRITE,3110) RA,BSTO,NIX C uncomment for v3 IF(NIX.GT.0) THEN READ(ITAPE) (IHX,I=1,NIX),(IRX,I=1,NIX) IPTS=2*IRX READ(ITAPE) (POVALU,I=1,IPTS) DO 13 LP=1,LRANG1 NBT=MAXNLG(LP)-LP+1 IF(NBT.GT.0) THEN DO 12 K=1,NBT READ(ITAPE) (POVALU,I=1,IPTS) 12 CONTINUE ENDIF 13 CONTINUE ENDIF READ (ITAPE) ((COEFF(I,L),I=1,3),L=1,LRANG2) READ (ITAPE) NAST WRITE(IWRITE,*)' NAST = ',NAST NAST = ABS(NAST) IF (NAST.GT.MZTAR) CALL RECOV2('TAPERD',' MZTAR',MZTAR,NAST) READ (ITAPE) (ENAT(I),I=1,NAST), (LAT(I),I=1,NAST), A (ISAT(I),I=1,NAST), (LPTY(I),I=1,NAST) C C THE FOLLOWING DATA DEPENDS ON LRGL, NSPN ETC. C 20 CONTINUE READ (ITAPE,END=210) LRGL2,NSPN2,NPTY2,NCFGP2,IPOL2 WRITE (IWRITE,3130) LRGL2,NSPN2,NPTY2,NCFGP2 READ (ITAPE) MNP1,NCONHP,NCHAN IF (NCHAN.GT.MZCHF) CALL RECOV2('TAPERD',' MZCHF',MZCHF,NCHAN) IF (MNP1.GT.MXMAT) CALL RECOV2('TAPERD','MXMAT ',MXMAT,MNP1) IF (NCFGP2.GT.MZNC2) CALL RECOV2('TAPERD',' MZNC2',MZNC2,NCFGP2) READ (ITAPE) (NCONAT(I),I=1,NAST) IF(NSPN2.NE.0)READ (ITAPE) (L2P(I),I=1,NCHAN) IF(NSPN2.EQ.0)READ (ITAPE,ERR=30) (L2P(I),I=1,NCHAN), X (KJ(I),I=1,NCHAN) ! NRB 30 READ (ITAPE) MORE2 C C READ THE UPPER TRIANGLE OF THE HAMILTONIAN MATRIX INTO HNP1 C MN2 = 2*MNP1 DO 110 M = 1,NCHAN MP = (M* (M+1))/2 MO = MP - M + 1 NQ = (M-1)*NRANG2 NBT = 1 DO 100 L = MO,MP C C ---> READ CONTINUUM-CONTINUUM BLOCKS INTO STORE IN TURN IF (ICODE.NE.21) GOTO 50 IF (L.EQ.MP) GOTO 40 IF (L2P(M).EQ.L2P(L-MO+1)) GOTO 50 40 CONTINUE READ (ITAPE) ((STORE(I,J),I=1,NRANG2),J=1,NRANG2) IF (L.GE.MP) GOTO 70 GOTO 60 C 50 CONTINUE READ (ITAPE) ((STORE(J,I),I=1,NRANG2),J=1,NRANG2) 60 CONTINUE KP = (L-MO)*NRANG2 KN = ((MN2-KP+1)*KP)/2 + NQ - KP C C ---> TRANSFER CONTINUUM-CONTINUUM BLOCKS INTO HNP1 FROM STORE 70 CONTINUE C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (MSKIP.LT.IABS(IPWINIT)) GO TO 100 C DO 90 J = 1,NRANG2 IF (L.NE.MP) THEN JJ = (J-1)* (MNP1-KP) - ((J-1)*J)/2 + KN C ELSE NBT = J JJ = ((NQ+J-1)* (MN2-NQ-J+2))/2 - J + 1 ENDIF C DO 80 I = NBT,NRANG2 HNP1(I+JJ) = STORE(I,J) 80 CONTINUE 90 CONTINUE 100 CONTINUE 110 CONTINUE C C READ BOUND-CONTINUUM BLOCKS INTO STORE C IF (NCFGP2.EQ.0) GOTO 160 IF (NCFGP2.GT.MZNC2) CALL RECOV2('TAPERD',' MZNC2',MZNC2,NCFGP2) DO 140 M = 1,NCHAN READ (ITAPE) ((STORE(I,J),I=1,NCFGP2),J=1,NRANG2) C C SKIP IF THIS SYMMETRY IS NOT NEEDED C IF (MSKIP.LT.IABS(IPWINIT)) GO TO 140 C KP = (M-1)*NRANG2 KQ = (KP* (MN2-KP+1))/2 + NCONHP DO 130 J = 1,NRANG2 JJ = (J-1)*MNP1 - (J* (J-1))/2 - J*KP + KQ DO 120 I = 1,NCFGP2 HNP1(I+JJ) = STORE(I,J) 120 CONTINUE 130 CONTINUE 140 CONTINUE C C READ BOUND-BOUND ELEMENTS INTO HNP1 C MP = (NCONHP* (MN2-NCONHP+1))/2 DO 150 I = 1,NCFGP2 MO = MP + 1 MP = MO + NCFGP2 - I READ (ITAPE) (HNP1(J),J=MO,MP) 150 CONTINUE C C READ THE ASYMPTOTIC COEFFICIENTS (ARRAY CF EXPANDED IN WRITOP) C 160 CONTINUE IF (LAMAX.EQ.0) GOTO 180 DO 170 I = 1,NCHAN READ (ITAPE) ((CF(I,J,K),K=1,LAMAX),J=I,NCHAN) 170 CONTINUE 180 CONTINUE c c if HDCORR.DAT exists, subtract corrections from diagonal c IF(EXH.AND.NCFGP2.GT.0)THEN READ(31,*)LRGL,NSPN,NPTY,mdiag IF(LRGL.NE.LRGL2.OR.NSPN.NE.NSPN2.OR.NPTY.NE.NPTY2)THEN WRITE(IWRITE,*) 'MIS-MATCH ON HDCORR.DAT FOR L S P:', X LRGL2,NSPN2,NPTY2,LRGL,NSPN,NPTY STOP 'MIS-MATCH ON HDCORR.DAT FOR L S P' ENDIF IF(MNP1.NE.mdiag)THEN WRITE(IWRITE,*) 'MIS-MATCH ON HDCORR.DAT FOR MNP1:', X MNP1,mdiag STOP 'MIS-MATCH ON HDCORR.DAT FOR MNP1' ENDIF do j=1,MNP1 iii = 1 + ((j-1)*(2*MNP1-j+2))/2 read(31,*)jjj,hcorr c write(IWRITE,*)j,HNP1(iii),HNP1(iii)-hcorr,hcorr HNP1(iii) = HNP1(iii) - hcorr enddo ENDIF C C CHECK THAT LRGL,NSPN,NPTY,NCFGP HAVE THE REQUIRED VALUES. C IF NOT, READ IN MORE DATA FROM FILE OR USER. C IF (MORE2.EQ.0) THEN MORE = 0 END IF IF (INAST.NE.0) GOTO 190 LRGL = LRGL2 NSPN = NSPN2 NPTY = NPTY2 190 CONTINUE IF (LRGL.EQ.LRGL2 .AND. NSPN.EQ.NSPN2 .AND. A NPTY.EQ.NPTY2) GOTO 230 IF (IPOLPH.LE.1) GOTO 200 WRITE (IWRITE,3140) MORE = 0 GOTO 220 C 200 CONTINUE IF (MORE2.GE.1) GOTO 20 210 CONTINUE C IF(iabs(N2HDAT).EQ.IFL) THEN !NO MORE FILES IF(INAST.EQ.0) THEN MORE=0 ELSE WRITE (IWRITE,3020) !! parallel IF(MSKIP.LT.INAST) MORE=1 END IF ELSE C C OPEN NEW FILE STG2HXXX.DAT IFL=IFL+1 IFLNM='STG2H00'//IFNM(IFL)//'.DAT' INQUIRE (FILE=IFLNM,EXIST=EX) IF(.NOT.EX) THEN WRITE(IWRITE,3200) IFLNM 3200 FORMAT('THE FILE ',A12,' DOES NOT EXIST .... stopping') STOP 'STG2HXXX.DAT DOES NOT EXIST!' ELSE CLOSE(ITAPE) OPEN (UNIT=ITAPE,FILE=IFLNM,STATUS='OLD',FORM='UNFORMATTED') WRITE(IWRITE,3210) IFLNM 3210 FORMAT(' OPENING FILE ',A12/) REWIND (ITAPE) END IF IFLNM='HDCORR00'//IFNM(IFL) INQUIRE (FILE=IFLNM,EXIST=EXH) IF(EXH)THEN CLOSE(31) OPEN (31,FILE=IFLNM,STATUS='OLD',FORM='FORMATTED') WRITE(IWRITE,3210) IFLNM ENDIF GO TO 20 END IF 220 CONTINUE IPLACE = -1 GOTO 390 230 CONTINUE C C SKIP IF THIS SYMMETRY IS NOT NEEDED C IF (MSKIP.LT.IABS(IPWINIT)) GO TO 235 C C TRANSFORM THE HAMILTONIAN TO REDUCE NUMBER OF N+1-CFGS - TWG 13/4/95 C if(ipolph.eq.2)ixsym(mskip)=0 IF(NCFGP2.GT.0)THEN CALL HAMNEW(HNP1,MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) MN2=2*MNP1 ENDIF NCFGP = NCFGP2 C C SHIFT THE N+1 CFGS - TWG 18/11/95. C IF(NCFGP2.GT.0)THEN CALL HAMNEW2(HNP1,MNP1,NCFGP2,LRGL2,NSPN2,NPTY2) ENDIF NCFGP = NCFGP2 C C WRITE OUT THE TARGET STATES. C 235 IF (MSKIP.GT.1) GOTO 250 WRITE (IWRITE,3030) CALL ORDER(ENAT,NORDER,NAST,1) IF (NCORR+ABS(NEST).EQ.0) THEN GROUND = ENAT(NORDER(1)) DO 240 II = 1,NAST I = NORDER(II) T = (ENAT(I)-GROUND)*2 IS = MIN(ISAT(I),9) WRITE (IWRITE,3040) SPIN(IS),LAT(I),PARITY(LPTY(I)),NCONAT(I), A ENAT(I),ZERO,T,II 240 CONTINUE ENDIF C C IF NEST=NAST, ADD ON THE EXTRA ENERGY TO EACH STATE AND TO THE C ASSOCIATED DIAGONAL ELEMENTS OF THE HAMILTONIAN MATRIX. C C EXTENSION '90JUL1/5: OBS.ENERGY INPUT IN RYD ABOVE 1 OR IN CM**-1 C EST(I)=0 ENSURES UNCORRECTED TERM ENERGY IF NO OBS.DATA AVAILABLE C 250 CONTINUE IF (NCORR+ABS(NEST).EQ.0) GOTO 310 MP = 0 E0 = ENAT(1) IF (MSKIP.LE.1) THEN C C IF NEST .LT. 0 GROUP ENERGIES AROUND ARITHMETIC MEANS C IF (NEST.LT.0) THEN DO II = 1,NAST EST(II) = ENAT(II) ENDDO ISTART = 0 DO 258 I = 1,-NEST IF(ISTART.GE.NAST.OR.MERGE(I).LE.0) GOTO 258 IF(MERGE(I).EQ.1) GOTO 256 EMEAN = 0.0D0 IEND = MIN(ISTART+MERGE(I),NAST) DO II=ISTART+1,IEND EMEAN = EMEAN+ENAT(II) ENDDO EMEAN = EMEAN / (IEND-ISTART) DO II=ISTART+1,IEND EST(II) = EMEAN ENDDO 256 ISTART = ISTART + MERGE(I) 258 CONTINUE C IF (ISTART.NE.NAST) THEN WRITE(IWRITE,*) '**** WARNING **** ', A 'Target State Energy MERGE only treats', B ISTART,' states' ENDIF ENDIF C C IF NCORR,ABS(NEST) .GT. 0 CORRECT ENERGIES C IF(EST(1).EQ.ZERO)E0=E0+ESHIFT DO 270 II = 1,NAST I = NORDER(II) I0=I IF(ISORT.GT.0)I0=II C AST(I) = ZERO IF(NEST.NE.0)AST(I) = EST(I0) - ENAT(I) IF(NCORR.GT.0)AST(I)=AST(I)+ECORR(I0) C IF (EST(1).EQ.ZERO) THEN AST(I) = ZERO IF(NEST.NE.0)AST(I) = EST(I0) IF(NCORR.GT.0)THEN IF(ECORR(I0).NE.0)THEN AST(I)=AST(I)+ECORR(I0) IF(EST(I0).EQ.0)AST(I)=AST(I)+(ENAT(I)-E0)*TWO ENDIF ENDIF AST(I)=AST(I)/TWO C IF (AST(I).EQ.ZERO)THEN ENAT(I) = ENAT(I) + ESHIFT GOTO 260 ENDIF C IF (EST(I0).LT.ZERO) AST(I) = -AST(I)/109737.3D0 AST(I) = AST(I) - ENAT(I) + E0 ENDIF C ENAT(I) = ENAT(I) + AST(I) C 260 CONTINUE T = (ENAT(I)-ENAT(1))*TWO IS = MIN(ISAT(I),9) WRITE(IWRITE,3040)SPIN(IS),LAT(I),PARITY(LPTY(I)),NCONAT(I), A ENAT(I),AST(I),T,II 270 CONTINUE ELSE IF(INAST.GT.0)THEN C ENAT HAS BEEN RE-READ AND SO MUST BE RE-ADJUSTED : NRB 4/11/94 CALL ORDER(ENAT,NORDER,NAST,1) DO 275 II=1,NAST I=NORDER(II) ENAT(I)=ENAT(I)+AST(I) 275 CONTINUE ENDIF C END MODS ENDIF C C SKIP IF THIS SYMMETRY IS NOT NEEDED IF (MSKIP.LT.IABS(IPWINIT)) GO TO 390 C C C LOOP OVER THE CHANNELS ASSOCIATED WITH EACH STATE IN THE H-MATRIX C DO 300 I = 1,NAST IF (NCONAT(I).LE.0) GOTO 300 MO = MP + 1 MP = MP + NCONAT(I) DO 290 M = MO,MP NQ = (M-1)*NRANG2 DO 280 J = 1,NRANG2 JJ = ((NQ+J-1)* (MN2-NQ-J+2))/2 + 1 HNP1(JJ) = HNP1(JJ) + AST(I) 280 CONTINUE 290 CONTINUE 300 CONTINUE C C WRITE THE FULL HAMILTONIAN MATRIX TO SCRATCH DISC C 310 CONTINUE IF (NOT1.GE.NRANG2) GOTO 330 REWIND IDISC1 DO 320 M = 1,NCHAN MO = ((M-1)*NRANG2* (MN2- (M-1)*NRANG2+1))/2 + 1 MP = (M*NRANG2* (MN2-M*NRANG2+1))/2 WRITE (IDISC1) (HNP1(I),I=MO,MP) 320 CONTINUE MO = (NCONHP* (MN2-NCONHP+1))/2 + 1 MP = (MNP1* (MNP1+1))/2 WRITE (IDISC1) (HNP1(I),I=MO,MP) C C WRITE OUT THE HAMILTONIAN MATRIX IF NBUG6=1 C 330 CONTINUE WRITE (IWRITE,3090) MNP1,NCHAN IF (NBUG6.NE.1) GOTO 370 M = 0 340 CONTINUE K = M L = M + 1 M = M + 8 JJ = MIN(M,MNP1) FO8(1) = FL8(1) DO 360 I = 1,JJ IF (I.LE.K) GOTO 350 FO8(1) = FL8(I-K) L = I 350 CONTINUE MO = ((I-1)* (MN2-I))/2 + L MP = MO + JJ - L WRITE (IWRITE,FO8) (HNP1(J),J=MO,MP) 360 CONTINUE WRITE (IWRITE,3040) IF (M.LT.MNP1) GOTO 340 C 370 CONTINUE WRITE (IWRITE,3100) IF (ABS(IPOLPH).EQ.2) GOTO 380 IF (MORE.LT.2 .AND. INAST.NE.0)THEN REWIND (ITAPE) REWIND (77) ENDIF GOTO 390 C 380 CONTINUE IF (MSKIP.GT.MZSLP) CALL RECOV2('TAPERD',' MZSLP',MZSLP,MSKIP) C 390 CONTINUE IF (MSKIP.LT.IABS(IPWINIT)) WRITE (IWRITE,3105)MSKIP,IABS(IPWINIT) C 3000 FORMAT (/35X,'SUBROUTINE TAPERD'/35X,17 ('-')) 3020 FORMAT (' CANNOT FIND ON INPUT FILE ANY DATA FOR THIS CASE'/) 3030 FORMAT (' TARGET OR CORE STATE DATA'/8X,' L', A 'PARITY CHANNELS ENERGY(AU) CORRECTION ENERGY(RY)') 3040 FORMAT (A8,I3,1X,A4,I7,F16.7,F12.7,F12.5,I8) 3050 FORMAT (' MAXNHF =',20I3) 3060 FORMAT (' MAXNLG =',20I3) 3070 FORMAT (' NELC =',I3,' NZ =',I3,' LRANG1 =',I3,' LRANG2 =',I3, A ' NRANG2 =',I3,' LAMAX =',I2,' LAM=', B I1/' MASS-CORRECTION(',I1,'), DARWIN-TERM(',I1, C '), SPIN-ORBIT(',I2,');',7X,'IZESP =',I3) 3090 FORMAT (/' HAMILTONIAN MATRIX - SIZE =',I5,I10,' CHANNELS'/) 3100 FORMAT (/' CORRECT DATA READ FROM STG2') 3105 FORMAT (/' (CALCULATION SKIPPED - SYMMETRY=: ',I3,' .LT. ',I3) 3110 FORMAT (' RA =',F10.5,', BSTO =',E12.5,'; NIX =',I2/) 3130 FORMAT (' LRGL =',I3,' NSPN =',I3,' NPTY =',I3,' NCFGP =',I5) 3140 FORMAT ( A ' SPECIFIED CASE OUT OF STEP WITH THIS CASE WHILE IPOLPH.GT.1' B /' REMAINING CASES SKIPPED') END SUBROUTINE VMUL(A,B,C,L,M,N) IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C MATRIX MULTIPLICATION A*B = C, C OUTPUTS RESULTS FOR C(1:L,1:N) C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C PARAMETER (MXMAT=MZCHF*MZNR2+MZNC2) C DIMENSION A(MZNR2,MXMAT),B(MXMAT,N),C(MZNR2,N) C----------------------------------------------------------------------- DO 40 J = 1,N DO 10 I = 1,L C(I,J) = 0.0D0 10 CONTINUE DO 30 K = 1,M DO 20 I = 1,L C(I,J) = A(I,K)*B(K,J) + C(I,J) 20 CONTINUE 30 CONTINUE 40 CONTINUE C END SUBROUTINE WRITOP IMPLICIT REAL*8 (A-H,O-Z) C C C C----------------------------------------------------------------------- C C WRITES THE BASIC QUANTITIES TO THE STG3 FILE C C----------------------------------------------------------------------- C INCLUDE 'PARAM' C PARAMETER (MXN21=MZNR2+1) C DIMENSION NORDER(MZTAR),IORDER(MZCHF) C COMMON /BASIC1/BSTO,RA,NELC,NZ,LRANG2,NRANG2,MAXNHF(MZLR2), A MAXNLG(MZLR2),MAXNC(MZLR2),LRANG1 COMMON /BASIC2/CF(MZCHF,MZCHF,MZLMX),ET(MZCHF),KJ(MZCHF), A L2P(MZCHF),LSTARG(MZCHF),NCONAT(MZTAR),LRGL,NSPN,NPTY COMMON /BASIN/EIGENS(MZNR2,MZLR2),ENDS(MXN21,MZLR2),DELTA,ETA COMMON /BUTT/COEFF(3,MZLR2),EK2MAX,EK2MIN,MAXNCB(MZLR2),NELCOR COMMON /CASES/MORE,MSKIP,IPWINIT,IPWFINAL,IPOLPH,INAST COMMON /CUPINT/MNP1,NCONHP,NCHAN,N2HDAT COMMON /DISTAP/IDISC1,IDISC2,IDISC3,IDISC4,ITAPE1,ITAPE2,ITAPE3, A ITAPE4 COMMON /MORDER/LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /RCASES/NOT1,NOT2,NOT3,NPROG,NCFGP COMMON /STATE/ENAT(MZTAR),LAT(MZTAR),ISAT(MZTAR),LPTY(MZTAR),NAST COMMON /NRBCOD/ICODE C COMMON /NRBDIP/MAXLD COMMON /NRBMNP/MNP2MX,IPRCENT C----------------------------------------------------------------------- C C IF MSKIP.GT.1, STG3 IS BEING REPEATED WITH NEW LRGL, NSPN, ETC. C SO IT IS UNNECESSARY TO WRITE OUT THIS FIRST SET OF DATA AGAIN. C CALL ORDER(ENAT,NORDER,NAST,1) C ITAPE = ITAPE3 NOTERM=MIN(NOT1,NRANG2) C IF (MSKIP.GT.IPWINIT) GOTO 20 C REWIND ITAPE LRNG2=LRANG2 !NRB IF(ICODE.EQ.25)LRNG2=-LRANG2 !FLAG FROM DARC NASTY=NAST IF(IPRCENT.NE.100)NASTY=-NAST !FLAG PARTITIONED WRITE (ITAPE) NELC,NZ,LRNG2,LAMAX,NASTY,RA,BSTO C ,MAXLD WRITE (ITAPE) (ENAT(NORDER(I)),I=1,NAST) WRITE (ITAPE) (LAT(NORDER(I)),I=1,NAST) WRITE (ITAPE) (ISAT(NORDER(I)),I=1,NAST) WRITE (ITAPE) ((COEFF(I,L),I=1,3),L=1,LRANG2) C IF(IPRCENT.NE.100)THEN !FOR PARTITIONED WRITE(ITAPE)IPRCENT,NOTERM DO L = 1,LRANG2 WRITE (ITAPE) (EIGENS(N,L),N=1,NOTERM) WRITE (ITAPE) (ENDS(N,L),N=1,NOTERM) ENDDO ENDIF C 20 CONTINUE C C THE FOLLOWING DATA DEPENDS ON LRGL, NSPN, NPTY C NF = 0 DO 40 i0 = 1,NAST i=norder(i0) IF (NCONAT(I).LE.0) GOTO 40 k0=0 do j=1,i-1 k0=k0+nconat(j) enddo NS = NF + 1 NF = NF + NCONAT(I) DO 30 K = NS,NF k0=k0+1 iorder(k)=k0 !i.e. replaces call order(et,... NRB et(k0)=-k !but still need for rsct, so... 30 CONTINUE 40 CONTINUE C C MNP2 = NOTERM*NCHAN + NCFGP MNP2N = (NOTERM*IPRCENT)/100 MNP2N = MNP2N*NCHAN + NCFGP WRITE (ITAPE) LRGL,NSPN,NPTY,NCHAN,MNP2N,MORE C WRITE (ITAPE) (NCONAT(NORDER(I)),I=1,NAST) IF(NSPN.NE.0)WRITE (ITAPE) (L2P(IORDER(I)),I=1,NCHAN) IF(NSPN.EQ.0)WRITE (ITAPE) (L2P(IORDER(I)),I=1,NCHAN) X ,(KJ(IORDER(I)),I=1,NCHAN) ! NRB C IF (LAMAX.EQ.0) RETURN DO 70 K = 1,LAMAX DO 60 I = 1,NCHAN DO 50 J = 1,NCHAN CF(J,I,K) = CF(I,J,K) 50 CONTINUE 60 CONTINUE 70 CONTINUE WRITE (ITAPE) (((CF(IORDER(I),IORDER(J),K),I=1,NCHAN),J=1,NCHAN), A K=1,LAMAX) C END