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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
|
subroutine ggauc(rs,zet,t,uu,vv,ww,ec,ecrs,eczet,dvcup,dvcdn)
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 uu : (grad d)*grad(abs(grad d))/(d**2 * (2*ks*g)**3)
c input vv : (laplacian d)/(d * (2*ks*g)**2)
c input ww : (grad d)*(grad zet)/(d * (2*ks*g)**2
c input ec : correlation energy
c input ecrs : derivative of ec with respect to rs
c input eczet : derivative of ec with respect to zet
c output dvcup,dvcdn : nonlocal parts of correlation potentials
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 thrdm,thrd2/-0.333333333333d0,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
t6 = t4*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
c energy done. now the potential:
ccrs = (c2+2.*c3*rs)/q7 - q6*(c4+2.*c5*rs+3.*c6*rs2)/q7**2
rsthrd = rs/3.d0
r4 = rsthrd*ccrs/coeff
c=========================================================
c fix made by Tomas Holmquist
if ((zet.eq.1.d0).or.(zet.eq.-1.d0)) then
gz = 0.d0
else
gz = ((1.d0+zet)**thrdm - (1.d0-zet)**thrdm)/3.d0
endif
c=========================================================
fac = delt/b+1.d0
bg = -3.d0*b2*ec*fac/(bet*g4)
bec = b2*fac/(bet*g3)
q8 = q5*q5+delt*q4*q5*t2
q9 = 1.d0+2.d0*b*t2
h0b = -bet*g3*b*t6*(2.d0+b*t2)/q8
h0rs = -rsthrd*h0b*bec*ecrs
fact0 = 2.d0*delt-6.d0*b
fact1 = q5*q9+q4*q9*q9
h0bt = 2.d0*bet*g3*t4*((q4*q5*fact0-delt*fact1)/q8)/q8
h0rst = rsthrd*t2*h0bt*bec*ecrs
h0z = 3.d0*gz*h0/g + h0b*(bg*gz+bec*eczet)
h0t = 2.*bet*g3*q9/q8
h0zt = 3.d0*gz*h0t/g+h0bt*(bg*gz+bec*eczet)
fact2 = q4*q5+b*t2*(q4*q9+q5)
fact3 = 2.d0*b*q5*q9+delt*fact2
h0tt = 4.d0*bet*g3*t*(2.d0*b/q8-(q9*fact3/q8)/q8)
h1rs = r3*r2*t2*(-r4+r1*t2/3.d0)
fact4 = 2.d0-r1*t2
h1rst = r3*r2*t2*(2.d0*r4*(1.d0-r1*t2)-thrd2*r1*t2*fact4)
h1z = gz*r3*r2*t2*(3.d0-4.d0*r1*t2)/g
h1t = 2.d0*r3*r2*(1.d0-r1*t2)
h1zt = 2.d0*gz*r3*r2*(3.d0-11.d0*r1*t2+4.d0*r1*r1*t4)/g
h1tt = 4.d0*r3*r2*r1*t*(-2.d0+r1*t2)
hrs = h0rs+h1rs
hrst = h0rst+h1rst
ht = h0t+h1t
htt = h0tt+h1tt
hz = h0z+h1z
hzt = h0zt+h1zt
comm = h+hrs+hrst+t2*ht/6.d0+7.d0*t2*t*htt/6.d0
pref = hz-gz*t2*ht/g
fact5 = gz*(2.d0*ht+t*htt)/g
comm = comm-pref*zet-uu*htt-vv*ht-ww*(hzt-fact5)
dvcup = comm + pref
dvcdn = comm - pref
return
end
|
