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

          Line data    Source code
       1             : MODULE m_corg91
       2             : !.....-----------------------------------------------------------------
       3             : !     gga91 correlation
       4             : !.....-----------------------------------------------------------------
       5             : CONTAINS
       6           0 :    SUBROUTINE corg91(fk,sk,gz,ec,ecrs,eczta,rs,zta,t,uu,vv,ww,h, &
       7             :                      dvcup,dvcdn)
       8             : !.....-----------------------------------------------------------------
       9             : !     input
      10             : !           rs: seitz radius
      11             : !         zta: relative spin polarization
      12             : !            t: abs(grad d)/(d*2.*ks*gz)
      13             : !           uu: (grad d)*grad(abs(grad d))/(d**2 * (2*ks*gz)**3)
      14             : !           vv: (laplacian d)/(d * (2*ks*gz)**2)
      15             : !           ww: (grad d)*(gradzta)/(d * (2*ks*gz)**2
      16             : !      output
      17             : !                  h: nonlocal part of correlation energy per electron
      18             : !          dvcup,-dn: nonlocal parts of correlation potentials.
      19             : 
      20             : !      with ks=sqrt(4*kf/pai), gz=[(1+zta)**(2/3)+(1-zta)**(2/3)]/2, &
      21             : !           kf=cbrt(3*pai**2*d).
      22             : !.....-----------------------------------------------------------------
      23             : !.....-----------------------------------------------------------------
      24             : !     .. previously untyped names ..
      25             :       IMPLICIT NONE
      26             : 
      27             :       REAL :: dvcdn,dvcup,ec,ecrs,eczta,fk,gz,h,rs,sk,t,uu,vv,ww,zta,a4, &
      28             :               alf,argmx,b,b2,b2fac,bec,bet,bg,c1,c2,c3,c4,c5,c6,cc,cc0, &
      29             :               ccrs,coeff,comm,cx,delt,fact0,fact1,fact2,fact3,fact4,fact5, &
      30             :               gm,gz3,gz4,h0,h0b,h0bt,h0rs,h0rst,h0t,h0tt,h0z,h0zt,h1,h1rs, &
      31             :               h1rst,h1t,h1tt,h1z,h1zt,hrs,hrst,ht,htt,hz,hzt,pon,pref,q4, &
      32             :               q5,q6,q7,q8,q9,r0,r1,r2,r3,r4,rs2,rs3,rsthrd,t2,t4,t6,thrd2, &
      33             :               thrdm,xnu,r1t2,expsm
      34             : !     ..
      35             :       DATA xnu,cc0,cx,alf/15.75592e0,0.00423500,-0.001667212e0,0.0900/
      36             :       DATA c1,c2,c3,c4/0.00256800,0.02326600,7.389e-6,8.723e0/
      37             :       DATA c5,c6,a4/0.472e0,7.389e-2,100.00/
      38             :       DATA thrdm,thrd2/-0.333333333333e0,0.666666666667e0/
      39             : !.....-----------------------------------------------------------------
      40             : !     expsm: argument of exponential-smallest.
      41           0 :       expsm=minexponent(expsm)/1.5
      42           0 :       argmx = 174.0
      43           0 :       bet = xnu*cc0
      44           0 :       delt = 2.e0*alf/bet
      45           0 :       gz3 = gz**3
      46           0 :       gz4 = gz3*gz
      47           0 :       pon = -delt*ec/ (gz3*bet)
      48           0 :       IF (pon > argmx) THEN
      49             :          b = 0.
      50             :       ELSE
      51           0 :          b = delt/ (exp(pon)-1.00)
      52             :       ENDIF
      53           0 :       b2 = b*b
      54           0 :       t2 = t*t
      55           0 :       t4 = t2*t2
      56           0 :       t6 = t4*t2
      57           0 :       rs2 = rs*rs
      58           0 :       rs3 = rs2*rs
      59           0 :       q4 = 1.00 + b*t2
      60           0 :       q5 = 1.00 + b*t2 + b2*t4
      61           0 :       q6 = c1 + c2*rs + c3*rs2
      62           0 :       q7 = 1.00 + c4*rs + c5*rs2 + c6*rs3
      63           0 :       cc = -cx + q6/q7
      64           0 :       r0 = (sk/fk)**2
      65           0 :       r1 = a4*r0*gz4
      66           0 :       coeff = cc - cc0 - 3.e0*cx/7.00
      67           0 :       r2 = xnu*coeff*gz3
      68             : ! tagu
      69           0 :       r1t2=max(-r1*t2,expsm)
      70           0 :       r3 = exp(r1t2)
      71             : ! tagu
      72           0 :       h0 = gz3* (bet/delt)*log(1.00+delt*q4*t2/q5)
      73           0 :       h1 = r3*r2*t2
      74           0 :       h = h0 + h1
      75             : !  local correlation option:
      76             : !     h=0.0e0
      77             : 
      78             : !  energy done. now the potential:
      79             : 
      80           0 :       ccrs = (c2+2.*c3*rs)/q7 - q6* (c4+2.*c5*rs+3.*c6*rs2)/q7**2
      81           0 :       rsthrd = rs/3.e0
      82           0 :       r4 = rsthrd*ccrs/coeff
      83           0 :       gm = ((1.00+zta)**thrdm- (1.00-zta)**thrdm)/3.00
      84           0 :       IF (pon > argmx) THEN
      85             :          b2fac = 0.
      86             :       ELSE
      87           0 :          b2fac = b2* (delt/b+1.00)
      88             :       ENDIF
      89           0 :       bg = -3.e0*ec*b2fac/ (bet*gz4)
      90           0 :       bec = b2fac/ (bet*gz3)
      91           0 :       q8 = q5*q5 + delt*q4*q5*t2
      92           0 :       q9 = 1.00 + 2.00*b*t2
      93           0 :       h0b = -bet*gz3*b*t6* (2.e0+b*t2)/q8
      94           0 :       h0rs = -rsthrd*h0b*bec*ecrs
      95           0 :       fact0 = 2.e0*delt - 6.00*b
      96           0 :       fact1 = q5*q9 + q4*q9*q9
      97           0 :       h0bt = 2.e0*bet*gz3*t4* ((q4*q5*fact0-delt*fact1)/q8)/q8
      98           0 :       h0rst = rsthrd*t2*h0bt*bec*ecrs
      99           0 :       h0z = 3.e0*gm*h0/gz + h0b* (bg*gm+bec*eczta)
     100           0 :       h0t = 2.*bet*gz3*q9/q8
     101           0 :       h0zt = 3.e0*gm*h0t/gz + h0bt* (bg*gm+bec*eczta)
     102           0 :       fact2 = q4*q5 + b*t2* (q4*q9+q5)
     103           0 :       fact3 = 2.e0*b*q5*q9 + delt*fact2
     104           0 :       h0tt = 4.e0*bet*gz3*t* (2.e0*b/q8- (q9* (fact3/q8))/q8)
     105           0 :       h1rs = r3*r2*t2* (-r4+r1*t2/3.e0)
     106           0 :       fact4 = 2.e0 - r1*t2
     107           0 :       h1rst = r3*r2*t2* (2.e0*r4* (1.00-r1*t2)-thrd2*r1*t2*fact4)
     108           0 :       h1z = gm*r3*r2*t2* (3.e0-4.e0*r1*t2)/gz
     109           0 :       h1t = 2.e0*r3*r2* (1.00-r1*t2)
     110           0 :       h1zt = 2.e0*gm*r3*r2* (3.e0-11.00*r1*t2+4.e0*r1*r1*t4)/gz
     111           0 :       h1tt = 4.e0*r3*r2*r1*t* (-2.e0+r1*t2)
     112           0 :       hrs = h0rs + h1rs
     113           0 :       hrst = h0rst + h1rst
     114           0 :       ht = h0t + h1t
     115           0 :       htt = h0tt + h1tt
     116           0 :       hz = h0z + h1z
     117           0 :       hzt = h0zt + h1zt
     118           0 :       comm = h + hrs + hrst + t2*ht/6.00 + 7.00*t2*t*htt/6.00
     119           0 :       pref = hz - gm*t2*ht/gz
     120           0 :       fact5 = gm* (2.e0*ht+t*htt)/gz
     121           0 :       comm = comm - pref*zta - uu*htt - vv*ht - ww* (hzt-fact5)
     122           0 :       dvcup = comm + pref
     123           0 :       dvcdn = comm - pref
     124             : 
     125             : !  local correlation option:
     126             : !     dvcup=0.0e0
     127             : !     dvcdn=0.0e0
     128             : 
     129           0 :    END SUBROUTINE corg91
     130             : END MODULE m_corg91

Generated by: LCOV version 1.13