c Hedin-Lundqvist version 2000 April, EHCC,
Kyoto University
c -----non-relativistic LDA atomic structure
calculation---------
c writen by
Masahiro Yamamoto,
c
Institute of Atomic Enenrgy, Kyoto University
c
March 1994
c ======= INPUT
==================================
c file:atomlda.dat
c
z : atomic number
c
no: # of orbitals considered
c
npq(no): pricipal quantum number
c
elu(no): angular momentum quantum number
c
spin(no):
dummy
c
wei(no) : # of electrons
c
e(no) : orbital energy in
a.u.
c
rv(651) : all electron potential r*v(r)
c
note:r(i)=dexp(x(i)),x(i)=-9.0d0+dble(i-1)*dx
c
dx=1.0d0/50.0d0
c iskip :if iskip
= 1, wave function and eigen value are
c
calculated for wei=0.0
c ============ output =================
c FILE='RSCFEN.DAT' :eigenvalue
c FILE='RSCFPO.DAT' :potential r*v(r)
c FILE='ROSUM.DAT' :electron
density
c FILE='KINETIC.DAT':kinetic energy
c file='wavef.dat' :wavefunction
c
=========================================
implicit real*8 (a-g,o-z)
dimension
npq(20),elu(20),spin(20),wei(20),e(20)
dimension
rv(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(651),p2(651)
dimension
p3(651),fo(651),go(651),fw(20,651),rosum(651),ek(20)
c------------------- data input ------------------
open (1,file='atomlda.dat')
read (1,*) z
read (1,*) no, iskip
write (6,*)
'z',z,'no',no
do i =1 , no
read (1,*)
npq(i),elu(i),spin(i),wei(i),e(i)
write (6,*) npq(i),elu(i),
spin(i),wei(i),e(i)
enddo
read (1,20) (rv(i), i=1, 651)
20 format (5(1pd15.7))
close (1)
write (6,20) (rv(i), i=1, 30)
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
p1(i)=rv(i)
p2(i)=rv(i)
p3(i)=rv(i)
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(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(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(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(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)
go to 2010
160 e(inls)=0.75*e(inls)
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(i)=rosum(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)/2.0
scrnc=seki1
DO 505 J=2, 650, 2
SEKI1=SEKI1+dx/3.0*(DEXP(3.0*X(J-1))*ROSUM(J-1)+4.0*
&
DEXP(3.0*X(J))*ROSUM(J)+DEXP(3.0*X(J+1))*ROSUM(J+1))
505 CONTINUE
WRITE (6,*) 'atomic number=', z,
'rosum=',SEKI1
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))*ROSUM(J-1)+4.0*DEXP(3.0*
1
X(J))*ROSUM(J)+DEXP(3.0*X(J+1))*ROSUM(J+1))
540 CONTINUE
530 SEKI1=SEKI1+dx/2.0*(DEXP(3.0*X(I-1))*ROSUM(I-1)+
1
DEXP(3.0*X(I))*ROSUM(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))*ROSUM(J-1)+4.0*DEXP(2.0*
1
X(J))*ROSUM(J)+DEXP(2.0*X(J+1))*ROSUM(J+1))
550 CONTINUE
560 SEKI2=SEKI2+dx/2.0*(DEXP(2.0*X(650))*ROSUM(650)+
1
DEXP(2.0*X(651))*ROSUM(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))*ROSUM(J-1)+4.0*DEXP(3.0*
1
X(J))*ROSUM(J)+DEXP(3.0*X(J+1))*ROSUM(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))*ROSUM(J-1)+4.0*DEXP(2.0*
1
X(J))*ROSUM(J)+DEXP(2.0*X(J+1))*ROSUM(J+1))
650 CONTINUE
640 SEKI2=R(I)*SEKI2
C ************ EXCHANGE-CORRELATION
POTENTIAL ***********
c pewdew zunger prb
1981 23 p5048
pai=acos(-1.0)
680 IF (ROSUM(I)) 690, 690, 700
690 PX=0.0D0
PC=0.0D0
GO TO 750
c ------- rosum=(4*pai)*ro(r)
-----------------
700 RS=(3.0/ROSUM(I))**(1.0/3.0)
ex=-(3.0/4.0/pai)*(9.0*pai/4.0)**(1.0/3.0)/rs
px=(4.0/3.0)*ex
c --- Hedin-Lundqvist ----
ahl=21.0d0
chl=0.0225d0
xhl=rs/ahl
aln=log(1.0d0+1.0d0/xhl)
pc=-chl*aln
c ----- CA --------
c if (rs-1.0) 701,702,702
c 701 pc=0.0311*log(rs)+(-0.048-0.0311/3.0)+
c &
(2.0/3.0)*0.0020*rs*log(rs)+(2.0*(-0.0116)-0.0020)/3.0*rs
c go to 750
c 702 pc=-0.1423*(1.0+(7.0/6.0)*1.0529*sqrt(rs)+
c &
(4.0/3.0)*0.3334*rs)/(1.0+1.0529*sqrt(rs)+0.3334*rs)**2.0
c -------------------------
750 F(I)=-Z+SEKI1+SEKI2+R(I)*(PC+PX)
P1(I)=rv(I)
C ----------- CONSTRUCT THE NEW
POTENTIAL ---------
DENOM=P3(I)+F(I)-P1(I)-P2(I)
IF (DENOM .LE. 1.0D-10) GO
TO 800
IF (P1(I)-F(I) .LE. 1.0D-10)
GO TO 800
AJUD=((P3(I)*F(I)-P1(I)*P2(I))/DENOM-F(I))/(P1(I)-F(I))
IF (AJUD .GE. 0.5) GO TO 800
rv(I)=(P3(I)*F(I)-P1(I)*P2(I))/DENOM
IF (AJUD .LE. 0.0) rv(I)=F(I)
GO TO 820
800 rv(I)=(P1(I)+F(I))/2.0
C -------------- POTENTIAL
PARAMETER-----------------
820 IF (DABS(rv(I)-P1(I)) .GE. PJUD)
PJUD=DABS(rv(I)-P1(I))
500 CONTINUE
C
##########################################################
DO 420 I=1, 651
P3(I)=P1(I)
P2(I)=F(I)
F(I)=0.0D0
420 CONTINUE
C -------- POTENTIAL CONVERGENCE
---------------------
write (6,*) 'pjud=',pjud
IF (PJUD. LE. 1.0D-4) GO TO 3000
if (iscf .gt. 990) 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(J-1)-P1(J-1))
FF2=(F(J)**2.0)*(rv(J)-P1(J))
FF3=(F(J+1)**2.0)*(rv(J+1)-P1(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(I)=0.0D0
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, 20
if (wei(inls) .le. 1.0d-4) go to
3002
do 3004 i=1, 651
rhs(i)=(e(inls)-rv(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, 20)
WRITE (6,3010) (E(I), I=1, 20)
3010 FORMAT (5(1PD15.7))
CLOSE (1)
OPEN(1, FILE='RSCFPO.DAT')
WRITE (1,3010) (rv(I), I=1, 651)
CLOSE (1)
OPEN(1, FILE='ROSUM.DAT')
WRITE (1,3010) (ROSUM(I), I=1, 651)
CLOSE (1)
OPEN(1, FILE='KINETIC.DAT')
WRITE (1,3010) (ek(i), i=1, 20)
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