subroutine exchcor(n,h,den,rvxc,etxc,eexc,a,b,c,itype)
c  calculate exchange correlation potentials and energys
c    for density functional theory from electron density
c  den(n) is electron density * (4*pi*r**2)
c  rvxc(n) is returned as vxc * r
c    to the total energy
c  eexc is the total exchange energy (int(den*exc))
c  etxc is eexc - int(den*vxc)
c  a(n), b(n), c(n), d(n) are work arrays
c  if itype=15 (gga case) need to calculate additional functions:
c       gd=abs(grad(density))/density
c       ggd=laplacian(density)/density
c       gggd=grad(density)*grad(abs(grad(density)))/density**2
c       for r<r0 (fixed in parameter statement), replace density by Gaussian
c         form so that Laplacian does not diverge
c
       implicit real*8 (a-h,o-z)
       parameter (r0=0.01d0)
       dimension den(n),rvxc(n),a(n),b(n),c(n)
       fpi=16*datan(1.d0)
       i0=r0/h
       i0=max0(i0,2)
       r00=i0*h
c
c  Exchange-correlation contribution  exc-vxc
c
      if (itype.ne.15) then
        do i=1,n
        r=i*h
        r2=r*r
        rho=den(i)/(fpi*r2)
        call excor(rho,vxc,exc,.false.,itype)
        a(i)=den(i)*(exc-vxc)
        b(i)=den(i)*exc
        rvxc(i)=r*vxc
        enddo
      endif
      if (itype.eq.15) then
        do i=1,n
        r=i*h
        r2=r*r
        a(i)=den(i)/(fpi*r2)
        b(i)=ddlog(a(i))
        enddo
        call nderiv(h,b,c,n)
        call nderiv(h,c,rvxc,n)
c replace values for r<r0 (i<i0) by Gaussian form
        t1=c(i0)/r00
        t2=rvxc(i0)
        bb=-0.75d0*t1+0.25*t2
        cr2=0.125d0*(t2-t1)
        br=bb*r00
        cr3=cr2*r00
        br2=br*r00
        cr4=cr3*r00
        do i=1,i0
        x=i
        x=x/i0
        x2=x*x
        x3=x*x2
        x4=x2*x2
        a(i)=a(i0)*dexp(-br2*(x2-1.d0)+cr4*(x4-1.d0))
        c(i)=-2*br*x+4*cr3*x3
        rvxc(i)=-2*bb+12*cr2*x2
        enddo
c
        do i=1,n
        r=i*h
c Laplacian term:
        ggd=2*c(i)/r+rvxc(i)+c(i)**2
c double derivative term:
        gggd=c(i)*(rvxc(i)+c(i)**2)*dsign(1.d0,c(i))
c gradient term:
        gd=dabs(c(i))
        call excorp(a(i),gd,ggd,gggd,vxc,exc,.false.,itype)
        a(i)=den(i)*(exc-vxc)
        b(i)=den(i)*exc
        rvxc(i)=r*vxc
        enddo
        endif
c  calculate exchange-correlation contribution to the potential
c       etxc=dsum(n,a,1)*h
c       eexc=dsum(n,b,1)*h
        etxc=overint(n,h,a)
        eexc=overint(n,h,b)
        write(6,*) 'etxc,eexc = ',etxc,eexc
        return
        end