xc.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"
 
#if defined(HAVE_CUDA)
#define HAVE_ACCEL 1
#endif
 
module xc_oct_m
  use accel_oct_m
  use allocate_hardware_aware_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use distributed_oct_m
  use electron_space_oct_m
  use exchange_operator_oct_m
  use global_oct_m
  use grid_oct_m
  use io_oct_m
  use io_function_oct_m
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lalg_basic_oct_m
  use lalg_adv_oct_m
  use libvdwxc_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use types_oct_m
  use unit_system_oct_m
  use xc_cam_oct_m
  use xc_f03_lib_m
#ifdef HAVE_LIBXC_FUNCS
  use xc_f03_funcs_m
#endif
  use xc_fbe_oct_m
  use xc_functional_oct_m
  use xc_interaction_oct_m
  use xc_lrc_oct_m
 
  implicit none
 
  private
  public ::             &
    xc_t,               &
    internal_quantities_t, &
    xc_init,            &
    xc_init_device_support, &
    xc_is_using_device, &
    xc_end,             &
    xc_write_info,      &
    xc_write_fxc_info,  &
    xc_is_orbital_dependent, &
    xc_is_not_size_consistent, &
    xc_is_energy_functional, &
    family_is_gga,      &
    family_is_supported, &
    family_is_mgga,     &
    family_is_nc_mgga,  &
    family_is_mgga_with_exc, &
    family_is_hybrid,   &
    in_family
 
  ! A Structure that contains the quantities needed to compute the functionals
  type internal_quantities_t
    real(real64), pointer     :: rho(:,:)       ! A pointer to the full density
 
    real(real64), allocatable :: dens(:,:)      ! Density (in the local frame of the magnetization)
    real(real64), allocatable :: gdens(:,:,:)   ! Gradient of the density
    real(real64), allocatable :: ldens(:,:)     ! Laplacian of the density
    real(real64), allocatable :: tau(:,:)       ! Kinetic energy density
 
    type(accel_mem_t) :: dens_buff              ! GPU buffer: extracted density in libxc layout (spin_channels x np)
    type(accel_mem_t) :: zk_buff                ! GPU buffer: energy per unit volume
    type(accel_mem_t) :: dedd_buff              ! GPU buffer: dE/dn
    type(accel_mem_t) :: sigma_buff             ! GPU buffer: contracted density gradient (1 or 3 x np)
    type(accel_mem_t) :: vsigma_buff            ! GPU buffer: dE/dsigma
    type(accel_mem_t) :: tau_buff               ! GPU buffer: kinetic energy density (spin_channels x np)
    type(accel_mem_t) :: ldens_buff             ! GPU buffer: Laplacian of the density (spin_channels x np)
    type(accel_mem_t) :: dedtau_buff            ! GPU buffer: dE/dtau
    type(accel_mem_t) :: dedldens_buff          ! GPU buffer: dE/dlaplacian
 
    ! Persistent pinned (page-locked) host staging buffers in libxc layout.
    real(real64), pointer, contiguous :: hbuf_dens(:,:)     => null()
    real(real64), pointer, contiguous :: hbuf_sigma(:,:)    => null()
    real(real64), pointer, contiguous :: hbuf_ldens(:,:)    => null()
    real(real64), pointer, contiguous :: hbuf_tau(:,:)      => null()
    real(real64), pointer, contiguous :: hbuf_zk(:)         => null()
    real(real64), pointer, contiguous :: hbuf_dedd(:,:)     => null()
    real(real64), pointer, contiguous :: hbuf_vsigma(:,:)   => null()
    real(real64), pointer, contiguous :: hbuf_dedldens(:,:) => null()
    real(real64), pointer, contiguous :: hbuf_dedtau(:,:)   => null()
 
    ! Density correction buffers: populated before XC if corrections are needed on GPU.
    type(accel_mem_t) :: buff_rho_core          ! GPU copy of st%rho_core (np)
    type(accel_mem_t) :: buff_frozen_rho        ! GPU copy of st%frozen_rho (np x nspin)
    integer           :: frozen_rho_np = 0      ! leading dimension of frozen_rho on CPU (= gr%np)
    real(real64)      :: amaldi_factor = m_one  ! Amaldi scaling factor applied to dens_buff (1.0 if inactive)
    ! True when dens_buff was pre-populated from st%buff_density by xc_update_internal_quantities.
    ! False when xc_compute_vxc must upload CPU l_dens to dens_buff each block.
    logical           :: dens_buff_prefilled = .false.
  end type internal_quantities_t
 
 
  type xc_t
    private
    integer,               public :: family
    integer,               public :: flags
    integer,               public :: kernel_family
    type(xc_functional_t), public :: functional(2,2)
 
    type(xc_functional_t), public :: kernel(2,2)
    type(xc_cam_t),        public :: cam
    type(xc_lrc_t),        public :: lrc
 
    logical,               public :: use_gi_ked
    integer,               public :: xc_density_correction
    logical,               public :: xcd_optimize_cutoff
    real(real64),          public :: xcd_ncutoff
    logical,               public :: xcd_minimum
    logical,               public :: xcd_normalize
    logical,               public :: parallel
    logical,               public :: xc_on_host
 
    type(internal_quantities_t), public :: quantities
 
  contains
    procedure :: compute_exchange => xc_compute_exchange
 
  end type xc_t
 
  real(real64), public, parameter :: xc_tiny = 1.0e-12_real64
 
  integer, public, parameter :: &
    LR_NONE = 0,        &
    lr_x    = 1
 
