PROGRAM graphatom  
  !*************************************************************
  !  program to calculate the self-consistent density functional
  !    atom ground state for atom with atomic number nz
  !************************************************************
  USE GlobalMath
  USE atomdata
  USE aeatom
  IMPLICIT NONE
  INTEGER, PARAMETER :: iscr1=30,iscr2=31, ifden=7,ifwfn=8  
  CHARACTER (len=4) :: flnm  
  CHARACTER (len=20) :: nm  
  CHARACTER (len=2) :: sym  
  CHARACTER (len=1) :: syml  
  REAL(8), POINTER :: r(:),den(:),rv(:),wfn(:,:)
  INTEGER, POINTER :: n,norbit,nps,npp,npd,npf,npg

  INTEGER :: i,j,io,many,l,istart
  TYPE (GridInfo), TARGET :: AEGrid
  TYPE (PotentialInfo), TARGET :: AEPot
  TYPE (OrbitInfo), TARGET :: AEOrbit
  TYPE (SCFInfo), TARGET :: AESCF

  CALL Init_GlobalConstants
  CALL SCFatom(AEGrid,AEPot,AEOrbit,AESCF)
  n=>AEGrid%n
  r=>AEGrid%r
  den=>AEPot%den
  rv=>AEPot%rv
  norbit=>AEOrbit%norbit
  nps=>AEOrbit%nps
  npp=>AEOrbit%npp
  npd=>AEOrbit%npd
  npf=>AEOrbit%npf
  npg=>AEOrbit%npg
  wfn=>AEOrbit%wfn

  !
  ! write density and wavefunctions
  !

  OPEN (unit = ifden, file=TRIM(AEPot%sym), form='formatted')
  WRITE(ifden,*) 'Completed calculations for ',TRIM(AEPOT%sym)
  WRITE(ifden,*) '  aeatom converged in',AESCF%iter,' iterations'
  WRITE(ifden,*) '     for nz = ',AEPot%nz
  WRITE(ifden,*) '    delta(density)  = ', AESCF%delta
  WRITE(ifden,*) '  Orbital energies:         '
  WRITE(ifden,*) ' n  l     occupancy       energy'
  DO io=1,AEorbit%norbit
     WRITE(ifden,'(i2,1x,i2,4x,1p2e15.7)') &
          AEOrbit%np(io),AEOrbit%l(io),&
          AEOrbit%occ(io),AEOrbit%eig(io)
  ENDDO
  WRITE(ifden,*)
  WRITE(ifden,*) ' Total energies'
  WRITE(ifden,*) '    One-electron contribution:  ',AESCF%eone
  WRITE(ifden,*) '    Coulomb contribution     :  ',AESCF%ecoul
  WRITE(ifden,*) '    Exch-correl contribution :  ',AESCF%eexc
  WRITE(ifden,*) '    Total                    :  ',AESCF%etot
  CLOSE(ifden)

  OPEN (unit = ifden, file = 'density', form = 'formatted')  
  DO i = 1, n  
     IF (r (i) .LE.6.d0) THEN  
        WRITE (ifden, '(1p6e12.4)') r (i), den (i)  
     ENDIF
  ENDDO
  CLOSE (ifden)  

  OPEN (unit = ifden, file = 'potential', form = 'formatted')  
  DO i = 1, n  
     IF (r (i) .LE.6.d0) THEN  
        WRITE (ifden, '(1p6e12.4)') r (i), rv (i)  
     ENDIF
  ENDDO
  CLOSE (ifden)  

  OPEN (unit = ifden, file = 'plotdensity', form = 'formatted')  
  nm = 'density'  
  WRITE (ifden, '("gplot -t ""Radial density for ",a2, &
       &""" -tx ""r (bohr)"" -f ",a8, " 1 2 lines")') AEpot%sym, nm
  CLOSE (ifden)  
  OPEN (unit = ifden, file = 'plotpotential', form = 'formatted')  
  nm = 'potential'  
  WRITE (ifden, '("gplot -t ""rxV(r) for ",a2, &
       &""" -tx ""r (bohr)"" -f ",a9, " 1 2 lines")') AEpot%sym, 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 (l.EQ.0) syml = 's'  
     IF (l.EQ.1) syml = 'p'  
     IF (l.EQ.2) syml = 'd'  
     IF (l.EQ.3) syml = 'f'  
     IF (l.EQ.4) syml = 'g'  
     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, '("gplot -t ""Radial ",a1, &
             &   " wavefunctions for ",a2, &
             &""" -tx ""r (bohr)"" -f ",a8, " 1 2 lines")') syml, AEpot%sym, nm
        IF (many.EQ.2) WRITE (ifwfn, '("gplot -t ""Radial ",a1, &
             &   " wavefunctions for ",a2, &
             &""" -tx ""r (bohr)"" -f ",a8, " 1 2 lines -f ",a8," 1 3 lines")') syml &
             , AEpot%sym, nm, nm
        IF (many.EQ.3) WRITE (ifwfn, '("gplot -t ""Radial ",a1, &
             &   " wavefunctions for ",a2, &
             &""" -tx ""r (bohr)"" -f ",a8, " 1 2 lines -f ",a8," 1 3 lines -f ", &
             &a8," 1 4 lines")') syml, AEpot%sym, nm, nm, nm
        IF (many.EQ.4) WRITE (ifwfn, '("gplot -t ""Radial ",a1, &
             &   " wavefunctions for ",a2, &
             &""" -tx ""r (bohr)"" -f ",a8, " 1 2 lines -f ",a8," 1 3 lines -f ", &
             &a8," 1 4 lines -f ",a8," 1 5 lines")') syml, AEpot%sym, nm, nm, nm, nm
        IF (many.EQ.5) WRITE (ifwfn, '("gplot -t ""Radial ",a1, &
             &   " wavefunctions for ",a2, &
             &""" -tx ""r (bohr)"" -f ",a8, " 1 2 lines -f ",a8," 1 3 lines -f ", &
             &a8," 1 4 lines -f ",a8," 1 5 lines -f ",a8," 1 6 lines")') syml,  &
             &AEpot%sym, nm, nm, nm, nm, nm
        IF (many.EQ.6) WRITE (ifwfn, '("gplot -t ""Radial ",a1, &
             &   " wavefunctions for ",a2, &
             &""" -tx ""r (bohr)"" -f ",a8, " 1 2 lines -f ",a8," 1 3 lines -f ", &
             &a8," 1 4 lines -f ",a8," 1 5 lines -f ",a8," 1 6 lines -f " &
             &,a8," 1 7 lines")') syml, AEpot%sym, nm, nm, nm, nm, nm, nm
        IF (many.EQ.7) WRITE (ifwfn, '("gplot -t ""Radial ",a1, &
             &   " wavefunctions for ",a2, &
             &""" -tx ""r (bohr)"" -f ",a8, " 1 2 lines -f ",a8," 1 3 lines -f ", &
             &a8," 1 4 lines -f ",a8," 1 5 lines -f ",a8," 1 6 lines -f " &
             &,a8," 1 7 lines -f ",a8," 1 8 lines" )') syml, AEpot%sym, nm, nm, nm, nm &
             &, nm, nm, nm
        CLOSE (ifwfn)  
     ENDIF
     istart = istart + many  

  ENDDO
END PROGRAM graphatom