exchpbe.f90 6.49 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
MODULE m_exchpbe
!----------------------------------------------------------------------
!     pbe exchange for a spin-unpolarized electronic system
!     k burke's modification of pw91 codes, may 14, 1996
!     modified again by k. burke, june 29, 1996, with simpler fx(s)
!     inclusion of HSE function by M. Schlipf 2009
!----------------------------------------------------------------------
!     references:
!     [a]j.p.~perdew, k.~burke, and m.~ernzerhof, submiited to prl, may96
!     [b]j.p. perdew and y. wang, phys. rev.  b {\bf 33},  8800  (1986);
!     {\bf 40},  3399  (1989) (e).
!     [c] B.~Hammer, L.~B.~Hansen and J.~K.~Norskov PRB 59 7413 (1999)
!     [d] J.~Heyd, G.~E.~Scuseria, M.~Ernzerhof, J. Chem. Phys. {\bf 118},
!     8207 (2003)
!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE exchpbe(xcpot,rho,s,u,v,lgga,lpot, &
                      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 24
      IMPLICIT NONE

25
!     .. Arguments
26
      TYPE(t_xcpot_data), INTENT (IN)  :: xcpot
27 28
      INTEGER, INTENT (IN)  :: lgga ! =0=>don't put in gradient corrections, just lda
      INTEGER, INTENT (IN)  :: lpot ! =0=>don't get potential and don't need u and v
29 30 31 32
      REAL,    INTENT (IN)  :: rho ! density
      REAL,    INTENT (IN)  :: s ! abs(grad rho)/(2*kf*rho), where kf=(3 pi^2 rho)^(1/3)
      REAL,    INTENT (IN)  :: u ! (grad rho)*grad(abs(grad rho))/(rho**2 * (2*kf)**3)
      REAL,    INTENT (IN)  :: v ! (laplacian rho)/(rho*(2*kf)**2) [pw86(24)]
33 34
      REAL,    INTENT (OUT) :: ex,vx ! exchange energy per electron (ex) and potential (vx)

35
!     new variables for the HSE functional
36 37 38 39
      REAL,    INTENT (OUT) :: vx_sr ! short ranged part of potential
      REAL :: kF,factor,fxhse
      REAL :: dFx_ds,d2Fx_ds2,dFx_dkF,d2Fx_dsdkF

40
!     .. local variables ..
41
      REAL :: ul,exunif,fs,fss,fxpbe,p0,s2
42 43 44
      REAL :: xwu,css,dxwu,ddx  ! for wu-cohen
      REAL, PARAMETER :: teo = 10.e0/81.e0 ! for wu-cohen
      REAL, PARAMETER :: cwu = 0.0079325 ! for wu-cohen
45 46
      REAL, PARAMETER :: thrd=1.e0/3.e0
      REAL, PARAMETER :: thrd4=4.e0/3.e0
47
      REAL, PARAMETER :: ax=-0.738558766382022405884230032680836e0 ! -0.75*(3/pi)^(1/3)
48

Matthias Redies's avatar
Matthias Redies committed
49 50 51
!----------------------------------------------------------------------
!     uk, ul defined after [a](13)  (for icorr==7)
!----------------------------------------------------------------------
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
!      IF (xcpot%is_name("pbe").OR.xcpot%is_name("wc").OR.
!     +     xcpot%is_name("PBEs")) THEN
!         uk=0.8040
!      ELSEIF (xcpot%is_name("rpbe")) THEN
!         uk=1.2450
!      ELSEIF (xcpot%is_name("Rpbe")) THEN  ! changed to [c]
!         uk=0.8040
!      ELSEIF (xcpot%is_name("pbe0").OR.xcpot%is_name("hse")
!     +        .OR. xcpot%is_name("lhse").OR.xcpot%is_name("vhse")) THEN
!         uk=0.8040
!      ELSE
!         CALL judft_error("BUG:Wrong xcpot",calledby="exchpbe.F")
!      ENDIF
!      IF (xcpot%is_name("PBEs")) THEN     ! pbe_sol
!         um=0.123456790123456d0
!      ELSE
!         um=0.2195149727645171e0
!      ENDIF
70

71
      ul=xcpot%um/xcpot%uk
Matthias Redies's avatar
Matthias Redies committed
72 73 74 75
!----------------------------------------------------------------------
!     construct lda exchange energy density:
!     e_x[unif] = -0.75 * (3/pi)^(1/3) * rho^(4/3)
!----------------------------------------------------------------------
76
      exunif = ax*rho**thrd
Matthias Redies's avatar
Matthias Redies committed
77
      IF (lgga == 0) THEN
78 79 80
         ex = exunif
         vx = ex*thrd4
         RETURN
81
      ENDIF
Matthias Redies's avatar
Matthias Redies committed
82 83 84 85 86
!----------------------------------------------------------------------
!     construct pbe enhancement factor
!     e_x[pbe]=e_x[unif]*fxpbe(s)
!     fxpbe(s)=1+xcpot%uk-xcpot%uk/(1+ul*s*s)                 [a](13)
!----------------------------------------------------------------------
87
      s2 = s*s
Matthias Redies's avatar
Matthias Redies committed
88

89
!     calculate fxpbe
90 91 92
      p0 = 1.e0 + ul*s2
      fxpbe = 1e0 + xcpot%uk - xcpot%uk/p0
      IF (xcpot%is_Rpbe) THEN
93
         p0 = exp( - ul*s2 )
94 95
         fxpbe = 1e0 + xcpot%uk * ( 1e0 - p0 )
      ELSEIF (xcpot%is_wc) THEN
96
         css = 1+cwu*s2*s2
97 98 99
         xwu = teo*s2 + (xcpot%um-teo)*s2*exp(-s2) + log(css)
         p0 = 1.e0 + xwu/xcpot%uk
         fxpbe = 1e0 + xcpot%uk - xcpot%uk/p0
100
      ENDIF
Matthias Redies's avatar
Matthias Redies committed
101
!     -gu
102
!     Mixing of short and long range components
103
      IF (xcpot%is_hse) THEN
Matthias Redies's avatar
Matthias Redies committed
104 105 106 107 108
         !     ex = (1-a)ex,SR + ex,LR
         !     = (1-a)ex,SR + (ex,PBE - ex,SR)
         !     = (fxpbe - a*fxhse)*exunif
         !     Calculate the enhancement factor fxhse and its derivatives
         !     as integral over the exchange hole (cf. [d])
109
         kF = (3.0 * pi_const**2 * rho)**thrd
110
         CALL judft_error("HSE not implemented",calledby="exchpbe")
Matthias Redies's avatar
Matthias Redies committed
111 112 113
         ! his creates a depency loop
         !     CALL calculateEnhancementFactor(kF, s, fxhse, dFx_ds, d2Fx_ds2,
         !     &        dFx_dkF, d2Fx_dsdkF)
114
         ex = exunif * (fxpbe - xcpot%exchange_weight * fxhse )
115
      ELSE
116
         ex = exunif*fxpbe
117
      END IF
Matthias Redies's avatar
Matthias Redies committed
118 119 120 121 122 123 124
      IF (lpot == 0) RETURN
!----------------------------------------------------------------------
!     energy done. now the potential:
!     find first and second derivatives of fx w.r.t s.
!     fs=(1/s)*d fxpbe/ ds
!     fss=d fs/ds
!----------------------------------------------------------------------
125
!     derivatives for the pbe part
126 127 128
      fs = 2.e0*xcpot%uk*ul/ (p0*p0)
      fss = -4.e0*ul*s*fs/p0
      IF (xcpot%is_Rpbe) THEN
129 130
         fs = 2.e0*ul*p0
         fss = -2.e0*ul*s*fs
131 132
      ELSEIF (xcpot%is_wc) THEN
         dxwu = 2*teo + 2*(xcpot%um-teo)*exp(-s2)*(1-s2) + 4*cwu*s2/css
133
         fs = dxwu / (p0*p0)
Matthias Redies's avatar
Matthias Redies committed
134 135
         ddx = 4*s*((xcpot%um-teo)*exp(-s2)*(s2-2)+2*cwu* &
                    (1-cwu*s2*s2)/css**2)
136
         fss = ( ddx - 2*s*dxwu*dxwu/(p0*xcpot%uk) ) / (p0*p0)
137
      ENDIF
Matthias Redies's avatar
Matthias Redies committed
138 139 140 141
!     -gu
!----------------------------------------------------------------------
!     calculate potential from [b](24)
!----------------------------------------------------------------------
142 143
      vx = exunif* (thrd4*fxpbe- (u-thrd4*s2*s)*fss-v*fs)

Matthias Redies's avatar
Matthias Redies committed
144
      IF ( .NOT. (xcpot%is_hse)) RETURN
145

Matthias Redies's avatar
Matthias Redies committed
146 147 148
!----------------------------------------------------------------------
!     short ranged potential (HSE functional)
!----------------------------------------------------------------------
149 150
!     calculate fs and fss for the HSE functional
!     where the 1st and 2nd derivative of Fx are known
151 152
      fs    = dFx_ds / s
      fss   = (d2Fx_ds2 - fs) / s
Matthias Redies's avatar
Matthias Redies committed
153 154
      vx_sr = exunif * ( thrd4*fxhse - (u-thrd4*s2*s)*fss - v*fs &
                        + thrd*kF * ( dFx_dkF - d2Fx_dsdkF*s ) )
155

Matthias Redies's avatar
Matthias Redies committed
156 157
   END SUBROUTINE exchpbe
END MODULE m_exchpbe