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
|
subroutine ggauxc(r,db,zet,i,imin,imax,donela,donegr,donez,
j uxcup,uxcdn)
c ggauxc calculates the exchange-correlation potential by taking
c the functional derivative of the energy with respect to the density
c input r : position coordinate (array)
c input db : spin up and down density. (array)
c input zet : relative spin polarization (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 donela : true if r^2*grad(density) is splined
c in/output donegr : true if abs(grad(density)) is splined
c in/output donez : true if zet is splined
c output uxcup : spin-up exchange-correlation potential
c output uxcdn : spin-down exchange-correlation potential
implicit logical (a-z)
double precision r,db,zet,uxcup,uxcdn,gradd,gradup,graddn
double precision fkup,fkdn,sk,g,sup,sdn,t,tempz,z2,gradz
double precision onethi,pi,d,rs,lapld,laplup,lapldn,deriv
double precision gagd,gagup,gagdn,uu,vv,ww,uup,vup,udn,vdn
double precision uclcup,uclcdn,ucgaup,ucgadn,uxup,uxdn
double precision ucup,ucdn,ec,ecrs,eczet
integer i,imin,imax,j
logical donela,donegr,donez
common/gga/gradd,gradup,graddn,fkup,fkdn,sk,g,sup,sdn,t
common/ggaz/tempz,z2
dimension r(561),db(2,561),zet(561),tempz(2,561),z2(2,561)
dimension gradd(561),gradup(561),graddn(561),fkup(561)
dimension fkdn(561),sk(561),g(561),sup(561),sdn(561),t(561)
onethi = 1.d0/3.d0
pi = 4.d0*datan(1.d0)
call laplac(r,gradup,graddn,i,imin,imax,donela,lapld,laplup,
j lapldn)
call grabgr(r,gradup,graddn,i,imin,imax,donegr,gagd,gagup,gagdn)
c calculate the exchange potential
if (db(1,i).gt.1.d-100) then
uup = gradup(i)*gagup/(db(1,i)**2*(2.d0*fkup(i))**3)
vup = laplup/(db(1,i)*(2.d0*fkup(i))**2)
call exchpt(2.d0*db(1,i),sup(i),uup,vup,uxup)
else
uxup = 0.d0
endif
if (db(2,i).gt.1.d-100) then
udn = graddn(i)*gagdn/(db(2,i)**2*(2.d0*fkdn(i))**3)
vdn = lapldn/(db(2,i)*(2.d0*fkdn(i))**2)
call exchpt(2.d0*db(2,i),sdn(i),udn,vdn,uxdn)
else
uxdn = 0.d0
endif
c exchange potential done, now calculate the correlation potential
if (.not.donez) then
do 10 j = imin,imax
tempz(1,j) = zet(j)
10 continue
call spline(r,imin+1,imax,tempz,z2)
donez = .true.
endif
gradz = deriv(r,tempz,z2,1,i,imin,imax)
d = db(1,i) + db(2,i)
if (d.gt.1.d-18) then
rs = (3.d0/(4.d0*pi*d))**onethi
uu = gradd(i)*gagd/(d**2*(2.d0*sk(i)*g(i))**3)
vv = lapld/(d*(2.d0*sk(i)*g(i))**2)
ww = gradd(i)*gradz/(d*(2.d0*sk(i)*g(i))**2)
call ldauc(rs,zet(i),ec,ecrs,eczet,uclcup,uclcdn)
call ggauc(rs,zet(i),t(i),uu,vv,ww,ec,ecrs,eczet,ucgaup,
j ucgadn)
ucup = uclcup + ucgaup
ucdn = uclcdn + ucgadn
else
ucdn = 0.d0
ucup = 0.d0
endif
uxcup = uxup + ucup
uxcdn = uxdn + ucdn
return
end
|
