subroutine energyanalysis(nz,n,h,mxo,occ,norbit,qcore,qvale,
     x nlst,llst,iscore,nmx,r,wfn,itype,coreden,valeden,den,tc,tv,
     x ccoul,vcoul,cext,vext,cexc,vexc,cvcoul,texc,ec,ev,ecv,etotal,
     x evale,mxwk,wk)
c  program to calculate the core and valence contributions to the
c   total energy
c   nz = nuclear charge
c   n = number of mesh points
c   h = step size between mesh points ( in bohr)
c   mxo = maximum number of orbitals
c   norbit = actual number of orbital
c   iscore(mxo) = 0 for core state , 1 for valence state
c   llst(mxo) = l for each orbital
c   nlst(mxo) = principal quantum number for each orbital
c   occ(mxo) = occupation of each orbital
c   itype = exchange-correlation type index
c   coreden = core density
c   valeden = valence density
c   den = coreden + valeden
c   tc= core kinetic energy
c   tv= valence kinetic energy
c   ccoul = core coulomb energy
c   vcoul = valence coulomb energy
c   cext = core external energy
c   vext = valence external energy
c   cexc = core exchange-correlation energy
c   vexc = valence exchange-correlation energy
c   cvcoul = core-valence coulomb interaction energy
c   texc = total exchange-correlation energy
c   ec = total core energy
c   ev = total valence energy
c   ecv = total core-valence interaction energy
c   etotal = total energy
c   evale = tv + vcoul + vext + cvcoul + texc - cexc
c   reference: U. von Barth and C. D. Gelatt, Phys. Rev. B 21, 2222(1980)
c 
      implicit real*8 (a-h,o-z)
      dimension occ(mxo),wfn(nmx,mxo),den(nmx)
      dimension r(nmx),coreden(nmx),valeden(nmx),wk(nmx,mxwk)
      dimension nlst(mxo),llst(mxo),iscore(mxo)
      noz = 0
      if (mxwk.lt.7) then
       write(6,*) 'error in energyanalysis - mxwk must be gt 7',mxwk
       stop
      endif
      !call dcopy (n,0.d0,0,coreden,1)
      !call dcopy (n,0.d0,0,valeden,1)
       coreden=0.d0
       valeden=0.d0
      do io=1,norbit
      if (iscore(io).eq.0) then
         do i=1,n
         coreden(i)=coreden(i)+occ(io)*(wfn(i,io)**2)
         enddo
      endif
      if (iscore(io).eq.1) then
         do i=1,n
         valeden(i)=valeden(i)+occ(io)*(wfn(i,io)**2)
         enddo
      endif
      enddo
      do i=1,n
       den(i)=coreden(i)+valeden(i)
      enddo
c
      etotal=0
      ec=0
      ev=0
      ecv=0
c kinetic energy
      tc=0
      tv=0
      do io=1,norbit
      if (iscore(io).eq.0) then
      call kinetic(n,h,llst(io),r,wfn(1,io),wfn(1,io),wk,wk(1,2),
     x   wk(1,3),ekin)
        tc=tc+occ(io)*ekin
        endif
      if (iscore(io).eq.1) then
      call kinetic(n,h,llst(io),r,wfn(1,io),wfn(1,io),wk,wk(1,2),
     x   wk(1,3),ekin)
        tv=tv+occ(io)*ekin
        endif
      enddo
      write(6,*) 'core and valence kinetic energies', tc,tv
c external potential interaction
      do i=1,n
       wk(i,6)=coreden(i)/r(i)
       wk(i,7)=valeden(i)/r(i)
      enddo
      cext=-2*nz*overint(n,h,wk(1,6))
      vext=-2*nz*overint(n,h,wk(1,7))
      write(6,*) 'core and valence external potential energies',cext,vext
c 
      qcore=overint(n,h,coreden)
      qvale=overint(n,h,valeden)
      qtot=qcore+qvale
      qchk=overint(n,h,den)
      write(6,*) 'qcore,qvale,qtot,qchk=', qcore,qvale,qtot,qchk
      call poisson(n,h,qcore,coreden,wk,wk(1,2),wk(1,3),ccoul)
      chkcvcoul1=overlap(n,h,wk(1,7),wk)
      call exchcor(n,h,coreden,wk,dum,cexc,wk(1,2),wk(1,3),wk(1,4),itype)
      call poisson(n,h,qvale,valeden,wk,wk(1,2),wk(1,3),vcoul)
      chkcvcoul2=overlap(n,h,wk(1,6),wk)
      call exchcor(n,h,valeden,wk,dum,vexc,wk(1,2),wk(1,3),wk(1,4),itype)
      call poisson(n,h,qtot,den,wk,wk(1,2),wk(1,3),cvcoul)
      call exchcor(n,h,den,wk,dum,texc,wk(1,2),wk(1,3),wk(1,4),itype)
      cvcoul=cvcoul-vcoul-ccoul
      write(6,*) 'core , valence and interaction coulomb energies'
      write(6,*) ccoul,vcoul,cvcoul,chkcvcoul1,chkcvcoul2
      write(6,*) 'core , valence and interaction exc energies'
      write(6,*) cexc,vexc,texc
      ec=tc+ccoul+cext+cexc
      ev=tv+vcoul+vext+vexc
      ecv=cvcoul+texc-cexc-vexc
      etotal=ec+ev+ecv
      evale=tv+vcoul+vext+cvcoul+texc-cexc
      write(6,*) 'ec = ',ec
      write(6,*) 'ev = ',ev
      write(6,*) 'ecv = ', ecv
      write(6,*) 'etotal = ',etotal
      write(6,*) 'evale = ',evale
      return
      end