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

c      with ks=sqrt(4*kf/pai), gz=[(1+zta)**(2/3)+(1-zta)**(2/3)]/2, &
c           kf=cbrt(3*pai**2*d).
c.....-----------------------------------------------------------------
c.....-----------------------------------------------------------------
c     .. previously untyped names ..
      IMPLICIT NONE

      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
c     ..
      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/
c.....-----------------------------------------------------------------
c     expsm: argument of exponential-smallest.
41
      expsm=minexponent(expsm)/1.5
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
      argmx = 174.0
      bet = xnu*cc0
      delt = 2.e0*alf/bet
      gz3 = gz**3
      gz4 = gz3*gz
      pon = -delt*ec/ (gz3*bet)
      IF (pon.gt.argmx) THEN
          b = 0.
      ELSE
          b = delt/ (exp(pon)-1.00)
      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
c+tagu
      r1t2=max(-r1*t2,expsm)
      r3 = exp(r1t2)
c-tagu
      h0 = gz3* (bet/delt)*log(1.00+delt*q4*t2/q5)
      h1 = r3*r2*t2
      h = h0 + h1
c  local correlation option:
c     h=0.0e0

c  energy done. now the potential:

      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
      IF (pon.gt.argmx) THEN
          b2fac = 0.
      ELSE
          b2fac = b2* (delt/b+1.00)
      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

c  local correlation option:
c     dvcup=0.0e0
c     dvcdn=0.0e0

      END SUBROUTINE corg91
      END MODULE m_corg91