xc_functional.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 xc_functional_oct_m
  use debug_oct_m
  use global_oct_m
  use libvdwxc_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use pseudo_oct_m
  use xc_f03_lib_m
#ifdef HAVE_LIBXC_FUNCS
  use xc_f03_funcs_m
#endif
 
  implicit none
 
  ! Although the following file only contain comments, we include it here to make sure it exists.
  ! Otherwise the code might compile, but not run properly, as the variables documentations
  ! will be incomplete.
#include "functionals_list.F90"
 
  private
  public ::                     &
    xc_functional_t,                &
    xc_functional_init,             &
    xc_functional_end,              &
    xc_functional_write_info,       &
    xc_functional_is_not_size_consistent, &
    xc_functional_is_energy_functional, &
    xc_get_default_functional
 
 
  integer, public, parameter ::   &
    XC_OEP_X = 901,               &  !< Exact exchange
    xc_oep_x_slater = 902,        &  
    xc_oep_x_fbe    = 903,        &  
    xc_ks_inversion = 904,        &  
    xc_rdmft_xc_m   = 905,        &  
    xc_oep_x_fbe_sl = 906,        &  
    xc_lda_c_fbe    = 907,        &  
    xc_lda_c_fbe_sl = 908,        &  
    xc_half_hartree = 917,        &  
    xc_vdw_c_vdwdf = 918,         &  
    xc_vdw_c_vdwdf2 = 919,        &  
    xc_vdw_c_vdwdfcx = 920,       &  
    xc_hyb_gga_xc_mvorb_hse06 = 921, &  
    xc_hyb_gga_xc_mvorb_pbeh = 922,  &  
    xc_mgga_x_nc_br = 923,           &  
    xc_mgga_x_nc_br_1 = 924,         &  
    xc_mgga_c_nc_cs = 925,           &  
    xc_mgga_x_nc_br_explicit = 926      
  !                                   !! with an explicit inversion of x(y), gamma = 0.8
 
  integer, public, parameter :: &
    XC_FAMILY_KS_INVERSION = 1024, &
    xc_family_rdmft = 2048, &
    xc_family_libvdwxc = 4096, &
    xc_family_nc_lda = 8192, &
    xc_family_nc_mgga = 16384
 
  type xc_functional_t
    ! Components are public by default
    integer         :: family = xc_family_unknown 
    integer         :: type   = 0
    integer         :: id     = 0
 
    integer         :: spin_channels = 0
    integer         :: flags         = 0
 
    logical, private :: from_libxc    = .false.
 
    type(xc_f03_func_t)               :: conf
    type(xc_f03_func_info_t), private :: info
    type(libvdwxc_t)                  :: libvdwxc
  end type xc_functional_t
 
  integer, public, parameter :: LIBXC_C_INDEX = 1000
 
