program atom
c  program to calculate the self-consistent density functional
c    atom ground state fo atom with atomic number nz
      implicit real*8 (a-h,o-z)
      parameter (mxnl=10,mxtl=mxnl*mxnl,nmx=20000,mxo=30,mxloop=200)
      parameter (iscr1=30,iscr2=31)
      parameter (ifden=7,ifwfn=8,ifsum=9)
      dimension occ(mxo),wfn(nmx,mxo),eig(mxo),w(nmx,mxo),e(mxo)
      dimension r(nmx),den(nmx),rv(nmx),denout(nmx)
      dimension nl(mxnl,mxnl),nlst(mxo),llst(mxo),no(mxnl)
      character*4 flnm
      character*20 nm
      character*2 atomname
      data nl/mxtl*0/
      open(unit=iscr1,status='scratch',form='unformatted')
      open(unit=iscr2,status='scratch',form='unformatted')
      write(6,*) 'Enter atomic symbol'
      read(5,'(a2)') atomname
      open (unit=ifsum,file=TRIM(atomname),form='formatted')
      call dfatom(nz,nps,npp,npd,npf,npg,n,h,mxnl,nl,mxo,occ,norbit,
     x nlst,llst,no,qf,nmx,r,wfn,eig,den,denout,rv,itype,mxloop,iscr1,
     x iscr2,mxo,w,e,small,ecoul,etxc,etot)
c
 1111 write(6,*) 'enter 0 if wish to recalculate atomic configuration'
      read(5,*) itis
      if (itis.eq.0) then
       write(6,*) 'for each of the following states enter new occ'
       do io=1,norbit
       write(6,600) io,nlst(io),llst(io),occ(io),eig(io)
  600  format(3i5,1p2e15.7)
       read(5,*) occ(io)
       enddo
      call rcdfatom(nz,nps,npp,npd,npf,npg,n,h,mxnl,nl,mxo,occ,norbit,
     x nlst,llst,no,qf,nmx,r,wfn,eig,den,denout,rv,itype,mxloop,iscr1,
     x iscr2,mxo,w,e,small,ecoul,etxc,etot)
       goto 1111
       endif
c
c write density and wavefunctions
c
       open(unit=ifden,file='density',form='formatted')
       do i=1,n
        if (r(i).le.6.d0) then
         write(ifden,610) r(i),den(i)
  610    format(1p6e12.4)
         endif
       enddo
         close (ifden)
         open(unit=ifden,file='plotdensity',form='formatted')
         nm='density'
         write(ifden,620) atomname,nm
         close (ifden)
         do i=1,n
           do io=1,norbit
           if(dabs(wfn(i,io)).lt.1.d-8) wfn(i,io)=0.d0
           enddo
         enddo
       istart=0
       do l=0,4
        if (l.eq.0) many=nps
        if (l.eq.1) many=npp-1
        if (l.eq.2) many=npd-2
        if (l.eq.3) many=npf-3
        if (l.eq.4) many=npg-4
        if (many.gt.0) then
         call mkname(l,flnm)
         open(unit=ifwfn,file='wfn'//flnm,form='formatted')
         do  i=1,n
          if (r(i).le.6.d0) then
           write(ifwfn,'(1p8e10.2)') r(i),(wfn(i,j),j=istart+1,istart+many)
           endif
           enddo
         close (ifwfn)
         nm='wfn'//flnm
         open(unit=ifwfn,file='plotwfn'//flnm,form='formatted')
         if (many.eq.1) write(ifwfn,620) atomname,nm
         if (many.eq.2) write(ifwfn,621) atomname,nm,nm
         if (many.eq.3) write(ifwfn,622) atomname,nm,nm,nm
         if (many.eq.4) write(ifwfn,623) atomname,nm,nm,nm,nm
         if (many.eq.5) write(ifwfn,624) atomname,nm,nm,nm,nm,nm
         if (many.eq.6) write(ifwfn,625) atomname,nm,nm,nm,nm,nm,nm
         if (many.eq.7) write(ifwfn,626) atomname,nm,nm,nm,nm,nm,nm,nm
  620    format('gplot -t ',a2,' -tx "r (bohr)" -f ',a8,' 1 2 lines')
  621    format('gplot -t ',a2,' -tx "r (bohr)" -f ',
     x     a8,' 1 2 lines -f ',a8,' 1 3 lines')
  622    format('gplot -t ',a2,' -tx "r (bohr)" -f ',
     x     a8,' 1 2 lines -f ',a8,' 1 3 lines -f ',
     x     a8,' 1 4 lines')
  623    format('gplot -t ',a2,' -tx "r (bohr)" -f ',
     x     a8,' 1 2 lines -f ',a8,' 1 3 lines -f ',
     x     a8,' 1 4 lines -f ',a8,' 1 5 lines')
  624    format('gplot -t ',a2,' -tx "r (bohr)" -f ',
     x     a8,' 1 2 lines -f ',a8,' 1 3 lines -f ',
     x     a8,' 1 4 lines -f ',a8,' 1 5 lines -f ',a8,' 1 6 lines')
  625    format('gplot -t ',a2,' -tx "r (bohr)" -f ',
     x     a8,' 1 2 lines -f ',a8,' 1 3 lines -f ',
     x     a8,' 1 4 lines -f ',a8,' 1 5 lines -f ',a8,' 1 6 lines -f ',
     x     a8,' 1 7 lines')
  626    format('gplot -t ',a2,' -tx "r (bohr)" -f ',
     x     a8,' 1 2 lines -f ',a8,' 1 3 lines -f ',
     x     a8,' 1 4 lines -f ',a8,' 1 5 lines -f ',a8,' 1 6 lines -f ',
     x     a8,' 1 7 lines -f ', a8,' 1 8 lines')
         close (ifwfn)
       endif
       istart=istart+many
       enddo

c write summary information
       write(ifsum,*) '****** Results for atom ',atomname,' ********'
       write(ifsum,*) ' n  l     occupancy       energy'
       eone=0.d0
       do io=1,norbit
       write(ifsum,630) nlst(io),llst(io),occ(io),eig(io)
  630  format(i2,1x,i2,4x,1p2e15.7)
       eone=eone+occ(io)*eig(io)
       enddo
       write(ifsum,*)
       write(ifsum,*) ' Total energies'
       write(ifsum,*) '    One-electron contribution:  ',eone
       write(ifsum,*) '    Coulomb contribution     :  ',ecoul
       write(ifsum,*) '    Exch-correl contribution :  ',etxc
       etot=eone-ecoul+etxc
       write(ifsum,*) '    Total                    :  ',etot
       
      stop
      end