C ****** TOTAL ENERGY CALCULATION FOR atomlsd.f
***********
c
non-relativistic Local Spin Density Approximation
c writen by
Masahiro Yamamoto,
c
Institute of Atomic Enenrgy, Kyoto University
c
March 1994
C NO:NUMBER OF OCCUPIED ORBITALS
C Wei(I):NUMBER OF ELECTRONS AT I LEVEL
C E(I):ENGENVALUES OF THE I ORBITAL
C ROSUM(I):4*PAI*DENTISY OF ELCTRONS AT THE I
POINT
C R(I):RADIUS OF I POINTS
C X(I)=LOG(R(I))
implicit real*8 (a-h,o-z)
dimension
F(651),G(2,651),E(40),ROSUM(2,651),X(651),R(651),
1
Wei(40),ro(651),npq(40),elu(40),nspin(40),ek(40),h(651)
c --------------- input --------------------
open (1,file='atomlsd.dat')
read (1,*) z
read (1,*) no
write (6,*)
'z',z,'no',no
do i =1 , no
read (1,*)
npq(i),elu(i),nspin(i),wei(i),ons
write (6,*)
npq(i),elu(i),nspin(i),wei(i)
enddo
close (1)
10 FORMAT (5(1PD15.7))
15 FORMAT (5(1PD10.3))
20 FORMAT (1PD15.7)
OPEN (1,FILE='RSCFEN.DAT')
READ (1,10) (E(I), I=1, 40)
CLOSE (1)
OPEN (1,FILE='ROSUM.DAT')
READ (1,10) (RO(I), I=1, 651)
read (1,10) (rosum(1,i), i=1, 651)
read (1,10) (rosum(2,i), i=1, 651)
CLOSE (1)
open (1,file='KINETIC.DAT')
read (1,10) (ek(i), i=1, 40)
close (1)
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 ************ SUMMATION OF ENERGY EIGENVALUE
****************
c --------- and orbital-dependent kinetic
energy ------
esum=0.0
eksum=0.0
DO 30 I=1, 20
ESUM=ESUM+Wei(I)*E(I)
eksum=eksum+wei(i)*ek(i)
30 CONTINUE
C ************ CHECK OF ROSUM + SCRNC
**********************
sum1=0.0+(r(1)**3)*rosum(1,1)/2.0
sum2=0.0+(r(1)**3)*rosum(2,1)/2.0
sum3=0.0+(r(1)**3)*ro(1)/2.0
scrnc=sum3
DO 700 I=2, 650, 2
SUM1=SUM1+dx/3.0*(DEXP(3.0*X(I-1))*ROSUM(1,I-1)+4.0*
1
DEXP(3.0*X(I))*ROSUM(1,I)+DEXP(3.0*X(I+1))*ROSUM(1,I+1))
SUM2=SUM2+dx/3.0*(DEXP(3.0*X(I-1))*ROSUM(2,I-1)+4.0*
1
DEXP(3.0*X(I))*ROSUM(2,I)+DEXP(3.0*X(I+1))*ROSUM(2,I+1))
SUM3=SUM3+dx/3.0*(DEXP(3.0*X(I-1))*RO(I-1)+4.0*
1
DEXP(3.0*X(I))*RO(I)+DEXP(3.0*X(I+1))*RO(I+1))
700 CONTINUE
WRITE (6,710) SUM1,sum2,sum3
710 FORMAT (1H ,'up, down-spin, total CHARGE=',3(1PD15.7))
C ************ HATREE POTENTIAL CORRECTION
*********************
DO 100 I=1, 651
SUM1=0.0D0
SUM2=0.0D0
IF (MOD(I,2)) 110,110,120
C ------ FOR EVEN I -------
110 IF (I .EQ. 2) GO TO 130
DO 140 J=2, I-2, 2
SUM1=SUM1+dx/3.0*(DEXP(3.0*X(J-1))*RO(J-1)+4.0*
1
DEXP(3.0*X(J))*RO(J)+DEXP(3.0*X(J+1))*RO(J+1))
140 CONTINUE
130
SUM1=SUM1+dx/2.0*(DEXP(3.0*X(I-1))*RO(I-1)+DEXP(3.0*X(I))*
1 RO(I))
SUM1=(SUM1+SCRNC)/R(I)
IF (I .EQ. 650) GO TO 150
DO 160 J=I+1, 649, 2
SUM2=SUM2+dx/3.0*(DEXP(2.0*X(J-1))*RO(J-1)+4.0*
1
DEXP(2.0*X(J))*RO(J)+DEXP(2.0*X(J+1))*RO(J+1))
160 CONTINUE
150
SUM2=SUM2+dx/2.0*(DEXP(2.0*X(650))*RO(650)+DEXP(2.0*X(651))
1 *RO(651))
GO TO 105
C ------ FOR ODD I -------
120 IF (I .EQ. 1) GO TO 170
DO 180 J=2, I-1, 2
SUM1=SUM1+dx/3.0*(DEXP(3.0*X(J-1))*RO(J-1)+4.0*
1
DEXP(3.0*X(J))*RO(J)+DEXP(3.0*X(J+1))*RO(J+1))
180 CONTINUE
170 SUM1=(SUM1+SCRNC)/R(I)
IF (I .EQ. 651) GO TO 195
DO 190 J=I+1, 650, 2
SUM2=SUM2+dx/3.0*(DEXP(2.0*X(J-1))*RO(J-1)+4.0*
1
DEXP(2.0*X(J))*RO(J)+DEXP(2.0*X(J+1))*RO(J+1))
190 CONTINUE
GO TO 105
195 SUM2=0.0D0
105 F(I)=-0.5D0*(SUM1+SUM2)
100 CONTINUE
C ************ EXCHANGE-CORRELATION
CORRECTION ****************
C ******** spin-dependent
EXCHANGE-CORRELATION POTENTIAL ******
c Pewdew Zunger PRB
1981 23 p5048
c Barth Hedin J.Phys.C
1972 5 p1629
c Vosko Wilk Nusair Can
J Phys 1980 58 p1200
pai=acos(-1.0)
do 400 i=1, 651
680 IF (RO(I)) 690, 690, 707
690 PXup=0.0
pxdown=0.0
PCup=0.0
pcdown=0.0
ex=0.0
ec=0.0
GO TO 751
c ------- rosum=(4*pai)*ro(r)
-----------------
707 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
ec=ecunpo+(ecpo-ecunpo)*fspol
c --------- spin-dependent
potential v(r)sigma -----
751 g(1,i)=ex+ec-(PCup+PXup)
g(2,i)=ex+ec-(PCdown+PXdown)
h(i)=ex+ec
c write (6,1000)
i,ro(i),rosum(1,i),rosum(2,i),spol
c write (6,1000)
i,ex,ec,pxup,pxdown,pcup,pcdown
1000 format(i5,6d15.7)
400 CONTINUE
C *********** INTEGRAL OVER THE WHOLE
SPHERE ****************
eint=dx/2.0*dexp(3.0*x(1))*f(1)*ro(1)
eint1=dx/2.0*dexp(3.0*x(1))*g(1,1)*rosum(1,1)
eint2=dx/2.0*dexp(3.0*x(1))*g(2,1)*rosum(2,1)
eint3=dx/2.0*dexp(3.0*x(1))*(-f(1)+h(1)-z/r(1))*ro(1)
DO 500 I=2, 650, 2
EINT=EINT+dx/3.0*(DEXP(3.0*X(I-1))*F(I-1)*RO(I-1)+
1
4.0*DEXP(3.0*X(I))*F(I)*RO(I)+DEXP(3.0*X(I+1))
2
*F(I+1)*RO(I+1))
EINT1=EINT1+dx/3.0*(DEXP(3.0*X(I-1))*g(1,I-1)*ROsum(1,I-1)+
1
4.0*DEXP(3.0*X(I))*g(1,I)*ROsum(1,I)+DEXP(3.0*X(I+1))
2
*g(1,I+1)*ROsum(1,I+1))
EINT2=EINT2+dx/3.0*(DEXP(3.0*X(I-1))*g(2,I-1)*ROsum(2,I-1)+
1
4.0*DEXP(3.0*X(I))*g(2,I)*ROsum(2,I)+DEXP(3.0*X(I+1))
2
*g(2,I+1)*ROsum(2,I+1))
EINT3=EINT3+dx/3.0*(DEXP(3.0*X(I-1))*(-f(i-1)+h(i-1)-z/r(i-1))
&
*RO(I-1)+4.0*DEXP(3.0*X(I))*(-f(i)+h(i)-z/r(i))*RO(I)
&
+DEXP(3.0*X(I+1))*(-f(i+1)+h(i+1)-z/r(i+1))*RO(I+1))
500 CONTINUE
C **************** TOTAL ENERGY
********************
ESUM=ESUM+EINT+eint1+eint2
eksum=eksum+eint3
OPEN (1,FILE='TOTE.DAT')
WRITE (1,610) ESUM,eksum
CLOSE (1)
WRITE (6,600) ESUM
600 FORMAT (1H ,'TOTAL ENERGY(AU) from eigenvalue=',1PD17.9)
WRITE (6,620) eksum
620 FORMAT (1H ,'TOTAL ENERGY(AU) from kinetic E.=',1PD17.9)
610 FORMAT (2PD17.9)
write (6,'(f15.7)')
(esum+eksum)/2.0
STOP
END