1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
subroutine ggaexc(r,db,i,imin,imax,doned,exc)
c ggaexc calculates the exchange-correlation energy
c input r : position coordinate (array)
c input db : spin up and down density. (array)
c input i : counter, i.e. present position = r(i)
c input imin : min value of i
c input imax : max value of i
c in/output doned : true if density is splined
c output exc : exchange-correlation energy
implicit logical (a-z)
logical doned
double precision r,db,gradd,gradup,graddn,fk,fkup,fkdn,sk,g
double precision sup,sdn,t,exc,ec
double precision onethi,twothi,pi,conkf,d,zet,rs,exup,exdn
double precision ex,eclda,ecgga
integer i,imin,imax
dimension r(561),db(2,561)
dimension gradd(561),gradup(561),graddn(561),fk(561),fkup(561)
dimension fkdn(561),sk(561),g(561),sup(561),sdn(561),t(561)
common/gga/gradd,gradup,graddn,fkup,fkdn,sk,g,sup,sdn,t
onethi = 1.d0/3.d0
twothi = 2.d0/3.d0
pi = 4.d0*datan(1.d0)
conkf = (3.d0*pi**2)**onethi
d = db(1,i) + db(2,i)
call grad(r,db,i,imin,imax,doned,gradd(i),gradup(i),graddn(i))
c calculate the exchange energy, first spin up, then spin down
if (db(1,i).gt.1.d-100) then
fkup(i) = conkf*(2.d0*db(1,i))**onethi
sup(i) = dabs(gradup(i))/(2.d0*fkup(i)*db(1,i))
call exchen(2.d0*db(1,i),sup(i),exup)
else
fkup(i) = 0.d0
sup(i) = 0.d0
exup = 0.d0
endif
if (db(2,i).gt.1.d-100) then
fkdn(i) = conkf*(2.d0*db(2,i))**onethi
sdn(i) = dabs(graddn(i))/(2.d0*fkdn(i)*db(2,i))
call exchen(2.d0*db(2,i),sdn(i),exdn)
else
fkdn(i) = 0.d0
sdn(i) = 0.d0
exdn = 0.d0
endif
if (d.gt.1d-100) then
ex = (exup*db(1,i) + exdn*db(2,i))/d
else
ex = 0.d0
endif
c exchange energy done, now calculate the correlation energy
fk(i) = conkf*d**onethi
sk(i) = dsqrt(4.d0*fk(i)/pi)
if (d.gt.1.d-100) then
zet = (db(1,i) - db(2,i))/d
else
zet = 0.d0
endif
g(i) = ((1.d0+zet)**twothi+(1.d0-zet)**twothi)/2.d0
if (d.gt.1.d-18) then
rs = (3.d0/(4.d0*pi*d))**onethi
t(i) = dabs(gradd(i))/(d*2*sk(i)*g(i))
call ldaec(rs,zet,eclda)
call ggaec(rs,zet,t(i),eclda,ecgga)
ec = eclda + ecgga
else
ec = 0.d0
endif
exc = ex + ec
return
end
|
