Line data Source code
1 : MODULE m_corpbe
2 :
3 : !----------------------------------------------------------------------
4 : ! official pbe correlation code. k. burke, may 14, 1996.
5 : !----------------------------------------------------------------------
6 : ! references:
7 : ! [a] j.p.~perdew, k.~burke, and m.~ernzerhof,
8 : ! {\sl generalized gradient approximation made simple}, sub.
9 : ! to phys. rev.lett. may 1996.
10 : ! [b] j. p. perdew, k. burke, and y. wang, {\sl real-space cutoff
11 : ! construction of a generalized gradient approximation: the pw91
12 : ! density functional}, submitted to phys. rev. b, feb. 1996.
13 : ! [c] j. p. perdew and y. wang, phys. rev. b {\bf 45}, 13244 (1992).
14 : !----------------------------------------------------------------------
15 : CONTAINS
16 153598753 : SUBROUTINE corpbe( &
17 : l_pbes,rs,zet,t,uu,vv,ww,lgga,lpot, &
18 : ec,vcup,vcdn,h,dvcup,dvcdn)
19 : !----------------------------------------------------------------------
20 : ! input: rs=seitz radius=(3/4pi rho)^(1/3)
21 : ! : zet=relative spin polarization = (rhoup-rhodn)/rho
22 : ! : t=abs(grad rho)/(rho*2.*ks*g) -- only needed for pbe
23 : ! : uu=(grad rho)*grad(abs(grad rho))/(rho**2 * (2*ks*g)**3)
24 : ! : vv=(laplacian rho)/(rho * (2*ks*g)**2)
25 : ! : ww=(grad rho)*(grad zet)/(rho * (2*ks*g)**2
26 : ! : uu,vv,ww, only needed for pbe potential
27 : ! : lgga=flag to do gga (0=>lsd only)
28 : ! : lpot=flag to do potential (0=>energy only)
29 : ! output: ec=lsd correlation energy from [a]
30 : ! : vcup=lsd up correlation potential
31 : ! : vcdn=lsd dn correlation potential
32 : ! : h=nonlocal part of correlation energy per electron
33 : ! : dvcup=nonlocal correction to vcup
34 : ! : dvcdn=nonlocal correction to vcdn
35 : !----------------------------------------------------------------------
36 :
37 : USE m_pbecor2
38 : IMPLICIT NONE
39 : LOGICAL,INTENT(IN) :: l_pbes
40 : INTEGER, INTENT (IN) :: lgga,lpot
41 : REAL, INTENT (IN) :: rs,zet,t,uu,vv,ww
42 : REAL, INTENT (OUT) :: dvcdn,dvcup,ec,h,vcdn,vcup
43 :
44 : ! thrd*=various multiples of 1/3
45 : ! numbers for use in lsd energy spin-interpolation formula, [c](9).
46 : ! gam= 2^(4/3)-2
47 : ! fzz=f''(0)= 8/(9*gam)
48 : ! numbers for construction of pbe
49 : ! gamma=(1-log(2))/pi^2
50 : ! bet=coefficient in gradient expansion for correlation, [a](4).
51 : ! eta=small number to stop d phi/ dzeta from blowing up at
52 : ! |zeta|=1.
53 :
54 : REAL, PARAMETER :: thrd=1.e0/3.e0
55 : REAL, PARAMETER :: thrdm=-thrd
56 : REAL, PARAMETER :: thrd2=2.e0*thrd
57 : REAL, PARAMETER :: sixthm=thrdm/2.e0
58 : REAL, PARAMETER :: thrd4=4.e0*thrd
59 : REAL, PARAMETER :: gam=0.5198420997897463295344212145565e0
60 : REAL, PARAMETER :: fzz=8.e0/ (9.e0*gam)
61 : REAL, PARAMETER :: gamma=0.03109069086965489503494086371273e0
62 : REAL, PARAMETER :: eta=1.e-12
63 : !----------------------------------------------------------------------
64 : !----------------------------------------------------------------------
65 : ! find lsd energy contributions, using [c](10) and table i[c].
66 : ! eu=unpolarized lsd correlation energy
67 : ! eurs=deu/drs
68 : ! ep=fully polarized lsd correlation energy
69 : ! eprs=dep/drs
70 : ! alfm=-spin stiffness, [c](3).
71 : ! alfrsm=-dalpha/drs
72 : ! f=spin-scaling factor from [c](9).
73 : ! construct ec, using [c](8)
74 :
75 : REAL :: alfm,alfrsm,b,b2,bec,bg,comm,ecrs,eczet,ep,eprs,eu,eurs, &
76 : f,fac,fact0,fact1,fact2,fact3,fact5,fz,g,g3,g4,gz,hb,hbt,hrs, &
77 : hrst,ht,htt,hz,hzt,pon,pref,q4,q5,q8,q9,rsthrd,rtrs, &
78 : t2,t4,t6,z4,delt,bet
79 : ! ..
80 :
81 153598753 : IF (l_pbes) THEN ! PBE_sol
82 : bet=0.046e0
83 : ELSE
84 153598753 : bet=0.06672455060314922e0
85 : ENDIF
86 153598753 : delt=bet/gamma
87 :
88 153598753 : rtrs = sqrt(rs)
89 : CALL pbecor2(0.0310907,0.21370,7.5957,3.5876,1.6382, &
90 153598753 : & 0.49294,rtrs,eu,eurs)
91 : CALL pbecor2(0.01554535,0.20548,14.1189,6.1977,3.3662, &
92 153598753 : & 0.62517,rtrs,ep,eprs)
93 : CALL pbecor2(0.0168869,0.11125,10.357,3.6231,0.88026, &
94 153598753 : & 0.49671,rtrs,alfm,alfrsm)
95 153598753 : z4 = zet**4
96 153598753 : f = ((1.e0+zet)**thrd4+ (1.e0-zet)**thrd4-2.e0)/gam
97 153598753 : ec = eu* (1.e0-f*z4) + ep*f*z4 - alfm*f* (1.e0-z4)/fzz
98 : !----------------------------------------------------------------------
99 : !----------------------------------------------------------------------
100 : ! lsd potential from [c](a1)
101 : ! ecrs = dec/drs [c](a2)
102 : ! eczet=dec/dzeta [c](a3)
103 : ! fz = df/dzeta [c](a4)
104 153598753 : ecrs = eurs* (1.e0-f*z4) + eprs*f*z4 - alfrsm*f* (1.e0-z4)/fzz
105 153598753 : fz = thrd4* ((1.e0+zet)**thrd- (1.e0-zet)**thrd)/gam
106 : eczet = 4.e0* (zet**3)*f* (ep-eu+alfm/fzz) + &
107 153598753 : fz* (z4*ep-z4*eu- (1.e0-z4)*alfm/fzz)
108 153598753 : comm = ec - rs*ecrs/3.e0 - zet*eczet
109 153598753 : vcup = comm + eczet
110 153598753 : vcdn = comm - eczet
111 153598753 : IF (lgga == 0) RETURN
112 : !----------------------------------------------------------------------
113 : !----------------------------------------------------------------------
114 : ! pbe correlation energy
115 : ! g=phi(zeta), given after [a](3)
116 : ! delt=bet/gamma
117 : ! b=a of [a](8)
118 153598753 : g = ((1.e0+zet)**thrd2+ (1.e0-zet)**thrd2)/2.e0
119 153598753 : g3 = g**3
120 153598753 : pon = -ec/ (g3*gamma)
121 153598753 : b = delt/ (exp(pon)-1.e0)
122 153598753 : b2 = b*b
123 153598753 : t2 = t*t
124 153598753 : t4 = t2*t2
125 153598753 : q4 = 1.e0 + b*t2
126 153598753 : q5 = 1.e0 + b*t2 + b2*t4
127 153598753 : h = g3* (bet/delt)*log(1.e0+delt*q4*t2/q5)
128 153598753 : IF (lpot == 0) RETURN
129 : !----------------------------------------------------------------------
130 : !----------------------------------------------------------------------
131 : ! energy done. now the potential, using appendix e of [b].
132 153598753 : g4 = g3*g
133 153598753 : t6 = t4*t2
134 153598753 : rsthrd = rs/3.e0
135 : gz = (((1.e0+zet)**2+eta)**sixthm- ((1.e0-zet)**2+eta)**sixthm)/ &
136 153598753 : & 3.e0
137 153598753 : fac = delt/b + 1.e0
138 153598753 : bg = -3.e0*b2*ec*fac/ (bet*g4)
139 153598753 : bec = b2*fac/ (bet*g3)
140 153598753 : q8 = q5*q5 + delt*q4*q5*t2
141 153598753 : q9 = 1.e0 + 2.e0*b*t2
142 153598753 : hb = -bet*g3*b*t6* (2.e0+b*t2)/q8
143 153598753 : hrs = -rsthrd*hb*bec*ecrs
144 153598753 : fact0 = 2.e0*delt - 6.e0*b
145 153598753 : fact1 = q5*q9 + q4*q9*q9
146 153598753 : hbt = 2.e0*bet*g3*t4* ((q4*q5*fact0-delt*fact1)/q8)/q8
147 153598753 : hrst = rsthrd*t2*hbt*bec*ecrs
148 153598753 : hz = 3.e0*gz*h/g + hb* (bg*gz+bec*eczet)
149 153598753 : ht = 2.e0*bet*g3*q9/q8
150 153598753 : hzt = 3.e0*gz*ht/g + hbt* (bg*gz+bec*eczet)
151 153598753 : fact2 = q4*q5 + b*t2* (q4*q9+q5)
152 153598753 : fact3 = 2.e0*b*q5*q9 + delt*fact2
153 153598753 : htt = 4.e0*bet*g3*t* (2.e0*b/q8- (q9*fact3/q8)/q8)
154 153598753 : comm = h + hrs + hrst + t2*ht/6.e0 + 7.e0*t2*t*htt/6.e0
155 153598753 : pref = hz - gz*t2*ht/g
156 153598753 : fact5 = gz* (2.e0*ht+t*htt)/g
157 153598753 : comm = comm - pref*zet - uu*htt - vv*ht - ww* (hzt-fact5)
158 153598753 : dvcup = comm + pref
159 153598753 : dvcdn = comm - pref
160 :
161 : END SUBROUTINE corpbe
162 : END MODULE m_corpbe
|
