ps.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2012 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module ps_oct_m
  use atomic_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use lalg_adv_oct_m
  use loct_math_oct_m
  use parser_oct_m
  use logrid_oct_m
  use math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use ps_cpi_oct_m
  use ps_fhi_oct_m
  use ps_hgh_oct_m
  use ps_xml_oct_m
  use ps_in_grid_oct_m
#ifdef HAVE_PSPIO
  use fpspio_m
#endif
  use ps_psf_oct_m
  use pseudo_oct_m
  use sort_oct_m
  use splines_oct_m
  use spline_filter_oct_m
  implicit none
 
  private
  public ::                     &
    ps_t,                       &
    ps_init,                    &
    ps_pspio_init,              &
    ps_separate,                &
    ps_filter,                  &
    ps_getradius,               &
    ps_derivatives,             &
    ps_debug,                   &
    ps_niwfs,                   &
    ps_bound_niwfs,             &
    ps_end,                     &
    ps_has_density,             &
    ps_has_nlcc,                &
    ps_density_volume,          &
    ps_guess_atomic_occupations
 
  integer, parameter, public :: &
    PS_FILTER_NONE = 0,         &
    ps_filter_ts   = 2,         &
    ps_filter_bsb  = 3
 
  integer, public, parameter ::  &
    PROJ_NONE = 0,  &
    proj_hgh  = 1,  &
    proj_kb   = 2,  &
    proj_rkb  = 3
 
  integer, public, parameter :: &
    PROJ_J_INDEPENDENT  = 0, & !< Non-relativisitic or scalar-relativistic pseudopotentials
    proj_j_dependent    = 1, & 
    proj_j_average      = 2    
 
  integer, parameter, public :: INVALID_L = 333
 
  character(len=4), parameter  :: ps_name(PSEUDO_FORMAT_UPF1:PSEUDO_FORMAT_PSP8) = &
    (/"upf1", "upf2", "qso ", "psml", "psf ", "cpi ", "fhi ", "hgh ", "psp8"/)
 
  type ps_t
    ! Components are public by default
    integer :: projector_type
    integer :: relativistic_treatment
 
    character(len=10), private :: label
 
    integer, private  :: ispin
 
    real(real64), private    :: z
    real(real64)             :: z_val
    type(valconf_t)   :: conf
 
    type(logrid_t), private :: g
 
    type(spline_t), allocatable :: ur(:, :)
    type(spline_t), allocatable, private :: ur_sq(:, :)
    logical, allocatable    :: bound(:, :)
 
    ! Kleinman-Bylander projectors stuff
    integer  :: lmax
    integer  :: llocal
 
    type(spline_t) :: vl
    logical        :: no_vl = .false.  
 
    real(real64) :: projectors_sphere_threshold
    real(real64) :: rc_max
 
    integer  :: kbc
    integer  :: projectors_per_l(5)
    real(real64), allocatable :: h(:,:,:), k(:, :, :)
    type(spline_t), allocatable :: kb(:, :)
    type(spline_t), allocatable :: dkb(:, :)
 
    logical :: nlcc
    type(spline_t) :: core
    type(spline_t) :: core_der
 
 
    !LONG-RANGE PART OF THE LOCAL POTENTIAL
 
    logical, private :: has_long_range
    logical, private :: is_separated
 
    type(spline_t), private :: vlr
    type(spline_t) :: vlr_sq
    type(spline_t) :: nlr
 
    real(real64) :: sigma_erf
 
    logical,        private     :: has_density
    type(spline_t), allocatable :: density(:)
    type(spline_t), allocatable :: density_der(:)
 
    logical          :: local
    integer          :: file_format
    integer, private :: pseudo_type
    integer          :: exchange_functional
    integer          :: correlation_functional
  end type ps_t
 
  real(real64), parameter :: eps = 1.0e-8_real64 
 
