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

c----------------------------------------------------------------------
50
c     uk, ul defined after [a](13)  (for icorr==7)
51
c----------------------------------------------------------------------
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 72
      
      ul=xcpot%um/xcpot%uk
73
c----------------------------------------------------------------------
74 75
c     construct lda exchange energy density: 
!     e_x[unif] = -0.75 * (3/pi)^(1/3) * rho^(4/3)  
76 77 78
c----------------------------------------------------------------------
      exunif = ax*rho**thrd
      IF (lgga.EQ.0) THEN
79 80 81
         ex = exunif
         vx = ex*thrd4
         RETURN
82 83
      ENDIF
c----------------------------------------------------------------------
84 85
c     construct pbe enhancement factor
c     e_x[pbe]=e_x[unif]*fxpbe(s)
86
c     fxpbe(s)=1+xcpot%uk-xcpot%uk/(1+ul*s*s)                 [a](13)
87 88
c----------------------------------------------------------------------
      s2 = s*s
89
      
90
!     calculate fxpbe
91 92 93
      p0 = 1.e0 + ul*s2
      fxpbe = 1e0 + xcpot%uk - xcpot%uk/p0
      IF (xcpot%is_Rpbe) THEN
94
         p0 = exp( - ul*s2 )
95 96
         fxpbe = 1e0 + xcpot%uk * ( 1e0 - p0 )
      ELSEIF (xcpot%is_wc) THEN
97
         css = 1+cwu*s2*s2
98 99 100
         xwu = teo*s2 + (xcpot%um-teo)*s2*exp(-s2) + log(css)
         p0 = 1.e0 + xwu/xcpot%uk
         fxpbe = 1e0 + xcpot%uk - xcpot%uk/p0
101
      ENDIF
102 103
c     -gu
!     Mixing of short and long range components
104
      IF (xcpot%is_hse) THEN
105 106 107 108 109 110
!     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])
         kF = (3.0 * pi_const**2 * rho)**thrd
111 112 113 114 115
         CALL judft_error("HSE not implemented",calledby="exchpbe")
         !This creates a depency loop
!     CALL calculateEnhancementFactor(kF, s, fxhse, dFx_ds, d2Fx_ds2,
!     &        dFx_dkF, d2Fx_dsdkF)
         ex = exunif * (fxpbe - xcpot%exchange_weight * fxhse )
116
      ELSE
117
         ex = exunif*fxpbe
118 119 120
      END IF
      IF (lpot.EQ.0) RETURN
c----------------------------------------------------------------------
121 122 123 124
c     energy done. now the potential:
c     find first and second derivatives of fx w.r.t s.
c     fs=(1/s)*d fxpbe/ ds
c     fss=d fs/ds
125
c----------------------------------------------------------------------
126
!     derivatives for the pbe part
127 128 129
      fs = 2.e0*xcpot%uk*ul/ (p0*p0)
      fss = -4.e0*ul*s*fs/p0
      IF (xcpot%is_Rpbe) THEN
130 131
         fs = 2.e0*ul*p0
         fss = -2.e0*ul*s*fs
132 133
      ELSEIF (xcpot%is_wc) THEN
         dxwu = 2*teo + 2*(xcpot%um-teo)*exp(-s2)*(1-s2) + 4*cwu*s2/css
134
         fs = dxwu / (p0*p0)
135 136 137
         ddx = 4*s*((xcpot%um-teo)*exp(-s2)*(s2-2)+2*cwu*
     &        (1-cwu*s2*s2)/css**2)
         fss = ( ddx - 2*s*dxwu*dxwu/(p0*xcpot%uk) ) / (p0*p0)
138
      ENDIF
139
c     -gu
140
c----------------------------------------------------------------------
141
c     calculate potential from [b](24)
142 143 144
c----------------------------------------------------------------------
      vx = exunif* (thrd4*fxpbe- (u-thrd4*s2*s)*fss-v*fs)

145
      IF (.not.(xcpot%is_hse)) RETURN
146 147 148 149

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

      END SUBROUTINE exchpbe
158
      END