1 files changed, 76 insertions, 0 deletions

diff --git a/src/labat/ggaexc.f b/src/labat/ggaexc.f
new file mode 100644
index 0000000..135b42e
--- /dev/null
+++ b/src/labat/ggaexc.f
@@ -0,0 +1,76 @@
+      subroutine ggaexc(r,db,i,imin,imax,doned,exc)
+c  ggaexc calculates the exchange-correlation energy
+c  input r         : position coordinate (array)
+c  input db        : spin up and down density. (array)
+c  input i         : counter, i.e. present position = r(i)
+c  input imin      : min value of i
+c  input imax      : max value of i
+c  in/output doned : true if density is splined
+c  output exc      : exchange-correlation energy
+      implicit logical (a-z)
+      logical doned
+      double precision r,db,gradd,gradup,graddn,fk,fkup,fkdn,sk,g
+      double precision sup,sdn,t,exc,ec
+      double precision onethi,twothi,pi,conkf,d,zet,rs,exup,exdn
+      double precision ex,eclda,ecgga
+      integer i,imin,imax
+      dimension r(561),db(2,561)
+      dimension gradd(561),gradup(561),graddn(561),fk(561),fkup(561)
+      dimension fkdn(561),sk(561),g(561),sup(561),sdn(561),t(561)
+      common/gga/gradd,gradup,graddn,fkup,fkdn,sk,g,sup,sdn,t
+      onethi = 1.d0/3.d0
+      twothi = 2.d0/3.d0
+      pi = 4.d0*datan(1.d0)
+      conkf = (3.d0*pi**2)**onethi
+      d = db(1,i) + db(2,i)
+      call grad(r,db,i,imin,imax,doned,gradd(i),gradup(i),graddn(i))
+
+c  calculate the exchange energy, first spin up, then spin down
+
+      if (db(1,i).gt.1.d-100) then
+         fkup(i) = conkf*(2.d0*db(1,i))**onethi
+         sup(i) = dabs(gradup(i))/(2.d0*fkup(i)*db(1,i))
+         call exchen(2.d0*db(1,i),sup(i),exup)
+      else
+         fkup(i) = 0.d0
+         sup(i) = 0.d0
+         exup = 0.d0
+      endif
+      if (db(2,i).gt.1.d-100) then
+         fkdn(i) = conkf*(2.d0*db(2,i))**onethi
+         sdn(i) = dabs(graddn(i))/(2.d0*fkdn(i)*db(2,i))
+         call exchen(2.d0*db(2,i),sdn(i),exdn)
+      else
+         fkdn(i) = 0.d0
+         sdn(i) = 0.d0
+         exdn = 0.d0
+      endif
+      if (d.gt.1d-100) then
+         ex = (exup*db(1,i) + exdn*db(2,i))/d
+      else
+         ex = 0.d0
+      endif
+
+c  exchange energy done, now calculate the correlation energy
+
+      fk(i) = conkf*d**onethi
+      sk(i) = dsqrt(4.d0*fk(i)/pi)
+      if (d.gt.1.d-100) then
+         zet = (db(1,i) - db(2,i))/d
+      else
+         zet = 0.d0
+      endif
+      g(i) = ((1.d0+zet)**twothi+(1.d0-zet)**twothi)/2.d0
+      if (d.gt.1.d-18) then
+         rs = (3.d0/(4.d0*pi*d))**onethi
+         t(i) = dabs(gradd(i))/(d*2*sk(i)*g(i))
+         call ldaec(rs,zet,eclda)
+         call ggaec(rs,zet,t(i),eclda,ecgga)
+         ec = eclda + ecgga
+      else
+         ec = 0.d0
+      endif
+      exc = ex + ec
+
+      return
+      end