Line data    Source code

       1              : MODULE m_xch91
       2              : !.....-----------------------------------------------------------------
       3              : !     gga91 exchange for a spin-unpolarized electronic system
       4              : !.....-----------------------------------------------------------------
       5              : !     input d: density
       6              : !           s:  abs(grad d)/(2*kf*d)
       7              : !           u:  (grad d)*grad(abs(grad d))/(d**2 * (2*kf)**3)
       8              : !           v: (laplacian d)/(d*(2*kf)**2)
       9              : !     output:  exchange energy (ex) and potential (vx) in ry.
      10              : 
      11              : !       kf=cbrt(3*pai**2*d).
      12              : !.....-----------------------------------------------------------------
      13              : CONTAINS
      14            0 :    SUBROUTINE xch91( &
      15              :       d,s,u,v, &
      16              :       exl,exg,vxl,vxg)
      17              : 
      18              :       IMPLICIT NONE
      19              : 
      20              :       REAL, INTENT (IN)  :: d,s,u,v
      21              :       REAL, INTENT (OUT) :: exl,exg,vxl,vxg
      22              : 
      23              :       REAL :: a,a1,a2,a3,a4,ax,b1,ex,f,fac,fs,fss, &
      24              :               p0,p1,p10,p11,p2,p3,p4,p5,p6,p7,p8,p9,s2,s3,s4, &
      25              :               thrd,thrd4,vx,a4s2,expsm
      26              : !.....-----------------------------------------------------------------
      27              : !     ..
      28              :       DATA a1,a2,a3,a4/0.19645,0.27430e0,0.15084e0,100.e0/
      29              :       DATA ax,a,b1/-0.7385588e0,7.7956e0,0.004e0/
      30              :       DATA thrd,thrd4/0.333333333333e0,1.33333333333e0/
      31              : !.....-----------------------------------------------------------------
      32              : !     expsm: argument of exponential-smallest.
      33            0 :       expsm=minexponent(expsm)/1.5
      34              : 
      35            0 :       fac = ax*d**thrd
      36            0 :       s2 = s*s
      37            0 :       s3 = s2*s
      38            0 :       s4 = s2*s2
      39            0 :       p0 = 1.e0/sqrt(1.e0+a*a*s2)
      40            0 :       p1 = log(a*s+1.e0/p0)
      41            0 :       a4s2=max(-a4*s2,expsm)
      42            0 :       p2 = exp(a4s2)
      43            0 :       p3 = 1.e0/ (1.e0+a1*s*p1+b1*s4)
      44            0 :       p4 = 1.e0 + a1*s*p1 + (a2-a3*p2)*s2
      45            0 :       f = p3*p4
      46            0 :       ex = fac*f
      47              : !  local exchange exl
      48              :       exl = fac
      49            0 :       exg = ex - exl
      50              : 
      51              : !  energy done. now the potential:
      52            0 :       p5 = b1*s2 - (a2-a3*p2)
      53            0 :       p6 = a1*s* (p1+a*s*p0)
      54            0 :       p7 = 2.e0* (a2-a3*p2) + 2.e0*a3*a4*s2*p2 - 4.e0*b1*s2*f
      55            0 :       fs = p3* (p3*p5*p6+p7)
      56            0 :       p8 = 2.e0*s* (b1-a3*a4*p2)
      57            0 :       p9 = a1*p1 + a*a1*s*p0* (3.e0-a*a*s2*p0*p0)
      58            0 :       p10 = 4.e0*a3*a4*s*p2* (2.e0-a4*s2) - 8.e0*b1*s*f - 4.e0*b1*s3*fs
      59            0 :       p11 = -p3*p3* (a1*p1+a*a1*s*p0+4.e0*b1*s3)
      60            0 :       fss = p3*p3* (p5*p9+p6*p8) + 2.e0*p3*p5*p6*p11 + p3*p10 + p7*p11
      61            0 :       vx = fac* (thrd4*f- (u-thrd4*s3)*fss-v*fs)
      62              : !  local exchange vxl:
      63            0 :       vxl = fac*thrd4
      64            0 :       vxg = vx - vxl
      65              : 
      66              : ! in ry and energy density.
      67            0 :       exl = exl*2.e0*d
      68            0 :       exg = exg*2.e0*d
      69            0 :       vxl = vxl*2.
      70            0 :       vxg = vxg*2.
      71              : 
      72            0 :    END SUBROUTINE xch91
      73              : END MODULE m_xch91