corg91.f90 4.47 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
MODULE m_corg91
!.....-----------------------------------------------------------------
!     gga91 correlation
!.....-----------------------------------------------------------------
CONTAINS
   SUBROUTINE corg91(fk,sk,gz,ec,ecrs,eczta,rs,zta,t,uu,vv,ww,h, &
                     dvcup,dvcdn)
!.....-----------------------------------------------------------------
!     input
!           rs: seitz radius
!         zta: relative spin polarization
!            t: abs(grad d)/(d*2.*ks*gz)
!           uu: (grad d)*grad(abs(grad d))/(d**2 * (2*ks*gz)**3)
!           vv: (laplacian d)/(d * (2*ks*gz)**2)
!           ww: (grad d)*(gradzta)/(d * (2*ks*gz)**2
!      output
!                  h: nonlocal part of correlation energy per electron
!          dvcup,-dn: nonlocal parts of correlation potentials.
19

Matthias Redies's avatar
Matthias Redies committed
20 21 22 23 24
!      with ks=sqrt(4*kf/pai), gz=[(1+zta)**(2/3)+(1-zta)**(2/3)]/2, &
!           kf=cbrt(3*pai**2*d).
!.....-----------------------------------------------------------------
!.....-----------------------------------------------------------------
!     .. previously untyped names ..
25 26
      IMPLICIT NONE

Matthias Redies's avatar
Matthias Redies committed
27 28 29 30 31 32 33 34
      REAL :: dvcdn,dvcup,ec,ecrs,eczta,fk,gz,h,rs,sk,t,uu,vv,ww,zta,a4, &
              alf,argmx,b,b2,b2fac,bec,bet,bg,c1,c2,c3,c4,c5,c6,cc,cc0, &
              ccrs,coeff,comm,cx,delt,fact0,fact1,fact2,fact3,fact4,fact5, &
              gm,gz3,gz4,h0,h0b,h0bt,h0rs,h0rst,h0t,h0tt,h0z,h0zt,h1,h1rs, &
              h1rst,h1t,h1tt,h1z,h1zt,hrs,hrst,ht,htt,hz,hzt,pon,pref,q4, &
              q5,q6,q7,q8,q9,r0,r1,r2,r3,r4,rs2,rs3,rsthrd,t2,t4,t6,thrd2, &
              thrdm,xnu,r1t2,expsm
!     ..
35 36 37 38
      DATA xnu,cc0,cx,alf/15.75592e0,0.00423500,-0.001667212e0,0.0900/
      DATA c1,c2,c3,c4/0.00256800,0.02326600,7.389e-6,8.723e0/
      DATA c5,c6,a4/0.472e0,7.389e-2,100.00/
      DATA thrdm,thrd2/-0.333333333333e0,0.666666666667e0/
Matthias Redies's avatar
Matthias Redies committed
39 40
!.....-----------------------------------------------------------------
!     expsm: argument of exponential-smallest.
41
      expsm=minexponent(expsm)/1.5
42 43 44 45 46 47
      argmx = 174.0
      bet = xnu*cc0
      delt = 2.e0*alf/bet
      gz3 = gz**3
      gz4 = gz3*gz
      pon = -delt*ec/ (gz3*bet)
Matthias Redies's avatar
Matthias Redies committed
48 49
      IF (pon > argmx) THEN
         b = 0.
50
      ELSE
Matthias Redies's avatar
Matthias Redies committed
51
         b = delt/ (exp(pon)-1.00)
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
      ENDIF
      b2 = b*b
      t2 = t*t
      t4 = t2*t2
      t6 = t4*t2
      rs2 = rs*rs
      rs3 = rs2*rs
      q4 = 1.00 + b*t2
      q5 = 1.00 + b*t2 + b2*t4
      q6 = c1 + c2*rs + c3*rs2
      q7 = 1.00 + c4*rs + c5*rs2 + c6*rs3
      cc = -cx + q6/q7
      r0 = (sk/fk)**2
      r1 = a4*r0*gz4
      coeff = cc - cc0 - 3.e0*cx/7.00
      r2 = xnu*coeff*gz3
Matthias Redies's avatar
Matthias Redies committed
68
! tagu
69 70
      r1t2=max(-r1*t2,expsm)
      r3 = exp(r1t2)
Matthias Redies's avatar
Matthias Redies committed
71
! tagu
72 73 74
      h0 = gz3* (bet/delt)*log(1.00+delt*q4*t2/q5)
      h1 = r3*r2*t2
      h = h0 + h1
Matthias Redies's avatar
Matthias Redies committed
75 76
!  local correlation option:
!     h=0.0e0
77

Matthias Redies's avatar
Matthias Redies committed
78
!  energy done. now the potential:
79 80 81 82 83

      ccrs = (c2+2.*c3*rs)/q7 - q6* (c4+2.*c5*rs+3.*c6*rs2)/q7**2
      rsthrd = rs/3.e0
      r4 = rsthrd*ccrs/coeff
      gm = ((1.00+zta)**thrdm- (1.00-zta)**thrdm)/3.00
Matthias Redies's avatar
Matthias Redies committed
84 85
      IF (pon > argmx) THEN
         b2fac = 0.
86
      ELSE
Matthias Redies's avatar
Matthias Redies committed
87
         b2fac = b2* (delt/b+1.00)
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
      ENDIF
      bg = -3.e0*ec*b2fac/ (bet*gz4)
      bec = b2fac/ (bet*gz3)
      q8 = q5*q5 + delt*q4*q5*t2
      q9 = 1.00 + 2.00*b*t2
      h0b = -bet*gz3*b*t6* (2.e0+b*t2)/q8
      h0rs = -rsthrd*h0b*bec*ecrs
      fact0 = 2.e0*delt - 6.00*b
      fact1 = q5*q9 + q4*q9*q9
      h0bt = 2.e0*bet*gz3*t4* ((q4*q5*fact0-delt*fact1)/q8)/q8
      h0rst = rsthrd*t2*h0bt*bec*ecrs
      h0z = 3.e0*gm*h0/gz + h0b* (bg*gm+bec*eczta)
      h0t = 2.*bet*gz3*q9/q8
      h0zt = 3.e0*gm*h0t/gz + h0bt* (bg*gm+bec*eczta)
      fact2 = q4*q5 + b*t2* (q4*q9+q5)
      fact3 = 2.e0*b*q5*q9 + delt*fact2
      h0tt = 4.e0*bet*gz3*t* (2.e0*b/q8- (q9* (fact3/q8))/q8)
      h1rs = r3*r2*t2* (-r4+r1*t2/3.e0)
      fact4 = 2.e0 - r1*t2
      h1rst = r3*r2*t2* (2.e0*r4* (1.00-r1*t2)-thrd2*r1*t2*fact4)
      h1z = gm*r3*r2*t2* (3.e0-4.e0*r1*t2)/gz
      h1t = 2.e0*r3*r2* (1.00-r1*t2)
      h1zt = 2.e0*gm*r3*r2* (3.e0-11.00*r1*t2+4.e0*r1*r1*t4)/gz
      h1tt = 4.e0*r3*r2*r1*t* (-2.e0+r1*t2)
      hrs = h0rs + h1rs
      hrst = h0rst + h1rst
      ht = h0t + h1t
      htt = h0tt + h1tt
      hz = h0z + h1z
      hzt = h0zt + h1zt
      comm = h + hrs + hrst + t2*ht/6.00 + 7.00*t2*t*htt/6.00
      pref = hz - gm*t2*ht/gz
      fact5 = gm* (2.e0*ht+t*htt)/gz
      comm = comm - pref*zta - uu*htt - vv*ht - ww* (hzt-fact5)
      dvcup = comm + pref
      dvcdn = comm - pref

Matthias Redies's avatar
Matthias Redies committed
125 126 127
!  local correlation option:
!     dvcup=0.0e0
!     dvcdn=0.0e0
128

Matthias Redies's avatar
Matthias Redies committed
129 130
   END SUBROUTINE corg91
END MODULE m_corg91