easypbe.f90 6.15 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
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
MODULE m_easypbe
!----------------------------------------------------------------------
!     easypbe---exchpbe
!             --corpbe  --- pbecor2
!----------------------------------------------------------------------
! easypbe is a driver for the pbe subroutines, using simple inputs
! k. burke, may 14, 1996.
!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE easypbe(xcpot, &
                      up,agrup,delgrup,uplap, &
                      dn,agrdn,delgrdn,dnlap, &
                      agr,delgr,lcor,lpot, &
                      exlsd,vxuplsd,vxdnlsd, &
                      eclsd,vcuplsd,vcdnlsd, &
                      expbe,vxuppbe,vxdnpbe, &
                      ecpbe,vcuppbe,vcdnpbe, &
                      vxupsr,vxdnsr)

      USE m_exchpbe
      USE m_corpbe
      USE m_types_xcpot_data
      IMPLICIT NONE

! .. Arguments ..
      TYPE(t_xcpot_data),INTENT(IN)::xcpot
      INTEGER, INTENT (IN) :: lcor ! flag to do correlation(=0=>don't)
      INTEGER, INTENT (IN) :: lpot     ! flag to do potential  (=0=>don't)
      REAL,    INTENT (IN) :: up,dn    ! density (spin up & down)
      REAL,    INTENT (IN) :: agrup,agrdn     ! |grad up|, |grad down|
      REAL,    INTENT (IN) :: delgrup,delgrdn ! delgrup=(grad up).(grad |grad up|)  (in pw91: gggru)
      REAL,    INTENT (IN) :: uplap,dnlap     ! grad^2 up=laplacian of up           (in pw91: g2ru)
      REAL,    INTENT (IN) :: agr,delgr       ! agr=|grad rho|,
! delgr=(grad rho).(grad |grad rho|)  (in pw91: gggr)
      REAL,    INTENT (OUT) :: vxuplsd,vxdnlsd ! up/down lsd exchange potential
      REAL,    INTENT (OUT) :: vcuplsd,vcdnlsd ! up/down lsd correlation potential
      REAL,    INTENT (OUT) :: exlsd,eclsd     ! lsd exchange / correlation energy density
      REAL,    INTENT (OUT) :: vxuppbe,vxdnpbe ! as above, but pbe quantities
      REAL,    INTENT (OUT) :: vcuppbe,vcdnpbe
      REAL,    INTENT (OUT) :: expbe,ecpbe     ! note that : exlsd=int d^3r rho(r) exlsd(r)
      REAL,    INTENT (OUT) :: vxupsr,vxdnsr

! .. local variables ..
! for exchange:
      REAL :: fk ! local fermi wavevector for 2*up=(3 pi^2 (2up))^(1/3)
      REAL :: s  ! dimensionless density gradient=|grad rho|/ (2*fk*rho)_(rho=2*up)
      REAL :: u  ! delgrad/(rho^2*(2*fk)**3)_(rho=2*up)
      REAL :: v  ! laplacian/(rho*(2*fk)**2)_(rho=2*up)
! for correlation:
      REAL :: zet    ! (up-dn)/rho
      REAL :: g      ! phi(zeta)
      REAL :: rs     ! (3/(4pi*rho))^(1/3)=local seitz radius=alpha/fk
      REAL :: sk     ! ks=thomas-fermi screening wavevector=sqrt(4fk/pi)
      REAL :: twoksg ! 2*ks*phi
      REAL :: t      ! correlation dimensionless gradient=|grad rho|/(2*ks*phi*rho)
      REAL :: uu     ! delgrad/(rho^2*twoksg^3)
      REAL :: rholap ! laplacian
      REAL :: vv     ! laplacian/(rho*twoksg^2)
      REAL :: ww     ! (|grad up|^2-|grad dn|^2-zet*|grad rho|^2)/(rho*twoksg)^2
      REAL :: ec     ! lsd correlation energy
      REAL :: vcup   ! lsd up correlation potential
      REAL :: vcdn   ! lsd down correlation potential
      REAL :: h      ! gradient correction to correlation energy
      REAL :: dvcup  ! gradient correction to up correlation potential
      REAL :: dvcdn  ! gradient correction to down correlation potential

      REAL :: rdum,eta,exdnlsd,exdnpbe,exuplsd,exuppbe,rho,rho2
!     ..
      REAL, PARAMETER :: thrd = 1.e0/3.e0
      REAL, PARAMETER :: thrd2 = 2.e0*thrd
      REAL, PARAMETER :: pi = 3.1415926535897932384626433832795e0
      REAL, PARAMETER :: pi32 = 29.608813203268075856503472999628e0 ! 3 pi**2
      REAL, PARAMETER :: alpha=1.91915829267751300662482032624669e0 ! (9pi/4)**thrd

      rho2 = 2.e0*up

!----------------------------------------------------------------------
! pbe exchange
! use  ex[up,dn]=0.5*(ex[2*up]+ex[2*dn]) (i.e., exact spin-scaling)
! do up exchange
!----------------------------------------------------------------------

      IF (rho2 > 1e-18) THEN
         fk = (pi32*rho2)**thrd
         s = 2.e0*agrup/ (2.e0*fk*rho2)
         u = 4.e0*delgrup/ (rho2*rho2* (2.e0*fk)**3)
         v = 2.e0*uplap/ (rho2* (2.e0*fk)**2)

         CALL exchpbe(xcpot,rho2,s,u,v,0,lpot,exuplsd,vxuplsd,rdum)
         CALL exchpbe(xcpot,rho2,s,u,v,1,lpot,exuppbe,vxuppbe,vxupsr)

      ELSE

         exuplsd = 0.e0
         vxuplsd = 0.e0
         exuppbe = 0.e0
         vxuppbe = 0.e0
         vxupsr  = 0.e0

      ENDIF

! repeat for down
      rho2 = 2.e0*dn

      IF (rho2 > 1e-18) THEN

         fk = (pi32*rho2)**thrd
         s = 2.e0*agrdn/ (2.e0*fk*rho2)
         u = 4.e0*delgrdn/ (rho2*rho2* (2.e0*fk)**3)
         v = 2.e0*dnlap/ (rho2* (2.e0*fk)**2)

         CALL exchpbe(xcpot,rho2,s,u,v,0,lpot,exdnlsd,vxdnlsd,rdum)
         CALL exchpbe(xcpot,rho2,s,u,v,1,lpot,exdnpbe,vxdnpbe,vxdnsr)

      ELSE

         exdnlsd = 0.e0
         vxdnlsd = 0.e0
         exdnpbe = 0.e0
         vxdnpbe = 0.e0
         vxdnsr  = 0.e0

      ENDIF

! construct total density and contribution to ex
      rho = up + dn
      exlsd = (exuplsd*up+exdnlsd*dn)/rho
      expbe = (exuppbe*up+exdnpbe*dn)/rho

      IF (lcor == 0) RETURN
!----------------------------------------------------------------------
! now do correlation
!----------------------------------------------------------------------

      IF (rho < 1.e-18) RETURN

      zet = (up-dn)/rho

! 999+
!     eta: eta should not be smaller than 1.e-16.
!          otherwise will cause floating invalid.
!          if bigger, the last digit may differ
!          from the run by aix without this zet-guard.
      eta = 1.e-16
      zet = min(zet,1.0-eta)
      zet = max(zet,-1.0+eta)
! 999-

      g = ((1.e0+zet)**thrd2+ (1.e0-zet)**thrd2)/2.e0
      fk = (pi32*rho)**thrd
      rs = alpha/fk
      sk = sqrt(4.e0*fk/pi)
      twoksg = 2.e0*sk*g
      t = agr/ (twoksg*rho)
      uu = delgr/ (rho*rho*twoksg**3)
      rholap = uplap + dnlap
      vv = rholap/ (rho*twoksg**2)
      ww = (agrup**2-agrdn**2-zet*agr**2)/ (rho*rho*twoksg**2)

      CALL corpbe( &
         xcpot%is_PBEs,rs,zet,t,uu,vv,ww,1,lpot, &
         ec,vcup,vcdn,h,dvcup,dvcdn)

      eclsd = ec
      ecpbe = ec + h
      vcuplsd = vcup
      vcdnlsd = vcdn
      vcuppbe = vcup + dvcup
      vcdnpbe = vcdn + dvcdn

   END SUBROUTINE easypbe
END MODULE m_easypbe