pseudopotential.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 pseudopotential_oct_m
  use debug_oct_m
  use element_oct_m
  use global_oct_m
  use io_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 pseudo_oct_m
  use pseudo_set_oct_m
  use species_oct_m
  use splines_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                        &
    pseudopotential_t,             &
    read_from_set,                 &
    read_from_default_file,        &
    get_set_directory,             &
    pseudopotential_real_nl_projector, &
    pseudopotential_nl_projector
 
  integer, public, parameter ::  &
    SPECIES_PSEUDO         = 7,             & !< pseudopotential
    species_pspio          = 110              
 
  type, extends(species_t) :: pseudopotential_t
    private
 
    integer                 :: type
    logical                 :: nlcc
 
    type(ps_t), allocatable, public :: ps
 
    integer :: user_lmax
    integer :: user_llocal
 
    integer, public :: pseudopotential_set_id
    logical, public :: pseudopotential_set_initialized
    type(pseudo_set_t), public :: pseudopotential_set
 
  contains
    procedure :: has_nlcc => pseudopotential_has_nlcc 
    procedure :: x_functional => pseudopotential_x_functional 
    procedure :: c_functional => pseudopotential_c_functional 
    procedure :: get_radius => pseudopotential_get_radius 
    procedure :: iwf_fix_qn => pseudopotential_iwf_fix_qn 
    procedure :: get_iwf_radius => pseudopotential_get_iwf_radius 
    procedure :: is_local => pseudopotential_is_local 
    procedure :: debug => pseudopotential_debug 
    procedure :: build => pseudopotential_build 
    procedure :: init_potential => pseudopotential_init_potential 
    procedure :: get_user_lmax => pseudopotential_get_user_lmax 
    procedure :: get_user_lloc => pseudopotential_get_user_lloc 
    procedure :: set_user_lmax => pseudopotential_set_user_lmax 
    procedure :: set_user_lloc => pseudopotential_set_user_lloc 
    procedure :: is_ps => pseudopotential_is_ps 
    procedure :: is_ps_with_nlcc => pseudopotential_is_ps_with_nlcc 
    procedure :: represents_real_atom => pseudopotential_represents_real_atom 
    final :: pseudopotential_finalize
  end type pseudopotential_t
 
  interface pseudopotential_t
    procedure pseudopotential_constructor
  end interface pseudopotential_t
 
