1 files changed, 82 insertions, 0 deletions

diff --git a/src/labat/ggauxc.f b/src/labat/ggauxc.f
new file mode 100644
index 0000000..d5fcb61
--- /dev/null
+++ b/src/labat/ggauxc.f
@@ -0,0 +1,82 @@
+      subroutine ggauxc(r,db,zet,i,imin,imax,donela,donegr,donez,
+     j     uxcup,uxcdn)
+c  ggauxc calculates the exchange-correlation potential by taking
+c  the functional derivative of the energy with respect to the density
+c  input r          : position coordinate (array)
+c  input db         : spin up and down density. (array)
+c  input zet        : relative spin polarization (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 donela : true if r^2*grad(density) is splined
+c  in/output donegr : true if abs(grad(density)) is splined
+c  in/output donez  : true if zet is splined
+c  output uxcup     : spin-up exchange-correlation potential
+c  output uxcdn     : spin-down exchange-correlation potential
+      implicit logical (a-z)
+      double precision r,db,zet,uxcup,uxcdn,gradd,gradup,graddn
+      double precision fkup,fkdn,sk,g,sup,sdn,t,tempz,z2,gradz
+      double precision onethi,pi,d,rs,lapld,laplup,lapldn,deriv
+      double precision gagd,gagup,gagdn,uu,vv,ww,uup,vup,udn,vdn
+      double precision uclcup,uclcdn,ucgaup,ucgadn,uxup,uxdn
+      double precision ucup,ucdn,ec,ecrs,eczet
+      integer i,imin,imax,j
+      logical donela,donegr,donez
+      common/gga/gradd,gradup,graddn,fkup,fkdn,sk,g,sup,sdn,t
+      common/ggaz/tempz,z2
+      dimension r(561),db(2,561),zet(561),tempz(2,561),z2(2,561)
+      dimension gradd(561),gradup(561),graddn(561),fkup(561)
+      dimension fkdn(561),sk(561),g(561),sup(561),sdn(561),t(561)
+      onethi = 1.d0/3.d0
+      pi = 4.d0*datan(1.d0)
+      call laplac(r,gradup,graddn,i,imin,imax,donela,lapld,laplup,
+     j     lapldn)
+      call grabgr(r,gradup,graddn,i,imin,imax,donegr,gagd,gagup,gagdn)
+
+c  calculate the exchange potential
+      if (db(1,i).gt.1.d-100) then
+         uup = gradup(i)*gagup/(db(1,i)**2*(2.d0*fkup(i))**3)
+         vup = laplup/(db(1,i)*(2.d0*fkup(i))**2)
+         call exchpt(2.d0*db(1,i),sup(i),uup,vup,uxup)
+      else
+         uxup = 0.d0
+      endif
+      if (db(2,i).gt.1.d-100) then
+         udn = graddn(i)*gagdn/(db(2,i)**2*(2.d0*fkdn(i))**3)
+         vdn = lapldn/(db(2,i)*(2.d0*fkdn(i))**2)
+         call exchpt(2.d0*db(2,i),sdn(i),udn,vdn,uxdn)
+      else
+         uxdn = 0.d0
+      endif
+
+c  exchange potential done, now calculate the correlation potential
+
+      if (.not.donez) then
+         do 10 j = imin,imax
+            tempz(1,j) = zet(j)
+ 10      continue
+         call spline(r,imin+1,imax,tempz,z2)
+         donez = .true.
+      endif 
+      gradz = deriv(r,tempz,z2,1,i,imin,imax)
+      d = db(1,i) + db(2,i)
+      if (d.gt.1.d-18) then
+         rs = (3.d0/(4.d0*pi*d))**onethi
+         uu = gradd(i)*gagd/(d**2*(2.d0*sk(i)*g(i))**3)
+         vv = lapld/(d*(2.d0*sk(i)*g(i))**2)
+         ww = gradd(i)*gradz/(d*(2.d0*sk(i)*g(i))**2)
+         call ldauc(rs,zet(i),ec,ecrs,eczet,uclcup,uclcdn)
+         call ggauc(rs,zet(i),t(i),uu,vv,ww,ec,ecrs,eczet,ucgaup,
+     j        ucgadn)
+         ucup = uclcup + ucgaup
+         ucdn = uclcdn + ucgadn
+      else
+         ucdn = 0.d0
+         ucup = 0.d0
+      endif
+      uxcup = uxup + ucup
+      uxcdn = uxdn + ucdn
+
+      return
+      end
+