c -----non-relativistic LSD atomic structure
calculation---------
c writen by
Masahiro Yamamoto,
c
Institute of Atomic Enenrgy, Kyoto University
c
March 1994
c ============== INPUT ===============
c file:atomlsd.dat
c
z : atomic number
c
no: # of orbitals considered
c
spin up and down are different orbitals
c
npq(no): pricipal quantum number
c
elu(no): angular momentum quantum number
c
nspin(no): spin quantum number up=1 down=2
c
wei(no) : # of electrons
c
e(no) : orbital energy in au
c
rv(2,651) : spin-dependent potential r*v(r)_sigma
c
iskip: if iskip=1,eigenvalue and eigenfunction are
c
calculated for wei=0.0
c ============ OUTPUT ===================
c FILE='RSCFEN.DAT':eigenvalue
c FILE='RSCFPO.DAT':potential for up spin and
down spin
c FILE='ROSUM.DAT':electron density,spin up and
down density
c FILE='KINETIC.DAT':kinetic energy
c file='wavef.dat':wavefunctions
implicit real*8 (a-g,o-z)
dimension
npq(40),elu(40),nspin(40),wei(40),e(40)
dimension
rv(2,651),r(651),x(651),f(651),g(651),rhs(651)
dimension
xx(4),y(4),s(4),cof(4),a(5),b(4),p1(2,651),p2(2,651)
dimension
p3(2,651),fo(651),go(651),fw(40,651),rosum(2,651)
dimension ro(651),rvsig(2,651),ek(40)
c------------------- data input ------------------
open (1,file='atomlsd.dat')
read (1,*) z
read (1,*) no,iskip
write (6,*) 'z=',z,'no=',no
do i =1 , no
read (1,*)
npq(i),elu(i),nspin(i),wei(i),e(i)
write (6,*)
npq(i),elu(i),nspin(i),wei(i),e(i)
enddo
read (1,20) (rv(1,i), i=1, 651)
read (1,20) (rv(2,i), i=1, 651)
20 format (5(1pd15.7))
close (1)
write (6,20) (rv(1,i), i=1, 20)
write (6,20) (rv(2,i), i=1, 20)
c -----coordinate----log. mesh
--------------------------
dx=1.0d0/50.0d0
do 40 i=1, 651
x(i)=-9.0d0+dble(i-1)*dx
r(i)=dexp(x(i))
40 continue
c
---------------potential-----------------------
do i=1 , 651
do j=1,2
p1(j,i)=rv(j,i)
p2(j,i)=rv(j,i)
p3(j,i)=rv(j,i)
enddo
enddo
c ------scf loop
1000-------------------------------
do 1000 iscf=1, 1000
c -------------- nlspin loop 2000----------
do 2000 inls=1, no
if ( iskip .eq. 1) go to 52
if (wei(inls) .le. 1.0d-4)
go to 2000
52 iflag=0
c --------- classical turning point
---------------
2010 iflag=iflag+1
do 55 i=1, 650
f(i)=0.0
g(i)=0.0
fo(i)=0.0
go(i)=0.0
55 continue
do 50 i= 1, 651
if
((rv(nspin(inls),i)/r(i)-e(inls)) .ge. 0.0) go to 60
50 continue
60 im=i
if (mod(im,2) .eq. 0) im=im+1
c -------- outermost point
---------
do 70 i=im, 651
if
((rv(nspin(inls),i)/r(i)-e(inls))*r(i)*r(i) .ge. 75.0)
& go to 80
70 continue
80 it=i
if (mod(it,2) .eq. 0) it=it+1
if (it .ge. 651) it=651
write (6,*) 'classical t.
p.', im,'outermost point',it
c --------------------- radial sch. eq.
r.h.s. ---------
do 110 j=1, 651
rhs(j)=2.0d0*rv(nspin(inls),j)-2.0d0*e(inls)*r(j)
if (nint(elu(inls))
.eq. 0) go to 110
rhs(j)=rhs(j)+(elu(inls)*(elu(inls)+1.0d0))/r(j)
110 continue
c ------------- r:small region / power
expansion approxi. -----
c --- potential rv is expanded by
sum[i=0to3]b(i)*r**i near r=0
do 112 i=1, 4
y(i)=rv(nspin(inls),i)
xx(i)=r(i)
112 continue
n=4
do 11 i=1, n
s(i)=0.
cof(i)=0.
11 continue
s(n)=-xx(1)
do 13 i=2, n
do 12 j=n+1-i, n-1
s(j)=s(j)-xx(i)*s(j+1)
12 continue
s(n)=s(n)-xx(i)
13 continue
do 16 j=1,n
phi=n
do 14 k=n-1, 1, -1
phi=k*s(k+1)+xx(j)*phi
14 continue
fff=y(j)/phi
bb=1.
do 15 k=n,1,-1
cof(k)=cof(k)+bb*fff
bb=s(k)+xx(j)*bb
15 continue
16 continue
c write (6,*) (i,cof(i), i=1, n)
do 113 i=1, 4
b(i)=cof(i)
113 continue
c --near r=0 wave function f is expanded by
sum[i=0to4]a(i)*r**i
eluu=elu(inls)
ee=e(inls)
a(1)=1.0d0
a(2)=b(1)*a(1)/(eluu+1.0d0)
a(3)=(b(1)*a(2)+(b(2)-ee)*a(1))/(2.0d0*eluu+3.0d0)
a(4)=(b(1)*a(3)+(b(2)-ee)*a(2)+b(3)*a(1))/(3.0d0*eluu+6.0d0)
a(5)=(b(1)*a(4)+(b(2)-ee)*a(3)+b(3)*a(2)+b(4)*a(1))
1 /(4.0d0*eluu+10.0d0)
do 116 j=1, 4
f(j)=0.0
g(j)=0.0
do 114 i=1, 5
iram=i-1
f(j)=f(j)+a(i)*r(j)**(nint(eluu)+1+iram)
g(j)=g(j)+a(i)*(eluu+1.0+dble(iram))*r(j)**(nint(eluu)+iram)
114 continue
116 continue
C ***** OUTWARD ADAMS-MOULTON METHOD
R(4)->R(IM-) ************
I=4
FF1=DEXP(X(I))*G(I)
FF2=DEXP(X(I-1))*G(I-1)
FF3=DEXP(X(I-2))*G(I-2)
FF4=DEXP(X(I-3))*G(I-3)
GG1=rhs(i)*F(I)
GG2=rhs(i-1)*F(I-1)
GG3=rhs(i-2)*F(I-2)
GG4=rhs(i-3)*F(I-3)
C ----------- OUTWARD LOOP 1001 ----------
DO 1001 I=4, im-1
C ------- PREDICT -------
FPR=F(I)+dx/24.0d0*(55.0d0*FF1-59.0d0*FF2+37.0d0*FF3-9.0d0*FF4)
GPR=G(I)+dx/24.0d0*(55.0d0*GG1-59.0d0*GG2+37.0d0*GG3-9.0d0*GG4)
C ------- CORRECT -------
FF0=DEXP(X(I+1))*GPR
GG0=rhs(i+1)*FPR
FCR=F(I)+dx/24.0d0*(9.0d0*FF0+19.0d0*FF1-5.0d0*FF2+FF3)
GCR=G(I)+dx/24.0d0*(9.0d0*GG0+19.0d0*GG1-5.0d0*GG2+GG3)
DO 111 J=1, 2
FF0=DEXP(X(I+1))*GCR
GG0=rhs(I+1)*FCR
FCR=F(I)+dx/24.0d0*(9.0d0*FF0+19.0d0*FF1-5.0d0*FF2+FF3)
GCR=G(I)+dx/24.0d0*(9.0d0*GG0+19.0d0*GG1-5.0d0*GG2+GG3)
111 CONTINUE
F(I+1)=FCR
G(I+1)=GCR
C ------------ PARAMETER EXCHANGE
FOR THE NEXT LOOP -------
FF4=FF3
FF3=FF2
FF2=FF1
FF1=DEXP(X(I+1))*G(I+1)
GG4=GG3
GG3=GG2
GG2=GG1
GG1=rhs(i+1)*F(I+1)
1001 CONTINUE
c ----check of the node
------------------------------
iinod=0
do 150 i=1, im-1
if (f(i)*f(i+1) .le. 0.0)
iinod=iinod+1
150 continue
inode=npq(inls)-nint(elu(inls))-1
if (inode-iinod) 155,170,160
155 e(inls)=1.25*e(inls)
write (6,*) '# of node',inode,'# of node
in this cal',iinod
go to 2010
160 e(inls)=0.75*e(inls)
write (6,*) '# of node',inode,'# of node
in this cal',iinod
go to 2010
c -------- inward adams-moulton method r(it)->r(im+) ----------
c -- input 4 parameters -----
170 continue
if (rhs(it) .lt. 0.0) stop 9999
dq=sqrt(rhs(it))
frf=1.0
do 120 ii=1,4
i=it+1-ii
fo(i)=frf*dexp(-dq*r(i))
go(i)=-dq*fo(i)
120 continue
i=it-3
FF1=DEXP(X(I))*Go(I)
FF2=DEXP(X(I+1))*Go(I+1)
FF3=DEXP(X(I+2))*Go(I+2)
FF4=DEXP(X(I+3))*Go(I+3)
GG1=rhs(i)*Fo(I)
GG2=rhs(i+1)*Fo(I+1)
GG3=rhs(i+2)*Fo(I+2)
GG4=rhs(i+3)*Fo(I+3)
c ********** inward loop 200 ******
do 200 ii=3, it-im-1
i=it-ii
C ------- PREDICT -------
FPR=Fo(I)-dx/24.0*(55.0*FF1-59.0*FF2+37.0*FF3-9.0*FF4)
GPR=Go(I)-dx/24.0*(55.0*GG1-59.0*GG2+37.0*GG3-9.0*GG4)
C ------- CORRECT -------
FF0=DEXP(X(I-1))*GPR
GG0=rhs(i-1)*FPR
FCR=Fo(I)-dx/24.0*(9.0*FF0+19.0*FF1-5.0*FF2+FF3)
GCR=Go(I)-dx/24.0*(9.0*GG0+19.0*GG1-5.0*GG2+GG3)
DO 1111 J=1, 2
FF0=DEXP(X(I-1))*GCR
GG0=rhs(I-1)*FCR
FCR=Fo(I)-dx/24.0*(9.0*FF0+19.0*FF1-5.0*FF2+FF3)
GCR=Go(I)-dx/24.0*(9.0*GG0+19.0*GG1-5.0*GG2+GG3)
1111 CONTINUE
Fo(I-1)=FCR
Go(I-1)=GCR
C ------------ PARAMETER EXCHANGE
FOR THE NEXT LOOP -------
FF4=FF3
FF3=FF2
FF2=FF1
FF1=DEXP(X(I-1))*Go(I-1)
GG4=GG3
GG3=GG2
GG2=GG1
GG1=rhs(i-1)*Fo(I-1)
200 CONTINUE
c ------------- integral of f r=0 to classical
t. p.--------------
sum=0.0d0+r(1)*f(1)*f(1)/2.0d0
do 220 i=2, im-1, 2
sum=sum+dx/3.0d0*(dexp(x(i-1))*f(i-1)**2+4.0d0*dexp(x(i))
1
*f(i)**2+dexp(x(i+1))*f(i+1)**2)
220 continue
c --------- integral of fo r=im to infinity
------
sumout=0.0+fo(it)*fo(it)/2.0/dq
do 225 i=im+1, it-1, 2
sumout=sumout+dx/3.0d0*(dexp(x(i-1))*fo(i-1)**2
&
+4.0d0*dexp(x(i))*fo(i)**2+dexp(x(i+1))*fo(i+1)**2)
225 continue
c ------- energy correction estimate see
hatree's book ---
dele=(-g(im)/f(im)+go(im)/fo(im))/(sum/f(im)/f(im)+sumout
&
/fo(im)/fo(im))
if (abs(dele/e(inls)) .gt. 0.4)
dele=0.8*abs(e(inls)/dele)*dele
e(inls)=e(inls)-dele/2.0
write (6,*) 'iscf/inls,iflag',
iscf,inls,iflag,e(inls),
$
dele/e(inls)
if (abs(dele/e(inls)) .ge. 1.0d-5) go to
2010
c --- if energy conv. is not satisfied go to
2010 --------
c ---- matching the f and fo at
classical t. p. ----
301 frfold=frf
if (abs((f(im)-fo(im))/f(im)) .le.
1.0d-6) go to 300
frf = frfold*f(im)/fo(im)
c --------- recal fo from calssical t. p. to
infinity
do 121 ii=1,4
i=it+1-ii
fo(i)=frf*dexp(-dq*r(i))
go(i)=-dq*fo(i)
121 continue
i=it-3
FF1=DEXP(X(I))*Go(I)
FF2=DEXP(X(I+1))*Go(I+1)
FF3=DEXP(X(I+2))*Go(I+2)
FF4=DEXP(X(I+3))*Go(I+3)
GG1=rhs(i)*Fo(I)
GG2=rhs(i+1)*Fo(I+1)
GG3=rhs(i+2)*Fo(I+2)
GG4=rhs(i+3)*Fo(I+3)
c ********** inward loop 201 ******
do 201 ii=3, it-im-1
i=it-ii
C ------- PREDICT -------
FPR=Fo(I)-dx/24.0*(55.0*FF1-59.0*FF2+37.0*FF3-9.0*FF4)
GPR=Go(I)-dx/24.0*(55.0*GG1-59.0*GG2+37.0*GG3-9.0*GG4)
C ------- CORRECT -------
FF0=DEXP(X(I-1))*GPR
GG0=rhs(i-1)*FPR
FCR=Fo(I)-dx/24.0*(9.0*FF0+19.0*FF1-5.0*FF2+FF3)
GCR=Go(I)-dx/24.0*(9.0*GG0+19.0*GG1-5.0*GG2+GG3)
DO 1112 J=1, 2
FF0=DEXP(X(I-1))*GCR
GG0=rhs(I-1)*FCR
FCR=Fo(I)-dx/24.0*(9.0*FF0+19.0*FF1-5.0*FF2+FF3)
GCR=Go(I)-dx/24.0*(9.0*GG0+19.0*GG1-5.0*GG2+GG3)
1112 CONTINUE
Fo(I-1)=FCR
Go(I-1)=GCR
C ------------ PARAMETER EXCHANGE
FOR THE NEXT LOOP -------
FF4=FF3
FF3=FF2
FF2=FF1
FF1=DEXP(X(I-1))*Go(I-1)
GG4=GG3
GG3=GG2
GG2=GG1
GG1=rhs(i-1)*Fo(I-1)
201 CONTINUE
go to 301
c ------------- integral of f*f r=0 to
infinity--------------
300 continue
sum=0.0d0+r(1)*f(1)*f(1)/2.0d0
do 221 i=2, im-1, 2
sum=sum+dx/3.0d0*(dexp(x(i-1))*f(i-1)**2+4.0d0*dexp(x(i))
1
*f(i)**2+dexp(x(i+1))*f(i+1)**2)
221 continue
sumout=0.0+fo(it)*fo(it)/2.0/dq
do 226 i=im+1, it-1, 2
sumout=sumout+dx/3.0d0*(dexp(x(i-1))*fo(i-1)**2
&
+4.0d0*dexp(x(i))*fo(i)**2+dexp(x(i+1))*fo(i+1)**2)
226 continue
sum=sum+sumout
c ------------------- normalization
-----------
do i=1, im
fw(inls,i)=f(i)/sqrt(sum)
enddo
do i=im+1,it
fw(inls,i)=fo(i)/sqrt(sum)
enddo
if (it .eq. 651) go to 227
do i=it+1, 651
fw(inls,i)=0.0
enddo
227 continue
c --------------- nomalization again
-----------------
229 sum=0.0+r(1)*fw(inls,1)*fw(inls,1)/2.0d0
do i=2, 650, 2
sum=sum+dx/3.0d0*(dexp(x(i-1))*fw(inls,i-1)**2+
&
4.0d0*dexp(x(i))*fw(inls,i)**2+dexp(x(i+1))*fw(inls,i+1)**2)
enddo
if (abs(sum-1.0d0) .le. 1.0d-10) go to
228
do i=1, 651
if (abs(fw(inls,i)) .gt. 1.0d-11)
&
fw(inls,i)=fw(inls,i)/sqrt(sum)
enddo
go to
229
228 continue
c ------------ summation of electron density
------
c ------------- rosum = 4*pai*ro(r)
------------------
do i=1, it
rosum(nspin(inls),i)=rosum(nspin(inls),i)
&
+wei(inls)*fw(inls,i)*fw(inls,i)/r(i)/r(i)
enddo
2000 continue
c ------ end of inls 2000 loop -------
c ------ rosum check
integ[r*r*rosum]=integ[4*pai*r*r*ro(r)]-----
SEKI1=0.0+r(1)*r(1)*r(1)*rosum(1,1)/2.0
seki2=0.0+r(1)*r(1)*r(1)*rosum(2,1)/2.0
scrnc=seki1+seki2
seki3=0.0+scrnc
do i=1, 651
ro(i)=rosum(1,i)+rosum(2,i)
enddo
DO 505 J=2, 650, 2
SEKI1=SEKI1+dx/3.0*(DEXP(3.0*X(J-1))*ROSUM(1,J-1)+4.0*
&
DEXP(3.0*X(J))*ROSUM(1,J)+DEXP(3.0*X(J+1))*ROSUM(1,J+1))
SEKI2=SEKI2+dx/3.0*(DEXP(3.0*X(J-1))*ROSUM(2,J-1)+4.0*
&
DEXP(3.0*X(J))*ROSUM(2,J)+DEXP(3.0*X(J+1))*ROSUM(2,J+1))
SEKI3=SEKI3+dx/3.0*(DEXP(3.0*X(J-1))*RO(J-1)+4.0*
&
DEXP(3.0*X(J))*RO(J)+DEXP(3.0*X(J+1))*RO(J+1))
505 CONTINUE
WRITE (6,*) 'atomic number=', z,
'rosum spin up=',SEKI1
write (6,*) 'rosum spin
down=',seki2, 'rosum=',seki3
C ################# MODIFYING THE
POTENTIAL ################
C ******** HARTREE POTENTIAL
****************
PAI4=4.0*3.141592654
DO 500 I=1, 651
SEKI1=0.0D0
SEKI2=0.0D0
IF (MOD(I,2)) 520,520,600
C ------------ FOR EVEN I ----------
520 IF (I .EQ. 2) GO TO 530
DO 540 J=2, I-2, 2
SEKI1=SEKI1+dx/3.0*(DEXP(3.0*X(J-1))*RO(J-1)+4.0*DEXP(3.0*
1
X(J))*RO(J)+DEXP(3.0*X(J+1))*RO(J+1))
540 CONTINUE
530 SEKI1=SEKI1+dx/2.0*(DEXP(3.0*X(I-1))*RO(I-1)+
1
DEXP(3.0*X(I))*RO(I))
SEKI1=SEKI1+SCRNC
IF (I .EQ. 650) GO TO 560
DO 550 J=I+1, 649, 2
SEKI2=SEKI2+dx/3.0*(DEXP(2.0*X(J-1))*RO(J-1)+4.0*DEXP(2.0*
1
X(J))*RO(J)+DEXP(2.0*X(J+1))*RO(J+1))
550 CONTINUE
560 SEKI2=SEKI2+dx/2.0*(DEXP(2.0*X(650))*RO(650)+
1
DEXP(2.0*X(651))*RO(651))
SEKI2=R(I)*SEKI2
GO TO 680
C ------------- FOR ODD I
-----------
600 IF (I. EQ. 1) GO TO 620
DO 610 J=2, I-1, 2
SEKI1=SEKI1+dx/3.0*(DEXP(3.0*X(J-1))*RO(J-1)+4.0*DEXP(3.0*
1
X(J))*RO(J)+DEXP(3.0*X(J+1))*RO(J+1))
610 CONTINUE
620 SEKI1=SEKI1+SCRNC
IF (I .EQ. 651) GO TO 640
DO 650 J=I+1, 650, 2
SEKI2=SEKI2+dx/3.0*(DEXP(2.0*X(J-1))*RO(J-1)+4.0*DEXP(2.0*
1
X(J))*RO(J)+DEXP(2.0*X(J+1))*RO(J+1))
650 CONTINUE
640 SEKI2=R(I)*SEKI2
C ******** spin-dependent
EXCHANGE-CORRELATION POTENTIAL ******
c pewdew zunger
parameterization PRB 1981 23 p5048
c von barth & hedin
form: J.Phys.c 1972 5 p1629
c vosko wilk nusair Can
J Phys 1980 58 p1200
pai=acos(-1.0)
680 IF (RO(I)) 690, 690, 700
690 PXup=0.0
pxdown=0.0
PCup=0.0
pcdown=0.0
GO TO 751
c ------- rosum=(4*pai)*ro(r)
-----------------
700 RS=(3.0/RO(I))**(1.0/3.0)
c ------------ spin polarization
spol--------------
if
(abs(rosum(1,i)-rosum(2,i)) .le. 1.0d-12) then
spol=0.0
fspol=0.0
dfspol=0.0
go to 704
endif
spol=(rosum(1,i)-rosum(2,i))/ro(i)
fspol=((1.0+spol)**(4.0/3.0)+(1.0-spol)**(4.0/3.0)-2.0)
& /(2.0**(4.0/3.0)-2.0)
dfspol=(4.0/3.0)*((1.0+spol)**(1.0/3.0)-(1.0-spol)**(1.0/3.0))
&
/(2.0**(4.0/3.0)-2.0)
c ------------ exchange
----------------------------------
704 exunpo=-(3.0/4.0/pai)*(9.0*pai/4.0)**(1.0/3.0)/rs
expo=exunpo*(2.0)**(1.0/3.0)
pxunpo=(4.0/3.0)*exunpo
pxpo=(4.0/3.0)*expo
c ----- unpo:unpolalized
,po:polarized,ex:exchange enrgy
c ----- px:exchange potential
ex=exunpo+(expo-exunpo)*fspol
pxup=pxunpo+fspol*(pxpo-pxunpo)+(expo-exunpo)*
&
(1.0-spol)*dfspol
pxdown=pxunpo+fspol*(pxpo-pxunpo)+(expo-exunpo)*
&
(-1.0-spol)*dfspol
c ------------- correlation
--------------------
if (rs-1.0) 701,702,702
701 pcunpo=0.0311*log(rs)+(-0.048-0.0311/3.0)+
&
(2.0/3.0)*0.0020*rs*log(rs)+(2.0*(-0.0116)-0.0020)/3.0*rs
pcpo=0.01555*log(rs)+(-0.0269-0.01555/3.0)+
&
(2.0/3.0)*0.0007*rs*log(rs)+(2.0*(-0.0048)-0.0007)/3.0*rs
ecunpo=0.0311*log(rs)-0.048+0.0020*rs*log(rs)-0.0116*rs
ecpo=0.01555*log(rs)-0.0269+0.0007*rs*log(rs)-0.0048*rs
go to 750
702 ecunpo=-0.1423/(1.0+1.0529*sqrt(rs)+0.3334*rs)
ecpo=-0.0843/(1.0+1.3981*sqrt(rs)+0.2611*rs)
pcunpo=ecunpo*(1.0+(7.0/6.0)*1.0529*sqrt(rs)+
&
(4.0/3.0)*0.3334*rs)/(1.0+1.0529*sqrt(rs)+0.3334*rs)
pcpo=ecpo*(1.0+(7.0/6.0)*1.3981*sqrt(rs)+
&
(4.0/3.0)*0.2611*rs)/(1.0+1.3981*sqrt(rs)+0.2611*rs)
750 pcup=pcunpo+fspol*(pcpo-pcunpo)+(ecpo-ecunpo)
&
*(1.0-spol)*dfspol
pcdown=pcunpo+fspol*(pcpo-pcunpo)+(ecpo-ecunpo)
&
*(-1.0-spol)*dfspol
c --------- spin-dependent
potential r*v(r)sigma -----
751 rvsig(1,i)=-Z+SEKI1+SEKI2+R(I)*(PCup+PXup)
rvsig(2,i)=-Z+SEKI1+SEKI2+R(I)*(PCdown+PXdown)
P1(1,I)=rv(1,I)
p1(2,i)=rv(2,i)
C ----------- CONSTRUCT THE NEW
POTENTIAL ---------
do jj=1, 2
DENOM=P3(jj,I)+rvsig(jj,i)-P1(jj,I)-P2(jj,I)
IF (DENOM .LE.
1.0D-10) GO TO 800
IF
(P1(jj,I)-rvsig(jj,i) .LE. 1.0D-10) GO TO 800
AJUD=((P3(jj,I)*rvsig(jj,i)-P1(jj,I)*P2(jj,I))/DENOM
&
-rvsig(jj,i))/(P1(jj,I)-rvsig(jj,i))
IF (AJUD .GE. 0.5) GO
TO 800
rv(jj,I)=(P3(jj,I)*rvsig(jj,i)-P1(jj,I)*P2(jj,I))/DENOM
IF (AJUD .LE. 0.0)
rv(jj,I)=rvsig(jj,i)
GO TO 820
800 rv(jj,I)=(P1(jj,I)+rvsig(jj,i))/2.0
C -------------- POTENTIAL
PARAMETER-----------------
820 IF (DABS(rv(jj,I)-P1(jj,I)) .GE. PJUD)
&
PJUD=DABS(rv(jj,I)-P1(jj,I))
enddo
500 CONTINUE
C
##########################################################
DO 420 I=1, 651
do j=1, 2
P3(j,I)=P1(j,I)
P2(j,I)=rvsig(j,i)
rvsig(j,i)=0.0
enddo
420 CONTINUE
C -------- POTENTIAL CONVERGENCE
---------------------
IF (PJUD. LE. 5.0D-5) GO TO 3000
C -------- ESTIMATE OF ENERGY
CHANGE BY PURTURBATION -----
DO 900 I=1, NO
IF (wei(I) .EQ. 0.0) GO TO
900
FF0=0.0D0
do j=1, 651
f(j)=fw(i,j)
enddo
DO 920 J=2, 650, 2
FF1=(F(J-1)**2.0)*(rv(nspin(i),J-1)-P1(nspin(i),J-1))
FF2=(F(J)**2.0)*(rv(nspin(i),J)-P1(nspin(i),J))
FF3=(F(J+1)**2.0)*(rv(nspin(i),J+1)-P1(nspin(i),J+1))
FF0=FF0+dx/3.0*(FF1+4.0*FF2+FF3)
920 CONTINUE
E(I)=E(I)+FF0/2.0
900 CONTINUE
C --------- CLEAR THE SUMMING PARAMETERS
---------------
DO 410 I=1, 651
ROSUM(1,I)=0.0
rosum(2,i)=0.0
410 CONTINUE
SCRNC=0.0D0
PJUD=0.0D0
c OPEN(1, FILE='RSCFEN.DAT')
c WRITE (1,3010) (E(I), I=1, NO)
c CLOSE (1)
c OPEN(1, FILE='RSCFPO.DAT')
c WRITE (1,3010) (P(I), I=1, 651)
c CLOSE (1)
1000 CONTINUE
C ************* SCF LOOP END
**********************************
c --------- kinetic energy
---------------------------
3000 continue
do 3002 inls=1, 40
ek(inls)=0.0
if (wei(inls) .le. 1.0d-4) go to
3002
do 3004 i=1, 651
rhs(i)=(e(inls)-rv(nspin(inls),i)/r(i))*fw(inls,i)*fw(inls,i)
3004 continue
ek(inls)=0.0+r(1)*rhs(1)/2.0
do i=2, 650, 2
ek(inls)=ek(inls)+dx/3.0d0*(dexp(x(i-1))*rhs(i-1)+
&
4.0d0*dexp(x(i))*rhs(i)+dexp(x(i+1))*rhs(i+1))
enddo
3002 continue
C ************** OUTPUT THE RESULTS CALCULATED
****************
DO 3030 I=1, 5
WRITE (6,3020)
3030 CONTINUE
3020 FORMAT (1H ,'!!!!!!!!!!!!! CALCULATION DONE
!!!!!!!!!!!!!!!')
OPEN(1, FILE='RSCFEN.DAT')
WRITE (1,3010) (E(I), I=1, 40)
WRITE (6,3010) (E(I), I=1, 40)
3010 FORMAT (5(1PD15.7))
CLOSE (1)
OPEN(1, FILE='RSCFPO.DAT')
WRITE (1,3010) (rv(1,I), I=1, 651)
write (1,3010) (rv(2,i), i=1, 651)
CLOSE (1)
OPEN(1, FILE='ROSUM.DAT')
write (1,3010) (ro(i), i =1, 651)
WRITE (1,3010) (ROSUM(1,I), I=1, 651)
write (1,3010) (rosum(2,i), i=1,651)
CLOSE (1)
OPEN(1, FILE='KINETIC.DAT')
WRITE (1,3010) (ek(i), i=1, 40)
WRITE (6,* ) PJUD
CLOSE (1)
open (1, file='wavef.dat')
do inls=1, no
write (1,3010) (fw(inls,i),
i=1, 651)
enddo
close (1)
stop
end