jellium.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2023-2024 N. Tancogne-Dejean
!!
!! 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 jellium_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use species_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                        &
    jellium_t,                     &
    jellium_sphere_t,              &
    jellium_slab_t,                &
    jellium_charge_t,              &
    species_from_file_t,           &
    species_user_defined_t,        &
    species_charge_density_t
 
  integer, public, parameter ::  &
    SPECIES_JELLIUM        = 3,             & !< jellium sphere.
    species_jellium_slab   = 4,             & 
    species_jellium_charge_density = 129,   & 
    species_usdef          = 123,           & 
    species_charge_density = 125,           & 
    species_from_file      = 126
 
  type, abstract, extends(species_t) :: jellium_t
    private
 
    real(real64) :: omega
 
  contains
    procedure :: iwf_fix_qn => jellium_iwf_fix_qn
    procedure :: get_iwf_radius => jellium_get_iwf_radius
    procedure :: is_local => jellium_is_local
    procedure :: init_potential => jellium_init_potential
    procedure :: debug => jellium_debug
    procedure :: build => jellium_build
    procedure :: get_omega => jellium_get_omega
    procedure :: is_user_defined => jellium_user_defined
  end type jellium_t
 
  type, extends(jellium_t) :: jellium_sphere_t
    private
 
    real(real64) :: jradius
 
  contains
    procedure :: radius => jellium_radius
    procedure :: set_radius => jellium_set_radius
    final :: jellium_sphere_finalize
  end type jellium_sphere_t
 
  type, extends(jellium_t) :: jellium_slab_t
    private
 
    real(real64) :: jthick
 
  contains
    procedure :: thickness => jellium_thick
    procedure :: set_thickness => jellium_set_thickness
    procedure :: get_density => jellium_slab_density
    final :: jellium_slab_finalize
  end type jellium_slab_t
 
  type, extends(jellium_t) :: jellium_charge_t
    private
 
    character(len=200), public :: density_formula
 
  contains
    procedure :: rho_string => jellium_rho_string
    final :: jellium_charge_finalize
  end type jellium_charge_t
 
  type, extends (jellium_t) :: species_from_file_t
    private
 
  contains
    final :: species_from_file_finalize
  end type species_from_file_t
 
  type, extends (jellium_t) :: species_user_defined_t
    private
 
    character(len=1024), public :: potential_formula
  contains
    procedure :: user_pot => jellium_userdef_pot
    final :: species_user_defined_finalize
  end type species_user_defined_t
 
  type, extends (jellium_t) :: species_charge_density_t
    private
 
    character(len=200), public :: density_formula
  contains
    procedure :: rho_string => species_rho_string
    final :: species_charge_density_finalize
  end type species_charge_density_t
 
  interface jellium_sphere_t
    procedure jellium_sphere_constructor
  end interface jellium_sphere_t
 
  interface jellium_slab_t
    procedure jellium_slab_constructor
  end interface jellium_slab_t
 
  interface jellium_charge_t
    procedure jellium_charge_constructor
  end interface jellium_charge_t
 
  interface species_from_file_t
    procedure species_from_file_constructor
  end interface species_from_file_t
 
  interface species_user_defined_t
    procedure species_user_defined_constructor
  end interface species_user_defined_t
 
  interface species_charge_density_t
    procedure species_charge_density_constructor
  end interface species_charge_density_t
 
 
