subroutine dfatom(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 iscr1,iscr2,mxwrk,w,e,small)
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
      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)
      dimension w(nmx,mxwrk),e(mxwrk)
      dimension nl(mxnl,mxnl),nlst(mxo),llst(mxo),no(mxnl)
      if (n.gt.nmx) then
       write(6,*) 'error in dfatom -- n,nmx = ',n,nmx
       stop
       endif
      rmix=rimix
      write(6,*) 'enter atomic number'
      read(5,*) nz
      write(6,*) 'Calculation for atomic number = ',nz
      write(6,*) 'enter maximum principle quantum numbers for s,p,d,f,g'
      read(5,*) nps,npp,npd,npf,npg
      write(6,600) nps,npp,npd,npf,npg
  600 format('ns np nd nf ng = ',5i4)
      icount=0
       if (nps.gt.0) then
       do is=1,nps
       icount=icount+1
       nl(is,1)=icount
       occ(icount)=2.d0
       nlst(icount)=is
       llst(icount)=0
       enddo
       endif
       if (npp.gt.1) then
       do ip=2,npp
       icount=icount+1
       nl(ip,2)=icount
       occ(icount)=6.d0
       nlst(icount)=ip
       llst(icount)=1
       enddo
       endif
       if (npd.gt.2) then
       do id=3,npd
       icount=icount+1
       nl(id,3)=icount
       occ(icount)=10.d0
       nlst(icount)=id
       llst(icount)=2
       enddo
       endif
       if (npf.gt.3) then
       do jf=4,npf
       icount=icount+1
       nl(jf,4)=icount
       occ(icount)=14.d0
       nlst(icount)=jf
       llst(icount)=3
       enddo
       endif
       if(npg.gt.4) then
       do ig=5,npg
       icount=icount+1
       nl(ig,5)=icount
       occ(icount)=18.d0
       nlst(icount)=ig
       llst(icount)=4
       enddo
       endif
       norbit=icount
       write(6,*) norbit, ' orbitals will be calculated'
       if (norbit.gt.mxo) then
        write(6,*) 'error in dfatom -- norbit,mxo=',norbit,mxo
        stop
        endif
c
       write(6,*)' Below are list the default occupations '
       write(6,610) 
  610 format(' n  l     occupancy')
       do io=1,norbit      
       write(6,620) nlst(io),llst(io),occ(io)
  620 format(i2,1x,i2,4x,1pe15.7)
       enddo
c
       write(6,*)' enter np l occ for all occupations for all revisions'
       write(6,*)'  enter 0 0 0. to end'
 1001 continue
      read(5,*) ip,l,xocc
      if (ip.le.0) go to 1011
      nfix=nl(ip,l+1)
      if (nfix.le.0.or.nfix.gt.norbit) then
       write(6,*) 'error in occupations -- ip,l,xocc',
     x     ip,l,xocc,nfix,norbit
       stop
      endif
      occ(nfix)=xocc
      go to 1001
 1011 continue
c
      write(6,*) ' Corrected occupations are: '
      write(6,610)
       electrons=0.d0
       do io=1,norbit      
       write(6,620) nlst(io),llst(io),occ(io)
       electrons=electrons+occ(io)
       enddo
      qf=nz-electrons
      write(6,*)
      write(6,*) 'nuclear charge    = ' , nz
      write(6,*) 'electronic charge = ', electrons
      write(6,*) 'net charge        = ', qf
c
c  read in # of mesh points and range of calculation
c
      n=5000
      range=50.d0
      small=1.d-10
      write(6,*) ' Calculational parameters used:'
      write(6,*) ' Number of mesh points:   ',n
      write(6,*) ' Integration range:       ',range
      h=range/n
      write(6,*) ' Integration step:        ',h
      write(6,*) ' Error tolerance:         ',small
      !read(5,*) n,range,small
      !write(6,*) 'n range h = ', n,range,h
      do ir=1,n
      r(ir)=ir*h
      enddo
      if (small.lt.0.d0.or.small.gt.1.d0) small=small0
      
c
c  calculate initial charge density from hydrogen-like functions
c    also initial energies
c
      !call dcopy(n,0.d0,0,den,1)
      den=0.d0
      zeff=nz
      do io=1,norbit
      np=nlst(io)
      l=llst(io)
      xocc=occ(io)
      zeff=zeff-0.5d0*xocc
      zeff=dmax1(zeff,1.d0)
      nzeff=zeff+0.1
      eig(io)=-(zeff/(np))**2
      write(6,*) io,np,l,xocc,eig(io)
      do ir=1,n
      wfn(ir,io)=hwfn(nzeff,np,l,r(ir))
      den(ir)=den(ir)+xocc*(wfn(ir,io)**2)
      enddo
      zeff=zeff-0.5*xocc
      enddo
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  choose form of exchange-correlation potential
      write(6,*) 'input exchange-correlation type -- 1 < itype < 15'
      !read(5,*) itype
      itype=14
      if (itype.eq.13) write(6,*) 'Wigner correlation'
      if (itype.eq.14) write(6,*) 'Perdew-Wang correlation'
      if (itype.lt.13) write(6,*) 'Ceperley-WVN correlation'
      if (itype.eq.15) write(6,*) 'Perdew''s gga-II ex & corr'
      if (itype.gt.15) then
       write(6,*) 'itype gt 15 -- resetting to Perdew-Wang form'
       itype=14
       endif
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    
      call potential(n,h,nz,electrons,den,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
      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 (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
      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 (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
      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 (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
      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 (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
      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 (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
      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
      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
      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

      function factorial(n)
      implicit real*8 (a-h,o-z)
      factorial=1
      if (n.lt.2) return
      do i=2,n
      factorial=factorial*i
      enddo
      return
      end
      function hwfn(nz,np,l,r)
c function to calculate the radial H wfn for nuclear charge nl
c                                            principal qn   np
c                                            orbital qn     l
c   r*(radial H wfn) is returned
      implicit real*8 (a-h,o-z)
      node=np-l-1
      scale=2.d0*nz/np
      rho=scale*r
      pref=scale*dsqrt(scale*factorial(np+l)/(2*np*factorial(node)))
      term=(rho**l)/factorial(2*l+1)
      sum=term
      if (node.gt.0) then
      do k=1,node
      term=-term*(node-k+1)*rho/(k*(2*l+1+k))
      sum=sum+term
      enddo
      endif
      hwfn=r*pref*dexp(-0.5d0*rho)*sum
c     write(6,*) 'r,hwfn=',r,hwfn
      return
      end