contains
 
 
  ! ---------------------------------------------------------
  subroutine ps_init(ps, namespace, label, z, user_lmax, user_llocal, ispin, filename)
    type(ps_t),        intent(out)   :: ps
    type(namespace_t), intent(in)    :: namespace
    character(len=10), intent(in)    :: label
    integer,           intent(in)    :: user_lmax
    integer,           intent(in)    :: user_llocal
    integer,           intent(in)    :: ispin
    real(real64),      intent(in)    :: z
    character(len=*),  intent(in)    :: filename
 
    integer :: l, ii, ll, is, ierr
    type(ps_psf_t) :: ps_psf
    type(ps_cpi_t) :: ps_cpi
    type(ps_fhi_t) :: ps_fhi
    type(ps_hgh_t) :: ps_hgh
    type(ps_xml_t) :: ps_xml
    real(real64), allocatable :: eigen(:, :)
 
    push_sub(ps_init)
 
    ps%exchange_functional = pseudo_exchange_unknown
    ps%correlation_functional = pseudo_correlation_unknown
 
    ! Fix the threshold to calculate the radius of the projector-function localization spheres:
 
    call messages_obsolete_variable(namespace, 'SpecieProjectorSphereThreshold', 'SpeciesProjectorSphereThreshold')
 
    !%Variable SpeciesProjectorSphereThreshold
    !%Type float
    !%Default 0.001
    !%Section System::Species
    !%Description
    !% The pseudopotentials may be composed of a local part, and a linear combination of nonlocal
    !% operators. These nonlocal projectors have "projector" form, <math> \left| v \right> \left< v \right| </math>
    !% (or, more generally speaking, <math> \left| u \right> \left< v \right| </math>).
    !% These projectors are localized in real space -- that is, the function <math>v</math>
    !% has a finite support around the nucleus. This region where the projectors are localized should
    !% be small or else the computation time required to operate with them will be very large.
    !%
    !% In practice, this localization is fixed by requiring the definition of the projectors to be
    !% contained in a sphere of a certain radius. This radius is computed by making sure that the
    !% absolute value of the projector functions, at points outside the localization sphere, is
    !% below a certain threshold. This threshold is set by <tt>SpeciesProjectorSphereThreshold</tt>.
    !%End
    call parse_variable(namespace, 'SpeciesProjectorSphereThreshold', 0.001_real64, ps%projectors_sphere_threshold)
    if (ps%projectors_sphere_threshold <= m_zero) call messages_input_error(namespace, 'SpeciesProjectorSphereThreshold')
 
    ps%file_format = pseudo_detect_format(filename)
 
    if (ps%file_format == pseudo_format_file_not_found) then
      call messages_write("Cannot open pseudopotential file '"//trim(filename)//"'.")
      call messages_fatal(namespace=namespace)
    end if
 
    if (ps%file_format == pseudo_format_unknown) then
      call messages_write("Cannot determine the pseudopotential type for species '"//trim(label)//"' from", &
        new_line = .true.)
      call messages_write("file '"//trim(filename)//"'.")
      call messages_fatal(namespace=namespace)
    end if
 
    ps%label   = label
    ps%ispin   = ispin
    ps%relativistic_treatment  = proj_j_independent
    ps%projector_type = proj_kb
    ps%sigma_erf = 0.625_real64 ! This is hard-coded to a reasonable value
 
    ps%projectors_per_l = 0
 
    select case (ps%file_format)
    case (pseudo_format_psf, pseudo_format_hgh)
      ps%has_density = .true.
    case default
      ps%has_density = .false.
    end select
 
    select case (ps%file_format)
    case (pseudo_format_psf)
      ps%pseudo_type   = pseudo_type_semilocal
 
      call ps_psf_init(ps_psf, ispin, filename, namespace)
 
      call valconf_copy(ps%conf, ps_psf%conf)
      ps%z      = z
      ps%conf%z = nint(z) ! atomic number
      ps%kbc    = 1     ! only one projector per angular momentum
 
      ps%lmax = ps_psf%ps_grid%no_l_channels - 1
 
      if (user_lmax /= invalid_l) then
        ps%lmax = min(ps%lmax, user_lmax) ! Maybe the file does not have enough components.
        if (user_lmax /= ps%lmax) then
          message(1) = "lmax in Species block for " // trim(label) // &
            " is larger than number available in pseudopotential."
          call messages_fatal(1, namespace=namespace)
        end if
      end if
 
      ps%conf%p = ps_psf%ps_grid%no_l_channels
      if (ps%lmax == 0) ps%llocal = 0 ! Vanderbilt is not acceptable if ps%lmax == 0.
 
      ! the local part of the pseudo
      if (user_llocal == invalid_l) then
        ps%llocal = 0
      else
        ps%llocal = user_llocal
      end if
 
      ps%projectors_per_l(1:ps%lmax+1) = 1
      if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
 
      call ps_psf_process(ps_psf, namespace, ps%lmax, ps%llocal)
      call logrid_copy(ps_psf%ps_grid%g, ps%g)
 
    case (pseudo_format_cpi, pseudo_format_fhi)
      ps%pseudo_type   = pseudo_type_semilocal
 
      if (ps%file_format == pseudo_format_cpi) then
        call ps_cpi_init(ps_cpi, trim(filename), namespace)
        ps%conf%p      = ps_cpi%ps_grid%no_l_channels
      else
        call ps_fhi_init(ps_fhi, trim(filename), namespace)
        ps%conf%p      = ps_fhi%ps_grid%no_l_channels
      end if
 
      ps%conf%z      = nint(z)
      ps%conf%symbol = label(1:2)
      ps%conf%type   = 1
      do l = 1, ps%conf%p
        ps%conf%l(l) = l - 1
      end do
 
      ps%z      = z
      ps%kbc    = 1     ! only one projector per angular momentum
 
      ps%lmax  = ps%conf%p - 1
 
      if (user_lmax /= invalid_l) then
        ps%lmax = min(ps%lmax, user_lmax) ! Maybe the file does not have enough components.
        if (user_lmax /= ps%lmax) then
          message(1) = "lmax in Species block for " // trim(label) // &
            " is larger than number available in pseudopotential."
          call messages_fatal(1, namespace=namespace)
        end if
      end if
 
      if (ps%lmax == 0) ps%llocal = 0 ! Vanderbilt is not acceptable if ps%lmax == 0.
 
      ! the local part of the pseudo
      if (user_llocal == invalid_l) then
        ps%llocal = 0
      else
        ps%llocal = user_llocal
      end if
 
      ps%projectors_per_l(1:ps%lmax+1) = 1
      if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
 
      if (ps%file_format == pseudo_format_cpi) then
        call ps_cpi_process(ps_cpi, ps%llocal, namespace)
        call logrid_copy(ps_cpi%ps_grid%g, ps%g)
      else
        call ps_fhi_process(ps_fhi, ps%lmax, ps%llocal, namespace)
        call logrid_copy(ps_fhi%ps_grid%g, ps%g)
      end if
 
    case (pseudo_format_hgh)
      ps%pseudo_type   = pseudo_type_kleinman_bylander
      ps%projector_type = proj_hgh
 
      call hgh_init(ps_hgh, trim(filename), namespace)
      call valconf_copy(ps%conf, ps_hgh%conf)
 
      ps%z        = z
      ps%z_val    = ps_hgh%z_val
      ps%kbc      = 3
      ps%llocal   = -1
      ps%lmax     = ps_hgh%l_max
      ps%conf%symbol = label(1:2)
      ps%sigma_erf = ps_hgh%rlocal ! We use the correct value
 
      ps%projectors_per_l(1:ps%lmax+1) = 1
 
      ! Get the occupations from the valence charge of the atom
      ps%conf%occ = m_zero
      ! We impose here non-spin-polarized occupations, to preserve the behavior of the code
      ! We might want to change this to get a better LCAO guess
      call ps_guess_atomic_occupations(namespace, ps%z, ps%z_val, 1, ps%conf)
      ! We need the information to solve the Schrodinder equation
      call valconf_copy(ps_hgh%conf, ps%conf)
 
      call hgh_process(ps_hgh, namespace)
      call logrid_copy(ps_hgh%g, ps%g)
 
      ! In case of spin-polarized calculations, we properly distribute the electrons
      if (ispin == 2) then
        call valconf_unpolarized_to_polarized(ps%conf)
      end if
 
    case (pseudo_format_qso, pseudo_format_upf1, pseudo_format_upf2, pseudo_format_psml, pseudo_format_psp8)
 
      call ps_xml_init(ps_xml, namespace, trim(filename), ps%file_format, ierr)
 
      ps%pseudo_type   = pseudo_type(ps_xml%pseudo)
      ps%exchange_functional = pseudo_exchange(ps_xml%pseudo)
      ps%correlation_functional = pseudo_correlation(ps_xml%pseudo)
 
      ps%z      = z
      ps%conf%z = nint(z)
 
      if (ps_xml%kleinman_bylander) then
        ps%conf%p = ps_xml%nwavefunctions
      else
        ps%conf%p = ps_xml%lmax + 1
      end if
 
      do ll = 0, ps_xml%lmax
        ps%conf%l(ll + 1) = ll
      end do
 
      ps%kbc   = ps_xml%nchannels
      ps%lmax  = ps_xml%lmax
 
      ps%projectors_per_l = 0
      do ll = 0, ps_xml%lmax
        ps%projectors_per_l(ll+1) = pseudo_nprojectors_per_l(ps_xml%pseudo, ll)
      end do
 
      if (ps_xml%kleinman_bylander) then
        ps%llocal = ps_xml%llocal
      else
        ! we have several options
        ps%llocal = 0                                     ! the default
        if (ps_xml%llocal >= 0) ps%llocal = ps_xml%llocal  ! the one given in the pseudopotential file
        if (user_llocal /= invalid_l) ps%llocal = user_llocal ! user supplied local component
        assert(ps%llocal >= 0)
        assert(ps%llocal <= ps%lmax)
      end if
 
      ps%g%nrval = ps_xml%grid_size
 
      safe_allocate(ps%g%rofi(1:ps%g%nrval))
      safe_allocate(ps%g%r2ofi(1:ps%g%nrval))
 
      do ii = 1, ps%g%nrval
        ps%g%rofi(ii) = ps_xml%grid(ii)
        ps%g%r2ofi(ii) = ps_xml%grid(ii)**2
      end do
 
    end select
 
    ps%local = (ps%lmax == 0 .and. ps%llocal == 0 ) .or. (ps%lmax == -1 .and. ps%llocal == -1)
 
    ! We allocate all the stuff
    safe_allocate(ps%kb   (0:ps%lmax, 1:ps%kbc))
    safe_allocate(ps%dkb  (0:ps%lmax, 1:ps%kbc))
    safe_allocate(ps%ur   (1:ps%conf%p, 1:ps%ispin))
    safe_allocate(ps%ur_sq(1:ps%conf%p, 1:ps%ispin))
    safe_allocate(ps%bound(1:ps%conf%p, 1:ps%ispin))
    safe_allocate(ps%h    (0:ps%lmax, 1:ps%kbc, 1:ps%kbc))
    safe_allocate(ps%density(1:ps%ispin))
    safe_allocate(ps%density_der(1:ps%ispin))
 
    call spline_init(ps%kb)
    call spline_init(ps%dkb)
    call spline_init(ps%vl)
    call spline_init(ps%core)
    call spline_init(ps%core_der)
    call spline_init(ps%density)
    call spline_init(ps%density_der)
 
    safe_allocate(eigen(1:ps%conf%p, 1:ps%ispin))
    eigen = m_zero
 
    ! Now we load the necessary information.
    select case (ps%file_format)
    case (pseudo_format_psf)
      call ps_psf_get_eigen(ps_psf, eigen)
      call ps_grid_load(ps, ps_psf%ps_grid)
      call ps_psf_end(ps_psf)
    case (pseudo_format_cpi)
      call ps_grid_load(ps, ps_cpi%ps_grid)
      call ps_cpi_end(ps_cpi)
    case (pseudo_format_fhi)
      call ps_grid_load(ps, ps_fhi%ps_grid)
      call ps_fhi_end(ps_fhi)
    case (pseudo_format_hgh)
      call hgh_get_eigen(ps_hgh, eigen)
      safe_allocate(ps%k    (0:ps%lmax, 1:ps%kbc, 1:ps%kbc))
      call hgh_load(ps, ps_hgh)
      call hgh_end(ps_hgh)
    case (pseudo_format_qso, pseudo_format_upf1, pseudo_format_upf2, pseudo_format_psml, pseudo_format_psp8)
      call ps_xml_load(ps, ps_xml, namespace)
      call ps_xml_end(ps_xml)
    end select
 
    if (ps_has_density(ps)) then
      do is = 1, ps%ispin
        call spline_der(ps%density(is), ps%density_der(is), ps%projectors_sphere_threshold)
      end do
    end if
 
    if (ps_has_nlcc(ps)) then
      call spline_der(ps%core, ps%core_der, ps%projectors_sphere_threshold)
    end if
 
    call ps_check_bound(ps, eigen)
 
    ps%has_long_range = .true.
    ps%is_separated = .false.
 
    call ps_info(ps, filename, namespace)
 
    safe_deallocate_a(eigen)
 
    pop_sub(ps_init)
  end subroutine ps_init
 
  !------------------------------------------------------------------------
 
  subroutine ps_info(ps, filename, namespace)
    type(ps_t),        intent(in) :: ps
    character(len=*),  intent(in) :: filename
    type(namespace_t), intent(in) :: namespace
 
    call messages_write("  Species '"//trim(ps%label)//"'", new_line = .true.)
    call messages_write("    type             : pseudopotential", new_line = .true.)
    call messages_write("    file             : '"//trim(filename)//"'")
    call messages_info(namespace=namespace)
 
    call messages_write("    file format      :")
    select case (ps%file_format)
    case (pseudo_format_upf1)
      call messages_write(" UPF1")
    case (pseudo_format_upf2)
      call messages_write(" UPF2")
    case (pseudo_format_qso)
      call messages_write(" QSO")
    case (pseudo_format_psml)
      call messages_write(" PSML")
    case (pseudo_format_psp8)
      call messages_write(" PSP8")
    case (pseudo_format_psf)
      call messages_write(" PSF")
    case (pseudo_format_cpi)
      call messages_write(" CPI")
    case (pseudo_format_fhi)
      call messages_write(" FHI")
    case (pseudo_format_hgh)
      call messages_write(" HGH")
    end select
    call messages_new_line()
 
    call messages_write("    valence charge   :")
    call messages_write(ps%z_val, align_left = .true., fmt = '(f4.1)')
    call messages_info(namespace=namespace)
 
    call messages_write("    atomic number    :")
    call messages_write(nint(ps%z), fmt = '(i4)')
    call messages_info(namespace=namespace)
 
    call messages_write("    form on file     :")
    select case (ps%pseudo_type)
    case (pseudo_type_ultrasoft)
      call messages_write(" ultrasoft")
    case (pseudo_type_semilocal)
      call messages_write(" semilocal")
    case (pseudo_type_kleinman_bylander)
      call messages_write(" kleinman-bylander")
    case (pseudo_type_paw)
      call messages_write(" paw")
    end select
    call messages_info(namespace=namespace)
 
    if (ps%pseudo_type == pseudo_type_semilocal) then
      call messages_write("    orbital origin   :")
      select case (ps%file_format)
      case (pseudo_format_psf)
        call messages_write(" calculated")
      case default
        call messages_write(" from file")
      end select
      call messages_info(namespace=namespace)
    end if
 
    call messages_write("    lmax             :")
    call messages_write(ps%lmax, fmt = '(i2)')
    call messages_info(namespace=namespace)
 
    call messages_write("    llocal           :")
    if (ps%llocal >= 0) then
      call messages_write(ps%llocal, fmt = '(i2)')
    else
      call messages_write(ps%llocal, fmt = '(i3)')
    end if
    call messages_info(namespace=namespace)
 
    call messages_write("    projectors per l :")
    call messages_write(ps%kbc, fmt = '(i2)')
    call messages_info(namespace=namespace)
 
    call messages_write("    total projectors :")
    if (ps%llocal < 0) then
      call messages_write(ps%kbc*(ps%lmax + 1), fmt = '(i2)')
    else
      call messages_write(ps%kbc*ps%lmax, fmt = '(i2)')
    end if
    call messages_info(namespace=namespace)
 
    if (ps%local) then
      call messages_write("    application form : local")
    else
      call messages_write("    application form : kleinman-bylander")
    end if
    call messages_info(namespace=namespace)
 
    call messages_write("    orbitals         :")
    call messages_write(ps_niwfs(ps), fmt='(i3)')
    call messages_info(namespace=namespace)
    call messages_write("    bound orbitals   :")
    call messages_write(ps_bound_niwfs(ps), fmt='(i3)')
    call messages_info(namespace=namespace)
 
    call messages_info(namespace=namespace)
 
  end subroutine ps_info
 
 
  ! ---------------------------------------------------------
  subroutine ps_separate(ps)
    type(ps_t),        intent(inout) :: ps
 
    real(real64), allocatable :: vsr(:), vlr(:), nlr(:)
    real(real64) :: r, exp_arg, arg, max_val_vsr
    integer :: ii
 
    push_sub(ps_separate)
 
    assert(ps%g%nrval > 0)
 
    safe_allocate(vsr(1:ps%g%nrval))
    safe_allocate(vlr(1:ps%g%nrval))
    safe_allocate(nlr(1:ps%g%nrval))
 
    ps%no_vl = .false.
 
    max_val_vsr = m_zero
 
    do ii = 1, ps%g%nrval
      r = ps%g%rofi(ii)
      arg = r/(ps%sigma_erf*sqrt(m_two))
      ! Note that the threshold is taken such that the Taylor expansion produces an error smaller than \epsilon
      ! given a third order Taylor expansion
      if(arg > 1e-3_real64) then
        vlr(ii)  = -ps%z_val*loct_erf(arg)/r
      else
        vlr(ii)  = -ps%z_val*m_two/(ps%sigma_erf*sqrt(m_two*m_pi)) * (m_one - arg**2 / m_three)
      end if
 
      vsr(ii) = spline_eval(ps%vl, r) - vlr(ii)
      if(abs(vsr(ii)) < 1.0e-14_real64) vsr(ii) = m_zero
      max_val_vsr = max(max_val_vsr, abs(vsr(ii)))
 
      exp_arg = -m_half*r**2/ps%sigma_erf**2
      if (exp_arg > m_min_exp_arg) then
        nlr(ii) = -ps%z_val/(ps%sigma_erf*sqrt(m_two*m_pi))**3*exp(exp_arg)
      else
        nlr(ii) = m_zero
      end if
    end do
 
    call spline_init(ps%vlr)
    call spline_fit(ps%g%nrval, ps%g%rofi, vlr, ps%vlr, ps%projectors_sphere_threshold)
 
    call spline_init(ps%vlr_sq)
    call spline_fit(ps%g%nrval, ps%g%r2ofi, vlr, ps%vlr_sq, ps%projectors_sphere_threshold)
 
    call spline_init(ps%nlr)
    call spline_fit(ps%g%nrval, ps%g%rofi, nlr, ps%nlr, ps%projectors_sphere_threshold)
 
    !overwrite vl
    call spline_end(ps%vl)
    call spline_init(ps%vl)
    call spline_fit(ps%g%nrval, ps%g%rofi, vsr, ps%vl, ps%projectors_sphere_threshold)
    if (max_val_vsr < 1.0e-12_real64) ps%no_vl = .true.
 
    safe_deallocate_a(vsr)
    safe_deallocate_a(vlr)
    safe_deallocate_a(nlr)
 
    ps%is_separated = .true.
 
    pop_sub(ps_separate)
  end subroutine ps_separate
 
 
  ! ---------------------------------------------------------
  subroutine ps_getradius(ps)
    type(ps_t), intent(inout) :: ps
    integer :: l, j
 
    push_sub(ps_getradius)
 
    ps%rc_max = m_zero
 
    do l = 0, ps%lmax
      if (l == ps%llocal) cycle
      do j = 1, ps%kbc
        ps%rc_max = max(ps%rc_max, ps%kb(l, j)%x_threshold)
      end do
    end do
 
    pop_sub(ps_getradius)
  end subroutine ps_getradius
 
 
  ! ---------------------------------------------------------
  subroutine ps_derivatives(ps)
    type(ps_t), intent(inout) :: ps
    integer :: l, j
 
    push_sub(ps_derivatives)
 
    do l = 0, ps%lmax
      do j = 1, ps%kbc
        call spline_der(ps%kb(l, j), ps%dkb(l, j), ps%projectors_sphere_threshold)
      end do
    end do
 
    pop_sub(ps_derivatives)
  end subroutine ps_derivatives
 
 
  ! ---------------------------------------------------------
  subroutine ps_filter(ps, filter, gmax)
    type(ps_t), intent(inout) :: ps
    integer,    intent(in)    :: filter
    real(real64),      intent(in)    :: gmax
 
    integer :: l, k, ispin
 
    real(real64) :: alpha, beta_fs, rmax, rcut, gamma, beta_rs
 
    push_sub(ps_filter)
    call profiling_in("PS_FILTER")
 
    select case (filter)
    case (ps_filter_none)
 
    case (ps_filter_ts)
      alpha = 1.1_real64
      gamma = m_two
 
      if(.not. ps%no_vl) then
        rmax = ps%vl%x_threshold
        ! For the special case l=-1, we use l=0 when doing the filtering, but it is not clear this is correct
        call spline_filter_mask(ps%vl, max(0, ps%llocal), rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
      end if
 
      do l = 0, ps%lmax
        if (l == ps%llocal) cycle
        do k = 1, ps%kbc
          call spline_filter_mask(ps%kb(l, k), l, ps%rc_max, gmax, alpha, gamma, ps%projectors_sphere_threshold)
        end do
      end do
 
      if (ps_has_nlcc(ps)) then
        rmax = ps%core%x_threshold
        call spline_filter_mask(ps%core, 0, rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
      end if
 
      if (ps_has_density(ps)) then
        do ispin = 1, ps%ispin
          if (abs(spline_integral(ps%density(ispin))) > 1.0e-12_real64) then
            rmax = ps%density(ispin)%x_threshold
            call spline_filter_mask(ps%density(ispin), 0, rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
            call spline_force_pos(ps%density(ispin), ps%projectors_sphere_threshold)
          end if
 
          if (abs(spline_integral(ps%density_der(ispin))) > 1.0e-12_real64) then
            rmax = ps%density_der(ispin)%x_threshold
            call spline_filter_mask(ps%density_der(ispin), 0, rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
          end if
        end do
      end if
 
    case (ps_filter_bsb)
      alpha   = 0.7_real64 ! The original was M_FOUR/7.0_real64
      beta_fs = 18.0_real64
      rcut    = 2.5_real64
      beta_rs = 0.4_real64
 
      ! For the special case l=-1, we use l=0 when doing the filtering, but it is not clear this is correct
      call spline_filter_bessel(ps%vl, max(0, ps%llocal), gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
      do l = 0, ps%lmax
        if (l == ps%llocal) cycle
        do k = 1, ps%kbc
          call spline_filter_bessel(ps%kb(l, k), l, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
        end do
      end do
 
      if (ps_has_nlcc(ps)) then
        call spline_filter_bessel(ps%core, 0, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
      end if
 
      if (ps_has_density(ps)) then
        do ispin = 1, ps%ispin
          call spline_filter_bessel(ps%density(ispin), 0, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
          call spline_force_pos(ps%density(ispin), ps%projectors_sphere_threshold)
          call spline_filter_bessel(ps%density_der(ispin), 0, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
        end do
      end if
 
    end select
 
    call profiling_out("PS_FILTER")
    pop_sub(ps_filter)
  end subroutine ps_filter
 
  ! ---------------------------------------------------------
  subroutine ps_check_bound(ps, eigen)
    type(ps_t), intent(inout) :: ps
    real(real64),      intent(in)    :: eigen(:,:)
 
    integer :: i, is, ir
    real(real64) :: ur1, ur2
 
    push_sub(ps_check_bound)
 
    ! Unbound states have positive eigenvalues
    where(eigen > m_zero)
      ps%bound = .false.
    elsewhere
      ps%bound = .true.
    end where
 
    ! We might not have information about the eigenvalues, so we need to check the wavefunctions
    do i = 1, ps%conf%p
      do is = 1, ps%ispin
        if (.not. ps%bound(i, is)) cycle
 
        do ir = ps%g%nrval, 3, -1
          ! First we look for the outmost value that is not zero
          if (abs(spline_eval(ps%ur(i, is), ps%g%rofi(ir))*ps%g%rofi(ir)) > m_zero) then
            ! Usually bound states have exponentially decaying wavefunctions,
            ! while unbound states have exponentially diverging
            ! wavefunctions. Therefore we check if the wavefunctions
            ! value is increasing with increasing radius. The fact
            ! that we do not use the wavefunctions outmost value that
            ! is not zero is on purpose, as some pseudopotential
            ! generators do funny things with that point.
            ur1 = spline_eval(ps%ur(i, is), ps%g%rofi(ir-2))*ps%g%rofi(ir-2)
            ur2 = spline_eval(ps%ur(i, is), ps%g%rofi(ir-1))*ps%g%rofi(ir-1)
            if ((ur1*ur2 > m_zero) .and. (abs(ur2) > abs(ur1))) ps%bound(i, is) = .false.
            exit
          end if
        end do
      end do
    end do
 
    pop_sub(ps_check_bound)
  end subroutine ps_check_bound
 
 
  ! ---------------------------------------------------------
  subroutine ps_debug(ps, dir, namespace, gmax)
    type(ps_t),        intent(in) :: ps
    character(len=*),  intent(in) :: dir
    type(namespace_t), intent(in) :: namespace
    real(real64),      intent(in) :: gmax
 
    ! We will plot also some Fourier transforms.
    type(spline_t), allocatable :: fw(:, :)
 
    integer  :: iunit
    integer  :: j, k, l
 
    push_sub(ps_debug)
 
    ! A text file with some basic data.
    iunit = io_open(trim(dir)//'/pseudo-info', namespace, action='write')
    write(iunit,'(a,/)')      ps%label
    write(iunit,'(a,a,/)')    'Format  : ', ps_name(ps%file_format)
    write(iunit,'(a,f6.3)')   'z       : ', ps%z
    write(iunit,'(a,f6.3,/)') 'zval    : ', ps%z_val
    write(iunit,'(a,i4)')     'lmax    : ', ps%lmax
    write(iunit,'(a,i4)')     'lloc    : ', ps%llocal
    write(iunit,'(a,i4,/)')   'kbc     : ', ps%kbc
    write(iunit,'(a,f9.5,/)') 'rcmax   : ', ps%rc_max
    write(iunit,'(a,/)')    'h matrix:'
    do l = 0, ps%lmax
      do k = 1, ps%kbc
        write(iunit,'(10f9.5)') (ps%h(l, k, j), j = 1, ps%kbc)
      end do
    end do
    if (allocated(ps%k)) then
      write(iunit,'(/,a,/)')    'k matrix:'
      do l = 0, ps%lmax
        do k = 1, ps%kbc
          write(iunit,'(10f9.5)') (ps%k(l, k, j), j = 1, ps%kbc)
        end do
      end do
    end if
 
    write(iunit,'(/,a)')    'orbitals:'
    do j = 1, ps%conf%p
      if (ps%ispin == 2) then
        write(iunit,'(1x,a,i2,3x,a,i2,3x,a,f5.1,3x,a,l1,3x,a,f5.1,3x,a,f5.1)') 'n = ', ps%conf%n(j), 'l = ', ps%conf%l(j), &
          'j = ', ps%conf%j(j), 'bound = ', all(ps%bound(j,:)), &
          'occ(1) = ', ps%conf%occ(j, 1), 'occ(2) = ', ps%conf%occ(j, 2)
      else
        write(iunit,'(1x,a,i2,3x,a,i2,3x,a,f5.1,3x,a,l1,3x,a,f5.1)') 'n = ', ps%conf%n(j), 'l = ', ps%conf%l(j), &
          'j = ', ps%conf%j(j), 'bound = ', all(ps%bound(j,:)), &
          'occ = ', ps%conf%occ(j, 1)
      end if
    end do
 
 
    call io_close(iunit)
 
    ! Local part of the pseudopotential
    iunit  = io_open(trim(dir)//'/local', namespace, action='write')
    call spline_print(ps%vl, iunit)
    call io_close(iunit)
 
    ! Local part of the pseudopotential
    iunit  = io_open(trim(dir)//'/local_long_range', namespace, action='write')
    call spline_print(ps%vlr, iunit)
    call io_close(iunit)
 
    ! Local part of the pseudopotential
    iunit  = io_open(trim(dir)//'/local_long_range_density', namespace, action='write')
    call spline_print(ps%nlr, iunit)
    call io_close(iunit)
 
    ! Fourier transform of the local part
    iunit = io_open(trim(dir)//'/local_ft', namespace, action='write')
    safe_allocate(fw(1, 1))
    call spline_init(fw(1, 1))
    call spline_3dft(ps%vl, fw(1, 1), ps%projectors_sphere_threshold, gmax = gmax)
    call spline_print(fw(1, 1), iunit)
    call spline_end(fw(1, 1))
    safe_deallocate_a(fw)
    call io_close(iunit)
 
    ! Kleinman-Bylander projectors
    iunit = io_open(trim(dir)//'/nonlocal', namespace, action='write')
    call spline_print(ps%kb, iunit)
    call io_close(iunit)
 
    iunit = io_open(trim(dir)//'/nonlocal_derivative', namespace, action='write')
    call spline_print(ps%dkb, iunit)
    call io_close(iunit)
 
    iunit = io_open(trim(dir)//'/nonlocal_ft', namespace, action='write')
    safe_allocate(fw(0:ps%lmax, 1:ps%kbc))
    call spline_init(fw)
    do k = 0, ps%lmax
      do j = 1, ps%kbc
        call spline_besselft(ps%kb(k, j), fw(k, j), k, threshold=ps%projectors_sphere_threshold, &
          gmax=m_four*gmax)
      end do
    end do
    call spline_print(fw, iunit)
    call spline_end(fw)
    safe_deallocate_a(fw)
    call io_close(iunit)
 
    ! Pseudo-wavefunctions
    iunit = io_open(trim(dir)//'/wavefunctions', namespace, action='write')
    call spline_print(ps%ur, iunit)
    call io_close(iunit)
 
    ! Density
    if (ps%has_density) then
      iunit = io_open(trim(dir)//'/density', namespace, action='write')
      call spline_print(ps%density, iunit)
      call io_close(iunit)
 
      iunit = io_open(trim(dir)//'/density_derivative', namespace, action='write')
      call spline_print(ps%density_der, iunit)
      call io_close(iunit)
    end if
 
    ! Non-linear core-corrections
    if (ps_has_nlcc(ps)) then
      iunit = io_open(trim(dir)//'/nlcc', namespace, action='write')
      call spline_print(ps%core, iunit)
      call io_close(iunit)
    end if
 
    pop_sub(ps_debug)
  end subroutine ps_debug
 
 
  ! ---------------------------------------------------------
  subroutine ps_end(ps)
    type(ps_t), intent(inout) :: ps
 
    if (.not. allocated(ps%kb)) return
 
    push_sub(ps_end)
 
    if (ps%is_separated) then
      call spline_end(ps%vlr)
      call spline_end(ps%vlr_sq)
      call spline_end(ps%nlr)
    end if
 
    call spline_end(ps%kb)
    call spline_end(ps%dkb)
    call spline_end(ps%ur)
    call spline_end(ps%ur_sq)
 
    call spline_end(ps%vl)
    call spline_end(ps%core)
    call spline_end(ps%core_der)
 
    if (allocated(ps%density)) call spline_end(ps%density)
    if (allocated(ps%density_der)) call spline_end(ps%density_der)
 
    call logrid_end(ps%g)
 
    safe_deallocate_a(ps%kb)
    safe_deallocate_a(ps%dkb)
    safe_deallocate_a(ps%ur)
    safe_deallocate_a(ps%ur_sq)
    safe_deallocate_a(ps%bound)
    safe_deallocate_a(ps%h)
    safe_deallocate_a(ps%k)
    safe_deallocate_a(ps%density)
    safe_deallocate_a(ps%density_der)
 
    pop_sub(ps_end)
  end subroutine ps_end
 
 
  ! ---------------------------------------------------------
  subroutine hgh_load(ps, ps_hgh)
    type(ps_t),     intent(inout) :: ps
    type(ps_hgh_t), intent(inout) :: ps_hgh
 
    push_sub(hgh_load)
 
    ! Fixes some components of ps
    ps%nlcc = .false.
    if (ps%lmax >= 0) then
      ps%rc_max = 1.1_real64 * maxval(ps_hgh%kbr(0:ps%lmax)) ! Increase a little.
    else
      ps%rc_max = m_zero
    end if
    ps%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc)
    ps%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc)
 
    ! now we fit the splines
    call get_splines()
 
    pop_sub(hgh_load)
 
  contains
 
    ! ---------------------------------------------------------
    subroutine get_splines()
      integer :: l, is, j, ip
      real(real64), allocatable :: hato(:), dens(:)
 
      push_sub(hgh_load.get_splines)
 
      safe_allocate(hato(1:ps_hgh%g%nrval))
      safe_allocate(dens(1:ps_hgh%g%nrval))
 
      ! Interpolate the KB-projection functions
      do l = 0, ps_hgh%l_max
        do j = 1, 3
          hato = ps_hgh%kb(:, l, j)
          call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%kb(l, j), ps%projectors_sphere_threshold)
        end do
      end do
 
      ! Now the part corresponding to the local pseudopotential
      ! where the asymptotic part is subtracted
      call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, ps_hgh%vlocal, ps%vl, ps%projectors_sphere_threshold)
 
      ! Define the table for the pseudo-wavefunction components (using splines)
      ! with a correct normalization function
      do is = 1, ps%ispin
        dens = m_zero
        do l = 1, ps%conf%p
          hato(2:ps_hgh%g%nrval) = ps_hgh%rphi(2:ps_hgh%g%nrval, l)/ps_hgh%g%rofi(2:ps_hgh%g%nrval)
          hato(1) = hato(2)
 
          do ip = 1, ps_hgh%g%nrval
            dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(m_four*m_pi)
          end do
 
          call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%ur(l, is), ps%projectors_sphere_threshold)
          call spline_fit(ps_hgh%g%nrval, ps_hgh%g%r2ofi, hato, ps%ur_sq(l, is), ps%projectors_sphere_threshold)
        end do
        call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
      end do
 
      safe_deallocate_a(hato)
      safe_deallocate_a(dens)
 
      pop_sub(hgh_load.get_splines)
    end subroutine get_splines
  end subroutine hgh_load
 
 
  ! ---------------------------------------------------------
  subroutine ps_grid_load(ps, ps_grid)
    type(ps_t),         intent(inout) :: ps
    type(ps_in_grid_t), intent(in)  :: ps_grid
 
    push_sub(ps_grid_load)
 
    ! Fixes some components of ps, read in ps_grid
    ps%z_val = ps_grid%zval
 
    ps%nlcc = ps_grid%core_corrections
 
    ps%h(0:ps%lmax, 1, 1) = ps_grid%dkbcos(1:ps%lmax+1)
 
    ! Increasing radius a little, just in case.
    ! I have hard-coded a larger increase of the cutoff for the filtering.
    ps%rc_max = maxval(ps_grid%kb_radius(1:ps%lmax+1)) * 1.5_real64
 
    ! now we fit the splines
    call get_splines(ps_grid%g)
 
    ! Passes from Rydbergs to Hartrees.
    ps%h(0:ps%lmax,:,:)    = ps%h(0:ps%lmax,:,:)    / m_two
 
    pop_sub(ps_grid_load)
 
  contains
 
    subroutine get_splines(g)
      type(logrid_t), intent(in) :: g
 
      real(real64), allocatable :: hato(:), dens(:)
      integer :: is, l, ir, nrc, ip
 
      push_sub(ps_grid_load.get_splines)
 
      safe_allocate(hato(1:g%nrval))
      safe_allocate(dens(1:g%nrval))
 
      ! the wavefunctions
      do is = 1, ps%ispin
 
        dens = m_zero
 
        do l = 1, ps_grid%no_l_channels
          hato(2:) = ps_grid%rphi(2:, l, 1+is)/g%rofi(2:)
          hato(1)  = first_point_extrapolate(g%rofi, hato)
 
          if(ps%conf%occ(l, is) > m_epsilon) then
            do ip = 1, g%nrval
              dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(m_four*m_pi)
            end do
          end if
 
          call spline_fit(g%nrval, g%rofi, hato, ps%ur(l, is), ps%projectors_sphere_threshold)
          call spline_fit(g%nrval, g%r2ofi, hato, ps%ur_sq(l, is), ps%projectors_sphere_threshold)
 
        end do
        call spline_fit(g%nrval, g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
      end do
 
 
      ! the Kleinman-Bylander projectors
      do l = 1, ps%lmax+1
        nrc = logrid_index(g, ps_grid%kb_radius(l)) + 1
        hato(1:nrc)         = ps_grid%KB(1:nrc, l)
        hato(nrc+1:g%nrval) = m_zero
 
        call spline_fit(g%nrval, g%rofi, hato, ps%kb(l-1, 1), ps%projectors_sphere_threshold)
      end do
 
      ! Now the part corresponding to the local pseudopotential
      ! where the asymptotic part is subtracted
      hato(:) = ps_grid%vlocal(:)/m_two
      call spline_fit(g%nrval, g%rofi, hato, ps%vl, ps%projectors_sphere_threshold)
 
      if (ps_grid%core_corrections) then
        ! find cutoff radius
        hato(2:) = ps_grid%chcore(2:)/(m_four*m_pi*g%r2ofi(2:))
 
        do ir = g%nrval-1, 2, -1
          if (hato(ir) > eps) then
            nrc = ir + 1
            exit
          end if
        end do
 
        hato(nrc:g%nrval) = m_zero
        hato(1) = first_point_extrapolate(g%rofi, hato)
 
        call spline_fit(g%nrval, g%rofi, hato, ps%core, ps%projectors_sphere_threshold)
      end if
 
      safe_deallocate_a(hato)
      safe_deallocate_a(dens)
 
      pop_sub(ps_grid_load.get_splines)
    end subroutine get_splines
  end subroutine ps_grid_load
 
  ! ---------------------------------------------------------
  subroutine ps_xml_load(ps, ps_xml, namespace)
    type(ps_t),        intent(inout) :: ps
    type(ps_xml_t),    intent(in)    :: ps_xml
    type(namespace_t), intent(in)    :: namespace
 
    integer :: ll, ip, is, ic, jc, ir, nrc, ii
    real(real64) :: rr, kbcos, kbnorm, dnrm, avgv, volume_element
    real(real64), allocatable :: vlocal(:), kbprojector(:), wavefunction(:), nlcc_density(:), dens(:,:)
    integer, allocatable :: cmap(:, :)
    real(real64), allocatable :: matrix(:, :), eigenvalues(:)
    logical :: density_is_known
    integer, allocatable :: order(:)
    real(real64), allocatable :: occ_tmp(:,:)
    real(real64), parameter :: tol_diagonal=1.0e-10_real64 ! tolerance for taking d_ij as a diagonal matrix
 
    push_sub(ps_xml_load)
 
    ! Vanderbilt-Kleinman-Bylander, i.e. more than one projector per l. Hamann uses 2 for ONCVPSP.
    if (pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
      if (ps%file_format == pseudo_format_psp8) then
        ps%relativistic_treatment = proj_j_average
        message(1) = "SOC from PSP8 is not currently supported."
        message(2) = "Only scalar relativistic effects will be considered."
        call messages_warning(2, namespace=namespace)
      else
        ps%relativistic_treatment = proj_j_dependent
      end if
    end if
 
    ps%nlcc = ps_xml%nlcc
 
    ps%z_val = ps_xml%valence_charge
 
    ps%has_density = ps_xml%has_density
 
    ! We start with the local potential
 
    ! Outside of the grid, we use the analytical long-range tail
    ! TODO: check usefulness: ps%g%nrval is initialized as ps_xml%grid_size
    ! This makes no sense, ps_xml%grid can only go up to ps_xml%grid_size
    ! We should just fit ps_xml%potential directly
    ! Also, if ip > ps_xml%grid_size, most likely ps_xml%grid is not defined
    safe_allocate(vlocal(1:ps%g%nrval))
    do ip = 1, ps_xml%grid_size
      rr = ps_xml%grid(ip)
      if (ip <= ps_xml%grid_size) then
        vlocal(ip) = ps_xml%potential(ip, ps%llocal)
      else
        vlocal(ip) = -ps_xml%valence_charge/rr
      end if
    end do
 
    ! We then fit this by a spline
    call spline_fit(ps%g%nrval, ps%g%rofi, vlocal, ps%vl, ps%projectors_sphere_threshold)
 
    safe_deallocate_a(vlocal)
 
    safe_allocate(kbprojector(1:ps%g%nrval))
    safe_allocate(wavefunction(1:ps%g%nrval))
 
    kbprojector = m_zero
    wavefunction = m_zero
 
    density_is_known = .false.
 
    ! We then proceed with the nonlocal projectors and the orbitals
    if (ps_xml%kleinman_bylander) then
 
      safe_allocate(cmap(0:ps_xml%lmax, 1:ps_xml%nchannels))
 
      ! We order the different projectors and create some mappings
      ! the order of the channels is determined by spin orbit and the j value
      do ll = 0, ps_xml%lmax
        do ic = 1, ps_xml%nchannels
          cmap(ll, ic) = ic
 
          if (ll == 0) cycle
          if (ll == ps_xml%llocal) cycle
          if (ps%relativistic_treatment /= proj_j_dependent) cycle
          if (ic > pseudo_nprojectors_per_l(ps_xml%pseudo, ll)) cycle ! This occurs for O and F for pseudodojo FR PBE for instance
 
          ! This is Octopus convention:
          ! We treat l+1/2 first and l-1/2 after, so we order the projectors this way
 
          assert(mod(ps_xml%nchannels, 2)==0)
          if (pseudo_projector_2j(ps_xml%pseudo, ll, ic) == 2*ll - 1) then
            cmap(ll, ic) = int((ic-1)/2)*2 + 2
          else
            assert(pseudo_projector_2j(ps_xml%pseudo, ll, ic) == 2*ll + 1)
            cmap(ll, ic) = int((ic-1)/2)*2 + 1
          end if
 
        end do
 
        ! check that all numbers are present for each l
        assert(sum(cmap(ll, 1:ps_xml%nchannels)) == (ps_xml%nchannels + 1)*ps_xml%nchannels/2)
      end do
 
      assert(all(cmap >= 0 .and. cmap <= ps_xml%nchannels))
 
 
      ! Weight of the projectors
      ! This is a diagonal matrix after the following treatment
      ps%h = m_zero
 
      ! We now take the matrix dij and we diagonalize it
      ! The nonlocal projectors are changed accordingly, and then fitted by splines
      if (pseudo_nprojectors(ps_xml%pseudo) > 0) then
        safe_allocate(matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels))
        safe_allocate(eigenvalues(1:ps_xml%nchannels))
 
        do ll = 0, ps_xml%lmax
 
          if (is_diagonal(ps_xml%nchannels, ps_xml%dij(ll, :, :), tol_diagonal) .or. &
            pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
            matrix = m_zero
            do ic = 1, ps_xml%nchannels
              eigenvalues(ic) = ps_xml%dij(ll, ic, ic)
              matrix(ic, ic) = m_one
            end do
          else
            ! diagonalize the coefficient matrix
            matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels) = ps_xml%dij(ll, 1:ps_xml%nchannels, 1:ps_xml%nchannels)
            call lalg_eigensolve(ps_xml%nchannels, matrix, eigenvalues)
          end if
 
          do ic = 1, ps_xml%nchannels
            kbprojector = m_zero
            do jc = 1, ps_xml%nchannels
              call lalg_axpy(ps_xml%grid_size, matrix(jc, ic), ps_xml%projector(:, ll, jc), kbprojector)
            end do
 
            call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, cmap(ll, ic)), ps%projectors_sphere_threshold)
 
            ps%h(ll, cmap(ll, ic), cmap(ll, ic)) = eigenvalues(ic)
 
          end do
        end do
 
        safe_deallocate_a(matrix)
        safe_deallocate_a(eigenvalues)
      end if
 
      ps%conf%p = ps_xml%nwavefunctions
 
      ! If we do not have a density but we have wavefunctions, we compute the
      ! pseudo-atomic density, from the pseudo-atomic wavefunctions
      ! In case we are doing a spin-polarized calculation, it is better to use the
      ! wavefunctions with spin-dependent occupations, instead the spin unresolved density
      ! given in the pseudopotential
      if ((.not. ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0) then
        safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
        dens = m_zero
      end if
 
      do ii = 1, ps_xml%nwavefunctions
 
        ps%conf%n(ii) = ps_xml%wf_n(ii)
        ps%conf%l(ii) = ps_xml%wf_l(ii)
 
        if (ps%ispin == 2) then
          ps%conf%occ(ii, 1) = min(ps_xml%wf_occ(ii), m_two*ps_xml%wf_l(ii) + m_one)
          ps%conf%occ(ii, 2) = ps_xml%wf_occ(ii) - ps%conf%occ(ii, 1)
        else
          ps%conf%occ(ii, 1) = ps_xml%wf_occ(ii)
        end if
 
        ps%conf%j(ii) = m_zero
        if (pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
          ps%conf%j(ii) = m_half*pseudo_wavefunction_2j(ps_xml%pseudo, ii)
        end if
 
        do ip = 1, ps%g%nrval
          if (ip <= ps_xml%grid_size) then
            wavefunction(ip) = ps_xml%wavefunction(ip, ii)
          else
            wavefunction(ip) = m_zero
          end if
        end do
 
        do is = 1, ps%ispin
          call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ii, is), ps%projectors_sphere_threshold)
          call spline_fit(ps%g%nrval, ps%g%r2ofi, wavefunction, ps%ur_sq(ii, is), ps%projectors_sphere_threshold)
        end do
 
        if (.not. ps_has_density(ps) .or. ps%ispin == 2) then
          do is = 1, ps%ispin
            do ip = 1, ps_xml%grid_size
              dens(ip, is) = dens(ip, is) + ps%conf%occ(ii, is)*wavefunction(ip)**2/(m_four*m_pi)
            end do
          end do
        end if
 
      end do
 
      !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
      if ((.not. ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0) then
        do is = 1, ps%ispin
          call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
        end do
        safe_deallocate_a(dens)
        density_is_known = .true.
        ps%has_density = .true.
      end if
 
      safe_deallocate_a(cmap)
 
    else !Not ps_xml%kleinman_bylander
 
      !Get the occupations from the valence charge of the atom
      ps%conf%occ = m_zero
      ps%conf%symbol = ps%label(1:2)
      call ps_guess_atomic_occupations(namespace, ps%z, ps%z_val, ps%ispin, ps%conf)
 
      ! In order to work in the following, we need to sort the occupations by angular momentum
      safe_allocate(order(1:ps_xml%lmax+1))
      safe_allocate(occ_tmp(1:ps_xml%lmax+1, 1:2))
      occ_tmp = ps%conf%occ
      call sort(ps%conf%l(1:ps_xml%lmax+1), order)
      do ll = 0, ps_xml%lmax
        ps%conf%occ(ll+1, 1:2) = occ_tmp(order(ll+1), 1:2)
      end do
      safe_deallocate_a(order)
      safe_deallocate_a(occ_tmp)
 
      !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
      if (abs(sum(ps%conf%occ) - ps%z_val ) < m_epsilon) then
        safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
        dens = m_zero
      end if
 
      do ll = 0, ps_xml%lmax
        ! we need to build the KB projectors
        ! the procedure was copied from ps_in_grid.F90 (r12967)
        dnrm = m_zero
        avgv = m_zero
        do ip = 1, ps_xml%grid_size
          rr = ps_xml%grid(ip)
          volume_element = rr**2*ps_xml%weights(ip)
          kbprojector(ip) = (ps_xml%potential(ip, ll) - ps_xml%potential(ip, ps%llocal))*ps_xml%wavefunction(ip, ll)
          dnrm = dnrm + kbprojector(ip)**2*volume_element
          avgv = avgv + kbprojector(ip)*ps_xml%wavefunction(ip, ll)*volume_element
        end do
 
        kbcos = dnrm/(safe_tol(avgv,1.0e-20_real64))
        kbnorm = m_one/(safe_tol(sqrt(dnrm),1.0e-20_real64))
 
        if (ll /= ps%llocal) then
          ps%h(ll, 1, 1) = kbcos
          kbprojector = kbprojector*kbnorm
        else
          ps%h(ll, 1, 1) = m_zero
        end if
 
        call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, 1), ps%projectors_sphere_threshold)
 
        ! wavefunctions, for the moment we pad them with zero
        call lalg_copy(ps_xml%grid_size, ps_xml%wavefunction(:, ll), wavefunction)
        wavefunction(ps_xml%grid_size+1:ps%g%nrval) = m_zero
 
        do is = 1, ps%ispin
          call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ll + 1, is), ps%projectors_sphere_threshold)
          call spline_fit(ps%g%nrval, ps%g%r2ofi, wavefunction, ps%ur_sq(ll + 1, is), ps%projectors_sphere_threshold)
        end do
 
        !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
        if (abs(sum(ps%conf%occ) - ps%z_val) < m_epsilon) then
          do is = 1, ps%ispin
            do ip = 1, ps_xml%grid_size
              dens(ip, is) = dens(ip, is) + ps%conf%occ(ll+1, is)*wavefunction(ip)**2/(m_four*m_pi)
            end do
          end do
        end if
      end do
 
      !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
      if (abs(sum(ps%conf%occ) - ps%z_val) < m_epsilon) then
        do is = 1, ps%ispin
          call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
        end do
        safe_deallocate_a(dens)
        density_is_known = .true.
        ps%has_density = .true.
      end if
    end if ! ps_xml%kleinman_bylander
 
 
    if (ps_has_density(ps) .and. .not. density_is_known) then
 
      safe_allocate(dens(1:ps%g%nrval, 1))
 
      dens(1:ps_xml%grid_size, 1) = ps_xml%density(1:ps_xml%grid_size)/ps%ispin
      dens(ps_xml%grid_size + 1:ps%g%nrval, 1) = m_zero
 
      do is = 1, ps%ispin
        call spline_fit(ps%g%nrval, ps%g%rofi, dens(:, 1), ps%density(is), ps%projectors_sphere_threshold)
      end do
 
      safe_deallocate_a(dens)
    end if
 
    ! Non-linear core-corrections
    ! We truncate the NLCC when below eps
    if (ps_xml%nlcc) then
 
      safe_allocate(nlcc_density(1:ps%g%nrval))
 
      call lalg_copy(ps_xml%grid_size, ps_xml%nlcc_density, nlcc_density)
 
      ! find cutoff radius
      do ir = ps_xml%grid_size - 1, 1, -1
        if (nlcc_density(ir) > eps) then
          nrc = ir + 1
          exit
        end if
      end do
 
      nlcc_density(nrc:ps%g%nrval) = m_zero
 
      call spline_fit(ps%g%nrval, ps%g%rofi, nlcc_density, ps%core, ps%projectors_sphere_threshold)
 
      safe_deallocate_a(nlcc_density)
    end if
 
    call ps_getradius(ps)
 
    !To be consistent with the other pseudopotentials, we are increasing the radius here
    ps%rc_max = ps%rc_max*1.05_real64
 
    safe_deallocate_a(kbprojector)
    safe_deallocate_a(wavefunction)
 
    pop_sub(ps_xml_load)
  end subroutine ps_xml_load
 
 
  ! ---------------------------------------------------------
  pure integer function ps_niwfs(ps)
    type(ps_t), intent(in) :: ps
 
    integer :: i, l
 
    ps_niwfs = 0
    do i = 1, ps%conf%p
      l = ps%conf%l(i)
      ps_niwfs = ps_niwfs + (2*l+1)
    end do
 
  end function ps_niwfs
 
  ! ---------------------------------------------------------
  pure integer function ps_bound_niwfs(ps)
    type(ps_t), intent(in) :: ps
 
    integer :: i, l
 
    ps_bound_niwfs = 0
    do i = 1, ps%conf%p
      l = ps%conf%l(i)
      if (any(.not. ps%bound(i,:))) cycle
      ps_bound_niwfs = ps_bound_niwfs + (2*l+1)
    end do
 
  end function ps_bound_niwfs
 
  !---------------------------------------
 
  pure logical function ps_has_density(ps) result(has_density)
    type(ps_t), intent(in) :: ps
 
    has_density = ps%has_density
 
  end function ps_has_density
 
  !---------------------------------------
 
  pure logical function ps_has_nlcc(ps) result(has_nlcc)
    type(ps_t), intent(in) :: ps
 
    has_nlcc = ps%nlcc
 
  end function ps_has_nlcc
 
  !---------------------------------------
  real(real64) function ps_density_volume(ps, namespace) result(volume)
    type(ps_t),        intent(in) :: ps
    type(namespace_t), intent(in) :: namespace
 
    integer :: ip, ispin
    real(real64) :: rr
    real(real64), allocatable ::vol(:)
    type(spline_t) :: volspl
 
    push_sub(ps_density_volume)
 
    if (.not. ps_has_density(ps)) then
      message(1) = "The pseudopotential does not contain an atomic density"
      call messages_fatal(1, namespace=namespace)
    end if
 
    safe_allocate(vol(1:ps%g%nrval))
 
    do ip = 1, ps%g%nrval
      rr = ps%g%rofi(ip)
      vol(ip) = m_zero
      do ispin = 1, ps%ispin
        vol(ip) = vol(ip) + spline_eval(ps%density(ispin), rr)*m_four*m_pi*rr**5
      end do
    end do
 
    call spline_init(volspl)
    call spline_fit(ps%g%nrval, ps%g%rofi, vol, volspl, ps%projectors_sphere_threshold)
    volume = spline_integral(volspl)
    call spline_end(volspl)
 
    safe_deallocate_a(vol)
 
    pop_sub(ps_density_volume)
  end function ps_density_volume
 
  subroutine ps_guess_atomic_occupations(namespace, zz, valcharge, ispin, conf)
    type(namespace_t), intent(in)     :: namespace
    real(real64),      intent(in)     :: zz
    real(real64),      intent(in)     :: valcharge
    integer,           intent(in)     :: ispin
    type(valconf_t),   intent(inout)  :: conf
 
    real(real64) :: val
 
    push_sub(ps_guess_atomic_occupations)
 
    val = valcharge
    conf%p = 0
    conf%l = 0
    conf%j = m_zero
    conf%n = 0
 
    write(message(1), '(a,a,a)') 'Debug: Guessing the atomic occupations for ', trim(conf%symbol), "."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    assert(valcharge <= zz)
 
    ! Here we populate the core states
    ! 1s state for all atoms after He
    if(int(zz) > 2 .and. val > zz - 2) then
      call fill_s_orbs(val, 2, 1)
    end if
    ! 2s state for all atoms after Be
    if(int(zz) > 4 .and. val > zz - 4) then
      call fill_s_orbs(val, 2, 2)
    end if
    ! 2p state for all atoms after Ne
    ! For pseudopotentials Al-Ar, we fill the 2s but not the 2p
    if(int(zz) > 18 .and. val > zz - 10) then
      call fill_p_orbs(val, 6, 2)
    end if
    ! 3s state for all atoms after Mg
    if(int(zz) > 12 .and. val > zz - 12) then
      call fill_s_orbs(val, 2, 3)
    end if
    ! 3p state for all atoms after Ar
    if(int(zz) > 18 .and. val > zz - 18) then
      call fill_p_orbs(val, 6, 3)
    end if
    ! 3d states for all atoms after Ni
    if(int(zz) > 28 .and. val > zz - 28) then
      call fill_d_orbs(val, 10, 3)
    end if
    ! 4s states for all atoms after Zn
    if(int(zz) > 30 .and. val > zz - 30) then
      call fill_s_orbs(val, 2, 4)
    end if
    ! 4p states for all atoms after Kr
    if(int(zz) > 36 .and. val > zz - 36) then
      call fill_p_orbs(val, 6, 4)
    end if
    ! 4d states for all atoms after Pd
    if(int(zz) > 46 .and. val > zz - 46) then
      call fill_d_orbs(val, 10, 4)
    end if
    ! 5s states for all atoms after Cd
    ! For Z=71 ot Z=80, the 4f is filled before the 5s/5p
    if((int(zz) > 48 .and. val > zz - 48) .or. &
      (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 62)) then
      call fill_s_orbs(val, 2, 5)
    end if
    ! 5p states for all atoms after Xe
    ! For Z=71 ot Z=80, the 4f is filled before the 5s/5p
    if((int(zz) > 54 .and. val > zz - 54) .or. &
      (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 68) ) then
      call fill_p_orbs(val, 6, 5)
    end if
    ! 4f states for all atoms after Yb
    ! Only after Z=80 for pseudopotentials
    if(int(zz) > 80 .and. val > zz - 68) then
      call fill_f_orbs(val, 14, 4)
    end if
 
 
    ! We leave here the valence states
    select case (int(zz))
    case (1)
      call fill_s_orbs(val, 1, 1) ! H 1s^1
    case (2)
      call fill_s_orbs(val, 2, 1) ! He 1s^2
    case (3)
      call fill_s_orbs(val, 1, 2) ! Li 1s^2 2s^1
    case (4)
      call fill_s_orbs(val, 2, 2) ! Be 1s^2 ; 2s^2
    case (5)
      call fill_p_orbs(val, 1, 2) ! B  1s^2 ; 2s^2 2p^1
    case (6)
      call fill_p_orbs(val, 2, 2) ! C       1s^2 ; 2s^2 2p^2
    case (7)
      call fill_p_orbs(val, 3, 2) ! N       1s^2 ; 2s^2 2
    case (8)
      call fill_p_orbs(val, 4, 2) ! O  1s^2 ; 2s^2 2p^4
    case (9)
      call fill_p_orbs(val, 5, 2) ! F       1s^2 ; 2s^2 2p^5
    case (10)
      call fill_p_orbs(val, 6, 2) ! Ne       1s^2 ; 2s^2 2p^6
    case (11)
      if(val > 6) call fill_p_orbs(val, 6, 2)
      call fill_s_orbs(val, 1, 3) ! Na      1s^2 2s^2 2p^6 ; 3s^1
    case (12)
      if(val > 6) call fill_p_orbs(val, 6, 2)
      call fill_s_orbs(val, 2, 3) ! Mg      1s^2 2s^2 2p^6 ; 3s^2
    case (13)
      if(val > 6) call fill_p_orbs(val, 6, 2)
      call fill_p_orbs(val, 1, 3) ! Al      1s^2 2s^2 2p^6 ; 3s^2 3p^1
    case (14)
      if(val > 6) call fill_p_orbs(val, 6, 2)
      call fill_p_orbs(val, 2, 3) ! Si      1s^2 2s^2 2p^6 ; 3s^2 3p^2
    case (15)
      if(val > 6) call fill_p_orbs(val, 6, 2)
      call fill_p_orbs(val, 3, 3) ! P       1s^2 2s^2 2p^6 ; 3s^2 3p^3
    case (16)
      if(val > 6) call fill_p_orbs(val, 6, 2)
      call fill_p_orbs(val, 4, 3) ! S       1s^2 2s^2 2p^6 ; 3s^2 3p^4
    case (17)
      if(val > 6) call fill_p_orbs(val, 6, 2)
      call fill_p_orbs(val, 5, 3) ! Cl      1s^2 2s^2 2p^6 ; 3s^2 3p^5
    case (18)
      if(val > 6) call fill_p_orbs(val, 6, 2)
      call fill_p_orbs(val, 6, 3) ! Ar      1s^2 2s^2 2p^6 ; 3s^2 3p^6
    case (19)
      call fill_s_orbs(val, 1, 4) ! K       1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^1
    case (20)
      call fill_s_orbs(val, 2, 4) ! Ca      1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2
    case (21)
      if (val > 1) call fill_s_orbs(val, 2, 4)
      call fill_d_orbs(val, 1, 3) ! Sc      1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^1
    case (22)
      if (val > 2) call fill_s_orbs(val, 2, 4)
      call fill_d_orbs(val, 2, 3) ! Ti      1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^2
    case (23)
      if (val > 3) call fill_s_orbs(val, 2, 4)
      call fill_d_orbs(val, 3, 3) ! V       1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^3
    case (24)
      if (val > 4) call fill_s_orbs(val, 2, 4)
      call fill_d_orbs(val, 4, 3) ! Cr      1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^4
    case (25)
      if (val > 5) call fill_s_orbs(val, 2, 4)
      call fill_d_orbs(val, 5, 3) ! Mn      1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^5
    case (26)
      if (val > 6) call fill_s_orbs(val, 2, 4)
      call fill_d_orbs(val, 6, 3) ! Fe      1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^6
    case (27)
      if (val > 7) call fill_s_orbs(val, 2, 4)
      call fill_d_orbs(val, 7, 3) ! Co      1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^7
    case (28)
      if (val > 8 ) call fill_s_orbs(val, 2, 4)
      call fill_d_orbs(val, 8, 3) ! Ni      1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^8
    case (29)
      call fill_s_orbs(val, 1, 4) ! Cu      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^1
    case (30)
      call fill_s_orbs(val, 2,  4) ! Zn      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2
    case (31)
      call fill_p_orbs(val, 1, 4) ! Ga      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^1
    case (32)
      call fill_p_orbs(val, 2, 4) ! Ge      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^2
    case (33)
      call fill_p_orbs(val, 3, 4) ! As      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^3
    case (34)
      call fill_p_orbs(val, 4, 4) ! Se      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^4
    case (35)
      call fill_p_orbs(val, 5, 4) ! Br      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^5
    case (36)
      call fill_p_orbs(val, 6, 4) ! Kr      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^6
    case (37)
      call fill_s_orbs(val, 1, 5) ! Rb      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 5s^1
    case (38)
      call fill_s_orbs(val, 2, 5) ! Sr      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 5s^2
    case (39)
      if (val > 2) call fill_d_orbs(val, 1, 4)
      call fill_s_orbs(val, 2, 5) ! Y       1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^1 5s^2
    case (40)
      if (val > 2) call fill_d_orbs(val, 2, 4)
      call fill_s_orbs(val, 2, 5) ! Zr      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^2 5s^2
    case (41)
      if (val > 1) call fill_d_orbs(val, 4, 4)
      call fill_s_orbs(val, 1, 5) ! Nb      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^4 5s^1
    case (42)
      if (val > 1) call fill_d_orbs(val, 5, 4)
      call fill_s_orbs(val, 1, 5) ! Mo      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^5 5s^1
    case (43)
      if (val > 2) call fill_d_orbs(val, 5, 4)
      call fill_s_orbs(val, 2, 5) ! Tc      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^5 5s^2
    case (44)
      if (val > 1) call fill_d_orbs(val, 7, 4)
      call fill_s_orbs(val, 1, 5) ! Ru      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^7 5s^1
    case (45)
      if (val > 1) call fill_d_orbs(val, 8, 4)
      call fill_s_orbs(val, 1, 5) ! Rh      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^8 5s^1
    case (46)
      call fill_d_orbs(val, 10, 4) ! Pd      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^10
    case (47)
      call fill_s_orbs(val, 1, 5) ! Ag      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^1
    case (48)
      call fill_s_orbs(val, 2, 5) ! Cd      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2
    case (49)
      call fill_p_orbs(val, 1, 5) ! In      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^1
    case (50)
      call fill_p_orbs(val, 2, 5) ! Sn      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^2
    case (51)
      call fill_p_orbs(val, 3, 5) ! Sb      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^3
    case (52)
      call fill_p_orbs(val, 4, 5) ! Te      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^4
    case (53)
      call fill_p_orbs(val, 5, 5) ! I       1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^5
    case (54)
      call fill_p_orbs(val, 6, 5) ! Xe      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^6
    case (55)
      call fill_s_orbs(val, 1, 6) ! Cs      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 6s^1
    case (56)
      call fill_s_orbs(val, 2, 6) ! Ba      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 6s^2
    case (57)
      if (val > 2) call fill_d_orbs(val, 1, 5)
      call fill_s_orbs(val, 2, 6) ! La      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 5d^1 6s^2
    case (58)
      if (val > 3) call fill_f_orbs(val, 1, 4)
      if (val > 2) call fill_d_orbs(val, 1, 5)
      call fill_s_orbs(val, 2, 6) ! Ce      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^1 5d^1 6s^2
    case (59)
      if (val > 2) call fill_f_orbs(val, 3, 4)
      call fill_s_orbs(val, 2, 6) ! Pr      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^3 6s^2
    case (60)
      if (val > 2) call fill_f_orbs(val, 4, 4)
      call fill_s_orbs(val, 2, 6) ! Nd      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^4 6s^2
    case (61)
      if (val > 2) call fill_f_orbs(val, 5, 4)
      call fill_s_orbs(val, 2, 6) ! Pm      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^5 6s^2
    case (62)
      if (val > 2) call fill_f_orbs(val, 6, 4)
      call fill_s_orbs(val, 2, 6) ! Sm      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^6 6s^2
    case (63)
      if (val > 2) call fill_f_orbs(val, 7, 4)
      call fill_s_orbs(val, 2, 6) ! Eu      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^7 6s^2
    case (64)
      if (val > 3) call fill_f_orbs(val, 7, 4)
      if (val > 2) call fill_d_orbs(val, 1, 5)
      call fill_s_orbs(val, 2, 6) ! Gd      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^7 5d^1 6s^2
    case (65)
      if (val > 2) call fill_f_orbs(val, 9, 4)
      call fill_s_orbs(val, 2, 6) ! Tb      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^9 6s^2
    case (66)
      if (val > 2) call fill_f_orbs(val, 10, 4)
      call fill_s_orbs(val, 2, 6) ! Dy      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^10 6s^2
    case (67)
      if (val > 2) call fill_f_orbs(val, 11, 4)
      call fill_s_orbs(val, 2, 6) ! Ho      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^11 6s^2
    case (68)
      if (val > 2) call fill_f_orbs(val, 12, 4)
      call fill_s_orbs(val, 2, 6) ! Er      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^12 6s^2
    case (69)
      if (val > 2) call fill_f_orbs(val, 13, 4)
      call fill_s_orbs(val, 2, 6) ! Tm      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^13 6s^2
    case (70)
      if (val > 2) call fill_f_orbs(val, 14, 4)
      call fill_s_orbs(val, 2, 6) ! Yb      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^14 6s^2
    case (71)
      if (val > 3) call fill_f_orbs(val, 14, 4)
      if (val > 2) call fill_d_orbs(val, 1, 5)
      call fill_s_orbs(val, 2, 6) ! Lu      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^1 6s^2
    case (72)
      if (val > 4) call fill_f_orbs(val, 14, 4)
      if (val > 2) call fill_d_orbs(val, 2, 5)
      call fill_s_orbs(val, 2, 6) ! Hf      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^2 6s^2
    case (73)
      if (val > 5) call fill_f_orbs(val, 14, 4)
      if (val > 2) call fill_d_orbs(val, 3, 5)
      call fill_s_orbs(val, 2, 6) ! Ta      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^3 6s^2
    case (74)
      if (val > 6) call fill_f_orbs(val, 14, 4)
      if (val > 2) call fill_d_orbs(val, 4, 5)
      call fill_s_orbs(val, 2, 6) ! W       1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^4 6s^2
    case (75)
      if (val > 7) call fill_f_orbs(val, 14, 4)
      if (val > 2) call fill_d_orbs(val, 5, 5)
      call fill_s_orbs(val, 2, 6) ! Re      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^5 6s^2
    case (76)
      if (val > 8) call fill_f_orbs(val, 14, 4)
      if (val > 2) call fill_d_orbs(val, 6, 5)
      call fill_s_orbs(val, 2, 6) ! Os      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^6 6s^2
    case (77)
      if (val > 9) call fill_f_orbs(val, 14, 4)
      if (val > 2) call fill_d_orbs(val, 7, 5)
      call fill_s_orbs(val, 2, 6) ! Ir      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^7 6s^2
    case (78)
      if (val > 10) call fill_f_orbs(val, 14, 4)
      if (val > 1) call fill_d_orbs(val, 9,   5)
      call fill_s_orbs(val, 1, 6) ! Pt      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^9 6s^1
    case (79)
      if (val > 21) call fill_f_orbs(val, 14, 4)
      if (val > 1) call fill_d_orbs(val, 10,  5)
      call fill_s_orbs(val, 1, 6) ! Au      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^1
    case (80)
      if (val > 12) call fill_f_orbs(val, 14, 4)
      if (val > 2) call fill_d_orbs(val, 10,  5)
      call fill_s_orbs(val, 2, 6) ! Hg      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2
    case (81)
      if (val > 10) call fill_d_orbs(val, 10, 5)
      if (val > 1) call fill_s_orbs(val, 2, 6)
      call fill_p_orbs(val, 1, 6) ! Tl      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^1
    case (82)
      if (val > 10) call fill_d_orbs(val, 10, 5)
      if (val > 2) call fill_s_orbs(val, 2, 6)
      call fill_p_orbs(val, 2, 6) ! Pb      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^2
    case (83)
      if (val > 10) call fill_d_orbs(val, 10, 5)
      if (val > 3) call fill_s_orbs(val, 2, 6)
      call fill_p_orbs(val, 3, 6) ! Bi      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^3
    case (84)
      if (val > 10) call fill_d_orbs(val, 10, 5)
      if (val > 4) call fill_s_orbs(val, 2, 6)
      call fill_p_orbs(val, 4, 6) ! Po      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^4
    case (85)
      if (val > 10) call fill_d_orbs(val, 10, 5)
      if (val > 5) call fill_s_orbs(val, 2, 6)
      call fill_p_orbs(val, 5, 6) ! At      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^5
    case (86)
      if (val > 10) call fill_d_orbs(val, 10, 5)
      if (val > 6) call fill_s_orbs(val, 2, 6)
      call fill_p_orbs(val, 6, 6) ! Rn      1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^6
    end select
 
    !If we attributed all the electrons, everything went fine
    if (val < m_epsilon) then
      !In case of spin-polarized calculations, we properly distribute the electrons
      if (ispin == 2) then
        call valconf_unpolarized_to_polarized(conf)
      end if
    else
      conf%occ = m_zero
      message(1) = "Error in attributing atomic occupations"
      call messages_warning(1, namespace=namespace)
    end if
 
    pop_sub(ps_guess_atomic_occupations)
 
  contains
    subroutine fill_s_orbs(val, max_occ, nn)
      real(real64), intent(inout) :: val
      integer, intent(in)  :: max_occ, nn
 
      conf%p = conf%p + 1
      conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
      val = val - conf%occ(conf%p,1)
      conf%l(conf%p) = 0
      conf%n(conf%p) = nn
    end subroutine fill_s_orbs
 
    subroutine fill_p_orbs(val, max_occ, nn)
      real(real64), intent(inout) :: val
      integer, intent(in)  :: max_occ, nn
 
      conf%p = conf%p + 1
      conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
      val = val - conf%occ(conf%p,1)
      conf%l(conf%p) = 1
      conf%n(conf%p) = nn
    end subroutine fill_p_orbs
 
    subroutine fill_d_orbs(val, max_occ, nn)
      real(real64), intent(inout) :: val
      integer, intent(in)  :: max_occ, nn
 
      conf%p = conf%p + 1
      conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
      val = val - conf%occ(conf%p,1)
      conf%l(conf%p) = 2
      conf%n(conf%p) = nn
    end subroutine fill_d_orbs
 
    subroutine fill_f_orbs(val, max_occ, nn)
      real(real64), intent(inout) :: val
      integer, intent(in)  :: max_occ, nn
 
      conf%p = conf%p + 1
      conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
      val = val - conf%occ(conf%p,1)
      conf%l(conf%p) = 3
      conf%n(conf%p) = nn
    end subroutine fill_f_orbs
 
  end subroutine ps_guess_atomic_occupations
 
#include "ps_pspio_inc.F90"
 
end module ps_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: