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

          Line data    Source code
       1             : MODULE m_easypbe
       2             : !----------------------------------------------------------------------
       3             : !     easypbe---exchpbe
       4             : !             --corpbe  --- pbecor2
       5             : !----------------------------------------------------------------------
       6             : ! easypbe is a driver for the pbe subroutines, using simple inputs
       7             : ! k. burke, may 14, 1996.
       8             : !----------------------------------------------------------------------
       9             : CONTAINS
      10    61863127 :    SUBROUTINE easypbe(xcpot, &
      11             :                       up,agrup,delgrup,uplap, &
      12             :                       dn,agrdn,delgrdn,dnlap, &
      13             :                       agr,delgr,lcor,lpot, &
      14             :                       exlsd,vxuplsd,vxdnlsd, &
      15             :                       eclsd,vcuplsd,vcdnlsd, &
      16             :                       expbe,vxuppbe,vxdnpbe, &
      17             :                       ecpbe,vcuppbe,vcdnpbe, &
      18             :                       vxupsr,vxdnsr)
      19             : 
      20             :       USE m_exchpbe
      21             :       USE m_corpbe
      22             :       USE m_types_xcpot_data
      23             :       IMPLICIT NONE
      24             : 
      25             : ! .. Arguments ..
      26             :       TYPE(t_xcpot_data),INTENT(IN)::xcpot
      27             :       INTEGER, INTENT (IN) :: lcor ! flag to do correlation(=0=>don't)
      28             :       INTEGER, INTENT (IN) :: lpot     ! flag to do potential  (=0=>don't)
      29             :       REAL,    INTENT (IN) :: up,dn    ! density (spin up & down)
      30             :       REAL,    INTENT (IN) :: agrup,agrdn     ! |grad up|, |grad down|
      31             :       REAL,    INTENT (IN) :: delgrup,delgrdn ! delgrup=(grad up).(grad |grad up|)  (in pw91: gggru)
      32             :       REAL,    INTENT (IN) :: uplap,dnlap     ! grad^2 up=laplacian of up           (in pw91: g2ru)
      33             :       REAL,    INTENT (IN) :: agr,delgr       ! agr=|grad rho|,
      34             : ! delgr=(grad rho).(grad |grad rho|)  (in pw91: gggr)
      35             :       REAL,    INTENT (OUT) :: vxuplsd,vxdnlsd ! up/down lsd exchange potential
      36             :       REAL,    INTENT (OUT) :: vcuplsd,vcdnlsd ! up/down lsd correlation potential
      37             :       REAL,    INTENT (OUT) :: exlsd,eclsd     ! lsd exchange / correlation energy density
      38             :       REAL,    INTENT (OUT) :: vxuppbe,vxdnpbe ! as above, but pbe quantities
      39             :       REAL,    INTENT (OUT) :: vcuppbe,vcdnpbe
      40             :       REAL,    INTENT (OUT) :: expbe,ecpbe     ! note that : exlsd=int d^3r rho(r) exlsd(r)
      41             :       REAL,    INTENT (OUT) :: vxupsr,vxdnsr
      42             : 
      43             : ! .. local variables ..
      44             : ! for exchange:
      45             :       REAL :: fk ! local fermi wavevector for 2*up=(3 pi^2 (2up))^(1/3)
      46             :       REAL :: s  ! dimensionless density gradient=|grad rho|/ (2*fk*rho)_(rho=2*up)
      47             :       REAL :: u  ! delgrad/(rho^2*(2*fk)**3)_(rho=2*up)
      48             :       REAL :: v  ! laplacian/(rho*(2*fk)**2)_(rho=2*up)
      49             : ! for correlation:
      50             :       REAL :: zet    ! (up-dn)/rho
      51             :       REAL :: g      ! phi(zeta)
      52             :       REAL :: rs     ! (3/(4pi*rho))^(1/3)=local seitz radius=alpha/fk
      53             :       REAL :: sk     ! ks=thomas-fermi screening wavevector=sqrt(4fk/pi)
      54             :       REAL :: twoksg ! 2*ks*phi
      55             :       REAL :: t      ! correlation dimensionless gradient=|grad rho|/(2*ks*phi*rho)
      56             :       REAL :: uu     ! delgrad/(rho^2*twoksg^3)
      57             :       REAL :: rholap ! laplacian
      58             :       REAL :: vv     ! laplacian/(rho*twoksg^2)
      59             :       REAL :: ww     ! (|grad up|^2-|grad dn|^2-zet*|grad rho|^2)/(rho*twoksg)^2
      60             :       REAL :: ec     ! lsd correlation energy
      61             :       REAL :: vcup   ! lsd up correlation potential
      62             :       REAL :: vcdn   ! lsd down correlation potential
      63             :       REAL :: h      ! gradient correction to correlation energy
      64             :       REAL :: dvcup  ! gradient correction to up correlation potential
      65             :       REAL :: dvcdn  ! gradient correction to down correlation potential
      66             : 
      67             :       REAL :: rdum,eta,exdnlsd,exdnpbe,exuplsd,exuppbe,rho,rho2
      68             : !     ..
      69             :       REAL, PARAMETER :: thrd = 1.e0/3.e0
      70             :       REAL, PARAMETER :: thrd2 = 2.e0*thrd
      71             :       REAL, PARAMETER :: pi = 3.1415926535897932384626433832795e0
      72             :       REAL, PARAMETER :: pi32 = 29.608813203268075856503472999628e0 ! 3 pi**2
      73             :       REAL, PARAMETER :: alpha=1.91915829267751300662482032624669e0 ! (9pi/4)**thrd
      74             : 
      75    61863127 :       rho2 = 2.e0*up
      76             : 
      77             : !----------------------------------------------------------------------
      78             : ! pbe exchange
      79             : ! use  ex[up,dn]=0.5*(ex[2*up]+ex[2*dn]) (i.e., exact spin-scaling)
      80             : ! do up exchange
      81             : !----------------------------------------------------------------------
      82             : 
      83    61863127 :       IF (rho2 > 1e-18) THEN
      84    61863127 :          fk = (pi32*rho2)**thrd
      85    61863127 :          s = 2.e0*agrup/ (2.e0*fk*rho2)
      86    61863127 :          u = 4.e0*delgrup/ (rho2*rho2* (2.e0*fk)**3)
      87    61863127 :          v = 2.e0*uplap/ (rho2* (2.e0*fk)**2)
      88             : 
      89    61863127 :          CALL exchpbe(xcpot,rho2,s,u,v,0,lpot,exuplsd,vxuplsd,rdum)
      90    61863127 :          CALL exchpbe(xcpot,rho2,s,u,v,1,lpot,exuppbe,vxuppbe,vxupsr)
      91             : 
      92             :       ELSE
      93             : 
      94           0 :          exuplsd = 0.e0
      95           0 :          vxuplsd = 0.e0
      96           0 :          exuppbe = 0.e0
      97           0 :          vxuppbe = 0.e0
      98           0 :          vxupsr  = 0.e0
      99             : 
     100             :       ENDIF
     101             : 
     102             : ! repeat for down
     103    61863127 :       rho2 = 2.e0*dn
     104             : 
     105    61863127 :       IF (rho2 > 1e-18) THEN
     106             : 
     107    61863127 :          fk = (pi32*rho2)**thrd
     108    61863127 :          s = 2.e0*agrdn/ (2.e0*fk*rho2)
     109    61863127 :          u = 4.e0*delgrdn/ (rho2*rho2* (2.e0*fk)**3)
     110    61863127 :          v = 2.e0*dnlap/ (rho2* (2.e0*fk)**2)
     111             : 
     112    61863127 :          CALL exchpbe(xcpot,rho2,s,u,v,0,lpot,exdnlsd,vxdnlsd,rdum)
     113    61863127 :          CALL exchpbe(xcpot,rho2,s,u,v,1,lpot,exdnpbe,vxdnpbe,vxdnsr)
     114             : 
     115             :       ELSE
     116             : 
     117           0 :          exdnlsd = 0.e0
     118           0 :          vxdnlsd = 0.e0
     119           0 :          exdnpbe = 0.e0
     120           0 :          vxdnpbe = 0.e0
     121           0 :          vxdnsr  = 0.e0
     122             : 
     123             :       ENDIF
     124             : 
     125             : ! construct total density and contribution to ex
     126    61863127 :       rho = up + dn
     127    61863127 :       exlsd = (exuplsd*up+exdnlsd*dn)/rho
     128    61863127 :       expbe = (exuppbe*up+exdnpbe*dn)/rho
     129             : 
     130    61863127 :       IF (lcor == 0) RETURN
     131             : !----------------------------------------------------------------------
     132             : ! now do correlation
     133             : !----------------------------------------------------------------------
     134             : 
     135    61863127 :       IF (rho < 1.e-18) RETURN
     136             : 
     137    61863127 :       zet = (up-dn)/rho
     138             : 
     139             : ! 999+
     140             : !     eta: eta should not be smaller than 1.e-16.
     141             : !          otherwise will cause floating invalid.
     142             : !          if bigger, the last digit may differ
     143             : !          from the run by aix without this zet-guard.
     144    61863127 :       eta = 1.e-16
     145    61863127 :       zet = min(zet,1.0-eta)
     146    61863127 :       zet = max(zet,-1.0+eta)
     147             : ! 999-
     148             : 
     149    61863127 :       g = ((1.e0+zet)**thrd2+ (1.e0-zet)**thrd2)/2.e0
     150    61863127 :       fk = (pi32*rho)**thrd
     151    61863127 :       rs = alpha/fk
     152    61863127 :       sk = sqrt(4.e0*fk/pi)
     153    61863127 :       twoksg = 2.e0*sk*g
     154    61863127 :       t = agr/ (twoksg*rho)
     155    61863127 :       uu = delgr/ (rho*rho*twoksg**3)
     156    61863127 :       rholap = uplap + dnlap
     157    61863127 :       vv = rholap/ (rho*twoksg**2)
     158    61863127 :       ww = (agrup**2-agrdn**2-zet*agr**2)/ (rho*rho*twoksg**2)
     159             : 
     160             :       CALL corpbe( &
     161             :          xcpot%is_PBEs,rs,zet,t,uu,vv,ww,1,lpot, &
     162    61863127 :          ec,vcup,vcdn,h,dvcup,dvcdn)
     163             : 
     164    61863127 :       eclsd = ec
     165    61863127 :       ecpbe = ec + h
     166    61863127 :       vcuplsd = vcup
     167    61863127 :       vcdnlsd = vcdn
     168    61863127 :       vcuppbe = vcup + dvcup
     169    61863127 :       vcdnpbe = vcdn + dvcdn
     170             : 
     171             :    END SUBROUTINE easypbe
     172             : END MODULE m_easypbe

Generated by: LCOV version 1.13