      subroutine excor(deni,vxc,exc,relat,itype)
      implicit real*8 (a-h,o-z)
c relativistic:
      parameter (balpha=137.03602d0,pi=3.14159265359d0)
c VWN:
      parameter (x0=-0.10498d0,a=0.0621814d0,b=3.72744d0,c=12.9352d0)
      parameter (xx0=x0*x0+b*x0+c,q=6.152029831d0)
      parameter (c1=2.d0*b/q,c2=b*x0/xx0,c3=2.d0*(b+2.d0*x0)/q)
c Perdew-Wang
      parameter (pa=0.031091d0,palp1=0.21370d0,pb1=7.5957d0,
     x pb2=3.5876d0,pb3=1.6382d0,pb4=0.49294d0,tpa=2.d0*pa)
      logical relat
      data init/0/
c input density in units of electrons/bohr**3
c output vxc potential in ryd
c output exc total exchange-correlation energy/den
c uses Kohn-Sham exchange
c adds relativistic exchange iff relat=.true.
c    Rajagopal, J.Phys.C 11, L943 (1978), 
c    MacDonald & Vosko, J.Phys.C 12, 2977 (1979)
c if itype eq 14 use Perdew & Wang, PRB, June 15, 1992
c if itype lt 13 use Vosko, Wilk, Nusair correlation,Can.J.Phys.58,1200(1980) 
c if itype eq 13 use Wigner exchange correlation, Pickett/Bylander, Kleinman 
c     B-K : Phys.Rev.B  vol.27 no.6 15 Mar 83
c     Pickett : Comp. Phys. Reports 9 (1989) p.115-198
c
c     q=dsqrt(4.d0*c-b*b)
c     q=dsqrt(4.d0*c-b*b)
c     c1=2.d0*b/q
c     c2=b*x0/xx0
c     c3=2.d0*(b+2.d0*x0)/q
      if (init.ne.itype) then
      if (itype.eq.13) write(6,*) 'Wigner correlation'
      if (itype.gt.13) write(6,*) 'Perdew-Wang correlation'
      if (itype.lt.13) write(6,*) 'Ceperley-WVN correlation'
      init=itype
      endif
c  inelegant fix for small densities
      den=dmax1(deni,1.d-68)
      fa=(3.d0*pi*pi*den)**0.333333333333333333d0
      vxnr=-2.d0*fa/pi
      if (itype.ge.13) then
      rs=(3.d0/(4.d0*pi*den))**0.33333333333333333d0
      if (itype.eq.13) then
c -- Wigner correlation --
      ec=-0.88d0/(7.79d0+rs)
      vc=-0.88d0*(4.d0*rs/3.d0+7.79d0)/((7.79d0+rs)**2)
      else
c -- Perdew Wang correlation --
      rsh=dsqrt(rs)
      ff=pb1*rsh+rs*(pb2+pb3*rsh+pb4*rs)
      ffp=0.5d0*pb1*rsh+rs*(pb2+1.5d0*pb3*rsh+2.d0*pb4*rs)
      xln=dlog(1.d0+1.d0/(tpa*ff))
      ec=-tpa*(1.d0+palp1*rs)*xln
      vc=-tpa*(1.d0+0.6666666666666d0*palp1*rs)*xln
     x-tpa*(1.d0+palp1*rs)*ffp/(3.d0*ff*(1.d0+tpa*ff))
      ec=2.d0*ec
      vc=2.d0*vc
      endif
      else
      x=(0.75d0/(den*pi))**0.1666666666666666667d0
      xx=x*x+b*x+c
      arg=q/(2.d0*x+b)
      ec=a*(dlog((x*x)/xx)+c1*datan(arg)
     x   -c2*(dlog(((x-x0)*(x-x0))/xx)+c3*datan(arg)))
      vc=ec-0.333333333333d0*a*(c/xx-(b*x0*x)/((x-x0)*xx))
      endif
      vxc=vxnr+vc
      exnr=0.75d0*vxnr
      exc=exnr+ec
      if(relat) then
      beta=fa/balpha
      bgam=dsqrt(1.d0+beta*beta)
      if(beta.gt.1.d-3) then
      rfac=(1.d0-dlog(beta+bgam)/(beta*bgam))   
      rrfac=rfac*(bgam/beta)
      else
      rfac=0.666666666666666666666666d0*(beta*beta)
      rrfac=0.666666666666666666666666d0*beta
      endif
      vxr=-1.5d0*vxnr*rfac
      exr=-1.5d0*exnr*(rrfac*rrfac)
      vxc=vxc+vxr
      exc=exc+exr
      endif
      return
      end
