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