atom_energy.F

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

MODULE atom_energy
   USE atom_admm_methods,               ONLY: atom_admm
   USE atom_electronic_structure,       ONLY: calculate_atom
   USE atom_fit,                        ONLY: atom_fit_density,&
                                              atom_fit_kgpot
   USE atom_grb,                        ONLY: atom_grb_construction
   USE atom_operators,                  ONLY: atom_int_release,&
                                              atom_int_setup,&
                                              atom_ppint_release,&
                                              atom_ppint_setup,&
                                              atom_relint_release,&
                                              atom_relint_setup
   USE atom_output,                     ONLY: atom_print_basis,&
                                              atom_print_info,&
                                              atom_print_method,&
                                              atom_print_orbitals,&
                                              atom_print_potential,&
                                              atom_write_pseudo_param
   USE atom_sgp,                        ONLY: atom_sgp_construction
   USE atom_types,                      ONLY: &
        atom_basis_type, atom_gthpot_type, atom_integrals, atom_optimization_type, atom_orbitals, &
        atom_p_type, atom_potential_type, atom_state, atom_type, create_atom_orbs, &
        create_atom_type, gth_pseudo, init_atom_basis, init_atom_potential, lmat, &
        read_atom_opt_section, release_atom_basis, release_atom_potential, release_atom_type, &
        set_atom
   USE atom_utils,                      ONLY: &
        atom_completeness, atom_condnumber, atom_consistent_method, atom_core_density, &
        atom_density, atom_local_potential, atom_read_external_density, atom_read_external_vxc, &
        atom_set_occupation, atom_write_zmp_restart, get_maxl_occ, get_maxn_occ
   USE atom_xc,                         ONLY: calculate_atom_ext_vxc,&
                                              calculate_atom_zmp
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
                                              cp_print_key_unit_nr
   USE input_constants,                 ONLY: &
        do_analytic, xc_funct_b3lyp, xc_funct_beefvdw, xc_funct_blyp, xc_funct_bp, &
        xc_funct_hcth120, xc_funct_no_shortcut, xc_funct_olyp, xc_funct_pade, xc_funct_pbe, &
        xc_funct_pbe0, xc_funct_tpss, xc_none
   USE input_section_types,             ONLY: section_vals_get,&
                                              section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: default_string_length,&
                                              dp
   USE mathconstants,                   ONLY: dfac,&
                                              fourpi,&
                                              gamma1,&
                                              rootpi
   USE periodic_table,                  ONLY: nelem,&
                                              ptable
   USE physcon,                         ONLY: bohr
#include "./base/base_uses.f90"

   IMPLICIT NONE
   PRIVATE
   PUBLIC  :: atom_energy_opt

   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_energy'

CONTAINS