contains
 
 
  ! ---------------------------------------------------------
  function jellium_slab_constructor(label, index) result(spec)
    class(jellium_slab_t),   pointer    :: spec
    character(len=*),        intent(in) :: label
    integer,              intent(in)    :: index
 
    push_sub(jellium_slab_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%omega = m_zero
    spec%jthick = -m_one
 
    pop_sub(jellium_slab_constructor)
  end function jellium_slab_constructor
 
  ! ---------------------------------------------------------
  subroutine jellium_slab_finalize(this)
    type(jellium_slab_t), intent(inout) :: this
 
    push_sub(jellium_slab_finalize)
 
    call species_end(this)
 
    pop_sub(jellium_slab_finalize)
  end subroutine jellium_slab_finalize
 
  ! ---------------------------------------------------------
  function jellium_sphere_constructor(label, index) result(spec)
    class(jellium_sphere_t), pointer    :: spec
    character(len=*),        intent(in) :: label
    integer,              intent(in)    :: index
 
    push_sub(jellium_sphere_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%omega = m_zero
    spec%jradius = m_half
 
    pop_sub(jellium_sphere_constructor)
  end function jellium_sphere_constructor
 
 
  ! ---------------------------------------------------------
  subroutine jellium_sphere_finalize(this)
    type(jellium_sphere_t), intent(inout) :: this
 
    push_sub(jellium_sphere_finalize)
 
    call species_end(this)
 
    pop_sub(jellium_sphere_finalize)
  end subroutine jellium_sphere_finalize
 
  ! ---------------------------------------------------------
  function jellium_charge_constructor(label, index) result(spec)
    class(jellium_charge_t), pointer    :: spec
    character(len=*),        intent(in) :: label
    integer,              intent(in)    :: index
 
    push_sub(jellium_charge_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%omega = m_zero
    spec%density_formula = ""
 
    pop_sub(jellium_charge_constructor)
  end function jellium_charge_constructor
 
 
  ! ---------------------------------------------------------
  subroutine jellium_charge_finalize(this)
    type(jellium_charge_t), intent(inout) :: this
 
    push_sub(jellium_charge_finalize)
 
    call species_end(this)
 
    pop_sub(jellium_charge_finalize)
  end subroutine jellium_charge_finalize
 
 
  ! ---------------------------------------------------------
  function species_from_file_constructor(label, index) result(spec)
    class(species_from_file_t), pointer :: spec
    character(len=*),        intent(in) :: label
    integer,              intent(in)    :: index
 
    push_sub(species_from_file_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%omega = m_zero
 
    pop_sub(species_from_file_constructor)
  end function species_from_file_constructor
 
  ! ---------------------------------------------------------
  subroutine species_from_file_finalize(this)
    type(species_from_file_t), intent(inout) :: this
 
    push_sub(species_from_file_finalize)
 
    call species_end(this)
 
    pop_sub(species_from_file_finalize)
  end subroutine species_from_file_finalize
 
 
  ! ---------------------------------------------------------
  function species_user_defined_constructor(label, index) result(spec)
    class(species_user_defined_t), pointer :: spec
    character(len=*),           intent(in) :: label
    integer,                 intent(in)    :: index
 
    push_sub(species_user_defined_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%potential_formula = ""
    spec%omega = m_zero
 
    pop_sub(species_user_defined_constructor)
  end function species_user_defined_constructor
 
 
  ! ---------------------------------------------------------
  subroutine species_user_defined_finalize(this)
    type(species_user_defined_t), intent(inout) :: this
 
    push_sub(species_user_defined_finalize)
 
    call species_end(this)
 
    pop_sub(species_user_defined_finalize)
  end subroutine species_user_defined_finalize
 
 
  ! ---------------------------------------------------------
  function species_charge_density_constructor(label, index) result(spec)
    class(species_charge_density_t), pointer :: spec
    character(len=*),             intent(in) :: label
    integer,                      intent(in) :: index
 
    push_sub(species_charge_density_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%omega = m_zero
    spec%density_formula = ""
 
    pop_sub(species_charge_density_constructor)
  end function species_charge_density_constructor
 
 
  ! ---------------------------------------------------------
  subroutine species_charge_density_finalize(this)
    type(species_charge_density_t), intent(inout) :: this
 
    push_sub(species_charge_density_finalize)
 
    call species_end(this)
 
    pop_sub(species_charge_density_finalize)
  end subroutine species_charge_density_finalize
 
  ! ---------------------------------------------------------
  real(real64) pure function jellium_get_omega(spec)
    class(jellium_t), intent(in) :: spec
    jellium_get_omega = spec%omega
  end function jellium_get_omega
 
 
  ! ---------------------------------------------------------
  real(real64) pure function jellium_radius(spec)
    class(jellium_sphere_t), intent(in) :: spec
    jellium_radius = spec%jradius
  end function jellium_radius
 
  ! ---------------------------------------------------------
  pure subroutine jellium_set_radius(spec, radius)
    class(jellium_sphere_t), intent(inout) :: spec
    real(real64),            intent(in)    :: radius
    spec%jradius = radius
  end subroutine jellium_set_radius
 
  ! ---------------------------------------------------------
  real(real64) pure function jellium_thick(spec)
    class(jellium_slab_t), intent(in) :: spec
    jellium_thick = spec%jthick
  end function jellium_thick
 
  ! ---------------------------------------------------------
  pure subroutine jellium_set_thickness(spec, thick)
    class(jellium_slab_t),   intent(inout) :: spec
    real(real64),            intent(in)    :: thick
    spec%jthick = thick
  end subroutine jellium_set_thickness
 
  ! ---------------------------------------------------------
  character(len=200) pure function jellium_rho_string(spec)
    class(jellium_charge_t), intent(in) :: spec
    jellium_rho_string = trim(spec%density_formula)
  end function jellium_rho_string
 
  ! ---------------------------------------------------------
  character(len=200) pure function species_rho_string(spec)
    class(species_charge_density_t), intent(in) :: spec
    species_rho_string = trim(spec%density_formula)
  end function species_rho_string
 
 
  ! ---------------------------------------------------------
  complex(real64) function jellium_userdef_pot(spec, dim, xx, r)
    class(species_user_defined_t),  intent(in) :: spec
    integer,                        intent(in) :: dim
    real(real64),                   intent(in) :: xx(:)
    real(real64),                   intent(in) :: r
 
    real(real64) :: pot_re, pot_im
 
    ! no push_sub, called too many times
 
    call parse_expression(pot_re, pot_im, dim, xx, r, m_zero, spec%potential_formula)
    jellium_userdef_pot = cmplx(pot_re, pot_im, real64)
 
  end function jellium_userdef_pot
 
  ! ---------------------------------------------------------
  real(real64) function jellium_userdef_omega(spec, dim) result(omega)
    class(species_user_defined_t), intent(in) :: spec
    integer,                       intent(in) :: dim
 
    real(real64), parameter :: delta = 0.01_real64        
    real(real64), parameter :: omega_default = 0.1_real64 
    real(real64) :: xx(dim), vv(3), curvature, dv2
    integer :: idir, ii, ndir
 
    push_sub(jellium_userdef_omega)
 
    curvature = m_zero
    ndir = 0
    do idir = 1, dim
      do ii = 1, 3
        xx = m_zero
        xx(idir) = ii * delta
        vv(ii) = real(jellium_userdef_pot(spec, dim, xx, norm2(xx)), real64)
      end do
      ! Second derivative of the potential at 2*delta along the direction idir.
      ! For a harmonic potential, V = omega^2 * x^2 / 2, this gives exactly omega^2.
      dv2 = (vv(3) - m_two*vv(2) + vv(1)) / delta**2
      if (dv2 > m_zero) then
        curvature = curvature + dv2
        ndir = ndir + 1
      end if
    end do
 
    if (ndir > 0) then
      omega = sqrt(curvature / ndir)
    else
      ! The potential is not confining close to the origin (e.g. a constant potential or a
      ! double well), so no length scale can be extracted from it and we use a default value.
      omega = omega_default
    end if
 
    pop_sub(jellium_userdef_omega)
  end function jellium_userdef_omega
 
  ! ---------------------------------------------------------
  subroutine jellium_iwf_fix_qn(spec, namespace, nspin, dim)
    class(jellium_t),  intent(inout) :: spec
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: nspin
    integer,           intent(in)    :: dim
 
    integer :: is, i, n1, n2, n3
 
    push_sub(jellium_iwf_fix_qn)
 
    select case (dim)
    case (1)
      do is = 1, nspin
        do i = 1, spec%niwfs
          spec%iwf_i(i, is) = i
          spec%iwf_n(i, is) = 0
          spec%iwf_l(i, is) = 0
          spec%iwf_m(i, is) = 0
          spec%iwf_j(i) = m_zero
        end do
      end do
 
    case (2)
      do is = 1, nspin
        i = 1
        n1 = 1
        n2 = 1
        do
          spec%iwf_i(i, is) = n1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2
          spec%iwf_m(i, is) = 0
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          spec%iwf_i(i, is) = n1+1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2
          spec%iwf_m(i, is) = 0
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          spec%iwf_i(i, is) = n1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2+1
          spec%iwf_m(i, is) = 0
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          n1 = n1 + 1; n2 = n2 + 1
        end do
      end do
 
    case (3)
      do is = 1, nspin
        i = 1
        n1 = 1
        n2 = 1
        n3 = 1
        do
          spec%iwf_i(i, is) = n1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2
          spec%iwf_m(i, is) = n3
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          spec%iwf_i(i, is) = n1+1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2
          spec%iwf_m(i, is) = n3
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          spec%iwf_i(i, is) = n1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2+1
          spec%iwf_m(i, is) = 0
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          spec%iwf_i(i, is) = n1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2
          spec%iwf_m(i, is) = n3+1
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          spec%iwf_i(i, is) = n1+1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2+1
          spec%iwf_m(i, is) = n3
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          spec%iwf_i(i, is) = n1+1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2
          spec%iwf_m(i, is) = n3+1
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          spec%iwf_i(i, is) = n1
          spec%iwf_n(i, is) = 1
          spec%iwf_l(i, is) = n2+1
          spec%iwf_m(i, is) = n3+1
          spec%iwf_j(i) = m_zero
          i = i + 1
          if (i>spec%niwfs) exit
 
          n1 = n1 + 1
          n2 = n2 + 1
          n3 = n3 + 1
        end do
      end do
    case default
      ! Not doing anything to allow for N-D simulations
    end select
 
    pop_sub(jellium_iwf_fix_qn)
  end subroutine jellium_iwf_fix_qn
 
  ! ---------------------------------------------------------
  real(real64) pure function jellium_get_iwf_radius(spec, ii, is, threshold) result(radius)
    class(jellium_t),  intent(in) :: spec
    integer,           intent(in) :: ii
    integer,           intent(in) :: is
    real(real64), optional,   intent(in) :: threshold
 
    integer :: nmax
 
    real(real64) threshold_
 
    threshold_ = optional_default(threshold, 0.001_real64)
 
    ! We scale the radius with the quantum number. Here we assume 1D, which is an upper bound for higher dimensions.
    ! The increase in radius is given by the Hermite polynomial
    ! Note that for spinless particles, this is not correct
    nmax = max(floor(spec%get_zval()/ m_two - m_one), 1)
    radius = sqrt(-m_two*log(threshold_)/spec%omega) * sqrt(nmax + m_half)
 
    ! The values for hydrogenic and harmonic-oscillator wavefunctions
    ! come from taking the exponential part (i.e. the one that controls
    ! the asymptotic behavior at large r), and setting it equal to
    ! the threshold.
  end function jellium_get_iwf_radius
 
  ! ---------------------------------------------------------
  logical pure function jellium_is_local(spec) result(is_local)
    class(jellium_t), intent(in) :: spec
 
    is_local = .true.
  end function jellium_is_local
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine jellium_init_potential(this, namespace, grid_cutoff, filter)
    class(jellium_t), intent(inout) :: this
    type(namespace_t),    intent(in)    :: namespace
    real(real64),         intent(in)    :: grid_cutoff
    integer,              intent(in)    :: filter
 
    real(real64) :: spacing
 
    push_sub(jellium_init_potential)
 
    ! The guess orbitals are a Gaussian of width sigma = 1/sqrt(omega) times Hermite
    ! polynomials. A width narrower than the spacing cannot be represented on the mesh
    spacing = m_pi / grid_cutoff
    this%omega = min(this%omega, m_one/spacing**2)
 
    pop_sub(jellium_init_potential)
  end subroutine jellium_init_potential
 
  ! ---------------------------------------------------------
  subroutine jellium_debug(spec, dir, namespace, gmax)
    class(jellium_t),  intent(in) :: spec
    character(len=*),  intent(in) :: dir
    type(namespace_t), intent(in) :: namespace
    real(real64),      intent(in) :: gmax
 
    character(len=256) :: dirname
    integer :: iunit
 
    if (.not. mpi_world%is_root()) then
      return
    end if
 
    push_sub(jellium_debug)
 
    dirname = trim(dir)//'/'//trim(spec%get_label())
 
    call io_mkdir(dirname, namespace)
 
    iunit = io_open(trim(dirname)//'/info', namespace, action='write')
 
    write(iunit, '(a,i3)')    'Index  = ', spec%get_index()
    write(iunit, '(2a)')      'Label  = ', trim(spec%get_label())
    write(iunit, '(a,f15.2)') 'z_val  = ', spec%get_zval()
    write(iunit, '(a,f15.2)') 'mass = ', spec%get_mass()
    write(iunit, '(a,f15.2)') 'vdw_radius = ', spec%get_vdw_radius()
    write(iunit, '(a,l1)')    'local  = ', spec%is_local()
 
    select type(spec)
    type is(species_from_file_t)
      write(iunit, '(a,f15.2)') 'z      = ', spec%get_z()
      write(iunit,'(a)')      'Species read from file "'//trim(spec%get_filename())//'".'
    type is(jellium_sphere_t)
      write(iunit, '(a,f15.2)') 'z      = ', spec%get_z()
      write(iunit, '(a,f15.2)') 'jradius= ', spec%radius()
    type is(jellium_slab_t)
      write(iunit, '(a,f15.2)') 'z      = ', spec%get_z()
      write(iunit, '(a,f15.2)') 'jthick= ', spec%thickness()
    type is(species_user_defined_t)
      write(iunit, '(2a)')      'usdef  = ', trim(spec%potential_formula)
    end select
 
    write(iunit, '(a,i3)')    'hubbard_l = ', spec%get_hubbard_l()
    write(iunit, '(a,f15.2)') 'hubbard_U = ', spec%get_hubbard_U()
    write(iunit, '(a,f15.2)') 'hubbard_j = ', spec%get_hubbard_j()
    write(iunit, '(a,f15.2)') 'hubbard_alpha = ', spec%get_hubbard_alpha()
 
    call io_close(iunit)
    pop_sub(jellium_debug)
  end subroutine jellium_debug
 
  ! ---------------------------------------------------------
  subroutine jellium_build(spec, namespace, ispin, dim, print_info)
    class(jellium_t),  intent(inout) :: spec
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: ispin
    integer,           intent(in)    :: dim
    logical, optional, intent(in)    :: print_info
 
    logical :: print_info_
    integer :: i
 
    push_sub(jellium_build)
 
    print_info_ = optional_default(print_info, .true.)
 
    ! masses are always in amu, so convert them to a.u.
    call spec%set_mass(units_to_atomic(unit_amu, spec%get_mass()))
 
    spec%has_density = .false.
 
    select type (spec)
    type is(species_user_defined_t)
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "',trim(spec%get_label()),'" is a user-defined potential.'
        i = min(237, len_trim(spec%potential_formula)-1) ! I subtract 1 to avoid the non-printable C "end-of-string" character.
        write(message(2),'(a,a)')      '   Potential = ', trim(spec%potential_formula(1:i))
        if (len(trim(spec%potential_formula)) > 237) then
          message(2) = trim(message(2))//'...'
        end if
        call messages_info(2, namespace=namespace)
      end if
      spec%niwfs = int(max(2*spec%get_zval(), m_one))
 
      ! spec%omega is used to build the guess orbitals, given by a Gaussian radial enveloppe
      ! of width \sigma = 1/sqrt(omega) times Hermite polynomials. We get it from the
      ! curvature of the potential, which sets the length scale of the problem.
      spec%omega = jellium_userdef_omega(spec, dim)
 
    type is(species_from_file_t)
      if (print_info_) then
        write(message(1),'(a)') 'Species read from file "'//trim(spec%get_filename())//'".'
        call messages_info(1, namespace=namespace)
      end if
      spec%niwfs = 2*nint(spec%get_zval()+m_half)
      spec%omega = 0.1_real64
 
    type is(jellium_sphere_t)
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "', trim(spec%get_label()), &
          '" is a jellium sphere / approximated point particle.'
        write(message(2),'(a,f11.6)')  '   Valence charge = ', spec%get_zval()
        write(message(3),'(a,f11.6)')  '   Radius [a.u]   = ', spec%jradius
        write(message(4),'(a,f11.6)')  '   Rs [a.u]       = ', spec%jradius * spec%get_zval() ** (-m_one/m_three)
        call messages_info(4, namespace=namespace)
      end if
      spec%niwfs = species_closed_shell_size(2*nint(spec%get_zval()+m_half))
      spec%omega = 0.1_real64
 
    type is(jellium_slab_t)
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "',trim(spec%get_label()),'" is a jellium slab.'
        write(message(2),'(a,f11.6)')  '   Valence charge  = ', spec%get_zval()
        write(message(3),'(a,f11.6)')  '   Thickness [a.u] = ', spec%jthick
        !write(message(4),'(a,f11.6)')  '   Rs [a.u]       = ', ( M_THREE /( M_FOUR *M_PI ) &
        !& *spec%get_zval() /( *sb%lsize(1) *sb%lsize(2) ) )**(1.0/3.0)
        call messages_info(3, namespace=namespace)
      end if
      spec%niwfs = 2*nint(spec%get_zval()+m_half)
      spec%omega = 0.1_real64
 
    type is(jellium_charge_t)
      spec%niwfs = int(max(2*spec%get_zval(), m_one))
      spec%omega = spec%get_zval()
      spec%has_density = .true.
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "', trim(spec%get_label()), '" is a distribution of charge:'
        write(message(2),'(a,a,a)') '   rho is enclosed in volume defined by the "', &
          trim(spec%density_formula), '" block'
        write(message(3),'(a,f11.6)')  '   Z = ', spec%get_zval()
        call messages_info(3, namespace=namespace)
      end if
 
    type is(species_charge_density_t)
      spec%niwfs = int(max(2*spec%get_zval(), m_one))
      spec%omega = spec%get_zval()
      spec%has_density = .true.
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "', trim(spec%get_label()), '" is a distribution of charge:'
        write(message(2),'(a,a)')   '   rho = ', trim(spec%density_formula)
        write(message(3),'(a,f11.6)')  '   Z = ', spec%get_zval()
        call messages_info(3, namespace=namespace)
      end if
    class default
      call messages_input_error(namespace, 'Species', 'Unknown species type')
    end select
 
    ! since there is no real cap, make sure there are at least a few available
    spec%niwfs = max(5, spec%niwfs)
 
    safe_allocate(spec%iwf_n(1:spec%niwfs, 1:ispin))
    safe_allocate(spec%iwf_l(1:spec%niwfs, 1:ispin))
    safe_allocate(spec%iwf_m(1:spec%niwfs, 1:ispin))
    safe_allocate(spec%iwf_i(1:spec%niwfs, 1:ispin))
    safe_allocate(spec%iwf_j(1:spec%niwfs))
 
    call spec%iwf_fix_qn(namespace, ispin, dim)
 
    write(message(1),'(a,i6,a,i6)') 'Number of orbitals: ', spec%niwfs
    if (print_info_) call messages_info(1, namespace=namespace)
 
    pop_sub(jellium_build)
  end subroutine jellium_build
 
  ! ---------------------------------------------------------
  logical pure function jellium_user_defined(spec)
    class(jellium_t), intent(in) :: spec
 
    select type(spec)
    class is(species_user_defined_t)
      jellium_user_defined = .true.
    class is(species_charge_density_t)
      jellium_user_defined = .true.
    class is(species_from_file_t)
      jellium_user_defined = .true.
    class default
      jellium_user_defined = .false.
    end select
 
  end function jellium_user_defined
 
  ! ---------------------------------------------------------
  real(real64) pure function jellium_slab_density(slab, box_dim) result(density)
    class(jellium_slab_t), intent(in) :: slab
    real(real64),           intent(in) :: box_dim(:)
 
    ! Note the factor of 4 as box_dim is half of the box size
    density = slab%get_zval() / (box_dim(1) * box_dim(2) * m_four * slab%jthick)
 
  end function jellium_slab_density
end module jellium_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: