species.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 species_oct_m
  use debug_oct_m
  use global_oct_m
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use string_oct_m
 
  implicit none
 
  private
  public ::                        &
    species_t,                     &
    species_wrapper_t,             &
    species_init,                  &
    species_end,                   &
    species_closed_shell_size,     &
    species_label_short,           &
    operator(==),                  &
    get_symbol
 
  interface operator(==)
    module procedure species_is_same_species
  end interface operator(==)
 
 
  integer, public, parameter :: LABEL_LEN=15
 
  type, abstract :: species_t
    private
    integer :: index
 
    character(len=LABEL_LEN) :: label
    real(real64)   :: z
    real(real64)   :: z_val
    !                                 !< minus the core charge in the case of the pseudopotentials
    real(real64)   :: mass
    real(real64)   :: vdw_radius
 
    logical, public :: has_density
 
    character(len=MAX_PATH_LEN) :: filename
 
    ! Derived classes need to access these quantities
    integer, public :: niwfs
    integer, allocatable, public :: iwf_l(:, :), iwf_m(:, :), iwf_i(:, :), iwf_n(:, :)
    real(real64), allocatable, public :: iwf_j(:)
 
    integer :: hubbard_l
    real(real64)   :: hubbard_U
    real(real64)   :: hubbard_j
    real(real64)   :: hubbard_alpha
 
  contains
    procedure(species_is_local),       deferred :: is_local
    procedure(species_get_iwf_radius), deferred :: get_iwf_radius
    procedure(species_iwf_fix_qn),     deferred :: iwf_fix_qn
    procedure(species_init_potential), deferred :: init_potential
    procedure(species_build),          deferred :: build
!    procedure, deferred :: get_average_energy !< G=0 contribution to the energy
!    procedure, deferred :: get_average_stress !< G=0 contribution to the stress tensor
    procedure(species_debug),          deferred :: debug
    procedure :: get_label => species_label
    procedure :: get_index => species_index
    procedure :: get_zval => species_zval
    procedure :: get_z => species_z
    procedure :: get_vdw_radius => species_vdw_radius
    procedure :: get_mass => species_mass
    procedure :: get_hubbard_l => species_hubbard_l
    procedure :: get_hubbard_u => species_hubbard_u
    procedure :: get_hubbard_j => species_hubbard_j
    procedure :: get_hubbard_alpha => species_hubbard_alpha
    procedure :: set_hubbard_l => species_set_hubbard_l
    procedure :: set_hubbard_u => species_set_hubbard_u
    procedure :: set_hubbard_j => species_set_hubbard_j
    procedure :: set_hubbard_alpha => species_set_hubbard_alpha
    procedure :: get_niwfs => species_niwfs
    procedure :: get_iwf_ilm => species_iwf_ilm
    procedure :: get_iwf_n => species_iwf_n
    procedure :: get_iwf_j => species_iwf_j
    procedure :: set_mass => species_set_mass
    procedure :: set_vdw_radius => species_set_vdw_radius
    procedure :: set_z => species_set_z
    procedure :: set_zval => species_set_zval
    procedure :: get_filename => species_filename
    procedure :: set_filename => species_set_filename
    procedure :: is_full => species_is_full
    procedure :: is_ps => species_is_ps
    procedure :: is_ps_with_nlcc => species_is_ps_with_nlcc
    procedure :: represents_real_atom => species_represents_real_atom
    procedure :: is_user_defined => species_user_defined
  end type species_t
 
  type species_wrapper_t
    class(species_t), pointer :: s
  end type
 
  abstract interface
    ! ---------------------------------------------------------
    logical function species_is_local(spec) result(is_local)
      import species_t
      class(species_t), intent(in) :: spec
    end function species_is_local
 
    ! ---------------------------------------------------------
    real(real64) function species_get_iwf_radius(spec, ii, is, threshold) result(radius)
      import real64
      import species_t
      class(species_t),         intent(in) :: spec
      integer,                  intent(in) :: ii
      integer,                  intent(in) :: is
      real(real64), optional,   intent(in) :: threshold
    end function species_get_iwf_radius
 
    ! ---------------------------------------------------------
    subroutine species_iwf_fix_qn(spec, namespace, nspin, dim)
      import species_t
      import namespace_t
      class(species_t),  intent(inout) :: spec
      type(namespace_t), intent(in)    :: namespace
      integer,           intent(in)    :: nspin
      integer,           intent(in)    :: dim
    end subroutine species_iwf_fix_qn
 
    ! ---------------------------------------------------------
    subroutine species_init_potential(this, namespace, grid_cutoff, filter)
      import species_t
      import namespace_t
      import real64
      class(species_t),    intent(inout) :: this
      type(namespace_t),   intent(in)    :: namespace
      real(real64),        intent(in)    :: grid_cutoff
      integer,             intent(in)    :: filter
    end subroutine species_init_potential
 
    ! ---------------------------------------------------------
    subroutine species_build(spec, namespace, ispin, dim, print_info)
      import species_t
      import namespace_t
      class(species_t),          intent(inout) :: spec
      type(namespace_t),         intent(in)    :: namespace
      integer,                   intent(in)    :: ispin
      integer,                   intent(in)    :: dim
      logical, optional,         intent(in)    :: print_info
    end subroutine species_build
 
    ! ---------------------------------------------------------
    subroutine species_debug(spec, dir, namespace, gmax)
      import species_t
      import namespace_t
      import real64
      class(species_t),        intent(in) :: spec
      character(len=*),        intent(in) :: dir
      type(namespace_t),       intent(in) :: namespace
      real(real64),            intent(in) :: gmax
    end subroutine species_debug
 
  end interface
 
  integer, private, parameter :: libxc_c_index = 1000
 
