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"
 
module xc_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 unit_system_oct_m
  use xc_f03_lib_m
  use xc_fbe_oct_m
  use xc_functional_oct_m
  use xc_interaction_oct_m
 
  implicit none
 
  private
  public ::             &
    xc_t,               &
    xc_init,            &
    xc_end,             &
    xc_write_info,      &
    xc_get_vxc,         &
    xc_get_nc_vxc,      &
    xc_get_fxc,         &
    xc_get_kxc,         &
    xc_is_orbital_dependent, &
    family_is_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
  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)
    real(real64),          public :: kernel_lrc_alpha
 
    real(real64), public   :: cam_omega
    real(real64), public   :: cam_alpha
    real(real64), public   :: cam_beta
    real(real64), public   :: cam_ext(3)
 
    logical         :: use_gi_ked
 
    integer :: xc_density_correction
    logical :: xcd_optimize_cutoff
    real(real64)   :: xcd_ncutoff
    logical :: xcd_minimum
    logical :: xcd_normalize
    logical :: parallel
 
    type(internal_quantities_t) :: quantities
  end type xc_t
 
  real(real64), parameter :: tiny      = 1.0e-12_real64
 
  integer, 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 + 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_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
 
    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()
 
    !we also need XC functionals that do not depend on the current
    !get both spin-polarized and unpolarized
    do isp = 1, 2
 
      call xc_functional_init(xcs%functional(func_x, isp), namespace, x_id, ndim, nel, isp)
      call xc_functional_init(xcs%functional(func_c, isp), namespace, c_id, ndim, nel, isp)
 
      call xc_functional_init(xcs%kernel(func_x, isp), namespace, xk_id, ndim, nel, isp)
      call xc_functional_init(xcs%kernel(func_c, isp), namespace, ck_id, ndim, nel, isp)
 
    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)
 
    ! Take care of hybrid functionals (they appear in the correlation functional)
    xcs%cam_omega = m_zero
    xcs%cam_alpha = m_zero
    xcs%cam_beta = m_zero
 
    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( any(abs(xcs%cam_ext) > m_epsilon) ) call set_hybrid_params(xcs,namespace)
 
        call xc_f03_hyb_cam_coef(xcs%functional(func_c,1)%conf, xcs%cam_omega, &
          xcs%cam_alpha, xcs%cam_beta)
      else
        ! we are doing Hartree-Fock plus possibly a correlation functional
        xcs%cam_omega = m_zero
        xcs%cam_alpha = m_one
        xcs%cam_beta = m_zero
      end if
 
      ! reset certain variables
      xcs%functional(func_x,1)%family = xc_family_oep
      xcs%functional(func_x,1)%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 XCKernelLRCAlpha
      !%Type float
      !%Default 0.0
      !%Section Hamiltonian::XC
      !%Description
      !% Set to a non-zero value to add a long-range correction for solids to the kernel.
      !% This is the <math>\alpha</math> parameter defined in S. Botti <i>et al.</i>, <i>Phys. Rev. B</i>
      !% 69, 155112 (2004). The <math>\Gamma = \Gamma` = 0</math> term <math>-\alpha/q^2</math> is taken
      !% into account by introducing an additional pole to the polarizability (see R. Stubner
      !% <i>et al.</i>, <i>Phys. Rev. B</i> 70, 245119 (2004)). The rest of the terms are included by
      !% multiplying the Hartree term by <math>1 - \alpha / 4 \pi</math>. The use of non-zero
      !% <math>\alpha</math> in combination with <tt>HamiltonianVariation</tt> = <tt>V_ext_only</tt>
      !% corresponds to account of only the <math>\Gamma = \Gamma` = 0</math> term.
      !% Applicable only to isotropic systems. (Experimental)
      !%End
 
      call parse_variable(namespace, 'XCKernelLRCAlpha', m_zero, xcs%kernel_lrc_alpha)
      if (abs(xcs%kernel_lrc_alpha) > m_epsilon) then
        call messages_experimental("Long-range correction to kernel", namespace=namespace)
      end if
 
      !%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 the <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
      xcs%cam_ext = m_zero
 
      if(parse_block(namespace, 'HybridCamParameters', blk) == 0) then
        do isp = 1, size(xcs%cam_ext)
          call parse_block_float(blk, 0, isp - 1, xcs%cam_ext(isp))
        end do
        call parse_block_end(blk)
      end if
 
      if( any(abs(xcs%cam_ext) > m_epsilon) ) then
        write(message(1), '(1x)')
        write(message(2), '(a,f8.5)') "External cam_alpha = ", xcs%cam_ext(1)
        write(message(3), '(a,f8.5)') "External cam_beta  = ", xcs%cam_ext(2)
        write(message(4), '(a,f8.5)') "External cam_omega = ", xcs%cam_ext(3)
        call messages_info(4, namespace=namespace)
      end if
 
      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
 
    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, optional, 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
 
  ! ---------------------------------------------------------
  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. (bitand(xcs%functional(ixc, 1)%flags, xc_flags_have_exc) /= 0 )) 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
 
  !  ---------------------------------------------------------
  subroutine copy_global_to_local(global, local, n_block, nspin, ip)
    real(real64),  intent(in)  :: global(:,:)
    real(real64),  intent(out) :: local(:,:)
    integer,       intent(in)  :: n_block
    integer,       intent(in)  :: nspin
    integer,       intent(in)  :: ip
 
    integer :: ib, is
 
    push_sub(copy_global_to_local)
 
    do is = 1, nspin
      !$omp parallel do
      do ib = 1, n_block
        local(is, ib) = global(ib + ip - 1, is)
      end do
    end do
 
    pop_sub(copy_global_to_local)
  end subroutine copy_global_to_local
 
  ! ---------------------------------------------------------
  subroutine copy_local_to_global(local, global, n_block, spin_channels, ip)
    real(real64),   intent(in)    :: local(:,:)
    real(real64),   intent(inout) :: global(:,:)
    integer,        intent(in)    :: n_block
    integer,        intent(in)    :: spin_channels
    integer,        intent(in)    :: ip
 
    integer :: ib, is
 
    push_sub(xc_compute_vxc.copy_local_to_global)
 
    do is = 1, spin_channels
      !$omp parallel do
      do ib = 1, n_block
        global(ib + ip - 1, is) = global(ib + ip - 1, is) + local(is, ib)
      end do
    end do
 
    pop_sub(xc_compute_vxc.copy_local_to_global)
  end subroutine copy_local_to_global
 
  subroutine set_hybrid_params(xcs,namespace)
    type(namespace_t), intent(in)  :: namespace
    type(xc_t)  :: xcs
    real(real64) :: parameters(3)
 
    parameters =  xcs%cam_ext
 
    select case(xcs%functional(func_c,1)%id)
    case(xc_hyb_gga_xc_pbeh) ! original pbe0 in libxc
      ! If the value is negative, we use the default value of PBE0
      if(parameters(1) < m_zero) parameters(1) = 0.25_real64
 
      call xc_f03_func_set_ext_params(xcs%functional(func_c,1)%conf, parameters)
      xcs%cam_alpha = parameters(1)
 
      write(message(1), '(a,f6.3,a)') 'Info: Setting PBEH mixing parameter (' , parameters(1) ,').'
      call messages_info(1)
 
    case(xc_hyb_gga_xc_cam_pbeh)
      call xc_f03_func_set_ext_params(xcs%functional(func_c,1)%conf, parameters)
      xcs%cam_alpha = parameters(1)
      xcs%cam_beta = parameters(2)
      xcs%cam_omega = parameters(3)
      ! check parameters
      call xc_f03_hyb_cam_coef(xcs%functional(func_c,1)%conf, xcs%cam_omega, &
        xcs%cam_alpha, xcs%cam_beta)
      write(message(1), '(a,f6.3,a)') 'Info: Setting alpha (CAM_PBE) parameter (' , xcs%cam_alpha ,').'
      write(message(2), '(a,f6.3,a)') 'Info: Setting  beta (CAM_PBE) parameter (' , xcs%cam_beta ,').'
      write(message(3), '(a,f6.3,a)') 'Info: Setting omega (CAM_PBE) parameter (' , xcs%cam_omega ,').'
      call messages_info(3)
 
    case(xc_hyb_lda_xc_cam_lda0) ! CAM LDA
      call xc_f03_func_set_ext_params(xcs%functional(func_c,1)%conf, parameters)
      xcs%cam_alpha = parameters(1)
      xcs%cam_beta = parameters(2)
      xcs%cam_omega = parameters(3)
      ! check parameters
      call xc_f03_hyb_cam_coef(xcs%functional(func_c,1)%conf, xcs%cam_omega, &
        xcs%cam_alpha, xcs%cam_beta)
      write(message(1), '(a,f6.3,a)') 'Info: Setting alpha (CAM_LDA) parameter (' , xcs%cam_alpha ,').'
      write(message(2), '(a,f6.3,a)') 'Info: Setting  beta (CAM_LDA) parameter (' , xcs%cam_beta ,').'
      write(message(3), '(a,f6.3,a)') 'Info: Setting omega (CAM_LDA) parameter (' , xcs%cam_omega ,').'
      call messages_info(3)
 
    case(xc_hyb_lda_xc_lda0)
      ! If the value is negative, we use the default value of LDA0
      if(parameters(1) < m_zero) parameters(1) = 0.25_real64
 
      call xc_f03_func_set_ext_params(xcs%functional(func_c,1)%conf, parameters)
      xcs%cam_alpha = parameters(1)
 
      write(message(1), '(a,f6.3,a)') 'Info: Setting LDA0 mixing parameter (' , parameters(1) ,').'
      call messages_info(1)
 
    case default
      assert(.false.)
    end select
 
  end subroutine set_hybrid_params
 
#include "xc_vxc_inc.F90"
#include "xc_fxc_inc.F90"
#include "xc_kxc_inc.F90"
 
#include "xc_vxc_nc_inc.F90"
 
end module xc_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: