Line data Source code
1 : MODULE m_exchpbe
2 : !----------------------------------------------------------------------
3 : ! pbe exchange for a spin-unpolarized electronic system
4 : ! k burke's modification of pw91 codes, may 14, 1996
5 : ! modified again by k. burke, june 29, 1996, with simpler fx(s)
6 : ! inclusion of HSE function by M. Schlipf 2009
7 : !----------------------------------------------------------------------
8 : ! references:
9 : ! [a]j.p.~perdew, k.~burke, and m.~ernzerhof, submiited to prl, may96
10 : ! [b]j.p. perdew and y. wang, phys. rev. b {\bf 33}, 8800 (1986);
11 : ! {\bf 40}, 3399 (1989) (e).
12 : ! [c] B.~Hammer, L.~B.~Hansen and J.~K.~Norskov PRB 59 7413 (1999)
13 : ! [d] J.~Heyd, G.~E.~Scuseria, M.~Ernzerhof, J. Chem. Phys. {\bf 118},
14 : ! 8207 (2003)
15 : !----------------------------------------------------------------------
16 : CONTAINS
17 614387492 : SUBROUTINE exchpbe(xcpot,rho,s,u,v,lgga,lpot, &
18 : ex,vx,vx_sr)
19 : USE m_hsefunctional, ONLY: calculateEnhancementFactor
20 : USE m_constants, ONLY: pi_const
21 : USE m_types_xcpot_data
22 : USE m_judft
23 : IMPLICIT NONE
24 :
25 : ! .. Arguments
26 : TYPE(t_xcpot_data), INTENT (IN) :: xcpot
27 : INTEGER, INTENT (IN) :: lgga ! =0=>don't put in gradient corrections, just lda
28 : INTEGER, INTENT (IN) :: lpot ! =0=>don't get potential and don't need u and v
29 : REAL, INTENT (IN) :: rho ! density
30 : REAL, INTENT (IN) :: s ! abs(grad rho)/(2*kf*rho), where kf=(3 pi^2 rho)^(1/3)
31 : REAL, INTENT (IN) :: u ! (grad rho)*grad(abs(grad rho))/(rho**2 * (2*kf)**3)
32 : REAL, INTENT (IN) :: v ! (laplacian rho)/(rho*(2*kf)**2) [pw86(24)]
33 : REAL, INTENT (OUT) :: ex,vx ! exchange energy per electron (ex) and potential (vx)
34 :
35 : ! new variables for the HSE functional
36 : REAL, INTENT (OUT) :: vx_sr ! short ranged part of potential
37 : REAL :: kF,factor,fxhse
38 : REAL :: dFx_ds,d2Fx_ds2,dFx_dkF,d2Fx_dsdkF
39 :
40 : ! .. local variables ..
41 : REAL :: ul,exunif,fs,fss,fxpbe,p0,s2
42 : REAL :: xwu,css,dxwu,ddx ! for wu-cohen
43 : REAL, PARAMETER :: teo = 10.e0/81.e0 ! for wu-cohen
44 : REAL, PARAMETER :: cwu = 0.0079325 ! for wu-cohen
45 : REAL, PARAMETER :: thrd=1.e0/3.e0
46 : REAL, PARAMETER :: thrd4=4.e0/3.e0
47 : REAL, PARAMETER :: ax=-0.738558766382022405884230032680836e0 ! -0.75*(3/pi)^(1/3)
48 :
49 : !----------------------------------------------------------------------
50 : ! uk, ul defined after [a](13) (for icorr==7)
51 : !----------------------------------------------------------------------
52 : ! IF (xcpot%is_name("pbe").OR.xcpot%is_name("wc").OR.
53 : ! + xcpot%is_name("PBEs")) THEN
54 : ! uk=0.8040
55 : ! ELSEIF (xcpot%is_name("rpbe")) THEN
56 : ! uk=1.2450
57 : ! ELSEIF (xcpot%is_name("Rpbe")) THEN ! changed to [c]
58 : ! uk=0.8040
59 : ! ELSEIF (xcpot%is_name("pbe0").OR.xcpot%is_name("hse")
60 : ! + .OR. xcpot%is_name("lhse").OR.xcpot%is_name("vhse")) THEN
61 : ! uk=0.8040
62 : ! ELSE
63 : ! CALL judft_error("BUG:Wrong xcpot",calledby="exchpbe.F")
64 : ! ENDIF
65 : ! IF (xcpot%is_name("PBEs")) THEN ! pbe_sol
66 : ! um=0.123456790123456d0
67 : ! ELSE
68 : ! um=0.2195149727645171e0
69 : ! ENDIF
70 :
71 614387492 : ul=xcpot%um/xcpot%uk
72 : !----------------------------------------------------------------------
73 : ! construct lda exchange energy density:
74 : ! e_x[unif] = -0.75 * (3/pi)^(1/3) * rho^(4/3)
75 : !----------------------------------------------------------------------
76 614387492 : exunif = ax*rho**thrd
77 614387492 : IF (lgga == 0) THEN
78 307193746 : ex = exunif
79 307193746 : vx = ex*thrd4
80 614387492 : RETURN
81 : ENDIF
82 : !----------------------------------------------------------------------
83 : ! construct pbe enhancement factor
84 : ! e_x[pbe]=e_x[unif]*fxpbe(s)
85 : ! fxpbe(s)=1+xcpot%uk-xcpot%uk/(1+ul*s*s) [a](13)
86 : !----------------------------------------------------------------------
87 307193746 : s2 = s*s
88 :
89 : ! calculate fxpbe
90 307193746 : p0 = 1.e0 + ul*s2
91 307193746 : fxpbe = 1e0 + xcpot%uk - xcpot%uk/p0
92 307193746 : IF (xcpot%is_Rpbe) THEN
93 0 : p0 = exp( - ul*s2 )
94 0 : fxpbe = 1e0 + xcpot%uk * ( 1e0 - p0 )
95 307193746 : ELSEIF (xcpot%is_wc) THEN
96 5773860 : css = 1+cwu*s2*s2
97 5773860 : IF(s2.GT.100) THEN ! This is introduced because the PGI compiler has problems with calculating exp(-somethingLarge).
98 569796 : xwu = teo*s2 + log(css)
99 : ELSE
100 5204064 : xwu = teo*s2 + (xcpot%um-teo)*s2*exp(-s2) + log(css)
101 : END IF
102 5773860 : p0 = 1.e0 + xwu/xcpot%uk
103 5773860 : fxpbe = 1e0 + xcpot%uk - xcpot%uk/p0
104 : ENDIF
105 : ! -gu
106 : ! Mixing of short and long range components
107 307193746 : IF (xcpot%is_hse) THEN
108 : ! ex = (1-a)ex,SR + ex,LR
109 : ! = (1-a)ex,SR + (ex,PBE - ex,SR)
110 : ! = (fxpbe - a*fxhse)*exunif
111 : ! Calculate the enhancement factor fxhse and its derivatives
112 : ! as integral over the exchange hole (cf. [d])
113 0 : kF = (3.0 * pi_const**2 * rho)**thrd
114 : !! CALL judft_error("HSE not implemented",calledby="exchpbe")
115 : ! his creates a depency loop
116 : !CALL timestart("hse: calc enhancement factor")
117 : CALL calculateEnhancementFactor(kF, s, fxhse, dFx_ds, d2Fx_ds2, &
118 0 : dFx_dkF, d2Fx_dsdkF)
119 : !CALL timestop("hse: calc enhancement factor")
120 0 : ex = exunif * (fxpbe - xcpot%exchange_weight * fxhse )
121 : ELSE
122 307193746 : ex = exunif*fxpbe
123 : END IF
124 307193746 : IF (lpot == 0) RETURN
125 : !----------------------------------------------------------------------
126 : ! energy done. now the potential:
127 : ! find first and second derivatives of fx w.r.t s.
128 : ! fs=(1/s)*d fxpbe/ ds
129 : ! fss=d fs/ds
130 : !----------------------------------------------------------------------
131 : ! derivatives for the pbe part
132 307193746 : fs = 2.e0*xcpot%uk*ul/ (p0*p0)
133 307193746 : fss = -4.e0*ul*s*fs/p0
134 307193746 : IF (xcpot%is_Rpbe) THEN
135 0 : fs = 2.e0*ul*p0
136 0 : fss = -2.e0*ul*s*fs
137 307193746 : ELSEIF (xcpot%is_wc) THEN
138 5773860 : IF (s2.GT.100) THEN ! This is introduced because the PGI compiler has problems with calculating exp(-somethingLarge).
139 569796 : dxwu = 2*teo + 4*cwu*s2/css
140 569796 : fs = dxwu / (p0*p0)
141 569796 : ddx = 4*s*(2*cwu*(1-cwu*s2*s2)/css**2)
142 : ELSE
143 5204064 : dxwu = 2*teo + 2*(xcpot%um-teo)*exp(-s2)*(1-s2) + 4*cwu*s2/css
144 5204064 : fs = dxwu / (p0*p0)
145 : ddx = 4*s*((xcpot%um-teo)*exp(-s2)*(s2-2)+2*cwu* &
146 5204064 : (1-cwu*s2*s2)/css**2)
147 : END IF
148 5773860 : fss = ( ddx - 2*s*dxwu*dxwu/(p0*xcpot%uk) ) / (p0*p0)
149 : ENDIF
150 : ! -gu
151 : !----------------------------------------------------------------------
152 : ! calculate potential from [b](24)
153 : !----------------------------------------------------------------------
154 307193746 : vx = exunif* (thrd4*fxpbe- (u-thrd4*s2*s)*fss-v*fs)
155 :
156 307193746 : IF ( .NOT. (xcpot%is_hse)) RETURN
157 :
158 : !----------------------------------------------------------------------
159 : ! short ranged potential (HSE functional)
160 : !----------------------------------------------------------------------
161 : ! calculate fs and fss for the HSE functional
162 : ! where the 1st and 2nd derivative of Fx are known
163 0 : fs = dFx_ds / s
164 0 : fss = (d2Fx_ds2 - fs) / s
165 : vx_sr = exunif * ( thrd4*fxhse - (u-thrd4*s2*s)*fss - v*fs &
166 0 : + thrd*kF * ( dFx_dkF - d2Fx_dsdkF*s ) )
167 :
168 : END SUBROUTINE exchpbe
169 : END MODULE m_exchpbe
|