contains
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine species_init(this, label, index)
    class(species_t), intent(inout) :: this
    character(len=*), intent(in)    :: label
    integer,          intent(in)    :: index
 
    push_sub(species_init)
 
    this%label = trim(label)
    this%index = index
 
    this%z = -m_one
    this%z_val = -m_one
    this%mass = -m_one
    this%vdw_radius = -m_one
    this%has_density = .false.
    this%niwfs = -1
    this%hubbard_l = -1
    this%hubbard_U = m_zero
    this%hubbard_j = m_zero
    this%hubbard_alpha = m_zero
    this%filename = ""
 
    pop_sub(species_init)
  end subroutine species_init
 
  ! ---------------------------------------------------------
  integer pure function species_closed_shell_size(min_niwfs) result(size)
    integer, intent(in) :: min_niwfs
 
    integer :: nn
 
    size = 0
    do nn = 1, min_niwfs
      if (size >= min_niwfs) exit
      size = size + nn**2
    end do
 
  end function species_closed_shell_size
 
  ! ---------------------------------------------------------
  character(len=LABEL_LEN) pure function species_label(species)
    class(species_t), intent(in) :: species
    species_label = trim(species%label)
  end function species_label
 
  ! ---------------------------------------------------------
  character(len=2) pure function species_label_short(species)
    class(species_t), intent(in) :: species
    species_label_short = species%label(1:2)
  end function species_label_short
 
  ! ---------------------------------------------------------
  integer pure function species_index(species)
    class(species_t), intent(in) :: species
    species_index = species%index
  end function species_index
 
  ! ---------------------------------------------------------
  real(real64) pure function species_zval(species)
    class(species_t), intent(in) :: species
    species_zval = species%z_val
  end function species_zval
 
  ! ---------------------------------------------------------
  pure subroutine species_set_zval(species, zval)
    class(species_t), intent(inout) :: species
    real(real64),     intent(in)    :: zval
    species%z_val = zval
  end subroutine species_set_zval
 
  ! ---------------------------------------------------------
  real(real64) pure function species_z(species)
    class(species_t), intent(in) :: species
    species_z = species%z
  end function species_z
 
  ! ---------------------------------------------------------
  pure subroutine species_set_z(species, z)
    class(species_t), intent(inout) :: species
    real(real64),     intent(in)    :: z
    species%z = z
  end subroutine species_set_z
 
  ! ---------------------------------------------------------
  real(real64) pure function species_mass(species)
    class(species_t), intent(in) :: species
    species_mass = species%mass
  end function species_mass
 
  ! ---------------------------------------------------------
  pure subroutine species_set_mass(species, mass)
    class(species_t), intent(inout) :: species
    real(real64),     intent(in)    :: mass
    species%mass = mass
  end subroutine species_set_mass
 
  ! ---------------------------------------------------------
  real(real64) pure function species_vdw_radius(species)
    class(species_t), intent(in) :: species
    species_vdw_radius = species%vdw_radius
  end function species_vdw_radius
 
  ! ---------------------------------------------------------
  pure subroutine species_set_vdw_radius(species, radius)
    class(species_t), intent(inout) :: species
    real(real64),     intent(in)    :: radius
    species%vdw_radius = radius
  end subroutine species_set_vdw_radius
 
  ! ---------------------------------------------------------
  integer pure function species_niwfs(species)
    class(species_t), intent(in) :: species
    species_niwfs = species%niwfs
  end function species_niwfs
 
  ! ---------------------------------------------------------
  integer pure function species_hubbard_l(species)
    class(species_t), intent(in) :: species
    species_hubbard_l = species%hubbard_l
  end function species_hubbard_l
 
  ! ---------------------------------------------------------
  real(real64) pure function species_hubbard_u(species)
    class(species_t), intent(in) :: species
    species_hubbard_u = species%hubbard_u
  end function species_hubbard_u
 
  ! ---------------------------------------------------------
  real(real64) pure function species_hubbard_j(species)
    class(species_t), intent(in) :: species
    species_hubbard_j = species%hubbard_j
  end function species_hubbard_j
 
  ! ---------------------------------------------------------
  real(real64) pure function species_hubbard_alpha(species)
    class(species_t), intent(in) :: species
    species_hubbard_alpha = species%hubbard_alpha
  end function species_hubbard_alpha
 
  ! ---------------------------------------------------------
  pure subroutine species_set_hubbard_l(species, hubbard_l)
    class(species_t), intent(inout) :: species
    integer,          intent(in)    :: hubbard_l
    species%hubbard_l = hubbard_l
  end subroutine species_set_hubbard_l
 
  ! ---------------------------------------------------------
  pure subroutine species_set_hubbard_u(species, hubbard_u)
    class(species_t), intent(inout) :: species
    real(real64),     intent(in)    :: hubbard_u
    species%hubbard_u = hubbard_u
  end subroutine species_set_hubbard_u
 
  ! ---------------------------------------------------------
  pure subroutine species_set_hubbard_j(species, hubbard_j)
    class(species_t), intent(inout) :: species
    real(real64),     intent(in)    :: hubbard_j
    species%hubbard_j = hubbard_j
  end subroutine species_set_hubbard_j
 
  ! ---------------------------------------------------------
  pure subroutine species_set_hubbard_alpha(species, hubbard_alpha)
    class(species_t), intent(inout) :: species
    real(real64),     intent(in)    :: hubbard_alpha
    species%hubbard_alpha = hubbard_alpha
  end subroutine species_set_hubbard_alpha
 
  ! ---------------------------------------------------------
  character(len=200) pure function species_filename(species)
    class(species_t), intent(in) :: species
    species_filename = trim(species%filename)
  end function species_filename
 
  ! ---------------------------------------------------------
  pure subroutine species_set_filename(species, filename)
    class(species_t), intent(inout) :: species
    character(len=*), intent(in)    :: filename
    species%filename = trim(filename)
  end subroutine species_set_filename
 
  ! ---------------------------------------------------------
  pure subroutine species_iwf_ilm(species, j, is, i, l, m)
    class(species_t), intent(in) :: species
    integer, intent(in)          :: j, is
    integer, intent(out)         :: i, l, m
 
    i = species%iwf_i(j, is)
    l = species%iwf_l(j, is)
    m = species%iwf_m(j, is)
  end subroutine species_iwf_ilm
 
  ! ---------------------------------------------------------
  pure subroutine species_iwf_n(species, j, is, n)
    class(species_t), intent(in) :: species
    integer, intent(in)          :: j, is
    integer, intent(out)         :: n
 
    n = species%iwf_n(j, is)
  end subroutine species_iwf_n
 
  ! ---------------------------------------------------------
  pure subroutine species_iwf_j(species, iorb, j)
    class(species_t), intent(in)  :: species
    integer,          intent(in)  :: iorb
    real(real64),     intent(out) :: j
 
    j = species%iwf_j(iorb)
  end subroutine species_iwf_j
 
  ! ---------------------------------------------------------
  subroutine species_end(species)
    class(species_t), intent(inout) :: species
 
    push_sub(species_end)
 
    safe_deallocate_a(species%iwf_n)
    safe_deallocate_a(species%iwf_l)
    safe_deallocate_a(species%iwf_m)
    safe_deallocate_a(species%iwf_i)
    safe_deallocate_a(species%iwf_j)
 
    pop_sub(species_end)
  end subroutine species_end
  ! ---------------------------------------------------------
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Private procedures
 
  ! ---------------------------------------------------------
  function get_symbol(label) result(c_symbol)
    character(len=*), intent(in)    :: label
    character(kind=c_char)          :: c_symbol(label_len+1)
 
    character(len=LABEL_LEN) :: symbol
 
    integer :: ilend
 
    ! use only the first part of the label to determine the element
    do ilend = 1, len(label)
      if (iachar(label(ilend:ilend)) >= iachar('a') .and. iachar(label(ilend:ilend)) <= iachar('z')) cycle
      if (iachar(label(ilend:ilend)) >= iachar('A') .and. iachar(label(ilend:ilend)) <= iachar('Z')) cycle
      exit
    end do
    ilend = ilend - 1
 
    symbol = label(1:ilend)
 
    c_symbol = c_null_char
    c_symbol(1:c_str_len(symbol)) = string_f_to_c(symbol)
  end function get_symbol
 
  ! ---------------------------------------------------------
  logical function species_is_same_species(spec1, spec2) result(same)
    class(species_t), intent(in) :: spec1, spec2
 
    push_sub(species_is_same_species)
 
    same = same_type_as(spec1, spec2)
    if(abs(spec1%z - spec2%z) > m_epsilon) same = .false.
    if(abs(spec1%z_val - spec2%z_val) > m_epsilon) same = .false.
    if(abs(spec1%mass - spec2%mass) > m_epsilon) same = .false.
 
    pop_sub(species_is_same_species)
  end function species_is_same_species
 
  ! ---------------------------------------------------------
  logical pure function species_is_full(this)
    class(species_t),        intent(in) :: this
 
    species_is_full = .false.
  end function species_is_full
 
  ! ---------------------------------------------------------
  logical pure function species_is_ps(this)
    class(species_t),        intent(in) :: this
 
    species_is_ps = .false.
  end function species_is_ps
 
  ! ---------------------------------------------------------
  logical pure function species_is_ps_with_nlcc(this)
    class(species_t),        intent(in) :: this
 
    species_is_ps_with_nlcc = .false.
  end function species_is_ps_with_nlcc
 
  ! ---------------------------------------------------------
  logical pure function species_represents_real_atom(spec)
    class(species_t), intent(in) :: spec
 
    species_represents_real_atom = .false.
  end function species_represents_real_atom
 
  ! ---------------------------------------------------------
  logical pure function species_user_defined(spec)
    class(species_t), intent(in) :: spec
 
    species_user_defined = .false.
  end function species_user_defined
 
 
end module species_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: