subroutine potential(n,h,nz,q,den,rv,ecoul,etxc,eexc,
     x v0,v0p,a,b,c,d,itype)
c  calculate veff for density functional theory from electron density
c  nz is nuclear charge
c  den(n) is electron density * (4*pi*r**2)
c  rv(n) is returned as veff * r
c  q is total electron charge
c  ecoul and etxc are coulomb and exchange-correlation contributions
c    to the total energy
c  eexc is the total exchange energy (int(den*exc))
c  v0=v(0)
c  v0p=dv/dr(0)
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),rv(n),a(n),b(n),c(n),d(n)
       i0=r0/h
       i0=max0(i0,2)
       r00=i0*h
       fpi=16*datan(1.d0)
c      
c  Coulomb contribution
c
       call poisson(n,h,q,den,rv,a,b,ecoul)
c      do i=1,n
c      write(6,*) i, 'poisson',rv(i)
c      enddo
c
c
c  Exchange-correlation contribution  exc-vxc
c
        call exchcor(n,h,den,d,etxc,eexc,a,b,c,itype)
c
        th=1
        !call daxpy(n,th,d,1,rv,1)
        rv=rv+d
c
c  add in nuclear contribution and calculate v0 and v0p
c
        do i=1,15
        r=i*h
        a(i)=rv(i)/r
        enddo
        call nderiv(h,a,b,15)
        v0=a(3)+3*(a(1)-a(2))
        v0p=b(3)+3*(b(1)-b(2))
        write(6,*) 'v0, v0p = ', v0,v0p
c       rewind 8
        do i=1,n
        rv(i)=rv(i)-2*nz
c       r=i*h
c       write(8,800) i,r,den(i),rv(i)
c 800   format(i5,1p3e15.7)
        enddo
        return
        end