contains
 
 
  ! ---------------------------------------------------------
  function pseudopotential_constructor(label, index) result(spec)
    class(pseudopotential_t),    pointer    :: spec
    character(len=*),         intent(in)    :: label
    integer,                  intent(in)    :: index
 
    push_sub(pseudopotential_constructor)
 
    safe_allocate(spec)
 
    call species_init(spec, label, index)
 
    spec%nlcc = .false.
 
    spec%user_lmax   = invalid_l
    spec%user_llocal = invalid_l
 
    spec%type = 0
 
    spec%pseudopotential_set_id = option__pseudopotentialset__none
    spec%pseudopotential_set_initialized = .false.
    call pseudo_set_nullify(spec%pseudopotential_set)
 
    pop_sub(pseudopotential_constructor)
  end function pseudopotential_constructor
 
 
  ! ---------------------------------------------------------
  subroutine pseudopotential_finalize(spec)
    type(pseudopotential_t), intent(inout) :: spec
 
    push_sub(pseudopotential_finalize)
 
    if (spec%pseudopotential_set_initialized) then
      call pseudo_set_end(spec%pseudopotential_set)
      spec%pseudopotential_set_initialized = .false.
    end if
 
    if (allocated(spec%ps)) call ps_end(spec%ps)
    safe_deallocate_a(spec%ps)
 
    call species_end(spec)
 
    pop_sub(pseudopotential_finalize)
  end subroutine pseudopotential_finalize
 
  ! ---------------------------------------------------------
  logical pure function pseudopotential_has_nlcc(spec)
    class(pseudopotential_t), intent(in) :: spec
    pseudopotential_has_nlcc = spec%nlcc
  end function pseudopotential_has_nlcc
 
 
  ! ---------------------------------------------------------
  integer pure function pseudopotential_x_functional(spec)
    class(pseudopotential_t), intent(in) :: spec
 
    pseudopotential_x_functional = spec%ps%exchange_functional
 
    ! if we do not know, try the pseudpotential set
    if (pseudopotential_x_functional == pseudo_exchange_unknown) then
      select case (spec%pseudopotential_set_id)
      case (                                     &
        option__pseudopotentialset__standard,   &
        option__pseudopotentialset__hgh_lda,    &
        option__pseudopotentialset__hgh_lda_sc, &
        option__pseudopotentialset__hscv_lda)
 
        pseudopotential_x_functional = option__xcfunctional__lda_x
 
      case (option__pseudopotentialset__hscv_pbe)
        pseudopotential_x_functional = option__xcfunctional__gga_x_pbe
      end select
    end if
 
  end function pseudopotential_x_functional
 
  ! ---------------------------------------------------------
 
  integer pure function pseudopotential_c_functional(spec)
    class(pseudopotential_t), intent(in) :: spec
 
    pseudopotential_c_functional = spec%ps%correlation_functional
 
    ! if we do not know, try the pseudpotential set
    if (pseudopotential_c_functional == pseudo_correlation_unknown) then
      select case (spec%pseudopotential_set_id)
      case (                                     &
        option__pseudopotentialset__standard,   &
        option__pseudopotentialset__hgh_lda,    &
        option__pseudopotentialset__hgh_lda_sc, &
        option__pseudopotentialset__hscv_lda)
 
        pseudopotential_c_functional = option__xcfunctional__lda_c_pz_mod/1000
 
      case (option__pseudopotentialset__hscv_pbe)
        pseudopotential_c_functional = option__xcfunctional__gga_c_pbe/1000
      end select
    end if
 
  end function pseudopotential_c_functional
 
  ! ---------------------------------------------------------
  real(real64) pure function pseudopotential_get_radius(spec) result(radius)
    class(pseudopotential_t),   intent(in) :: spec
 
    radius = spec%ps%rc_max
  end function pseudopotential_get_radius
 
  ! ---------------------------------------------------------
  integer pure function pseudopotential_get_user_lloc(spec) result(lloc)
    class(pseudopotential_t),   intent(in) :: spec
    lloc = spec%user_llocal
  end function pseudopotential_get_user_lloc
 
  ! ---------------------------------------------------------
  integer pure function pseudopotential_get_user_lmax(spec) result(lmax)
    class(pseudopotential_t),   intent(in) :: spec
    lmax = spec%user_lmax
  end function pseudopotential_get_user_lmax
 
  ! ---------------------------------------------------------
  pure subroutine pseudopotential_set_user_lmax(spec, ll)
    class(pseudopotential_t), intent(inout) :: spec
    integer,                  intent(in)    :: ll
    spec%user_lmax = ll
  end subroutine pseudopotential_set_user_lmax
 
  ! ---------------------------------------------------------
  pure subroutine pseudopotential_set_user_lloc(spec, ll)
    class(pseudopotential_t), intent(inout) :: spec
    integer,                  intent(in)    :: ll
    spec%user_llocal = ll
  end subroutine pseudopotential_set_user_lloc
 
  ! ---------------------------------------------------------
 
  character(len=MAX_PATH_LEN) function get_set_directory(set_id) result(filename)
    integer,         intent(in)   :: set_id
 
    push_sub(get_set_directory)
 
    select case (set_id)
    case (option__pseudopotentialset__standard)
      filename = trim(conf%share)//'/pseudopotentials/PSF'
    case (option__pseudopotentialset__sg15)
      filename = trim(conf%share)//'/pseudopotentials/quantum-simulation.org/sg15/'
    case (option__pseudopotentialset__hgh_lda)
      filename = trim(conf%share)//'/pseudopotentials/HGH/lda/'
    case (option__pseudopotentialset__hgh_lda_sc)
      filename = trim(conf%share)//'/pseudopotentials/HGH/lda_sc/'
    case (option__pseudopotentialset__hscv_lda)
      filename = trim(conf%share)//'/pseudopotentials/quantum-simulation.org/hscv/lda/'
    case (option__pseudopotentialset__hscv_pbe)
      filename = trim(conf%share)//'/pseudopotentials/quantum-simulation.org/hscv/pbe/'
    case (option__pseudopotentialset__pseudodojo_lda)
      filename = trim(conf%share)//'/pseudopotentials/pseudo-dojo.org/nc-sr-04_pw_standard/'
    case (option__pseudopotentialset__pseudodojo_pbe)
      filename = trim(conf%share)//'/pseudopotentials/pseudo-dojo.org/nc-sr-04_pbe_standard/'
    case (option__pseudopotentialset__pseudodojo_pbesol)
      filename = trim(conf%share)//'/pseudopotentials/pseudo-dojo.org/nc-sr-04_pbesol_standard/'
    case (option__pseudopotentialset__none)
      filename = ''
    case default
      assert(.false.)
    end select
 
    pop_sub(get_set_directory)
  end function get_set_directory
 
  ! ---------------------------------------------------------
  subroutine read_from_set(spec, set_id, set, read_data)
    type(pseudopotential_t), intent(inout) :: spec
    integer,                 intent(in)    :: set_id
    type(pseudo_set_t),      intent(in)    :: set
    integer,                 intent(out)   :: read_data
 
    type(element_t) :: el
 
    push_sub(read_from_set)
 
    spec%pseudopotential_set_id = set_id
    spec%pseudopotential_set = set
 
    call element_init(el, get_symbol(spec%get_label()))
 
    if (spec%pseudopotential_set_id /= option__pseudopotentialset__none .and. pseudo_set_has(spec%pseudopotential_set, el)) then
      call spec%set_filename(pseudo_set_file_path(spec%pseudopotential_set, el))
 
      ! these might have been set before
      if (spec%get_z() < 0) call spec%set_z(real(element_atomic_number(el), real64))
      if (spec%user_lmax == invalid_l) spec%user_lmax = pseudo_set_lmax(spec%pseudopotential_set, el)
      if (spec%user_llocal == invalid_l) spec%user_llocal = pseudo_set_llocal(spec%pseudopotential_set, el)
      if (spec%get_mass() < 0) call spec%set_mass(element_mass(el))
      if (spec%get_vdw_radius() < 0) call spec%set_vdw_radius(element_vdw_radius(el))
      read_data = 8
    else
      read_data = 0
    end if
 
    call element_end(el)
 
    pop_sub(read_from_set)
  end subroutine read_from_set
 
 
  ! ---------------------------------------------------------
  subroutine read_from_default_file(iunit, read_data, spec)
    integer,                  intent(in)    :: iunit
    integer,                  intent(inout) :: read_data
    class(pseudopotential_t), intent(inout) :: spec
 
    character(len=LABEL_LEN) :: label
    character(len=MAX_PATH_LEN) :: filename
    type(element_t) :: element
    integer :: lmax, llocal
    real(real64) :: zz
 
    push_sub(read_from_default_file)
 
    backspace(iunit)
 
    read(iunit,*) label, filename, zz, lmax, llocal
 
    call spec%set_filename(trim(conf%share)//'/pseudopotentials/'//trim(filename))
 
    assert(trim(label) == trim(spec%get_label()))
 
    read_data = 8
 
    ! get the mass, vdw radius and atomic number for this element
    call element_init(element, get_symbol(label))
 
    assert(element_valid(element))
 
    ! these might have been set before
    if (spec%get_z() < 0) call spec%set_z(zz)
    if (spec%get_z() < 0) call spec%set_z(real(element_atomic_number(element), real64))
    if (spec%user_lmax == invalid_l) spec%user_lmax = lmax
    if (spec%user_llocal == invalid_l) spec%user_llocal = llocal
    if (spec%get_mass() < 0) call spec%set_mass(element_mass(element))
    if (spec%get_vdw_radius() < 0) call spec%set_vdw_radius(element_vdw_radius(element))
 
    call element_end(element)
 
    pop_sub(read_from_default_file)
  end subroutine read_from_default_file
 
  ! ---------------------------------------------------------
  subroutine pseudopotential_real_nl_projector(spec, np, x, r, l, lm, i, uV)
    class(pseudopotential_t), intent(in)  :: spec
    integer,                  intent(in)  :: np
    real(real64),             intent(in)  :: x(:,:)
    real(real64),             intent(in)  :: r(:)
    integer,                  intent(in)  :: l, lm, i
    real(real64),             intent(out) :: uv(:)
 
    integer :: ip
    real(real64) :: ylm
 
    push_sub(pseudopotential_real_nl_projector)
 
    if (np > 0) then
      uv(1:np) = r(1:np)
      call spline_eval_vec(spec%ps%kb(l, i), np, uv)
 
      do ip = 1, np
        call ylmr_real(x(1:3, ip), l, lm, ylm)
        uv(ip) = uv(ip) * ylm
      end do
    end if
 
    pop_sub(pseudopotential_real_nl_projector)
  end subroutine pseudopotential_real_nl_projector
 
  ! ---------------------------------------------------------
  subroutine pseudopotential_nl_projector(spec, np, x, r, l, lm, i, uV)
    class(pseudopotential_t), intent(in)  :: spec
    integer,                  intent(in)  :: np
    real(real64),             intent(in)  :: x(:,:)
    real(real64),             intent(in)  :: r(:)
    integer,                  intent(in)  :: l, lm, i
    complex(real64),          intent(out) :: uv(:)
 
    integer :: ip
    complex(real64) :: ylm
 
    push_sub(pseudopotential_nl_projector)
 
    if (np > 0) then
      uv(1:np) = r(1:np)
      call spline_eval_vec(spec%ps%kb(l, i), np, uv)
 
      do ip = 1, np
        call ylmr_cmplx(x(1:3, ip), l, lm, ylm)
        uv(ip) = uv(ip) * ylm
      end do
    end if
 
    pop_sub(pseudopotential_nl_projector)
  end subroutine pseudopotential_nl_projector
 
  ! ---------------------------------------------------------
  subroutine pseudopotential_iwf_fix_qn(spec, namespace, nspin, dim)
    class(pseudopotential_t), intent(inout) :: spec
    type(namespace_t),        intent(in)    :: namespace
    integer,                  intent(in)    :: nspin
    integer,                  intent(in)    :: dim
 
    integer :: is, n, i, l, m
 
    push_sub(pseudopotential_iwf_fix_qn)
 
    do is = 1, nspin
      n = 1
      do i = 1, spec%ps%conf%p
        if (n > spec%niwfs) exit
        l = spec%ps%conf%l(i)
 
        if (.not. spec%ps%bound(i,is)) cycle
 
        do m = -l, l
          spec%iwf_i(n, is) = i
          spec%iwf_n(n, is) = spec%ps%conf%n(i)
          spec%iwf_l(n, is) = l
          spec%iwf_m(n, is) = m
          spec%iwf_j(n) = spec%ps%conf%j(i)
          n = n + 1
        end do
 
      end do
    end do
 
    pop_sub(pseudopotential_iwf_fix_qn)
  end subroutine pseudopotential_iwf_fix_qn
 
  ! ---------------------------------------------------------
  real(real64) function pseudopotential_get_iwf_radius(spec, ii, is, threshold) result(radius)
    class(pseudopotential_t), intent(in) :: spec
    integer,                  intent(in) :: ii
    integer,                  intent(in) :: is
    real(real64), optional,          intent(in) :: threshold
 
    real(real64) threshold_
 
    push_sub(pseudopotential_get_iwf_radius)
 
    threshold_ = optional_default(threshold, spec%ps%projectors_sphere_threshold)
    assert(ii <= spec%ps%conf%p)
    radius = spline_x_threshold(spec%ps%ur(ii, is), threshold_)
 
    pop_sub(pseudopotential_get_iwf_radius)
  end function pseudopotential_get_iwf_radius
 
  ! ---------------------------------------------------------
  logical function pseudopotential_is_local(spec) result(is_local)
    class(pseudopotential_t), intent(in) :: spec
 
    push_sub(pseudopotential_is_local)
 
    is_local = .false.
    if (spec%ps%lmax == 0 .and. spec%ps%llocal == 0) is_local = .true.
 
    pop_sub(pseudopotential_is_local)
  end function pseudopotential_is_local
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine pseudopotential_init_potential(this, namespace, grid_cutoff, filter)
    class(pseudopotential_t),  intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
    real(real64),              intent(in)    :: grid_cutoff
    integer,                   intent(in)    :: filter
 
    character(len=256) :: dirname
    integer            :: iorb
    real(real64) :: local_radius, orbital_radius
 
    push_sub(pseudopotential_init_potential)
 
    call ps_separate(this%ps)
 
    call ps_getradius(this%ps)
 
    if (filter /= ps_filter_none .and. this%ps%projector_type /= proj_hgh) then
      call ps_filter(this%ps, filter, grid_cutoff)
      call ps_getradius(this%ps) ! radius may have changed
    end if
 
    call ps_derivatives(this%ps)
 
    local_radius = this%ps%vl%x_threshold
 
    orbital_radius = m_zero
    ! FIXME: should take max over spins too here.
    do iorb = 1, this%get_niwfs()
      orbital_radius = max(orbital_radius, this%get_iwf_radius(this%iwf_i(iorb, 1), is = 1))
    end do
 
    call messages_write('Info: Pseudopotential for '//trim(this%get_label()), new_line = .true.)
    call messages_write('  Radii for localized parts:', new_line = .true.)
    call messages_write('    local part     = ')
    call messages_write(local_radius, fmt = 'f5.1', units = units_out%length, new_line = .true.)
    call messages_write('    non-local part = ')
    call messages_write(this%ps%rc_max, fmt = 'f5.1', units = units_out%length, new_line = .true.)
    call messages_write('    orbitals       = ')
    call messages_write(orbital_radius, fmt = 'f5.1', units = units_out%length, new_line = .true.)
    call messages_info(namespace=namespace)
 
    if (max(local_radius, this%ps%rc_max) > 6.0_real64) then
      call messages_write("One of the radii of your pseudopotential's localized parts seems", new_line = .true.)
      call messages_write("unusually large; check that your pseudopotential is correct.")
      call messages_warning(namespace=namespace)
    end if
 
    if (orbital_radius > 20.0_real64) then
      call messages_write("The radius of the atomic orbitals given by your pseudopotential seems", new_line = .true.)
      call messages_write("unusually large; check that your pseudopotential is correct.")
      call messages_warning(namespace=namespace)
    end if
 
    if (debug%info) then
      write(dirname, '(a)') 'debug/geometry'
      call io_mkdir(dirname, namespace)
      call this%debug(trim(dirname), namespace, grid_cutoff)
    end if
 
    pop_sub(pseudopotential_init_potential)
  end subroutine pseudopotential_init_potential
 
  ! ---------------------------------------------------------
  subroutine pseudopotential_debug(spec, dir, namespace, gmax)
    class(pseudopotential_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(pseudopotential_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,i3)')    'Type   = ', spec%type
    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,l1)')    'nlcc   = ', spec%nlcc
    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()
 
    if (debug%info) call ps_debug(spec%ps, trim(dirname), namespace, gmax)
 
    call io_close(iunit)
    pop_sub(pseudopotential_debug)
  end subroutine pseudopotential_debug
 
  ! ---------------------------------------------------------
  subroutine pseudopotential_build(spec, namespace, ispin, dim, print_info)
    class(pseudopotential_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(pseudopotential_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.
 
    ! allocate structure
    safe_allocate(spec%ps)
    if (spec%type == species_pspio) then
      call ps_pspio_init(spec%ps, namespace, spec%get_label(), spec%get_z(), spec%user_lmax, &
        ispin, spec%get_filename())
    else
      call ps_init(spec%ps, namespace, spec%get_label(), spec%get_z(), spec%user_lmax, &
        spec%user_llocal, ispin, spec%get_filename())
    end if
    call spec%set_zval(spec%ps%z_val)
    spec%nlcc = spec%ps%nlcc
    spec%niwfs = ps_bound_niwfs(spec%ps)
 
    ! invalidate these variables as they should not be used after
    spec%user_lmax = invalid_l
    spec%user_llocal = invalid_l
 
    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)
 
    pop_sub(pseudopotential_build)
  end subroutine pseudopotential_build
 
  ! ---------------------------------------------------------
  logical pure function pseudopotential_is_ps(this)
    class(pseudopotential_t),  intent(in) :: this
 
    pseudopotential_is_ps = .true.
  end function pseudopotential_is_ps
 
  ! ---------------------------------------------------------
  logical pure function pseudopotential_is_ps_with_nlcc(this)
    class(pseudopotential_t),  intent(in) :: this
 
    pseudopotential_is_ps_with_nlcc = this%has_nlcc()
  end function pseudopotential_is_ps_with_nlcc
 
  ! ---------------------------------------------------------
  logical pure function pseudopotential_represents_real_atom(spec)
    class(pseudopotential_t), intent(in) :: spec
 
    pseudopotential_represents_real_atom = .true.
  end function pseudopotential_represents_real_atom
 
 
 
end module pseudopotential_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: