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