allelectron.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 allelectron_oct_m
  use atomic_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use loct_math_oct_m
  use logrid_oct_m
  use math_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use ps_oct_m
  use root_solver_oct_m
  use species_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                        &
    allelectron_t,                 &
    soft_coulomb_t,                &
    full_delta_t,                  &
    full_gaussian_t,               &
    full_anc_t
 
  type, abstract, extends(species_t) :: allelectron_t
    private
 
    real(real64) :: omega
    real(real64) :: sigma = -m_one       
 
    type(valconf_t), public  :: conf
 
  contains
    procedure :: get_sigma => allelectron_sigma 
    procedure :: get_omega => allelectron_omega 
    procedure :: set_sigma => allelectron_set_sigma 
    procedure :: iwf_fix_qn => allelectron_iwf_fix_qn 
    procedure :: get_iwf_radius => allelectron_get_iwf_radius 
    procedure :: is_local => allelectron_is_local 
    procedure :: init_potential => allelectron_init_potential 
    procedure :: debug => allelectron_debug 
    procedure :: build => allelectron_build 
    procedure :: is_full => allelectron_is_full 
    procedure :: represents_real_atom => allelectron_represents_real_atom 
  end type allelectron_t
 
  type, extends(allelectron_t) :: soft_coulomb_t
    private
 
    real(real64) :: softening
 
  contains
    procedure :: get_softening => soft_coulomb_softening 
    procedure :: set_softening => soft_coulomb_set_softening 
  end type soft_coulomb_t
 
  type, extends(allelectron_t) :: full_anc_t
    private
 
    real(real64) :: aa = -m_one          
    real(real64) :: bb = m_zero          
 
  contains
    procedure :: a => full_anc_a 
    procedure :: b => full_anc_b 
    procedure :: set_a => full_anc_set_a 
  end type full_anc_t
 
  type, extends(allelectron_t) :: full_delta_t
    private
 
  contains
  end type full_delta_t
 
 
  type, extends(allelectron_t) :: full_gaussian_t
    private
 
  contains
  end type full_gaussian_t
 
 
  interface soft_coulomb_t
    procedure soft_coulomb_constructor
  end interface soft_coulomb_t
 
  interface full_delta_t
    procedure full_delta_constructor
  end interface full_delta_t
 
  interface full_gaussian_t
    procedure full_gaussian_constructor
  end interface full_gaussian_t
 
  interface full_anc_t
    procedure full_anc_constructor
  end interface full_anc_t
 
  real(real64) :: alpha_p
  type(logrid_t), pointer :: grid_p
 
