exchpbe.F 6.88 KB
Newer Older
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 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
      MODULE m_exchpbe
c----------------------------------------------------------------------
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
c----------------------------------------------------------------------
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);
c     {\bf 40},  3399  (1989) (e).
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},
c     8207 (2003) 
c----------------------------------------------------------------------
      CONTAINS
      SUBROUTINE exchpbe(
     >                   icorr,rho,s,u,v,lgga,lpot,
     <                   ex,vx,vx_sr)
#ifndef CPP_GF
      USE m_hsefunctional, ONLY: calculateEnhancementFactor
      USE m_hybridmix
#endif
      USE m_constants,     ONLY: pi_const
      USE icorrkeys

      IMPLICIT NONE

      ! .. Arguments
      INTEGER, INTENT (IN)  :: icorr
      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
      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)]
      REAL,    INTENT (OUT) :: ex,vx ! exchange energy per electron (ex) and potential (vx)

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

      ! .. local variables ..
      REAL :: um,uk,ul,exunif,fs,fss,fxpbe,p0,s2
      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
      REAL, PARAMETER :: thrd=1.e0/3.e0
      REAL, PARAMETER :: thrd4=4.e0/3.e0
      REAL, PARAMETER :: ax=-0.738558766382022405884230032680836e0  ! -0.75*(3/pi)^(1/3)

c----------------------------------------------------------------------
c uk, ul defined after [a](13)  (for icorr==7)
c----------------------------------------------------------------------
      IF ((icorr == 7).OR.(icorr == 10).OR.(icorr == 11)) THEN
        uk=0.8040
      ELSEIF (icorr.EQ.8) THEN
        uk=1.2450
      ELSEIF (icorr.EQ.9) THEN    ! changed to [c]
        uk=0.8040
      ELSEIF (icorr.EQ.icorr_pbe0   .OR. icorr.EQ.icorr_hse
     +   .OR. icorr.EQ.icorr_hseloc .OR. icorr.EQ.icorr_vhse ) THEN
        uk=0.8040
      ELSE
        WRITE (6,'(//'' icorr is not correctly transferred. icorr='',i5)
     &    ') icorr
        STOP
      ENDIF
      IF (icorr == 11) THEN  ! pbe_sol
        um=0.123456790123456d0
      ELSE
        um=0.2195149727645171e0
      ENDIF

      ul=um/uk
c----------------------------------------------------------------------
c construct lda exchange energy density: 
                          ! e_x[unif] = -0.75 * (3/pi)^(1/3) * rho^(4/3)  
c----------------------------------------------------------------------
      exunif = ax*rho**thrd
      IF (lgga.EQ.0) THEN
          ex = exunif
          vx = ex*thrd4
          RETURN
      ENDIF
c----------------------------------------------------------------------
c construct pbe enhancement factor
c       e_x[pbe]=e_x[unif]*fxpbe(s)
c       fxpbe(s)=1+uk-uk/(1+ul*s*s)                 [a](13)
c----------------------------------------------------------------------
      s2 = s*s
c+gu
      IF (icorr.EQ.7 .OR. icorr.EQ. 8 .OR. icorr.EQ.icorr_pbe0 .OR.
     &    icorr.EQ.icorr_hse .OR. icorr.EQ.icorr_hseloc .OR.
     &    icorr.EQ.icorr_vhse ) THEN
        ! calculate fxpbe
        p0 = 1.e0 + ul*s2
        fxpbe = 1e0 + uk - uk/p0
      ELSEIF (icorr.EQ.9) THEN
        p0 = exp( - ul*s2 )
        fxpbe = 1e0 + uk * ( 1e0 - p0 )
      ELSEIF (icorr.EQ.10) THEN
        css = 1+cwu*s2*s2
        xwu = teo*s2 + (um-teo)*s2*exp(-s2) + log(css)
        p0 = 1.e0 + xwu/uk
        fxpbe = 1e0 + uk - uk/p0
      ENDIF
c-gu
      ! Mixing of short and long range components
      IF ( icorr == icorr_hse .OR. icorr == icorr_hseloc .OR.
     +     icorr == icorr_vhse ) THEN
        ! 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
#ifndef CPP_GF
        CALL calculateEnhancementFactor(kF, s, fxhse, dFx_ds, d2Fx_ds2,
     &                                  dFx_dkF, d2Fx_dsdkF)
        IF ( icorr == icorr_vhse ) THEN
          ex = exunif * (fxpbe - aMix_VHSE() * fxhse )
        ELSE
          ex = exunif * (fxpbe - aMix_HSE * fxhse )
        END IF
#endif
      ELSE
        ex = exunif*fxpbe
      END IF
      IF (lpot.EQ.0) RETURN
c----------------------------------------------------------------------
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
c----------------------------------------------------------------------
c+gu
      IF (icorr.EQ.7 .OR. icorr.EQ.8 .OR. icorr.EQ.icorr_pbe0 .OR.
     &    icorr.EQ.icorr_hse .OR. icorr.EQ.icorr_hseloc .OR.
     &    icorr.EQ.icorr_vhse ) THEN
        ! derivatives for the pbe part
        fs = 2.e0*uk*ul/ (p0*p0)
        fss = -4.e0*ul*s*fs/p0
      ELSEIF (icorr.EQ.9) THEN
        fs = 2.e0*ul*p0
        fss = -2.e0*ul*s*fs
      ELSEIF (icorr.EQ.10) THEN
        dxwu = 2*teo + 2*(um-teo)*exp(-s2)*(1-s2) + 4*cwu*s2/css
        fs = dxwu / (p0*p0)
        ddx = 4*s*((um-teo)*exp(-s2)*(s2-2)+2*cwu*(1-cwu*s2*s2)/css**2)
        fss = ( ddx - 2*s*dxwu*dxwu/(p0*uk) ) / (p0*p0)
      ENDIF
c-gu
c----------------------------------------------------------------------
c calculate potential from [b](24)
c----------------------------------------------------------------------
      vx = exunif* (thrd4*fxpbe- (u-thrd4*s2*s)*fss-v*fs)

      IF ( icorr /= icorr_hse .AND. icorr /= icorr_hseloc .AND.
     &     icorr /= icorr_vhse ) RETURN

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

      END SUBROUTINE exchpbe
      END MODULE m_exchpbe