contains
 
  ! ---------------------------------------------------------
  subroutine xc_write_info(xcs, iunit, namespace)
    type(xc_t),                  intent(in) :: xcs
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    integer :: ifunc
 
    push_sub(xc_write_info)
 
    write(message(1), '(a)') "Exchange-correlation:"
    call messages_info(1, iunit=iunit, namespace=namespace)
 
    do ifunc = func_x, func_c
      call xc_functional_write_info(xcs%functional(ifunc, 1), iunit, namespace)
    end do
 
    if (abs(xcs%cam%alpha) + abs(xcs%cam%beta) > m_epsilon) then
      write(message(1), '(1x)')
      write(message(2), '(a,f8.5)') "Exact exchange mixing = ", xcs%cam%alpha
      write(message(3), '(a,f8.5)') "Exact exchange for short-range beta = ", xcs%cam%beta
      write(message(4), '(a,f8.5)') "Exact exchange range-separate omega = ", xcs%cam%omega
      call messages_info(4, iunit=iunit, namespace=namespace)
    end if
 
 
    pop_sub(xc_write_info)
  end subroutine xc_write_info
 
  ! ---------------------------------------------------------
  subroutine xc_write_fxc_info(xcs, iunit, namespace)
    type(xc_t),                  intent(in) :: xcs
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    integer :: ifunc
 
    push_sub(xc_write_fxc_info)
 
    write(message(1), '(a)') "Exchange-correlation kernel:"
    call messages_info(1, iunit=iunit, namespace=namespace)
 
    do ifunc = func_x, func_c
      call xc_functional_write_info(xcs%kernel(ifunc, 1), iunit, namespace)
    end do
 
    pop_sub(xc_write_fxc_info)
  end subroutine xc_write_fxc_info
 
 
  ! ---------------------------------------------------------
  subroutine xc_init_device_support(xcs, namespace)
    type(xc_t),        intent(inout) :: xcs
    type(namespace_t), intent(in)    :: namespace
 
    push_sub(xc_init_device_support)
 
    xcs%xc_on_host = .not. xc_is_using_device(namespace)
 
    pop_sub(xc_init_device_support)
  end subroutine xc_init_device_support
 
  ! ---------------------------------------------------------
  logical function xc_is_using_device(namespace)
    type(namespace_t), intent(in) :: namespace
 
    logical :: xc_force_cpu
 
    push_sub(xc_is_using_device)
 
    !%Variable XCForceFunctionalsOnCPU
    !%Type logical
    !%Default no
    !%Section Hamiltonian::XC
    !%Description
    !% By setting this variable to <tt>yes</tt>, you force Octopus to evaluate XC functionals
    !% on the CPU, even if CUDA/HIP support is available.
    !%End
    call parse_variable(namespace, 'XCForceFunctionalsOnCPU', .false., xc_force_cpu)
 
