|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 30720 (0x7800)
Types: TextFile
Names: »per10«
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
└─⟦4334b4c0b⟧
└─⟦this⟧ »per10«
c program 10 acrz
c
c acrzgf values. oscillator strengths from nomerical mchf radial
c funktions, with an adaption to calculate electric quadrupole
c oscillator strengths. c.f.fischer, k.m.s.saxena, m.godefroid.
c ref. in comp. phys. commun. 9 (1975) 381 and
c comp. phys. commun. 15 (1978) 275 and
comp. phys. commun. 17 (1979) 426.
c
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C * *
C * JULY 1974 *
C * *
C * PROGRAM TO EVALUATE LENGTH AND VELOCITY FORM OF 'GF' VALUES *
C * USING NUMERICAL MCHF RADIAL FUNCTIONS. *
C * *
C * *
C * BY *
C * *
C * *
C * C. FROESE FISCHER AND K.M.S. SAXENA *
C * *
C * DEPARTMENT OF APPLIED MATHEMATICS *
C * UNIVERSITY OF WATERLOO, WATERLOO, ONTARIO N2L 3G1 *
C * C A N A D A. *
C * *
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C * *
C * DESCRIPTION OF THE COMMON BLOCKS:- *
C * ================================ *
C * *
C * (1) COMMON BLOCK /PARAM/. *
C * -------------------- *
C * *
C * THIS COMMON BLOCK STORES SEVERAL NUMERICAL CONSTANTS *
C * IN DOUBLE PRECISION AND VALUES OF CERTAIN PARAMETERS OF THE *
C * CALCULATIONS. *
C * *
C * THE VARIABLES IN THIS COMMON BLOCK ARE AS FOLLOWS. *
C * *
C * D0 : DOUBLE PRECISION CONSTANT 0 *
C * D1 : DOUBLE PRECISION CONSTANT 1 *
C * D2 : DOUBLE PRECISION CONSTANT 2 *
C * D3 : DOUBLE PRECISION CONSTANT 3 *
C * D4 : DOUBLE PRECISION CONSTANT 4 *
C * D5 : DOUBLE PRECISION CONSTANT 1/2 *
C * D6 : DOUBLE PRECISION CONSTANT 6 *
C * D10: DOUBLE PRECISION CONSTANT 10 *
c * d15: double precision constant 15 *
c * d30: double precision constant 30 *
C * H : 0.0625, THE STEP-SIZE DELTA(RHO) *
C * H1 : H/3 *
C * Z : ATOMIC NUMBER *
C * NO : MAXIMUM NUMBER OF POINTS IN THE RANGE OF THE FUNCTIONS *
C * =220 *
C * ND : NO-2 *
C * *
C * MOST OF THESE CONSTANTS ARE SET IN THE BLOCK DATA *
C * SEGMENT. *
C * *
C * *
C * (2) THE BLANK COMMON //. *
C * ------------------- *
C * *
C * THIS COMMON SECTION STORES THE INFORMATION ABOUT *
C * THE NUMERICAL MCHF WAVEFUNCTIONS. *
C * *
C * THE VARIABLES IN THIS COMMON BLOCK ARE AS FOLLOWS. *
C * *
C * R(220) : VALUES OF R(J)=EXP(RHO(J))/Z *
C * RR(220) : VALUES OF R(J)**2 *
C * R2(220) : VALUES OF DSQRT(R(J)) *
C * P(50,220) : VALUES OF P(R(J))/R2(J) FOR ORBITALS I=1,50 *
C * AZ(50) : A0=LIMIT OF ( P(R)/(R**(L+1)) ) AS R->0 FOR *
C * ORBITALS I=1,50 *
C * L(50) : L-QUANTUM NUMBERS FOR ORBITALS I=1,50 *
C * MAX(50) : NUMBER OF MESH-POINTS IN THE NUMERICAL *
C * ORBITALS I=1,50 *
C * *
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C * *
C * DESCRIPTION OF THE OTHER IMPORTANT VARIABLES:- *
C * ============================================ *
C * *
C * ASTER : AN INTEGER VARIABLE READ IN FORMAT(A1). *
C * BLANK : A BLANK CHARACTER ' ', DEFINED AS AN INTEGER *
C * AND SET BY A DATA STATEMENT IN THE MAIN PROGRAM. *
C * SYMM(5) : AN INTEGER VARIABLE STORING THE SPECTROSCOPIC *
C * SYMBOLS 'S','P','D','F', AND 'G' IN FORMAT(A1). *
C * SYMM(I) CORRESPONDS TO THE SYMMETRY WITH L=(I-1). *
C * JREAD : GENERAL INPUT CHANNEL FOR A SET OF CALCULATIONS. *
C * IREAD1 : INPUT CHANNEL FOR THE WAVEFUNCTIONS REQUIRED IN *
C * A CASE. *
C * IREAD2 : INPUT CHANNEL FOR THE TRANSITION INTEGRALS DATA *
C * REQUIRED IN A CASE. *
C * IWRITE : OUTPUT CHANNEL REQUIRED TO OUTPUT THE RESULTS OF *
C * A CASE. *
C * ATOM : ATOMIC SYMBOL/DESIGNATION IN FORMAT(A8). *
C * IZ : ATOMIC NUMBER Z. *
C * NWF : NUMBER OF ORBITALS INVOLVED. *
C * TRMI : TERM SYMBOL FOR THE INITIAL STATE IN FORMAT(A8). *
C * TRMF : TERM SYMBOL FOR THE FINAL STATE IN FORMAT(A8). *
C * NCFGI : NUMBER OF CONFIGURATIONS IN THE INITIAL STATE. *
C * NCFGF : NUMBER OF CONFIGURATIONS IN THE FINAL STATE. *
C * MULT : MULTIPLICITY (2S+1) OF THE TWO STATES. *
c * kctas : phase convention *
c * = 0 fano and racah *
c * = 1 condon and shortley *
C * EL(50,3) : AN INTEGER VARIABLE STORING THE THREE CHARACTER *
C * ALPHANUMERIC SYMBOL OF THE ORBITALS I=1,50 IN *
C * FORMAT(3A1) FOR EACH OF THEM. *
C * ETI : TOTAL ENERGY (IN A.U.) FOR THE INITIAL STATE. *
C * ETF : TOTAL ENERGY (IN A.U.) FOR THE FINAL STATE. *
C * CONFAI(30): CHARACTERS 01-08 (IN A8 FORMAT) OF THE ALPHA- *
C * NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF *
C * THE INITIAL STATE CONFIGURATIONS I=1,30. *
C * CONFBI(30): CHARACTERS 09-16 (IN A8 FORMAT) OF THE ALPHA- *
C * NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF *
C * THE INITIAL STATE CONFIGURATIONS I=1,30. *
C * CONFCI(30): CHARACTERS 17-24 (IN A8 FORMAT) OF THE ALPHA- *
C * NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF *
C * THE INITIAL STATE CONFIGURATIONS I=1,30. *
C * CONFAF(30): CHARACTERS 01-08 (IN A8 FORMAT) OF THE ALPHA- *
C * NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF *
C * THE FINAL STATE CONFIGURATIONS I=1,30. *
C * CONFBF(30): CHARACTERS 09-16 (IN A8 FORMAT) OF THE ALPHA- *
C * NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF *
C * THE FINAL STATE CONFIGURATIONS I=1,30. *
C * CONFCF(30): CHARACTERS 17-24 (IN A8 FORMAT) OF THE ALPHA- *
C * NUMERIC SYMBOL (IN TOTAL OF 24 CHARACTERS) OF *
C * THE FINAL STATE CONFIGURATIONS I=1,30. *
C * WTI(30) : WEIGHTS OF THE INITIAL STATE CONFIGURATIONS I=1,30*
C * WTF(30) : WEIGHTS OF THE FINAL STATE CONFIGURATIONS I=1,30*
C * COEF : MULTILPYING COEFFICIENT IN A TRANSITION INTEGRAL. *
c * ktens : the order of the tensor *
c * = 1 dipole transition *
c * = 2 quadropol transition *
C * KRHO : THE NUMBER INDICATING THE ORBITAL EL(KRHO) FROM *
C * THE INITIAL STATE CONFIGURATION APPEARING IN A *
C * TRANSITION INTEGRAL. *
C * KSIG : THE NUMBER INDICATING THE ORBITAL EL(KSIG) FROM *
C * THE FINAL STATE CONFIGURATION APPEARING IN A *
C * TRANSITION INTEGRAL. *
C * JI,JF : THE NUMBERS INDICATING THE INITIAL AND FINAL *
C * STATE CONFIGURATIONS FOR A TRANSITION INTEGRAL. *
C * II1,II2, : THE NUMBERS INDICATING THE ORBITALS EL(II1), *
C * II3 EL(II2), AND EL(II3) FROM THE INITIAL STATE *
C * CONFIGURATION JI APPEARING IN THE THREE OVERLAP *
C * INTEGRALS IN A TRANSITION INTEGRAL CONTRIBUTION. *
C * IF1,IF2, : THE NUMBERS INDICATING THE ORBITALS EL(IF1), *
C * IF3 EL(IF2), AND EL(IF3) FROM THE FINAL STATE *
C * CONFIGURATION JF APPEARING IN THE THREE OVERLAP *
C * INTEGRALS IN A TRANSITION INTEGRAL CONTRIBUTION. *
C * IP1,IP2, : THE NUMBERS INDICATING THE POWERS OF THE THREE *
C * IP3 OVERLAP INTEGRALS IN A TRANSITION INTEGRAL *
C * CONTRIBUTION. *
C * *
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C * *
C * NON STANDARD FORTRAN:- *
C * ==================== *
C * *
C * THIS PROGRAM IS WRITTEN AS A SYSTEM 360 DOUBLE PRE- *
C * CISION PROGRAM. HOWEVER, ON COMPUTERS WITH A WORD LENGTH OF *
C * 48 BITS OR MORE IT SHOULD BE USED IN SINGLE PRECISION ONLY. *
C * IN ORDER TO CONVERT THIS PROGRAM TO A SINGLE PRECISION PROGRAM *
C * THE FOLLOWING SHOULD BE DONE: *
C * *
C * (1) REMOVE ALL 'IMPLICIT REAL*8(A-H,O-Z)' CARDS. *
C * (2) REMOVE '*8' FROM THE FUNCTION DEFINITION CARDS, FOR *
C * EXAMPLE, 'REAL FUNCTION GRAD*8(I,J)' WILL BE REPLACED *
C * BY 'REAL FUNCTION GRAD(I,J)' *
C * (3) CHANGE SYSTEM FUNCTION NAMES 'DSQRT','DABS', ETC. TO *
C * 'SQRT','ABS', ETC. *
C * (4) CHANGE DOUBLE PRECISION CONSTANTS TO SINGLE PRECISION *
C * CONTANTS. *
C * (5) CHANGE D-FORMAT CODES TO E-FORMAT CODES. *
C * *
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C ------------------------------------------------------------------
C *** M A I N
C ------------------------------------------------------------------
C
c implicit real*8(a-h,o-z)
INTEGER ASTER,BLANK,SYMM(5),EL(50,3)
COMMON /PARAM/D0,D1,D2,D3,D4,D5,D6,D10,d15,d30,H,H1,Z,NO,ND
COMMON R(220),RR(220),R2(220),P(50,220),AZ(50),L(50),MAX(50)
DIMENSION CONFAI(30),CONFBI(30),CONFCI(30),WTI(30),
1 CONFAF(30),CONFBF(30),CONFCF(30),WTF(30)
DATA BLANK/1H /
DATA SYMM(1),SYMM(2),SYMM(3),SYMM(4),SYMM(5)/1HS,1HP,1HD,1HF,1HG/
C
C *** THE INPUT FORMATS. IT IS TO BE NOTED HERE THAT THE '102 FORMAT'
C *** USED TO INPUT THE NUMERICAL MCHF RADIAL FUNCTION DATA IS NOT THE
C *** SAME AS ACTUALLY USED IN THE MCHF PROGRAM TO PUNCH THESE RADIAL
C *** FUNCTIONS. THIS HAS BEEN DONE IN ORDER TO KEEP THE DATA-FIELD
C *** LIMITED TO COLUMNS 1-72 AND USE COLUMNS 73-80 OF THE CARDS FOR
C *** IDENTIFICATION AND SEQUENCE NUMBERS ONLY. HOWEVER, ONE CAN USE
C *** THE MCHF RADIAL FUNCTION DATA AS IT IS PUNCHED OUT IN THE MCHF
C *** PROGRAM BY CHANGING THE '102 FORMAT' TO:
C102 FORMAT(32X,3A1,I6,22X,D16.8/(7F11.7))
C
101 FORMAT(A1,3I4,A8,2I4,2A8,5I4)
102 FORMAT( 3A1,I6,22X,D16.8/(6F11.7))
103 FORMAT(F16.8)
104 FORMAT(3A8,F12.8)
105 FORMAT(A1,F12.8,3HRI(,2I2,1H,,2I2,1H),9I4)
106 format(i4)
C
C *** THE OUTPUT FORMATS.
C
200 FORMAT(5H1SET ,I3/5H CASE,I3/8H *******//)
201 FORMAT(23X37HOSCILLATOR STRENGTHS CALCULATION FOR , A8,4H (Z=,I3,2
1H) ,A8,2H->,A8,12H TRANSITION )
202 FORMAT(1H ,22X,83(1H*)/)
203 FORMAT(/9H JREAD =,I3,38H (GENERAL DATA INPUT UNIT FOR THE SET) /
19H IREAD1 =,I3,45H (WAVEFUNCTION DATA INPUT UNIT FOR THIS CASE) /
29H IREAD2 =,I3,53H (TRANSITION INTEGRALS DATA INPUT UNIT FOR THIS
2CASE) /
39H IWRITE =,I3,36H (PRINTED OUTPUT UNIT FOR THIS CASE) /
49H NWF =,I3,35H (NUMBER OF WAVEFUNCTIONS INVOLVED) /
59H NCFGI =,I3,48H (NUMBER OF CONFIGURATIONS IN THE INITIAL STATE)
6/9H NCFGF =,I3,46H (NUMBER OF CONFIGURATIONS IN THE FINAL STATE)/
79H MULT =,I3,40H (MULTIPLICITY (2S+1) OF THE TWO STATES)
8/9h ictas =,i3,20h (phase conventions) )
204 FORMAT(//5H THE ,I2,28H ORBITALS ARE THE FOLLOWING:/)
205 FORMAT(2X,9(I1,3X),21(I2,2X))
206 FORMAT(/2X,30(3A1,1X)/)
207 FORMAT(2X,30(I2,2X))
208 FORMAT(//38H CONFIGURATIONS/WEIGHTS OF THE INITIAL,A8,
124HSTATE WITH TOTAL ENERGY= ,F16.8)
209 FORMAT( 86H -----------------------------------------------------
1--------------------------------/)
210 FORMAT(I7,1H.,3A8,F12.7)
211 FORMAT(//38H CONFIGURATIONS/WEIGHTS OF THE FINAL,A8,
124HSTATE WITH TOTAL ENERGY= ,F16.8)
212 FORMAT(64H1TERM INITIAL STATE CONFIGURATION FINAL STATE CONFI
1GURATION11X12HCONTRIBUTION18X6HLENGTH8X8HVELOCITY)
213 FORMAT(128H ---- --------------------------- ----------------
1--------- ------------------------- ---------- ---
2-------//)
214 FORMAT(//91H ERROR ( L(KRHO)-L(KSIG).NE.+1/-1 ) FOUND IN THE TRAN
1SITION INTEGRAL PART OF THE FOLLOWING,I4,32H-TH INPUT CARD FROM IR
2EAD2 FILE./)
215 FORMAT(//70H ERROR (DIFFERENT L VALUES) FOUND IN THE OVERLAP PART
1OF THE FOLLOWING,I4,32H-TH INPUT CARD FROM IREAD2 FILE./)
216 FORMAT( 2H (,I2,2H) ,I4,1H.,3A8,4H -> ,I2,1H.,3A8,4H : <,3A1,3H/O/
1,3A1,1H>,17X,1H:,2X,F13.8,2X,F13.8)
217 FORMAT(1H+,79X,3H* <,3A1,1H/,3A1,3H>**,I1)
218 FORMAT(1H ,79X,3H* <,3A1,1H/,3A1,3H>**,I1)
219 FORMAT(/35X40HFINAL OSCILLATOR STRENGTHS (LENGTH) =,d14.8/
1 63X12H(VELOCITY) =, d14.8)
220 FORMAT(/42X35HENERGY DIFFERENCE OF THE STATES = ,D16.8,5H CM-1/
174X3H= ,D16.8,10H ANGSTROMS/74X3H= ,D16.8,6H A. U.)
221 FORMAT(//91H ERROR ( L(KRHO)-L(KSIG).NE.+2/-2/0) FOUND IN THE TRAN
1SITION INTEGRAL PART OF THE FOLLOWING,I4,32H-1H INPUT CARD FROM IR
1EAD2 FILE./)
222 FORMAT(//37H THE ORDER OF THE TENSOR OPERATOR IS ,I4,//)
C
C *** SET JREAD EQUAL TO THE CARD INPUT CHANNEL OF YOUR INSTALLATION
C *** & NSET=0, AND START THE DATA READING FOR THE FISRT SET OF THE
C *** CALCULATIONS.
C
JREAD=1
NSET=0
c* set input and output channels
zone readf(1200,1,stderror)
zone readf1(200,1,stderror)
zone readf2(200,1,stderror)
zone writef(200,1,stderror)
call zassign(readf,1)
call zassign(readf1,2)
call zassign(readf2,3)
call zassign(writef,7)
call open(readf,4,'readfile',0)
call open(readf1,4,'readfile1',0)
call open(readf2,4,'readfile2',0)
C
C
C *** START THE DATA READING FOR A NEW SET OF CALCULATIONS.
C
1 READ(JREAD,101) ASTER, JREAD
C
C *** IF 'ASTER' IS NOT A BLANK CHARACTER, ALL SETS OF CALCULATIONS ARE
C *** FINISHED AND THE PROGRAM STOPS.
C *** IF 'ASTER' IS A BLANK CHARACTER, SET NSET=NSET+1 & NCASE=0,
C *** AND START THE DATA READING FOR THE FIRST CASE OF THIS SET.
C
IF(ASTER.NE.BLANK) GO TO 15
NSET=NSET+1
NCASE=0
C
C *** START THE DATA READING FOR A NEW CASE OF THE NSET-TH SET.
C
2 READ(JREAD,101) ASTER,IREAD1,IREAD2,IWRITE,ATOM,IZ,NWF,TRMI,TRMF,
1 NCFGI,NCFGF,MULT,ISO,icstas
C
C *** IF 'ASTER' IS NOT A BLANK CHARACTER, ALL CASES OF CALCULATIONS ARE
C *** FINISHED AND THE PROGRAM PROCEEDS TO TAKE THE NEXT SET OF CASES.
C *** IF 'ASTER' IS A BLANK CHARACTER, SET NCASE=NCASE+1, AND PROCEED
C *** FOR READING ADDITIONAL DATA, AND CALCULATIONS ETC. FOR THIS CASE.
C
IF(ASTER.NE.BLANK) GO TO 1
NCASE=NCASE+1
READ(JREAD,103) ETI
READ(JREAD,104) (CONFAI(NC),CONFBI(NC),
1 CONFCI(NC),WTI(NC),NC=1,NCFGI)
READ(JREAD,103) ETF
READ(JREAD,104) (CONFAF(NC),CONFBF(NC),
1 CONFCF(NC),WTF(NC),NC=1,NCFGF)
WRITE(IWRITE,200) NSET,NCASE
Z=IZ
WRITE(IWRITE,201) ATOM,IZ,TRMI,TRMF
WRITE(IWRITE,202)
WRITE(IWRITE,203) JREAD,IREAD1,IREAD2,IWRITE,NWF,NCFGI,NCFGF,MULT
C
C *** CALCULATE THE 'R' VARIABLE-MESH FOR THE NUMERICAL RADIAL ORBITALS
C *** P(R).
C
RHO=-4.0
DO 3 J=1,NO
R(J)=eXP(RHO)/Z
RR(J)=R(J)*R(J)
R2(J)=SQRT(R(J))
3 RHO=RHO+H
C
C *** READ THE RADIAL ORBITALS AND CALCULATE THEIR L-QUANTUM NUMBERS.
C
DO 5 I=1,NWF
READ(IREAD1,102) (EL(I,J),J=1,3),M,AZ(I),(P(I,J),J=1,M)
J=3
IF(EL(I,1).NE.BLANK) J=2
IF(EL(I,J).EQ.SYMM(1)) L(I)=0
IF(EL(I,J).EQ.SYMM(2)) L(I)=1
IF(EL(I,J).EQ.SYMM(3)) L(I)=2
IF(EL(I,J).EQ.SYMM(4)) L(I)=3
IF(EL(I,J).EQ.SYMM(5)) L(I)=4
MAX(I)=M
IF(M.EQ.NO) GO TO 5
M=M+1
DO 4 J=M,NO
4 P(I,J)=D0
5 CONTINUE
read(iread2,106) ktens
write(iwrite,222) ktens
WRITE(IWRITE,204) NWF
N30=MIN0(NWF,30)
WRITE(IWRITE,205) (I,I=1,N30)
WRITE(IWRITE,206) ((EL(I,J),J=1,3),I=1,N30)
IF(NWF.LE.30) GO TO 6
WRITE(IWRITE,207) (I,I=31,NWF)
WRITE(IWRITE,206) ((EL(I,J),J=1,3),I=31,NWF)
6 WRITE(IWRITE,208) TRMI,ETI
WRITE(IWRITE,209)
WRITE(IWRITE,210) (NC,CONFAI(NC),CONFBI(NC),
1 CONFCI(NC),WTI(NC),NC=1,NCFGI)
WRITE(IWRITE,211) TRMF,ETF
WRITE(IWRITE,209)
WRITE(IWRITE,210) (NC,CONFAF(NC),CONFBF(NC),
1 CONFCF(NC),WTF(NC),NC=1,NCFGF)
WRITE(IWRITE,212)
WRITE(IWRITE,213)
C
C *** INITIALIZE THE TOTAL LENGTH (SL) AND THE TOTAL VELOCITY (SV)
C *** 'GF' VALUES TO ZERO, SET NTERM=0, AND START READING THE TRANSITION
C *** INTEGRALS DATA FOR THE CASE IN QUESTION.
C
SL = D0
SV = D0
NTERM=0
C
C *** READ DATA FOR THE NEXT TRANSITION INTEGRAL TERM.
C
7 READ(IREAD2,105) ASTER,COEF,KRHO,JI,KSIG,JF,II1,IF1,IP1,
1 II2,IF2,IP2,
2 II3,IF3,IP3
C
C *** IF 'ASTER' IS NOT A BLANK CHARACTER, ALL TRANSITION INTEGRAL
C *** TERMS HAVE BEEN READ AND USED. THEN THE PROGRAM PROCEEDS TO DO
C *** THE CALCULATIONS OF THE FINAL RESULTS OF OSCILLATOR STRENGTHS
C *** FOR THE CASE IN QUESTION.
C *** IF 'ASTER' IS A BLANK CHARACTER, SET NTERM=NTERM+1, AND PROCEED
C *** TO EVALUATE THE CONTRIBUTIONS TL AND TV OF THIS TERM TO THE
C *** TOTAL 'GF' VALUES SL AND SV.
C
IF(ASTER.NE.BLANK) GO TO 14
NTERM=NTERM+1
C
C *** BEFORE PROCEEDING TO USE THIS TRANSITION INTEGRAL INPUT DATA,
C *** CHECK FOR ANY OBVIOUS ERROR IN IT. IF ERROR IS FOUND, IT IS
C *** PRINTED AND THE PROGRAM STOPS.
C
go to (25,26),ktens
25 IF(IABS(L(KRHO)-L(KSIG)).EQ.1) GO TO 8
C
C *** ERROR FOUND IN AN INPUT DATA FROM FILE IREAD2.
C
WRITE(IWRITE,214) NTERM
GO TO 12
26 if(iabs(l(krho))-l(ksig)).eq.2.or.(l(krho)-l(ksig)).eq.0) go to 8
write(iwrite,221) nterm
go to 12
8 IF(IP1.EQ.0) GO TO 9
IF(L(II1).EQ.L(IF1)) GO TO 9
C
C *** ERROR FOUND IN AN INPUT DATA FROM FILE IREAD2.
C
GO TO 11
9 IF(IP2.EQ.0) GO TO 10
IF(L(II2).EQ.L(IF2)) GO TO 10
C
C *** ERROR FOUND IN AN INPUT DATA FROM FILE IREAD2.
C
GO TO 11
10 IF(IP3.EQ.0) GO TO 13
IF(L(II3).EQ.L(IF3)) GO TO 13
C
C *** ERROR FOUND IN AN INPUT DATA FROM FILE IREAD2.
C
11 WRITE(IWRITE,215) NTERM
12 WRITE(IWRITE,105) ASTER,COEF,KRHO,JI,KSIG,JF,II1,IF1,IP1,
1 II2,IF2,IP2,
2 II3,IF3,IP3
CALL EXIT
C
C *** NO ERROR IN THE TRANSITION INTEGRAL DATA IS FOUND AND HENCE
C *** PROCEED TO USE IT FOR EVALUATING TL AND TV VALUES.
C
13 goto(27,28) ktens
27 k=2
l1=l(krho)
L2 = L(KSIG)
LL = MAX0(L1,L2)
if(ictas.eq.1) k=ll-l1
d=(-1)**k*sqrt(ll*d1)*wti(ji)*wtf(jf)*coef
goto 29
l1=l(krho)
L2=L(KSIG)
IF(L1.EQ.L2) GOTO 30
LL=MAX0(L1,L2)
XX=D3*LL*(LL-D1)
XX=XX/(D2*(D2*LL-D1))
d=sqrt(xx)*wti(ji)*wtf(jf)*coef
GOTO 29
30 k=2
1 if(icstas.eq.1) k=1
xx=d1*l1*(l1+d1)*(d2*l1+d1)
XX=XX/(D2*L1+D3)*(D2*L1-D1)
D=(-1)**K*sqrt(XX)*WTI(JI)*WTF(JF)*COEF
C
c
c*** multiply by the overlaps integral if any
c
29 IF (IP1 .NE. 0) D = D*QUADR(II1,IF1,0)**IP1
IF (IP2 .NE. 0) D = D*QUADR(II2,IF2,0)**IP2
IF (IP3 .NE. 0) D = D*QUADR(II3,IF3,0)**IP3
goto(31,32),ktens
31 TL = D*QUADR(KRHO,KSIG,1)
TV = D*GRAD(KRHO,KSIG)
goto 33
32 tl=d*quadr(krho,ksig,2)
tv=d*grad2(krho,ksig)
C
C *** OUTPUT THE CALCULATED TERM CONTRIBUTIONS TL AND TV.
C
33 WRITE(IWRITE,216) NTERM,JI,CONFAI(JI),CONFBI(JI),CONFCI(JI),
1 JF,CONFAF(JF),CONFBF(JF),CONFCF(JF),
2 (EL(KRHO,J),J=1,3),(EL(KSIG,J),J=1,3),TL,TV
IF(IP1.NE.0) WRITE(IWRITE,217) (EL(II1,J),J=1,3),(EL(IF1,J),J=1,3)
1 ,IP1
IF(IP2.NE.0) WRITE(IWRITE,218) (EL(II2,J),J=1,3),(EL(IF2,J),J=1,3)
1 ,IP2
IF(IP3.NE.0) WRITE(IWRITE,218) (EL(II3,J),J=1,3),(EL(IF3,J),J=1,3)
1 ,IP3
C
C *** ADD THE EVALUATED TERM CONTRIBUTIONS TL AND TV TO THE TOTAL
C *** 'GF' VALUES SL AND SV.
C
SL = SL + TL
SV = SV + TV
C
C *** PROCEED TO TAKE THE NEXT TRANSITION INTEGRAL DATA.
C
GO TO 7
C
C *** ALL TRANSITION INTEGRALS DATA HAVE BEEN USED AND HENCE PROCEED
C *** TO DO THE FINAL CALCULATIONS AND OUTPUT THEM.
C
14 D = abs(ETI-ETF)
goto(34,35),ktens
34 CL = D2*D*abs(DFLOAT(MULT))/D3
CV = D2*abs(DFLOAT(MULT))/(D*D3)
goto 36
35 cl=d1*d**3*abs(float(mult))/d30
cv=d2*d*abs(float(mult))/d15
36 gl=cl*sl**2
GV = CV*SV**2
if(ktens.ne.2) goto 37
gl=gl*5.325132272e-05
gv=gv*5.325132272e-05
37 write(iwrite,219) gl,gv
DD = D*D2*109737.3e0
ANGS = D10**8/DD
WRITE(IWRITE,220) DD,ANGS,D
C
C *** CALCULATIONS FOR THE PRESENT CASE ARE FINISHED AND THE FINAL
C *** RESULTS HAVE ALSO BEEN OUTPUT.
C *** IF THE NEXT CASE TO BE DONE IS A CASE ISO-ELECTRONIC TO THE ONE
C *** JUST FINISHED, REWIND THE DEVICE 'IREAD2'.
C *** THIS WILL ENABLE US TO REUSE THE TRANSITION INTEGRALS DATA.
C *** PROCEED TO TAKE THE NEXT CASE.
C
if(iso.eq.1) call setposition(readf2,0,0)
GO TO 2
15 STOP
END
C
C ------------------------------------------------------------------
C *** B L O C K D A T A
C ------------------------------------------------------------------
C
subroutine bldata
c IMPLICIT REAL*8(A-H,O-Z)
COMMON /PARAM/D0,D1,D2,D3,D4,D5,D6,D10,d15,d30,H,H1,Z,NO,ND
data d0,d1,d2,d3,d4,d5/0.e0,1.e0,2.e0,3.e0,4.e0,.5e0/
data d6,d10,d15,d30/6.e0,10.e0,15.e0,30.e0/
DATA H,H1/.0625e0,.4166666666666667e-01/
DATA NO,ND/220,218/
END
C
C ------------------------------------------------------------------
C *** Q U A D R
C ------------------------------------------------------------------
C
C *** QUADR INTEGRATES P(I)*P(J)*R**KK BY SIMPSON'S RULE
C
REAL FUNCTION QUADR(I,J,KK)
c IMPLICIT REAL*8(A-H,O-Z)
COMMON /PARAM/D0,D1,D2,D3,D4,D5,D6,D10,d15,d30,H,H1,Z,NO,ND
COMMON R(220),RR(220),R2(220),P(50,220),AZ(50),L(50),MAX(50)
K = KK + 2
LI = L(I)
LJ = L(J)
DEN = LI + LJ + 1 + K
ZR = Z*R(4)
BI = (P(I,4)/(AZ(I)*R2(4)*R(4)**LI) - D1+ZR/(LI+1) )/ZR**2
BJ = (P(J,4)/(AZ(J)*R2(4)*R(4)**LJ) - D1+ZR/(LJ+1) )/ZR**2
ALPHA= (D1/(LI + 1) + D1/(LJ + 1))/(DEN + D1)
ZR = Z*R(1)
BETA = (DEN+D1)*ALPHA**2 - D2*(BI+BJ+D1/((LI+1)*(LJ+1)))/(DEN+D2)
D = P(I,1)*P(J,1)*R(1)**K*(((BETA*ZR+ALPHA)*ZR+D1)/(DEN*H1)+D5)
M = MIN0(MAX(I),MAX(J)) - 1
DO 1 JJ = 2,M,2
JP = JJ + 1
1 D = D + D2*P(I,JJ)*P(J,JJ)*R(JJ)**K+P(I,JP)*P(J,JP)*R(JP)**K
QUADR = D*H1
RETURN
END
C
C ------------------------------------------------------------------
C *** G R A D
C ------------------------------------------------------------------
C
C *** THE GRAD FUNCTION SUBPROGRAM COMPUTES THE FOLLOWING DIRECTLY
C *** <P(J)!D + L(I)/R !P(I)> WITH L(I) > L(J)
C
REAL FUNCTION GRAD(I,J)
c IMPLICIT REAL*8(A-H,O-Z)
COMMON /PARAM/D0,D1,D2,D3,D4,D5,D6,D10,d15,d30,H,H1,Z,NO,ND
COMMON R(220),RR(220),R2(220),P(50,220),AZ(50),L(50),MAX(50)
LL = MAX0(L(I),L(J))
II = I
JJ = J
IF ( L(I) .GT. L(J) ) GO TO 1
II = J
JJ = I
1 A1 = (LL+D5)/(LL*(LL+1)*(2*LL+1))
GRAD = R(1)*P(I,1)*P(J,1)*(D1 + A1*Z*R(1))
DL = D5*P(I,1)*P(J,1)*R(1)
MM = MIN0(MAX(I)+1,MAX(J)+1,ND)
K = 2
F1 = D5*(P(II,K+1) - P(II,K-1))
F2 = P(II,K+1) - D2*P(II,K) + P(II,K-1)
G0 = P(JJ,K)*R(K)
G1 = D5*(P(JJ,K+1)*R(K+1) - P(JJ,K-1)*R(K-1))
G2 = P(JJ,K+1)*R(K+1) - D2*P(JJ,K)*R(K) + P(JJ,K-1)*R(K-1)
GRAD = GRAD + D2*F1*G0 +(D2*F2*G1 + F1*G2)/D3
DL = DL + D2*P(II,K)*P(JJ,K)*R(K) + P(II,K+1)*P(JJ,K+1)*R(K+1)
DO 2 K = 4,MM,2
F1 = D5*(P(II,K+1) - P(II,K-1))
F2 = P(II,K+1) - D2*P(II,K) + P(II,K-1)
F3 = D5*(P(II,K+2) - P(II,K-2)) - D2*F1
F4 = P(II,K+2) + P(II,K-2) - D4*(P(II,K+1) + P(II,K-1))
1 + D6*P(II,K)
G0 = P(JJ,K)*R(K)
G1 = D5*(P(JJ,K+1)*R(K+1) - P(JJ,K-1)*R(K-1))
G2 = P(JJ,K+1)*R(K+1) - D2*P(JJ,K)*R(K) + P(JJ,K-1)*R(K-1)
G3 = D5*(P(JJ,K+2)*R(K+2) - P(JJ,K-2)*R(K-2)) -D2*G1
G4 = P(JJ,K+2)*R(K+2) + P(JJ,K-2)*R(K-2) - D4*(P(JJ,K+1)*R(K+1)
1 + P(JJ,K-1)*R(K-1)) + D6*P(JJ,K)*R(K)
GRAD = GRAD + D2*F1*G0 +(D2*F2*G1 + F1*G2)/D3
1 - (F1*G4-F4*G1 + D4*(F2*G3-F3*G2))/90.D0
2 DL = DL + D2*P(II,K)*P(JJ,K)*R(K) + P(II,K+1)*P(JJ,K+1)*R(K+1)
GRAD = GRAD + (LL+D5)*DL*H1
IF (II .EQ. I) GRAD = - GRAD
RETURN
END
C
C ------------------------------------------------------------------
C *** GRAD2
C ------------------------------------------------------------------
C
C *** THE GRAD2 FUNCTION SUBPROGRAM COMPUTES THE FOLLOWING DIRECTLY
C *** <P(I)!R.D + F(DELTA)! P(J)>
C
REAL FUNCTION GRAD2(I,J)
c IMPLICIT REAL*8(A-H,O-Z)
COMMON/PARAM/D0,D1,D2,D3,D4,D5,D6,D10,D15,D30,H,H1,Z,NO,ND
COMMON R(220),RR(220),R2(220),P(50,220),AZ(50),L(50),MAX(50)
DIMENSION Q(220)
JJ=I
II=J
DO1K=1,NO
1 Q(K)=P(JJ,K)*R(K)
LI=L(I)
LJ=L(J)
A1=(LI+LJ+D2)/((LI+D1)*(LJ+D1))
A2=((LJ+D5)*(LJ+D1)+(LJ+D1+D5)*(LI+D1))/((LI+D1)*(LJ+D1))
A=A1-A2*(LI+LJ+D3)/((LI+LJ+D4)*(LJ+D5))
FACT=(LJ+D5)/(LI+LJ+D3)
G=R(1)**2*P(I,1)*P(J,1)*FACT*(1.+A*Z*R(1))
MM=MIN0(MAX(I)+1,MAX(J)+1,ND)
K=2
F1=D5*(P(II,K+1)-P(II,K-1))
F2=P(II,K+1)-D2*P(II,K)+P(II,K-1)
G0=Q(K)*R(K)
G1=D5*(Q(K+1)*R(K+1)-Q(K-1)*R(K-1))
G2=Q(K+1)*R(K+1)-D2*Q(K)*R(K)+Q(K-1)*R(K-1)
G=G+D2*F1*G0+(D2*F2*G1+F1*G2)/D3
DO2K=4,MM,2
F1=D5*(P(II,K+1)-P(II,K-1))
F2=P(II,K+1)-D2*P(II,K)+P(II,K-1)
F3=D5*(P(II,K+2)-P(II,K-2))-D2*F1
F4=P(II,K+2)+P(II,K-2)-D4*(P(II,K+1)+P(II,K-1))
1 +D6*P(II,K)
G0=Q(K)*R(K)
G1=D5*(Q(K+1)*R(K+1)-Q(K-1)*R(K-1))
G2=Q(K+1)*R(K+1)-D2*Q(K)*R(K)+Q(K-1)*R(K-1)
G3=D5*(Q(K+2)*R(K+2)-Q(K-2)*R(K-2))-D2*G1
G4=Q(K+2)*R(K+2)+Q(K-2)*R(K-2)-D4*(Q(K+1)*R(K+1)
1 +Q(K-1)*R(K-1))+D6*Q(K)*R(K)
G=G+D2*F1*G0+(D2*F2*G1+F1*G2)/D3
1 -(F1*G4-F4*G1+D4*(F2*G3-F3*G2))/90.D0
2 CONTINUE
U=QUADR(JJ,II,0)
G=G-D5*U
DELTA=LJ-LI
IF(DELTA)10,11,12
10 GRAD2=G-(LI-D2)*U
RETURN
11 GRAD2=G+(D1+D5)*U
RETURN
12 GRAD2=G+(LI+D3)*U
RETURN
END
▶EOF◀