! **************************************************************************************************
!> \brief Compute the atomic energy.
!> \param atom_section  ATOM input section
!> \par History
!>    * 11.2016 geometrical response basis set [Juerg Hutter]
!>    * 07.2016 ADMM analysis [Juerg Hutter]
!>    * 02.2014 non-additive kinetic energy term [Juerg Hutter]
!>    * 11.2013 Zhao, Morrison, and Parr (ZMP) potential [Daniele Varsano]
!>    * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
!>    * 05.2009 electronic density fit [Juerg Hutter]
!>    * 04.2009 refactored and renamed to atom_energy_opt() [Juerg Hutter]
!>    * 08.2008 created as atom_energy() [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_energy_opt(atom_section)
      TYPE(section_vals_type), POINTER                   :: atom_section

      CHARACTER(len=*), PARAMETER                        :: routineN = 'atom_energy_opt'

      CHARACTER(LEN=2)                                   :: elem
      CHARACTER(LEN=default_string_length)               :: filename, xcfstr
      CHARACTER(LEN=default_string_length), &
         DIMENSION(:), POINTER                           :: tmpstringlist
      INTEGER :: do_eric, do_erie, handle, i, im, in, iw, k, maxl, mb, method, mo, n_meth, n_rep, &
         nder, nr_gh, num_gau, num_gto, num_pol, reltyp, shortcut, zcore, zval, zz
      INTEGER, DIMENSION(0:lmat)                         :: maxn
      INTEGER, DIMENSION(:), POINTER                     :: cn
      LOGICAL                                            :: dm, do_gh, do_zmp, doread, eri_c, eri_e, &
                                                            explicit, graph, had_ae, had_pp, &
                                                            lcomp, lcond, pp_calc, read_vxc
      REAL(KIND=dp)                                      :: crad, delta, lambda
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ext_density, ext_vxc
      REAL(KIND=dp), DIMENSION(0:lmat, 10)               :: pocc
      TYPE(atom_basis_type), POINTER                     :: ae_basis, pp_basis
      TYPE(atom_integrals), POINTER                      :: ae_int, pp_int
      TYPE(atom_optimization_type)                       :: optimization
      TYPE(atom_orbitals), POINTER                       :: orbitals
      TYPE(atom_p_type), DIMENSION(:, :), POINTER        :: atom_info
      TYPE(atom_potential_type), POINTER                 :: ae_pot, p_pot
      TYPE(atom_state), POINTER                          :: state
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(section_vals_type), POINTER :: admm_section, basis_section, external_vxc_section, &
         method_section, opt_section, potential_section, powell_section, sgp_section, xc_func, &
         xc_section, zmp_restart_section, zmp_section

      CALL timeset(routineN, handle)

      ! What atom do we calculate
      CALL section_vals_val_get(atom_section, "ATOMIC_NUMBER", i_val=zval)
      CALL section_vals_val_get(atom_section, "ELEMENT", c_val=elem)
      zz = 0
      DO i = 1, nelem
         IF (ptable(i)%symbol == elem) THEN
            zz = i
            EXIT
         END IF
      END DO
      IF (zz /= 1) zval = zz

      ! read and set up inofrmation on the basis sets
      ALLOCATE (ae_basis, pp_basis)
      basis_section => section_vals_get_subs_vals(atom_section, "AE_BASIS")
      NULLIFY (ae_basis%grid)
      CALL init_atom_basis(ae_basis, basis_section, zval, "AE")
      NULLIFY (pp_basis%grid)
      basis_section => section_vals_get_subs_vals(atom_section, "PP_BASIS")
      CALL init_atom_basis(pp_basis, basis_section, zval, "PP")

      ! print general and basis set information
      logger => cp_get_default_logger()
      iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log")
      IF (iw > 0) CALL atom_print_info(zval, "Atomic Energy Calculation", iw)
      CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
      iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log")
      IF (iw > 0) THEN
         CALL atom_print_basis(ae_basis, iw, " All Electron Basis")
         CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis")
      END IF
      CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET")

      ! read and setup information on the pseudopotential
      NULLIFY (potential_section)
      potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL")
      ALLOCATE (ae_pot, p_pot)
      CALL init_atom_potential(p_pot, potential_section, zval)
      CALL init_atom_potential(ae_pot, potential_section, -1)

      ! if the ERI's are calculated analytically, we have to precalculate them
      eri_c = .FALSE.
      CALL section_vals_val_get(atom_section, "COULOMB_INTEGRALS", i_val=do_eric)
      IF (do_eric == do_analytic) eri_c = .TRUE.
      eri_e = .FALSE.
      CALL section_vals_val_get(atom_section, "EXCHANGE_INTEGRALS", i_val=do_erie)
      IF (do_erie == do_analytic) eri_e = .TRUE.
      CALL section_vals_val_get(atom_section, "USE_GAUSS_HERMITE", l_val=do_gh)
      CALL section_vals_val_get(atom_section, "GRID_POINTS_GH", i_val=nr_gh)

      ! information on the states to be calculated
      CALL section_vals_val_get(atom_section, "MAX_ANGULAR_MOMENTUM", i_val=maxl)
      maxn = 0
      CALL section_vals_val_get(atom_section, "CALCULATE_STATES", i_vals=cn)
      DO in = 1, MIN(SIZE(cn), 4)
         maxn(in - 1) = cn(in)
      END DO
      DO in = 0, lmat
         maxn(in) = MIN(maxn(in), ae_basis%nbas(in))
         maxn(in) = MIN(maxn(in), pp_basis%nbas(in))
      END DO

      ! read optimization section
      opt_section => section_vals_get_subs_vals(atom_section, "OPTIMIZATION")
      CALL read_atom_opt_section(optimization, opt_section)

      had_ae = .FALSE.
      had_pp = .FALSE.

      ! Check for the total number of electron configurations to be calculated
      CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", n_rep_val=n_rep)
      ! Check for the total number of method types to be calculated
      method_section => section_vals_get_subs_vals(atom_section, "METHOD")
      CALL section_vals_get(method_section, n_repetition=n_meth)

      ! integrals
      ALLOCATE (ae_int, pp_int)

      ALLOCATE (atom_info(n_rep, n_meth))

      DO in = 1, n_rep
         DO im = 1, n_meth

            NULLIFY (atom_info(in, im)%atom)
            CALL create_atom_type(atom_info(in, im)%atom)

            atom_info(in, im)%atom%optimization = optimization

            atom_info(in, im)%atom%z = zval

            xc_section => section_vals_get_subs_vals(method_section, "XC", i_rep_section=im)
            atom_info(in, im)%atom%xc_section => xc_section

            ! Retrieve the XC functional string for upf output
            ! This does not work within the routine atom_write_upf below for unknown reasons
            xc_func => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
            CALL section_vals_get(xc_func, explicit=explicit)
            IF (explicit) THEN
               CALL section_vals_val_get(xc_func, "_SECTION_PARAMETERS_", i_val=shortcut)
               SELECT CASE (shortcut)
               CASE (xc_funct_no_shortcut, xc_none)
                  xcfstr = "NONE"
               CASE (xc_funct_pbe0)
                  xcfstr = "PBE0"
               CASE (xc_funct_beefvdw)
                  xcfstr = "BEEFVDW"
               CASE (xc_funct_b3lyp)
                  xcfstr = "B3LYP"
               CASE (xc_funct_blyp)
                  xcfstr = "BLYP"
               CASE (xc_funct_bp)
                  xcfstr = "BP"
               CASE (xc_funct_pade)
                  xcfstr = "PADE"
               CASE (xc_funct_pbe)
                  xcfstr = "PBE"
               CASE (xc_funct_tpss)
                  xcfstr = "TPSS"
               CASE (xc_funct_olyp)
                  xcfstr = "OLYP"
               CASE (xc_funct_hcth120)
                  xcfstr = "HCTH120"
               CASE DEFAULT
                  xcfstr = "DFT"
               END SELECT
            ELSE
               xcfstr = "DFT"
            END IF

            ! ZMP Reading input sections if they are found initialize everything
            do_zmp = .FALSE.
            doread = .FALSE.
            read_vxc = .FALSE.

            zmp_section => section_vals_get_subs_vals(method_section, "ZMP")
            CALL section_vals_get(zmp_section, explicit=do_zmp)
            atom_info(in, im)%atom%do_zmp = do_zmp
            CALL section_vals_val_get(zmp_section, "FILE_DENSITY", c_val=filename)
            atom_info(in, im)%atom%ext_file = filename
            CALL section_vals_val_get(zmp_section, "GRID_TOL", &
                                      r_val=atom_info(in, im)%atom%zmpgrid_tol)
            CALL section_vals_val_get(zmp_section, "LAMBDA", r_val=lambda)
            atom_info(in, im)%atom%lambda = lambda
            CALL section_vals_val_get(zmp_section, "DM", l_val=dm)
            atom_info(in, im)%atom%dm = dm

            zmp_restart_section => section_vals_get_subs_vals(zmp_section, "RESTART")
            CALL section_vals_get(zmp_restart_section, explicit=doread)
            atom_info(in, im)%atom%doread = doread
            CALL section_vals_val_get(zmp_restart_section, "FILE_RESTART", c_val=filename)
            atom_info(in, im)%atom%zmp_restart_file = filename

            ! ZMP Reading external vxc section, if found initialize
            external_vxc_section => section_vals_get_subs_vals(method_section, "EXTERNAL_VXC")
            CALL section_vals_get(external_vxc_section, explicit=read_vxc)
            atom_info(in, im)%atom%read_vxc = read_vxc
            CALL section_vals_val_get(external_vxc_section, "FILE_VXC", c_val=filename)
            atom_info(in, im)%atom%ext_vxc_file = filename
            CALL section_vals_val_get(external_vxc_section, "GRID_TOL", &
                                      r_val=atom_info(in, im)%atom%zmpvxcgrid_tol)

            ALLOCATE (state)

            ! get the electronic configuration
            CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, &
                                      c_vals=tmpstringlist)

            ! set occupations
            CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity)
            state%maxl_occ = get_maxl_occ(state%occ)
            state%maxn_occ = get_maxn_occ(state%occ)

            ! set number of states to be calculated
            state%maxl_calc = MAX(maxl, state%maxl_occ)
            state%maxl_calc = MIN(lmat, state%maxl_calc)
            state%maxn_calc = 0
            DO k = 0, state%maxl_calc
               state%maxn_calc(k) = MAX(maxn(k), state%maxn_occ(k))
            END DO

            ! is there a pseudo potential
            pp_calc = ANY(INDEX(tmpstringlist(1:), "CORE") /= 0)
            IF (pp_calc) THEN
               ! get and set the core occupations
               CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
               CALL atom_set_occupation(tmpstringlist, state%core, pocc)
               zcore = zval - NINT(SUM(state%core))
               CALL set_atom(atom_info(in, im)%atom, zcore=zcore, pp_calc=.TRUE.)
            ELSE
               state%core = 0._dp
               CALL set_atom(atom_info(in, im)%atom, zcore=zval, pp_calc=.FALSE.)
            END IF

            CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_section=im)
            CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im)
            CALL set_atom(atom_info(in, im)%atom, method_type=method, relativistic=reltyp)

            IF (atom_consistent_method(method, state%multiplicity)) THEN
               iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%METHOD_INFO", extension=".log")
               CALL atom_print_method(atom_info(in, im)%atom, iw)
               CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%METHOD_INFO")

               iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log")
               IF (pp_calc) THEN
                  IF (iw > 0) CALL atom_print_potential(p_pot, iw)
               ELSE
                  IF (iw > 0) CALL atom_print_potential(ae_pot, iw)
               END IF
               CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL")
            END IF

            ! calculate integrals
            IF (pp_calc) THEN
               ! general integrals
               CALL atom_int_setup(pp_int, pp_basis, &
                                   potential=p_pot, eri_coulomb=eri_c, eri_exchange=eri_e)
               ! potential
               CALL atom_ppint_setup(pp_int, pp_basis, potential=p_pot)
               !
               NULLIFY (pp_int%tzora, pp_int%hdkh)
               !
               CALL set_atom(atom_info(in, im)%atom, basis=pp_basis, integrals=pp_int, potential=p_pot)
               state%maxn_calc(:) = MIN(state%maxn_calc(:), pp_basis%nbas(:))
               CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
               had_pp = .TRUE.
            ELSE
               ! general integrals
               CALL atom_int_setup(ae_int, ae_basis, potential=ae_pot, &
                                   eri_coulomb=eri_c, eri_exchange=eri_e)
               ! potential
               CALL atom_ppint_setup(ae_int, ae_basis, potential=ae_pot)
               ! relativistic correction terms
               CALL atom_relint_setup(ae_int, ae_basis, reltyp, zcore=REAL(zval, dp))
               !
               CALL set_atom(atom_info(in, im)%atom, basis=ae_basis, integrals=ae_int, potential=ae_pot)
               state%maxn_calc(:) = MIN(state%maxn_calc(:), ae_basis%nbas(:))
               CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
               had_ae = .TRUE.
            END IF

            CALL set_atom(atom_info(in, im)%atom, state=state)

            CALL set_atom(atom_info(in, im)%atom, coulomb_integral_type=do_eric, &
                          exchange_integral_type=do_erie)
            atom_info(in, im)%atom%hfx_pot%do_gh = do_gh
            atom_info(in, im)%atom%hfx_pot%nr_gh = nr_gh

            NULLIFY (orbitals)
            mo = MAXVAL(state%maxn_calc)
            mb = MAXVAL(atom_info(in, im)%atom%basis%nbas)
            CALL create_atom_orbs(orbitals, mb, mo)
            CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals)

            IF (atom_consistent_method(method, state%multiplicity)) THEN
               !Calculate the electronic structure
               iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SCF_INFO", extension=".log")
               CALL calculate_atom(atom_info(in, im)%atom, iw)
               CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SCF_INFO")

               ! ZMP If we have the external density do zmp
               IF (atom_info(in, im)%atom%do_zmp) THEN
                  CALL atom_write_zmp_restart(atom_info(in, im)%atom)

                  ALLOCATE (ext_density(atom_info(in, im)%atom%basis%grid%nr))
                  ext_density = 0._dp
                  CALL atom_read_external_density(ext_density, atom_info(in, im)%atom, iw)
                  CALL calculate_atom_zmp(ext_density=ext_density, &
                                          atom=atom_info(in, im)%atom, &
                                          lprint=.TRUE.)
                  DEALLOCATE (ext_density)
               END IF
               ! ZMP If we have the external v_xc calculate KS quantities
               IF (atom_info(in, im)%atom%read_vxc) THEN
                  ALLOCATE (ext_vxc(atom_info(in, im)%atom%basis%grid%nr))
                  ext_vxc = 0._dp
                  CALL atom_read_external_vxc(ext_vxc, atom_info(in, im)%atom, iw)
                  CALL calculate_atom_ext_vxc(vxc=ext_vxc, &
                                              atom=atom_info(in, im)%atom, &
                                              lprint=.TRUE.)
                  DEALLOCATE (ext_vxc)
               END IF

               ! Print out the orbitals if requested
               iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ORBITALS", extension=".log")
               CALL section_vals_val_get(atom_section, "PRINT%ORBITALS%XMGRACE", l_val=graph)
               IF (iw > 0) THEN
                  CALL atom_print_orbitals(atom_info(in, im)%atom, iw, xmgrace=graph)
               END IF
               CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ORBITALS")

               ! generate a UPF file
               iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%UPF_FILE", extension=".upf", &
                                         file_position="REWIND")
               IF (iw > 0) THEN
                  CALL atom_write_upf(atom_info(in, im)%atom, iw, xcfstr)
               END IF
               CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%UPF_FILE")

               ! perform a fit of the total electronic density
               iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_DENSITY", extension=".log")
               IF (iw > 0) THEN
                  CALL section_vals_val_get(atom_section, "PRINT%FIT_DENSITY%NUM_GTO", i_val=num_gto)
                  powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
                  CALL atom_fit_density(atom_info(in, im)%atom, num_gto, 0, iw, powell_section=powell_section)
               END IF
               CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_DENSITY")

               ! Optimize a local potential for the non-additive kinetic energy term in KG
               iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_KGPOT", extension=".log")
               IF (iw > 0) THEN
                  CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_GAUSSIAN", i_val=num_gau)
                  CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_POLYNOM", i_val=num_pol)
                  powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
                  CALL atom_fit_kgpot(atom_info(in, im)%atom, num_gau, num_pol, iw, powell_section=powell_section)
               END IF
               CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_KGPOT")

               ! generate a response basis
               iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%RESPONSE_BASIS", extension=".log")
               IF (iw > 0) THEN
                  CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DELTA_CHARGE", r_val=delta)
                  CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DERIVATIVES", i_val=nder)
                  CALL atom_response_basis(atom_info(in, im)%atom, delta, nder, iw)
               END IF
               CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%RESPONSE_BASIS")

            END IF

         END DO
      END DO

      ! generate a geometrical response basis (GRB)
      iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS", extension=".log")
      IF (iw > 0) THEN
         CALL atom_grb_construction(atom_info, atom_section, iw)
      END IF
      CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")

      ! Analyze basis sets
      iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ANALYZE_BASIS", extension=".log")
      IF (iw > 0) THEN
         CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%OVERLAP_CONDITION_NUMBER", l_val=lcond)
         CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%COMPLETENESS", l_val=lcomp)
         crad = ptable(zval)%covalent_radius*bohr
         IF (had_ae) THEN
            IF (lcond) CALL atom_condnumber(ae_basis, crad, iw)
            IF (lcomp) CALL atom_completeness(ae_basis, zval, iw)
         END IF
         IF (had_pp) THEN
            IF (lcond) CALL atom_condnumber(pp_basis, crad, iw)
            IF (lcomp) CALL atom_completeness(pp_basis, zval, iw)
         END IF
      END IF
      CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ANALYZE_BASIS")

      ! Analyse ADMM approximation
      iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ADMM", extension=".log")
      admm_section => section_vals_get_subs_vals(atom_section, "PRINT%ADMM")
      IF (iw > 0) THEN
         CALL atom_admm(atom_info, admm_section, iw)
      END IF
      CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ADMM")

      ! Generate a separable form of the pseudopotential using Gaussian functions
      iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO", extension=".log")
      sgp_section => section_vals_get_subs_vals(atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
      IF (iw > 0) THEN
         CALL atom_sgp_construction(atom_info, sgp_section, iw)
      END IF
      CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")

      ! clean up
      IF (had_ae) THEN
         CALL atom_int_release(ae_int)
         CALL atom_ppint_release(ae_int)
         CALL atom_relint_release(ae_int)
      END IF
      IF (had_pp) THEN
         CALL atom_int_release(pp_int)
         CALL atom_ppint_release(pp_int)
         CALL atom_relint_release(pp_int)
      END IF
      CALL release_atom_basis(ae_basis)
      CALL release_atom_basis(pp_basis)

      CALL release_atom_potential(p_pot)
      CALL release_atom_potential(ae_pot)

      DO in = 1, n_rep
         DO im = 1, n_meth
            CALL release_atom_type(atom_info(in, im)%atom)
         END DO
      END DO
      DEALLOCATE (atom_info)

      DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)

      CALL timestop(handle)

   END SUBROUTINE atom_energy_opt

! **************************************************************************************************
!> \brief Calculate response basis contraction coefficients.
!> \param atom   information about the atomic kind
!> \param delta  variation of charge used in finite difference derivative calculation
!> \param nder   number of wavefunction derivatives to calculate
!> \param iw     output file unit
!> \par History
!>    * 10.2011 Gram-Schmidt orthogonalization [Joost VandeVondele]
!>    * 05.2009 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_response_basis(atom, delta, nder, iw)

      TYPE(atom_type), POINTER                           :: atom
      REAL(KIND=dp), INTENT(IN)                          :: delta
      INTEGER, INTENT(IN)                                :: nder, iw

      INTEGER                                            :: i, ider, j, k, l, lhomo, m, n, nhomo, &
                                                            s1, s2
      REAL(KIND=dp)                                      :: dene, emax, expzet, fhomo, o, prefac, &
                                                            zeta
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: amat
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: rbasis
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: wfn
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: ovlp
      TYPE(atom_state), POINTER                          :: state

      WRITE (iw, '(/," ",79("*"),/,T30,A,A/," ",79("*"))') "RESPONSE BASIS for ", ptable(atom%z)%symbol

      state => atom%state
      ovlp => atom%integrals%ovlp

      ! find HOMO
      lhomo = -1
      nhomo = -1
      emax = -HUGE(1._dp)
      DO l = 0, state%maxl_occ
         DO i = 1, state%maxn_occ(l)
            IF (atom%orbitals%ener(i, l) > emax) THEN
               lhomo = l
               nhomo = i
               emax = atom%orbitals%ener(i, l)
               fhomo = state%occupation(l, i)
            END IF
         END DO
      END DO

      s1 = SIZE(atom%orbitals%wfn, 1)
      s2 = SIZE(atom%orbitals%wfn, 2)
      ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
      s2 = MAXVAL(state%maxn_occ) + nder
      ALLOCATE (rbasis(s1, s2, 0:lmat))
      rbasis = 0._dp

      DO ider = -nder, nder
         dene = REAL(ider, KIND=dp)*delta
         CPASSERT(fhomo > ABS(dene))
         state%occupation(lhomo, nhomo) = fhomo + dene
         CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
         wfn(:, :, :, ider) = atom%orbitals%wfn
         state%occupation(lhomo, nhomo) = fhomo
      END DO
      ! reset
      CALL calculate_atom(atom, iw=0, noguess=.TRUE.)

      DO l = 0, state%maxl_occ
         ! occupied states
         DO i = 1, MAX(state%maxn_occ(l), 1)
            rbasis(:, i, l) = wfn(:, i, l, 0)
         END DO
         ! differentiation
         DO ider = 1, nder
            i = MAX(state%maxn_occ(l), 1)
            SELECT CASE (ider)
            CASE (1)
               rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
            CASE (2)
               rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
            CASE (3)
               rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
                                               + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
            CASE DEFAULT
               CPABORT("")
            END SELECT
         END DO

         ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
         n = state%maxn_occ(l) + nder
         m = atom%basis%nbas(l)
         DO i = 1, n
            DO j = 1, i - 1
               o = DOT_PRODUCT(rbasis(1:m, j, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
               rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
            END DO
            o = DOT_PRODUCT(rbasis(1:m, i, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
            rbasis(1:m, i, l) = rbasis(1:m, i, l)/SQRT(o)
         END DO

         ! check
         ALLOCATE (amat(n, n))
         amat(1:n, 1:n) = MATMUL(TRANSPOSE(rbasis(1:m, 1:n, l)), MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, 1:n, l)))
         DO i = 1, n
            amat(i, i) = amat(i, i) - 1._dp
         END DO
         IF (MAXVAL(ABS(amat)) > 1.e-12) THEN
            WRITE (iw, '(/,A,G20.10)') " Orthogonality error  ", MAXVAL(ABS(amat))
         END IF
         DEALLOCATE (amat)

         ! Quickstep normalization
         WRITE (iw, '(/,A,T30,I3)') " Angular momentum :", l

         WRITE (iw, '(/,A,I0,A,I0,A)') "             Exponent     Coef.(Quickstep Normalization), first ", &
            n - nder, " valence ", nder, " response"
         expzet = 0.25_dp*REAL(2*l + 3, dp)
         prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
         DO i = 1, m
            zeta = (2._dp*atom%basis%am(i, l))**expzet
            WRITE (iw, '(4X,F20.10,4X,15ES20.6)') atom%basis%am(i, l), ((prefac*rbasis(i, k, l)/zeta), k=1, n)
         END DO

      END DO

      DEALLOCATE (wfn, rbasis)

      WRITE (iw, '(" ",79("*"))')

   END SUBROUTINE atom_response_basis

! **************************************************************************************************
!> \brief Write pseudo-potential in Quantum Espresso UPF format.
!> \param atom  information about the atomic kind
!> \param iw    output file unit
!> \param xcfstr ...
!> \par History
!>    * 09.2012 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_write_upf(atom, iw, xcfstr)

      TYPE(atom_type), POINTER                           :: atom
      INTEGER, INTENT(IN)                                :: iw
      CHARACTER(LEN=*), INTENT(IN)                       :: xcfstr

      CHARACTER(LEN=default_string_length)               :: string
      INTEGER                                            :: i, ibeta, ii, j, k, kk, l, lmax, n0, &
                                                            nbeta, nr, nwfc, nwfn
      LOGICAL                                            :: soc, up
      REAL(KIND=dp)                                      :: occa, occb, pf, rhoatom, rl, rlp, rmax, &
                                                            vnlm, vnlp
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: beta, corden, dens, ef, locpot, rp
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dij
      TYPE(atom_gthpot_type), POINTER                    :: pot

      IF (.NOT. atom%pp_calc) RETURN
      IF (atom%potential%ppot_type /= gth_pseudo) RETURN
      pot => atom%potential%gth_pot
      CPASSERT(.NOT. pot%lsdpot)
      soc = pot%soc

      WRITE (iw, '(A)') '<UPF version="2.0.1">'

      WRITE (iw, '(T4,A)') '<PP_INFO>'
      WRITE (iw, '(T8,A)') 'Converted from CP2K GTH format'
      WRITE (iw, '(T8,A)') '<PP_INPUTFILE>'
      CALL atom_write_pseudo_param(pot, iw)
      WRITE (iw, '(T8,A)') '</PP_INPUTFILE>'
      WRITE (iw, '(T4,A)') '</PP_INFO>'
      WRITE (iw, '(T4,A)') '<PP_HEADER'
      WRITE (iw, '(T8,A)') 'generated="Generated in analytical, separable form"'
      WRITE (iw, '(T8,A)') 'author="Goedecker/Hartwigsen/Hutter/Teter"'
      WRITE (iw, '(T8,A)') 'date="Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"'
      WRITE (iw, '(T8,A)') 'comment="This file generated by CP2K ATOM code"'
      CALL compose(string, "element", cval=pot%symbol)
      WRITE (iw, '(T8,A)') TRIM(string)
      WRITE (iw, '(T8,A)') 'pseudo_type="NC"'
      WRITE (iw, '(T8,A)') 'relativistic="no"'
      WRITE (iw, '(T8,A)') 'is_ultrasoft="F"'
      WRITE (iw, '(T8,A)') 'is_paw="F"'
      WRITE (iw, '(T8,A)') 'is_coulomb="F"'
      IF (soc) THEN
         WRITE (iw, '(T8,A)') 'has_so="T"'
      ELSE
         WRITE (iw, '(T8,A)') 'has_so="F"'
      END IF
      WRITE (iw, '(T8,A)') 'has_wfc="F"'
      WRITE (iw, '(T8,A)') 'has_gipaw="F"'
      WRITE (iw, '(T8,A)') 'paw_as_gipaw="F"'
      IF (pot%nlcc) THEN
         WRITE (iw, '(T8,A)') 'core_correction="T"'
      ELSE
         WRITE (iw, '(T8,A)') 'core_correction="F"'
      END IF
      WRITE (iw, '(T8,A)') 'functional="'//TRIM(ADJUSTL(xcfstr))//'"'
      CALL compose(string, "z_valence", rval=pot%zion)
      WRITE (iw, '(T8,A)') TRIM(string)
      CALL compose(string, "total_psenergy", rval=2._dp*atom%energy%etot)
      WRITE (iw, '(T8,A)') TRIM(string)
      WRITE (iw, '(T8,A)') 'wfc_cutoff="0.0E+00"'
      WRITE (iw, '(T8,A)') 'rho_cutoff="0.0E+00"'
      lmax = -1
      DO l = 0, lmat
         IF (pot%nl(l) > 0) lmax = l
      END DO
      CALL compose(string, "l_max", ival=lmax)
      WRITE (iw, '(T8,A)') TRIM(string)
      WRITE (iw, '(T8,A)') 'l_max_rho="0"'
      WRITE (iw, '(T8,A)') 'l_local="-3"'
      nr = atom%basis%grid%nr
      CALL compose(string, "mesh_size", ival=nr)
      WRITE (iw, '(T8,A)') TRIM(string)
      nwfc = SUM(atom%state%maxn_occ)
      CALL compose(string, "number_of_wfc", ival=nwfc)
      WRITE (iw, '(T8,A)') TRIM(string)
      IF (soc) THEN
         nbeta = pot%nl(0) + 2*SUM(pot%nl(1:))
      ELSE
         nbeta = SUM(pot%nl)
      END IF
      CALL compose(string, "number_of_proj", ival=nbeta)
      WRITE (iw, '(T8,A)') TRIM(string)//'/>'

      ! Mesh
      up = atom%basis%grid%rad(1) < atom%basis%grid%rad(nr)
      WRITE (iw, '(T4,A)') '<PP_MESH'
      WRITE (iw, '(T8,A)') 'dx="1.E+00"'
      CALL compose(string, "mesh", ival=nr)
      WRITE (iw, '(T8,A)') TRIM(string)
      WRITE (iw, '(T8,A)') 'xmin="1.E+00"'
      rmax = MAXVAL(atom%basis%grid%rad)
      CALL compose(string, "rmax", rval=rmax)
      WRITE (iw, '(T8,A)') TRIM(string)
      WRITE (iw, '(T8,A)') 'zmesh="1.E+00">'
      WRITE (iw, '(T8,A)') '<PP_R type="real"'
      CALL compose(string, "size", ival=nr)
      WRITE (iw, '(T12,A)') TRIM(string)
      WRITE (iw, '(T12,A)') 'columns="4">'
      IF (up) THEN
         WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=1, nr)
      ELSE
         WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=nr, 1, -1)
      END IF
      WRITE (iw, '(T8,A)') '</PP_R>'
      WRITE (iw, '(T8,A)') '<PP_RAB type="real"'
      CALL compose(string, "size", ival=nr)
      WRITE (iw, '(T12,A)') TRIM(string)
      WRITE (iw, '(T12,A)') 'columns="4">'
      IF (up) THEN
         WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=1, nr)
      ELSE
         WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=nr, 1, -1)
      END IF
      WRITE (iw, '(T8,A)') '</PP_RAB>'
      WRITE (iw, '(T4,A)') '</PP_MESH>'

      ! NLCC
      IF (pot%nlcc) THEN
         WRITE (iw, '(T4,A)') '<PP_NLCC type="real"'
         CALL compose(string, "size", ival=nr)
         WRITE (iw, '(T8,A)') TRIM(string)
         WRITE (iw, '(T8,A)') 'columns="4">'
         ALLOCATE (corden(nr))
         corden(:) = 0.0_dp
         CALL atom_core_density(corden, atom%potential, "RHO", atom%basis%grid%rad)
         IF (up) THEN
            WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=1, nr)
         ELSE
            WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=nr, 1, -1)
         END IF
         DEALLOCATE (corden)
         WRITE (iw, '(T4,A)') '</PP_NLCC>'
      END IF

      ! local PP
      WRITE (iw, '(T4,A)') '<PP_LOCAL type="real"'
      CALL compose(string, "size", ival=nr)
      WRITE (iw, '(T8,A)') TRIM(string)
      WRITE (iw, '(T8,A)') 'columns="4">'
      ALLOCATE (locpot(nr))
      locpot(:) = 0.0_dp
      CALL atom_local_potential(locpot, pot, atom%basis%grid%rad)
      IF (up) THEN
         WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=1, nr)
      ELSE
         WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=nr, 1, -1)
      END IF
      DEALLOCATE (locpot)
      WRITE (iw, '(T4,A)') '</PP_LOCAL>'

      ! nonlocal PP
      WRITE (iw, '(T4,A)') '<PP_NONLOCAL>'
      ALLOCATE (rp(nr), ef(nr), beta(nr))
      ibeta = 0
      DO l = 0, lmat
         IF (pot%nl(l) == 0) CYCLE
         rl = pot%rcnl(l)
         rp(:) = atom%basis%grid%rad
         ef(:) = EXP(-0.5_dp*rp*rp/(rl*rl))
         DO i = 1, pot%nl(l)
            pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
            j = l + 2*i - 1
            pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
            beta(:) = pf*rp**(l + 2*i - 2)*ef
            beta(:) = beta*rp
            !
            ibeta = ibeta + 1
            CALL compose(string, "<PP_BETA", counter=ibeta)
            WRITE (iw, '(T8,A)') TRIM(string)
            CALL compose(string, "angular_momentum", ival=l)
            WRITE (iw, '(T12,A)') TRIM(string)
            WRITE (iw, '(T12,A)') 'type="real"'
            CALL compose(string, "size", ival=nr)
            WRITE (iw, '(T12,A)') TRIM(string)
            WRITE (iw, '(T12,A)') 'columns="4">'
            IF (up) THEN
               WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
            ELSE
               WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
            END IF
            CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.TRUE.)
            WRITE (iw, '(T8,A)') TRIM(string)
            IF (soc .AND. l /= 0) THEN
               ibeta = ibeta + 1
               CALL compose(string, "<PP_BETA", counter=ibeta)
               WRITE (iw, '(T8,A)') TRIM(string)
               CALL compose(string, "angular_momentum", ival=l)
               WRITE (iw, '(T12,A)') TRIM(string)
               WRITE (iw, '(T12,A)') 'type="real"'
               CALL compose(string, "size", ival=nr)
               WRITE (iw, '(T12,A)') TRIM(string)
               WRITE (iw, '(T12,A)') 'columns="4">'
               IF (up) THEN
                  WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
               ELSE
                  WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
               END IF
               CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.TRUE.)
               WRITE (iw, '(T8,A)') TRIM(string)
            END IF
         END DO
      END DO
      DEALLOCATE (rp, ef, beta)
      ! nonlocal PP matrix elements
      ALLOCATE (dij(nbeta, nbeta))
      dij = 0._dp
      IF (soc) THEN
         ! from reverse engineering we guess: hnl, knl, hnl, knl, ...
         n0 = pot%nl(0)
         DO l = 0, lmat
            IF (pot%nl(l) == 0) CYCLE
            IF (l == 0) THEN
               ! no SO term for l=0
               dij(1:n0, 1:n0) = pot%hnl(1:n0, 1:n0, 0)
            ELSE
               ibeta = n0 + 2*SUM(pot%nl(1:l - 1))
               DO i = 1, pot%nl(l)
                  ii = 2*(i - 1) + 1
                  DO k = 1, pot%nl(l)
                     kk = 2*(k - 1) + 1
                     vnlp = pot%hnl(i, k, l) + 0.5_dp*l*pot%knl(i, k, l)
                     vnlm = pot%hnl(i, k, l) - 0.5_dp*(l + 1)*pot%knl(i, k, l)
                     dij(ibeta + ii, ibeta + kk) = vnlm
                     dij(ibeta + ii + 1, ibeta + kk + 1) = vnlp
                  END DO
               END DO
            END IF
         END DO
      ELSE
         DO l = 0, lmat
            IF (pot%nl(l) == 0) CYCLE
            ibeta = SUM(pot%nl(0:l - 1)) + 1
            i = ibeta + pot%nl(l) - 1
            dij(ibeta:i, ibeta:i) = pot%hnl(1:pot%nl(l), 1:pot%nl(l), l)
         END DO
      END IF
      WRITE (iw, '(T8,A)') '<PP_DIJ type="real"'
      WRITE (iw, '(T12,A)') 'columns="4">'
      WRITE (iw, '(T12,4ES25.12E3)') ((0.5_dp*dij(i, j), j=1, nbeta), i=1, nbeta)
      WRITE (iw, '(T8,A)') '</PP_DIJ>'
      DEALLOCATE (dij)
      WRITE (iw, '(T4,A)') '</PP_NONLOCAL>'

      ! atomic wavefunctions
      ALLOCATE (beta(nr))
      WRITE (iw, '(T4,A)') '<PP_PSWFC>'
      nwfn = 0
      DO l = 0, lmat
         DO i = 1, 10
            IF (ABS(atom%state%occupation(l, i)) == 0._dp) CYCLE
            occa = atom%state%occupation(l, i)
            occb = 0.0_dp
            IF (soc .AND. l /= 0) THEN
               occb = FLOOR(occa/2)
               occa = occa - occb
            END IF
            nwfn = nwfn + 1
            CALL compose(string, "<PP_CHI", counter=nwfn)
            WRITE (iw, '(T8,A)') TRIM(string)
            CALL compose(string, "l", ival=l)
            WRITE (iw, '(T12,A)') TRIM(string)
            CALL compose(string, "occupation", rval=occa)
            WRITE (iw, '(T12,A)') TRIM(string)
            CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
            WRITE (iw, '(T12,A)') TRIM(string)
            WRITE (iw, '(T12,A)') 'columns="4">'
            beta = 0._dp
            DO k = 1, atom%basis%nbas(l)
               beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
            END DO
            beta(:) = beta*atom%basis%grid%rad
            IF (up) THEN
               WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=1, nr)
            ELSE
               WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=nr, 1, -1)
            END IF
            CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.TRUE.)
            WRITE (iw, '(T8,A)') TRIM(string)
            IF (soc .AND. l /= 0) THEN
               nwfn = nwfn + 1
               CALL compose(string, "<PP_CHI", counter=nwfn)
               WRITE (iw, '(T8,A)') TRIM(string)
               CALL compose(string, "l", ival=l)
               WRITE (iw, '(T12,A)') TRIM(string)
               CALL compose(string, "occupation", rval=occb)
               WRITE (iw, '(T12,A)') TRIM(string)
               CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
               WRITE (iw, '(T12,A)') TRIM(string)
               WRITE (iw, '(T12,A)') 'columns="4">'
               beta = 0._dp
               DO k = 1, atom%basis%nbas(l)
                  beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
               END DO
               beta(:) = beta*atom%basis%grid%rad
               IF (up) THEN
                  WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=1, nr)
               ELSE
                  WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=nr, 1, -1)
               END IF
               CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.TRUE.)
               WRITE (iw, '(T8,A)') TRIM(string)
            END IF
         END DO
      END DO
      WRITE (iw, '(T4,A)') '</PP_PSWFC>'
      DEALLOCATE (beta)

      ! atomic charge
      ALLOCATE (dens(nr))
      WRITE (iw, '(T4,A)') '<PP_RHOATOM type="real"'
      CALL compose(string, "size", ival=nr)
      WRITE (iw, '(T8,A)') TRIM(string)
      WRITE (iw, '(T8,A)') 'columns="4">'
      CALL atom_density(dens, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
                        "RHO", atom%basis%grid%rad)
      IF (up) THEN
         WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
      ELSE
         WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
      END IF
      WRITE (iw, '(T4,A)') '</PP_RHOATOM>'

      ! Check atomic electron count
      rhoatom = SUM(fourpi*atom%basis%grid%wr(:)*dens(:))
      IF (ABS(pot%zion - rhoatom) > 0.01_dp) THEN
         CALL cp_warn(__LOCATION__, "Valence charge and electron count differ by more than 0.01")
      END IF

      DEALLOCATE (dens)
      ! PP SOC information
      IF (soc) THEN
         WRITE (iw, '(T4,A)') '<PP_SPIN_ORB>'
!        <PP_RELWFC.1 index="1" els="6S" nn="1" lchi="0" jchi="5.000000000000e-1" oc="1.000000000000e0"/>
         nwfn = 0
         DO l = 0, lmat
            DO i = 1, 10
               IF (ABS(atom%state%occupation(l, i)) == 0._dp) CYCLE
               occa = atom%state%occupation(l, i)
               occb = 0.0_dp
               IF (l /= 0) THEN
                  occb = FLOOR(occa/2)
                  occa = occa - occb
               END IF
               nwfn = nwfn + 1
               CALL compose(string, "<PP_RELWFC", counter=nwfn)
               WRITE (iw, '(T8,A)') TRIM(string)
               CALL compose(string, "index", ival=nwfn)
               WRITE (iw, '(T12,A)') TRIM(string)
               !!CALL compose(string, "els", cval=xxxx)
               !!WRITE (iw, '(T12,A)') TRIM(string)
               CALL compose(string, "nn", ival=nwfn)
               WRITE (iw, '(T12,A)') TRIM(string)
               CALL compose(string, "lchi", ival=l)
               WRITE (iw, '(T12,A)') TRIM(string)
               IF (l == 0) THEN
                  rlp = 0.5_dp
               ELSE
                  rlp = l - 0.5_dp
               END IF
               CALL compose(string, "jchi", rval=rlp)
               WRITE (iw, '(T12,A)') TRIM(string)
               CALL compose(string, "oc", rval=occa)
               WRITE (iw, '(T12,A)') TRIM(string)
               WRITE (iw, '(T12,A)') '/>'
               IF (l /= 0) THEN
                  nwfn = nwfn + 1
                  CALL compose(string, "<PP_RELWFC", counter=nwfn)
                  WRITE (iw, '(T8,A)') TRIM(string)
                  CALL compose(string, "index", ival=nwfn)
                  WRITE (iw, '(T12,A)') TRIM(string)
                  !!CALL compose(string, "els", cval=xxxx)
                  !!WRITE (iw, '(T12,A)') TRIM(string)
                  CALL compose(string, "nn", ival=nwfn)
                  WRITE (iw, '(T12,A)') TRIM(string)
                  CALL compose(string, "lchi", ival=l)
                  WRITE (iw, '(T12,A)') TRIM(string)
                  rlp = l + 0.5_dp
                  CALL compose(string, "jchi", rval=rlp)
                  WRITE (iw, '(T12,A)') TRIM(string)
                  CALL compose(string, "oc", rval=occa)
                  WRITE (iw, '(T12,A)') TRIM(string)
                  WRITE (iw, '(T12,A)') '/>'
               END IF
            END DO
         END DO
         ibeta = 0
         DO l = 0, lmat
            IF (pot%nl(l) == 0) CYCLE
            DO i = 1, pot%nl(l)
               ibeta = ibeta + 1
               CALL compose(string, "<PP_RELBETA", counter=ibeta)
               WRITE (iw, '(T8,A)') TRIM(string)
               CALL compose(string, "index", ival=ibeta)
               WRITE (iw, '(T12,A)') TRIM(string)
               CALL compose(string, "lll", ival=l)
               WRITE (iw, '(T12,A)') TRIM(string)
               IF (l == 0) THEN
                  rlp = 0.5_dp
               ELSE
                  rlp = l - 0.5_dp
               END IF
               CALL compose(string, "jjj", rval=rlp)
               WRITE (iw, '(T12,A)') TRIM(string)
               WRITE (iw, '(T12,A)') '/>'
               IF (l > 0) THEN
                  ibeta = ibeta + 1
                  CALL compose(string, "<PP_RELBETA", counter=ibeta)
                  WRITE (iw, '(T8,A)') TRIM(string)
                  CALL compose(string, "index", ival=ibeta)
                  WRITE (iw, '(T12,A)') TRIM(string)
                  CALL compose(string, "lll", ival=l)
                  WRITE (iw, '(T12,A)') TRIM(string)
                  rlp = l + 0.5_dp
                  CALL compose(string, "jjj", rval=rlp)
                  WRITE (iw, '(T12,A)') TRIM(string)
                  WRITE (iw, '(T12,A)') '/>'
               END IF
            END DO
         END DO
         WRITE (iw, '(T4,A)') '</PP_SPIN_ORB>'
      END IF

      WRITE (iw, '(A)') '</UPF>'

   END SUBROUTINE atom_write_upf