#if defined(HAVE_LIBXC_DEVICE) && defined(HAVE_ACCEL)
    xc_is_using_device = .false.
    if (accel_is_enabled()) then
      xc_is_using_device = .not. xc_force_cpu
    end if
#elif defined(HAVE_LIBXC_DEVICE) && !defined(HAVE_ACCEL)
    if(.not. xc_force_cpu) then
      call messages_warning("Libxc has device support but no accelerator support is enabled in Octopus.", namespace=namespace)
    end if
    xc_is_using_device = .false.
#else
    xc_is_using_device = .false.
#endif
 
    pop_sub(xc_is_using_device)
  end function xc_is_using_device
 
  ! ---------------------------------------------------------
  subroutine xc_init(xcs, namespace, ndim, periodic_dim, nel, x_id, c_id, xk_id, ck_id, hartree_fock, ispin)
    type(xc_t),        intent(out) :: xcs
    type(namespace_t), intent(in)  :: namespace
    integer,           intent(in)  :: ndim
    integer,           intent(in)  :: periodic_dim
    real(real64),      intent(in)  :: nel
    integer,           intent(in)  :: x_id
    integer,           intent(in)  :: c_id
    integer,           intent(in)  :: xk_id
    integer,           intent(in)  :: ck_id
    logical,           intent(in)  :: hartree_fock
    integer,           intent(in)  :: ispin
 
    integer :: isp, xc_major, xc_minor, xc_micro
    logical :: ll
    type(block_t)   :: blk
    type(xc_cam_t)  :: cam_ext
 
    push_sub(xc_init)
 
    call xc_f03_version(xc_major, xc_minor, xc_micro)
 
    xcs%family = 0
    xcs%flags  = 0
    xcs%kernel_family = 0
 
    call parse()
 
    call xc_lrc_init(xcs%lrc, namespace, ndim, periodic_dim)
 
    call xc_init_device_support(xcs, namespace)
 
    ! 1D/2D functionals (e.g. lda_x_1d_soft) call xc_integrate() which uses malloc() inside
    ! GPU kernels. With large meshes this exhausts the CUDA device heap (cudaErrorIllegalAddress).
    if (ndim /= 3) xcs%xc_on_host = .true.
 
    !we also need XC functionals that do not depend on the current
    !get both spin-polarized and unpolarized
 
    ! TODO: check that nel should not be spin polarized here
    do isp = 1, 2
 
      call xc_functional_init(xcs%functional(func_x, isp), namespace, x_id, ndim, nel, isp, xcs%xc_on_host)
      call xc_functional_init(xcs%functional(func_c, isp), namespace, c_id, ndim, nel, isp, xcs%xc_on_host)
 
      call xc_functional_init(xcs%kernel(func_x, isp), namespace, xk_id, ndim, nel, isp, xcs%xc_on_host)
      call xc_functional_init(xcs%kernel(func_c, isp), namespace, ck_id, ndim, nel, isp, xcs%xc_on_host)
 
    end do
 
    xcs%family = ior(xcs%family, xcs%functional(func_x,1)%family)
    xcs%family = ior(xcs%family, xcs%functional(func_c,1)%family)
 
    xcs%flags = ior(xcs%flags, xcs%functional(func_x,1)%flags)
    xcs%flags = ior(xcs%flags, xcs%functional(func_c,1)%flags)
 
    xcs%kernel_family = ior(xcs%kernel_family, xcs%kernel(func_x,1)%family)
    xcs%kernel_family = ior(xcs%kernel_family, xcs%kernel(func_c,1)%family)
 
 
    if (xc_is_not_size_consistent(xcs, namespace) .and. periodic_dim > 0) then
      message(1) = "Cannot perform a periodic calculation with a functional"
      message(2) = "that depends on the number of electrons."
      call messages_fatal(2, namespace=namespace)
    end if
 
    ! Take care of hybrid functionals (they appear in the correlation functional)
    xcs%cam = cam_null
 
    ll =  (hartree_fock) &
      .or.(xcs%functional(func_x,1)%id == xc_oep_x) &
      .or. family_is_hybrid(xcs)
    if (ll) then
      if ((xcs%functional(func_x,1)%id /= 0).and.(xcs%functional(func_x,1)%id /= xc_oep_x)) then
        message(1) = "You cannot use an exchange functional when performing"
        message(2) = "a Hartree-Fock calculation or using a hybrid functional."
        call messages_fatal(2, namespace=namespace)
      end if
 
      if (periodic_dim == ndim) then
        call messages_experimental("Fock operator (Hartree-Fock, OEP, hybrids) in fully periodic systems", namespace=namespace)
      end if
 
      ! get the mixing coefficient for hybrids
      if (family_is_hybrid(xcs)) then
        if( .not. cam_ext%is_null() ) call set_hybrid_params(xcs, namespace, cam_ext)
        call xc_f03_hyb_cam_coef(xcs%functional(func_c,1)%conf, xcs%cam%omega, &
          xcs%cam%alpha, xcs%cam%beta)
        call xc_f03_hyb_cam_coef(xcs%functional(func_c,2)%conf, xcs%cam%omega, &
          xcs%cam%alpha, xcs%cam%beta)
      else
        ! we are doing Hartree-Fock plus possibly a correlation functional
        xcs%cam = cam_exact_exchange
      end if
 
      ! reset certain variables
      xcs%functional(func_x,1)%family = xc_family_oep
      xcs%functional(func_x,1)%id     = xc_oep_x
      xcs%functional(func_x,2)%family = xc_family_oep
      xcs%functional(func_x,2)%id     = xc_oep_x
      if (.not. hartree_fock) then
        xcs%family             = ior(xcs%family, xc_family_oep)
      end if
    end if
 
    if (in_family(xcs%family, [xc_family_lca])) then
      call messages_not_implemented("LCA current functionals", namespace) ! not even in libxc!
    end if
 
    call messages_obsolete_variable(namespace, 'MGGAimplementation')
    call messages_obsolete_variable(namespace, 'CurrentInTau', 'XCUseGaugeIndependentKED')
 
    if (xcs%functional(func_x, 1)%id == xc_mgga_x_tb09 .and. periodic_dim /= 3) then
      message(1) = "mgga_x_tb09 functional can only be used for 3D periodic systems"
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (family_is_mgga(xcs%family) .or. family_is_nc_mgga(xcs%family)) then
      !%Variable XCUseGaugeIndependentKED
      !%Type logical
      !%Default yes
      !%Section Hamiltonian::XC
      !%Description
      !% If true, when evaluating the XC functional, a term including the (paramagnetic or total) current
      !% is added to the kinetic-energy density such as to make it gauge-independent.
      !% Applies only to meta-GGA (and hybrid meta-GGA) functionals.
      !%End
      call parse_variable(namespace, 'XCUseGaugeIndependentKED', .true., xcs%use_gi_ked)
    end if
 
    pop_sub(xc_init)
 
  contains
 
    subroutine parse()
 
      push_sub(xc_init.parse)
 
      ! the values of x_id,  c_id, xk_id, and c_id are read outside the routine
 
      !%Variable XCDensityCorrection
      !%Type integer
      !%Default none
      !%Section Hamiltonian::XC::DensityCorrection
      !%Description
      !% This variable controls the long-range correction of the XC
      !% potential using the <a href=http://arxiv.org/abs/1107.4339>XC density representation</a>.
      !%Option none 0
      !% No correction is applied.
      !%Option long_range_x 1
      !% The correction is applied to the exchange potential.
      !%End
      call parse_variable(namespace, 'XCDensityCorrection', lr_none, xcs%xc_density_correction)
 
      if (xcs%xc_density_correction /= lr_none) then
        call messages_experimental('XC density correction', namespace=namespace)
 
        if(ispin /= unpolarized) then
          call messages_not_implemented('XCDensityCorrection with SpinComponents /= unpolarized', namespace)
        end if
 
        !%Variable XCDensityCorrectionOptimize
        !%Type logical
        !%Default true
        !%Section Hamiltonian::XC::DensityCorrection
        !%Description
        !% When enabled, the density cutoff will be
        !% optimized to replicate the boundary conditions of the exact
        !% XC potential. If the variable is set to no, the value of
        !% the cutoff must be given by <tt>XCDensityCorrectionCutoff</tt>
        !% variable.
        !%End
        call parse_variable(namespace, 'XCDensityCorrectionOptimize', .true., xcs%xcd_optimize_cutoff)
 
        !%Variable XCDensityCorrectionCutoff
        !%Type float
        !%Default 0.0
        !%Section Hamiltonian::XC::DensityCorrection
        !%Description
        !% The value of the cutoff applied to the XC density.
        !%End
        call parse_variable(namespace, 'XCDensityCorrectionCutoff', m_zero, xcs%xcd_ncutoff)
 
        !%Variable XCDensityCorrectionMinimum
        !%Type logical
        !%Default true
        !%Section Hamiltonian::XC::DensityCorrection
        !%Description
        !% When enabled, the cutoff optimization will
        !% return the first minimum of the <math>q_{xc}</math> function if it does
        !% not find a value of -1 (<a href=http://arxiv.org/abs/1107.4339>details</a>).
        !% This is required for atoms or small
        !% molecules, but may cause numerical problems.
        !%End
        call parse_variable(namespace, 'XCDensityCorrectionMinimum', .true., xcs%xcd_minimum)
 
        !%Variable XCDensityCorrectionNormalize
        !%Type logical
        !%Default true
        !%Section Hamiltonian::XC::DensityCorrection
        !%Description
        !% When enabled, the correction will be
        !% normalized to reproduce the exact boundary conditions of
        !% the XC potential.
        !%End
        call parse_variable(namespace, 'XCDensityCorrectionNormalize', .true., xcs%xcd_normalize)
 
      end if
 
      !%Variable ParallelXC
      !%Type logical
      !%Default true
      !%Section Execution::Parallelization
      !%Description
      !% When enabled, additional parallelization
      !% will be used for the calculation of the XC functional.
      !%End
      call messages_obsolete_variable(namespace, 'XCParallel', 'ParallelXC')
      call parse_variable(namespace, 'ParallelXC', .true., xcs%parallel)
 
 
      !%Variable HybridCAMParameters
      !%Type block
      !%Section Hamiltonian::XC
      !%Description
      !% This variable specifies the <math>\alpha, \beta, \omega</math> for CAM-type
      !% hybrid functionals. Defaults are zero.
      !%End
      cam_ext = cam_null
 
      if(parse_block(namespace, 'HybridCamParameters', blk) == 0) then
        call parse_block_float(blk, 0, 0, cam_ext%alpha)
        call parse_block_float(blk, 0, 1, cam_ext%beta)
        call parse_block_float(blk, 0, 2, cam_ext%omega)
        call parse_block_end(blk)
      end if
 
      if(.not. cam_ext%is_null()) then
        call cam_ext%print(namespace, msg="Info: Setting external CAM parameters")
      endif
 
      pop_sub(xc_init.parse)
    end subroutine parse
 
  end subroutine xc_init
 
 
  ! ---------------------------------------------------------
  subroutine xc_end(xcs)
    type(xc_t), intent(inout) :: xcs
 
    integer :: isp
 
    push_sub(xc_end)
 
    do isp = 1, 2
      call xc_functional_end(xcs%functional(func_x, isp))
      call xc_functional_end(xcs%functional(func_c, isp))
      call xc_functional_end(xcs%kernel(func_x, isp))
      call xc_functional_end(xcs%kernel(func_c, isp))
    end do
    xcs%family = 0
    xcs%flags  = 0
 
    ! Free the persistent host staging buffers (see internal_quantities_t), allocated
    ! with the hardware-aware allocator in xc_compute_vxc.
    if (associated(xcs%quantities%hbuf_dens)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_dens(1,1)), size(xcs%quantities%hbuf_dens, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_dens)
    end if
    if (associated(xcs%quantities%hbuf_sigma)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_sigma(1,1)), size(xcs%quantities%hbuf_sigma, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_sigma)
    end if
    if (associated(xcs%quantities%hbuf_ldens)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_ldens(1,1)), size(xcs%quantities%hbuf_ldens, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_ldens)
    end if
    if (associated(xcs%quantities%hbuf_tau)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_tau(1,1)), size(xcs%quantities%hbuf_tau, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_tau)
    end if
    if (associated(xcs%quantities%hbuf_zk)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_zk(1)), size(xcs%quantities%hbuf_zk, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_zk)
    end if
    if (associated(xcs%quantities%hbuf_dedd)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_dedd(1,1)), size(xcs%quantities%hbuf_dedd, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_dedd)
    end if
    if (associated(xcs%quantities%hbuf_vsigma)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_vsigma(1,1)), size(xcs%quantities%hbuf_vsigma, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_vsigma)
    end if
    if (associated(xcs%quantities%hbuf_dedldens)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_dedldens(1,1)), &
        size(xcs%quantities%hbuf_dedldens, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_dedldens)
    end if
    if (associated(xcs%quantities%hbuf_dedtau)) then
      call deallocate_hardware_aware(c_loc(xcs%quantities%hbuf_dedtau(1,1)), size(xcs%quantities%hbuf_dedtau, kind=int64)*8_int64)
      nullify(xcs%quantities%hbuf_dedtau)
    end if
 
    pop_sub(xc_end)
  end subroutine xc_end
 
  ! ---------------------------------------------------------
  logical pure function xc_is_orbital_dependent(xcs)
    type(xc_t), intent(in) :: xcs
 
    xc_is_orbital_dependent = family_is_hybrid(xcs) .or. &
      in_family(xcs%functional(func_x,1)%family, [xc_family_oep]) .or. &
      in_family(xcs%family, [xc_family_mgga, xc_family_nc_mgga])
 
  end function xc_is_orbital_dependent
 
  ! ---------------------------------------------------------
  pure logical function family_is_gga(family, only_collinear)
    integer,           intent(in) :: family
    logical, optional, intent(in) :: only_collinear
 
    if(optional_default(only_collinear, .false.)) then
      family_is_gga = in_family(family, [xc_family_gga, xc_family_hyb_gga, &
        xc_family_mgga, xc_family_hyb_mgga, xc_family_libvdwxc])
    else
      family_is_gga = in_family(family, [xc_family_gga, xc_family_hyb_gga, &
        xc_family_mgga, xc_family_hyb_mgga, xc_family_libvdwxc, xc_family_nc_mgga])
    end if
  end function  family_is_gga
 
  !----------------------------------------------------------------------
  pure logical function family_is_supported(family)
    integer, intent(in) :: family
 
    family_is_supported = in_family(family, [xc_family_lda, xc_family_hyb_lda, xc_family_gga, xc_family_hyb_gga, &
      xc_family_mgga, xc_family_hyb_mgga, xc_family_libvdwxc])
  end function  family_is_supported
 
  ! ---------------------------------------------------------
  pure logical function family_is_mgga(family, only_collinear)
    integer,           intent(in) :: family
    logical, optional, intent(in) :: only_collinear
 
    if(optional_default(only_collinear, .false.)) then
      family_is_mgga = in_family(family, [xc_family_mgga, xc_family_hyb_mgga])
    else
      family_is_mgga = in_family(family, [xc_family_mgga, xc_family_hyb_mgga, xc_family_nc_mgga])
    end if
  end function family_is_mgga
 
  pure logical function family_is_nc_mgga(family)
    integer, intent(in) :: family
 
    family_is_nc_mgga = bitand(family, xc_family_nc_mgga) /= 0
  end function family_is_nc_mgga
 
  ! ---------------------------------------------------------
  logical pure function family_is_mgga_with_exc(xcs)
    type(xc_t), intent(in) :: xcs
 
    integer :: ixc
 
    family_is_mgga_with_exc = .false.
    do ixc = 1, 2
      if (in_family(xcs%functional(ixc, 1)%family, [xc_family_mgga, xc_family_hyb_mgga, xc_family_nc_mgga]) &
        .and. xc_functional_is_energy_functional(xcs%functional(ixc, 1))) then
        family_is_mgga_with_exc = .true.
      end if
    end do
  end function family_is_mgga_with_exc
 
  logical pure function family_is_hybrid(xcs)
    type(xc_t), intent(in) :: xcs
 
    integer :: ixc
 
    family_is_hybrid = .false.
    do ixc = 1, 2
      if (in_family(xcs%functional(ixc, 1)%family, [xc_family_hyb_lda, xc_family_hyb_gga, xc_family_hyb_mgga])) then
        family_is_hybrid = .true.
      end if
    end do
  end function family_is_hybrid
 
  pure logical function in_family(family, xc_families)
    integer, intent(in) :: family
    integer, intent(in) :: xc_families(:)
 
    in_family = bitand(family, sum(xc_families)) /= 0
  end function in_family
 
  logical function xc_compute_exchange(xc, theory_level)
    class(xc_t), intent(in) :: xc
    integer,     intent(in) :: theory_level
 
    integer, parameter  :: exchange_theory_level(3) = [hartree, hartree_fock, rdmft]
 
    push_sub(xc_compute_exchange)
 
    xc_compute_exchange = any(exchange_theory_level == theory_level) .or. &
      (theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(xc)) .or. &
      (xc%functional(func_x, 1)%id == xc_oep_x_slater) .or. &
      (bitand(xc%family, xc_family_oep) /= 0)
 
    pop_sub(xc_compute_exchange)
 
  end function xc_compute_exchange
 
  ! ---------------------------------------------------------
  subroutine set_hybrid_params(xcs, namespace, cam_ext)
    type(xc_t),        intent(inout)  :: xcs
    type(namespace_t), intent(in)     :: namespace
    type(xc_cam_t),    intent(in)     :: cam_ext
 
    real(real64), parameter :: default_alpha_pbe0 = 0.25_real64 
    real(real64)            :: parameters(3)
 
    push_sub(set_hybrid_params)
 
    ! LibXC expects an array of reals, ordered [alpha, beta, omega]
    parameters = cam_ext%as_array()
    xcs%cam%alpha = parameters(1)
 
    select case(xcs%functional(func_c, 1)%id)
 
    case(xc_hyb_gga_xc_pbeh, xc_hyb_lda_xc_lda0) ! original PBE0/LDA0 in libxc
      if(parameters(1) < m_zero) parameters(1) = default_alpha_pbe0
      call xc_f03_func_set_ext_params(xcs%functional(func_c, 1)%conf, parameters)
      call xc_f03_func_set_ext_params(xcs%functional(func_c, 2)%conf, parameters)
      write(message(1), '(a,f6.3,a)') 'Info: Setting mixing parameter (' , parameters(1) ,').'
      call messages_info(1)
 
    case(xc_hyb_gga_xc_cam_pbeh, xc_hyb_lda_xc_cam_lda0)
      xcs%cam%beta = parameters(2)
      xcs%cam%omega = parameters(3)
      call xc_f03_func_set_ext_params(xcs%functional(func_c, 1)%conf, parameters)
      call xc_f03_func_set_ext_params(xcs%functional(func_c, 2)%conf, parameters)
      ! check parameters
      call xc_f03_hyb_cam_coef(xcs%functional(func_c,1)%conf, xcs%cam%omega, &
        xcs%cam%alpha, xcs%cam%beta)
      call xc_f03_hyb_cam_coef(xcs%functional(func_c,2)%conf, xcs%cam%omega, &
        xcs%cam%alpha, xcs%cam%beta)
      call xcs%cam%print(namespace, msg="Setting CAM parameters:")
 
    case default
      assert(.false.)
 
    end select
 
    pop_sub(set_hybrid_params)
  end subroutine set_hybrid_params
 
  ! ---------------------------------------------------------
  logical function xc_is_not_size_consistent(xcs, namespace)
    type(xc_t),        intent(in) :: xcs
    type(namespace_t), intent(in) :: namespace
 
    xc_is_not_size_consistent = xc_functional_is_not_size_consistent(xcs%functional(func_x,1), namespace) &
      .or. xc_functional_is_not_size_consistent(xcs%functional(func_c,1), namespace)
  end function xc_is_not_size_consistent
 
  ! ---------------------------------------------------------
  logical pure function xc_is_energy_functional(xcs)
    type(xc_t), intent(in)  :: xcs
    xc_is_energy_functional = xc_functional_is_energy_functional(xcs%functional(func_x,1)) &
      .or. xc_functional_is_energy_functional(xcs%functional(func_c,1))
  end function xc_is_energy_functional
 
end module xc_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: