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
|
subroutine ggaec(rs,zet,t,ec,h)
c called by subroutine ggaexc
c gga91 correlation
c input rs : seitz radius
c input zet : relative spin polarization
c input t : abs(grad d)/(d*2.*ks*g)
c input : correlation energy per electron (ec)
c output h : nonlocal part of correlation energy per electron
implicit double precision (a-h,o-z)
data xnu,cc0,cx,alf/15.75592d0,0.004235d0,-0.001667212d0,0.09d0/
data c1,c2,c3,c4/0.002568d0,0.023266d0,7.389d-6,8.723d0/
data c5,c6,a4/0.472d0,7.389d-2,100.d0/
data thrd2/0.666666666667d0/
pi = 4.d0*datan(1.d0)
fk = 1.91915829d0/rs
sk = dsqrt(4.d0*fk/pi)
g = ((1.d0+zet)**thrd2+(1.d0-zet)**thrd2)/2.d0
bet = xnu*cc0
delt = 2.d0*alf/bet
g3 = g**3
g4 = g3*g
pon = -delt*ec/(g3*bet)
b = delt/(dexp(pon)-1.d0)
b2 = b*b
t2 = t*t
t4 = t2*t2
rs2 = rs*rs
rs3 = rs2*rs
q4 = 1.d0+b*t2
q5 = 1.d0+b*t2+b2*t4
q6 = c1+c2*rs+c3*rs2
q7 = 1.d0+c4*rs+c5*rs2+c6*rs3
cc = -cx + q6/q7
r0 = (sk/fk)**2
r1 = a4*r0*g4
coeff = cc-cc0-3.d0*cx/7.d0
r2 = xnu*coeff*g3
r3 = dexp(-r1*t2)
h0 = g3*(bet/delt)*dlog(1.d0+delt*q4*t2/q5)
h1 = r3*r2*t2
h = h0 + h1
return
end
|
