C ************ TOTAL ENERGY CALCULATION FOR
atomlda.f ***********
c
for non-relativistic LDA cal.
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(651),E(29),ROSUM(651),X(651),R(651),
1
Wei(29),rv(651),ek(29),elu(29),h(651)
c --------------- input --------------------
open (1,file='atomlda.dat')
read (1,*) z
read (1,*) no
write (6,*)
'z',z,'no',no
do i =1 , no
read (1,*) npq
,elu(i) ,spin ,wei(i),ons
write (6,*) npq
,elu(i) , spin ,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, 20)
CLOSE (1)
OPEN (1,FILE='ROSUM.DAT')
READ (1,10) (ROSUM(I), I=1, 651)
CLOSE (1)
open (1,file='RSCFPO.DAT')
read (1,10) (rv(i), i=1, 651)
close (1)
open (1,file='KINETIC.DAT')
read (1,10) (ek(i), i=1, 20)
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 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)/2.0
scrnc=sum1
DO 700 I=2, 650, 2
SUM1=SUM1+dx/3.0*(DEXP(3.0*X(I-1))*ROSUM(I-1)+4.0*
1
DEXP(3.0*X(I))*ROSUM(I)+DEXP(3.0*X(I+1))*ROSUM(I+1))
700 CONTINUE
WRITE (6,710) SUM1
710 FORMAT (1H ,'TOTAL CHARGE=',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))*ROSUM(J-1)+4.0*
1
DEXP(3.0*X(J))*ROSUM(J)+DEXP(3.0*X(J+1))*ROSUM(J+1))
140 CONTINUE
130
SUM1=SUM1+dx/2.0*(DEXP(3.0*X(I-1))*ROSUM(I-1)+DEXP(3.0*X(I))*
1 ROSUM(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))*ROSUM(J-1)+4.0*
1
DEXP(2.0*X(J))*ROSUM(J)+DEXP(2.0*X(J+1))*ROSUM(J+1))
160 CONTINUE
150
SUM2=SUM2+dx/2.0*(DEXP(2.0*X(650))*ROSUM(650)+DEXP(2.0*X(651))
1 *ROSUM(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))*ROSUM(J-1)+4.0*
1
DEXP(3.0*X(J))*ROSUM(J)+DEXP(3.0*X(J+1))*ROSUM(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))*ROSUM(J-1)+4.0*
1
DEXP(2.0*X(J))*ROSUM(J)+DEXP(2.0*X(J+1))*ROSUM(J+1))
190 CONTINUE
GO TO 105
195 SUM2=0.0D0
105 F(I)=-0.5D0*(SUM1+SUM2)
100 CONTINUE
C ************ EXCHANGE-CORRELATION
CORRECTION ****************
pai=acos(-1.0)
DO 400 I=1, 651
IF (ROSUM(I)) 410,410,420
410 G(I)=0.0D0
GO TO 400
c ------- rosum=(4*pai)*ro(r)
-----------------
420 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
if (rs-1.0) 701,702,702
701 pc=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
ec=0.0311*log(rs)-0.048+0.0020*rs*log(rs)-0.0116*rs
go to 750
702 pc=-0.1423*(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)**2.0
ec=-0.1423/(1.0+1.0529*sqrt(rs)+0.3334*rs)
750 g(I)=ex+ec-(PC+PX)
h(i)=ex+ec
c write (6,1000) i,ex,ec,px,pc
1000 format(i5,4d15.7)
400 CONTINUE
C *********** INTEGRAL OVER THE WHOLE
SPHERE ****************
EINT=dx/2.0*DEXP(3.0*X(1))*(F(1)+G(1))*ROSUM(1)
EINT1=dx/2.0*DEXP(3.0*X(1))*(-F(1)+h(1)-z/r(1))*ROSUM(1)
DO 500 I=2, 650, 2
EINT=EINT+dx/3.0*(DEXP(3.0*X(I-1))*(F(I-1)+G(I-1))*ROSUM(I-1)+
1
4.0*DEXP(3.0*X(I))*(F(I)+G(I))*ROSUM(I)+DEXP(3.0*X(I+1))
2
*(F(I+1)+G(I+1))*ROSUM(I+1))
EINT1=EINT1+dx/3.0*(DEXP(3.0*X(I-1))*(-F(I-1)+h(I-1)-z/r(i-1))
&
*ROSUM(I-1)+4.0*DEXP(3.0*X(I))*(-F(I)+h(I)-z/r(i))*ROSUM(I)
&
+DEXP(3.0*X(I+1))*(-F(I+1)+h(I+1)-z/r(i+1))*ROSUM(I+1))
500 CONTINUE
C **************** TOTAL ENERGY
********************
ESUM=ESUM+EINT
eksum=eksum+eint1
OPEN (1,FILE='TOTE.DAT')
WRITE (1,610) ESUM,eksum
CLOSE (1)
WRITE (6,600) ESUM
600 FORMAT (1H ,'TOTAL E(AU) from eigenvalue=',1PD17.9)
write (6,620) eksum
620 format (1h ,'total E(au) from kinetic energy',1pd17.9)
610 FORMAT (2PD17.9)
write (6,'(f15.7)')
(esum+eksum)/2.0
STOP
END