atom_output.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Routines that print various information about an atomic kind.
! **************************************************************************************************
MODULE atom_output
   USE atom_types,                      ONLY: &
        atom_basis_type, atom_gthpot_type, atom_potential_type, atom_state, atom_type, cgto_basis, &
        ecp_pseudo, gth_pseudo, gto_basis, lmat, no_pseudo, num_basis, sgp_pseudo, sto_basis, &
        upf_pseudo
   USE atom_utils,                      ONLY: get_maxl_occ,&
                                              get_maxn_occ,&
                                              get_rho0
   USE cp_files,                        ONLY: close_file,&
                                              open_file
   USE input_constants,                 ONLY: &
        barrier_conf, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
        do_rhf_atom, do_rks_atom, do_rohf_atom, do_sczoramp_atom, do_uhf_atom, do_uks_atom, &
        do_zoramp_atom, poly_conf, xc_none
   USE input_cp2k_check,                ONLY: xc_functionals_expand
   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                              section_vals_get_subs_vals2,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: default_string_length,&
                                              dp
   USE mathconstants,                   ONLY: dfac,&
                                              pi,&
                                              rootpi
   USE periodic_table,                  ONLY: ptable
   USE physcon,                         ONLY: evolt
   USE xc_derivatives,                  ONLY: xc_functional_get_info
   USE xc_libxc,                        ONLY: libxc_check_existence_in_libxc,&
                                              libxc_get_reference_length
   USE xmgrace,                         ONLY: xm_graph_data,&
                                              xm_graph_info,&
                                              xm_write_defaults,&
                                              xm_write_frame,&
                                              xm_write_frameport
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_output'
 
   PUBLIC :: atom_print_state, atom_print_energies, atom_print_iteration, &
             atom_print_basis, atom_print_method, atom_print_info, atom_print_potential, &
             atom_print_basis_file, atom_write_pseudo_param, atom_print_orbitals, &
             atom_print_zmp_iteration
 
CONTAINS
 
! **************************************************************************************************
!> \brief Print an information string related to the atomic kind.
!> \param zval  atomic number
!> \param info  information string
!> \param iw    output file unit
!> \par History
!>    * 09.2008 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_info(zval, info, iw)
      INTEGER, INTENT(IN)                                :: zval
      CHARACTER(len=*), INTENT(IN)                       :: info
      INTEGER, INTENT(IN)                                :: iw
 
      WRITE (iw, '(/," ",A,T40,A," [",A,"]",T62,"Atomic number:",T78,I3,/)') &
         adjustl(trim(info)), trim(ptable(zval)%name), trim(ptable(zval)%symbol), zval
 
   END SUBROUTINE atom_print_info
 
! **************************************************************************************************
!> \brief Print information about electronic state.
!> \param state  electronic state
!> \param iw     output file unit
!> \par History
!>    * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
!>    * 11.2009 print multiplicity [Juerg Hutter]
!>    * 08.2008 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_state(state, iw)
      TYPE(atom_state)                                   :: state
      INTEGER, INTENT(IN)                                :: iw
 
      CHARACTER(LEN=1), DIMENSION(0:7), PARAMETER :: &
         label = (/"S", "P", "D", "F", "G", "H", "I", "K"/)
 
      INTEGER                                            :: j, l, mc, mlc, mlo, mm(0:lmat), mo
 
      cpassert(lmat <= 7)
      WRITE (iw, '(/,T2,A)') "Electronic structure"
      WRITE (iw, '(T5,A,T71,F10.2)') "Total number of core electrons", sum(state%core)
      WRITE (iw, '(T5,A,T71,F10.2)') "Total number of valence electrons", sum(state%occ)
      WRITE (iw, '(T5,A,T71,F10.2)') "Total number of electrons", sum(state%occ + state%core)
      SELECT CASE (state%multiplicity)
      CASE (-1)
         WRITE (iw, '(T5,A,T68,A)') "Multiplicity", "not specified"
      CASE (-2)
         WRITE (iw, '(T5,A,T72,A)') "Multiplicity", "high spin"
      CASE (-3)
         WRITE (iw, '(T5,A,T73,A)') "Multiplicity", "low spin"
      CASE (1)
         WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "singlet"
      CASE (2)
         WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "doublet"
      CASE (3)
         WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "triplet"
      CASE (4)
         WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "quartet"
      CASE (5)
         WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "quintet"
      CASE (6)
         WRITE (iw, '(T5,A,T75,A)') "Multiplicity", "sextet"
      CASE (7)
         WRITE (iw, '(T5,A,T75,A)') "Multiplicity", "septet"
      CASE DEFAULT
      END SELECT
 
      mlo = get_maxl_occ(state%occ)
      mlc = get_maxl_occ(state%core)
      mm = get_maxn_occ(state%core)
 
      IF (state%multiplicity == -1) THEN
         DO l = 0, max(mlo, mlc)
            mo = state%maxn_occ(l)
            IF (sum(state%core(l, :)) == 0) THEN
               WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occ(l, j), j=1, mo)
            ELSE
               mc = mm(l)
               cpassert(sum(state%occ(l, 1:mc)) == 0)
               WRITE (iw, advance="no", fmt='(A5,T9,A1,10F6.2)') label(l), "[", (state%core(l, j), j=1, mc)
               WRITE (iw, fmt='(A1,F5.2,10F6.2)') "]", (state%occ(l, j), j=mc + 1, mc + mo)
            END IF
         END DO
      ELSE
         WRITE (iw, '(T5,A)') "Alpha Electrons"
         DO l = 0, max(mlo, mlc)
            mo = state%maxn_occ(l)
            IF (sum(state%core(l, :)) == 0) THEN
               WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occa(l, j), j=1, mo)
            ELSE
               mc = mm(l)
               WRITE (iw, advance="no", fmt='(A5,T9,A1,10F6.2)') label(l), "[", (0.5_dp*state%core(l, j), j=1, mc)
               WRITE (iw, fmt='(A1,F5.2,10F6.2)') "]", (state%occa(l, j), j=1, mo)
            END IF
         END DO
         WRITE (iw, '(T5,A)') "Beta Electrons"
         DO l = 0, max(mlo, mlc)
            mo = state%maxn_occ(l)
            IF (sum(state%core(l, :)) == 0) THEN
               WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occb(l, j), j=1, mo)
            ELSE
               mc = mm(l)
               WRITE (iw, advance="no", fmt='(A5,T9,A1,10F6.2)') label(l), "[", (0.5_dp*state%core(l, j), j=1, mc)
               WRITE (iw, fmt='(A1,F5.2,10F6.2)') "]", (state%occb(l, j), j=1, mo)
            END IF
         END DO
      END IF
      WRITE (iw, *)
 
   END SUBROUTINE atom_print_state
 
! **************************************************************************************************
!> \brief Print energy components.
!> \param atom  information about the atomic kind
!> \param iw    output file unit
!> \par History
!>    * 05.2010 print virial coefficient [Juerg Hutter]
!>    * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
!>    * 09.2008 print orbital energies [Juerg Hutter]
!>    * 08.2008 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_energies(atom, iw)
      TYPE(atom_type)                                    :: atom
      INTEGER, INTENT(IN)                                :: iw
 
      INTEGER                                            :: i, l, n
      REAL(kind=dp)                                      :: drho
 
      WRITE (iw, '(/,A,T36,A,T61,F20.12)') " Energy components [Hartree]", &
         "    Total Energy ::", atom%energy%etot
      WRITE (iw, '(T36,A,T61,F20.12)') "     Band Energy ::", atom%energy%eband
      WRITE (iw, '(T36,A,T61,F20.12)') "  Kinetic Energy ::", atom%energy%ekin
      WRITE (iw, '(T36,A,T61,F20.12)') "Potential Energy ::", atom%energy%epot
      IF (atom%energy%ekin /= 0.0_dp) THEN
         WRITE (iw, '(T36,A,T61,F20.12)') "   Virial (-V/T) ::", -atom%energy%epot/atom%energy%ekin
      END IF
      WRITE (iw, '(T36,A,T61,F20.12)') "     Core Energy ::", atom%energy%ecore
      IF (atom%energy%exc /= 0._dp) &
         WRITE (iw, '(T36,A,T61,F20.12)') "       XC Energy ::", atom%energy%exc
      WRITE (iw, '(T36,A,T61,F20.12)') "  Coulomb Energy ::", atom%energy%ecoulomb
      IF (atom%energy%eexchange /= 0._dp) &
         WRITE (iw, '(T34,A,T61,F20.12)') "HF Exchange Energy ::", atom%energy%eexchange
      IF (atom%potential%ppot_type /= no_pseudo) THEN
         WRITE (iw, '(T20,A,T61,F20.12)') "    Total Pseudopotential Energy ::", atom%energy%epseudo
         WRITE (iw, '(T20,A,T61,F20.12)') "    Local Pseudopotential Energy ::", atom%energy%eploc
         IF (atom%energy%elsd /= 0._dp) &
            WRITE (iw, '(T20,A,T61,F20.12)') "     Local Spin-potential Energy ::", atom%energy%elsd
         WRITE (iw, '(T20,A,T61,F20.12)') " Nonlocal Pseudopotential Energy ::", atom%energy%epnl
      END IF
      IF (atom%potential%confinement) THEN
         WRITE (iw, '(T36,A,T61,F20.12)') "     Confinement ::", atom%energy%econfinement
      END IF
 
      IF (atom%state%multiplicity == -1) THEN
         WRITE (iw, '(/,A,T20,A,T30,A,T36,A,T49,A,T71,A,/)') " Orbital energies", &
            "State", "L", "Occupation", "Energy[a.u.]", "Energy[eV]"
         DO l = 0, atom%state%maxl_calc
            n = atom%state%maxn_calc(l)
            DO i = 1, n
               WRITE (iw, '(T23,I2,T30,I1,T36,F10.3,T46,F15.6,T66,F15.6)') &
                  i, l, atom%state%occupation(l, i), atom%orbitals%ener(i, l), atom%orbitals%ener(i, l)*evolt
            END DO
            IF (n > 0) WRITE (iw, *)
         END DO
      ELSE
         WRITE (iw, '(/,A,T20,A,T30,A,T36,A,T42,A,T55,A,T71,A,/)') " Orbital energies", &
            "State", "Spin", "L", "Occupation", "Energy[a.u.]", "Energy[eV]"
         DO l = 0, atom%state%maxl_calc
            n = atom%state%maxn_calc(l)
            DO i = 1, n
               WRITE (iw, '(T23,I2,T29,A,T36,I1,T42,F10.3,T52,F15.6,T68,F13.6)') &
                  i, "alpha", l, atom%state%occa(l, i), atom%orbitals%enera(i, l), atom%orbitals%enera(i, l)*evolt
            END DO
            DO i = 1, n
               WRITE (iw, '(T23,I2,T29,A,T36,I1,T42,F10.3,T52,F15.6,T68,F13.6)') &
                  i, " beta", l, atom%state%occb(l, i), atom%orbitals%enerb(i, l), atom%orbitals%enerb(i, l)*evolt
            END DO
            IF (n > 0) WRITE (iw, *)
         END DO
      END IF
 
      CALL get_rho0(atom, drho)
      WRITE (iw, '(/,A,T66,F15.6)') " Total Electron Density at R=0: ", drho
 
   END SUBROUTINE atom_print_energies
 
! **************************************************************************************************
!> \brief Printing of the atomic iterations when ZMP is active.
!> \param iter  current iteration number
!> \param deps  convergence
!> \param atom  intormation about the atomic kind
!> \param iw    output file unit
!> \author D. Varsano [daniele.varsano@nano.cnr.it]
! **************************************************************************************************
   SUBROUTINE atom_print_zmp_iteration(iter, deps, atom, iw)
      INTEGER, INTENT(IN)                                :: iter
      REAL(dp), INTENT(IN)                               :: deps
      TYPE(atom_type), INTENT(IN)                        :: atom
      INTEGER, INTENT(IN)                                :: iw
 
      IF (iter == 1) THEN
         WRITE (iw, '(/," ",79("*"),/,T33,"Integral",T48,"Integral",/,T3,A,T16,A,T33,A,T46,A,T69,A/," ",79("*"))') &
            "Iteration", "Convergence", "rho diff.", "rho*v_xc[au]", "Energy[au]"
      END IF
      WRITE (iw, '(T3,I9,T15,G13.6,T30,G13.6,T46,G13.6,T61,F20.12)') iter, deps, atom%rho_diff_integral, &
         atom%energy%exc, atom%energy%etot
 
   END SUBROUTINE atom_print_zmp_iteration
 
! **************************************************************************************************
!> \brief Print convergence information.
!> \param iter  current iteration number
!> \param deps  convergency
!> \param etot  total energy
!> \param iw    output file unit
!> \par History
!>    * 08.2008 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_iteration(iter, deps, etot, iw)
      INTEGER, INTENT(IN)                                :: iter
      REAL(dp), INTENT(IN)                               :: deps, etot
      INTEGER, INTENT(IN)                                :: iw
 
      IF (iter == 1) THEN
         WRITE (iw, '(/," ",79("*"),/,T19,A,T38,A,T70,A,/," ",79("*"))') &
            "Iteration", "Convergence", "Energy [au]"
      END IF
      WRITE (iw, '(T20,i8,T34,G14.6,T61,F20.12)') iter, deps, etot
 
   END SUBROUTINE atom_print_iteration
 
! **************************************************************************************************
!> \brief Print atomic basis set.
!> \param atom_basis  atomic basis set
!> \param iw          output file unit
!> \param title       header to print on top of the basis set
!> \par History
!>    * 09.2008 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_basis(atom_basis, iw, title)
      TYPE(atom_basis_type)                              :: atom_basis
      INTEGER, INTENT(IN)                                :: iw
      CHARACTER(len=*)                                   :: title
 
      INTEGER                                            :: i, j, l
 
      WRITE (iw, '(/,A)') trim(title)
      SELECT CASE (atom_basis%basis_type)
      CASE (gto_basis)
         IF (atom_basis%geometrical) THEN
            WRITE (iw, '(/," ",21("*"),A,22("*"))') " Geometrical Gaussian Type Orbitals "
            WRITE (iw, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", atom_basis%aval, &
               " Proportionality factor: ", atom_basis%cval
         ELSE
            WRITE (iw, '(/," ",21("*"),A,21("*"))') " Uncontracted Gaussian Type Orbitals "
         END IF
         DO l = 0, lmat
            IF (atom_basis%nbas(l) > 0) THEN
               SELECT CASE (l)
               CASE DEFAULT
                  WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
                     "X Exponents: ", (i, atom_basis%am(i, l), i=1, atom_basis%nbas(l))
               CASE (0)
                  WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
                     "s Exponents: ", (i, atom_basis%am(i, 0), i=1, atom_basis%nbas(0))
               CASE (1)
                  WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
                     "p Exponents: ", (i, atom_basis%am(i, 1), i=1, atom_basis%nbas(1))
               CASE (2)
                  WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
                     "d Exponents: ", (i, atom_basis%am(i, 2), i=1, atom_basis%nbas(2))
               CASE (3)
                  WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
                     "f Exponents: ", (i, atom_basis%am(i, 3), i=1, atom_basis%nbas(3))
               END SELECT
            END IF
         END DO
         WRITE (iw, '(" ",79("*"))')
      CASE (cgto_basis)
         WRITE (iw, '(/," ",22("*"),A,22("*"))') " Contracted Gaussian Type Orbitals "
         DO l = 0, lmat
            IF (atom_basis%nbas(l) > 0) THEN
               IF (l == 0) WRITE (iw, '(A)') " s Functions"
               IF (l == 1) WRITE (iw, '(A)') " p Functions"
               IF (l == 2) WRITE (iw, '(A)') " d Functions"
               IF (l == 3) WRITE (iw, '(A)') " f Functions"
               IF (l >= 3) WRITE (iw, '(A)') " x Functions"
               DO i = 1, atom_basis%nprim(l)
                  WRITE (iw, '(F15.6,5(T21,6F10.6,/))') &
                     atom_basis%am(i, l), (atom_basis%cm(i, j, l), j=1, atom_basis%nbas(l))
               END DO
            END IF
         END DO
         WRITE (iw, '(" ",79("*"))')
      CASE (sto_basis)
         WRITE (iw, '(/," ",28("*"),A,29("*"))') " Slater Type Orbitals "
         DO l = 0, lmat
            DO i = 1, atom_basis%nbas(l)
               SELECT CASE (l)
               CASE DEFAULT
                  WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, l), "X Exponent :", atom_basis%as(i, l)
               CASE (0)
                  WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 0), "S Exponent :", atom_basis%as(i, 0)
               CASE (1)
                  WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 1), "P Exponent :", atom_basis%as(i, 1)
               CASE (2)
                  WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 2), "D Exponent :", atom_basis%as(i, 2)
               CASE (3)
                  WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 3), "F Exponent :", atom_basis%as(i, 3)
               END SELECT
            END DO
         END DO
         WRITE (iw, '(" ",79("*"))')
      CASE (num_basis)
         cpabort("")
      CASE DEFAULT
         cpabort("")
      END SELECT
 
   END SUBROUTINE atom_print_basis
 
! **************************************************************************************************
!> \brief Print the optimized atomic basis set into a file.
!> \param atom_basis  atomic basis set
!> \param wfn ...
!> \par History
!>    * 11.2016 revised output format [Matthias Krack]
!>    * 11.2011 Slater basis functions [Juerg Hutter]
!>    * 03.2011 created [Juerg Hutter]
!> \note  The basis set is stored as the file 'OPT_BASIS' inside the current working directory.
!>        It may be a good idea, however, to specify the name of this file via some input section.
! **************************************************************************************************
   SUBROUTINE atom_print_basis_file(atom_basis, wfn)
      TYPE(atom_basis_type)                              :: atom_basis
      REAL(kind=dp), DIMENSION(:, :, 0:), OPTIONAL       :: wfn
 
      INTEGER                                            :: i, im, iw, l
      REAL(kind=dp)                                      :: expzet, prefac, zeta
 
      CALL open_file(file_name="OPT_BASIS", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
      SELECT CASE (atom_basis%basis_type)
      CASE (gto_basis)
         IF (atom_basis%geometrical) THEN
            WRITE (iw, '(/," ",21("*"),A,22("*"))') " Geometrical Gaussian Type Orbitals "
            WRITE (iw, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", atom_basis%aval, &
               " Proportionality factor: ", atom_basis%cval
         ELSE
            WRITE (iw, '(T3,A)') "BASIS_TYPE GAUSSIAN"
         END IF
         DO l = 0, lmat
            IF (atom_basis%nbas(l) > 0) THEN
               SELECT CASE (l)
               CASE DEFAULT
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "X_EXPONENTS ", (atom_basis%am(i, l), i=1, atom_basis%nbas(l))
               CASE (0)
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "S_EXPONENTS ", (atom_basis%am(i, 0), i=1, atom_basis%nbas(0))
               CASE (1)
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "P_EXPONENTS ", (atom_basis%am(i, 1), i=1, atom_basis%nbas(1))
               CASE (2)
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "D_EXPONENTS ", (atom_basis%am(i, 2), i=1, atom_basis%nbas(2))
               CASE (3)
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "F_EXPONENTS ", (atom_basis%am(i, 3), i=1, atom_basis%nbas(3))
               END SELECT
            END IF
         END DO
      CASE (cgto_basis)
         cpabort("")
      CASE (sto_basis)
         WRITE (iw, '(T3,A)') "BASIS_TYPE SLATER"
         DO l = 0, lmat
            IF (atom_basis%nbas(l) > 0) THEN
               SELECT CASE (l)
               CASE DEFAULT
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "X_EXPONENTS ", (atom_basis%as(i, l), i=1, atom_basis%nbas(l))
                  WRITE (iw, '(T3,A,60I3)') &
                     "X_QUANTUM_NUMBERS ", (atom_basis%ns(i, l), i=1, atom_basis%nbas(l))
               CASE (0)
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "S_EXPONENTS ", (atom_basis%as(i, 0), i=1, atom_basis%nbas(0))
                  WRITE (iw, '(T3,A,60I3)') &
                     "S_QUANTUM_NUMBERS ", (atom_basis%ns(i, 0), i=1, atom_basis%nbas(0))
               CASE (1)
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "P_EXPONENTS ", (atom_basis%as(i, 1), i=1, atom_basis%nbas(1))
                  WRITE (iw, '(T3,A,60I3)') &
                     "P_QUANTUM_NUMBERS ", (atom_basis%ns(i, 1), i=1, atom_basis%nbas(1))
               CASE (2)
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "D_EXPONENTS ", (atom_basis%as(i, 2), i=1, atom_basis%nbas(2))
                  WRITE (iw, '(T3,A,60I3)') &
                     "D_QUANTUM_NUMBERS ", (atom_basis%ns(i, 2), i=1, atom_basis%nbas(2))
               CASE (3)
                  WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
                     "F_EXPONENTS ", (atom_basis%as(i, 3), i=1, atom_basis%nbas(3))
                  WRITE (iw, '(T3,A,60I3)') &
                     "F_QUANTUM_NUMBERS ", (atom_basis%ns(i, 3), i=1, atom_basis%nbas(3))
               END SELECT
            END IF
         END DO
      CASE (num_basis)
         cpabort("")
      CASE DEFAULT
         cpabort("")
      END SELECT
 
      IF (PRESENT(wfn)) THEN
         SELECT CASE (atom_basis%basis_type)
         CASE DEFAULT
         CASE (gto_basis)
            IF (.NOT. atom_basis%geometrical) THEN
               WRITE (iw, '(/,T3,A)') "ORBITAL COEFFICENTS (Quickstep normalization)"
               im = min(6, SIZE(wfn, 2))
               DO l = 0, lmat
                  IF (atom_basis%nbas(l) > 0) THEN
                     WRITE (iw, '(T3,A,I3)') "L Quantum Number:", l
                     ! Quickstep normalization
                     expzet = 0.25_dp*real(2*l + 3, dp)
                     prefac = sqrt(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
                     DO i = 1, atom_basis%nbas(l)
                        zeta = (2._dp*atom_basis%am(i, l))**expzet
                        WRITE (iw, '(T5,F14.8,2x,6F12.8)') atom_basis%am(i, l), wfn(i, 1:im, l)*prefac/zeta
                     END DO
                  END IF
               END DO
            END IF
         END SELECT
      END IF
 
      CALL close_file(unit_number=iw)
 
   END SUBROUTINE atom_print_basis_file
 
! **************************************************************************************************
!> \brief Print information about the electronic structure method in use.
!> \param atom  information about the atomic kind
!> \param iw    output file unit
!> \par History
!>    * 09.2015 direct use of the LibXC Fortran interface [Andreas Gloess]
!>    * 10.2012 LibXC interface [Fabien Tran]
!>    * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
!>    * 04.2009 print geometrical Gaussian type orbitals [Juerg Hutter]
!>    * 09.2008 new subroutine's prototype; print relativistic methods [Juerg Hutter]
!>    * 09.2008 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_method(atom, iw)
      TYPE(atom_type)                                    :: atom
      INTEGER, INTENT(IN)                                :: iw
 
      CHARACTER(len=160)                                 :: shortform
      CHARACTER(LEN=20)                                  :: tmpstr
      CHARACTER(len=:), ALLOCATABLE                      :: reference
      INTEGER                                            :: ifun, il, meth, myfun, reltyp
      LOGICAL                                            :: lsd
      TYPE(section_vals_type), POINTER                   :: xc_fun, xc_fun_section, xc_section
 
      NULLIFY (xc_fun, xc_fun_section, xc_section)
 
      meth = atom%method_type
 
      xc_section => atom%xc_section
      xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
      SELECT CASE (meth)
      CASE DEFAULT
         cpabort("")
      CASE (do_rks_atom)
         CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
      CASE (do_uks_atom)
         CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
      CASE (do_rhf_atom)
         myfun = xc_none
      CASE (do_uhf_atom)
         myfun = xc_none
      CASE (do_rohf_atom)
         myfun = xc_none
      END SELECT
 
      SELECT CASE (meth)
      CASE DEFAULT
         cpabort("")
      CASE (do_rks_atom)
         IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Restricted Kohn-Sham Calculation')")
      CASE (do_uks_atom)
         IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Unrestricted Kohn-Sham Calculation')")
      CASE (do_rhf_atom)
         IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Restricted Hartree-Fock Calculation')")
      CASE (do_uhf_atom)
         IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Unrestricted Hartree-Fock Calculation')")
      CASE (do_rohf_atom)
         IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Restricted Open-Shell Kohn-Sham Calculation')")
      END SELECT
 
      ! zmp
      IF (atom%do_zmp) THEN
         IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Method on atomic radial density')")
         IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Lambda : ',F5.1)") atom%lambda
         IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Reading external density : ',A20)") atom%ext_file
         IF (atom%dm) THEN
            IF (iw > 0) WRITE (iw, fmt="(' ZMP       | The file is in the form of a density matrix')")
         ELSE
            IF (iw > 0) WRITE (iw, fmt="(' ZMP       | The file is in the form of a linear density')")
         END IF
         IF (atom%doread) THEN
            IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Restarting calculation from ',A20,' file if present')") atom%zmp_restart_file
         END IF
      ELSE IF (atom%read_vxc) THEN
         IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Calculating density from external V_xc')")
         IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Reading external v_xc file : ',A20)") atom%ext_vxc_file
      END IF
 
      IF (atom%pp_calc) THEN
         IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Nonrelativistic Calculation')")
      ELSE
         reltyp = atom%relativistic
 
         SELECT CASE (reltyp)
         CASE DEFAULT
            cpabort("")
         CASE (do_nonrel_atom)
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Nonrelativistic Calculation')")
         CASE (do_zoramp_atom)
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using ZORA(MP)')")
         CASE (do_sczoramp_atom)
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using scaled ZORA(MP)')")
         CASE (do_dkh0_atom)
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Douglas-Kroll 0th order')")
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using kietic energy scaling')")
         CASE (do_dkh1_atom)
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Douglas-Kroll 1st order')")
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Foldy-Wouthuysen transformation')")
         CASE (do_dkh2_atom)
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Douglas-Kroll 2nd order')")
         CASE (do_dkh3_atom)
            IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Douglas-Kroll 3rd order')")
         END SELECT
      END IF
 
      lsd = (meth == do_uks_atom)
 
      IF (myfun /= xc_none) THEN
         CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", c_val=tmpstr)
         IF (iw > 0) WRITE (iw, fmt="(' FUNCTIONAL| ROUTINE=',a)") trim(tmpstr)
         CALL xc_functionals_expand(xc_fun_section, xc_section)
         IF (iw > 0) THEN
            ifun = 0
            DO
               ifun = ifun + 1
               xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
               IF (.NOT. ASSOCIATED(xc_fun)) EXIT
               IF (libxc_check_existence_in_libxc(xc_fun)) THEN
                  ALLOCATE (CHARACTER(LEN=libxc_get_reference_length(xc_fun, lsd)) :: reference)
               ELSE
                  ALLOCATE (CHARACTER(LEN=20*default_string_length) :: reference)
               END IF
               CALL xc_functional_get_info(xc_fun, lsd=lsd, reference=reference, shortform=shortform)
               WRITE (iw, fmt="(' FUNCTIONAL| ',a,':')") &
                  trim(xc_fun%section%name)
               DO il = 1, len_trim(reference), 67
                  WRITE (iw, fmt="(' FUNCTIONAL| ',a67)") reference(il:)
               END DO
               DEALLOCATE (reference)
            END DO
         END IF
      ELSE
         IF (iw > 0) WRITE (iw, fmt="(' FUNCTIONAL| NO EXCHANGE-CORRELATION FUNCTIONAL USED.')")
      END IF
 
   END SUBROUTINE atom_print_method
 
! **************************************************************************************************
!> \brief Print information about the pseudo-potential.
!> \param potential pseudo-potential
!> \param iw        output file unit
!> \par History
!>    * 05.2017 SGP pseudo-potentials [Juerg Hutter]
!>    * 02.2016 pseudo-potential in Quantum Espresso UPF format [Juerg Hutter]
!>    * 01.2016 new confinement potential form [Juerg Hutter]
!>    * 03.2010 extension of GTH pseudo-potential definition [Juerg Hutter]
!>    * 05.2009 GTH pseudo-potential [Juerg Hutter]
!>    * 09.2008 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_potential(potential, iw)
      TYPE(atom_potential_type)                          :: potential
      INTEGER, INTENT(IN)                                :: iw
 
      CHARACTER(len=60)                                  :: pline
      INTEGER                                            :: i, j, k, l
 
      SELECT CASE (potential%ppot_type)
      CASE (no_pseudo)
         WRITE (iw, '(/," ",28("*"),A,27("*"))') " All Electron Potential "
      CASE (gth_pseudo)
         WRITE (iw, '(/," ",29("*"),A,29("*"))') " GTH Pseudopotential "
         WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%gth_pot%zion
         WRITE (iw, '(T10,A,T66,F15.6)') " Rc ", potential%gth_pot%rc
         WRITE (pline, '(5F12.6)') (potential%gth_pot%cl(i), i=1, potential%gth_pot%ncl)
         WRITE (iw, '(T10,A,T21,A60)') " C1 C2 ... ", adjustr(pline)
         IF (potential%gth_pot%lpotextended) THEN
            DO k = 1, potential%gth_pot%nexp_lpot
               WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LPot: rc=", potential%gth_pot%alpha_lpot(k), &
                  "CX=", (potential%gth_pot%cval_lpot(i, k), i=1, potential%gth_pot%nct_lpot(k))
            END DO
         END IF
         IF (potential%gth_pot%nlcc) THEN
            DO k = 1, potential%gth_pot%nexp_nlcc
               WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LSDPot: rc=", potential%gth_pot%alpha_nlcc(k), &
                  "CX=", (potential%gth_pot%cval_nlcc(i, k)*4.0_dp*pi, i=1, potential%gth_pot%nct_nlcc(k))
            END DO
         END IF
         IF (potential%gth_pot%lsdpot) THEN
            DO k = 1, potential%gth_pot%nexp_lsd
               WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LSDPot: rc=", potential%gth_pot%alpha_lsd(k), &
                  "CX=", (potential%gth_pot%cval_lsd(i, k), i=1, potential%gth_pot%nct_lsd(k))
            END DO
         END IF
         DO l = 0, lmat
            IF (potential%gth_pot%nl(l) > 0) THEN
               WRITE (iw, '(T10,A,T76,I5)') " Angular momentum ", l
               WRITE (iw, '(T10,A,T66,F15.6)') " Rcnl ", potential%gth_pot%rcnl(l)
               WRITE (iw, '(T10,A,T76,I5)') " Nl ", potential%gth_pot%nl(l)
               WRITE (pline, '(5F12.6)') (potential%gth_pot%hnl(1, j, l), j=1, potential%gth_pot%nl(l))
               WRITE (iw, '(T10,A,T21,A60)') " Hnl ", adjustr(pline)
               DO i = 2, potential%gth_pot%nl(l)
                  WRITE (pline, '(T21,5F12.6)') (potential%gth_pot%hnl(i, j, l), j=i, potential%gth_pot%nl(l))
                  WRITE (iw, '(T21,A60)') adjustr(pline)
               END DO
            END IF
         END DO
      CASE (upf_pseudo)
         WRITE (iw, '(/," ",29("*"),A,29("*"))') " UPF Pseudopotential "
         DO k = 1, potential%upf_pot%maxinfo
            WRITE (iw, '(A80)') potential%upf_pot%info(k)
         END DO
      CASE (sgp_pseudo)
         WRITE (iw, '(/," ",29("*"),A,29("*"))') " SGP Pseudopotential "
         WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%sgp_pot%zion
      CASE (ecp_pseudo)
         WRITE (iw, '(/," ",26("*"),A,27("*"))') " Effective Core Potential "
         WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%ecp_pot%zion
         DO k = 1, potential%ecp_pot%nloc
            IF (k == 1) THEN
               WRITE (iw, '(T10,A,T40,I3,T49,2F16.8)') " Local Potential ", potential%ecp_pot%nrloc(k), &
                  potential%ecp_pot%bloc(k), potential%ecp_pot%aloc(k)
            ELSE
               WRITE (iw, '(T40,I3,T49,2F16.8)') potential%ecp_pot%nrloc(k), &
                  potential%ecp_pot%bloc(k), potential%ecp_pot%aloc(k)
            END IF
         END DO
         DO l = 0, potential%ecp_pot%lmax
            WRITE (iw, '(T10,A,I3)') " ECP l-value ", l
            DO k = 1, potential%ecp_pot%npot(l)
               WRITE (iw, '(T40,I3,T49,2F16.8)') potential%ecp_pot%nrpot(k, l), &
                  potential%ecp_pot%bpot(k, l), potential%ecp_pot%apot(k, l)
            END DO
         END DO
      CASE DEFAULT
         cpabort("")
      END SELECT
      IF (potential%confinement) THEN
         IF (potential%conf_type == poly_conf) THEN
            WRITE (iw, '(/,T10,A,T51,F12.6," * (R /",F6.2,")**",F6.2)') &
               " Confinement Potential ", potential%acon, potential%rcon, potential%scon
         ELSE IF (potential%conf_type == barrier_conf) THEN
            WRITE (iw, '(/,T10,A)') " Confinement Potential s*F[(r-ron)/w] "
            WRITE (iw, '(T57,A,F12.6,A)') "s     =", potential%acon, "   Ha"
            WRITE (iw, '(T57,A,F12.6,A)') "w     =", potential%rcon, " Bohr"
            WRITE (iw, '(T57,A,F12.6,A)') "ron   =", potential%scon, " Bohr"
         ELSE
            cpabort("")
         END IF
      ELSE
         WRITE (iw, '(/,T10,A)') " No Confinement Potential is applied "
      END IF
      WRITE (iw, '(" ",79("*"))')
 
   END SUBROUTINE atom_print_potential
 
! **************************************************************************************************
!> \brief Print GTH pseudo-potential parameters.
!> \param gthpot  pseudo-potential
!> \param iunit   output file unit
!> \par History
!>    * 09.2012 created [Juerg Hutter]
!> \note  The pseudo-potential is written into the 'iunit' file unit or as the file 'GTH-PARAMETER'
!>        inside the current working directory if the I/O unit is not given explicitly.
! **************************************************************************************************
   SUBROUTINE atom_write_pseudo_param(gthpot, iunit)
      TYPE(atom_gthpot_type), INTENT(INOUT)              :: gthpot
      INTEGER, INTENT(IN), OPTIONAL                      :: iunit
 
      INTEGER                                            :: i, iw, j, k, n
 
      IF (PRESENT(iunit)) THEN
         iw = iunit
      ELSE
         CALL open_file(file_name="GTH-PARAMETER", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
      END IF
      WRITE (iw, '(A2,A)') gthpot%symbol, adjustl(trim(gthpot%pname))
      WRITE (iw, '(4I5)') gthpot%econf(0:3)
      WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rc, gthpot%ncl, (gthpot%cl(i), i=1, gthpot%ncl)
      IF (gthpot%lpotextended) THEN
         WRITE (iw, '(A,I5)') "  LPOT", gthpot%nexp_lpot
         DO i = 1, gthpot%nexp_lpot
            WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lpot(i), gthpot%nct_lpot(i), &
               (gthpot%cval_lpot(j, i), j=1, gthpot%nct_lpot(i))
         END DO
      END IF
      IF (gthpot%lsdpot) THEN
         WRITE (iw, '(A,I5)') "  LSD ", gthpot%nexp_lsd
         DO i = 1, gthpot%nexp_lsd
            WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lsd(i), gthpot%nct_lsd(i), &
               (gthpot%cval_lsd(j, i), j=1, gthpot%nct_lsd(i))
         END DO
      END IF
      IF (gthpot%nlcc) THEN
         WRITE (iw, '(A,I5)') " NLCC ", gthpot%nexp_nlcc
         DO i = 1, gthpot%nexp_nlcc
            WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_nlcc(i), gthpot%nct_nlcc(i), &
               (gthpot%cval_nlcc(j, i)*4.0_dp*pi, j=1, gthpot%nct_nlcc(i))
         END DO
      END IF
      n = 0
      DO i = lmat, 0, -1
         IF (gthpot%nl(i) > 0) THEN
            n = i + 1
            EXIT
         END IF
      END DO
      WRITE (iw, '(I8)') n
      DO i = 0, n - 1
         WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rcnl(i), gthpot%nl(i), (gthpot%hnl(1, k, i), k=1, gthpot%nl(i))
         SELECT CASE (gthpot%nl(i))
         CASE (2)
            WRITE (iw, '(T49,F20.14)') gthpot%hnl(2, 2, i)
         CASE (3)
            WRITE (iw, '(T49,2F20.14)') gthpot%hnl(2, 2, i), gthpot%hnl(2, 3, i)
            WRITE (iw, '(T69,F20.14)') gthpot%hnl(3, 3, i)
         CASE DEFAULT
            DO j = 2, gthpot%nl(i)
               WRITE (iw, '(T29,5F20.14)') (gthpot%hnl(j, k, i), k=j, gthpot%nl(i))
            END DO
         END SELECT
      END DO
      IF (.NOT. PRESENT(iunit)) CALL close_file(unit_number=iw)
 
   END SUBROUTINE atom_write_pseudo_param
 
! **************************************************************************************************
!> \brief Print atomic orbitals.
!> \param atom  information about the atomic kind
!> \param iw    output file unit
!> \param xmgrace ...
!> \par History
!>    * 04.2013 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_orbitals(atom, iw, xmgrace)
      TYPE(atom_type), POINTER                           :: atom
      INTEGER, INTENT(IN)                                :: iw
      LOGICAL, INTENT(IN), OPTIONAL                      :: xmgrace
 
      CHARACTER(LEN=40)                                  :: fnbody
      INTEGER                                            :: z
      LOGICAL                                            :: graph
 
      SELECT CASE (atom%method_type)
      CASE DEFAULT
         cpabort("")
      CASE (do_rks_atom)
         CALL atom_print_orbitals_helper(atom, atom%orbitals%wfn, "", iw)
      CASE (do_uks_atom)
         CALL atom_print_orbitals_helper(atom, atom%orbitals%wfna, "Alpha", iw)
         CALL atom_print_orbitals_helper(atom, atom%orbitals%wfnb, "Beta", iw)
      CASE (do_rhf_atom)
         CALL atom_print_orbitals_helper(atom, atom%orbitals%wfn, "", iw)
      CASE (do_uhf_atom)
         CALL atom_print_orbitals_helper(atom, atom%orbitals%wfna, "Alpha", iw)
         CALL atom_print_orbitals_helper(atom, atom%orbitals%wfnb, "Beta", iw)
      CASE (do_rohf_atom)
         cpabort("")
      END SELECT
 
      graph = .false.
      IF (PRESENT(xmgrace)) graph = xmgrace
      IF (graph .AND. iw > 0) THEN
         z = atom%z
         fnbody = trim(ptable(z)%symbol)//"_PPorbital"
         SELECT CASE (atom%method_type)
         CASE DEFAULT
            cpabort("")
         CASE (do_rks_atom)
            CALL atom_orbitals_grace(atom, atom%orbitals%wfn, fnbody)
         CASE (do_uks_atom)
            CALL atom_orbitals_grace(atom, atom%orbitals%wfna, trim(fnbody)//"alpha")
            CALL atom_orbitals_grace(atom, atom%orbitals%wfnb, trim(fnbody)//"beta")
         CASE (do_rhf_atom)
            CALL atom_orbitals_grace(atom, atom%orbitals%wfn, fnbody)
         CASE (do_uhf_atom)
            CALL atom_orbitals_grace(atom, atom%orbitals%wfna, trim(fnbody)//"alpha")
            CALL atom_orbitals_grace(atom, atom%orbitals%wfnb, trim(fnbody)//"beta")
         CASE (do_rohf_atom)
            cpabort("")
         END SELECT
      END IF
 
   END SUBROUTINE atom_print_orbitals
 
! **************************************************************************************************
!> \brief Print atomic orbitals of the given spin.
!> \param atom         information about the atomic kind
!> \param wfn          atomic orbitals
!> \param description  description string
!> \param iw           output file unit
!> \par History
!>    * 04.2013 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_print_orbitals_helper(atom, wfn, description, iw)
      TYPE(atom_type), POINTER                           :: atom
      REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT)  :: wfn
      CHARACTER(len=*), INTENT(IN)                       :: description
      INTEGER, INTENT(IN)                                :: iw
 
      INTEGER                                            :: b, l, maxl, nb, nv, v
 
      WRITE (iw, '(/,A,A,A)') " Atomic orbital expansion coefficients [", description, "]"
 
      maxl = atom%state%maxl_calc
      DO l = 0, maxl
 
         nb = atom%basis%nbas(l)
         nv = atom%state%maxn_calc(l)
         IF (nb > 0 .AND. nv > 0) THEN
            nv = min(nv, SIZE(wfn, 2))
            DO v = 1, nv
               WRITE (iw, '(/,"    ORBITAL      L = ",I1,"      State = ",I3)') l, v
               DO b = 1, nb
                  WRITE (iw, '("      ",ES23.15)') wfn(b, v, l)
               END DO
            END DO
         END IF
      END DO
   END SUBROUTINE atom_print_orbitals_helper
 
! **************************************************************************************************
!> \brief Print atomic orbitals of the given spin.
!> \param atom         information about the atomic kind
!> \param wfn          atomic orbitals
!> \param fnbody       body of file name
!> \par History
!>    * 02.2025 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_orbitals_grace(atom, wfn, fnbody)
      TYPE(atom_type), POINTER                           :: atom
      REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT)  :: wfn
      CHARACTER(len=*), INTENT(IN)                       :: fnbody
 
      CHARACTER(LEN=1), DIMENSION(0:8)                   :: lname
      CHARACTER(LEN=1), DIMENSION(1:9)                   :: wnum
      CHARACTER(LEN=40)                                  :: fname, legend
      INTEGER                                            :: b, i, iw, l, m, maxl, nb, nv, v
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: gdata, wfnr
      REAL(kind=dp), DIMENSION(4)                        :: world_coord
 
      lname = ['s', 'p', 'd', 'f', 'g', 'h', 'j', 'k', 'l']
      wnum = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
      m = atom%basis%grid%nr
      maxl = atom%state%maxl_calc
      DO l = 0, maxl
         fname = trim(fnbody)//"_"//lname(l)//".agr"
         nb = atom%basis%nbas(l)
         nv = atom%state%maxn_calc(l)
         IF (nb > 0 .AND. nv > 0) THEN
            CALL open_file(file_name=fname, file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
            nv = min(nv, SIZE(wfn, 2))
            ALLOCATE (wfnr(m, nv))
            wfnr = 0.0_dp
            DO v = 1, nv
               DO b = 1, nb
                  wfnr(:, v) = wfnr(:, v) + wfn(b, v, l)*atom%basis%bf(:, b, l)
               END DO
            END DO
            world_coord(1) = 0.0_dp
            world_coord(2) = minval(wfnr) - 0.5_dp
            world_coord(3) = 15.0_dp
            world_coord(4) = maxval(wfnr) + 0.5_dp
            !
            CALL xm_write_defaults(iw)
            CALL xm_write_frameport(iw)
            CALL xm_write_frame(iw, world_coord, &
                                title="PP Radial Wavefunction", &
                                subtitle=lname(l)//"-Quantum Number", &
                                xlabel="Radius [Bohr]", &
                                ylabel="")
            DO i = 0, nv - 1
               legend = "WFN "//wnum(i + 1)
               CALL xm_graph_info(iw, i, 2.5_dp, legend)
            END DO
            ALLOCATE (gdata(m, 2))
            gdata(1:m, 1) = atom%basis%grid%rad(1:m)
            DO i = 0, nv - 1
               gdata(1:m, 2) = wfnr(1:m, i + 1)
               CALL xm_graph_data(iw, i, gdata)
            END DO
            DEALLOCATE (gdata, wfnr)
            CALL close_file(iw)
         END IF
      END DO
   END SUBROUTINE atom_orbitals_grace
 
END MODULE atom_output