Created: November 02, 1998 --- Last Updated: August 06, 1999
By Mark D. Pauli
C HARTREE-FOCK-SLATER SELF-CONSISTENT ATOMIC FIELD PROGRAM
C ORIGINALLY WRITTEN BY SHERWOOD SKILLMAN
C RCA LABORATORIES, PRINCETON, NEW JERSEY, SPRING 1961
C MODIFIED BY FRANK HERMAN, SUMMER 1961
C FURTHER MODIFIED BY RICHARD KORTUM, LOCKHEED RESEARCH
C LABORATORIES, PALO ALTO, CALIFORNIA, SUMMER 1962
C MODIFIED BY L F MATTHEISS FOR CCL MONITOR SYSTEM
C MODIFIED BY S.CRAMPIN AND D.K.SALDIN FOR MUFPOT USE
C MODIFIED BY M.PAULI FOR LKKR USE
C ----------------------------------------------------------------------
C ITERATION NUMBER,MEASURE OF SELF-CONSISTENCY,AND ATOMIC NUMBER Z
C ARE ALWAYS PRINTED ON-LINE.
C OUTPUT TAPE B-5 CONTAINS SELF-CONSISTENT ATOMIC POTENTIAL,
C ENERGY EIGENVALUES, AND RADIAL WAVE FUNCTIONS
C SUCCESSIVE SOLUTIONS SEPARATED BY END OF FILE
C Z = ATOMIC NUMBER. XION = IONICITY. ZZZ = XION+1
C IF KEY = 0 NORMALIZED NUMERICAL ATOMIC POTENTIAL IS READ IN
C FOR EVERY 4TH MESH POINT FROM 1 TO 437.
C IF KEY = 1 NUMERICAL ATOMIC POTENTIAL IS TO BE READ IN.
C IF KEY = 2 EXTRAPOLATE FOR STARTING POTENTIAL
C IPRATT IS THE NUMBER OF CONSECUTIVE ITERATIONS TO USE THE PRATT
C IMPROVEMENT SCHEME FOLLOWING EACH APPLICATION OF THE ARITHMETIC
C AVERAGE SCHEME.
C MAXIT = MAXIMUM NO. OF ITERATIONS UNLESS MAXIT IS READ AS 0 ,
C IN WHICH CASE MAXIT IS SET TO 20.
C IF KUT = 0, MODIFIED H-F-S POTENTIAL IS USED IN WAVE EQUATION.
C IF KUT = 1, IDEAL H-F-S POTENTIAL IS USED IN WAVE EQUATION.
C ----------------------------------------------------------------------
CHARACTER*4 Q001,Q002,Q003,Q004,Q005,Q006,Q007,Q008,
* REC1,REC2,REC3,REC4
CHARACTER*8 TITLE(9)
DIMENSION X(521),DENM(521),RSATOR(521),RSCORE(521)
COMMON XNUM(521),RSVALE(521),RSATOM(521),V(521),
1 SNLO(521),R(521),RU(521),RUEXCH(521),XI(521),XJ(521),
2 SNL(14,521),RUINL1(521),RUFNL1(521),RUINL2(521),RUFNL2(521)
3 ,RU2(521),RU3(521),NNLZ(24),WWNL(24),NKKK(24),EE(24),A(4,5)
EQUIVALENCE (RSCORE,XNUM),(RSVALE,DENM) ,(QQ,RSATOM)
LOGICAL BRANCH,NPTRAD,BRIEF,CHG
1 FORMAT(F8.5,9F7.5)
3 FORMAT(F4.0,2I4,F10.6)
5 FORMAT(6H WWW= ,F4.0,6H ZZZ= ,F4.0,6H Z= ,F4.0,11H NCORES= ,I4
1,11H NVALES= ,I4,11H NCSPVS= ,I4/25H CONTROL CARDS INCORRECT.)
7 FORMAT(I4,F10.6,F11.4)
8 FORMAT (I5,0PF7.1,1PE16.8,I8,0PF11.4,I9,0PF11.4//)
10 FORMAT (1PE15.7,1P4E14.7)
11 FORMAT (//1X,'NCARDS=',I4/1X,'MCARDS=',I4)
100 FORMAT (I4,2F8.6,4I4,3F14.9)
2002 FORMAT (I4,I7,4X,1P4E14.7/)
2012 FORMAT(1P5E14.7)
2010 FORMAT (1PE15.7,1P4E14.7,1X,'Z',I3,I4)
C2011 FORMAT (1PE15.7,57X,'Z',I3,I4)
2011 FORMAT(1PE14.7,57X,'Z',I3,I4)
3001 FORMAT(4(A4,4X),I4)
DATA Q001,Q002,Q003,Q004,Q005,Q006,Q007,Q008/4HCONT,4HPUNC
U ,4H POT,4H RAD,4H-BRF,4H-CHG,4HNRAD,4HHYDR/
OPEN(5,FILE='input.dat',STATUS='OLD')
OPEN(7,FILE='hsall.raw.out',STATUS='NEW')
OPEN(8,FILE='hsall.header.out',STATUS='NEW')
OPEN(9,FILE='hsinfo.out',STATUS='NEW')
OPEN(10,FILE='hswf.out',STATUS='NEW')
OPEN(11,FILE='hspot.out',STATUS='NEW')
2005 READ (5,3001) REC1,REC2,REC3,REC4,NA
IF(REC1.EQ.Q008) CALL HYDROG
IF(REC1.EQ.Q008) GO TO 2005
IF(NA.EQ.0) NA=6
IF(REC1.NE.Q001.AND.REC1.NE.Q005) GO TO 3000
BRIEF=.FALSE.
CHG=.FALSE.
NPTRAD=.FALSE.
IF(REC1.EQ.Q005) BRIEF = .TRUE.
IF(REC2.EQ.Q006) CHG = .TRUE.
IF(REC4.EQ.Q007) NPTRAD = .TRUE.
3002 WRITE (9,9000)
9000 FORMAT(1H1)
15 CONTINUE
DO 6833 I=1,521
X(I)=0.
RSCORE(I)=0.
RUINL1(I)=0.
RU2(I)=0.
RU(I)=0.
R(I)=0.
RSVALE(I)=0.
RUFNL1(I)=0.
RU3(I)=0.
XI(I)=0.
V(I)=0.
RSATOM(I)=0.
RUINL2(I)=0.
XNUM(I)=0.
XJ(I)=0.
RUEXCH(I)=0.
RUFNL2(I)=0.
DENM(I)=0.
SNLO(I)=0.
DO 6833 J=1,14
SNL(J,I)=0.
6833 CONTINUE
DO 6844 I=1,24
NNLZ(I)=0
WWNL(I)=0.
NKKK(I)=0
EE(I)=0.
6844 CONTINUE
DO 6845 I=1,4
DO 6845 J=1,5
A(I,J)=0.
6845 CONTINUE
C ----------------------------------------------------------------------
C Start of Original Herman-Skillman program
C ----------------------------------------------------------------------
C READ HEADING CARD.
READ (5,9001) TITLE
WRITE (9,9002) TITLE
WRITE(10,9003)"HERMAN-S"
WRITE(11,9003)"HERMAN-S MODIFIED POTENTIAL RV(R) ON 441 pt GRID,
1IN HARTREES IDEALLY = -ZU"
WRITE(11,9003)"441 GRID STARTS WITH dx=0.0025 AND DOUBLES AFTER
1EVERY 40 POINTS."
WRITE(10,9004)TITLE
WRITE(11,9004)TITLE
9001 FORMAT (10A8)
9002 FORMAT(1X,9A8///)
9003 FORMAT(A)
9004 FORMAT(9A8)
C ----------------------------------------------------------------------
C READ CONTROL CARDS AND INPUT POTENTIALS. CALCULATE TRIAL POTENTIAL
READ (5,100) KEY,TOL,THRESH,MESH,IPRATT,MAXIT,KUT,RADION,RATIO,
1ALPHA
WRITE (9,9005) KEY,TOL,THRESH,MESH,IPRATT,MAXIT,KUT,RADION,RATIO,
1 ALPHA
9005 FORMAT (6X,'KEY TOL THRESH MESH IPRATT MAXIT KUT',
1 ' RADION RATIO ALPHA'/
2 6X,I2,F9.4,F10.6,I6,2I7,I6,F14.7,2F12.7//)
IF(RADION) 8003, 8003, 8004
8003 ASSIGN 8503 TO NT
GO TO 970
8004 IF(RATIO)8011,8010,8012
8011 WRITE (9,9006)
9006 FORMAT(16H RATIO INCORRECT)
STOP
8010 RATIO=1.0
8012 ASSIGN 8512 TO NT
C ----------------------------------------------------------------------
970 NCARDS=90
IF(MAXIT) 104,104,105
104 MAXIT=20
105 CONTINUE
107 NBLOCK = (MESH)/40
MMXY=NBLOCK
C CONSTRUCT X MESH AND R MESH
I=1
X(I)=0.0
R(I)=0.0
DELTAX=0.0025E0
DO 25 J=1,NBLOCK
DO 24 K=1,40
I=I+1
24 X(I)=X(I-1)+DELTAX
25 DELTAX=DELTAX+DELTAX
C READ IN NORMALIZED ATOMIC POTENTIAL; FORM BASED ON IF KEY = 0,1,2
IF(KEY-1)1203,208,14
14 READ (5,10)(RU2(M),M=1,441)
READ (5,10)(RU3(M),M=1,441)
ZE2=-RU2(1)/2.0
ZE3=-RU3(1)/2.0
GO TO 205
1203 READ (5,1)(RU2(M),M=1,437,4)
GO TO 205
208 READ (5,10)(RU3(M),M=1,441)
ZE3=-RU3(1)/2.0
205 READ (5,3) Z,NCORES,NVALES,XION
206 NFILES =NFILES+1
IZ=Z
C NCSPVS=Number of CoreS Plus ValenceS
NCSPVS=NCORES+NVALES
WRITE(10,9007)Z,NCSPVS
WRITE(11,9007)Z,NCSPVS
9007 FORMAT(F9.4/I4)
C=0.88534138E0/Z**(1.E0/3.E0)
ZION=XION
TWOION=XION+XION
ZZZ=XION+1.
TWOZZZ = ZZZ+ZZZ
DO 224 I=2,MESH
224 R(I)=C*X(I)
210 READ (5,7)(NNLZ(I),WWNL(I),EE(I),I=1,NCSPVS)
WRITE (9,9008) (NNLZ(I),WWNL(I),EE(I),I=1,NCSPVS)
9008 FORMAT(I4,F10.6,F14.3)
WWW=0.0
DO 112 I=1,NCSPVS
112 WWW=WWW+WWNL(I)
IF(ABS(Z+1.0-WWW-ZZZ)-0.001E0) 114,113,113
113 WRITE (9,5)WWW,ZZZ,Z,NCORES,NVALES,NCSPVS
STOP
114 CONTINUE
IF (KEY-1) 204,207,16
C CONSTRUCT ATOMIC POTENIAL
16 IF(ABS(ZE3-ZE2-Z+ZE3)-0.001E0) 18,18,1818
1818 WRITE (9,1817)Z,ZE2,ZE3
1817 FORMAT (27H0STARTING POTENTIALS AND Z=,F6.2,9H IN ERROR,2F6.2)
STOP
18 DO 21 I=1,441
RU(I) = RU3(I)+RU3(I)-RU2(I)
21 CONTINUE
GO TO 211
204 TWOZ=Z+Z
DO 110 I=1,437,4
110 RU(I)=-RU2(I)*TWOZ
RU(441)=RU(437)
RU(445)=RU(437)
M=9
DO111 I=1,437,4
M=M-1
IF(M) 1100,1101,1101
1100 RU(I+1)=(22.0*RU(I)+11.0*RU(I+4)-RU(I+8))/32.0
RU(I+2)=(10.0*RU(I)+15.0*RU(I+4)-RU(I+8))/24.0
RU(I+3)=( 6.0*RU(I)+27.0*RU(I+4)-RU(I+8))/32.0
M=9
GO TO 111
1101 RU(I+1)=(21.0*RU(I)+14.0*RU(I+4)-3.0*RU(I+8))/32.0
RU(I+2)=( 3.0*RU(I)+ 6.0*RU(I+4)- RU(I+8))/ 8.0
RU(I+3)=( 5.0*RU(I)+30.0*RU(I+4)-3.0*RU(I+8))/32.0
111 CONTINUE
GO TO 211
207 IF (ABS(ZE3-Z)-0.001E0) 215,215,220
215 DO 213 I=1,441
RU(I)=RU3(I)
213 CONTINUE
GO TO 211
220 ZOZ=Z/ZE3
DO 223 I=1,441
RU(I)=RU3(I)*ZOZ
223 CONTINUE
211 V(1) = -9.9D35
M=MIN0(441,MESH)
IF(KUT) 1212,1260,1212
1212 DO 1213 I=2,M
1213 V(I)= RU(I)/R(I)
C IF MESH > M, THEN FILL IN POTENTIAL USING THEORETICAL LIMIT VALUES
IF(MESH-M)1216,1216,1214
1214 DO 1215 I=442, MESH
1215 V(I)=-TWOION/R(I)
1216 LIMIT=M
ICUT= MESH
IC=MESH
GO TO 268
1260 CONTINUE
ICUT=0
DO 26 I=2,M
IF (ICUT) 261,261,263
261 IF( TWOZZZ+RU(I)) 264,264,262
262 ICUT =I
263 V(I)=-TWOZZZ/R(I)
GO TO 26
264 V(I)=RU(I)/R(I)
26 CONTINUE
IF(ICUT) 265,265,266
265 ICUT=M
266 LIMIT=ICUT
IF(MESH-M) 268,268,267
267 CONTINUE
DO 115 I=442,MESH
115 V(I)=-TWOZZZ/R(I)
268 CONTINUE
DELTA =1000000.0
NITER=0
NONMON =3
IPRSW=0
C START ITERATION
27 MCARDS=90
IF(MAXIT-NITER) 405,28,28
405 CONTINUE
WRITE (9,450)
DO 410 I=1,MESH,5
WRITE (9,451)I,X (I),RU3(I),RUINL1(I),RUFNL1(I), RUINL2(I),R
1UFNL2(I),RU(I)
410 CONTINUE
450 FORMAT ('1'//8X,'I',6X,'X',11X,'RUINL1',10X,'RUFNL1',10X,
1 'RUINL2',10X,'RUFNL2',12X,'RU'/)
451 FORMAT (1X,I8,F12.7,1P6E16.7)
GO TO 120
28 DO 29 I=1,MESH
RSCORE(I)=0.0
29 RSVALE(I)=0.0
C ----------------------------------------------------------------------
C SOLVE SCHROEDINGER EQUATION FOR EACH STATE IN TURN. ALSO CALCULATE
C CORE AND VALENCE CHARGE DENSITY.
DO 35 M=1,NCSPVS
C E = inputted Binding Energy of the quantum level nl.
E=EE(M)
C Truncate by division to get NN = n of nlm from NNLZ
NN=NNLZ(M)/100
C Truncate by division and use n to get LAM = l of nlm
LAM =NNLZ(M)/10 -10*NN
XL= LAM
CALL SCHEQ(Z,E,LAM,NN,KKK,MESH,C,THRESH)
IF(M-NCORES)30,30,32
30 DO 31 I=1,KKK
31 RSCORE(I)=RSCORE(I)+WWNL(M)*SNLO(I)**2
GO TO 34
32 DO 33 I=1,KKK
33 RSVALE(I)=RSVALE(I)+WWNL(M)*SNLO(I)**2
34 DO 201 I=1,KKK
201 SNL(M,I)=SNLO(I)
NKKK(M)=KKK
MCARDS=MCARDS+2+((KKK-1)/40)*8
35 EE(M)=E
C ----------------------------------------------------------------------
C CALCULATE TOTAL CHARGE DENSITY AND ATOMIC EXCHANGE POTENTIAL
DO 37 I=1,MESH
36 RSATOM(I)=RSCORE(I)+RSVALE(I)
37 RUEXCH(I)=-6.0*ALPHA*((3.0*R(I)*RSATOM(I))/315.82734E0)**
1 (1.E0/3.E0)
C ----------------------------------------------------------------------
C CALCULATE ATOMIC COULOMB POTENTIAL
364 RSATOR(1)=0.0
DO 39 I=2,MESH
39 RSATOR(I)=RSATOM(I)/R(I)
CALL INTEGR(RSATOM,XI,R,MMXY)
CALL INTEGR(RSATOR,XJ,R,MMXY)
40 CONTINUE
41 DO 44 I=1,MESH
XI(I)=-2.0*Z+2.0*(XI(I)+R(I)*(XJ(MESH)-XJ(I)))
C SET RADIUS OF CHARGED SPHERE EQUAL TO IONIC RADIUS. THIS IS READ
C IN AS SPECIFIED ON THE DATA CARD REFERED BY PROGRAM CARD 065.
C IF RADION IS 0.0 OR BLANK, THIS ROUTINE IS DELETED FROM MAIN H-S
C PROGRAM.
GO TO NT,(8503,8512)
8512 XION=ZION*RATIO
IF(R(I)-RADION) 8001, 8002, 8002
8001 XI(I)=XI(I)+2.0*XION/RADION*R(I)
GO TO 8503
8002 XI(I)=XI(I)+2.0*XION
C ----------------------------------------------------------------------
8503 XJ(I)=XI(I)+RUEXCH(I)
44 CONTINUE
74 DO 75 I=1,MESH
RUINL1(I)=RUINL2(I)
RUFNL1(I)=RUFNL2(I)
RUINL2(I)=RU(I)
75 RUFNL2(I)=XJ(I)
76 NITER=NITER+1
PDELTA =DELTA
DELTA =0.0
DO 48 I=1,LIMIT
45 SNLO(I)=RU(I)-XJ(I)
46 XI(I)=ABS(SNLO(I))
IF (XI(I)-DELTA)48,48,47
47 DELTA=XI(I)
IDELTA=I
48 CONTINUE
WRITE (9,908)
908 FORMAT (/2X,'ITER',4X,'Z',7X,'DELTA',8X,'I(DEL) X(DEL) ',
1 ' I(CUT) X(CUT)')
WRITE (9,8)NITER,Z,DELTA,IDELTA,X(IDELTA),ICUT,X(ICUT)
C TEST SELF-CONSISTENCY OF ATOMIC POTENTIAL.
50 IF(DELTA-TOL) 58,51,51
C IF SCF CRITERION NOT SATISFIED, CALCULATE NEXT TRIAL POTENTIAL.
51 IF(IPRSW) 80,80,82
80 DO 81 I=2,LIMIT
81 RU(I)=0.5E0*(RU(I)+XJ(I))
IF (MESH-LIMIT) 280,280,270
270 RUZM=XJ(MESH)
RATIO =(RUZM-RU(LIMIT))/(RUZM-XJ(LIMIT))
DO 271 I=LIMIT,MESH
271 RU(I)=RUZM-RATIO*(RUZM-XJ(I))
LIMIT =MESH
280 IPRSW=IPRATT
GO TO 87
82 CONTINUE
282 IF(NONMON) 80,80,183
183 IF(PDELTA-DELTA)184,184,185
C IF DELTA IS NOT MONOTONIC DECREASING FOUR TIMES, BYPASS
C PRATT IMPROVEMENT SCHEME
184 NONMON =NONMON-1
IF(NONMON)80,80,185
185 ALPH=0.5E0
C ----------------------------------------------------------------------
C PRATT IMPROVEMENT SCHEME
DO 849 I=2,ICUT
XNUM(I)=RUINL1(I)*RUFNL2(I)-RUINL2(I)*RUFNL1(I)
DENM(I)=RUFNL2(I)-RUFNL1(I)-RUINL2(I)+RUINL1(I)
182 IF(ABS(DENM(I)/RUINL2(I))-0.0001) 83,83,84
83 CONTINUE
832 ALPH=0.5E0
GO TO 846
84 ALPH=(XNUM(I)/DENM(I)-RUFNL2(I))/SNLO(I)
IF(ALPH) 841,845,842
841 ALPH=0.0
GO TO 846
842 IF (0.5E0-ALPH) 843,845,845
843 ALPH=0.5E0
845 CONTINUE
846 XI(I)=ALPH
849 CONTINUE
IPRSW=IPRSW-1
IF(KUT) 1873,1849,1873
1849 CONTINUE
IC=ICUT+20
IC1=ICUT+1
ADEL=XI(ICUT)/20.0
DO 873 I=IC1,IC
XI(I)=XI(I-1)-ADEL
873 XJ(I)=XI(I)
1873 CONTINUE
XJ(1)=0.5E0
XJ(2)=XI(2)
ASUM=XI(2)+XI(3)+XI(4)+XI(5)
DO 874 I=3,ICUT
XJ(I)=ASUM*0.2E0
874 ASUM=ASUM-XI(I-2)+XI(I+3)
IF(KUT) 1875,1874,1875
1874 CONTINUE
IC1=IC+1
DO 875 I=IC1,MESH
XJ(I)=0.0
875 RU(I)=RUFNL2(I)
1875 CONTINUE
DO 876 I=2,IC
876 RU(I)= RUFNL2(I)+XJ(I)*SNLO(I)
87 CONTINUE
IF(KUT) 2870,1870,2870
2870 ICUT=MESH
LIMIT=MESH
DO 2871 I=2,MESH
VLAST= V(I)
V(I)=RU(I)/R(I)
2871 XI(I)=V(I)-VLAST
GO TO 1861
1870 CONTINUE
ICUT =0
DO 86 I=2,MESH
VLAST=V(I)
IF(ICUT) 851,851,853
851 IF(TWOZZZ+RU(I)) 854,854,852
852 ICUT=I
853 V(I)=-TWOZZZ/R(I)
GO TO 86
854 V(I)=RU(I)/R(I)
86 XI(I)=V(I)-VLAST
1861 CONTINUE
XI(1)=0.0
C ----------------------------------------------------------------------
C NEXT TRIAL EIGENVALUES PREDICTED BY PERTURBATION THEORY
NCARDS=90
DO 55 M=1,NCSPVS
K=(NKKK(M)-1)/40
H=0.0025E0*C
ASUM=0.0
A1=0.0
I=1
DO 54 J=1,K
DO 53 L=1,40
I=I+1
A2=XI(I)*SNL(M,I)**2
53 A1=A1+A2*H
ASUM=ASUM+A1-(A2/2.0)*H
H=H+H
54 A1=(A2/2.0)*H
EE(M)=EE(M)+ASUM
55 NCARDS=NCARDS+8*K+2
WRITE (9,9008) (NNLZ(I),WWNL(I),EE(I),I=1,NCSPVS)
GO TO 27
58 CONTINUE
580 WRITE (9,11)NCARDS,MCARDS
WRITE (9,9000)
IF(REC2.NE.Q002) GO TO 1901
1900 WRITE (7,3) Z,NCORES,NVALES,XION
1901 IF(BRIEF) GO TO 586
WRITE (8,9009)
9009 FORMAT (' Z CORE VAL ION')
586 CONTINUE
WRITE (8,3) Z,NCORES,NVALES,XION
IF(KUT) 5810,5810,584
5810 CONTINUE
DO 583 I=1,441
IF(TWOZZZ+RUINL2(I)) 583,583,581
581 DO 582 M=I,441
582 RUINL2(M)=-TWOZZZ
GO TO 584
583 CONTINUE
584 CONTINUE
IF(BRIEF) GO TO 585
WRITE (8,9010)
9010 FORMAT (/// 109HTHE FOLLOWING DATA IS THE MODIFIED HFS POTENTIAL,
1RV(R), GIVEN AT THE 441 POINTS OF THE MESH IDEALLY = -2Z*U/)
585 CONTINUE
BRANCH = REC3.NE.Q003
NC=0
DO 588 MIN=1,440,5
MAX= MIN+4
NC=NC+1
IF(BRANCH) GO TO 3100
1902 WRITE (7,2010) (RUINL2(M),M=MIN,MAX),IZ,NC
WRITE (11,2015)(RUINL2(M)/2.0,M=MIN,MAX)
2015 FORMAT (1P5E14.7)
3100 IF(BRIEF) GO TO 588
1903 WRITE (8,2010)(RUINL2(M),M=MIN,MAX),IZ,NC
588 CONTINUE
NC=NC+1
IF(BRANCH) GO TO 3101
1904 WRITE (7,2011) RUINL2(441),IZ,NC
WRITE (11,2016) RUINL2(441)/2.0
2016 FORMAT (1PE14.7)
3101 IF(BRIEF) GO TO 3102
1905 WRITE (8,2011)RUINL2(441),IZ,NC
3102 BRANCH=REC4.NE.Q004
WRITE (8,9000)
C ----------------------------------------------------------------------
C LOOP OVER THE ORBITALS
DO 59 M=1,NCSPVS
NLZ=NNLZ(M)
KKK=NKKK(M)
XL=NLZ/10-10*(NLZ/100)
LP=XL+1.0
C COMPUTE FIRST TERM OF SERIES. (SNL(R)/R**(LAM+1) AT R=0)
DO 591 I=1,4
A(I,1)=1.0
A(I,2)=R(I+1)
A(I,3)=R(I+1)*R(I+1)
A(I,4)=R(I+1)*A(I,3)
591 A(I,5)=SNL (M,I+1)/R(I+1)**LP
CALL CROSYM (4)
IF (NPTRAD) GO TO 4100
IF(BRIEF) GO TO 3203
WRITE (8,9300) NLZ, KKK
WRITE(10,9126)INT(XL),KKK,WWNL(M)/(4*XL+2)
9126 FORMAT(I4/I4/F9.4)
9300 FORMAT (/// 71H THE FOLLOWING DATA IS THE NORMALIZED RADIAL WAVE
1FUNCTION FOR ORBITAL , I3, 21H , OUT TO MESH POINT ,I3/)
3203 IF(BRANCH) GO TO 3103
1906 WRITE (7,2002) NLZ,KKK,XL,EE(M),WWNL(M),A(1,5)
3103 IF(BRIEF) GO TO 3104
1907 WRITE (8,2002) NLZ,KKK,XL,EE(M),WWNL(M),A(1,5)
3104 K1=KKK-1
NC=0
DO598 MIN=1,K1 ,5
MAX= MIN+4
NC=NC+1
IF(BRANCH) GO TO 3105
WRITE (10,2013) (SNL(M,I),I=MIN,MAX)
1908 WRITE (7,2010) (SNL(M,I),I=MIN,MAX),IZ,NC
3105 IF(BRIEF) GO TO 598
1909 WRITE (8,2010)(SNL(M,I),I=MIN,MAX),IZ,NC
2013 FORMAT(1PE14.7,1P4E14.7)
598 CONTINUE
NC=NC+1
IF(BRANCH) GO TO 3106
WRITE (10,2011) SNL(M,KKK),IZ,NC
1910 WRITE (7,2011) SNL(M,KKK),IZ,NC
3106 IF(BRIEF) GO TO 59
1911 WRITE (8,2011)SNL(M,KKK),IZ,NC
59 CONTINUE
C ----------------------------------------------------------------------
4100 BRANCH=REC2 .NE. Q002 .AND. (.NOT.CHG)
IF(BRIEF) GO TO 3204
WRITE (8,9500)
9500 FORMAT ('1'///1X,'TOTAL CHARGE DENSITY OUTPUT'//)
3204 CONTINUE
NC=0
IF(REC2.NE.Q006) GO TO 3205
WRITE(7,3202) IZ
3202 FORMAT(I4,4H 1,4H 0,2(4H 1),4H 0)
3205 CONTINUE
DO 590 MIN=1,436,5
MAX=MIN+4
NC=NC+1
IF(BRANCH) GO TO 3107
WRITE (7,4001) (RSATOM(I),I=MIN,MAX),IZ,NC
3107 IF(BRIEF) GO TO 590
WRITE (8,4001) (RSATOM(I),I=MIN,MAX),IZ,NC
590 CONTINUE
4001 FORMAT (1PE15.7,1P4E14.7,1X,'T',I3,I4)
NC=NC+1
IF(BRANCH) GO TO 3108
WRITE (7,4002) RSATOM(441),IZ,NC
3108 IF(BRIEF) GO TO 3200
WRITE (8,4002) RSATOM(441),IZ,NC
3200 CONTINUE
4002 FORMAT (1PE15.7,57X,'T',I3,I4)
IF(NCORES.EQ.0) GO TO 956
IF(BRIEF) GO TO 3207
WRITE (8,9700)
9700 FORMAT ('1'///1X,'CORE CHARGE DENSITY OUTPUT'//)
3207 NC=0
DO 595 MIN=1,436,5
MAX=MIN+4
NC=NC+1
IF(BRANCH.OR.CHG) GO TO 3109
WRITE (7,4003) (RSCORE(I),I=MIN,MAX),IZ,NC
3109 IF(BRIEF) GO TO 595
596 WRITE (8,4003) (RSCORE(I),I=MIN,MAX),IZ,NC
595 CONTINUE
4003 FORMAT (1PE15.7,1P4E14.7,1X,'C',I3,I4)
NC=NC+1
IF(BRANCH.OR.CHG) GO TO 3110
WRITE (7,4004) RSCORE(441),IZ,NC
3110 IF(BRIEF) GO TO 3201
594 WRITE (8,4004) RSCORE(441),IZ,NC
3201 CONTINUE
4004 FORMAT (1PE15.7,57X,'C',I3,I4)
C ----------------------------------------------------------------------
956 CONTINUE
IF(KEY-1) 958,960,965
958 KEY=1
960 DO 961 I=1,MESH
961 RU3(I)=RU(I)
ZE3=Z
GO TO 120
965 DO 966 I=1,MESH
RU2(I)=RU3(I)
966 RU3(I)=RU(I)
ZE2=ZE3
ZE3=Z
120 GO TO 2005
3000 STOP
END
C ----------------------------------------------------------------------
C ----------------------------------------------------------------------
SUBROUTINE SCHEQ(ZZ,EN,LAMBDA,NOFL,KKK,MESS,SCF,THRESH)
C SCHEQ: SCHROEDINGER EQUATION SUBROUTINE
C COMPUTE ENERGY EIGENVALUE AND WAVE FUNCTION
C ORIGINALLY WRITTEN BY SHERWOOD SKILLMAN
C RCA LABORATORIES, PRINCETON, NEW JERSEY, SPRING 1961
C MODIFIED BY FRANK HERMAN, SUMMER 1961
C FURTHER MODIFIED BY RICHARD KORTUM AND PAUL KELLY, LOCKHEED
C RESEARCH LABORATORIES, PALO ALTO, CALIFORNIA, SUMMER 1962
C --------------------------------------------------------
C R(521) - H-S mesh
C ---------------------------------------------------------
DIMENSION DENM(521),QQ(521),P(5),Q(5),T(5),D(5),RSCORE(521)
COMMON XNUM(521),RSVALE(521),RSATOM(521),V(521),
1 SNLO(521),R(521),RU(521),RUEXCH(521),XI(521),XJ(521),
2 SNL(14,521),RUINL1(521),RUFNL1(521),RUINL2(521),RUFNL2(521)
3 ,RU2(521),RU3(521),NNLZ(24),WWNL(24),NKKK(24),EE(24),A(4,5)
EQUIVALENCE (RSCORE,XNUM),(RSVALE,DENM) ,(QQ,RSATOM)
C SET UP CONSTANTS AND INITIALIZE
Z=ZZ
DO 7643 JK=1,5
P(JK)=0.0
T(JK)=0.0
D(JK)=0.0
Q(JK)=0.0
7643 CONTINUE
C LAMBDA = l and NN = n of nlm.
LAM=LAMBDA
NN=NOFL
MESH=MESS
C C is the radial mesh constant; x=r/u u=C/Z^1/3
C=SCF
MANY = 200
73 E=EN
EG=0.0
MOREV=0
LESSV=0
EMORE=0.0
ELESS=0.0
MORE=0
LESS=0
DE=0.0
NPRINT=0
LAMM=LAM-1
LAMP=LAM+1
XLP=LAMP
NDCR=NN-LAMP
C B = l(l+1)
B=LAM*LAMP
OC=R(2)
H=OC
HSQ=H*H
B3=(V(3)-V(2))/H-Z/HSQ
Y=H+H
FLPS=4*LAM+6
SLPT=6*LAM+12
ELPT=8*LAM+20
A1=-Z/XLP
YSQ=Y*Y
B1=-Z-Z
AB1=A1*B1
AB3=A1*B3
C RAISE H AND Y TO LAM+1
HTL=H
YTL=Y
IF(LAM)77,1102,1100
1100 DO 1101 I=1,LAM
HTL=HTL*H
1101 YTL=YTL*Y
1102 H1=HSQ
BOHS=B/HSQ
BOH=B1/H
BTH=B3*H
BQ3=BOHS+BOH+BTH
BQ4=BOHS/4.0+BOH/2.0+BTH+BTH
EPL=8+LAM
FPL=5+LAM**2
XIFC=C*0.2170139E-4
C START OUTWARD INTEGRATION
10 NPRINT=NPRINT+1
EPS =E-EG
EG =E
IF(MANY-NPRINT) 900,75,75
900 WRITE (9,901)NN,LAM ,Z
901 FORMAT (' NO CONVERGENCE ON NN=',I4,5X,'LAM=',I2,5X,'Z=',F6.2)
CALL DUMP
77 NSTOP=77
777 WRITE (9,770)NSTOP,E
WRITE (*,*)'Humbug, we have an problem'
770 FORMAT('!!!ERROR at NSTOP=',I4,' IN SUBROUTINE SCHEQ',1P1E16.8)
CALL DUMP
75 DO 11 I=1,MESH
11 SNLO(I)=0.0
IF(NPRINT-1) 77,500,529
500 CONTINUE
DO 502 I= 4,MESH
QQ(I) = V(I)+B/(R(I)*R(I))-E
502 CONTINUE
516 M= MESH
DO 520 I=4,MESH
IF(QQ(M)) 519,520,520
519 IK=M+1
GO TO 525
520 M=M-1
521 NSTOP =521
C Q(charge density) IS EVERYWHERE POSITIVE
GO TO 777
525 IF(MESH-IK) 526,526,535
526 EPS = QQ(MESH-40)
E = E+EPS
529 CONTINUE
DO 530 I=4,MESH
530 QQ(I) = QQ(I)-EPS
GO TO 516
535 CONTINUE
14 NCROSS=0
SIGN=1.0
H=OC
Y=H+H
C B= LAM*(LAM+1)
C B1= -2.*Z
B2=3.0*Z/H-E+2.0*V(2)-V(3)
C B3=(V(3)-V(2))/H -Z/HSQ
C A1= -Z/(LAM+1)
A2=(AB1+B2)/FLPS
C A2=(A1*B1+B2)/(4*LAM+6)
A3=(A2*B1+A1*B2+B3)/SLPT
C A3=(A2*B1+A1*B2+B3)/(6*LAM+12)
A4=(A3*B1+A2*B2+AB3)/ELPT
C A4=(A3*B1+A2*B2+A1*B3)/(8*LAM+20)
P(3)=(1.0+H*(A1+H*(A2+H*(A3+H*A4))))*HTL
C P(3)=(1.0+A1*H+A2*H**2+A3*H**3+A4*H**4)*H**(XL+1.0)
P(4)=(1.0+Y*(A1+Y*(A2+Y*(A3+Y*A4))))*YTL
C P(4)=(1.0+A1*Y+A2*Y**2+A3*Y**3+A4*Y**4)*Y**(XL+1.0)
Q(3)=BQ3+B2
C Q(3)=(B+B1*H+B2*H**2+B3*H**3)/H**2
Q(4)=BQ4+B2
C Q(4)=(B+B1*Y+B2*Y**2+B3*Y**3)/Y**2
SNLO(2)=P(3)
SNLO(3)=P(4)
I=3
DX=OC
H1=H**2
H2=H1/12.0
T(3)=P(3)*(1.0-H2*Q(3))
T(4)=P(4)*(1.0-H2*Q(4))
D(4)=T(4)-T(3)
NCOUNT=3
NINT=2
15 I=I+1
C IF END OF MESH IS REACHED, MODIFY TRIAL EIGENVALUE
IF(I-MESH) 16,3,3
3 IF(NDCR-NCROSS) 32,33,33
C RETURN TO BEGINNING OF OUTWARD INTEGRATION IF NECESSARY
16 Q(5) =QQ(I)
IF(IK-I) 29,29,21
21 D(5)=D(4)+H1*Q(4)*P(4)
T(5)=D(5)+T(4)
IF(1.0-ABS(H2*Q(5))) 3,3,501
501 P(5)=T(5)/(1.0-H2*Q(5))
SNLO(I) = P(5)
IF(SIGN) 211,77,212
211 IF(P(5)) 23,23,22
212 IF(P(5)) 22,23,23
22 NCROSS=NCROSS+1
C COUNT CHANGES IN SIGN
SIGN=-SIGN
23 NCOUNT=NCOUNT+1
IF(7-NCOUNT)77,24,25
24 NCOUNT=2
25 NINT=NINT+1
IF(40-NINT)77,26,27
26 DX=DX+DX
H=DX
H1=H**2
H2=H1/12.0
NINT=0
T(5)=P(5)*(1.0-H2*Q(5))
T(3)=P(3)*(1.0-H2*Q(3))
D(5)=T(5)-T(3)
27 DO 28 K=1,4
P(K)=P(K+1)
T(K)=T(K+1)
D(K)=D(K+1)
28 Q(K)=Q(K+1)
GO TO 15
29 IF(NCOUNT-2)77,30,21
30 IF(NINT-4)21,21,31
C MATCHING RADIUS HAS BEEN REACHED GOING OUT
C IF NDCR NOT EQUAL TO NCROSS, MODIFY TRIAL EIGENVALUE
31 EIGEN=E
IF(NDCR-NCROSS) 32,35,33
32 MORE=1
C TOO MANY CROSSINGS, INCREASE ABSF(E)
MOREV=MOREV+1
IF(MOREV-1) 50,53,52
50 NSTOP=50
GO TO 777
52 IF(E -EMORE) 53,54,54
53 EMORE=E
54 IF (LESS) 55,56,64
55 NSTOP=55
GO TO 777
56 E=1.25E0*EG
GO TO 10
33 LESS=1
C TOO FEW CROSSINGS, DECREASE ABSF(E)
LESSV=LESSV+1
IF(LESSV-1) 57,60,59
57 NSTOP=57
GO TO 777
59 IF(ELESS- E) 60,61,61
60 ELESS=E
61 IF(MORE) 62,63,64
62 NSTOP=62
GO TO 777
63 E=0.75E0*EG
GO TO 10
64 E=0.5E0*(EMORE+ELESS)
GO TO 10
35 IF(ABS(SNLO(I-1))-ABS(SNLO(I-2))) 351,354,354
C CHECK TO SEE THAT WAVE IS IN THE DAMPED REGION (ABSOLUTE VALUE
C DECREASING AND SIGNS ALIKE)
351 IF (P(5)) 352,21 ,353
352 IF(SNLO(I-2)) 401,21,21
353 IF(SNLO(I-2)) 21,21,401
C LARGE ABSOLUTE VALUE OF P IN WHAT SHOULD BE THE DAMPED REGION
354 IF (1.0D25-ABS(P(5))) 33,33,21
C INDICATES TOO FEW PEAKS, DECREASE ABSF(E)
C NOW NDCR = NCROSS AND MATCHING RADIUS LIES IN DAMPED REGION
401 IMATCH=I-2
XMATCH=R(I-2)
PPOUT=(T(4)-T(2)-0.5E0*(P(4)-P(2)))/H
S2=PPOUT/P(3)
C INTEGRATION IS BY 8 APPLICATIONS OF NEWTON-COTES CLOSED
C QUADRATURE FOR FIVE INTERVALS ON EACH BLOCK
C XIFC =(5*H(BLOCK 1)/288)/2 ,H(1) =0.0025*SCALE FACTOR
SUM1=0.0
XIF=XIFC
I=1
VALUE=0.0
36 MM=8
SUM2=0.0
XIF=XIF+XIF
37 Y=VALUE
VALUE=SNLO(I+5)**2
SUM2=SUM2+19.0*(VALUE+Y)+75.0*(SNLO(I+4)**2+SNLO(I+1)**2)
1+50.0*(SNLO(I+2)**2+SNLO(I+3)**2)
I=I+5
IF (IMATCH-I) 77,39,371
371 MM=MM-1
IF(MM)77,38,37
38 SUM1=SUM2*XIF+SUM1
GO TO 36
39 SUM1= SUM1+SUM2*XIF
40 S1=SUM1/P(3)**2
PMATCH=P(3)
IF(NN-1)77,41,42
41 XINW=EPL*XMATCH
C FOR N =1, START INWARD INTEGRATION AT(8+LAM)*XMATCH OR X MAX
GO TO 421
42 XINW=FPL*XMATCH
C FOR N NOT=1, START AT (5+LAM)*XMATCH OR X MAX (END OF MESH)
421 DO 44 I= 41,MESH,40
IF(XINW-R(I)) 43,43,44
43 KKK =I
GO TO 45
44 CONTINUE
KKK =MESH
45 I =KKK
DX =R(I-1)-R(I)
H =DX
XIF=0.1736111E-1*DX
HSQ=H*H
HSQ12=HSQ/12.0
Q(3)= QQ(I)
P(3)= EXP(-R(I)*SQRT(Q(3)))
402 SUM3=P(3)/Q(3)
I=I-1
Q(4) =QQ(I)
404 P(4)= EXP(-R(I)*SQRT(Q(4)))
IF (ABS(P(4))-1.0E-35) 4041,4041,405
4041 KKK=KKK -40
IF(KKK -IMATCH) 4042,4042,45
4042 WRITE (9,4043)Z ,NN,LAM,KKK
4043 FORMAT (6H0AT Z=,F6.0,6H NL =,I3,I1,7H KKK =,I5,22H IS LESS THA
1N IMATCH =,I5,43H INWARD INTEGRATION WILL BE TRIED AT KKK+40)
KKK =KKK +40
P(4)=1.5E-35
P(3)=1.0E-35
405 IF(PMATCH)102,77,103
102 P(3)=-P(3)
P(4)=-P(4)
103 SNLO(I+1)=P(3)
SNLO(I)=P(4)
T(3)=P(3)*(1.0-HSQ12 *Q(3))
T(4)=P(4)*(1.0-HSQ12*Q(4))
D(4)=T(4)-T(3)
104 DO 106 M=2,40
I=I-1
Q(5) =QQ(I)
D(5)=HSQ*Q(4)*P(4)+D(4)
T(5)=D(5)+T(4)
P(5)=T(5)/(1.0-HSQ12*Q(5))
IF(I-IMATCH+1)77,200,105
105 SNLO(I)=P(5)
DO 106 K=1,4
P(K)=P(K+1)
T(K)=T(K+1)
D(K)=D(K+1)
106 Q(K)=Q(K+1)
Q(5) =QQ(I-2)
D(5)=HSQ*Q(4)*P(4)+D(4)
T(5)=D(5)+T(4)
P(5)=T(5)/(1.0-HSQ12*Q(5))
P(5)=1.09375E0*P(4)+0.2734375E0*P(5)-0.546875E0*P(3)
1 +0.21875E0*P(2)-0.0390625E0*P(1)
I=I-1
DX=DX/2.0
Q(5) =QQ(I)
H=DX
HSQ=H*H
HSQ12=HSQ/12.0
T(5)=P(5)*(1.0-HSQ12*Q(5))
T(4)=P(4)*(1.0-HSQ12*Q(4))
D(5)=T(5)-T(4)
SNLO(I)=P(5)
DO 107 L=1,4
P(L)=P(L+1)
T(L)=T(L+1)
D(L)=D(L+1)
107 Q(L)=Q(L+1)
GO TO 104
C MATCHING RADIUS HAS BEEN REACHED COMING IN
200 K=KKK
VALUE=SNLO(K)**2
GO TO 202
2001 CONTINUE
201 SUM3=SUM3+XIF*SUM4
XIF=XIF*0.5E0
202 MM=8
SUM4 =0.0
203 Y=VALUE
VALUE=SNLO(K-5)**2
SUM4=SUM4+19.0*(VALUE+Y)+75.0*(SNLO(K-1)**2+SNLO(K-4)**2)
1+50.0*(SNLO(K-2)**2+SNLO(K-3)**2)
K=K-5
IF(K-IMATCH) 77,2031,2030
2030 MM=MM-1
IF(MM) 77,2001,203
2031 SUM3=SUM3+XIF*SUM4
204 S3=SUM3/P(4)**2
PPIN=(T(5)-T(3)-0.5E0*(P(5)-P(3)))/H
S4=PPIN/P(4)
DE=(S2-S4)/(S1-S3)
IF(ABS(DE/E)-THRESH) 300,205,205
205 E=E+DE
IF(E)10,206,206
206 E=E-DE
DE=DE/2.0
GO TO 205
C IMPROVE TRIAL EIGENVALUE BY PERTURBATION THEORY IF NECESSARY
C CALCULATE THE NORMALIZED WAVE FUNCTIONS
300 POP=PMATCH/P(4)
DO 302 J=IMATCH,KKK
302 SNLO(J)=SNLO(J)*POP
SUM1=0.0
J=1
XIF=XIFC
VALUE=0.0
303 MM=8
XIF=XIF+XIF
SUM2=0.0
304 Y=VALUE
VALUE=SNLO(J+5)**2
SUM2=SUM2+19.0*(VALUE+Y)+75.0*(SNLO(J+4)**2+SNLO(J+1)**2)
1+50.0*(SNLO(J+2)**2+SNLO(J+3)**2)
J=J+5
MM=MM-1
IF(MM)77,305,304
305 SUM1=SUM1+XIF*SUM2
IF(KKK-J)77,307,303
307 C1=SQRT(SUM1)
IF(SNLO(3))308,77,310
308 C1=-C1
310 DO 311 I=1,KKK
311 SNLO(I)=SNLO(I)/C1
EN=E
RETURN
END
C ----------------------------------------------------------------------
C ----------------------------------------------------------------------
SUBROUTINE INTEGR(X,Y,R,N)
C SIMPSON'S RULE INTEGRATION SCHEME
DIMENSION X(521),Y(521),R(521)
H=R(2)
Y(1)=0.0
Y(2)=H*(5.0*X(1)+8.0*X(2)-X(3))/12.0
DO 20 J=1,N
DO 10 K=1,40
I=40*(J-1)+K
IF(I-40*N) 5,10,10
5 Y(I+2)=Y(I)+H*(X(I)+4.0*X(I+1)+X(I+2))/3.0
10 CONTINUE
H=H+H
IF(I-40*N) 15,20,15
15 Y(I+2)=Y(I+1)+H*(5.0*X(I+1)+8.0*X(I+2)-X(I+3))/12.0
20 CONTINUE
RETURN
END
C ----------------------------------------------------------------------
C ----------------------------------------------------------------------
SUBROUTINE CROSYM(M)
C SIMULTANEOUS EQUATION SOLVER
C WRITTEN BY I.C. HANSON, SCIENTIFIC COMPUTATION DEPARTMENT,
C LOCKHEED MISSILES AND SPACE COMPANY, SUNNYVALE, CALIFORNIA
C SOLVE M SIMULTANEOUS EQUATIONS BY THE METHOD OF CROUT
DIMENSION RSCORE(521)
COMMON XNUM(521),RSVALE(521),RSATOM(521),V(521),
1 SNLO(521),R(521),RU(521),RUEXCH(521),XI(521),XJ(521),
2 SNL(14,521),RUINL1(521),RUFNL1(521),RUINL2(521),RUFNL2(521)
3 ,RU2(521),RU3(521),NNLZ(24),WWNL(24),NKKK(24),EE(24),A(4,5)
N=M+1
I1=1
1 I3=I1
SUM=ABS(A(I1,I1))
DO3I=I1,M
IF(SUM-ABS(A(I,I1)))2,3,3
2 I3=I
SUM=ABS(A(I,I1))
3 CONTINUE
IF(I3-I1)4,6,4
4 DO 5J=1,N
SUM=-A(I1,J)
A(I1,J)=A(I3,J)
5 A(I3,J)=SUM
6 I3=I1+1
DO7I=I3,M
SUMP=ABS(A(I1,I1))
IF(SUMP.EQ.0.0) GO TO 65
A(I,I1)=A(I,I1)/A(I1,I1)
GO TO 7
65 A(I,I1)=0.0
7 CONTINUE
14 J2=I1-1
I3=I1+1
IF(J2)8,11,8
8 DO9J=I3,N
DO9I=1,J2
9 A(I1,J)=A(I1,J)-A(I1,I)*A(I,J)
IF(I1-M)11,15,11
11 J2=I1
I1=I1+1
DO12I=I1,M
DO12J=1,J2
12 A(I,I1)=A(I,I1)-A(I,J)*A(J,I1)
IF(I1-M)1,14,1
15 DO17I=1,M
J2=M-I
I3=J2+1
SUMQ=ABS(A(I3,I3))
IF(SUMQ.EQ.0.0) GO TO 155
A(I3,N)=A(I3,N)/A(I3,I3)
GO TO 158
155 A(I3,N)=0.0
158 IF(J2)16,18,16
16 DO17J=1,J2
17 A(J,N)=A(J,N)-A(I3,N)*A(J,I3)
18 RETURN
END
C ----------------------------------------------------------------------
C ----------------------------------------------------------------------
SUBROUTINE HYDROG
DIMENSION X(441),R(441),RSATOM(441)
MESH=441
NBLOCK=MESH/40
I=1
X(I)=0.0E0
R(1)=0.0E0
DELTAX=0.0025E0
DO 20 J=1,NBLOCK
DO 10 K=1,40
I=I+1
10 X(I)=X(I-1)+DELTAX
20 DELTAX=DELTAX+DELTAX
DO 30 I=2,MESH
30 R(I)=X(I)*0.88534138E0
DO 45 I=1,MESH
IF(R(I).GE.40.0E0) GO TO 40
RSATOM(I)=4.0E0*R(I)*R(I)*EXP(-2.0E0*R(I))
GO TO 45
40 RSATOM(I)=0.0E0
45 CONTINUE
WRITE(7,100)
100 FORMAT(4H 1,4H 1,4H 0,2(4H 1),4H 0)
NC=0
DO 50 MIN=1,441,5
MAX=MIN+4
MAX=MIN0(MAX,441)
NC=NC+1
WRITE(7,200) (RSATOM(I),I=MIN,MAX),NC
50 CONTINUE
200 FORMAT(1PE15.7,1P4E14.7,1X,"T",3H 1,I4)
RETURN
END
C ----------------------------------------------------------------------
C ----------------------------------------------------------------------
SUBROUTINE DUMP
C SUBROUTINE DUMP (DUMMY PROGRAM FOR HERMAN'S SCF)
STOP
END
__________________________________________________________________________