C PROGRAM TO COMPUTE TOTAL, DEGENERACY AVERAGED CROSS SECTIONS, C SIG(J2,J1), FROM S-MATRICES SAVED ON TAPE BY SCATTERING PROG. C 14 OCT 94 (SG) MOLSCAT VERSION 14 COMPATIBLE C C NAMELIST (SAVE) VARIABLES: C JTOTL - SKIP JTOT.LT.JTOTL C JTOTU - SKIP JTOT.GT.JTOTU C JSTEP - SKIP JTOT NOT EQUAL TO JTOTL(JSTEP)JTOTU C JTOUT - IF .LT. 999999 (DEFAULT) NO ECHO FOR JTOT.LT.JTOUT C IT - INPUT TAPE UNIT (DEFAULT 10) C MXSIG - IF .GT. 0 WILL LIMIT THE NUMBER OF ACCUM SIG C LFMT - FORMATTED/UNFORMATTED DATA SET C C *** ***** ***** ***** ***** ***** ***** ***** ***** ***** ***** *** C *** ***** LFMT CONTROLS FORMATTED (VERSION 10 AND EARLIER) * *** C *** ***** UNFORMATTED (VERSION 11 AND LATER) * *** C *** ***** ***** ***** ***** ***** ***** ***** ***** ***** ***** *** C C *** NOTE. WORKS FOR ITYPE=1,2,5 (AND EP AND CS EQUIVALENTS) C *** CHANGES 4/20/92 FOR ITYPE=3 C C SIG(JF,JI) OUTPUT ON CARDS IN FORMAT FOR STATISTICAL EQUILIBRIUM C PROGRAM. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MXCH=300,MXJLEV=500,MXS=MXCH*MXCH,MXSG=50000) C DIMS OF ELEVEL (IN COMMON/CMBASE/) AND ENERGY FIXED BY MOLSCAT PARAMETER (MXEL=1000,MXNRG=100) DIMENSION NBASIS(MXCH),J(MXCH),L(MXCH),WVEC(MXCH),JLEV(MXJLEV) DIMENSION ENERGY(MXNRG),MINJT(MXNRG),MAXJT(MXNRG),ISTSIG(MXNRG) DIMENSION SREAL(MXS),SIMAG(MXS),SIG(MXSG) DIMENSION ELEVEL(MXEL) INTEGER ITYPE,NLEV,NQN,NNRG,NOPEN CHARACTER*4 LABEL(20) CHARACTER*1 COUT,CC,CS,BL,ST,XST LOGICAL LOGI,LOGJ,LFMT,LFORM,LPRINT C C LFORM AND LPRINT CONTROL CROSS SECTION OUTPUT DATA LFORM/.TRUE./, LPRINT/.FALSE./ C NAMELIST /SAVE/ JTOTL,JTOTU,JSTEP,JTOUT,IT,MXSIG,LFMT C C DEFAULT VALUES FOR NAMELIST VARIABLES DATA IT/10/, MXSIG/ 0/ DATA JTOTU/999999/,JTOTL/0/,JSTEP/1/,JTOUT/999999/ DATA LFMT/.FALSE./ C C CRITERION FOR PRINTING SIGMA DATA EPS/1.D-8/ DATA JSTEPT/-1/ DATA COUT/'X'/, CC/'C'/, CS/'S'/, BL/' '/,ST/'*'/ DATA PI/3.14159 26535 89793 D0/ C C FORMATS FOR SAVE TAPE (FOR VERSION 11 AND EARLIER) 100 FORMAT(20A4/3I4,F8.4,I4) 101 FORMAT(20I4) 102 FORMAT(I4/(5E16.8)) 103 FORMAT(2I4,E16.8,I4,E16.8,I4) 104 FORMAT(I4/(2I4,E16.8) ) 105 FORMAT(5E16.8) C C SUPPRESS UNDERFLOWS C CALL XUFLOW(0) C C INITIALIZATION . . . DO 10 I=1,MXNRG MINJT(I)=-1 10 MAXJT(I)=-1 DO 11 I=1,MXEL 11 ENERGY(I)=0.D0 C C GET (POSSIBLE) NAMELIST INPUT DATA READ(5,SAVE,END=8000) 8000 WRITE(6,608) IT,JTOTL,JTOTU,JSTEP,JTOUT 608 FORMAT(' *** PROGRAM TO COMPUTE STATE-TO-STATE CROSS SECTIONS'/ 1 ' FROM S-MATRICES SAVED ON TAPE UNIT',I3/ 2 ' NAMELIST PARAMETERS ARE JTOTL =',I8/ 3 30X,'JTOTU =',I8/30X,'JSTEP =',I8/30X,'JTOUT =',I8) IF (MXSIG.GT.0) WRITE(6,643) MXSIG 643 FORMAT(' *** NOTE POSSIBLE LIMIT BY MXSIG =',I6) IF (JTOUT.LT.999999) THEN WRITE(6,688) JTOUT 688 FORMAT(' *** INPUT ECHO SUPPRESSED FOR JTOT.LE.',I7) ENDIF C C GET AND PROCESS HEADER RECORD FROM SAVE TAPE IF (LFMT) THEN READ(IT,100) LABEL,ITYPE,NLEV,NQN,URED,IP ELSE READ(IT) LABEL,ITYPE,NLEV,NQN,URED,IP ENDIF WRITE(6,600) LABEL,ITYPE,NLEV,NQN,URED,IP 600 FORMAT(/' SAVE TAPE LABEL =',20A4/ 2 ' ITYPE =',I3,5X,'NLEV, NQN =',2I4,5X,'URED =',F8.4, 3 5X,'IPROGM =',I3) IF (LFMT.AND.IP.GE.11) THEN WRITE(6,*) ' *** IPROGM INCONSISTENT WITH FORMATTED TAPE' STOP ENDIF COUT=BL IF (ITYPE.LT.10) COUT=CC IF (ITYPE.GE.20 .AND. ITYPE.LT.30) COUT=CS ITYPE=ITYPE-10*(ITYPE/10) IF (.NOT.(ITYPE.EQ.1.OR.ITYPE.EQ.2.OR.ITYPE.EQ.5.OR.ITYPE.EQ.6 1 .OR.ITYPE.EQ.7.OR.ITYPE.EQ.3)) THEN WRITE(6,*) ' *** CODE NOT IMPLEMENTED FOR ITYPE',ITYPE STOP ENDIF C C GET JLEV INFORMATION FROM TAPE NSQ=NLEV*NQN IF (NSQ.GT.MXJLEV) THEN WRITE(6,*) ' *** NLEV*NQN.GT.MXJLEV',MXJLEV STOP ENDIF IF (LFMT) THEN READ(IT,101) (JLEV(I),I=1,NSQ) ELSE READ(IT) (JLEV(I),I=1,NSQ) ENDIF IDENT=0 IF (ITYPE.EQ.3) CALL IDCHK(IDENT,NLEV,JLEV) C C GET NLEVEL/ELEVEL INFORMATION IF (LFMT) THEN IF (IP.GE.3) READ(IT,102) NLEVEL,(ELEVEL(I),I=1,NLEVEL) ELSE READ(IT) NLEVEL,(ELEVEL(I),I=1,NLEVEL) ENDIF DO 1002 I=1,NLEV 1002 WRITE(6,601) I,ELEVEL(ABS(JLEV(NLEV*(NQN-1)+I))), 1 (JLEV(NLEV*(K-1)+I),K=1,NQN) 601 FORMAT(' LEVEL',I4,' ENERGY',F10.3,' QUANTUM NOS. ',10I4) C C GET COLLISION ENERGIES FROM SAVE TAPE IF (LFMT) THEN READ(IT,102) NNRG,(ENERGY(I),I=1,NNRG) ELSE READ(IT) NNRG,(ENERGY(I),I=1,NNRG) ENDIF IF (NNRG.GT.MXNRG) THEN WRITE(6,*) ' *** NNRG.GT.MXNRG',NNRG,MXNRG STOP ENDIF WRITE(6,602) NNRG,(ENERGY(I),I=1,NNRG) 602 FORMAT(/' NUMBER OF COLLISION ENERGIES IS',I4/(' ',1P,5D14.4)) C C**** GET TOP 'LEVEL NO.' FROM JLEV(I,NQN) ILOW=NLEV*(NQN-1)+1 ITOP=NLEV*NQN NTOP=0 DO 1101 I=ILOW,ITOP 1101 NTOP=MAX(NTOP,ABS(JLEV(I))) WRITE(6,681) NTOP,NLEV 681 FORMAT(/' *** NO. OF LEVELS FOR SIGMA(,) IS',I4,' NLEV =',I4) IF (MXSIG.GT.0) NTOP=MIN(NTOP,MXSIG) NSQ=NTOP*NTOP ISTSIG(1)=0 DO 1001 I=2,NNRG 1001 ISTSIG(I)=ISTSIG(I-1)+NSQ C CLEAR SIG() ISTOP=NTOP*NTOP*NNRG IF (ISTOP.GT.MXSG) THEN WRITE(6,*) ' *** STORAGE FOR SIGMA() EXCEEDED',ISTOP,MXSG STOP ENDIF DO 1100 I=1,ISTOP 1100 SIG(I)=0. C C START READ SAVE TAPE LOOP OVER JTOT / ENERGY VALUES . . . C JSTEPX=-1 JTOLD=-1 C C FIRST, GET JTOT,INRG,M RECORD ----- 2000 IF (LFMT) THEN READ(IT,103,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M ELSE IF (IP.LT.14) THEN READ(IT,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M ELSE READ(IT,END=9000) JTOT,INRG,ECHK,IEXCH,WT,M,NOPEN ENDIF ENDIF C N.B. JTOT .LT. 0 INDICATES A OBSOLETE 'CHECKPOINT' RECORD. IF (JTOT.LT.0) THEN WRITE(6,698) JTOT 698 FORMAT('0 * * * NOTE. ROLLOUT MARKER ENCOUNTERED. NROLL =',I6) GO TO 2000 ENDIF C ECHO INPUT IF (JTOUT.GE.999999.OR.(JTOUT.LT.999999.AND.JTOT.GT.JTOUT)) THEN IF (IP.LT.14) THEN WRITE(6,662) IT,JTOT,M,INRG,ECHK,IEXCH,WT 662 FORMAT(' UNIT(',I2,'): JTOT ',I4,'.',I3,' E(',I3,')=', 1 F10.3,' IEX,WT=',I2,F6.2) ELSE WRITE(6,663) IT,JTOT,M,INRG,ECHK,IEXCH,WT,NOPEN 663 FORMAT(' UNIT(',I2,'): JTOT ',I4,'.',I3,' E(',I3,')=', 1 F10.3,' IEX,WT,NOP=',I2,F7.2,I5) ENDIF ENDIF C CALCULATE JSTEP VALUE FROM TAPE (JSTEPT) IF (JTOT.EQ.JTOLD .OR. JTOLD.EQ.-1) GO TO 2011 JSTEPT=JTOT-JTOLD IF (JSTEPT.EQ.JSTEPX) GO TO 2011 IF (JSTEPX.EQ.-1) GO TO 2912 WRITE(6,649) JSTEPT,JSTEPX 649 FORMAT('0 *** ERROR JSTEPT,JSTEPX =',2I6) STOP 2912 JSTEPX=JSTEPT 2011 JTOLD=JTOT C CHECK THAT ENERGY(INRG) CORRESPONDS WITH HEADER RECORD. . . IF (ABS(ECHK-ENERGY(INRG)).LT.1.D-6) GO TO 2002 WRITE(6,697) INRG,ECHK 697 FORMAT('0 * * * WARNING. FOR ',I4,'-TH ENERGY, ECHECK =',D16.8) C C GET J,L,WVEC VALUES ----- 2002 IF (LFMT) THEN READ(IT,104,END=9090) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN) IF (NOPEN.GT.MXCH) THEN WRITE(6,696) MXCH,JTOT,INRG 696 FORMAT('0 *** ERROR. NO. OF OPEN CHANNELS EXCEEDS',I6, 1 ' JTOT,INRG =',2I6) STOP ENDIF ELSE IF (IP.GE.14) THEN IF (NOPEN.GT.MXCH) THEN WRITE(6,696) MXCH,JTOT,INRG STOP ENDIF READ(IT,END=9090) (J(I),L(I),WVEC(I),I=1,NOPEN) ELSE READ(IT,END=9090) NOPEN,(J(I),L(I),WVEC(I),I=1,NOPEN) IF (NOPEN.GT.MXCH) THEN WRITE(6,696) MXCH,JTOT,INRG STOP ENDIF ENDIF ENDIF C C FINALLY GET S-MATRICES --- NSQ=NOPEN*NOPEN IF (LFMT) THEN READ(IT,105,END=9090) (SREAL(I),I=1,NSQ) READ(IT,105,END=9090) (SIMAG(I),I=1,NSQ) ELSE CALL SREAD(IT,NOPEN,SREAL,IEND) IF (IEND.GT.0) GO TO 9090 CALL SREAD(IT,NOPEN,SIMAG,IEND) IF (IEND.GT.0) GO TO 9090 ENDIF C C ALLOW FOR SKIPPING S-MATRICES ON JTOTL,JTOTU,JSTEP INPUT. JDIF=JTOT-JTOTL IF (JTOT.LT.JTOTL .OR. JTOT.GT.JTOTU .OR. 1 JDIF-JSTEP*(JDIF/JSTEP).NE.0) THEN WRITE(6,682) JTOT 682 FORMAT(' --- SKIPPING JTOT',I6,' NOT IN JTOTL (JSTEP) JTOTU') GO TO 2000 ENDIF C C ACCUMULATE TO CROSS SECTION; BOOKKEEPING FOR MINJT, MAXJT IF (MINJT(INRG).LT.0) MINJT(INRG)=JTOT MAXJT(INRG)=MAX0(MAXJT(INRG),JTOT) IJ=0 C LET II MEASURE INITIAL AND JJ MEASURE FINAL QUANT. STATES XJT=(2*JTOT+1)*PI IF (IP.LT.14) THEN IF (IEXCH.NE.0) XJT=XJT*WT ELSE IF (WT.NE.0.D0) XJT=XJT*WT ENDIF DMX=0.D0 OMX=0.D0 DO 3000 II=1,NOPEN C LEVI/LEVJ ARE 'LEVEL NO.'; IXI/IXJ 'SIG INDX' FROM JLEV(LEV,NQN) C BEGINNING IN V14, NEGATIVE SIG INDX FOR "INCOMPLETE" CS VALUES LEVI=J(II) IXI=JLEV(NLEV*(NQN-1)+LEVI) LOGI=IXI.LT.0 IXI=ABS(IXI) DO 3000 JJ=1,NOPEN LEVJ=J(JJ) IXJ=JLEV(NLEV*(NQN-1)+LEVJ) LOGJ=IXJ.LT.0 IXJ=ABS(IXJ) FACT=XJT/(WVEC(JJ)*WVEC(JJ)) IJ=IJ+1 IF (II.EQ.JJ) SREAL(IJ)=1.D0-SREAL(IJ) C INDEX OF SIG(JJ,II) IS ISTSIG(INRG)+NTOP*(IXJ-1)+IXI IX=ISTSIG(INRG)+NTOP*(IXJ-1)+IXI ADD = (SREAL(IJ)*SREAL(IJ)+SIMAG(IJ)*SIMAG(IJ))*FACT C GET DEGENERACY DENOMINATOR IF (ITYPE.EQ.3) THEN C BELOW IS ITYPE=3 W/ TAKAYANAGI IDENTICAL PARTICLE COUNTING JI1=JLEV(LEVJ) JI2=JLEV(LEVJ+NLEV) DEGEN=DBLE((2*JI1+1)*(2*JI2+1)) IF (IDENT.NE.0) THEN JF1=JLEV(LEVI) JF2=JLEV(LEVI+NLEV) IF (JI1.EQ.JI2) DEGEN=DEGEN/2.D0 IF (JF1.EQ.JF2) DEGEN=DEGEN/2.D0 ENDIF ELSE C FOLLOWING APPLIES TO LINEARS AND TO SYM/ASYM TOPS DEGEN=DBLE(2*JLEV(LEVJ)+1) ENDIF ADD=ADD/DEGEN FACT=1.D0 IF (LOGI.AND.LOGJ) FACT=-FACT SIG(IX)=SIG(IX)+ADD*FACT IF (II.EQ.JJ) OMX=MAX(OMX,ADD) IF (II.NE.JJ) DMX=MAX(DMX,ADD) 3000 CONTINUE IF (JTOUT.GE.999999.OR.(JTOUT.LT.999999.AND.JTOT.GT.JTOUT)) 1 WRITE(6,653) OMX,DMX 653 FORMAT(10X,'MAX CHANGE IN DIAG, OFF-DIAG SIGMA',1P,2D12.4) C GO BACK FOR MORE JTOT, INRG SETS . . . GO TO 2000 C C END OF FILE ON (10) 9000 WRITE(6,695) IT 695 FORMAT(' *** NOTE. NORMAL END OF FILE REACHED ON UNIT (',I3,')') GO TO 9123 9090 WRITE(6,690) IT 690 FORMAT(' ***** ***** ABNORMAL EOF REACHED ON UNIT (',I3,')') C C PRINT CROSS SECTIONS. CODE REMAINS FOR TOTALS OUT OF LEVEL 9123 JSTOUT=JSTEP*JSTEPT XJSTEP=JSTOUT WRITE(6,123) JSTEPT,JSTEP,XJSTEP 123 FORMAT(/' *****'/' JSTEP ON TAPE =',I4/ 1 ' JSTEP REQUESTED IN &SAVE INPUT =',I4/ 2 ' CROSS SECTIONS MULTIPLIED BY',F9.2/' *****') IF (LFORM) WRITE(6,702) 702 FORMAT(/12X,'ENERGY JTOTL JSTEP JTOTU',9X, 1 'F I',10X,'SIG(F,I)'/) DO 4000 I=1,NNRG DO 4100 II=1,NTOP EREL=ENERGY(I)-ELEVEL(II) SUMII=0. DO 4101 JJ=1,NTOP IST=ISTSIG(I)+(II-1)*NTOP+JJ SIG(IST)=SIG(IST)*XJSTEP C ***OFF-DIAGONAL**** TOTAL OUT OF LEVEL II IF (II.NE.JJ.AND.SIG(IST).GT.0.D0) SUMII=SUMII+SIG(IST) XST=BL IF (SIG(IST).LT.0.D0) XST=ST IF (ABS(SIG(IST)).LT.EPS) GO TO 4101 IF (LFORM) THEN C FORMAT BELOW MATCHES MOLSCAT VERSION 14 OUTPUT ROUTINE WRITE(6,703)BL,ENERGY(I),MINJT(I),JSTOUT,MAXJT(I), 1 JJ,II,ABS(SIG(IST)),XST 703 FORMAT(A1,F19.6,I5,2I7,5X,2I5,1P,D20.6,1X,A1) ELSE WRITE(6,700) ENERGY(I),MINJT(I),JSTOUT,MAXJT(I),JJ,II,SIG(IST) 2 ,COUT,NLEV 700 FORMAT(' ',F19.6,I5,'(',I3,')',I5,5X,2I5,E20.6,5X,A1,I2) ENDIF 4101 CONTINUE IF (LPRINT) WRITE(6,692) EREL,II,SUMII 692 FORMAT(' *** REL E =',F10.2,' SUM OUT OF LEVEL',I4,' =',F8.3) 4100 CONTINUE 4000 CONTINUE STOP END SUBROUTINE IDCHK(IDENT,NLEV,JLEV) DIMENSION JLEV(NLEV) LOGICAL J1GEJ2,J1LEJ2,J1AND2 WRITE(6,600) 600 FORMAT('0 ***'/' *** FOR MOD(ITYPE,10) = 3 ATTEMP TO DETERMINE' 1 ,' IDENT FROM JLEV DATA'/' ***') J1GEJ2=.TRUE. J1LEJ2=.TRUE. DO 1000 I=1,NLEV J1=JLEV(I) J2=JLEV(NLEV+I) J1GEJ2=J1GEJ2.AND.J1.GE.J2 1000 J1LEJ2=J1LEJ2.AND.J1.LE.J2 IF (J1LEJ2.AND.J1GEJ2) THEN WRITE(6,601) 601 FORMAT('0 J1=J2 FOR ALL BASIS FUNCTIONS ', 1 ' -- ASSUME NON-IDENTICAL') IDENT=0 GO TO 9000 ENDIF IF ((.NOT.J1LEJ2).AND.(.NOT.J1GEJ2)) THEN WRITE(6,602) 602 FORMAT('0 BASIS NOT ORDERED J1.GT.J2 OR J1.LE.J2', 1 ' -- ASSUME NON-IDENTICAL') IDENT=0 GO TO 9000 ENDIF C IF WE REACH BELOW, EITHER J1.GE.J2 **OR** J1.LE.J2 FOR ALL BASIS WRITE(6,603) 603 FORMAT('0 BASIS IS ORDERED EITHER J1.GT.J2 OR J1.LE.J2', 1 ' -- ASSUME IDENTICAL PARTICLES') IDENT=1 9000 WRITE(6,604) IDENT 604 FORMAT('0 *** IDENT HAS BEEN SET TO',I3) RETURN END SUBROUTINE SWRITE(IU,N,S) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION S(N,N) C WRITE(IU) ((S(I,J),J=1,I),I=1,N) RETURN C ENTRY SREAD(IU,N,S,IEND) IEND=0 READ(IU,END=9999) ((S(I,J),J=1,I),I=1,N) DO 1000 I=1,N DO 1000 J=1,I-1 1000 S(J,I)=S(I,J) RETURN 9999 IEND=1 RETURN END