contains
 
  ! ---------------------------------------------------------
  subroutine xc_functional_init(functl, namespace, id, ndim, nel, spin_channels)
    type(xc_functional_t), intent(inout) :: functl
    type(namespace_t),     intent(in)    :: namespace
    integer,               intent(in)    :: id
    integer,               intent(in)    :: ndim
    real(real64),          intent(in)    :: nel
    integer,               intent(in)    :: spin_channels
 
    integer :: interact_1d
    real(real64)   :: alpha, parameters(2)
 
    push_sub(xc_functional_init)
 
    ! initialize structure
    functl%id = id
    functl%spin_channels = spin_channels
 
    if (functl%id == 0) then
      functl%family = xc_family_none
    else
      ! get the family of the functional
      functl%family = xc_f03_family_from_id(functl%id)
      ! this also ensures it is actually a functional defined by the linked version of libxc
 
      if (functl%family == xc_family_unknown) then
 
        select case (functl%id)
        case (xc_oep_x, xc_oep_x_slater, xc_oep_x_fbe, xc_oep_x_fbe_sl, xc_lda_c_fbe_sl)
          functl%family = xc_family_oep
 
        case (xc_ks_inversion)
          functl%family = xc_family_ks_inversion
 
        case (xc_half_hartree)
          call messages_experimental("half-Hartree exchange", namespace=namespace)
          functl%family = xc_family_lda ! XXX not really
 
        case (xc_vdw_c_vdwdf, xc_vdw_c_vdwdf2, xc_vdw_c_vdwdfcx)
          functl%family = xc_family_libvdwxc
          !functl%flags = functl%flags + XC_FLAGS_HAVE_VXC + XC_FLAGS_HAVE_EXC
 
        case (xc_rdmft_xc_m)
          functl%family = xc_family_rdmft
 
        case (xc_hyb_gga_xc_mvorb_hse06, xc_hyb_gga_xc_mvorb_pbeh)
          functl%family = xc_family_hyb_gga
 
        case (xc_mgga_x_nc_br, xc_mgga_x_nc_br_1, xc_mgga_x_nc_br_explicit, xc_mgga_c_nc_cs)
          functl%family = xc_family_nc_mgga
 
        case default
          call messages_input_error(namespace, 'XCFunctional', 'Unknown functional')
 
        end select
      end if
    end if
 
    if (functl%family == xc_family_oep) then
      functl%type = xc_exchange
      functl%flags = xc_flags_1d + xc_flags_2d + xc_flags_3d
 
    else if (functl%family == xc_family_ks_inversion .or. functl%family == xc_family_rdmft) then
      functl%type = xc_exchange_correlation
      functl%flags = xc_flags_1d + xc_flags_2d + xc_flags_3d
 
    else if (functl%family == xc_family_libvdwxc) then
      call xc_f03_func_init(functl%conf, xc_lda_c_pw, spin_channels)
      functl%info = xc_f03_func_get_info(functl%conf)
      functl%type = xc_f03_func_info_get_kind(functl%info)
      functl%flags = xc_f03_func_info_get_flags(functl%info)
      ! Convert Octopus code for functional into corresponding libvdwxc code:
      call libvdwxc_init(functl%libvdwxc, namespace, functl%id - xc_vdw_c_vdwdf + 1)
 
    else if (functl%id == xc_half_hartree) then
      functl%type = xc_exchange_correlation
      functl%flags = xc_flags_1d + xc_flags_2d + xc_flags_3d
 
    else if(functl%id == xc_lda_c_fbe) then
      functl%family = xc_family_lda
      functl%type = xc_correlation
      functl%flags = xc_flags_have_exc + xc_flags_have_vxc + xc_flags_3d
 
    else if (functl%id == xc_mgga_x_nc_br .or. functl%id == xc_mgga_x_nc_br_1 .or. functl%id == xc_mgga_x_nc_br_explicit) then
      functl%type = xc_exchange
      functl%flags = xc_flags_have_vxc + xc_flags_have_exc + xc_flags_3d
 
    else if(functl%id == xc_mgga_c_nc_cs) then
      functl%type = xc_correlation
      functl%flags = xc_flags_have_vxc + xc_flags_have_exc + xc_flags_3d
 
    else if (functl%family == xc_family_none) then
      functl%type = -1
      functl%flags = xc_flags_1d + xc_flags_2d + xc_flags_3d
 
    else ! handled by libxc
      ! initialize
      functl%from_libxc = .true.
 
      !For the two MVORB functionals, we initialize libxc with the non-MVORB functionals
      select case (functl%id)
      case (xc_hyb_gga_xc_mvorb_hse06)
        call xc_f03_func_init(functl%conf, xc_hyb_gga_xc_hse06, spin_channels)
 
      case (xc_hyb_gga_xc_mvorb_pbeh)
        call xc_f03_func_init(functl%conf, xc_hyb_gga_xc_pbeh, spin_channels)
 
      case default
        call xc_f03_func_init(functl%conf, functl%id, spin_channels)
      end select
      functl%info     = xc_f03_func_get_info(functl%conf)
      functl%type     = xc_f03_func_info_get_kind(functl%info)
      functl%flags    = xc_f03_func_info_get_flags(functl%info)
 
      ! FIXME: no need to say this for kernel
      if (bitand(functl%flags, xc_flags_have_exc) == 0) then
        message(1) = 'Specified functional does not have total energy available.'
        message(2) = 'Corresponding component of energy will just be left as zero.'
        call messages_warning(2, namespace=namespace)
      end if
 
      if (bitand(functl%flags, xc_flags_have_vxc) == 0) then
        message(1) = 'Specified functional does not have XC potential available.'
        message(2) = 'Cannot run calculations. Choose another XCFunctional.'
        call messages_fatal(2, namespace=namespace)
      end if
    end if
 
    ! Octopus-defined functionals need to provide the proper flags
    if (functl%family /= xc_family_none) then
      call xc_check_dimension(functl, ndim, namespace)
    end if
 
    ! Yukawa hybrids not supported yet
    if (bitand(functl%flags, xc_flags_hyb_camy) /= 0) then
      call messages_not_implemented("Yukawa-range separated hybrids", namespace=namespace)
    end if
 
    ! FIXME: aren`t there other parameters that can or should be set?
    ! special parameters that have to be configured
    select case (functl%id)
      ! FIXME: aren`t there other Xalpha functionals?
    case (xc_lda_c_xalpha)
 
      !%Variable Xalpha
      !%Type float
      !%Default 1.0
      !%Section Hamiltonian::XC
      !%Description
      !% The parameter of the Slater X<math>\alpha</math> functional. Applies only for
      !% <tt>XCFunctional = xc_lda_c_xalpha</tt>.
      !%End
      call parse_variable(namespace, 'Xalpha', m_one, parameters(1))
 
      call xc_f03_func_set_ext_params(functl%conf, parameters(1))
 
      ! FIXME: doesn`t this apply to other 1D functionals?
    case (xc_lda_x_1d_soft, xc_lda_c_1d_csc)
      !%Variable Interaction1D
      !%Type integer
      !%Default interaction_soft_coulomb
      !%Section Hamiltonian::XC
      !%Description
      !% When running in 1D, one has to soften the Coulomb interaction. This softening
      !% is not unique, and several possibilities exist in the literature.
      !%Option interaction_exp_screened 0
      !% Exponentially screened Coulomb interaction.
      !% See, <i>e.g.</i>, M Casula, S Sorella, and G Senatore, <i>Phys. Rev. B</i> <b>74</b>, 245427 (2006).
      !%Option interaction_soft_coulomb 1
      !% Soft Coulomb interaction of the form <math>1/\sqrt{x^2 + \alpha^2}</math>.
      !%End
      call messages_obsolete_variable(namespace, 'SoftInteraction1D_alpha', 'Interaction1D')
      call parse_variable(namespace, 'Interaction1D', option__interaction1d__interaction_soft_coulomb, interact_1d)
 
      !%Variable Interaction1DScreening
      !%Type float
      !%Default 1.0
      !%Section Hamiltonian::XC
      !%Description
      !% Defines the screening parameter <math>\alpha</math> of the softened Coulomb interaction
      !% when running in 1D.
      !%End
      call messages_obsolete_variable(namespace, 'SoftInteraction1D_alpha', 'Interaction1DScreening')
      call parse_variable(namespace, 'Interaction1DScreening', m_one, alpha)
      parameters(1) = real(interact_1d, real64)
      parameters(2) = alpha
      call xc_f03_func_set_ext_params(functl%conf, parameters(1))
 
    case (xc_gga_x_lb)
      if (parse_is_defined(namespace, 'LB94_modified')) then
        call messages_obsolete_variable(namespace, 'LB94_modified')
      end if
 
      if (parse_is_defined(namespace, 'LB94_threshold')) then
        call messages_obsolete_variable(namespace, 'LB94_threshold')
      end if
    end select
 
 
    ! For functionals that depend on the number of electrons, we set the number of electrons
    ! This is currently the case of GGA_K_TFLW, LDA_X_2D_PRM, and LDA_X_RAE
    if (xc_functional_is_not_size_consistent(functl, namespace)) then
      assert(xc_f03_func_info_get_n_ext_params(functl%info) == 1)
      parameters(1) = nel
      call xc_f03_func_set_ext_params(functl%conf, parameters(1))
      write(message(1), '(a,i1)') "Info: Setting the number of electrons for the functional for spin ", spin_channels
      call messages_info(1, namespace=namespace)
    end if
 
    pop_sub(xc_functional_init)
  end subroutine xc_functional_init
 
 
  ! ---------------------------------------------------------
  subroutine xc_functional_end(functl)
    type(xc_functional_t), intent(inout) :: functl
 
    push_sub(xc_functional_end)
 
    if (functl%family /= xc_family_none .and. functl%family /= xc_family_oep .and. &
      functl%family /= xc_family_ks_inversion .and. functl%id /= xc_half_hartree &
      .and. functl%id /= xc_lda_c_fbe &
      .and. functl%family /= xc_family_nc_lda .and. functl%family /= xc_family_nc_mgga) then
      call xc_f03_func_end(functl%conf)
    end if
 
    if (functl%family == xc_family_libvdwxc) then
      call libvdwxc_end(functl%libvdwxc)
    end if
 
    pop_sub(xc_functional_end)
  end subroutine xc_functional_end
 
 
  ! ---------------------------------------------------------
  subroutine xc_functional_write_info(functl, iunit, namespace)
    type(xc_functional_t),       intent(in) :: functl
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    character(len=120) :: family
    integer :: ii
 
    push_sub(xc_functional_write_info)
 
    if (functl%family == xc_family_oep) then
      ! this is handled separately
 
      select case (functl%id)
      case (xc_oep_x)
        write(message(1), '(2x,a)') 'Exchange'
        write(message(2), '(4x,a)') 'Exact exchange'
        call messages_info(2, iunit=iunit, namespace=namespace)
 
      case(xc_oep_x_slater)
        write(message(1), '(2x,a)') 'Exchange'
        write(message(2), '(4x,a)') 'Slater exchange'
        call messages_info(2, iunit=iunit, namespace=namespace)
 
      case (xc_oep_x_fbe)
        write(message(1), '(2x,a)') 'Exchange'
        write(message(2), '(4x,a)') 'Force-based local exchange'
        write(message(3), '(4x,a)') '[1] Tancogne-Dejean et al., J. Chem. Phys. 160, 024103 (2024)'
        call messages_info(3, iunit=iunit, namespace=namespace)
 
      case (xc_oep_x_fbe_sl)
        write(message(1), '(2x,a)') 'Exchange'
        write(message(2), '(4x,a)') 'Force-based local exchange - Sturm-Liouville'
        call messages_info(2, iunit=iunit, namespace=namespace)
 
      case (xc_lda_c_fbe_sl)
        write(message(1), '(2x,a)') 'Correlation'
        write(message(2), '(4x,a)') 'Force-based local-density correlation - Sturm-Liouville'
        call messages_info(2, iunit=iunit, namespace=namespace)
 
      case default
        assert(.false.)
      end select
 
    else if (functl%family == xc_family_ks_inversion) then
      ! this is handled separately
      select case (functl%id)
      case (xc_ks_inversion)
        write(message(1), '(2x,a)') 'Exchange-Correlation:'
        write(message(2), '(4x,a)') '  KS Inversion'
        call messages_info(2, iunit=iunit, namespace=namespace)
      end select
 
    else if (functl%family == xc_family_libvdwxc) then
      call libvdwxc_write_info(functl%libvdwxc, iunit=iunit)
 
    else if (functl%id == xc_mgga_x_nc_br) then
      write(message(1), '(2x,a)') 'Exchange'
      write(message(2), '(4x,a)') 'Noncollinear Becke-Roussel (MGGA)'
      write(message(3), '(4x,a)') '[1] N. Tancogne-Dejean, A. Rubio, and C. A. Ullrich, Phys. Rev. B 107, 165111 (2023)'
      call messages_info(3, iunit)
 
    else if (functl%id == xc_mgga_x_nc_br_1) then
      write(message(1), '(2x,a)') 'Exchange'
      write(message(2), '(4x,a)') 'Noncollinear Becke-Roussel, gamma = 1.0 (MGGA)'
      write(message(3), '(4x,a)') '[1] N. Tancogne-Dejean, A. Rubio, and C. A. Ullrich, Phys. Rev. B 107, 165111 (2023)'
      call messages_info(3, iunit)
 
    else if (functl%id == xc_mgga_x_nc_br_explicit) then
      write(message(1), '(2x,a)') 'Exchange'
      write(message(2), '(4x,a)') 'Noncollinear Becke-Roussel (MGGA) with explicit inversion of x(y)'
      write(message(3), '(4x,a)') '[1] N. Tancogne-Dejean, A. Rubio, and C. A. Ullrich, Phys. Rev. B 107, 165111 (2023)'
      write(message(4), '(4x,a)') '[2] E. Proynov et al., Chemical Physics Letters 455, 103 (2008)'
      call messages_info(4, iunit)
 
    else if (functl%id == xc_mgga_c_nc_cs) then
      write(message(1), '(2x,a)') 'Correlation'
      write(message(2), '(4x,a)') 'Noncollinear Colle-Salvetti (MGGA)'
      write(message(3), '(4x,a)') '[1] N. Tancogne-Dejean, A. Rubio, and C. A. Ullrich, Phys. Rev. B 107, 165111 (2023)'
      call messages_info(3, iunit)
 
 
    else if (functl%id == xc_half_hartree) then
      write(message(1), '(2x,a)') 'Exchange-Correlation:'
      write(message(2), '(4x,a)') 'Half-Hartree two-electron exchange'
      call messages_info(2, iunit=iunit, namespace=namespace)
 
    else if(functl%id == xc_lda_c_fbe) then
      write(message(1), '(2x,a)') 'Correlation'
      write(message(2), '(4x,a)') 'Force-based LDA correlation'
      write(message(3), '(4x,a)') '[1] N. Tancogne-Dejean, M. Penz, M. Ruggenthaler, A. Rubio, J. Phys. Chem. A 129, 3132 (2025)'
      call messages_info(3, iunit=iunit, namespace=namespace)
 
    else if (functl%family /= xc_family_none) then ! all the other families
      select case (functl%type)
      case (xc_exchange)
        write(message(1), '(2x,a)') 'Exchange'
      case (xc_correlation)
        write(message(1), '(2x,a)') 'Correlation'
      case (xc_exchange_correlation)
        write(message(1), '(2x,a)') 'Exchange-correlation'
      case (xc_kinetic)
        call messages_not_implemented("kinetic-energy functionals", namespace=namespace)
      case default
        write(message(1), '(a,i6,a,i6)') "Unknown functional type ", functl%type, ' for functional ', functl%id
        call messages_fatal(1, namespace=namespace)
      end select
 
      select case (functl%family)
      case (xc_family_lda)
        write(family,'(a)') "LDA"
      case (xc_family_gga)
        write(family,'(a)') "GGA"
      case (xc_family_mgga)
        write(family,'(a)') "MGGA"
      case (xc_family_lca)
        call messages_not_implemented("Support of LCA functional", namespace=namespace)
      case (xc_family_hyb_lda)
        write(family,'(a)') "Hybrid LDA"
      case (xc_family_hyb_gga)
        write(family,'(a)') "Hybrid GGA"
      case (xc_family_hyb_mgga)
        write(family,'(a)') "Hybrid MGGA"
      end select
      write(message(2), '(4x,4a)') trim(xc_f03_func_info_get_name(functl%info)), ' (', trim(family), ')'
      call messages_info(2, iunit=iunit, namespace=namespace)
 
      ii = 0
      do while(ii >= 0)
        write(message(1), '(4x,a,i1,2a)') '[', ii + 1, '] ', &
          trim(xc_f03_func_reference_get_ref(xc_f03_func_info_get_references(functl%info, ii)))
        call messages_info(1, iunit=iunit, namespace=namespace)
      end do
    end if
 
    pop_sub(xc_functional_write_info)
  end subroutine xc_functional_write_info
 
 
  integer function xc_get_default_functional(dim, pseudo_x_functional, pseudo_c_functional) result(default)
    integer, intent(in) :: dim
    integer, intent(in) :: pseudo_x_functional, pseudo_c_functional
 
    push_sub(xc_get_default_functional)
 
    default = 0
 
    if (pseudo_x_functional /= pseudo_exchange_any) then
      default = pseudo_x_functional
    else
      select case (dim)
      case (3)
        default = xc_lda_x
      case (2)
        default = xc_lda_x_2d
      case (1)
        default = xc_lda_x_1d_soft
      end select
    end if
 
    assert(default >= 0)
 
    if (pseudo_c_functional /= pseudo_correlation_any) then
      default = default + libxc_c_index*pseudo_c_functional
    else
      select case (dim)
      case (3)
        default = default + libxc_c_index * xc_lda_c_pz_mod
      case (2)
        default = default + libxc_c_index * xc_lda_c_2d_amgb
      case (1)
        default = default + libxc_c_index * xc_lda_c_1d_csc
      end select
    end if
 
    assert(default >= 0)
 
    pop_sub(xc_get_default_functional)
  end function xc_get_default_functional
 
  subroutine xc_check_dimension(functl, ndim, namespace)
    type(xc_functional_t), intent(in) :: functl
    integer,               intent(in) :: ndim
    type(namespace_t),     intent(in) :: namespace
 
    logical :: ok
 
    push_sub(xc_check_dimension)
 
    ! Check that the functional is used for the correct space dimension
    ok = bitand(functl%flags, xc_flags_1d) /= 0
    if (ndim == 1 .and. (.not. ok)) then
      message(1) = 'Cannot use the specified functionals in 1D.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    ok = bitand(functl%flags, xc_flags_2d) /= 0
    if (ndim == 2 .and. (.not. ok)) then
      message(1) = 'Cannot use the specified functionals in 2D.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    ok = bitand(functl%flags, xc_flags_3d) /= 0
    if (ndim == 3 .and. (.not. ok)) then
      message(1) = 'Cannot use the specified functionals in 3D.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub(xc_check_dimension)
  end subroutine xc_check_dimension
 
  logical function xc_functional_is_not_size_consistent(functl, namespace)
    type(xc_functional_t), intent(in) :: functl
    type(namespace_t),     intent(in) :: namespace
 
    integer :: n_ext_params, ip
    character(len=128) :: ext_params_name
 
    xc_functional_is_not_size_consistent = .false.
    if (.not. functl%from_libxc) return
 
    push_sub(xc_functional_is_not_size_consistent)
 
    ! External parameters
    n_ext_params = xc_f03_func_info_get_n_ext_params(functl%info)
    do ip = 0, n_ext_params-1
      ext_params_name = xc_f03_func_info_get_ext_params_name(functl%info, ip)
      ! Internal parameters starts with '_'
      if (ext_params_name(1:1) == '_') cycle
      if (trim(xc_f03_func_info_get_ext_params_description(functl%info, ip)) == 'Number of electrons') then
        xc_functional_is_not_size_consistent = .true.
      end if
    end do
 
    ! We only support size inconsistent functionals with a single parameter at the moment
    if (xc_functional_is_not_size_consistent .and. n_ext_params > 1) then
      message(1) = 'The selected functional is currently not supported.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub(xc_functional_is_not_size_consistent)
  end function xc_functional_is_not_size_consistent
 
  logical pure function xc_functional_is_energy_functional(functl)
    type(xc_functional_t), intent(in) :: functl
 
    xc_functional_is_energy_functional = bitand(functl%flags, xc_flags_have_exc) /= 0
  end function xc_functional_is_energy_functional
 
end module xc_functional_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: