subroutine fcdfatom(nz,nps,npp,npd,npf,npg,n,h,mxnl,nl,mxo,occ,
     x norbit,nlst,llst,no,qf,nmx,r,wfn,eig,den,denout,rv,itype,mxloop,
     x qcore,coreden,totden,iscore,iscr1,iscr2,mxwrk,w,e,small)
c  frozen core version of dfatom
c  program to calculate the self-consistent density functional
c    atom ground state fo atom with atomic number nz
c  iscr1,iscr2 are scratch units
c  no(l+1) is returned as expected number of bound states for l 
c    for self-consistent potential rv
c    iscore(io) = 0 for core state, 1 for valence state
c  configuration set before calling this subroutine
c  den, denout are valence density only
      implicit real*8 (a-h,o-z)
      parameter (small0=1.d-15,rimix=0.2)
      parameter (niter=100,mxnz=10)
      dimension occ(mxo),wfn(nmx,mxo),eig(mxo)
      dimension r(nmx),den(nmx),rv(nmx),denout(nmx),coreden(nmx)
      dimension totden(nmx),w(nmx,mxwrk),e(mxwrk)
      dimension nl(mxnl,mxnl),nlst(mxo),llst(mxo),no(mxnl),iscore(mxo)
      if (n.gt.nmx) then
       write(6,*) 'error in dfatom -- n,nmx = ',n,nmx
       stop
       endif
      rmix=rimix
       if (norbit.gt.mxo) then
        write(6,*) 'error in fcdfatom -- norbit,mxo=',norbit,mxo
        stop
        endif
c
       electrons=0.d0
       do io=1,norbit      
       if(iscore(io).eq.1) then
       write(6,620) nlst(io),llst(io),occ(io)
  620 format(i2,1x,i2,4x,1pe15.7)
       electrons=electrons+occ(io)
       endif
       enddo
      qf=nz-electrons-qcore
      write(6,*)
      write(6,*) 'nuclear charge    = ' , nz
      write(6,*) 'core charge    = ' , qcore
      write(6,*) 'valence electronic charge = ', electrons
      write(6,*) 'net charge        = ', qf
c
c  check charge
c
      qcal=overint(n,h,den)
      write(6,*) 'qcal electrons = ',qcal, electrons
c  rescale density
      rescale=electrons/qcal
      !call dscal(n,rescale,den,1)
      den=rescale*den
c
c
       loop=-1
       cnvrg=small*(nz**2)
       write(6,*) 'Density convergence parameter set to ',cnvrg
       rewind iscr1
c
c    start iteration loop
c
 2001 continue
      loop=loop+1
c  calculate  potential * r == rv for given density
c    
      do i=1,n
       totden(i)=den(i)+coreden(i)
      enddo
      etot=electrons+qcore
      call potential(n,h,nz,etot,totden,rv,ecoul,etxc,eexc,v0,v0p,
     x       w,w(1,2),w(1,3),w(1,4),itype)
c
c  solve for bound states of Schroedinger equation
c
      icount=0
      qf=nz-electrons-qcore
      jierr=0
      it=0
c  s states :
      if (nps.gt.0) then
      it=it+1
      emin=-nz*nz
      l=0
      nroot=nps
      if (nroot+4.gt.mxwrk) then
       write(6,*) 'error in boundsch --increase mxwrk gt',nroot+4
       stop
       endif
      call boundsch(n,e,nmx,w,l,nz,nroot,rv,h,v0,v0p,qf,emin,
     x eig(icount+1),niter,ierr,ntroot,w(1,nroot+1),w(1,nroot+2),
     x w(1,nroot+3),w(1,nroot+4))
      no(l+1)=ntroot
      ii=ierr
      do io=1,nroot
      icount=icount+1
      jj=mod(ii,10)
      if (iscore(icount).eq.1) then
      if (jj.eq.1.or.jj.eq.2) then
      eig(icount)=e(io)
      !call dcopy (n,w(1,io),1,wfn(1,icount),1)
      wfn(:,icount)=w(:,io)
      else
      jierr=jierr+1
      endif
      endif
      ii=ii/10
