LCOV - code coverage report
Current view: top level - xc-pot - exchpbe.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 29 39 74.4 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          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   247452508 :    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   247452508 :       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   247452508 :       exunif = ax*rho**thrd
      77   247452508 :       IF (lgga == 0) THEN
      78   123726254 :          ex = exunif
      79   123726254 :          vx = ex*thrd4
      80   123726254 :          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   123726254 :       s2 = s*s
      88             : 
      89             : !     calculate fxpbe
      90   123726254 :       p0 = 1.e0 + ul*s2
      91   123726254 :       fxpbe = 1e0 + xcpot%uk - xcpot%uk/p0
      92   123726254 :       IF (xcpot%is_Rpbe) THEN
      93           0 :          p0 = exp( - ul*s2 )
      94           0 :          fxpbe = 1e0 + xcpot%uk * ( 1e0 - p0 )
      95   123726254 :       ELSEIF (xcpot%is_wc) THEN
      96    11565304 :          css = 1+cwu*s2*s2
      97    11565304 :          xwu = teo*s2 + (xcpot%um-teo)*s2*exp(-s2) + log(css)
      98    11565304 :          p0 = 1.e0 + xwu/xcpot%uk
      99    11565304 :          fxpbe = 1e0 + xcpot%uk - xcpot%uk/p0
     100             :       ENDIF
     101             : !     -gu
     102             : !     Mixing of short and long range components
     103   123726254 :       IF (xcpot%is_hse) THEN
     104             :          !     ex = (1-a)ex,SR + ex,LR
     105             :          !     = (1-a)ex,SR + (ex,PBE - ex,SR)
     106             :          !     = (fxpbe - a*fxhse)*exunif
     107             :          !     Calculate the enhancement factor fxhse and its derivatives
     108             :          !     as integral over the exchange hole (cf. [d])
     109           0 :          kF = (3.0 * pi_const**2 * rho)**thrd
     110           0 :          CALL judft_error("HSE not implemented",calledby="exchpbe")
     111             :          ! his creates a depency loop
     112             :          !     CALL calculateEnhancementFactor(kF, s, fxhse, dFx_ds, d2Fx_ds2,
     113             :          !     &        dFx_dkF, d2Fx_dsdkF)
     114           0 :          ex = exunif * (fxpbe - xcpot%exchange_weight * fxhse )
     115             :       ELSE
     116   123726254 :          ex = exunif*fxpbe
     117             :       END IF
     118   123726254 :       IF (lpot == 0) RETURN
     119             : !----------------------------------------------------------------------
     120             : !     energy done. now the potential:
     121             : !     find first and second derivatives of fx w.r.t s.
     122             : !     fs=(1/s)*d fxpbe/ ds
     123             : !     fss=d fs/ds
     124             : !----------------------------------------------------------------------
     125             : !     derivatives for the pbe part
     126   123726254 :       fs = 2.e0*xcpot%uk*ul/ (p0*p0)
     127   123726254 :       fss = -4.e0*ul*s*fs/p0
     128   123726254 :       IF (xcpot%is_Rpbe) THEN
     129           0 :          fs = 2.e0*ul*p0
     130           0 :          fss = -2.e0*ul*s*fs
     131   123726254 :       ELSEIF (xcpot%is_wc) THEN
     132    11565304 :          dxwu = 2*teo + 2*(xcpot%um-teo)*exp(-s2)*(1-s2) + 4*cwu*s2/css
     133    11565304 :          fs = dxwu / (p0*p0)
     134             :          ddx = 4*s*((xcpot%um-teo)*exp(-s2)*(s2-2)+2*cwu* &
     135    11565304 :                     (1-cwu*s2*s2)/css**2)
     136    11565304 :          fss = ( ddx - 2*s*dxwu*dxwu/(p0*xcpot%uk) ) / (p0*p0)
     137             :       ENDIF
     138             : !     -gu
     139             : !----------------------------------------------------------------------
     140             : !     calculate potential from [b](24)
     141             : !----------------------------------------------------------------------
     142   123726254 :       vx = exunif* (thrd4*fxpbe- (u-thrd4*s2*s)*fss-v*fs)
     143             : 
     144   123726254 :       IF ( .NOT. (xcpot%is_hse)) RETURN
     145             : 
     146             : !----------------------------------------------------------------------
     147             : !     short ranged potential (HSE functional)
     148             : !----------------------------------------------------------------------
     149             : !     calculate fs and fss for the HSE functional
     150             : !     where the 1st and 2nd derivative of Fx are known
     151           0 :       fs    = dFx_ds / s
     152           0 :       fss   = (d2Fx_ds2 - fs) / s
     153             :       vx_sr = exunif * ( thrd4*fxhse - (u-thrd4*s2*s)*fss - v*fs &
     154           0 :                         + thrd*kF * ( dFx_dkF - d2Fx_dsdkF*s ) )
     155             : 
     156             :    END SUBROUTINE exchpbe
     157             : END MODULE m_exchpbe

Generated by: LCOV version 1.13