! **************************************************************************************************
!> \brief Produce an XML attribute string 'tag="counter/rval/ival/cval"'
!> \param string   composed string
!> \param tag      attribute tag
!> \param counter  counter
!> \param rval     real variable
!> \param ival     integer variable
!> \param cval     string variable
!> \param isfinal  close the current XML element if this is the last attribute
!> \par History
!>    * 09.2012 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE compose(string, tag, counter, rval, ival, cval, isfinal)
      CHARACTER(len=*), INTENT(OUT)                      :: string
      CHARACTER(len=*), INTENT(IN)                       :: tag
      INTEGER, INTENT(IN), OPTIONAL                      :: counter
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: rval
      INTEGER, INTENT(IN), OPTIONAL                      :: ival
      CHARACTER(len=*), INTENT(IN), OPTIONAL             :: cval
      LOGICAL, INTENT(IN), OPTIONAL                      :: isfinal

      CHARACTER(LEN=default_string_length)               :: str
      LOGICAL                                            :: fin

      IF (PRESENT(counter)) THEN
         WRITE (str, "(I12)") counter
      ELSEIF (PRESENT(rval)) THEN
         WRITE (str, "(G18.8)") rval
      ELSEIF (PRESENT(ival)) THEN
         WRITE (str, "(I12)") ival
      ELSEIF (PRESENT(cval)) THEN
         WRITE (str, "(A)") TRIM(ADJUSTL(cval))
      ELSE
         WRITE (str, "(A)") ""
      END IF
      fin = .FALSE.
      IF (PRESENT(isfinal)) fin = isfinal
      IF (PRESENT(counter)) THEN
         IF (fin) THEN
            WRITE (string, "(A,A1,A,A1)") TRIM(ADJUSTL(tag)), '.', TRIM(ADJUSTL(str)), '>'
         ELSE
            WRITE (string, "(A,A1,A)") TRIM(ADJUSTL(tag)), '.', TRIM(ADJUSTL(str))
         END IF
      ELSE
         IF (fin) THEN
            WRITE (string, "(A,A2,A,A2)") TRIM(ADJUSTL(tag)), '="', TRIM(ADJUSTL(str)), '>"'
         ELSE
            WRITE (string, "(A,A2,A,A1)") TRIM(ADJUSTL(tag)), '="', TRIM(ADJUSTL(str)), '"'
         END IF
      END IF
   END SUBROUTINE compose

END MODULE atom_energy