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

          Line data    Source code
       1             : MODULE m_corpbe
       2             : 
       3             : !----------------------------------------------------------------------
       4             : !  official pbe correlation code. k. burke, may 14, 1996.
       5             : !----------------------------------------------------------------------
       6             : ! references:
       7             : ! [a] j.p.~perdew, k.~burke, and m.~ernzerhof,
       8             : !     {\sl generalized gradient approximation made simple}, sub.
       9             : !     to phys. rev.lett. may 1996.
      10             : ! [b] j. p. perdew, k. burke, and y. wang, {\sl real-space cutoff
      11             : !     construction of a generalized gradient approximation:  the pw91
      12             : !     density functional}, submitted to phys. rev. b, feb. 1996.
      13             : ! [c] j. p. perdew and y. wang, phys. rev. b {\bf 45}, 13244 (1992).
      14             : !----------------------------------------------------------------------
      15             : CONTAINS
      16    61863127 :    SUBROUTINE corpbe( &
      17             :       l_pbes,rs,zet,t,uu,vv,ww,lgga,lpot, &
      18             :       ec,vcup,vcdn,h,dvcup,dvcdn)
      19             : !----------------------------------------------------------------------
      20             : !  input: rs=seitz radius=(3/4pi rho)^(1/3)
      21             : !       : zet=relative spin polarization = (rhoup-rhodn)/rho
      22             : !       : t=abs(grad rho)/(rho*2.*ks*g)  -- only needed for pbe
      23             : !       : uu=(grad rho)*grad(abs(grad rho))/(rho**2 * (2*ks*g)**3)
      24             : !       : vv=(laplacian rho)/(rho * (2*ks*g)**2)
      25             : !       : ww=(grad rho)*(grad zet)/(rho * (2*ks*g)**2
      26             : !       :  uu,vv,ww, only needed for pbe potential
      27             : !       : lgga=flag to do gga (0=>lsd only)
      28             : !       : lpot=flag to do potential (0=>energy only)
      29             : !  output: ec=lsd correlation energy from [a]
      30             : !        : vcup=lsd up correlation potential
      31             : !        : vcdn=lsd dn correlation potential
      32             : !        : h=nonlocal part of correlation energy per electron
      33             : !        : dvcup=nonlocal correction to vcup
      34             : !        : dvcdn=nonlocal correction to vcdn
      35             : !----------------------------------------------------------------------
      36             : 
      37             :       USE m_pbecor2
      38             :       IMPLICIT NONE
      39             :       LOGICAL,INTENT(IN)    :: l_pbes
      40             :       INTEGER, INTENT (IN)  :: lgga,lpot
      41             :       REAL,    INTENT (IN)  :: rs,zet,t,uu,vv,ww
      42             :       REAL,    INTENT (OUT) :: dvcdn,dvcup,ec,h,vcdn,vcup
      43             : 
      44             : ! thrd*=various multiples of 1/3
      45             : ! numbers for use in lsd energy spin-interpolation formula, [c](9).
      46             : !      gam= 2^(4/3)-2
      47             : !      fzz=f''(0)= 8/(9*gam)
      48             : ! numbers for construction of pbe
      49             : !      gamma=(1-log(2))/pi^2
      50             : !      bet=coefficient in gradient expansion for correlation, [a](4).
      51             : !      eta=small number to stop d phi/ dzeta from blowing up at
      52             : !          |zeta|=1.
      53             : 
      54             :       REAL, PARAMETER :: thrd=1.e0/3.e0
      55             :       REAL, PARAMETER :: thrdm=-thrd
      56             :       REAL, PARAMETER :: thrd2=2.e0*thrd
      57             :       REAL, PARAMETER :: sixthm=thrdm/2.e0
      58             :       REAL, PARAMETER :: thrd4=4.e0*thrd
      59             :       REAL, PARAMETER :: gam=0.5198420997897463295344212145565e0
      60             :       REAL, PARAMETER :: fzz=8.e0/ (9.e0*gam)
      61             :       REAL, PARAMETER :: gamma=0.03109069086965489503494086371273e0
      62             :       REAL, PARAMETER :: eta=1.e-12
      63             : !----------------------------------------------------------------------
      64             : !----------------------------------------------------------------------
      65             : ! find lsd energy contributions, using [c](10) and table i[c].
      66             : ! eu=unpolarized lsd correlation energy
      67             : ! eurs=deu/drs
      68             : ! ep=fully polarized lsd correlation energy
      69             : ! eprs=dep/drs
      70             : ! alfm=-spin stiffness, [c](3).
      71             : ! alfrsm=-dalpha/drs
      72             : ! f=spin-scaling factor from [c](9).
      73             : ! construct ec, using [c](8)
      74             : 
      75             :       REAL :: alfm,alfrsm,b,b2,bec,bg,comm,ecrs,eczet,ep,eprs,eu,eurs, &
      76             :               f,fac,fact0,fact1,fact2,fact3,fact5,fz,g,g3,g4,gz,hb,hbt,hrs, &
      77             :               hrst,ht,htt,hz,hzt,pon,pref,q4,q5,q8,q9,rsthrd,rtrs, &
      78             :               t2,t4,t6,z4,delt,bet
      79             : !     ..
      80             : 
      81    61863127 :       IF (l_pbes) THEN ! PBE_sol
      82             :          bet=0.046e0
      83             :       ELSE
      84    61863127 :          bet=0.06672455060314922e0
      85             :       ENDIF
      86    61863127 :       delt=bet/gamma
      87             : 
      88    61863127 :       rtrs = sqrt(rs)
      89             :       CALL pbecor2(0.0310907,0.21370,7.5957,3.5876,1.6382, &
      90    61863127 :       &            0.49294,rtrs,eu,eurs)
      91             :       CALL pbecor2(0.01554535,0.20548,14.1189,6.1977,3.3662, &
      92    61863127 :       &            0.62517,rtrs,ep,eprs)
      93             :       CALL pbecor2(0.0168869,0.11125,10.357,3.6231,0.88026, &
      94    61863127 :       &            0.49671,rtrs,alfm,alfrsm)
      95    61863127 :       z4 = zet**4
      96    61863127 :       f = ((1.e0+zet)**thrd4+ (1.e0-zet)**thrd4-2.e0)/gam
      97    61863127 :       ec = eu* (1.e0-f*z4) + ep*f*z4 - alfm*f* (1.e0-z4)/fzz
      98             : !----------------------------------------------------------------------
      99             : !----------------------------------------------------------------------
     100             : ! lsd potential from [c](a1)
     101             : ! ecrs = dec/drs [c](a2)
     102             : ! eczet=dec/dzeta [c](a3)
     103             : ! fz = df/dzeta [c](a4)
     104    61863127 :       ecrs = eurs* (1.e0-f*z4) + eprs*f*z4 - alfrsm*f* (1.e0-z4)/fzz
     105    61863127 :       fz = thrd4* ((1.e0+zet)**thrd- (1.e0-zet)**thrd)/gam
     106             :       eczet = 4.e0* (zet**3)*f* (ep-eu+alfm/fzz) + &
     107    61863127 :               fz* (z4*ep-z4*eu- (1.e0-z4)*alfm/fzz)
     108    61863127 :       comm = ec - rs*ecrs/3.e0 - zet*eczet
     109    61863127 :       vcup = comm + eczet
     110    61863127 :       vcdn = comm - eczet
     111    61863127 :       IF (lgga == 0) RETURN
     112             : !----------------------------------------------------------------------
     113             : !----------------------------------------------------------------------
     114             : ! pbe correlation energy
     115             : ! g=phi(zeta), given after [a](3)
     116             : ! delt=bet/gamma
     117             : ! b=a of [a](8)
     118    61863127 :       g = ((1.e0+zet)**thrd2+ (1.e0-zet)**thrd2)/2.e0
     119    61863127 :       g3 = g**3
     120    61863127 :       pon = -ec/ (g3*gamma)
     121    61863127 :       b = delt/ (exp(pon)-1.e0)
     122    61863127 :       b2 = b*b
     123    61863127 :       t2 = t*t
     124    61863127 :       t4 = t2*t2
     125    61863127 :       q4 = 1.e0 + b*t2
     126    61863127 :       q5 = 1.e0 + b*t2 + b2*t4
     127    61863127 :       h = g3* (bet/delt)*log(1.e0+delt*q4*t2/q5)
     128    61863127 :       IF (lpot == 0) RETURN
     129             : !----------------------------------------------------------------------
     130             : !----------------------------------------------------------------------
     131             : ! energy done. now the potential, using appendix e of [b].
     132    61863127 :       g4 = g3*g
     133    61863127 :       t6 = t4*t2
     134    61863127 :       rsthrd = rs/3.e0
     135             :       gz = (((1.e0+zet)**2+eta)**sixthm- ((1.e0-zet)**2+eta)**sixthm)/ &
     136    61863127 :       &      3.e0
     137    61863127 :       fac = delt/b + 1.e0
     138    61863127 :       bg = -3.e0*b2*ec*fac/ (bet*g4)
     139    61863127 :       bec = b2*fac/ (bet*g3)
     140    61863127 :       q8 = q5*q5 + delt*q4*q5*t2
     141    61863127 :       q9 = 1.e0 + 2.e0*b*t2
     142    61863127 :       hb = -bet*g3*b*t6* (2.e0+b*t2)/q8
     143    61863127 :       hrs = -rsthrd*hb*bec*ecrs
     144    61863127 :       fact0 = 2.e0*delt - 6.e0*b
     145    61863127 :       fact1 = q5*q9 + q4*q9*q9
     146    61863127 :       hbt = 2.e0*bet*g3*t4* ((q4*q5*fact0-delt*fact1)/q8)/q8
     147    61863127 :       hrst = rsthrd*t2*hbt*bec*ecrs
     148    61863127 :       hz = 3.e0*gz*h/g + hb* (bg*gz+bec*eczet)
     149    61863127 :       ht = 2.e0*bet*g3*q9/q8
     150    61863127 :       hzt = 3.e0*gz*ht/g + hbt* (bg*gz+bec*eczet)
     151    61863127 :       fact2 = q4*q5 + b*t2* (q4*q9+q5)
     152    61863127 :       fact3 = 2.e0*b*q5*q9 + delt*fact2
     153    61863127 :       htt = 4.e0*bet*g3*t* (2.e0*b/q8- (q9*fact3/q8)/q8)
     154    61863127 :       comm = h + hrs + hrst + t2*ht/6.e0 + 7.e0*t2*t*htt/6.e0
     155    61863127 :       pref = hz - gz*t2*ht/g
     156    61863127 :       fact5 = gz* (2.e0*ht+t*htt)/g
     157    61863127 :       comm = comm - pref*zet - uu*htt - vv*ht - ww* (hzt-fact5)
     158    61863127 :       dvcup = comm + pref
     159    61863127 :       dvcdn = comm - pref
     160             : 
     161             :    END SUBROUTINE corpbe
     162             : END MODULE m_corpbe

Generated by: LCOV version 1.13