|
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: 111360 (0x1b300)
Types: TextFile
Names: »tmchf77«
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
└─⟦4334b4c0b⟧
└─⟦this⟧ »tmchf77«
program mchf77
c version 1981-07-14
LONG ASTER,NEXT
LOGICAL FAIL,RSCAN(20),OMIT,PRINT,OUT,ORTHO,EZERO,LD,ALL
REAL EL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
zone oucz(128,1,stderror)
INTEGER OUC,OUD,OUF,OUH
COMMON /INOUT/ IUC,IUD,IUF,IUH,OUC,OUD,OUF,OUH
data z/0.d0/,tol/0.d0/,no/0/,nd/0/,nwf/0/,np/0/,ncfg/0/,
1 ib/0/,ic/0/,nscf/0/
data pde/220*0.d0/,ek/20*0.d0/,e/400*0.d0/,ed/0.d0/,azd/0.d0/,
1 az/20*0.d0/,sum/20*0.d0/,s/20*0.d0/,dpm/20*0.d0/,
2 acc/20*0.d0/,meth/20*0/,ipr/0/
data iuc/0/,iud/0/,iuf/0/,iuh/0/,ouc/0/,oud/0/,ouf/0/,ouh/0/
C
ASTER='*'
READ (5,3) IUC,IUD,IUF,IUH,OUC,OUD,OUF,OUH
3 FORMAT(4I3)
if (ouc.eq.6) goto 525
call zassign(oucz,ouc)
call open(oucz,4,'hfwf',0)
C
C ***** INITIALIZE
C
525 CALL INIT
do 600 i=1,20
rscan(i)=.false.
600 continue
ortho=.false.
omit=.false.
ezero=.false.
ezero=.false.
all=.false.
H3 = H/D3
H1 = D2*H3
CH = H*H/D12
EH = EXP(-H)
c the variables in 13 must be initialized
acfg=0.d0
id=0
cfgtol=0.d0
scftol=0.d0
ld=.false.
c end of new cards
1 FAIL = .FALSE.
DO 4 I=1,20
DPM(I) = D10
DO 4 J=1,20
IF (I .NE. J) E(I,J) = 1.E-10
4 CONTINUE
C
C ***** DETERMINE DATA ABOUT THE PROBLEM
C
CALL DATAIO(NEW)
13 READ (5,2) PRINT,NSCF,IC,ACFG,ID,CFGTOL,SCFTOL,LD
2 FORMAT(L3,2I3,F3.1,I3,2F6.1,L3)
550 format(' test output',l3,2i8,f8.1,i8,2f8.1,l3)
a write(6,550) print,nscf,ic,acfg,id,cfgtol,scftol,ld
C
C ***** PERFORM THE MCHF ITERATION
C
CALL SCF(ETOTAL,ACFG,SCFTOL,CFGTOL,LD,NEW)
IF (FAIL) GO TO 6
C
C ***** OUTPUT RESULTS IF PRINT = .TRUE.
C
IF (IB .GT. NWF ) GO TO 15
CALL OUTPUT(PRINT)
CALL SUMMRY(ETOTAL)
C
C ***** CHECK FOR ISOELECTRONIC SEQUENCE OR END OF CASE.
C
15 READ (5,10) NEXT,ATOM,ZZ,(ACC(I),I=1,NWF)
10 FORMAT(1X,I2,A6,F6.0,20F3.1)
IF (NEXT .EQ. 0) GO TO 1
C
C ***** SCALE RESULTS FOR ANOTHER MEMBER OF THE ISOELECTRONIC SEQUENCE
C
CALL SCALE(ZZ)
WRITE (6,14) ATOM,TERM
14 FORMAT(1H1,9X,2A6)
CALL ORTHOG
GO TO 13
C
C ***** SEARCH FOR AN * CARD DENOTING AN END OF CASE
C
6 READ (5,12) NEXT
12 FORMAT(A1)
IF (READERR.SHIFT.(-8).EQ.25) GOTO 5
IF (NEXT .EQ. ASTER) GO TO 1
GO TO 6
5 STOP
END
C
C ------------------------------------------------------------------
C 3-2 A R R A Y
C ------------------------------------------------------------------
C K
C FROM THE AVERAGE ENERGIES AS WELL AS FROM CONTRIBUTION OF F
C K
C AND G WHICH ARISE EITHER AS DEVIATIONS FROM AN AVERAGE ENERGY OR
C AS AN INTERACTION, TWO ARRAYS A(I,J,K) AND B(I,J,K) ARE
C 2(K-1)
C GENERATED. THE VALUE OF A IS THE COEFFICIENT OF Y (J,J)
C IJK
C IN THE POTENTIAL FUNCTION FOR THE I'TH RADIAL FUNCTION,
C WHEREAS THE VALUE
C !L -L !+2(K-1)
C OF B IS THE COEFFICIENT OF Y I J (I,J)P (R)
C IJK J
C IN THE EXCHANGE FUNCTION.
C
C IN AN MCHF CALCULATION, THESE COEFFICIENTS DEPEND ON THE
C MIXING COEFFICIENTS.
C
C
SUBROUTINE ARRAY
REAL QC
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
data r/220*0.d0/,rr/220*0.d0/,r2/220*0.d0/,p/4400*0.d0/,
1 y/4400*0.d0/,yk/220*0.d0/,yr/220*0.d0/,
2 x/220*0.d0/,n/20*0/,l/20*0/,max/20*0/
C
C ***** INITIALIZE ARRAYS
C
DO 1 I=1,NWF
SUM(I) = D0
DO 1 J=1,NWF
DO 1 K=1,5
A(I,J,K) = D0
1 B(I,J,K) = D0
C
C ***** DETERMINE COEFFICIENTS FOR AVERAGE ENERGIES OF CONFIGURATIONS
C
DO 2 K =1,NCFG
C = WT(K)**2
DO 3 I=1,NWF
IF (QC(I,K) .EQ. D0) GO TO 3
LI = L(I)+1
IF (LI .GT. 5) GO TO 3
DO 4 J=1,NWF
IF (QC(J,K) .EQ. D0) GO TO 4
LJ = L(J)+1
IF (LJ .GT. 5) GO TO 4
IF (I .NE. J) GO TO 5
CC = C*QC(I,K)*(QC(I,K) - D1)
A(I,I,1) = A(I,I,1) + CC
IF (LI .EQ. 1) GO TO 4
DO 8 KL = 2,LI
8 A(I,I,KL) = A(I,I,KL) - CC*CA(LI-1,KL-1)
GO TO 4
5 CC = C*QC(I,K)*QC(J,K)
A(I,J,1) = A(I,J,1) + CC
DO 9 KL =1,5
9 B(I,J,KL) = B(I,J,KL) - CC*CB(LI,LJ,KL)
4 CONTINUE
3 SUM(I) = SUM(I) + C*QC(I,K)
2 CONTINUE
C
C ***** CORRECT COEFFICIENTS FOR DEVIATIONS FROM THE AVERAGE FOR A
C ***** SPECIFIC SL TERM AS DETERMINED BY THE INPUT DATA
C
IF (NF .EQ. 0) GO TO 10
C
C ***** CORRECTIONS DUE TO 'FK' INTEGRALS
C
DO 11 I=1,NF
K = KFG(I)/2 + 1
NI = NCI(I)
NJ = NCJ(I)
CON = CFG(I)*WT(NI)*WT(NJ)
II = IFG(I)
IJ = JFG(I)
A(II,IJ,K) = A(II,IJ,K) + CON
11 A(IJ,II,K) = A(IJ,II,K) + CON
10 IF (NG .EQ. 0) GO TO 12
C
C ***** CORRECTIONS DUE TO 'GK' INTEGRALS
C
NFG = NF + NG
NFP = NF + 1
DO 13 I= NFP,NFG
II = IFG(I)
IJ = JFG(I)
NI = NCI(I)
NJ = NCJ(I)
LI = L(II)
LJ = L(IJ)
KL = (KFG(I) - IABS(LI - LJ))/2 + 1
CON = CFG(I)*WT(NI)*WT(NJ)
B(II,IJ,KL) = B(II,IJ,KL) + CON
13 B(IJ,II,KL) = B(IJ,II,KL) + CON
C
C ***** DIVIDE ALL ARRAYS FOR FUNCTION I BY THE EXPECTED OCCUPATION
C ***** NUMBER FOR FUNCTION I, IF DIFFERENT FROM ZERO
C
12 DO 15 I = 1,NWF
IF (SUM(I) .EQ. D0) GO TO 15
DO 14 J=1,NWF
DO 14 K=1,5
A(I,J,K) = A(I,J,K)/SUM(I)
14 B(I,J,K) = B(I,J,K)/SUM(I)
15 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-3 D A T A
C ------------------------------------------------------------------
C
C DATA CONCERNING THE NUMBER OF CONFIGURATIONS (NCFG), THE NUMBER
C AND TYPE OF ELECTRONS IN EACH CONFIGURATION, AS WELL AS DATA
C ASSOCIATED WITH THE ENERGY EXPRESSION ARE READ AND STORED.
C
SUBROUTINE DATAIO(NEW)
REAL QC
LONG BL,EL,EL1,EL2
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,LWT,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
COMMON /TEMP/ATM(20),TRM(20),ZZ(20),IND(20),IDUM(3181)
INTEGER OUC,OUD,OUF,OUH
COMMON /INOUT/ IUC,IUD,IUF,IUH,OUC,OUD,OUF,OUH
data atm/20*0.d0/,trm/20*0.d0/,zz/20*0.d0/,ind/20*0/,idum/3181*0/
C
C ***** READ 'ATOM' CARD
C
BL=' '
READ (5,1) ATOM,TERM,Z,NO,NWF,NIT,NCFG,NF,NG,NR,NL,ORTHO,
1 OMIT,NEW
1 FORMAT(2A6,F6.0,I6,7I3,2L3,I3)
IB = NWF - NIT + 1
ND = NO - 2
WRITE (6,2) ATOM,TERM,Z
2 FORMAT(1H1///9X,33HHARTREE-FOCK WAVE FUNCTIONS FOR ,2A6,4H Z =,
1 F5.1//27X,13HCONFIGURATION,15X,6HWEIGHT//)
OMIT = .NOT. OMIT
C
C ***** READ 'CONFIGURATION' CARDS AND NORMALIZE THE WEIGHTS
C
W = D0
DO 3 I=1,NCFG
READ (IUC,4) CONFIG(1,I),CONFIG(2,I),CONFIG(3,I),WTT,LWT
IF (.NOT. LWT) WT(I) = WTT
4 FORMAT(3A8,F10.8,L1)
3 W = W + WT(I)**2
C
C ***** IF INPUT WEIGHTS ARE ALL ZERO, SET EACH = 1/SQRT(NCFG)
C
IF (W .NE. D0) GO TO 21
W = NCFG
DO 22 I=1,NCFG
22 WT(I) = D1
21 W = SQRT(W)
DO 18 I = 1,NCFG
WT(I) = WT(I)/W
18 WRITE (6,5) I, CONFIG(1,I), CONFIG(2,I), CONFIG(3,I), WT(I)
5 FORMAT(14X,I2,6X,3A8,F19.8)
WRITE (6,6)(I,I=1,NCFG)
6 FORMAT(//9X,10HINPUT DATA/9X,10H----- ----//13X,13HWAVE FUNCTION,
1 11H PROCEDURE,6X,41HNUMBER OF ELECTRONS IN EACH CONFIGURATION/
2 17X,22HNL SIGMA METH ACC OPT,20I4/40X,20I4)
WRITE (6,66)
66 FORMAT(//)
C
C ***** READ 'ELECTRON ' CARD
C
DO 9 I = 1,NWF
READ (5,7) EL(I),N(I),L(I),S(I),METH(I),ACC(I),IND(I),
1 (QC(I,J),J=1,NCFG)
7 FORMAT(A3,2I3,F6.2,I3 ,F3.1,I3,15F3.0/24X,15F3.0/24X,15F3.0)
WRITE (6,8) I,EL(I),N(I),L(I),S(I),METH(I),ACC(I),IND(I)
1 ,(QC(I,J),J=1,NCFG)
8 FORMAT(I8, 2X,A3,2I3,F7.1,I4,F4.1,I4,1X,20F4.0/40X,20F4.0)
SUM(I) = D0
DO 36 J = 1,NCFG
IF (QC(I,J) .GT. D0 ) SUM(I)=SUM(I) + QC(I,J)*WT(J)**2
IF ( QC(I,J) .LT. D0 ) QC(I,J) = 1.E-15
36 CONTINUE
C
C ***** INITIALIZE ARRAYS, IF NECESSARY
C
IF (IND(I) .EQ. 1) GO TO 9
E(I,I) = D0
EK(I) = D0
AZ(I) = D0
9 RSCAN(I) = .TRUE.
WRITE (6,10)
10 FORMAT(//9X,9HENERGY = ,9X,12HE(AVERAGE) + //)
C
C ***** READ 'FK' AND 'GK' CARDS
C
NFG = NF + NG
IF (NFG .EQ. 0) GO TO 13
DO 20 I=1,NFG
READ (IUD,11) CFG(I),W,KFG(I),IFG(I),NCI(I),JFG(I),NCJ(I)
11 FORMAT(F12.8,A1,I1,1X,2I2,1X,2I2)
IF ( OUD .EQ. 0 ) GO TO 20
I1 = IFG(I)
I2 = JFG(I)
WRITE (6,12) CFG(I),NCI(I),NCJ(I),W,KFG(I),EL(I1),EL(I2)
12 FORMAT(26X,F14.8,3H*C(,I2,4H)*C(,I2,2H) ,A1,I1,1H(,A3,1H,,A3,1H))
20 CONTINUE
13 IF (NR .EQ. 0) GO TO 17
C
C ***** READ 'RK' CARDS AND SET RSCAN(I)=.FALSE. IF WAVEFUNCTION I
C ***** IS PRESENT IN THE RK INTEGRAL
C
DO 14 I=1,NR
READ (IUD,15) CR(I),KR(I),I1R(I),I2R(I),NCRI(I),J1R(I),J2R(I),
1 NCRJ(I),IO(I),JO(I),IQ(I)
15 FORMAT(F12.8,1X,I1,1X,3I2,1X,3I2,2X,I2,1X,I2,1X,I2)
I1 = I1R(I)
I2 = I2R(I)
J1 = J1R(I)
J2 = J2R(I)
EL1 = BL
EL2 = BL
IF (IQ(I) .EQ. 0) GO TO 30
II = IO(I)
JJ = JO(I)
RSCAN(II) = .FALSE.
RSCAN(JJ) = .FALSE.
EL1 = EL(II)
EL2 = EL(JJ)
30 IF ( OUD .EQ. 0 ) GO TO 23
WRITE (6,16) CR(I),NCRI(I),NCRJ(I),KR(I),EL(I1),EL(I2),EL(J1),
1 EL(J2),EL1,EL2,IQ(I)
16 FORMAT(26X,F14.8,3H*C(,I2,4H)*C(,I2,3H) R,I1,1H(,2A3,1H,,2A3,2H)<
1 ,A3,1H!,A3,1H>,I2)
23 RSCAN(I1) = .FALSE.
RSCAN(I2) = .FALSE.
RSCAN(J1) = .FALSE.
14 RSCAN(J2) = .FALSE.
17 IF (NL .EQ. 0) GO TO 31
C
C ***** READ 'L' CARDS AND SET RSCAN(I)=.FALSE. IF WAVEFUNCTION I
C ***** IS PRESENT IN THE L INTEGRAL
C
DO 32 I = 1,NL
READ (IUD,33) CL(I),ILI(I),NCLI(I),ILJ(I),NCLJ(I),ILO(I),JLO(I),
1 LQ(I)
33 FORMAT(F12.8,2X,2I2,1X,2I2,2X,I2,1X,I2,1X,I2)
I1 = ILI(I)
J1 = ILJ(I)
EL1 = BL
EL2 = BL
IF (LQ(I) .EQ. 0) GO TO 34
I2 = ILO(I)
J2 = JLO(I)
EL1 = EL(I2)
EL2 = EL(J2)
RSCAN(I2) = .FALSE.
RSCAN(J2) = .FALSE.
34 IF ( OUD .EQ. 0 ) GO TO 24
WRITE (6,35) CL(I),NCLI(I),NCLJ(I),EL(I1),EL(J1),EL1,EL2,LQ(I)
35 FORMAT(26X,F14.8,3H*C(,I2,4H)*C(,I2,4H) L(,A3,1H,,A3,2H)<,
1 A3,1H!,A3,1H>,I2)
24 RSCAN(I1) = .FALSE.
32 RSCAN(J2) = .FALSE.
31 IF ( OUD .EQ. 0 ) WRITE(6,25) NF, NG, NR, NL
25 FORMAT(19X,I3,17H FK INTEGRAL(S) +/19X,I3,17H GK INTEGRAL(S) +/
1 19X,I3,17H RK INTEGRAL(S) +/19X,I3,15H L INTEGRAL(S) )
C
C ***** COMPUTE THE INITIAL ARRAY AND INITIAL RADIAL FUNCTIONS
C
IF ( NIT .NE. 0) CALL ARRAY
CALL WAVEFN
RETURN
STOP
END
C
C ------------------------------------------------------------------
C 3-4 D E
C ------------------------------------------------------------------
C
C THIS ROUTINE CONTROLS THE SOLUTION OF THE DIFFERENTTIAL EQUATION
C FOR THE RADIAL FUNCTION P . ONE OF THREE METHODS IS SELECTED -
C I1
C M1, M2, OR M3 - FOR SOLVING THE EQUATIONS, THE INITIAL CHOICE
C BEING DETERMINED BY AN INPUT PARAMTER, METH(I), EXCEPT WHEN NO
C EXCHANGE IS PRESENT, IN WHICH CASE M2 IS SELECTED. (FOR FURTHER
C INFORMATION SEE SEC. 7-4)
C
C VALUE OF METH(I) METHOD
C --------------- ------
C < OR =1 M1 WITH SEARCH FOR AN ACCEPTABLE SOLUTION
C =2 M2 WITH SEARCH FOR AN ACCEPTABLE SOLUTION
C =3 M1 WITHOUT ANY CHECKING OF THE SOLUTION
C
C IF M1 FAILS TO FIND AN ACCEPTABLE SOLUTION, THE RADIAL FUNCTIONS
C ARE ORTHOGONALIZED, OFF-DIAGONAL ENERGY PARAMETERS RECOMPUTED,
C AND THE METHOD TRIED AGAIN. SHOULD IT CONTINUE TO FAIL, METH(I)
C IS SET TO 2.
C
SUBROUTINE DE(I1)
LONG EL
LONG ASTER(3)
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,OTHER,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO ,ALL
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /TEMP/P2(220),HQ(220),XX(220),AZZ,PP,FN,EM,FM,EU,FU,
1 DELTAE,M,NODE,MK,KK,NJ,IDUM(1980)
ASTER(1)=' '
ASTER(2)='* '
ASTER(3)='**'
I = I1
ED2 = E(I,I)
KK= MAX0(1,METH(I))
IF (NWF .EQ. 1) KK = 2
NODE = N(I) - L(I) - 1
C
C ***** CALL METHOD TO SOLVE THE DIFFERENTIAL EQUATION
C
CALL METHOD(I)
IF ( FAIL ) GO TO 25
C
12 PN = SQRT(QUAD(I,M,PDE,PDE))
DO 9 J = 1,M
9 PDE(J) = PDE(J)/PN
AZZ = AZZ/PN
C
C ***** SET THE ACCELERATING PARAMETER
C
C
IF (IPR .NE. I ) GO TO 33
ED2 = ED2 - E(I,I)
IF (ED1*ED2 .GT. D0) ACC(I) = .75*ACC(I)
IF (ED1*ED2 .LT. D0) ACC(I) = (D1 + D3*ACC(I))/D4
33 C = ACC(I)
CD = D1 - C
C
C ***** IMPROVE THE ESTIMATES
C
MAX(I) = M
DP = D0
DO 21 J = 1,M
DIFF = P(I,J)-PDE(J)
DP = AMAX1(DP ,ABS(DIFF)*R2(J))
21 P(I,J) = PDE(J) + C*DIFF
IF (M .EQ. NO) GO TO 26
M = M + 1
DO 24 J = M,NO
24 P(I,J) = D0
AZ(I) = CD*AZZ + C*AZ(I)
IF (C .EQ. D0) GO TO 28
26 PNN = SQRT(QUADR(I,I,0))
DO 30 J = 1,M
30 P(I,J) = P(I,J)/PNN
AZ(I) = AZ(I)/PNN
C
C ***** CHECK THE ORTHOGONALIZATION
C
28 DO 60 J = 1,NWF
IF (J .GE. IB .AND. OMIT ) GO TO 60
IF (N(I) .EQ. N(J) .OR. L(I) .NE. L(J)) GO TO 60
IF (A(I,J,1) .EQ. D0 .AND. .NOT. ORTHO) GO TO 60
J1 = J
J2 = I
IF ( J .LT. IB .OR. DPM(J) .LT. DPM(I)) GO TO 61
J1 = I
J2 = J
61 C = QUADR(J1,J2,0)
D = SQRT(D1 - C*C)
WRITE (6,63) EL(J1),EL(J2),C
63 FORMAT(6X,1H(,A3,1H/,A3,2H)=,1PE8.1)
IF (P(J2,1) - C*P(J1,1) .LT. D0) D = -D
DO 62 JJ = 1,NO
62 P(J2,JJ) = (P(J2,JJ) - C*P(J1,JJ))/D
AZZ = (AZ(J2) - C*AZ(J1))/D
IF (AZZ .GT. D0) AZ(J2) = AMAX1(AZZ,D5*AZ(J2))
IF (J2 .EQ. I) GO TO 60
CALL YKF(J2,J2,0)
DO 64 JJ = 1,NO
64 Y(J2,JJ) = YK(JJ)
60 CONTINUE
C
C ***** GENERATE THE IMPROVED Y0 ARRAY
C
CALL YKF(I,I,0)
DO 31 J = 1,NO
31 Y(I,J) = YK(J)
WRITE (6,17) EL(I),E(I,I),AZ(I),PN,ASTER(KK),DP
17 FORMAT(20X,A3,2F15.7,F12.7, A2,1PE10.2)
C
DPM(I) = DP
IF (IPR .EQ. I1) ED1 = ED2
IF (IPR .NE. I1) ED1 = ED2 - E(I1,I1)
IPR = I1
RETURN
C
C ***** IF METHOD FAILED TO FIND AN ACCEPTABLE SOLUTION, ORTHOGONALIZE
C ***** THE ESTIMATES AND TRY AGAIN
C
25 CALL ORTHOG
CALL GRANGE
27 CALL METHOD(I)
IF ( FAIL ) GO TO 23
GO TO 12
C
C ***** ERROR RETURN FROM SECOND TRY. IF M1 WAS USED,SWITCH TO
C M2 AND TRY ONCE MORE.
C
23 IF ( KK .EQ. 2) RETURN
KK = 2
GO TO 27
END
C
C ------------------------------------------------------------------
C 3-5 D I A G
C ------------------------------------------------------------------
C
C THE DIAG SUBROUTINE COMPUTES AN ENERGY MATRIX AND, IF THE NUMBER
C OF CONFIGURATIONS IS GREATER THAN 1, FINDS AN EIGENVALUE AND
C EIGENVECTOR OF THIS MATRIX.
C
C GIVEN THE ACCELERATING PARAMETER FOR CONFIGURATION MIXING,
C ACFG, AND THE RELATIVE CONVERGENCE CRITERION, CFGTOL, THE
C ROUTINE RETURNS THE NEW VALUE OF THE TOTAL ENERGY AND THE
C VALUE OF A LOGICAL VARIABLE, ECONV, WHICH WILL BE .TRUE. IF
C THE RELATIVE DIFFERENCE OF THE TWO ENERGIES IS LESS
C THAN CFGTOL. THE NEW MIXING COEFFICIENTS ARE STORED IN COMMON.
C
C
SUBROUTINE DIAG(ETOTAL,ECONV,ACFG,CFGTOL,NC,LAST)
REAL QC
LONG EL
LOGICAL ECONV,LAST
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /TEMP/W(40,40),GKIJ(5),FKIJ,C,CC,CCC,RATIO,ETL,
1WP(40),CONT,ETPREV,LIJ,LI,LJ,IM,N1,N2,NFP,NFG,NM,IP,JP,II,JJ,I,J
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /HMATRX/ET(40,40)
INTEGER OUC,OUD,OUF,OUH
COMMON /INOUT/ IUC,IUD,IUF,IUH,OUC,OUD,OUF,OUH
C
C ***** COMPUTE KINETIC ENERGY IF NECESSARY
C
c nc may be negative in the call. that is nonsense
if (nc.le.0) nc=1
IF (IB .GT. NWF) GO TO 51
DO 50 I = IB,NWF
50 EK(I) = -D5*HL(I,I)
C
C ***** SET UP THE ENERGY MATRIX FOR DIFFERENT CONFIGURATIONS LEAVING
C ***** AN NC BY NC PRINCIPAL SUBMATRIX UNCHAMGED
C
51 DO 52 I = 1,NCFG
52 WP(I) = WT(I)
DO 1 II = NC,NCFG
DO 1 JJ = 1,II
ET(JJ,II) = D0
1 ET(II,JJ) = D0
DO 2 I = 1,NWF
IF (SUM(I) .EQ. D0) GO TO 2
LI = L(I) + 1
IF (LI .GT. 5) GO TO 2
DO 3 J = 1,I
IF (SUM(J) .EQ. D0) GO TO 3
LJ = L(J) + 1
IF (LJ .GT. 5) GO TO 3
LIJ = IABS(LI - LJ)
KMAX = (LI + LJ - LIJ)/2
IF (I .NE. J) FKIJ = FK(I,J,0)
DO 4 K = 1,KMAX
4 GKIJ(K) = GK(I,J,LIJ+2*(K-1))
DO 5 II = NC,NCFG
C = QC(I,II)
IF (I .EQ. J) GO TO 6
CC = C*QC(J,II)
ET(II,II) = ET(II,II) + CC*FKIJ
DO 7 K = 1,KMAX
CCC = CB(LI,LJ,K)
7 ET(II,II) = ET(II,II) - CC*CCC*GKIJ(K)
GO TO 5
6 CC = D5*C*(C-D1)
ET(II,II) = ET(II,II) + C*EK(I) + CC*GKIJ(1)
IF (LI .EQ. 1) GO TO 5
DO 8 K = 2,KMAX
8 ET(II,II) = ET(II,II) - CA(LI-1,K-1)*CC*GKIJ(K)
5 CONTINUE
3 CONTINUE
2 CONTINUE
IF (NF .EQ. 0) GO TO 10
C
C ***** ADD CONTRIBUTIONS FROM 'FK' INTEGRALS
C
DO 11 I=1,NF
N1 = NCI(I)
N2 = NCJ(I)
CONT = D5*CFG(I)*FK(IFG(I),JFG(I),KFG(I))
ET(N1,N2) = ET(N1,N2) + CONT
11 ET(N2,N1) = ET(N2,N1) + CONT
10 IF (NG .EQ. 0) GO TO 12
C
C ***** ADD CONTRIBUTIONS FROM 'GK' INTEGRALS
C
NFP = NF+1
NFG = NF + NG
DO 13 I=NFP,NFG
N1 = NCI(I)
N2 = NCJ(I)
CONT = D5*CFG(I)*GK(IFG(I),JFG(I),KFG(I))
ET(N1,N2) = ET(N1,N2) + CONT
13 ET(N2,N1) = ET(N2,N1) + CONT
12 IF (NR .EQ. 0) GO TO 41
C
C ***** ADD CONTRIBUTIONS FROM 'RK' INTEGRALS
C
DO 15 I =1,NR
N1 = NCRI(I)
N2 = NCRJ(I)
CC = D1
IF (IQ(I) .NE. 0) CC = QUADR(IO(I),JO(I),0)**IQ(I)
CONT =D5*CR(I)*CC *RK(I1R(I),I2R(I),J1R(I),J2R(I),KR(I))
ET(N1,N2) = ET(N1,N2) + CONT
15 ET(N2,N1) = ET(N2,N1) + CONT
41 IF ( NL .EQ. 0 ) GO TO 14
C
C ***** ADD CONTRIBUTION FROM THE 'L' INTEGRALS
C
DO 42 I = 1,NL
N1 = NCLI(I)
N2 = NCLJ(I)
CC = D1
IF (LQ(I) .NE. 0) CC = QUADR(ILO(I),JLO(I),0)**LQ(I)
CONT = D5*CC*CL(I)*HL(ILI(I),ILJ(I))
ET(N1,N2) = ET(N1,N2) + CONT
42 ET(N2,N1) = ET(N2,N1) + CONT
14 IF (NCFG .EQ. 1) GO TO 37
IF (ID .GT. 0) GO TO 38
C
C ***** COMPUTE EIGENVALUE AND EIGENVECTOR BY THE METHOD OF
C ***** SEC. 6-8.2 THIS METHOD MAY CAUSE DIFFICULTIES WHEN NEAR
C ***** DEGENERACY EFFECTS ARE PRESENT SINCE IT MAY CONVERGE TO
C ***** THE WRONG EIGENVALUE. THE CODE UP TO, BUT NOT INCLUDING
C ***** STATEMENT NUMBER 31, MAY BE REPLACED BY A CALL TO A MORE
C ***** REFINED LIBRARY SUBROUTINE.
C
ETPREV = D0
DO 30 II=1,5
ETL = D0
C
C ***** DETERMINE ESTIMATE OF THE EIGENVALUE
C
DO 16 I=1,NCFG
CONT = D0
DO 17 J=1,NCFG
17 CONT = CONT + WT(J)*ET(I,J)
16 ETL = ETL + WT(I)*CONT
C
C ***** SOLVE SYSTEM OF EQUATIONS FOR EIGENVECTOR
C
DO 18 I=1,NCFG
DO 19 J=1,NCFG
19 W(I,J) = ET(I,J)
18 W(I,I) = W(I,I) - ETL
WT(1) = D1
IF (NCFG .NE. 2) GO TO 54
WT(2) = -W(2,1)/W(2,2)
GO TO 20
54 NM = NCFG -1
DO 21 I=2,NM
IP =I+1
DO 21 J=IP,NCFG
RATIO = W(J,I)/W(I,I)
DO 23 K=1,NCFG
23 W(J,K) = W(J,K) - RATIO*W(I,K)
21 CONTINUE
WT(NCFG) = -W(NCFG,1)/W(NCFG,NCFG)
DO 24 I = 2,NM
J = NCFG -I+1
JP = J+1
WT(J) = - W(J,1)
DO 25 K = JP,NCFG
25 WT(J) = WT(J) - W(J,K)*WT(K)
24 WT(J) = WT(J)/W(J,J)
20 DO 26 I=2,NCFG
26 WT(1) = WT(1) + WT(I)**2
WT(1) = D1/SQRT(WT(1))
DO 27 I=2,NCFG
27 WT(I) = WT(I)*WT(1)
C
C ***** ITERATE, IF NECESSARY, OTHERWISE OUTPUT RESULTS
C
IF (ABS((ETPREV-ETL)/ETL) .LT. 1.E-5) GO TO 31
30 ETPREV = ETL
WRITE (6,40)
40 FORMAT(///10X,47HMATRIX DIAGONALIZATION PROCEDURE NOT CONVERGING )
DELTAE = D0
GO TO 33
31 DELTAE = ETL -ETOTAL
33 ETOTAL = ETL
WRITE (6,32) ETOTAL
32 FORMAT(//10X,15HTOTAL ENERGY = ,F18.9 )
39 WRITE (6,34)
34 FORMAT(/20X,6HWEIGHT,13X,13HENERGY MATRIX )
CC = D0
DO 28 I = 1,NCFG
WT(I) = WT(I) + ACFG*(WP(I) - WT(I))
28 CC = CC + WT(I)**2
CC = D1/SQRT(CC)
DO 35 I=1,NCFG
35 WT(I) = WT(I)*CC
II = NC
IF (LAST) II = 1
DO 45 I = II,NCFG
45 WRITE (6,36) I,WT(I),(ET(I,J),J=1,I)
36 FORMAT(I10,F12.8,2X,7F15.7/(24X,7F15.7))
70 IF (.NOT. LAST .OR. OUC .EQ. 0 ) GO TO 49
C
C ***** PUNCH CONFIGURATIONS AND WEIGHTS ON UNIT OUC
C
WRITE(OUC,46) ATOM,TERM,ETOTAL
46 FORMAT(3X,2A6,F14.7)
DO 47 J = 1,NCFG
47 WRITE(OUC,48) CONFIG(1,J),CONFIG(2,J),CONFIG(3,J),WT(J)
48 FORMAT(3A8,F10.7)
IF ( OUH .EQ. 0 ) GO TO 49
DO 60 I = 1,NCFG
60 WRITE (OUH,61) (ET(I,J),J=1,I)
61 FORMAT(5F14.7)
49 ECONV = .FALSE.
IF (ABS(DELTAE/ETOTAL) .LE. CFGTOL) ECONV = .TRUE.
RETURN
37 ETOTAL = ET(1,1)
DELTAE = D0
WRITE (6,32) ETOTAL
GO TO 70
38 DELTAE = D0
ETOTAL = ET(1,1)
GO TO 39
END
C
C ------------------------------------------------------------------
C 3-6 D I F F
C ------------------------------------------------------------------
C
C
C STORES LP IN THE ARRAY YK. THE DIFFERENCE APPROXIMATION OF
C I
C EQ. (6-14) IS USED.
C
C
SUBROUTINE DIFF(I)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
C
C ***** FORM DD + 2Z/R -L(L+1)/RR!P(I)>
C
MM = MAX(I) - 3
FL = L(I)
TZ = Z + Z
C = (FL+D5)**2
HH = 180.E0*H*H
DO 11 K = 4,MM
11 YK(K) = (D2*(P(I,K+3)+P(I,K-3)) - 27.E0*(P(I,K+2)+P(I,K-2)) +
1 270.E0*(P(I,K+1)+P(I,K-1)) - 490.E0*P(I,K))/HH +
2 P(I,K)*(TZ*R(K) - C)
C
C ***** BECAUSE OF THE POSSIBILITY OF EXTENSIVE CANCELLATION NEAR THE
C ***** ORIGIN, SEARCH FOR THE POINT WHERE THE ASYMPTOTIC BEHAVIOUR
C ***** BEGINS AND SMOOTH THE ORIGIN.
C
LEXP = L(I) + 2
Y1 = YK(4)/R2(4)/R(4)**LEXP
Y2 = YK(5)/R2(5)/R(5)**LEXP
DO 1 K = 4,100
KP = K+2
Y3 = YK(KP)/R2(KP)/R(KP)**LEXP
IF (Y2 .EQ. D0) GO TO 1
IF (ABS(Y1/Y2 - D1) .LT. .05 .AND. ABS(Y3/Y2 - D1) .LT. .05)
1 GO TO 2
Y1 = Y2
Y2 = Y3
1 CONTINUE
WRITE (6,3) I
3 FORMAT(6X,47HASYMPTOTIC REGION NOT FOUND FOR FUNCTION NUMBER,I3)
STOP
C
C ***** ASYMPTOTIC REGION HAS BEEN FOUND
C
2 KP = K
KM = KP - 1
DO 4 K = 1,KM
4 YK(K) = Y1*R2(K)*R(K)**LEXP
MM = MM + 1
YK(MM) = (-(P(I,MM+2)+P(I,MM-2)) + D16*(P(I,MM+1)+P(I,MM-1))
1 -D30*P(I,MM))/(D12*H*H) + P(I,MM)*(TZ*R(MM) - C)
MM = MM + 1
YK(MM) = (P(I,MM+1) + P(I,MM-1) - D2*P(I,MM))/(H*H) +
1 P(I,MM)*(TZ*R(MM) - C)
MM = MM+1
DO 5 K =MM,NO
5 YK(K) = D0
RETURN
END
C
C ------------------------------------------------------------------
C 3-7 D Y K
C ------------------------------------------------------------------
C
C STORES IN YK THE VALUES OF THE INTEGRAL OF
C K
C P (S/R) (DP /DS - P /S) INTEGRATED OVER THE INTERVAL (0,R)
C I J J
C
C WHICH ENTER INTO THE SPIN-ORBIT CALCULATION.
C
C
SUBROUTINE DYK(I,J,K)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
DEN = L(I)+ L(J)+ 2 + K
FACT = L(J)
DO 1 JJ =1,2
1 YK(JJ) = FACT*P(I,JJ)*P(J,JJ)*R(JJ)/DEN
A = EH**K
AA = A*A
A = D4*A
F1 = FACT*P(I,1)*P(J,1)*R(1)
F2 = FACT*P(I,2)*P(J,2)*R(2)
DO 8 M =3,ND
F3 =D5*P(I,M)*((P(J,M+1)-P(J,M-1))/H - P(J,M))*R(M)
YK(M) = YK(M-2)*AA + H3*(F3+ A*F2 + AA*F1)
F1 = F2
8 F2 = F3
A = A*(EH)**3
AA = A*A/D16
C = 2*K+3
HH = C*H3
YK(NO)= YK(ND)
F1 = YK(NO)
F2 = F1
DO 9 MM = 3,NO
M = NO -MM+1
F3 =YK(M)
YK(M) = YK(M+2)*AA + HH*(F3 +A*F2 + AA*F1)
F1 = F2
9 F2 = F3
RETURN
END
C
C ------------------------------------------------------------------
C 3-8 E K I N
C ------------------------------------------------------------------
C
C RETURNS THE VALUE OF THE INTEGRAL OF
C
C (2/R)P (Y P + X )
C J I I I
C
C INTEGRATED WITH RESPECT TO R.
C
C
FUNCTION EKIN(I,II)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
CALL XCH(I,2)
CALL POTL(I)
DO 1 J=1,NO
YK(J) = YR(J)
1 YR(J) = P(II,J)
EKIN = D2*QUADS(I,II,1) + QUAD(II,NO ,YR,X)
RETURN
END
C
C ------------------------------------------------------------------
C 3-9 E L A G R
C ------------------------------------------------------------------
C
C THE ROUTINE FIRST DETERMINES WHETHER OR NOT THE OFF-DIAGONAL
C ENERGY PARAMETER CONNECTING THIS PAIR OF FUNCTIONS MAY BE SET TO
C ZERO; IF NOT, IT COMPUTES A LAGRANGE MULTIPLIER USING AN
C EQUATION LIKE EQ. (7-8) PROVIDED !Q - Q ! > 0.05. THE
C I J
C EQUATION IS MODIFIED TO INCLUDE CONTRIBUTIONS FROM THE
C INTERACTIONS BETWEEN CONFIGURATIONS.
C
C
FUNCTION ELAGR(I,J)
REAL QC
LONG EL
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /TEMP/QI,QJ,G,DI,DJ,DII,DJJ,DIJ,DJI,IDUM(3303)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
QI = SUM(I)
QJ = SUM(J)
EZERO = .TRUE.
G = D0
IF (E(I,J) .EQ. D0) GO TO 1
ALL = .TRUE.
CALL ROTATE(I,J)
C
C ***** IF EZERO IS 'TRUE', THEN THE OFF-DIAGONAL ENERGY PARAMETERS
C ***** SHOULD BE ZERO
C
IF (EZERO) GO TO 1
IF (.NOT. EZERO .AND. G .EQ. D0) GO TO 1
G = D0
IF (ABS(QI - QJ) .LT. 5.E-2) GO TO 1
C
C ***** COMPUTE THE OFF-DIAGONAL ENERGY PARAMETER WHEN BOTH
C ***** FUNCTIONS ARE IN THE SCF ITERATION AND QI NOT EQUAL TO QJ
C
DO 13 K = 1,5
C = A(I,I,K) - A(J,I,K) - B(J,I,K)
IF (ABS(C ) .GT. 1.E-10) G = G + C*RK(I,I,I,J,2*(K-1))
C = A(J,J,K) - A(I,J,K) - B(I,J,K)
IF (ABS(C ) .GT. 1.E-10) G = G - C*RK(J,J,J,I,2*(K-1))
13 CONTINUE
DO 14 M = 1,NWF
IF (M .EQ. I .OR. M.EQ. J) GO TO 14
DO 15 K = 1,5
C = A(I,M,K) - A(J,M,K)
IF (ABS(C) .GT. 1.E-10) G = G + C*RK(I,M,J,M,2*(K-1))
C = B(I,M,K) - B(J,M,K)
IF (ABS(C) .GT. 1.E-10) G = G + C*RK(I,J,M,M,IABS(L(I)-L(M))
1 + 2*(K-1))
15 CONTINUE
14 CONTINUE
C
C
C ***** ADD CONTRIBUTION FROM THE INTERACTION OF CONFIGURATIONS
C
ALL = .FALSE.
IF (NR .EQ. 0) GO TO 11
DO 20 II = 1,NR
N1 = NCRI(II)
N2 = NCRJ(II)
C = WT(N1)*WT(N2)*CR(II)
CI = C/SUM(I)
CJ = C/SUM(J)
CALL RPERT(I,J,II)
DD = CI*DI - CJ*DJ
G = G + D5*DD
20 CONTINUE
C
C
C ***** ADD CONTRIBUTION FROM THE INTERACTION OF CONFIGURATIONS - HL
C
11 IF (NL .EQ. 0) GO TO 1
DO 30 II = 1,NL
N1 = NCLI(II)
N2 = NCLJ(II)
C = WT(N1)*WT(N2)*CL(II)
CI = C/SUM(I)
CJ = C/SUM(J)
CALL LPERT(I,J,II)
DD = CI*DI - CJ*DJ
G = G + D5*DD
30 CONTINUE
1 ELAGR = G
RETURN
END
C
C ------------------------------------------------------------------
C 3-10 F K
C ------------------------------------------------------------------
C K
C RETURNS THE VALUE OF F (I,J)
C
C
FUNCTION FK(I,J,K)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
IF (K .NE. 0) GO TO 1
DO 3 JJ = 1,NO
3 YK(JJ) = Y(I,JJ)
GO TO 4
1 CALL YKF(I,I,K)
4 FK = QUADS(J,J,1)
RETURN
END
C
C ------------------------------------------------------------------
C 3-11 G K
C ------------------------------------------------------------------
C K
C RETURNS THE VALUE OF G (I,J).
C
C
FUNCTION GK(I,J,K)
CALL YKF(I,J,K)
GK = QUADS(I,J,1)
RETURN
END
C
C ------------------------------------------------------------------
C 3-12 G R A N G E
C ------------------------------------------------------------------
C
C CONTROLS THE CALCULATION OF OFF-DIAGONAL ENERGY PARAMETERS.
C IT SEARCHES FOR ALL PAIRS (I,J) WHICH ARE CONSTRAINED THROUGH AN
C ORTHOGONALITY REQUIREMENT. WHEN ONE OF THE PAIR , SAY P
C I
C MUST BE ORTHOGONAL NOT ONLY TO P BUT ALSO TO P WHERE N = N ,
C J K J K
C A SYSTEM OF EQUATIONS MUST BE SOLVED, THE EXACT FORM DEPENDING ON
C WHETHER OR NOT ANY OF THE FUNCTIONS ARE PART OF THE FROZEN CORE.
C WHEN ONLY ONE PAIR WITH A GIVEN SET OF PRINCIPAL QUANTUM NUMBERS
C IS PRESENT, ELAGR(I,J) IS USED TO DETERMINE THE OFF- DIAGONAL
C ENERGY PARAMETERS AS LONG AS !Q -Q ! > 0.05. OTHERWISE EQ.
C I J
C (7-10) MODIFIED FOR CONFIGURATION INTERACTION IS USED.
C
C
SUBROUTINE GRANGE
REAL EL
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,OTHER,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /TEMP/QI,QJ,D,DI,DJ,DL,G,FN,FL,RI(5),RJ(5),IL,M,MM,MX,IDUM(3279)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
DIMENSION KI(4),KJ(4),KT(4),AA(4,5),CC(4),XX(4)
EQUIVALENCE (I1,KI(1)),(I2,KI(2)),(J1,KJ(1)),(J2,KJ(2))
EQUIVALENCE (AA(1,5),CC(1))
IF (NWF .EQ. 1 .OR. IB .GT. NWF) RETURN
WRITE (6,26)
26 FORMAT(/)
II = MAX0(2,IB)
DO 2 I = II,NWF
C
C ***** WHEN ORTHO IS 'TRUE', CHECK IF ANOTHER ORBITAL WITH INDEX > I
C ***** HAS THE SAME NL VALUES
C
IF ( I .EQ. NWF .OR. .NOT. ORTHO ) GO TO 1
OTHER = .FALSE.
IP = I+1
DO 3 I2 = IP,NWF
IF (N(I2) .EQ. N(I) .AND. L(I2) .EQ. L(I) ) OTHER = .TRUE.
3 CONTINUE
IF (OTHER) GO TO 2
1 IL = I-1
DO 4 J = 1,IL
IF (N(I) .EQ. N(J) .OR. L(I) .NE. L(J)) GO TO 4
IF (A(I,J,1) .EQ. D0 .AND. .NOT. ORTHO) GO TO 4
KI(1) = I
KJ(1) = J
NI = 1
NJ = 1
IF (.NOT. ORTHO) GO TO 6
C
C ***** AN (I,J) PAIR CONSTRAINED BY AN ORTHOGONALITY REQUIREMENT HAS
C ***** BEEN FOUND. CHECK IF ANOTHER ORBITAL WITH INDEX > J HAS THE
C ***** SAME NL VALUES.
C
OTHER = .FALSE.
JP = J+1
DO 5 J2 = JP,NWF
IF (N(J2) .EQ. N(J) .AND. L(J2) .EQ. L(J)) OTHER = .TRUE.
5 CONTINUE
IF (OTHER) GO TO 4
C
C ***** THE (I,J) PAIR HAS THE HIGHEST POSSIBLE INDICES FOR THE GIVEN
C ***** SET OF NL VALUES
C
C ***** SEARCH FOR OTHER ORBITALS WITH THE SAME NL VALUES
C
DO 7 K = 1,NWF
IF (K .EQ. I .OR. N(K) .NE. N(I) .OR. L(K) .NE. L(I)) GO TO 8
NI = NI + 1
KI(NI) = K
GO TO 7
8 IF (K .EQ. J .OR. N(K) .NE. N(J) .OR. L(K) .NE. L(J)) GO TO 7
NJ = NJ + 1
KJ(NJ) = K
7 CONTINUE
IF ( NI*NJ .EQ. 1 ) GO TO 6
C
C ***** ONE OF THE NL VALUES HAS MORE THAN ONE ORBITAL ASSOCIATED WITH
C ***** IT
C
IF (NJ-2) 41,42,30
C
C ***** NOW NJ = 2
C
42 IF ( J1 .GE. IB .OR. J2 .GE. IB) GO TO 43
C
C ***** NOW NJ = 2 AND J1 AND J2 ARE BOTH PART OF THE FROZEN CORE
C
OV = QUADR(J1,J2,0)
DEN = D1 - OV**2
DO 31 K = 1,NI
KK = KI(K)
IF (KK .LT. IB ) GO TO 31
C1 = HL(KK,J1) - EKIN(KK,J1)
C2 = HL(KK,J2) - EKIN(KK,J2)
E(KK,J1) = (C1 - C2*OV)/DEN
E(KK,J2) = (C2 - C1*OV)/DEN
WRITE(6,24) EL(KK),EL(J1),E(KK,J1),EL(KK),EL(J2),E(KK,J2)
31 CONTINUE
GO TO 4
C
C ***** TOO MANY ORBITALS WITH THE SAME NL VALUES
C
30 WRITE(6,23) (KI(K),K=1,NI),(KJ(K),K=1,NJ)
23 FORMAT(/10X,50HTOO MANY ORBITALS WITH THE SAME NL VALUE..INDICES
1 ,3HARE/10X,20I3)
STOP
C
C ***** NJ = 2 BUT WITH EITHER J1 OR J2 NOT PART OF THE FROZEN CORE
C
43 IF (NI .GT. 2) GO TO 30
C
C ***** IF NECESSARY, INTERCHANGE THE TWO GROUPS SO THAT NI .LE. NJ
C
IF (NI .LE. NJ) GO TO 15
DO 10 K = 1,NI
10 KT(K) = KI(K)
NII = NI
DO 11 K = 1,NJ
11 KI(K) = KJ(K)
NI = NJ
DO 12 K = 1,NII
12 KJ(K) = KT(K)
NJ = NII
C
C ***** SET UP THE SYSTEM OF EQUATIONS FOR UP TO 4 OFF-DIAGONAL
C ***** ENERGY PARAMETERS
C
15 DO 48 K = 1,4
XX(K) = D0
DO 49 KK = 1,5
49 AA(K,KK) = D0
48 CONTINUE
IF (NI .EQ. 1) I2 = 0
K1 = 1
K2 = 2
K3 = 3
K4 = 4
IROW = 1
OVI = D0
IF (I2 .NE. 0) OVI = QUADR(I1,I2,0)
OVJ = QUADR(J1,J2,0)
DO 50 K = 1,2
DO 51 KK = 1,2
IF (I1 .EQ. 0) GO TO 55
IF (I1 .LT. IB) GO TO 52
C
C ***** NOW I1 .GE. IB
C
IF (J1 .LT. IB) GO TO 53
C
C
C ***** NOW I1 AND J1 .GE. IB
C
CC(IROW) = HL(I1,J1) - EKIN(I1,J1) + HL(J1,I1) - EKIN(J1,I1)
AA(IROW,K1) = D1 + SUM(I1)/SUM(J1)
IF (I2 .GT. 0) AA(IROW,K3) = SUM(I2)*OVI/SUM(J1)
GO TO 54
C
C ***** NOW I1 .GE. IB BUT J1 < IB
C
53 CC(IROW) = HL(I1,J1) - EKIN(I1,J1)
AA(IROW,K1) = D1
54 AA(IROW,K2) = OVJ
GO TO 55
C
C ***** I1 IS PART OF THE FROZEN CORE
C
52 IF (J1 .LT. IB) GO TO 55
C
C ***** I1 < IB BUT J1 .GE. IB
C
CC(IROW) = HL(J1,I1) - EKIN(J1,I1)
AA(IROW,K1) = SUM(I1)/SUM(J1)
IF (I2 .NE. 0) AA(IROW,K3) = SUM(I2)*OVI/SUM(J1)
C
C ***** NOW INTERCHANGE J1 AND J2 AND INCREMENT THE ROW
C
55 JT = J1
J1 = J2
J2 = JT
KTT = K1
K1 = K2
K2 = KTT
KTT = K3
K3 = K4
K4 = KTT
51 IROW = IROW + 1
C
C ***** NOW INTERCHANGE I1 AND I2 AND REPEAT
IT = I1
I1 = I2
I2 = IT
KTT = K1
K1 = K3
K3 = KTT
KTT = K2
K2 = K4
50 K4 = KTT
C
C ***** NOW SOLVE THE SYSTEM OF EQUATIONS
C
DO 60 K = 1,3
C
C ***** SEARCH FOR THE LARGEST PIVOT IN THE K'TH COLUMN
C
LM = K
DM = ABS(AA(K,K))
KP = K+1
DO 61 KK = KP,4
DD = ABS(AA(KK,K))
IF (DM .GE. DD) GO TO 61
DM = DD
LM = KK
61 CONTINUE
IF (LM .EQ. K) GO TO 65
C
C ***** INTERCHANGE THE ROWS
C
DO 62 KK = K,5
T = AA(K,KK)
AA(K,KK) = AA(LM,KK)
62 AA(LM,KK) = T
C
C ***** ELIMINATE THE K'TH VARIABLE
C
65 IF (AA(K,K) .EQ. D0) GO TO 60
DO 63 KK = KP,4
RATIO = AA(KK,K)/AA(K,K)
DO 64 KKK = KP,5
64 AA(KK,KKK) = AA(KK,KKK) - RATIO*AA(K,KKK)
63 CONTINUE
60 CONTINUE
C
C ***** BACKSUBSTITUTE ALLOWING FOR THE POSSIBILITY THAT THE RANK
C ***** OF THE SYSTEM MAY ONLY BE 2.
C
IF (AA(4,4) .NE. D0) XX(4) = CC(4)/AA(4,4)
IF (AA(3,3) .NE. D0) XX(3) = (CC(3) - AA(3,4)*XX(4))/AA(3,3)
XX(2) = (CC(2) - AA(2,3)*XX(3) - AA(2,4)*XX(4))/AA(2,2)
XX(1) = (CC(1) - AA(1,2)*XX(2) - AA(1,3)*XX(3) - AA( 1,4)*XX(4))/
1 AA(1,1)
C
C ***** COMPUTE AND PRINT THE OFF-DIAGONAL ENERGY PARAMETERS
C
IROW = 1
DO 70 K =1,2
IF (I1 .EQ. 0) GO TO 72
DO 71 KK = 1,2
E(I1,J1) = XX(IROW)
E(J1,I1) = SUM(I1)*XX(IROW)/SUM(J1)
WRITE(6,24) EL(I1),EL(J1),E(I1,J1),EL(J1),EL(I1),E(J1,I1)
24 FORMAT(10X,2HE(,2A3,4H) = ,F10.6,4X,2HE(,2A3,4H) = ,F10.6)
C
C ***** INTERCHANGE J1 AND J2
C
JT = J1
J1 = J2
J2 = JT
71 IROW = IROW + 1
C
C ***** INTERCHANGE I1 AND I2
C
72 IT = I1
I1 = I2
70 I2 = IT
GO TO 4
C
C ***** NJ = 1 BUT NI .GT. 1
C
41 IF (J1 .GE. IB) GO TO 43
C
C ***** NOW NJ = 1 AND J1 IS PART OF THE FROZEN CORE
C
DO 20 K = 1,NI
KK= KI(K)
E(KK,J1) = HL(KK,J1) - EKIN(KK,J1)
20 WRITE (6,25) EL(KK),EL(J1),E(KK,J1)
25 FORMAT(10X,2HE(,2A3,4H) = ,F10.6)
GO TO 4
C
C ***** NJ AND NI = 1
C
6 IF (J .GE. IB) GO TO 28
E(I,J) = HL(I,J) - EKIN(I,J)
GO TO 22
28 G = ELAGR(I,J)
IF ( G .NE. D0) GO TO 27
IF ( .NOT. EZERO ) GO TO 21
GO TO 4
27 G = D2*G/(QI - QJ)
E(I,J) = G*QJ
E(J,I) = G*QI
GO TO 22
21 E(I,J) = HL(I,J) - EKIN(I,J)
E(J,I) = HL(J,I) - EKIN(J,I)
E(I,J) = QJ*(E(I,J) + E(J,I))/(QI + QJ)
E(J,I) = QI*E(I,J)/QJ
22 WRITE(6,24) EL(I),EL(J),E(I,J),EL(J),EL(I),E(J,I)
4 CONTINUE
2 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-13 H L
C ------------------------------------------------------------------
C
C RETURNS THE VALUE OF <I!L!J>, USING A SPECIAL FORMULA TO
C PRESERVE SYMMETRY.
C
FUNCTION HL(I,J)
REAL EL
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
IF (IABS(L(I)-L(J)) .EQ. 0) GO TO 3
WRITE(6,4) EL(I),L(I),EL(J),L(J)
4 FORMAT(10X,44HUNALLOWED L VALUES OCCURRED IN HL SUBROUTINE/
1 2(10X,A3,9H HAS L = ,I3))
3 LI = L(I)
C = 2*LI + 1
A1 = -D2/(C*(LI+1))
A2 = A1/((C+D2)*(LI+1))
A3 = A2/((LI+2)*(LI+1))
ZR = Z*R(1)
HL = H*C*P(I,1)*P(J,1)*(D1+ZR*(A1+ZR*(A2+ZR*A3)))
MM = MIN0(MAX(I)+3,MAX(J)+3,ND-1)
K = 2
C = D4/D3
DI1 = P(I,K+1) - P(I,K-1)
DI2 = P(I,K+1) - D2*P(I,K) + P(I,K-1)
DJ1 = P(J,K+1) - P(J,K-1)
DJ2 = P(J,K+1) - D2*P(J,K) + P(J,K-1)
HL = HL + DI1*DJ1 + C*DI2*DJ2
DO 1 K = 4,MM,2
DI1 = P(I,K+1) - P(I,K-1)
DI2 = P(I,K+1) - D2*P(I,K) + P(I,K-1)
DI4 = P(I,K+2) - D4*(P(I,K+1)+P(I,K-1)) + D6*P(I,K) +P(I,K-2)
DI3 = P(I,K+2) - P(I,K-2) - D2*DI1
DI5 = P(I,K+3)-P(I,K-3) - D4*(P(I,K+2)-P(I,K-2))
1 + 5.E0*(P(I,K+1)-P(I,K-1))
DI6 = P(I,K+3)+P(I,K-3) - D6*(P(I,K+2)+P(I,K-2))
1 + 15.E0*(P(I,K+1)+P(I,K-1)) - 20.E0*P(I,K)
DJ1 = P(J,K+1) - P(J,K-1)
DJ2 = P(J,K+1) - D2*P(J,K) + P(J,K-1)
DJ4 = P(J,K+2) - D4*(P(J,K+1)+P(J,K-1)) + D6*P(J,K) +P(J,K-2)
DJ3 = P(J,K+2) - P(J,K-2) - D2*DJ1
DJ5 = P(J,K+3)-P(J,K-3) - D4*(P(J,K+2)-P(J,K-2))
1 + 5.E0*(P(J,K+1)-P(J,K-1))
DJ6 = P(J,K+3)+P(J,K-3) - D6*(P(J,K+2)+P(J,K-2))
1 + 15.E0*(P(J,K+1)+P(J,K-1)) - 20.E0*P(J,K)
1 HL = HL + DI1*DJ1 + C*DI2*DJ2 + (DI3*DJ3 + DI2*DJ4+DI4*DJ2)/45.E0
1 -(DI3*DJ5+DI5*DJ3)/252.E0 - (DI2*DJ6+DI6*DJ2-1.1*DI4*DJ4)/378.E0
TZ = Z + Z
C = (LI + D5)**2
HL2 = D5*(TZ*R(1) - C)*P(I,1)*P(J,1)
DO 2 K = 2,MM,2
2 HL2 = HL2 + D2*(TZ*R(K) - C)*P(I,K)*P(J,K)
1 + (TZ*R(K+1) - C)*P(I,K+1)*P(J,K+1)
HL = -HL/(D2*H) + HL2*H1
RETURN
END
C
C ------------------------------------------------------------------
C 3-14 H N O R M
C ------------------------------------------------------------------
C
C RETURNS THE VALUE OF THE NORMALIZATION CONSTANT FOR AN (NL)
C HYDROGENIC FUNCTION WITH NUCLEAR CHARGE ZZ.
C
C
FUNCTION HNORM(N,L,ZZ)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
M = L + L + 1
A = N + L
B = M
T = A
D = B
M = M - 1
IF (M .EQ. 0) GO TO 2
DO 1 I = 1,M
A = A - D1
B = B - D1
T = T*A
1 D = D*B
2 HNORM = SQRT(ZZ*T)/( N*D)
RETURN
END
C
C ------------------------------------------------------------------
C 3-15 H W F
C ------------------------------------------------------------------
C
C RETURNS THE VALUE OF AN UNNORMALIZED (NL) HYDROGENIC FUNCTION
C WITH NUCLEAR CHARGE ZZ AND RADIUS R.
C
C
FUNCTION HWF(N,L,ZZ,R)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
K = N-L-1
P = D1
A = D1
B = K
C = N+ L
X = -D2*ZZ*R/N
C
C ***** TEST IF UNDERFLOW MAY OCCUR, IF SO SET HWF = 0
C
IF ( X .LT. -180.E0 ) GO TO 5
IF (K) 1,2,3
3 DO 4 I = 1,K
P = D1 + A/B*P/C*X
A = A + D1
B = B - D1
4 C = C - D1
2 HWF = P*EXP(X/D2)*(-X)**(L+1)
RETURN
1 WRITE(6,7) N,L,ZZ,R
7 FORMAT(51H FORBIDDEN COMBINATION OF N AND L IN HWF SUBPROGRAM/
1 4H N = ,I4,6H L = ,I4,6H Z = ,F6.1,6H R = ,F8.4)
STOP
5 HWF = D0
RETURN
END
C
C ------------------------------------------------------------------
C 3-16 I N I T
C ------------------------------------------------------------------
C
C INITIALIZES BASIC CONSTANTS OF THE PROGRAM INCLUDING THOSE
C WHICH DEFINE THE AVERAGE ENERGY OF A CONFIGURATION.
C
C
SUBROUTINE INIT
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
C
C THE ARRAYS IN STATE COMMON AREA INITIALIZED
C A.L. 1981-07-13
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
DATA WT/40*0.E0/,QC/800*0.E0/,CFG/200*0.E0/,CR/300*0.E0/
DATA NCI/200*0.E0/,NCJ/200*0.E0/,NCRI/300*0.E0/,NCRJ/300*0.E0/
DATA CL/20*0.E0/,NCLI/20*0.E0/,NCLJ/20*0.E0/
DATA KFG/200*0/,NF/0/,NG/0/,NR/0/,NL/0/
DATA IFG/200*0/,JFG/200*0/,KR/300*0/,I1R/300*0/,I2R/300*0/
DATA J1R/300*0/,J2R/300*0/,IO/300*0/,JO/300*0/,IQ/300*0/
DATA ILI/20*0/,ILJ/20*0/,ILO/20*0/,JLO/20*0/,LQ/20*0/
C ***** SET THE COMMONLY USED CONSTANTS
C
D0 = 0.E0
D1 = 1.E0
D2 = 2.E0
D3 = 3.E0
D4 = 4.E0
D5 = 1./2.E0
D6 = 6.E0
D8 = 8.E0
D10 = 10.E0
D12 = 12.E0
D16 = 16.E0
D30 = 30.E0
C
C ***** SET THE STARTING POINT AND THE STEP SIZE
C
RHO = -4.E0
H = 1./16.E0
C
C ***** CLEAR THE ARRAYS FOR THE AVERAGE INTERACTIONS
C
DO 1 I = 1,4
DO 1 J = 1,4
1 CA(I,J) = 0.E0
DO 2 I = 1,5
DO 2 J = 1,5
DO 3 K = 1,5
3 CB(I,J,K) = 0.E0
2 CONTINUE
C
C ***** AVERAGE INTERACTIONS FOR EQUIVALENT ELECTRONS
C
C ***** P - P
C
CA(1,1) = 2.E0/25.E0
C
C ***** D - D
C
CA(2,1) = 2.E0/63.E0
CA(2,2) = 2.E0/63.E0
C
C ***** F - F
C
CA(3,1) = 4.E0/ 195.E0
CA(3,2) = 2.E0/ 143.E0
CA(3,3) = 100.E0/5577.E0
C
C ***** G - G
C
CA(4,1) = 20.E0/ 1309.E0
CA(4,2) = 162.E0/ 17017.E0
CA(4,3) = 20.E0/ 2431.E0
CA(4,4) = 4410.E0/371943.E0
C
C
C ***** AVERAGE INTERACTIONS FOR NON-EQUIVALENT ELECTRONS
C
C ***** S - ( S, P, D, F, G )
C
CB(1,1,1) = 1.E0/ 2.E0
CB(2,1,1) = 1.E0/ 6.E0
CB(3,1,1) = 1.E0/10.E0
CB(4,1,1) = 1.E0/14.E0
CB(5,1,1) = 1.E0/18.E0
C
C ***** P - ( P, D, F, G )
C
CB(2,2,1) = 1.E0/ 6.E0
CB(2,2,2) = 1.E0/ 15.E0
CB(3,2,1) = 1.E0/ 15.E0
CB(3,2,2) = 3.E0/ 70.E0
CB(4,2,1) = 3.E0/ 70.E0
CB(4,2,2) = 2.E0/ 63.E0
CB(5,2,1) = 2.E0/ 63.E0
CB(5,2,2) = 5.E0/198.E0
C
C ***** D - ( D, F, G )
C
CB(3,3,1) = 1.E0/ 10.E0
CB(3,3,2) = 1.E0/ 35.E0
CB(3,3,3) = 1.E0/ 35.E0
CB(4,3,1) = 3.E0/ 70.E0
CB(4,3,2) = 2.E0/105.E0
CB(4,3,3) = 5.E0/231.E0
CB(5,3,1) = 1.E0/ 35.E0
CB(5,3,2) = 10.E0/693.E0
CB(5,3,3) = 5.E0/286.E0
C
C ***** F - ( F, G )
C
CB(4,4,1) = 1.E0/ 14.E0
CB(4,4,2) = 2.E0/ 105.E0
CB(4,4,3) = 1.E0/ 77.E0
CB(4,4,4) = 50.E0/3003.E0
CB(5,4,1) = 2.E0/ 63.E0
CB(5,4,2) = 1.E0/ 77.E0
CB(5,4,3) = 10.E0/1001.E0
CB(5,4,4) = 35.E0/2574.E0
C
C ***** G - ( G )
C
CB(5,5,1) = 1.E0/ 18.E0
CB(5,5,2) = 10.E0/ 693.E0
CB(5,5,3) = 9.E0/ 1001.E0
CB(5,5,4) = 10.E0/ 1287.E0
CB(5,5,5) = 245.E0/21879.E0
C
C ***** SYMMETRIZE THE ARRAY
C
DO 4 I = 2,5
IM = I - 1
DO 5 J = 1,IM
DO 6 K = 1,4
6 CB(J,I,K) = CB(I,J,K)
5 CONTINUE
4 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-17 L P E R T
C -----------------------------------------------------------------
C
C LPERT DETERMINES THE EFFECT OF A PERTURBATION IN THE FORM OF A
C ROTATION OF THE I'TH AND J'TH ORBITAL BOTH ON THE ENERGY AND ON THE
C STATIONARY CONDITION, DUE TO THE PRESENCE OF A GIVEN L INTEGRAL.
C
C
SUBROUTINE LPERT(IIN,JIN,M)
REAL EL,QC
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /TEMP/QI,QJ,G,DI,DJ,DII,DJJ,DIJ,DJI,IDUM(3303)
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
DIMENSION IND(6)
EQUIVALENCE (I1,IND(1)),(I2,IND(2)),(I3,IND(3)),(I4,IND(4)),
1 (I5,IND(5)),(I6,IND(6))
I1 = ILI(M)
I2 = ILJ(M)
J = IIN
I = JIN
IPW = LQ(M)
IF (IPW .EQ. 0) GO TO 1
I5 = ILO(M)
I6 = JLO(M)
OVLAP = QUADR(I5,I6,0)**IPW
IF ((I5-I)*(I6-I)*(I5-J)*(I6-J) .NE. 0) GO TO 1
EZERO = .FALSE.
RETURN
1 DO 10 KP = 1,2
DI = D0
DII = D0
DIJ = D0
DO 2 K = 1,2
IF (IND(K) .NE. I) GO TO 2
IND(K) = J
DI = DI + HL(I1,I2)
IF (.NOT. ALL) GO TO 3
DO 4 K2 = 1,2
IF (IND(K2) .NE. J) GO TO 5
IND(K2) = I
DIJ = DIJ + HL(I1,I2)
IND(K2) = J
5 IF (IND(K2) .NE. I) GO TO 4
IND(K2) = J
DII = DII + HL(I1,I2)
IND(K2) = I
4 CONTINUE
3 IND(K) = I
2 CONTINUE
IF (IPW .EQ. 0) GO TO 6
DI = DI*OVLAP
DII = DII*OVLAP
DIJ = DIJ*OVLAP
6 IF (KP .EQ. 2) RETURN
DJ = DI
DJJ = DII
DJI = DIJ
I = IIN
J = JIN
10 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-18 M E T H O D
C ------------------------------------------------------------------
C
C USES M1, M2, OR M3 TO SOLVE THE RADIAL EQUATION. IF THE INPUT
C DATA INDICATED METH(I) = 3, THEN THIS SOLUTION IS RETURNED TO
C DE. OTHERWISE, THE ROUTINE SEARCHES FOR AN ACCEPTABLE SOLUTION
C WHICH IS BOTH POSITIVE NEAR THE ORIGIN AND HAS THE REQUIRED
C NUMBER OF NODES.
C
C
SUBROUTINE METHOD(I)
REAL EL
LOGICAL ONE,V1,V2,FIRST
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,OTHER ,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /TEMP/P2(220),HQ(220),XX(220),AZZ,PP,FN,EM,FM,EU,FU,
1 DELTAE,M,NODE,MK,KK,NJ,IDUM(1980)
DIMENSION P1(220)
EQUIVALENCE (PDE(1),P1(1))
C
C ***** 'FIRST' MUST BE 'TRUE' THE FIRST TIME SOLVE IS CALLED FOR
C ***** POTENTIAL AND EXCHANGE TO BE COMPUTED
C ***** 'ONE' IS SET 'TRUE' AS SOON AS AN ACCEPTABLE SOLUTION HAS BEEN
C ***** FOUND
C ***** 'EU' IS THE UPPER BOUND OF THE ENERGY PARAMETER
C ***** 'EM' IS THE MINIMUM VALUE OF THE ENERGY PARAMETER
C
FIRST = .TRUE.
FAIL = .FALSE.
EM = D0
EU = ((Z - AMIN1(D5*S(I),D2*S(I)))/N(I))**2
FU = EU
MK = 0
17 CALL SOLVE(I,FIRST)
C
C ***** IF KK NOT EQUAL 1, OMIT THE NODE CHECKING
C
IF (KK .EQ. 3) GO TO 51
C
C ***** COUNT THE NUMBER OF NODES
C
NC = NODEC(M)
C
C ***** IF NODE COUNT IS OFF BY NO MORE THAN 1 AND DELTAE IS STILL
C ***** QUITE LARGE, APPLY THE DELTAE CORRECTION
C
IF (IABS(NC-NODE) .EQ. 1 .AND. ABS(DELTAE/ED) .GT. 0.02D0)
1 GO TO 46
C
C ***** BRANCH ACCORDING TO WHETHER THE NODE COUNT IS TOO SMALL,
C ***** JUST RIGHT, OR TOO LARGE
C
IF (NC - NODE ) 8,9,10
C
C ***** THE SOLUTION HAS THE CORRECT NUMBER OF NODES
C
9 V2 = ABS(DELTAE)/ED .LT. 1.E-5
IF (PDE(1) .LT. D0 .AND. .NOT. V2) GO TO 46
IF (PDE(1) .GT. D0) GO TO 51
DO 52 J = 1,NO
52 PDE(J) = - PDE(J)
PP = -D2 - PP
51 AZZ = AZD*(D1 + PP)
E(I,I) = ED
RETURN
C
C ***** THE SOLUTION HAS TOO FEW NODES
C
8 IF (PDE(1) .LE. D0) GO TO 11
DEL = D1 - ED/EU
EU = ED
IF ( DEL .LT. .05E0) FU = FU*((L(I)+1+NC)/FN)**2.5
IF (DEL .GE. .05E0) FU = ED*((L(I)+1+NC)/FN)**2.5
IF (FU .LT. EM) FU = D5*(EU + EM)
IF (ABS(FU - ED) .LT. 0.001E0) GO TO 27
ED = FU
GO TO 33
C
C ***** TRY A NEW VALUE OF ED WHICH MUST LIE WITHIN THE UPPER AND
C ***** LOWER BOUND
C
11 EDP = ED
ED = ED*((L(I)+1+NC)/FN)**2.5
IF (ED .GE. EU ) ED = D5*(EU + EDP)
IF (ED .LE. EM ) ED = D5*(EM + EDP)
33 MK = MK + 1
IF ( EU .LE. EM ) WRITE(6,30) EM,EU,ED
30 FORMAT(6X,48HWARNING: DIFFICULTY WITH NODE COUNTING PROCEDURE/
1 6X,42HLOWER BOUND ON ED GREATER THAN UPPER BOUND/
2 6X,5HEL = ,F10.6,7H EU = ,F10.6,7H ED = ,F10.6)
FIRST = .FALSE.
IF ( MK .GT. 3*N(I) .OR. EU-EM .LT. FN**(-3)) GO TO 27
GO TO 17
C
C ***** THE SOLUTION HAS TOO MANY NODES
C
10 IF (PDE(1) .LT. D0) GO TO 11
DEL = D1 - EM/ED
EM = ED
IF (DEL .LT. 0.05E0) FM = FM*((L(I)+1+NC)/FN)**2.5
IF (DEL .GE. 0.05E0) FM = ED*((L(I)+1+NC)/FN)**2.5
IF (FM .GT. EU) FM = D5*(EU + EM)
IF (ABS(FM - ED) .LT. 0.001E0) GO TO 27
ED = FM
GO TO 33
C
C ***** ADJUST ENERGY TO LIE BETWEEN UPPER AND LOWER BOUND
C
46 ED = ED - DELTAE
IF (ED .LE. EM .OR. ED .GE. EU) ED = ED + DELTAE + DELTAE
IF (ED .GT. EM .AND. ED .LT. EU) GO TO 33
ED = ED - DELTAE
DELTAE = D5*DELTAE
GO TO 46
C
C ***** METHOD WAS UNABLE TO FIND AN ACCEPTABLE SOLUTION
C
27 WRITE (6,28) KK,EL(I),NC,NJ,ED,EM,EU
28 FORMAT(10X,6HMETHOD,I2,38H UNABLE TO SOLVE EQUATION FOR ELECTRON,
1 A3/10X,5HNC = ,I3,3X,5HNJ = ,I3,3X,5HED = ,F10.6,3X,5HEL = ,
2 F10.6,3X,5HEU = ,F10.6)
FAIL = .TRUE.
RETURN
END
C
C
C ------------------------------------------------------------------
C 3-19 N M R V S
C ------------------------------------------------------------------
C
C GIVEN TWO STARTING VALUES, PDE(1) AND PDE(2), VALUES OF PDE(J),
C J=3,4,...,NJ+1 ARE OBTAINED BY OUTWARD INTEGRATION OF
C Y" = YR Y + F
C USING THE DISCRETIZATION OF EQ. (6-27 ) WITH THE DIFFERENCE
C CORRECTION. WITH PDE(NJ) GIVEN, THE TAIL PROCEDURE IS APPLIED TO
C PDE(J),J=NJ+1, NJ+2,...,MM, WHERE MM IS DETERMINED AUTOMATICALLY
C AND DELTA IS THE DIFFERENCE BETWEEN PDE(NJ+1) FOR OUTWARD AND
C INWARD INTEGRATION. (SEE EQ 6-32, 6-33, AND 6-37 FOR FURTHER
C DETAILS.)
C
C
SUBROUTINE NMRVS(NJ,DELTA,MM,PDE,F)
DIMENSION PDE(220),F(220),A(150),D(150)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
EQUIVALENCE (G,G3)
C
C ***** INTEGRATE OUTWARD TO NJ+1
C
Y1 = PDE(1)
Y2= PDE(2)
G1 = YR(1)
G2 = YR(2)
M = NJ + 1
DO 1 I = 3,M
G3 = YR(I)
Y3 = (Y2+Y2-Y1 + (D10*G2*Y2 + G1*Y1) + F(I-1)) / (D1 - G3)
PDE(I) = Y3
Y1 = Y2
Y2 = Y3
G1 = G2
1 G2 = G3
DELTA = Y3
C
C ***** APPLY THE TAIL PROCEDURE
C
K = 1
PDE(M) = -(D1 - G1)*Y1 + F(M)
A(1) = D1 - G
D(1) = -(D2 + D10*G)
22 RATIO = A(K)/D(K)
C
C ***** THE INTEGER 149 IN THE NEXT STATEMENT IS THE DIMENSION OF A
C ***** MINUS 1
C
IF (K .GE. 149 .OR. M .EQ. ND) GO TO 23
K = K +1
M = M+1
G = YR(M)
A(K) = D1 - G
D(K) = -(D2 + D10*G) - A(K)*RATIO
PDE(M) = -PDE(M-1)*RATIO + F(M)
IF (ABS(PDE(M))+ABS(PDE(M-1)) .GT. TOL .OR. K .LT. 9) GO TO 22
20 CON =SQRT(EH)*EXP(-SQRT(ABS(G/CH-.25)/RR(M))*(R(M+1)-R(M)))
PDE(M) = PDE(M)/(D(K) + CON*(D1- YR(M+1)))
J = M+1
DO 2 I= J,NO
2 PDE(I) = D0
DO 3 J = 2,K
I = M-J+1
II = K-J+1
3 PDE(I) = (PDE(I)-A(II+1)*PDE(I+1))/D(II)
C
C ***** SET DELTA = DIFFERENCE OF THE TWO SOLUTIONS AT NJ+1
C ***** MM = NUMBER OF POINTS IN THE RANGE OF THE SOLUTION
C
DELTA = DELTA - PDE(I)
MM = M
RETURN
23 WRITE (6,24)
24 FORMAT(6X,52HWARNING: FUNCTIONS TRUNCATED BY NMRVS IN TAIL REGION)
GO TO 20
END
C
C ------------------------------------------------------------------
C 3-20 N O D E C
C ------------------------------------------------------------------
C
C COUNTS THE NUMBER OF NODES OF THE FUNCTION PDE(J) IN THE RANGE
C J=40,...,M-10. THE NODE COUNTING PROCEURE TERMINATES AS SOON AS
C A LOCAL MAX (OR MIN) IS ENCOUNTERED AFTER THE MAIN MAXIMUM FOR
C WHICH !PDE(J)! < 0.75 MAX !PDE(J)!
C J
C
FUNCTION NODEC(M)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
NCC= 0
DM = 0.E0
MM = M - 10
DO 36 J = 40,MM
X1 = ABS(PDE(J)*R2(J))
IF (X1 .LT. .75E0*DM) GO TO 12
DM = AMAX1(DM,X1)
NODEC = NCC
12 IF (PDE(J)*PDE(J+1) .LE. 0.E0) NCC= NCC+ 1
IF (X1 .EQ. 0.E0) NCC = NCC - 1
36 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-21 O R T H O G
C ------------------------------------------------------------------
C
C THIS ROUTINE ORTHOGONALIZES THE SET OF RADIAL FUNCTIONS WHEN AN
C ORTHOGONALITY CONSTRAIN APPLIES. A GRAM-SCHMIDT TYPE OF PROCESS
C IS USED. WHEN MORE THAN ONE RADIAL FUNCTION WITH A GIVEN (NL) IS
C PRESENT, IT MAY BE NECESSARY TO SOLVE A 2X2 SYSTEM OF EQUATIONS.
C
C
SUBROUTINE ORTHOG
REAL EL
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,OTHER,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /TEMP/QI,QJ,D,DI,DJ,DL,G,FN,FL,RI(5),RJ(5),IL,M,MM,MX,K,
2IDUM(3278)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
DIMENSION KI(4),KJ(4),KT(4)
EQUIVALENCE (I1,KI(1)),(I2,KI(2)),(J1,KJ(1)),(J2,KJ(2))
IF (NWF .EQ. 1 .OR. IB .GT. NWF) RETURN
WRITE (6,26)
26 FORMAT(/)
II = MAX0(2,IB)
DO 2 I = II,NWF
IL = I-1
DO 3 J = 1,IL
IF (N(I) .EQ. N(J) .OR. L(I) .NE. L(J)) GO TO 3
IF (A(I,J,1) .EQ. D0 .AND. .NOT. ORTHO) GO TO 3
KJ(1) = J
NJ = 1
IF (.NOT. ORTHO) GO TO 9
C
C ***** AN (I,J) PAIR CONSTRAINED BY AN ORTHOGONALITY REQUIREMENT HAS
C ***** BEEN FOUND. CHECK IF ANOTHER ORBITAL WITH INDEX > J HAS THE
C ***** SAME NL VALUES.
C
IF ( J .EQ. IL ) GO TO 1
OTHER = .FALSE.
JP = J+1
DO 5 J2 = JP,IL
IF (N(J2) .EQ. N(J) .AND. L(J2) .EQ. L(J)) OTHER = .TRUE.
5 CONTINUE
IF (OTHER) GO TO 3
C
C ***** THE (I,J) PAIR HAS THE HIGHEST POSSIBLE INDICES FOR THE GIVEN
C ***** SET OF NL VALUES
C ***** SEARCH FOR OTHER ORBITALS WITH THE SAME NL VALUES
C
1 IF ( J .EQ. 1 ) GO TO 8
JM = J - 1
C
DO 7 K = 1,JM
IF (N(K) .NE. N(J) .OR. L(K) .NE. L(J)) GO TO 7
NJ = NJ + 1
KJ(NJ) = K
7 CONTINUE
8 IF (NJ .EQ. 1) GO TO 9
OV = QUADR(J1,J2,0)
C1 = QUADR(I,J1,0)
C2 = QUADR(I,J2,0)
DEN = D1 - OV**2
A1 = -(C1 - C2*OV)/DEN
A2 = -(C2 - C1*OV)/DEN
WRITE (6,6) EL(I),EL(J1),C1,EL(I),EL(J2),C2,EL(J1),EL(J2),OV
D = SQRT(D1 + C1*A1 + C2*A2)
MX = MAX0(MAX(I),MAX(J1),MAX(J2))
DO 34 JJ = 1,MX
34 P(I,JJ) = (P(I,JJ) + A1*P(J1,JJ) + A2*P(J2,JJ))/D
AZ(I) = (AZ(I) + A1*AZ(J1) + A2*AZ(J2))/D
GO TO 35
9 D = QUADR(I,J,0)
IF ( ABS(D) .LT. 1.E-9 ) GO TO 3
DI = SQRT(D1 - D*D)
IF (P(I,1)-D*P(J,1) .LT. D0) DI = -DI
IMAX = MAX(I)
JMAX = MAX(J)
MX = MAX0(IMAX,JMAX)
DO 4 K=1,MX
4 P(I,K) = (P(I,K) - D *P(J,K))/DI
AZ(I) = (AZ(I) - D*AZ(J))/DI
WRITE(6,6) EL(I),EL(J),D
6 FORMAT(6X,1H(,A3,1H/,A3,4H) = ,1PE13.3)
C
C ***** COMPUTE AND STORE Y0(I,I)
C
35 CALL YKF(I,I,0)
DO 10 K = 1,NO
10 Y(I,K) = YK(K)
3 CONTINUE
2 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-22 O U T P U T
C ------------------------------------------------------------------
C
C THE RADIAL FUNCTIONS AND ORTHOGONALITY INTEGRALS ARE PRINTED,
C IF PRINT IS .TRUE. THE FUNCTIONS P/SQRT(R) WILL ALSO BE
C PUNCHED (OR STORED) ON UNIT OUF, IF OUF .NE. 0.
C
C
SUBROUTINE OUTPUT(PRINT)
REAL EL
LOGICAL PRINT
INTEGER OUC,OUD,OUF,OUH
COMMON /INOUT/ IUC,IUD,IUF,IUH,OUC,OUD,OUF,OUH
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
DIMENSION OUT(8)
IF ( .NOT. PRINT ) GO TO 30
C
C ***** PRINT RADIAL FUNCTIONS, 7 PER PAGE
C
ML = IB
2 MU = MIN0(ML+7,NWF)
I = MU - ML + 1
MX = 0
DO 1 J = ML,MU
1 MX = MAX0(MX,MAX(J))
WRITE (6,5) ATOM,TERM,(EL(J),J=ML,MU)
5 FORMAT(1H1,9X,19HWAVE FUNCTIONS FOR ,2A6//10X,1HR,8(10X,A3))
K= 0
KK = 0
DO 6 J = 1,MX
DO 9 JJ = ML,MU
IJ = JJ - ML + 1
9 OUT(IJ) = P(JJ,J)*R2(J)
K = K+1
IF (K .LE. 10) GO TO 6
K = 1
KK = KK+1
IF (KK .LT. 5) GO TO 21
KK = 0
WRITE (6,23)
23 FORMAT(1H1//)
GO TO 6
21 WRITE (6,8)
8 FORMAT(1X)
6 WRITE (6,10) R(J),(OUT(JJ),JJ=1,I)
10 FORMAT(F13.5,F15.6,7F13.6)
DO 15 J = ML,MU
IJ = J - ML + 1
15 OUT(IJ) = DPM(J)
WRITE (6,16) (OUT(J),J=1,I)
16 FORMAT(4X,10HMAX. DIFF. ,F15.7,7F13.7)
ML = ML+8
IF (ML .LE. NWF) GO TO 2
IF ( NWF .LE. 1) GO TO 30
C
C ***** PRINT ORTHOGONALITY INTEGRALS
C
WRITE (6,11) ATOM,TERM
11 FORMAT(////10X,33HORTHOGONALITY INTEGRALS FOR ATOM ,A6,6H TERM ,A6
1 //20X, 4H(NL),3X,4H(NL),7X,8HINTEGRAL //)
LM = IB
ML = MAX0(2,LM)
DO 12 I = ML,NWF
JF = I - 1
DO 13 J = 1,JF
IF (L(I) .NE. L(J)) GO TO 13
T = QUADR(I,J,0)
WRITE (6,17) EL(I),EL(J),T
17 FORMAT(21X,A3,4X,A3,F15.8)
13 CONTINUE
12 CONTINUE
30 IF ( OUF .EQ. 0) GO TO 14
C
C ***** OUTPUT FUNCTIONS ON UNIT OUF FOR FUTURE INPUT
C
DO 3 I = IB,NWF
MX = MAX(I)
WRITE (OUF,4) ATOM,TERM,EL(I),MX,Z,E(I,I),EK(I),AZ(I)
4 FORMAT( 6H ATOM ,A6,6H TERM ,A6,8HELECTRON ,A3,I6,F6.2/
1 4H E =,1PE14.7,4H I =,1PE14.7,4H AZ=,1PE14.7)
3 WRITE (OUF,7) (P(I,J),J =1,MX)
7 FORMAT(5(1PE14.7))
C
14 RETURN
END
C
C ------------------------------------------------------------------
C 3-23 P O T L
C ------------------------------------------------------------------
C
C COMPUTES AND STORES THE POTENTIAL FUNCTION
C 2(K-1)
C YR = SUM A Y (J,J;R)
C J,K IJK
C
SUBROUTINE POTL(I)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
DO 1 J=1,NO
1 YR(J) = D0
DO 2 J=1,NWF
C = A(I,J,1)
IF ( ABS(C) .LE. 1.E-10) GO TO 2
DO 3 JJ = 1,NO
3 YR(JJ) = C*Y(J,JJ) + YR(JJ)
DO 5 K = 2,5
C = A(I,J,K)
IF (ABS(C) .LE. 1.E-10) GO TO 5
CALL YKF(J,J,2*K-2)
DO 7 JJ=1,NO
7 YR(JJ) = YR(JJ) + C*YK(JJ)
5 CONTINUE
2 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-24 Q U A D
C ------------------------------------------------------------------
C
C EVALUATES THE INTEGRAL OF F(R)G(R) WITH RESPECT TO R , WHERE
C F(R) AND G(R) HAVE THE SAME ASYMPTOTIC PROPERTIES AS P (R). THE
C I
C COMPOSITE SIMPSON'S RULE IS USED. THE INTEGRAND IS ZERO FOR R >
C R .
C M
C
FUNCTION QUAD(I,M,F,G)
DIMENSION F(220),G(220)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
D = (D1 + D5*Z*R(1))/(H1*(2*L(I) + 3))
QUAD = RR(1)* F(1)*G(1)*( D -D5)
DO 1 J = 2,M,2
1 QUAD = QUAD + D2*RR(J)*F(J)*G(J)+ RR(J-1)*F(J-1)*G(J-1)
QUAD = QUAD*H1
RETURN
END
C
C ------------------------------------------------------------------
C 3-25 Q U A D R
C ------------------------------------------------------------------
C
C KK
C EVALUATES THE INTEGRAL OF R P (R) P (R) WITH RESPECT TO R
C I J
C
FUNCTION QUADR(I,J,KK)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
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 3-26 Q U A D S
C ------------------------------------------------------------------
C
C
C KK
C EVALUATES THE INTEGRAL OF (1/R) YK(R) P (R) P (R) WITH
C I J
C RESPECT TO R.
C
C
FUNCTION QUADS(I,J,KK)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
DEN = L(I) + L(J) + 3
K = 2 - KK
CD = D1 + Z*R(1)*(DEN-D1)/((DEN+D1)*((L(I)+1)*(L(J)+1)))
D = YK(1)*P(I,1)*P(J,1)*R(1)**K*( CD/(DEN*H1)+ D5)
MX = MIN0(MAX(I),MAX(J)) - 1
DO 1 M = 2,MX,2
1 D = D + D2*YK(M)*P(I,M)*P(J,M)*R(M)**K +
1 YK(M+1)*P(I,M+1)*P(J,M+1)*R(M+1)**K
QUADS = D*H1
RETURN
END
C
C ------------------------------------------------------------------
C 3-27 R K
C ------------------------------------------------------------------
C
C K
C EVALUATES R (I, J; II, JJ)
C
C
FUNCTION RK(I,J,II,JJ,K)
CALL YKF(I,II,K)
RK = QUADS(J,JJ,1)
RETURN
END
C
C ------------------------------------------------------------------
C 3-28 R O T A T E
C ------------------------------------------------------------------
C
C THIS ROUTINE ANALYSES THE ENERGY EXPRESSION TO DETERMINE THE
C STATIONARY CONDITION WITH RESPECT TO ROTATION OF ORBIALS I AND J.
C IF THE CONDITION IS ZERO, THE OFF-DIAGONAL ENERGY PARAMETERS MAY
C BE SET TO ZERO; OTHERWISE THE ORBITALS ARE ROTATED SO AS TO SATISFY
C THE STATIONAY CONDITION TO FIRST ORDER IN THE PERTURBATION.
C
C ROTATE SHOULD ONLY BE CALLED WHEN BOTH ORBITALS ARE BEING
C VARIED, AN ORTHOGONALITY CONSTRAINT APPLIES TO THE PAIR, AND NO
C OTHER ORBITALS ARE PRESENT WITH THE SAME PRINCIPAL QUANTUM NUMBERS.
C AN INCORRECT VALUE OF ORTHO = .FALSE. MAY CAUSE ROTATE TO BE CALLED
C WHEN THE ABOVE CONDITIONS DO NOT APPLY.
C
C
SUBROUTINE ROTATE(I,J)
REAL EL,QC
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /TEMP/QI,QJ,G,DI,DJ,DII,DJJ,DIJ,DJI,IDUM(3303)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
G = D0
DG = D0
IF (QI .EQ. D2*(2*L(I)+1) .AND. QJ .EQ. D2*(2*L(J)+1))
1 GO TO 44
IF (ABS(QI - QJ) .LT. 1.E-14) GO TO 16
C = D5*(QI - QJ)
G = G -C*HL(I,J)
DG = DG -C*(HL(I,I) - HL(J,J))
C
16 DO 13 K = 1,5
C = QI*(A(I,I,K) - A(I,J,K) - B(I,J,K))
IF (ABS(C) .LT. 1.E-10) GO TO 21
G = G + C*RK(I,I,I,J,2*(K-1))
FKII = FK(I,I,2*(K-1))
FKIJ = FK(I,J,2*(K-1))
GKIJ = GK(I,J,2*(K-1))
DG = DG +C*(FKII - FKIJ - D2*GKIJ)
21 CJ = QJ*(A(J,J,K) - A(J,I,K) - B(J,I,K))
IF (ABS(CJ) .LT. 1.E-10) GO TO 13
FKJJ = FK(J,J,2*(K-1))
IF (ABS(C) .GE. 1.E-10) GO TO 22
FKIJ = FK(I,J,2*(K-1))
GKIJ = GK(I,J,2*(K-1))
22 G = G - CJ*RK(J,J,J,I,2*(K-1))
DG = DG + CJ*(FKJJ -FKIJ - D2*GKIJ)
13 CONTINUE
DO 14 M = 1,NWF
IF (M .EQ. I .OR. M.EQ. J) GO TO 14
DO 15 K = 1,5
C = A(I,M,K)*QI - A(J,M,K)*QJ
IF (ABS(C) .LT. 1.E-10) GO TO 23
G = G + C*RK(I,M,J,M,2*(K-1))
DG = DG + C*(FK(I,M,2*(K-1)) - FK(J,M,2*(K-1)))
23 C = B(I,M,K)*QI - B(J,M,K)*QJ
IF (ABS(C) .LT. 1.E-10) GO TO 15
KK = IABS(L(I) - L(M)) + 2*(K-1)
G = G + C*RK(I,J,M,M,KK)
DG = DG + C*(GK(I,M,KK) - GK(J,M,KK))
15 CONTINUE
14 CONTINUE
C
C
C ***** ADD CONTRIBUTION FROM THE INTERACTION OF CONFIGURATIONS
C
IF (NR .EQ. 0) GO TO 11
DO 20 II = 1,NR
N1 = NCRI(II)
N2 = NCRJ(II)
C = WT(N1)*WT(N2)*CR(II)
CALL RPERT(I,J,II)
IF (.NOT. EZERO ) GO TO 45
G = G + D5*C*(DI - DJ)
DG = DG +D5*C*(DIJ - DJJ - DII + DJI)
20 CONTINUE
C
C
C ***** ADD CONTRIBUTION FROM THE INTERACTION OF CONFIGURATIONS - HL
C
11 IF (NL .EQ. 0) GO TO 101
DO 30 II = 1,NL
N1 = NCLI(II)
N2 = NCLJ(II)
C = WT(N1)*WT(N2)*CL(II)
CALL LPERT(I,J,II)
IF (.NOT. EZERO ) GO TO 45
G = G + D5*C*(DI - DJ)
DG = DG + D5*C*(DIJ - DJJ - DII + DJI)
30 CONTINUE
101 WRITE (6,100) G,DG
100 FORMAT(10X,12HCONDITION = ,F10.6,4X,12HVARIATION = ,F10.6)
IF (ABS(G) .GT. 1.E-9 .OR. ABS(E(I,J)) .GT. 2.E-10) GO TO 40
G = D0
44 E(I,J) = D0
RETURN
40 EPS = G/DG
EZERO = .FALSE.
DD = SQRT(D1 + EPS*EPS)
DO 41 JJ = 1,NO
PI = (P(I,JJ) + EPS*P(J,JJ))/DD
P(J,JJ) = (P(J,JJ) - EPS*P(I,JJ))/DD
41 P(I,JJ) = PI
CALL YKF(I,I,0)
DO 42 JJ=1,NO
42 Y(I,JJ) = YK(JJ)
CALL YKF(J,J,0)
DO 43 JJ = 1,NO
43 Y(J,JJ) = YK(JJ)
RETURN
C
C ***** ONE OF THE PAIR OF ORBITALS OCCURS IN AN OVERLAP INTEGRAL AND
C ***** THE PAIR CANNOT BE ROTATED.
C
45 G = D0
RETURN
END
C
C ------------------------------------------------------------------
C 3-29 R P E R T
C ------------------------------------------------------------------
C
C RPERT DETERMINES THE EFFECT OF A PERTURBATION IN THE FORM OF A
C ROTATION OF THE I'TH AND J'TH ORBITAL BOTH ON THE ENERGY AND ON THE
C STATIONARY CONDITION, DUE TO THE PRESENCE OF A GIVEN RK INTEGRAL.
C
C
SUBROUTINE RPERT(IIN,JIN,M)
REAL EL,QC
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /TEMP/QI,QJ,G,DI,DJ,DII,DJJ,DIJ,DJI,IDUM(3303)
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
DIMENSION IND(6)
EQUIVALENCE (I1,IND(1)),(I2,IND(2)),(I3,IND(3)),(I4,IND(4)),
1 (I5,IND(5)),(I6,IND(6))
I1 = I1R(M)
I2 = I2R(M)
I3 = J1R(M)
I4 = J2R(M)
KK = KR(M)
J = IIN
I = JIN
IPW = IQ(M)
IF (IPW .EQ. 0) GO TO 1
I5 = IO(M)
I6 = JO(M)
OVLAP = QUADR(I5,I6,0)**IPW
IF ((I5-I)*(I6-I)*(I5-J)*(I6-J) .NE. 0) GO TO 1
EZERO = .FALSE.
RETURN
1 DO 10 KP = 1,2
DI = D0
DII = D0
DIJ = D0
DO 2 K = 1,4
IF (IND(K) .NE. I) GO TO 2
IND(K) = J
DI = DI + RK(I1,I2,I3,I4,KK)
IF (.NOT. ALL) GO TO 3
DO 4 K2 = 1,4
IF (IND(K2) .NE. J) GO TO 5
IND(K2) = I
DIJ = DIJ + RK(I1,I2,I3,I4,KK)
IND(K2) = J
5 IF (IND(K2) .NE. I) GO TO 4
IND(K2) = J
DII = DII + RK(I1,I2,I3,I4,KK)
IND(K2) = I
4 CONTINUE
3 IND(K) = I
2 CONTINUE
IF (IPW .EQ. 0) GO TO 6
DI = DI*OVLAP
DII = DII*OVLAP
DIJ = DIJ*OVLAP
6 IF (KP .EQ. 2) RETURN
DJ = DI
DJJ = DII
DJI = DIJ
I = IIN
J = JIN
10 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-30 S C A L E
C ------------------------------------------------------------------
C
C THE CURRENT RADIAL FUNCTIONS ARE SCALED ACCORDING TO THE
C PROCEDURES OF SEC. 7-2 . VALUES OF AZ AND E(I,I), THE STARTING
C VALUES AND THE DIAGONAL ENERGY PARAMETERS ARE ALSO SCALED.
C
C
SUBROUTINE SCALE(ZZ)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /TEMP/RATIO,SR,SC,SS,F0,F1,F2,F3,PNORM,THETA,K,IDUM(3300)
DIMENSION RS(220),PS(220)
EQUIVALENCE (RS(1),YR(1)),(PS(1),X(1))
C
C ***** SCALE VALUES OF R=RS, P=PS AND ONE-ELECTRON PARAMETERS.
C ***** GENERATE NEW VALUES OF R, R*R, AND SQRT(R)
C
RATIO = Z/ZZ
SR = SQRT(RATIO)
DO 1 J = 1,NO
R(J) = R(J)*RATIO
RR(J) = R(J)*R(J)
1 R2(J) = R2(J)*SR
DO 2 I = 1,NWF
SC = (ZZ-S(I))/(Z-S(I))
SS = SC*RATIO
E(I,I) = E(I,I)*SC**2
a write(6,499) ss, sc, sr, ratio, zz, z
499 format(' ss, sc, sr, ratio, zz, z',6e10.3)
DO 3 J = 1,NO
RS(J) = R(J)/SS
3 PS(J) = P(I,J)*SC
SC = (ZZ - D5*S(I))/(Z - D5*S(I))
AZ(I) = AZ(I)*SC**(L(I)+1)*SQRT(SC)
a if (sc.lt.d0) write(6,500) sc,zz,z,s(i)
500 format(' test of sc line 2524 ',4f10.8)
K = 3
C
C ***** INTERPOLATE THE (RS,PS) FUNCTIONS FOR VALUES OF P AT THE SET
C ***** OF POINTS R
C
DO 4 J = 1,NO
C
C ***** SEARCH FOR THE NEAREST ENTRIES IN THE (RS,PS) TABLE
C
5 IF (K .EQ. ND) GO TO 7
IF (RS(K) .GT. R(J)) GO TO 6
K = K + 1
GO TO 5
C
C ***** INTERPOLATE
C
6 THETA = ALOG(R(J)/RS(K-1))/H
F0 = PS(K-2)
F1 = PS(K-1)
F2 = PS(K)
F3 = PS(K+1)
P(I,J) = D5*(F1+F2) + (THETA -D5)*(F2 - F1) +
1 THETA*(THETA - D1)*(F0 - F1 - F2 + F3)/D4
GO TO 4
7 P(I,J) = D0
4 CONTINUE
MAX(I) = NO
C
C *****NORMALIZE THE INTERPOLATED FUNCTION
C
PNORM = SQRT(QUADR(I,I,0))
a if (pnorm.le.0.d0) write(6,501) pnorm
501 format(' pnorm error ',f10.7)
DO 10 J = 1,NO
10 P(I,J) = P(I,J)/PNORM
C
C ***** COMPUTE NEW VALUES OF Y0(I,I) AND STORE
C
CALL YKF(I,I,0)
DO 11 J = 1,NO
11 Y(I,J) = YK(J)
2 CONTINUE
Z = ZZ
RETURN
END
C
C ------------------------------------------------------------------
C 3-31 S C F
C -----------------------------------------------------------------
C
C THIS ROUTINE CONTROLS THE SCF PROCEDURE DESCRIBED IN CHAPTER
C 7. IF CERTAIN INPUT PARAMETERS ARE ZERO (OR BLANK) THEY WILL BE
C SET TO THEIR DEFAULT VALUE.
C
C PARAMETER DEFAULT VALUE
C -------- -------------
C CFGTOL 1.E-10
C SCFTOL 1.E-7
C IC (NWF + 1 - IB)/4 + 3
C NSCF 12
C
C THE SELF-CONSISTENCY CONVERGENCE CRITERION IS
C
C Z2 = SQRT( SCFTOL*(Z*NWF/2) )
C
C IT IS INCREASED BY A FACTOR TWO AT THE END OF EACH ITERATION WHEREAS
C CFGTOL IS INCREASED BY SQRT(2).
C
C
SUBROUTINE SCF(ETOTAL,ACFG,SCFTOL,CFGTOL,LD,NEW)
REAL EL
LOGICAL FAIL,RSCAN(20),OMIT,ECONV,ORTHO,EZERO,LD,LAST,CONV,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
INTEGER OUC,OUD,OUF,OUH
COMMON /INOUT/ IUC,IUD,IUF,IUH,OUC,OUD,OUF,OUH
COMMON /HMATRX/ ET(40,40)
C
C ***** DETERMINE THE NUMBER OF NEW CONFIGURATION STATES
C
IF (NEW .EQ. 0) NEW = NCFG
NC = NCFG - NEW +1
C
C ***** READ HMATRIX FROM UNIT IUH , IF DIFFERENT FROM ZERO
C
IF ( IUH .EQ. 0 .OR. NC .LE. 1 ) GO TO 101
NCM = NC - 1
DO 102 I = 1,NCM
READ (IUH,103) (ET(I,J),J=1,I)
103 FORMAT(5F14.7)
C
C ***** SYMMETRIZE THE MATRIX
C
DO 104 J = 1,I
104 ET(J,I) = ET(I,J)
102 CONTINUE
C
C ***** SET THE SCF CONVERGENCE PARAMETER TO AN OPTIMISTIC VALUE
C
c eq in next two statements corrected to le
101 IF ( CFGTOL .LE. D0) CFGTOL = .1E-9
IF ( SCFTOL .LE. D0) SCFTOL = 1.E-7
IF ( IC .EQ. 0 .AND. NCFG .EQ. 1 ) IC = 3 + (NWF + 1 - IB)/4
TOL = SQRT(Z)*1.E-10
Z2 = SCFTOL*SQRT(Z*NWF)
WRITE (6,15)
15 FORMAT(//)
WRITE (6,16) ORTHO,OMIT,ACFG,SCFTOL,NO
16 FORMAT(10X,44HORTHOGONALITY BETWEEN CONFIGURATIONS =,L4/
1 10X,44HWEAK ORTHOGONALIZATION DURING THE SCF CYCLE=,L4/
2 10X,44HACCELERATING PARAMETER FOR MCHF ITERATION =,F5.2/
3 10X,44HSCF CONVERGENCE TOLERANCE (FUNCTIONS) =,1PE9.2
4 ,/10X,44HNUMBER OF POINTS IN THE MAXIMUM RANGE =,I4)
C
C ***** SET ITERATION PARAMETERS
C
IPR = 0
ECONV = .FALSE.
LAST = .FALSE.
IF (IB .GT. NWF ) LAST = .TRUE.
DP1 = D0
ETOTAL = D0
IF (NSCF .EQ. 0) NSCF = 12
IF ( .NOT. LD ) GO TO 9
CALL DIAG(ETOTAL,ECONV,ACFG,CFGTOL,NC,LAST)
IF (LAST) RETURN
IF ( NCFG .GT. 1 .AND. ID .NE. 1) CALL ARRAY
C
C ***** PERFORM NSCF SELF-CONSISTENT FIELD ITERATIONS
C
9 DO 100 I = 1,NSCF
WRITE(6,7) I,CFGTOL,Z2
7 FORMAT(//10X,17HITERATION NUMBER ,I2/10X,16H----------------//
1 10X,50HCONVERGENCE CRITERIA:ENERGY (CFGTOL) =,1PE9.1/
2 10X,50H :FUNCTION(SCFTOL*SQRT(Z*NWF))=,1PE9.1)
DP1 = D0
IF (IB .GT. NWF) GO TO 12
CALL GRANGE
C
C ***** SOLVE EACH DIFFERENTIAL EQUATION IN TURN
C
WRITE(6,14)
14 FORMAT(/20X,3H EL,9X,2HED,13X,2HAZ,11X,4HNORM,7X,3HDPM)
DO 2 J = IB,NWF
CALL DE(J)
IF ( FAIL ) GO TO 3
DP = DPM(J)*SQRT(SUM(J))
IF ( DP1 .GE. DP ) GO TO 2
DP1 = DP
JJ = J
2 CONTINUE
IF ((NCFG .EQ. 1 .OR. ID .EQ. 1) .AND. DP1 .LT. Z2) GO TO 6
IF ( IC .LE. 0) GO TO 6
C
C ***** SOLVE IC DIFFERENTIAL EQUATIONS EACH TIME SELECTING THE
C ***** ONE WITH THE LARGEST DPM
C
DO 4 II =1,IC
CALL DE(JJ)
IF ( FAIL ) GO TO 3
DP1 = D0
DO 5 J = IB,NWF
DP = SQRT(SUM(J))*DPM(J)
IF ( DP1 .GT. DP ) GO TO 5
JJ = J
DP1 = DP
5 CONTINUE
IF (DP1 .LT. Z2) GO TO 6
4 CONTINUE
6 CALL ORTHOG
IF (.NOT.(NCFG .EQ. 1 .OR. ID .EQ. 1)) GO TO 12
IF (DP1 .GT. Z2) GO TO 1
IF (DP1 .LE. Z2 .AND. LAST ) GO TO 12
LAST = .TRUE.
GO TO 1
12 CALL DIAG(ETOTAL,ECONV,ACFG,CFGTOL,NC,LAST)
IF ( NCFG .GT. 1 .AND. ID .NE. 1) CALL ARRAY
C
C ***** IF FUNCTIONS APPEAR TO HAVE CONVERGED,SOLVE EACH AGAIN, IN
C ***** TURN, AND TEST AGAIN
C
IF (IB .GT. NWF) RETURN
CONV = ECONV .AND. DP1 .LE. Z2
IF (CONV .AND. LAST) RETURN
IF (CONV) LAST =.TRUE.
C
C ***** INCREASE THE CONVERGENCE CRITERION FOR SELF-CONSISTENCY
C
1 Z2 = D2*Z2
WRITE(6,8) EL(JJ),DP1
8 FORMAT(/10X,34HLEAST SELF-CONSISTENT FUNCTION IS ,A3,
1 27H (WEIGHTED MAXIMUM CHANGE =,1PE10.2,1H))
100 CFGTOL = 1.4E0*CFGTOL
WRITE (6,13)
13 FORMAT(10X,46HSCF ITERATIONS HAVE NOT CONVERGED TO REQUESTED
1 ,9H ACCURACY)
FAIL = .TRUE.
3 RETURN
END
C
C ------------------------------------------------------------------
C 3-32 S E A R C H
C ------------------------------------------------------------------
C
C THIS ROUTINE SEARCHES FOR THE NJ>70 SUCH THAT YR(J) > 0 FOR ALL
C J > NJ.
C
C
SUBROUTINE SEARCH(NJ,I)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
IA = 70
IL = NO
4 IF (YR(IA) .LT. D0) GO TO 3
IA = IA + 2
IF (IA .LT. IL ) GO TO 4
NJ = MAX0(70,MAX(I)-100)
RETURN
3 NK = (IA + IL)/2
IF (YR(NK) .LT. D0) GO TO 1
IL = NK
GO TO 2
1 IA = NK
2 IF (IL - IA .GT. 1) GO TO 3
NJ = IL - 7
RETURN
END
C
C ------------------------------------------------------------------
C 3-33 S M
C ------------------------------------------------------------------
C
C 3 K
C EVALUATES THE INTEGRAL OF (1/R) P (R) P (R) Z (J, J; R) WITH
C I I
C RESPECT TO R.
C
C
FUNCTION SM(I,J,K)
CALL ZK(J,J,K)
SM = QUADS(I,I,3)*2.921787
RETURN
END
C
C ------------------------------------------------------------------
C 3-34 S N
C ------------------------------------------------------------------
C
C 3 K
C EVALUATES THE INTEGRAL OF (1/R) P (R) P (R) Z (I, J; R) WITH
C I J
C RESPECT TO R.
C
FUNCTION SN(I,J,K)
CALL ZK(I,J,K)
SN = QUADS(I,J,3)*2.921787
RETURN
END
C
C ------------------------------------------------------------------
C 3-35 S O L V E
C ------------------------------------------------------------------
C
C WHEN FIRST IS .TRUE., SOLVE COMPUTES THE POTENTIAL AND EXCHANGE
C FUNCTION AND INITIALIZES VARIABLES FOR THE I'TH RADIAL EQUATION.
C THE VECTOR P1 IS THE SOLUTION OF THE RADIAL EQUATION AND P2 THE
C VARIATION OF THE SOLUTION WITH RESPECT TO THE ENERGY PARAMETER
C E(I,I).
C
C
SUBROUTINE SOLVE(I,FIRST)
REAL EL
LOGICAL ONE,FIRST
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /TEMP/ P2(220),HQ(220),XX(220),AZZ,PP,FN,EM,FM,EU,FU,
1 DELTAE,M,NODE,MK,KK,NJ,IDUM(1980)
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
DIMENSION ZERO(220),P1(220)
EQUIVALENCE (ZERO(1),XX(1)),(PDE(1),P1(1))
C
C ***** IF FIRST IS 'TRUE', CALL POTL AND XCH AND SET UP ARRAYS
C
IF (.NOT. FIRST) GO TO 17
CALL POTL(I)
CALL XCH(I,3)
FN = N(I)
FL = L(I)
V = YR(1)/R(1)
B4 = Z*(FL+1.3333333333)/((FL+D1)*(FL+D2))
CN = (D2*Z/FN)**(L(I) +1)
C = D4*FL +D6
CD = (FL+D5)**2
XY = X(1)
XP = X(2)
ED = E(I,I)
X1 = X(1)
X2 = X(2)
X3 = X(3)
X4 = X(4)
DO 1 J = 3,ND
X5 = X(J+2)
X(J) =CH*(-X5+24.E0*(X4+X2) + 194.E0*X3 - X1)/20.E0
X1 = X2
X2= X3
X3 = X4
1 X4 = X5
X(NO-1) = CH*(X4 + D10*X3 + X2)
DO 4 J = 1,NO
4 YK(J) = -D2*(Z - YR(J))*R(J) + CD
X1 = CH*P(I,1)*(YK(1)+ED*RR(1))
X2 = CH*P(I,2)*(YK(2)+ED*RR(2))
X3 = CH*P(I,3)*(YK(3)+ED*RR(3))
X4 = CH*P(I,4)*(YK(4)+ED*RR(4))
DO 7 J = 3,ND
X5 = CH* P(I,J+2)*(YK(J+2)+ED*RR(J+2))
X(J) = X(J) - (X5 -D4*(X2 + X4) + D6*X3 +X1)/20.E0
X1 = X2
X2 = X3
X3 = X4
7 X4 = X5
RL = L(I) + 2.5
X(2) = R(2)**RL*(X(5)/R(5)**RL - D3*(X(4)/R(4)**RL -
1 X(3)/R(3)**RL))
C
C ***** DETERMINE LOWER BOUND ON THE ENERGY PARAMETER
C
DO 11 JJ = 15,ND
J = NO - JJ
IF (YK(J) .LT. D0) GO TO 63
11 CONTINUE
WRITE ( 6,12)
12 FORMAT(10X,49HPOTENTIAL FUNCTION TOO SMALL - 2R*(Z-Y)<(L+.5)**2)
STOP
63 EM = -YK(J)/RR(J)
FM = EM
C
C ***** DETERMINE DIAGONAL ENERGY PARAMETER
C
F1 = D0
C11 = D0
M = MIN0(MAX(I),NO-1)
DO 5 J = 2,M
FNUM = P(I,J+1) - P(I,J) - P(I,J) + P(I,J-1) -CH*( YK(J+1)*
1 P(I,J+1) + D10*YK(J)*P(I,J) + YK(J-1)*P(I,J-1))-X(J)
DEL1 = RR(J+1)*P(I,J+1) + D10*RR(J)*P(I,J) + RR(J-1)*P(I,J-1)
F1 = F1 +P(I,J)*FNUM
C11 = C11 + P(I,J)*DEL1
5 CONTINUE
ED = F1/(C11*CH)
IF (ED .GT. EM) GO TO 19
C
C ***** ERROR MESSAGE AND ENERGY ADJUSTMENT FOR AN ENERGY PARAMETER
C ***** TOO SMALL FOR THE RANGE OF THE FUNCTION
C
WRITE (6,24) ED
24 FORMAT(10X,5HED = ,F10.6,36H; ADJUSTED TO ALLOWED MINIMUM ENERGY )
ED = EM
IF ( ABS(FM - E(I,I)) .GT. 1.E-6 ) GO TO 19
WRITE (6,64) EL(I)
64 FORMAT(//10X,14HITERATION FOR ,A3,24H MAY BE CONVERGING TO A ,
1 18HCONTINUUM FUNCTION/10X,
2 39HCHECK THAT E(CORE)-E(TOTAL) IS POSITIVE )
C
C ***** CHECK IF UPPER BOUND IS CORRECT
C
19 IF ( D10*ED .LT. EU) GO TO 18
EU = D10*ED
FU = EU
18 AZD = AZ(I)
17 DO 26 J=1,NO
YR(J) = (YK(J) + ED*RR(J))*CH
26 ZERO(J) = D0
C
C ***** SEARCH FOR THE POINT AT WHICH YR BECOMES POSITIVE
C
CALL SEARCH(NJ,I)
C
C ***** COMPUTE STARTING VALUES FROM SERIES EXPANSION
C
B3 = (V + V + ED - (Z/FN)**2)/C
DO 6 J = 1,2
HW = HWF(N(I),L(I),Z,R(J))/CN
6 HQ(J) = AZD*(HW + R(J)**(L(I)+3)*B3*(D1-R(J)*B4))/R2(J)
C
C ***** OBTAIN HOMOGENEOUS SOLUTION
C
CALL NMRVS(NJ,DELH,MH,HQ,ZERO)
P1(1) = HQ(1) + XY/C
P1(2) = HQ(2) + XP/C
C
C ***** OBTAIN PARTICULAR SOLUTION
C
CALL NMRVS(NJ,DEL1,M1,P1,X)
C
C ***** DETERMINE THE ENERGY ADJUSTMENT REQUIRED FOR A SOLUTION WITH
C ***** GIVEN A0
C
M = MAX0(M1,MH)
PNORM = D0
DO 50 J = 1,M
50 PNORM = PNORM + RR(J)*HQ(J)*P1(J)
Y1 = P1(NJ-1)
Y2 = P1(NJ)
Y3 = P1(NJ+1)
DELTA = Y2 - Y1 + Y2 - Y3 +YR(NJ-1)*Y1 + D10*YR(NJ)*Y2
1 + YR(NJ+1)*Y3 + X(NJ)
DELTAE = HQ(NJ)*DELTA/(H*H*PNORM)
PP = -DEL1/DELH
C
C ***** MATCH AT THE JOIN FOR A SOLUTION OF THE DIFFERENTIAL EQUATION
C
DO 13 J = 1,NO
13 P1(J) = P1(J) + PP*HQ(J)
C
C ***** IF THE EQUATIONS APPEAR TO BE NEARLY
C **** SINGULAR, SOLVE THE VARIATIONAL EQUATIONS
C
IF (KK .NE. 2) RETURN
X1 = P(I,1)*RR(1)
X2 = P(I,2)*RR(2)
P2(1) = X1/C
P2(2) = X2/C
DO 8 J = 3,NO
X3 = P(I,J)*RR(J)
XX(J-1) = (D10*X2 + X1 + X3)*CH
X1 = X2
8 X2 = X3
CALL NMRVS(NJ,DEL2,M2,P2,XX)
AA = -DEL2/DELH
M = MAX0(M,M2)
DO 9 J = 1,NO
9 P2(J) = P2(J) + AA*HQ(J)
A11 = QUAD(I,M,P2,P2)
B11 = QUAD(I,M,P1,P2)
C11 = QUAD(I,M,P1,P1) - D1
DISC = B11*B11 - A11*C11
IF ( DISC .LT. D0 ) GO TO 70
DE1 = -(B11+SQRT(DISC))/A11
DE2 = C11/A11/DE1
IF ( P1(3)+DE1*P2(3) .LT. D0) DE1 = DE2
GO TO 71
70 DE1 = C11/A11
71 DO 301 J = 1,NO
P1(J) = P1(J) + DE1*P2(J)
301 CONTINUE
PP = PP + DE1*AA
RETURN
END
C
C ------------------------------------------------------------------
C 3-36 S U M M R Y
C ------------------------------------------------------------------
C
C THE RESULTS OF A CALCULATION ARE SUMMARIZED. THESE INCLUDE
C THE FOLLOWING FOR EACH ELECTRON:
C
C E(NL) - DIAGONAL ENERGY PARAMETER
C AZ(NL) - STARTING PARAMETER, P(R)/R**(L+1) AS R -> 0.
C SIGMA - SCREENING PARAMETER AS DEFINED BY EQ. (7-2) .
C 1/R**3 - EXPECTED VALUE OF <1/R**3>
C 1/R - EXPECTED VALUE OF <1/R>
C R - EXPECTED MEAN RADIUS
C R**2 - EXPECTED VALUE OF <R**2>
C I(NL) - -(1/2)<NL!L!NL>
C KE - I(NL) + Z <R>
C
C THESE RESULTS ARE FOLLOWED BY:
C
C TOTAL ENERGY (ET)
C KINETIC ENERGY (EK)
C POTENTIAL ENERGY (EP) = ET - EN
C RATIO = - EP/EN
C K K K
C THE VALUES OF ALL F , G , R AND <NL!L!N'L> INTEGRALS WHICH ENTER
C INTO THE CALCULATION ARE PRINTED, IF OUD > 0.
C IF OUD > 0, THE SPIN-ORBIT PARAMETERS ZETA, AS DEFINED BY BLUME
C AND WATSON (1962) ARE COMPUTED FOR EACH CONFIGURATION.
C THIS CALCULATION AGREES WITH THEIR FORMULAS FOR CONFIGURATIONS
C WITH ONLY ONE OPEN SHELL (THE ONLY CASE CONSIDERED IN THEIR
C PAPER). WHEN MORE THAN ONE OPEN SHELL IS PRESENT, THE
C "EXCHANGE PART" BETWEEN OPEN SHELLS IS IGNORED. THE SPIN-SPIN
C PARAMETERS ARE ALSO PRINTED.
C
SUBROUTINE SUMMRY(ET)
REAL EL,QC
INTEGER OUC,OUD,OUF,OUH
LONG D,F,G
C DATA D/4H /,F/1HF/,G/1HG/
COMMON /INOUT/ IUC,IUD,IUF,IUH,OUC,OUD,OUF,OUH
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /TEMP/R3(20),SS(3),R1,RM,RMM,RH,SC,QI,QJ,SP,C,CC,OUT(5),
1 EKINP,EN,EPOT,RATIO,LI,LJ,K,KF,I1,I2,J1,J2,I,J,MIN,IW(5),
2 IDUM(3221)
COMMON /PRIVATE/ CV(3,4,3),CN(3,4,4)
DATA CV/-1.0E0,-.6E0,0.E0,-.6E0,-.2E0,-.1714285714285714,.2E0,
1 2*-.08571428571428571,.2571428571428571,.08571428571428571,
2 -.028571428571428571,2*0.E0,-.4285714285714286,0.E0,
3 -.3428571428571428,-.2380952380952380,-.5142857142857143,
4 -.2857142857142857,-.08571428571428571,-.4761904761904762,
5 -.0571428571428571,-.0649350649350649,8*0.E0,-.1948051948051948
6 ,0.E0,-.2597402597402597,-.1748251748251748/
DATA CN/ 2.E0,1.8E0,0.E0,-1.E0,.8E0,.85714285714285714,.8E0,-.6E0,
1 .5142857142857143,1.285714285714286,.5142857142857143,
2 -.4285714285714286,-1.E0,-1.2E0,1.714285714285714,0.E0,
3 -.1428571428571428,-.04761904761904762,-1.6857142857142857,
4 .1714285714285714,-.2571428571428571,-2.476190476190476,
5 -.7142857142857142,-.1428571428571428,2*0.E0,
6 -1.285714285714286,0.E0,-.3428571428571428,-.4761904761904762,
7 .5142857142857143,0.E0,.03376623376623377,.9523809523809524,
8 -.3688311688311688,-.0649350649350649,8*0.E0,-.1948051948051948
9 ,0.E0,.2597402597402597,0.E0/
D=' '
F='F'
G='G'
WRITE (6,9) ATOM,TERM
9 FORMAT(/// 24X,5HATOM ,A6,3X,5HTERM ,A6//67X,13HMEAN VALUE OF,24X,
1 22HONE ELECTRON INTEGRALS /5X,2HNL,7X,5HE(NL),9X,6HAZ(NL),
2 5X,5HSIGMA,8X,6H1/R**3,8X,3H1/R,10X,1HR,13X,4HR**2,10X,
3 5HI(NL),9X,2HKE )
EN = D0
C
C ***** COMPUTE AND PRINT ONE-ELECTRON PARAMETERS
C
DO 10 I = 1,NWF
R1 = QUADR(I,I,-1)
RM = QUADR(I,I,1)
RMM = QUADR(I,I,2)
EKINP = EK(I) + Z*R1
EN = EN+ SUM(I)*EKINP
RH = 3*N(I)*N(I) - L(I)*(L(I) + 1)
SC = Z - D5*RH/RM
S(I) = SC
R3(I) = D0
IF (L(I) .NE. 0) R3(I) = QUADR(I,I,-3)
WRITE (6,15)EL(I),E(I,I),AZ(I),SC,R3(I),R1,RM,RMM,EK(I),EKINP
15 FORMAT( 4X,A3,F14.7,F15.7,F8.3,2F14.5,4F14.7)
10 CONTINUE
IF ( NL .EQ. 0 ) GO TO 31
C
C ***** ADD CONTRIBUTION FROM THE 'L' INTEGRALS
C
DO 32 I = 1,NL
N1 = NCLI(I)
N2 = NCLJ(I)
CC = CL(I)*WT(N1)*WT(N2)
IF (LQ(I) .NE. 0) CC = CC*QUADR(ILO(I),JLO(I),0)**LQ(I)
CONT = CC*(HL(ILI(I),ILJ(I)) - D2*Z*QUADR(ILI(I),ILJ(I),-1))
32 EN = EN + CONT
31 EPOTL = ET - EN
RATIO =-EPOTL/EN
WRITE (6,26) ET,EN,EPOTL,RATIO
26 FORMAT(/4X,15HTOTAL ENERGY = ,F14.7,4X,17HKINETIC ENERGY = ,F14.7,
1 4X,19HPOTENTIAL ENERGY = ,F14.7 ,4X,8HRATIO = ,F11.9)
IF ( OUD .EQ. 0 ) RETURN
WRITE (OUD,126)
126 FORMAT(/2X,27HVALUES OF F AND G INTEGRALS /)
C
C ***** PRINT TABLES OF 'FK' AND 'GK' INTEGRALS WHICH WERE USED IN
C ***** DETERMINING THE ENERGY
C
DO 17 J = IB,NWF
DO 20 I = 1,J
KF = 0
MIN = 0
DO 16 K = 1,5
IF (A(I,J,K) .EQ. D0) GO TO 18
IW(K) = MIN
OUT(K) = FK(I,J,MIN)
MIN = MIN + 2
16 KF = KF + 1
18 IF (KF .NE. 0) WRITE(OUD,19) (D,F,IW(K),EL(I),EL(J),OUT(K),K=1,KF)
19 FORMAT( 2(2X,2A1,I1,1H(,A3,1H,,A3,4H ) =, F10.7,2X))
MIN = IABS(L(I) - L(J))
KF = 0
DO 24 K = 1,5
IF (B(I,J,K) .EQ. D0) GO TO 25
IW(K) = MIN
OUT(K) = GK(I,J,MIN)
MIN = MIN + 2
24 KF = KF + 1
25 IF (KF .NE. 0) WRITE(OUD,19) (D,G,IW(K),EL(I),EL(J),OUT(K),K=1,KF)
20 CONTINUE
17 CONTINUE
IF (NR .EQ. 0) GO TO 27
C
C ***** PRINT TABLES OF 'RK' INTEGRALS
C
WRITE (OUD,21)
21 FORMAT(//2X,21HVALUES OF R INTEGRALS //)
DO 22 I = 1,NR
I1 = I1R(I)
I2 = I2R(I)
J1 = J1R(I)
J2 = J2R(I)
OUT(1) = RK(I1,I2,J1,J2,KR(I))
22 WRITE (OUD,23) KR(I),EL(I1),EL(I2),EL(J1),EL(J2),OUT(1)
23 FORMAT(2X,1HR,I1,2H (,2A3,1H,, 2A3,4H ) =, F10.7 )
27 IF (NL .EQ. 0) GO TO 13
C
C ***** PRINT TABLES OF 'L' INTEGRALS
C
WRITE (OUD,28)
28 FORMAT(2X,21HVALUES OF L INTEGRALS //)
DO 29 I = 1,NL
I1 = ILI(I)
J1 = ILJ(I)
OUT(1) = HL(I1,J1)
29 WRITE(OUD,30) EL(I1),EL(J1),OUT(1)
30 FORMAT(2X,2HL(,A3,1H,,A3,4H) = ,F16.7)
13 WRITE (OUD,1)
1 FORMAT(//3X,20HSPIN-ORBIT PARAMETER,10X,20HSPIN-SPIN PARAMETERS/)
C
C ***** COMPUTE AND PRINT SPIN-ORBIT AND SPIN-SPIN PARAMETERS
C
DO 2 K = 1,NCFG
DO 7 I = 1,NWF
QI = QC(I,K)
LI = L(I)
IF (QI .EQ. 4*LI+2 .OR. LI .EQ. 0 .OR. QI .EQ. D0) GO TO 7
IF (LI .EQ. 4) GO TO 7
SP = Z*R3(I)*5.843574
DO 4 J = 1,NWF
IF (L(J) .EQ. 4) GO TO 4
QJ = QC(J,K)
IF (QJ .EQ. D0 .OR. I .EQ. J) GO TO 4
LJ = L(J) + 1
C = QJ
SP = SP - D2*C*SM(I,J,0)
IF (QJ .NE. 4*LJ-2) GO TO 4
DO 5 KK = 1,3
CC = CV(LI,LJ ,KK)
IF (CC .NE. D0) SP = SP - C*CC*V(I,J,(LI+LJ)+2*(KK-1-(LI+LJ)/2))
5 CONTINUE
DO 6 KK = 1,4
CC = CN(LI,LJ ,KK)
IF (CC .NE. D0) SP = SP - C*CC*SN(I,J,(LI+LJ)-3+2*(KK -(LI+LJ)/2))
6 CONTINUE
4 CONTINUE
C= QI
IF (C .EQ. D1) GO TO 8
SS(1) = SM(I,I,0)
C = C + C - D3
SP = SP - C*SS(1)
GO TO (3,11,12), LI
11 SS(2) = SM(I,I,2)
SP = SP + .85714142857*SS(2)
GO TO 3
12 SS(2) = SM(I,I,2)
SS(3) = SM(I,I,4)
SP = SP + SS(2) + .45454545455*SS(3)
3 WRITE (OUD,14) K,EL(I),SP,(D,EL(I),EL(I),SS(II),II=1,LI)
14 FORMAT( 1X,I3,6H ZETA(,A3,3H) =,F11.3,3X,A4,3HM0(,A3,1H,,A3,
1 4H) = , F7.4,A4/34X,3HM2(,A3,1H,,A3,4H) = ,F7.4
2 ,A4/34X,3HM4(,A3,1H,,A3,4H) = ,F7.4)
GO TO 7
8 WRITE (OUD,14) K,EL(I),SP
7 CONTINUE
2 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-37 V
C ------------------------------------------------------------------
C
C K
C EVALUATES V (I,J) AS DEFINED BY BLUME AND WATSON (1962).
C
FUNCTION V(I,J,K)
CALL DYK(I,J,K)
V = QUADS(I,J,2)
CALL DYK(J,I,K)
V = (V - QUADS(I,J,2))*2.921787
RETURN
END
C
C ------------------------------------------------------------------
C 3-38 W A V E F N
C ------------------------------------------------------------------
C
C THIS ROUTINE INITIALIZES RADIAL FUNCTIONS BY THE PROCEDURE
C INDICATED BY IND(I).
C
C VALUE OF IND(I) METHOD
C --------------- ------
C -1 FUNCTIONS READ FROM UNIT IU2
C 0 SCREENED HYDROGENIC FUNCTIONS WITH ZZ=Z-S(I)
C 1 FUNCTIONS IN MEMORY LEFT UNCHANGED
C 0
C THE SET OF FUNCTIONS ARE THEN ORTHOGONALIZED, Y (I, I;R) AND THE
C DIAGONAL ENERGY PARAMETERS COMPUTED, WHEN NECESSARY.
C
C
SUBROUTINE WAVEFN
REAL EL
LOGICAL FAIL,RSCAN(20),OMIT,ECONV,ORTHO,EZERO,ALL
INTEGER OUC,OUD,OUF,OUH
COMMON /INOUT/ IUC,IUD,IUF,IUH,OUC,OUD,OUF,OUH
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO,ALL
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /TEMP/ATM(20),TRM(20),ZZ(20),IND(20),PN,Z2,FN,M,K,IDUM(3173)
COMMON /LABEL/CONFIG(3,40),EL(20),ATOM,TERM
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
LONG TITLE(3,3)
C
C ***** GENERATE ARRAYS FOR R,R*R AND SQRT(R) WITH A CONSTANT MESH
C ***** SIZE IN THE LOG(Z*R) VARIABLE
C
DO 1 I=1,NO
R(I)= EXP(RHO)/Z
RR(I) = R(I)*R(I)
R2(I) = SQRT(R(I))
1 RHO = RHO + H
RHO = RHO - NO*H
C
C ***** SET PARAMTERS FOR ELECTRONS AND INITIALIZE FUNCTIONS
C
DO 9 I = 1,NWF
MAX(I) = NO
PN = HNORM(N(I),L(I),Z-S(I))
IF (IND(I)) 7,8,9
C
C ***** READ DATA CARDS FOR ELECTRON I IF IND=-1 ON 'ELECTRON' CARD
C
7 READ(IUF,10) ATM(I),TRM(I),M,ZZ(I),E(I,I),EK(I),AZ(I),
1 (P(I,J),J=1,M)
10 FORMAT (6X,A6,6X,A6,11X,I6,F6.0/3(4X,E14.7)/(5(1PE14.7)))
C
C ***** SCALE RESULTS IF CARDS ARE FOR AN ATOM WITH A DIFFERENT Z
C
IF (Z .EQ. ZZ(I)) GO TO 20
Z2 = Z/ZZ(I)
DO 11 J=1,M
11 P(I,J)= P(I,J)*Z2
20 IF (M .EQ. NO) GO TO 4
MAX(I) = M
C
C ***** SET REMAINING VALUES IN THE RANGE = 0.
C
M = M +1
DO 13 J=M,NO
13 P(I,J) = D0
GO TO 4
C
C ***** DETERMINE ESTIMATES OF THE WAVE FUNCTIONS BY THE SCREENED
C ***** HYDROGENIC APPROXIMATION
C
8 DO 3 J=1,NO
P(I,J) = PN*HWF(N(I),L(I),Z-S(I),R(J))/R2(J)
3 CONTINUE
C
C ***** IF AZ(I) WAS NOT SPECIFIED ON INPUT DATA, COMPUTE
C
AZ(I) = PN*(D2*(Z - D5*S(I))/N(I))**(L(I) + 1)
C
C ***** ORTHOGONALIZE TO INNER FUNCTIONS
C
4 IF (I .EQ. 1 .OR. I .LT. IB) GO TO 5
IM = I - 1
DO 6 II =1,IM
IF (L(I) .NE. L(II) .OR. N(I) .EQ. N(II) ) GO TO 6
IF (.NOT. ORTHO .AND. A(I,II,1) .EQ. D0 ) GO TO 6
PN = QUADR(I,II,0)
PNN = SQRT(D1 - PN*PN)
IF (P(I,1) - PN*P(II,1) .LT. D0) PNN = -PNN
DO 25 J = 1,NO
25 P(I,J) =(P(I,J) - PN*P(II,J))/PNN
EK(I) = -D5*HL(I,I)
6 CONTINUE
C
C ***** COMPUTE Y0(I,I) AND STORE
C
5 CALL YKF(I,I,0)
DO 2 J= 1,NO
2 Y(I,J) = YK(J)
9 CONTINUE
WRITE(6,14)
14 FORMAT(/// 8X,18HINITIAL ESTIMATES //10X,2HNL,
1 4X,5HSIGMA,6X,5HE(NL),4X,6HEK(NL),4X,6HAZ(NL),4X,9HFUNCTIONS//)
C
C ***** COMPUTE ONE-ELECTRON ENERGY PARAMETERS IF THEY WERE NOT
C ***** SPECIFIED ON INPUT.
C
DO 15 I = 1,NWF
IF ( EK(I) .EQ. D0) EK(I) = -D5*HL(I,I)
C IF (E(I,I) .EQ. D0) E(I,I) = HL(I,I) - EKIN(I,I)
K = IND(I) + 2
IF (IND(I) .GE. 0 ) GO TO 17
TITLE(2,1) = ATM(I)
TITLE(3,1) = TRM(I)
17 WRITE(6,19) EL(I),S(I),E(I,I),EK(I),AZ(I)
19 FORMATO(9X,A3,F9.2,F11.3,F10.3,F10.3)
IF (K-1) 15,505,510
505 WRITE(6,500) ATM(I),TRM(I)
510 IF (K-2) 15,506,511
506 WRITE(6,501) ATM(I),TRM(I)
511 IF (K-3) 15,507,15
507 WRITE(6,502) ATM(I),TRM(I)
500 FORMAT('SCALED ',2A6)
501 FORMAT('SCREENED HYDROGENIC',2A6)
502 FORMAT('UNCHANGED ',2A6)
15 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C 3-39 X C H
C ------------------------------------------------------------------
C
C THIS ROUTINE COMPUTES FUNCTIONS ASSOCIATED WITH THE EXCHANGE
C FUNCTION FOR THE I'TH RADIAL EQUATION, INCLUDING CONTRIBUTIONS
C FROM THE INTERACTIONS. THE EXACT FORM OF THE FUNCTION
C DEPENDS ON THE VALUE OF IOPT.
C
C VALUE OF IOPT FUNCTION
C ------------- --------
C 1 SQRT(R) X(R)
C 2 X(R)/SQRT(R)
C 3 R SQRT(R) ( X(R) + SUM E P )
C IJ J
C
SUBROUTINE XCH(I,IOPT)
REAL QC
LOGICAL FAIL,RSCAN(20),OMIT,ORTHO,EZERO,ALL
COMMON /TEST/FAIL,RSCAN,ORTHO,OMIT,EZERO ,ALL
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
COMMON /WAVE/PDE(220),EK(20),E(20,20),ED,AZD,AZ(20),SUM(20),S(20),
1 DPM(20),ACC(20),METH(20),IPR
COMMON /COEFF/A(20,20,5),B(20,20,5),CA(4,4),CB(5,5,5)
COMMON /STATE/WT(40),QC(20,40),CFG(200),CR(300),KFG(200),NF,NG,NR,
1 NL,IFG(200),NCI(200),JFG(200),NCJ(200),KR(300),I1R(300),I2R(300)
2 ,NCRI(300),J1R(300),J2R(300),NCRJ(300),IO(300),JO(300),IQ(300)
3 ,CL(20),ILI(20),NCLI(20),ILJ(20),NCLJ(20),ILO(20),JLO(20),LQ(20)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
DO 1 J=1,NO
1 X(J) = D0
DO 2 J=1,NWF
DO 3 K=1,5
C = B(I,J,K)*D2
IF (ABS(C) .LE. 1.E-10) GO TO 3
CALL YKF(I,J,2*(K-1) + IABS(L(I) - L(J)))
DO 6 JJ =1,NO
6 X(JJ) = X(JJ)+ C*YK(JJ)*P(J,JJ)
3 CONTINUE
2 CONTINUE
IF (RSCAN(I)) GO TO 4
IF (NR .EQ. 0) GO TO 22
DO 10 II=1,NR
N1 = NCRI(II)
N2 = NCRJ(II)
C = WT(N1)*WT(N2)*CR(II)/SUM(I)
DO 11 I2=1,2
DO 12 I1=1,2
IF (I1R(II) .NE. I) GO TO 13
CC = C
IF (IQ(II) .NE. 0) CC = CC*QUADR(IO(II),JO(II),0)**IQ(II)
CALL YKF(I2R(II),J2R(II),KR(II))
M = J1R(II)
DO 14 J=1,NO
14 X(J) = X(J) +CC*P(M,J)*YK(J)
13 III = I1R(II)
I1R(II)= I2R(II)
I2R(II)= III
III = J1R(II)
J1R(II) = J2R(II)
12 J2R(II) = III
III = I1R(II)
I1R(II) = J1R(II)
J1R(II) = III
III = I2R(II)
I2R(II)= J2R(II)
11 J2R(II)= III
IJ1 = IO(II)
IJ2 = JO(II)
K = 1
21 IF (IJ1 .NE. I) GO TO 19
CC = C*IQ(II)*RK(I1R(II),I2R(II),J1R(II),J2R(II),KR(II))*
1 QUADR(IO(II),JO(II),0)**(IQ(II)-1)
DO 20 J = 1,NO
20 X(J) = X(J) + CC*R(J)*P(IJ2,J)
19 IJT = IJ1
IJ1 = IJ2
IJ2 = IJT
K = K + 1
IF ( K .EQ. 2 ) GO TO 21
10 CONTINUE
C
C ***** ADD CONTRIBUTION FORM THE L INTEGRALS
C
22 IF (NL .EQ. 0) GO TO 4
DO 30 II = 1,NL
N1 = NCLI(II)
N2 = NCLJ(II)
C = WT(N1)*WT(N2)*CL(II)/SUM(I)
DO 31 I1 = 1,2
IF (ILI(II) .NE. I) GO TO 32
CALL DIFF(ILJ(II))
CC = C
IF (LQ(II) .NE. 0) CC = CC*QUADR(ILO(II),JLO(II),0)**LQ(II)
DO 34 J = 1,NO
34 X(J) = X(J) + CC*YK(J)/R(J)
32 IF (ILO(II) .NE. I) GO TO 35
LJ = JLO(II)
CC = C*LQ(II)*HL(ILI(II),ILJ(II))
IF (LQ(II) .GT. 1) CC = CC*QUADR(ILO(II),JLO(II),0)**(LQ(II)-1)
DO 36 J = 1,NO
36 X(J) = X(J) + CC*R(J)*P(LJ,J)
35 III = ILI(II)
ILI(II) = ILJ(II)
ILJ(II) = III
III = ILO(II)
ILO(II) = JLO(II)
31 JLO(II) = III
30 CONTINUE
4 GO TO (15,16,17),IOPT
16 DO 18 J = 1,NO
18 X(J) = X(J)/R(J)
GO TO 15
17 DO 5 J =1,NO
5 X(J) = R(J)*X(J)
DO 7 J = 1,NWF
C = E(I,J)
IF (C .EQ. D0 .OR. (J .EQ. I)) GO TO 7
DO 9 JJ = 1,NO
9 X(JJ) = X(JJ) + C*P(J,JJ)*RR(JJ)
7 CONTINUE
C
C ***** CHECK IF EXCHANGE IS ZERO: IF SO, METHOD 2 SHOULD BE USED.
C
15 IF (METH(I) .EQ. 2) RETURN
IF ( ABS(X(1)) + ABS(X(2)) + ABS(X(3)) .EQ. D0 ) METH(I) = 2
RETURN
END
C
C ------------------------------------------------------------------
C 3-40 Y K F
C ------------------------------------------------------------------
C
C K
C STORES Y (I, J; R) IN THE ARRAY YK
C
C
SUBROUTINE YKF(I,J,K)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
CALL ZK(I,J,K)
A = EH**(K+1)
C = 2*K+1
A2 = A*A
H90 = C*H3/D30
A3 = A2*A*H90
AI = H90/A
AN = 114.E0*A*H90
A34 = 34.E0*H90
F1 = YK(NO)*EH**K
F2 = YK(NO)
F3 = YK(NO-1)
F4 = YK(ND)
DO 9 MM = 2,ND
M = NO -MM
F5 = YK(M-1)
YK(M) = YK(M+2)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)
F1 = F2
F2 = F3
F3 = F4
9 F4 = F5
YK(1) = YK(3)*A2+C*H3*(F4 + D4*A*F3 + A2*F2)
RETURN
END
C
C ------------------------------------------------------------------
C 3-41 Z K
C ------------------------------------------------------------------
C
C K
C STORES Z (I, J; R) IN THE ARRAY YK.
C
C
SUBROUTINE ZK(I,J,K)
COMMON R(220),RR(220),R2(220),P(20,220),Y(20,220),YK(220),YR(220),
1 X(220),N(20),L(20),MAX(20)
COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,NP,NCFG,IB,IC,ID
1 ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,NSCF
DEN = L(I) + L(J) + 3+ K
FACT = (D1/(L(I)+1) + D1/(L(J)+1))/(DEN + D1)
A = EH**K
A2 = A*A
H90 = H/90.E0
A3 = A2*A*H90
AI = H90/A
AN = 114.E0*A*H90
A34 = 34.E0*H90
F1 = RR(1)*P(I,1)*P(J,1)
F2 = RR(2)*P(I,2)*P(J,2)
F3 = RR(3)*P(I,3)*P(J,3)
F4 = RR(4)*P(I,4)*P(J,4)
YK(1) = F1*(D1 + Z*R(1)*FACT)/DEN
YK(2) = F2*(D1 + Z*R(2)*FACT)/DEN
YK(3) = YK(1)*A2 + H3*(F3 + D4*A*F2 + A2*F1)
DO 8 M = 5,NO
F5 = RR(M)*P(I,M)*P(J,M)
YK(M-1) = YK(M-3)*A2 + ( AN*F3 + A34*(F4+A2*F2)-F5*AI-F1*A3)
F1 = F2
F2 = F3
F3 = F4
8 F4 = F5
YK(NO) = A*YK(NO-1)
IF (IABS(I-J) + IABS(K) .NE. 0) GO TO 2
C
C ***** FOR Y0(I,I) SET THE LIMIT TO 1 AND REMOVE OSCILLATIONS
C ***** INTRODUCED BY THE USE OF SIMPSON'S RULE
C
M2 = (NO/2)*2 - 1
M1 = M2 - 1
C1 = D1 - YK(M1)
C2 = D1 - YK(M2)
DO 3 M = 1,M1,2
YK(M) = YK(M) + C1
3 YK(M+1) = YK(M+1) + C2
YK(NO) = D1
YK(NO-1) = D1
2 RETURN
END
C
C ------------------------------------------------------------------
C 4 T E S T D A T A
C ------------------------------------------------------------------
C
C
C
/ end of mchf77
▶EOF◀