c     write(6,*) 'io,ierr',io,ierr,jj,e(io)
      enddo
      endif
c  p states :
      if (npp.gt.1) then
      it=it+1
      emin=-nz*nz/4.d0
      l=1
      nroot=npp-1
      if (nroot+4.gt.mxwrk) then
       write(6,*) 'error in boundsch --increase mxwrk gt',nroot+4
       stop
       endif
      call boundsch(n,e,nmx,w,l,nz,nroot,rv,h,v0,v0p,qf,emin,
     x eig(icount+1),niter,ierr,ntroot,w(1,nroot+1),w(1,nroot+2),
     x w(1,nroot+3),w(1,nroot+4))
      no(l+1)=ntroot
      ii=ierr
      do io=1,nroot
      icount=icount+1
      jj=mod(ii,10)
      if (iscore(icount).eq.1) then
      if (jj.eq.1.or.jj.eq.2) then
      eig(icount)=e(io)
      !call dcopy (n,w(1,io),1,wfn(1,icount),1)
      wfn(:,icount)=w(:,io)
      else
      jierr=jierr+1
      endif
      endif
      ii=ii/10
c     write(6,*) 'io,ierr',io,ierr,jj,e(io)
      enddo
      endif
c  d states :
      if (npd.gt.2) then
      it=it+1
      emin=-nz*nz/9.d0
      l=2
      nroot=npd-2
      if (nroot+4.gt.mxwrk) then
       write(6,*) 'error in boundsch --increase mxwrk gt',nroot+4
       stop
       endif
      call boundsch(n,e,nmx,w,l,nz,nroot,rv,h,v0,v0p,qf,emin,
     x eig(icount+1),niter,ierr,ntroot,w(1,nroot+1),w(1,nroot+2),
     x w(1,nroot+3),w(1,nroot+4))
      no(l+1)=ntroot
      ii=ierr
      do io=1,nroot
      icount=icount+1
      jj=mod(ii,10)
      if (iscore(icount).eq.1) then
      if (jj.eq.1.or.jj.eq.2) then
      eig(icount)=e(io)
      !call dcopy (n,w(1,io),1,wfn(1,icount),1)
      wfn(:,icount)=w(:,io)
      else
      jierr=jierr+1
      endif
      endif
      ii=ii/10
c     write(6,*) 'io,ierr',io,ierr,jj,e(io)
      enddo
      endif
c  f states :
      if (npf.gt.3) then
      it=it+1
      emin=-nz*nz/16.d0
      l=3
      nroot=npf-3
      if (nroot+4.gt.mxwrk) then
       write(6,*) 'error in boundsch --increase mxwrk gt',nroot+4
       stop
       endif
      call boundsch(n,e,nmx,w,l,nz,nroot,rv,h,v0,v0p,qf,emin,
     x eig(icount+1),niter,ierr,ntroot,w(1,nroot+1),w(1,nroot+2),
     x w(1,nroot+3),w(1,nroot+4))
      no(l+1)=ntroot
      ii=ierr
      do io=1,nroot
      icount=icount+1
      jj=mod(ii,10)
      if (iscore(icount).eq.1) then
      if (jj.eq.1.or.jj.eq.2) then
      eig(icount)=e(io)
      !call dcopy (n,w(1,io),1,wfn(1,icount),1)
      wfn(:,icount)=w(:,io)
      else
      jierr=jierr+1
      endif
      endif
      ii=ii/10
c     write(6,*) 'io,ierr',io,ierr,jj,e(io)
      enddo
      endif
c  g states :
      if (npg.gt.4) then
      it=it+1
      emin=-nz*nz/25.d0
      l=4
      nroot=npg-4
      if (nroot+4.gt.mxwrk) then
       write(6,*) 'error in boundsch --increase mxwrk gt',nroot+4
       stop
       endif
      call boundsch(n,e,nmx,w,l,nz,nroot,rv,h,v0,v0p,qf,emin,
     x eig(icount+1),niter,ierr,ntroot,w(1,nroot+1),w(1,nroot+2),
     x w(1,nroot+3),w(1,nroot+4))
      no(l+1)=ntroot
      ii=ierr
      do io=1,nroot
      icount=icount+1
      jj=mod(ii,10)
      if (iscore(icount).eq.1) then
      if (jj.eq.1.or.jj.eq.2) then
      eig(icount)=e(io)
      !call dcopy (n,w(1,io),1,wfn(1,icount),1)
      wfn(:,icount)=w(:,io)
      else
      jierr=jierr+1
      endif
      endif
      ii=ii/10
c     write(6,*) 'io,ierr',io,ierr,jj,e(io)
      enddo
      endif
c
c  calculate new density
c
      !call dcopy(n,0.d0,0,denout,1)
      denout=0.d0
      write(6,*) '  results for loop = ',loop
      write(6,*) ' n  l     occupancy       energy'
      do io=1,norbit
      if (iscore(io).eq.1) then
      write(6,630) nlst(io),llst(io),occ(io),eig(io)
  630 format(i2,1x,i2,4x,1p2e15.7)
      if (occ(io).gt.small) then
        do i=1,n
        denout(i)=denout(i)+occ(io)*(wfn(i,io)**2)
        enddo
      endif
      endif
      enddo
      qcal=overint(n,h,denout)
      write(6,*) 'qcal electrons = ',qcal, electrons
c  rescale density
      rescale=electrons/qcal
      !call dscal(n,rescale,denout,1)
      denout=rescale*denout
c
c  call broyden iteration
c
      write(iscr1) (den(i),i=1,n)
      write(iscr1) (denout(i),i=1,n)
      call rhoiter(loop,n,den,denout,rmix,w,w(1,2),w(1,3),w(1,4),
     x w(1,5),w(1,6),iscr1,iscr2,mxloop,delta)
c  correct and recale density 
       do i=1,n
       if (den(i).lt.0.d0) den(i)=0.d0
       enddo
      qcal=overint(n,h,den)
      write(6,*) 'qcal electrons = ',qcal, electrons
c  rescale density
      rescale=electrons/qcal
      !call dscal(n,rescale,den,1)
      den=rescale*den
c
      write(6,*) '  results for loop ,delta = ',loop,delta
      write(6,*) ' n  l     occupancy       energy'
      eone=0.d0
      do io=1,norbit
      write(6,630) nlst(io),llst(io),occ(io),eig(io)
      eone=eone+occ(io)*eig(io)
      enddo
      write(6,*)
      write(6,*) ' Total energies'
      write(6,*) '    One-electron contribution:  ',eone
      write(6,*) '    Coulomb contribution     :  ',ecoul
      write(6,*) '    Exch-correl contribution :  ',etxc
      etot=eone-ecoul+etxc
      write(6,*) '    Total                    :  ',etot
c
      if (loop.lt.4) go to 2001
       if (delta.gt.cnvrg.and.loop.gt.mxloop) then
       write(6,*)'calculation terminating without density convergence'
       write(6,*) 'delta, cnvrg =', delta,cnvrg
       write(6,*) 'loop,mxloop = ',loop,mxloop
       stop
       endif
      if (delta.gt.cnvrg) go to 2001
c
c  converged result
c
      write(6,*) '  calculation converged in',loop,' iterations'
      write(6,*) '     for nz = ',nz
      write(6,*) '    delta(density)  = ', delta
      write(6,*) '  results for loop = ',loop
      write(6,*) ' n  l     occupancy       energy'
      eone=0.d0
      do io=1,norbit
      write(6,630) nlst(io),llst(io),occ(io),eig(io)
      eone=eone+occ(io)*eig(io)
      enddo
      write(6,*)
      write(6,*) ' Total energies'
      write(6,*) '    One-electron contribution:  ',eone
      write(6,*) '    Coulomb contribution     :  ',ecoul
      write(6,*) '    Exch-correl contribution :  ',etxc
      etot=eone-ecoul+etxc
      write(6,*) '    Total                    :  ',etot
      return
      end