contains
 
 
  ! ---------------------------------------------------------
  function soft_coulomb_constructor(label, index) result(spec)
    class(soft_coulomb_t), pointer  :: spec
    character(len=*), intent(in)    :: label
    integer,          intent(in)    :: index
 
    push_sub(soft_coulomb_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%sigma = -m_one
    spec%omega = m_zero
    spec%softening = -m_one
 
    pop_sub(soft_coulomb_constructor)
  end function soft_coulomb_constructor
 
  ! ---------------------------------------------------------
  function full_anc_constructor(label, index, a) result(spec)
    class(full_anc_t), pointer      :: spec
    character(len=*), intent(in)    :: label
    integer,          intent(in)    :: index
    real(real64),     intent(in)    :: a
 
    push_sub(full_anc_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%sigma = -m_one
    spec%omega = m_zero
    spec%aa = a
    spec%bb = m_zero ! determined later by the code
 
    pop_sub(full_anc_constructor)
  end function full_anc_constructor
 
  ! ---------------------------------------------------------
  function full_delta_constructor(label, index, sigma) result(spec)
    class(full_delta_t), pointer    :: spec
    character(len=*), intent(in)    :: label
    integer,          intent(in)    :: index
    real(real64),     intent(in)    :: sigma
 
    push_sub(full_delta_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
    spec%sigma = sigma
    spec%omega = m_zero
 
    pop_sub(full_delta_constructor)
  end function full_delta_constructor
 
 
  ! ---------------------------------------------------------
  function full_gaussian_constructor(label, index, sigma) result(spec)
    class(full_gaussian_t), pointer :: spec
    character(len=*), intent(in)    :: label
    integer,          intent(in)    :: index
    real(real64),     intent(in)    :: sigma
 
    push_sub(full_gaussian_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%sigma = sigma
    spec%omega = m_zero
 
    pop_sub(full_gaussian_constructor)
  end function full_gaussian_constructor
 
 
  ! ---------------------------------------------------------
  real(real64) pure function allelectron_sigma(spec)
    class(allelectron_t), intent(in) :: spec
    allelectron_sigma = spec%sigma
  end function allelectron_sigma
 
  ! ---------------------------------------------------------
  pure subroutine allelectron_set_sigma(spec, sigma)
    class(allelectron_t), intent(inout) :: spec
    real(real64),         intent(in)    :: sigma
    spec%sigma = sigma
  end subroutine allelectron_set_sigma
 
  ! ---------------------------------------------------------
  real(real64) pure function allelectron_omega(spec)
    class(allelectron_t), intent(in) :: spec
    allelectron_omega = spec%omega
  end function allelectron_omega
 
  ! ---------------------------------------------------------
  real(real64) pure function full_anc_a(spec)
    class(full_anc_t), intent(in) :: spec
    full_anc_a = spec%aa
  end function full_anc_a
 
  ! ---------------------------------------------------------
  real(real64) pure function full_anc_b(spec)
    class(full_anc_t), intent(in) :: spec
    full_anc_b = spec%bb
  end function full_anc_b
 
  ! ---------------------------------------------------------
  pure subroutine full_anc_set_a(spec, a)
    class(full_anc_t), intent(inout) :: spec
    real(real64),      intent(in)    :: a
    spec%aa = a
  end subroutine full_anc_set_a
 
  ! ---------------------------------------------------------
  real(real64) pure function soft_coulomb_softening(spec)
    class(soft_coulomb_t), intent(in) :: spec
    soft_coulomb_softening = spec%softening
  end function soft_coulomb_softening
 
  ! ---------------------------------------------------------
  pure subroutine soft_coulomb_set_softening(spec, soft)
    class(soft_coulomb_t), intent(inout) :: spec
    real(real64),          intent(in)    :: soft
    spec%softening = soft
  end subroutine soft_coulomb_set_softening
 
 
  ! ---------------------------------------------------------
  subroutine allelectron_iwf_fix_qn(spec, namespace, nspin, dim)
    class(allelectron_t),  intent(inout) :: spec
    type(namespace_t),     intent(in)    :: namespace
    integer,               intent(in)    :: nspin
    integer,               intent(in)    :: dim
 
    integer :: is, n, i, l, m, n1, n2, ii, nn
    push_sub(allelectron_iwf_fix_qn)
 
 
    select type(spec)
    type is (soft_coulomb_t)
      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
          n = 1
          nn = 1
          ii = 1
          ! just up to the highest principal quantum number, actually
          do i = 1, spec%niwfs
            if (n > spec%niwfs) exit
            do l = 0, i-1
              do m = -l, l
                spec%iwf_i(n, is) = ii
                spec%iwf_n(n, is) = nn
                spec%iwf_l(n, is) = l
                spec%iwf_m(n, is) = m
                spec%iwf_j(n) = m_zero
                n = n + 1
              end do
              ii = ii + 1
            end do
            nn = nn + 1
          end do
        end do
 
      end select
 
    class default
      ! We need to find the occupations of the atomic wavefunctions to be able to get the atomic density
      ! Note that here we store the configuration is spec%ps%conf, even if we do not have a pseudo
      ! Not so clean...
      spec%conf%symbol = species_label_short(spec)
      call ps_guess_atomic_occupations(namespace, spec%get_z(), spec%get_zval(), nspin, spec%conf)
 
      ! Now that we have the occupations and the correct quantum numbers, we can attribute them
      do is = 1, nspin
        n = 1
        do i = 1, spec%conf%p
          if (n > spec%niwfs) exit
          l = spec%conf%l(i)
 
          do m = -l, l
            spec%iwf_i(n, is) = i
            spec%iwf_n(n, is) = spec%conf%n(i)
            spec%iwf_l(n, is) = l
            spec%iwf_m(n, is) = m
            spec%iwf_j(n) = spec%conf%j(i)
            n = n + 1
          end do
        end do
      end do
 
      ! We now overwrite spec%niwfs because we know now the number of orbitals
      spec%niwfs = n-1
 
    end select
 
    pop_sub(allelectron_iwf_fix_qn)
  end subroutine allelectron_iwf_fix_qn
 
  ! ---------------------------------------------------------
  real(real64) function allelectron_get_iwf_radius(spec, ii, is, threshold) result(radius)
    class(allelectron_t), intent(in) :: spec
    integer,              intent(in) :: ii
    integer,              intent(in) :: is
    real(real64), optional,      intent(in) :: threshold
 
    real(real64) threshold_
 
    push_sub(allelectron_get_iwf_radius)
 
    threshold_ = optional_default(threshold, 0.001_real64)
 
    select type(spec)
    type is(soft_coulomb_t)
      radius = -ii*log(threshold_)/spec%get_zval()
    class default ! full species
      ! For n=1, we use the direct analytical expression
      radius = -log(threshold_*sqrt(m_pi/(spec%get_zval()**3)))/spec%get_zval()
 
      ! For higher orbitals, we could use the first moment of the hydrogenic wavefunction
      ! See Physics of atoms and molecules, Bransden and Joachain
      ! Eq. 3.68
      ! radius = (real(ii, real64) **2/spec%get_z())*(M_ONE + M_HALF*(M_ONE-real(ll*(ll+1), real64)/(spec%get_z()**2)))
      ! However, these values are too small and leads to unbound orbitals.
      ! Hence, we use the scaling of the first moment, that goes as n^2 to scale the n=1 radius
      radius = radius * ii**2
    end select
 
    ! 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.
 
    pop_sub(allelectron_get_iwf_radius)
  end function allelectron_get_iwf_radius
 
  ! ---------------------------------------------------------
  logical pure function allelectron_is_local(spec) result(is_local)
    class(allelectron_t), intent(in) :: spec
 
    is_local = .true.
  end function allelectron_is_local
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine allelectron_init_potential(this, namespace, grid_cutoff, filter)
    class(allelectron_t), intent(inout) :: this
    type(namespace_t),    intent(in)    :: namespace
    real(real64),         intent(in)    :: grid_cutoff
    integer,              intent(in)    :: filter
 
    real(real64) :: grid_aa, grid_bb, bb(1), startval(1)
    integer :: grid_np
    logical :: conv
    type(root_solver_t) :: rs
 
    push_sub(allelectron_init_potential)
 
    select type(this)
    type is(full_anc_t)
      ! We first construct a logarithmic grid
      ! Then we find the value of b on this grid
      ! Finally, we evaluate the scaled potential following Eq. 19
      call logrid_find_parameters(namespace, int(this%get_z()), grid_aa, grid_bb, grid_np)
      safe_allocate(grid_p)
      call logrid_init(grid_p, logrid_psf, grid_aa, grid_bb, grid_np)
 
      alpha_p = this%aa
      ! We start from b=-0.1, which is close to the solution for a=4
      startval(1) = -0.1_real64
 
      ! solve equation
      call root_solver_init(rs, namespace, 1, solver_type=root_newton, maxiter=500, abs_tolerance=1.0e-20_real64)
      call droot_solver_run(rs, func_anc, bb, conv, startval=startval)
 
      this%bb = bb(1)
 
      if (.not. conv) then
        write(message(1),'(a)') 'Root finding in species_pot did not converge/'
        call messages_fatal(1, namespace=namespace)
      end if
 
      if(debug%info) then
        write(message(1),'(a, f12.6)')  'Debug: Optimized value of b for the ANC potential = ', this%bb
        call messages_info(1, namespace=namespace)
      end if
      call logrid_end(grid_p)
      safe_deallocate_p(grid_p)
    end select
 
    pop_sub(allelectron_init_potential)
  contains
    ! ---------------------------------------------------------
    ! In order to use the root finder, we need an auxiliary function
    ! that returns the value of the function to minimize and its Jacobian matrix
    ! given a value of the parameters.
    ! This routine does it for the ANC potential, where xin corresponds to the value b
    subroutine func_anc(xin, ff, jacobian)
      real(real64), intent(in)  :: xin(:)
      real(real64), intent(out) :: ff(:), jacobian(:,:)
 
      real(real64), allocatable :: rho(:)
      real(real64) :: norm
      integer :: ip
 
      push_sub(func_anc)
 
      safe_allocate(rho(1:grid_p%nrval))
      norm = m_one/sqrt(m_pi)
      do ip = 1, grid_p%nrval
        rho(ip) = -grid_p%rofi(ip) * loct_erf(grid_p%rofi(ip)*alpha_p) + xin(1)*exp(-alpha_p**2*grid_p%r2ofi(ip))
        rho(ip) =  norm * exp(rho(ip))
      end do
 
      ! First, we calculate the function to be minimized
      ff(1) = sum(grid_p%drdi*rho**2*grid_p%r2ofi) - m_one/m_four/m_pi
 
      ! Now the jacobian.
      do ip = 1, grid_p%nrval
        rho(ip) = m_two*rho(ip)**2*exp(-alpha_p**2*grid_p%r2ofi(ip))
      end do
      jacobian(1,1) = sum(grid_p%drdi*rho*grid_p%r2ofi)
 
      safe_deallocate_a(rho)
      pop_sub(func_anc)
    end subroutine func_anc
 
  end subroutine allelectron_init_potential
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  subroutine allelectron_debug(spec, dir, namespace, gmax)
    class(allelectron_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_grp_is_root(mpi_world)) then
      return
    end if
 
    push_sub(allelectron_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      = ', spec%get_z()
    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()
    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(allelectron_debug)
  end subroutine allelectron_debug
 
  ! ---------------------------------------------------------
  subroutine allelectron_build(spec, namespace, ispin, dim, print_info)
    class(allelectron_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_
 
    push_sub(allelectron_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 (soft_coulomb_t)
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "',trim(spec%get_label()),'" is a soft-Coulomb potential.'
        call messages_info(1, namespace=namespace)
      end if
      spec%niwfs = species_closed_shell_size(2*nint(spec%get_zval()+m_half))
      spec%omega = 0.1_real64
 
      ! since there is no real cap, make sure there are at least a few available
      spec%niwfs = max(5, spec%niwfs)
 
    type is (full_delta_t)
      spec%has_density = .true.
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "',trim(spec%get_label()),'" is an all-electron atom.'
        write(message(2),'(a,f11.6)')  '   Z = ', spec%get_z()
        write(message(3),'(a)')  '   Potential will be calculated solving the Poisson equation'
        write(message(4),'(a)')  '   for a delta density distribution.'
        call messages_info(4, namespace=namespace)
      end if
      spec%niwfs = species_closed_shell_size(floor(m_half*spec%get_zval()+m_half))
      ! We add one more shell because heavier elements might have not complete shells
      ! An example is Ag with 5s1 electrons but no 4f electrons
      ! We correct this number later
      spec%niwfs = species_closed_shell_size(spec%niwfs+1)
      spec%omega = spec%get_z()
 
    type is (full_gaussian_t)
      spec%has_density = .true.
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "',trim(spec%get_label()),'" is an all-electron atom.'
        write(message(2),'(a,f11.6)')  '   Z = ', spec%get_z()
        write(message(3),'(a)')  '   Potential will be calculated solving the Poisson equation'
        write(message(4),'(a)')  '   for a Gaussian density distribution.'
        call messages_info(4, namespace=namespace)
      end if
      spec%niwfs = species_closed_shell_size(floor(m_half*spec%get_zval()+m_half))
      ! We add one more shell because heavier elements might have not complete shells
      ! An example is Ag with 5s1 electrons but no 4f electrons
      ! We correct this number later
      spec%niwfs = species_closed_shell_size(spec%niwfs+1)
      spec%omega = spec%get_z()
      assert(spec%sigma > 0)
 
    type is (full_anc_t)
      spec%has_density = .false.
      if (print_info_) then
        write(message(1),'(a,a,a)')    'Species "',trim(spec%get_label()),'" is an all-electron atom.'
        write(message(2),'(a,f11.6)')  '   Z = ', spec%get_z()
        write(message(3),'(a)')  '   Potential is an analytical norm-conserving regularized potential'
        call messages_info(3, namespace=namespace)
      end if
      spec%niwfs = species_closed_shell_size(floor(m_half*spec%get_zval()+m_half))
      ! We add one more shell because heavier elements might have not complete shells
      ! An example is Ag with 5s2 electrons but no 4f electrons
      ! We correct this number later
      spec%niwfs = species_closed_shell_size(spec%niwfs+1)
      spec%omega = spec%get_z()
      assert(spec%aa > 0)
 
    class default
      call messages_input_error(namespace, 'Species', 'Unknown species type')
    end select
 
    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(allelectron_build)
  end subroutine allelectron_build
 
  ! ---------------------------------------------------------
  logical pure function allelectron_is_full(this)
    class(allelectron_t),        intent(in) :: this
 
    select type (this)
    class is (soft_coulomb_t)
      allelectron_is_full = .false.
    class default
      allelectron_is_full = .true.
    end select
  end function allelectron_is_full
 
  ! ---------------------------------------------------------
  logical pure function allelectron_represents_real_atom(spec)
    class(allelectron_t), intent(in) :: spec
 
    allelectron_represents_real_atom = .true.
  end function allelectron_represents_real_atom
 
 
 
end module allelectron_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: