1 files changed, 43 insertions, 0 deletions

diff --git a/src/labat/ggaec.f b/src/labat/ggaec.f
new file mode 100644
index 0000000..3177859
--- /dev/null
+++ b/src/labat/ggaec.f
@@ -0,0 +1,43 @@
+      subroutine ggaec(rs,zet,t,ec,h)
+c  called by subroutine ggaexc
+c  gga91 correlation
+c  input rs  : seitz radius
+c  input zet : relative spin polarization
+c  input t   : abs(grad d)/(d*2.*ks*g)
+c  input     : correlation energy per electron (ec)
+c  output h  : nonlocal part of correlation energy per electron
+      implicit double precision (a-h,o-z)
+      data xnu,cc0,cx,alf/15.75592d0,0.004235d0,-0.001667212d0,0.09d0/
+      data c1,c2,c3,c4/0.002568d0,0.023266d0,7.389d-6,8.723d0/
+      data c5,c6,a4/0.472d0,7.389d-2,100.d0/
+      data thrd2/0.666666666667d0/
+      pi = 4.d0*datan(1.d0)
+      fk = 1.91915829d0/rs
+      sk = dsqrt(4.d0*fk/pi)
+      g = ((1.d0+zet)**thrd2+(1.d0-zet)**thrd2)/2.d0
+      bet = xnu*cc0
+      delt = 2.d0*alf/bet
+      g3 = g**3
+      g4 = g3*g
+      pon = -delt*ec/(g3*bet)
+      b = delt/(dexp(pon)-1.d0)
+      b2 = b*b
+      t2 = t*t
+      t4 = t2*t2
+      rs2 = rs*rs
+      rs3 = rs2*rs
+      q4 = 1.d0+b*t2
+      q5 = 1.d0+b*t2+b2*t4
+      q6 = c1+c2*rs+c3*rs2
+      q7 = 1.d0+c4*rs+c5*rs2+c6*rs3
+      cc = -cx + q6/q7
+      r0 = (sk/fk)**2
+      r1 = a4*r0*g4
+      coeff = cc-cc0-3.d0*cx/7.d0
+      r2 = xnu*coeff*g3
+      r3 = dexp(-r1*t2)
+      h0 = g3*(bet/delt)*dlog(1.d0+delt*q4*t2/q5)
+      h1 = r3*r2*t2
+      h = h0 + h1
+      return
+      end