xc_interaction.F90

Go to the documentation of this file.

!! Copyright (C) 2021 Nicolas 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 xc_interaction_oct_m
  use comm_oct_m
  use debug_oct_m
  use distributed_oct_m
  use electron_space_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_partner_oct_m
  use, intrinsic :: iso_fortran_env
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use density_interaction_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use space_oct_m
  use states_elec_dim_oct_m
  use xc_f03_lib_m
  use xc_functional_oct_m
 
  implicit none
 
  private
  public ::                &
    xc_interaction_t,      &
    xc_interaction_compute,&
    calc_tb09_c,           &
    calc_mvorb_alpha
 
  type, extends(density_interaction_t) :: xc_interaction_t
    private
 
  contains
    procedure :: init => xc_interaction_init
    procedure :: calculate => xc_interaction_calculate
    procedure :: calculate_energy => xc_interaction_calculate_energy
    procedure :: end => xc_interaction_end
    final :: xc_interaction_finalize
  end type xc_interaction_t
 
  interface xc_interaction_t
    module procedure xc_interaction_constructor
  end interface xc_interaction_t
 
  integer, public, parameter :: &
    FUNC_X = 1,         &
    func_c = 2
 
 
contains
 
  ! ---------------------------------------------------------
  function xc_interaction_constructor(partner) result(this)
    class(interaction_partner_t), target, intent(inout) :: partner
    class(xc_interaction_t),               pointer       :: this
 
    push_sub(xc_interaction_constructor)
 
    allocate(this)
 
    this%label = "exchange-correlation interaction"
 
    this%partner => partner
 
    pop_sub(xc_interaction_constructor)
  end function xc_interaction_constructor
 
  ! ---------------------------------------------------------
  subroutine xc_interaction_init(this)
    class(xc_interaction_t),    intent(inout) :: this
 
    push_sub(xc_interaction_init)
 
    pop_sub(xc_interaction_init)
  end subroutine xc_interaction_init
 
  ! ---------------------------------------------------------
  subroutine xc_interaction_calculate(this)
    class(xc_interaction_t),             intent(inout) :: this
 
    push_sub(xc_interaction_calculate)
 
    this%density = m_zero
 
    pop_sub(xc_interaction_calculate)
  end subroutine xc_interaction_calculate
 
  ! ---------------------------------------------------------
  subroutine xc_interaction_calculate_energy(this)
    class(xc_interaction_t),             intent(inout) :: this
 
    push_sub(xc_interaction_calculate_energy)
 
    this%energy = m_zero
 
    pop_sub(xc_interaction_calculate_energy)
  end subroutine xc_interaction_calculate_energy
 
 
  ! ---------------------------------------------------------
  subroutine xc_interaction_compute(this)
    class(xc_interaction_t),             intent(inout) :: this
    push_sub(xc_interaction_compute)
 
    this%density = m_zero
 
    pop_sub(xc_interaction_compute)
  end subroutine xc_interaction_compute
 
  ! ---------------------------------------------------------
  subroutine xc_interaction_end(this)
    class(xc_interaction_t), intent(inout) :: this
 
    push_sub(xc_interaction_end)
 
    safe_deallocate_a(this%density)
 
    call interaction_end(this)
 
    pop_sub(xc_interaction_end)
  end subroutine xc_interaction_end
 
 
  ! ---------------------------------------------------------
  subroutine xc_interaction_finalize(this)
    type(xc_interaction_t), intent(inout) :: this
 
    push_sub(xc_interaction_finalize)
 
    call xc_interaction_end(this)
 
    pop_sub(xc_interaction_finalize)
  end subroutine xc_interaction_finalize
 
  ! -----------------------------------------------------
  subroutine calc_tb09_c(mesh, space, functl, dens, gdens, ispin, rcell_volume)
    class(mesh_t),         intent(in) :: mesh
    class(space_t),        intent(in) :: space
    type(xc_functional_t), intent(inout) :: functl(:)
    real(real64),          intent(in) :: dens(:,:)
    real(real64),          intent(in) :: gdens(:,:,:)
    integer,               intent(in) :: ispin
    real(real64),          intent(in) :: rcell_volume
 
    real(real64), allocatable :: gnon(:)
    real(real64) :: gn(space%dim), n, parameters(1)
    integer :: ii
 
    push_sub(calc_tb09_c)
 
    safe_allocate(gnon(1:mesh%np))
 
    do ii = 1, mesh%np
      if (ispin == unpolarized) then
        n = dens(ii, 1)
        gn(1:space%dim) = gdens(ii, 1:space%dim, 1)
      else
        n = dens(ii, 1) + dens(ii, 2)
        gn(1:space%dim) = gdens(ii, 1:space%dim, 1) + gdens(ii, 1:space%dim, 2)
      end if
 
      if (n <= 1e-7_real64) then
        gnon(ii) = m_zero
      else
        gnon(ii) = sqrt(sum((gn(1:space%dim)/n)**2))
      end if
    end do
 
    parameters(1) =  -0.012_real64 + 1.023_real64*sqrt(dmf_integrate(mesh, gnon)/rcell_volume)
 
    call xc_f03_func_set_ext_params(functl(1)%conf, parameters)
 
    safe_deallocate_a(gnon)
 
    pop_sub(calc_tb09_c)
  end subroutine calc_tb09_c
 
  ! ---------------------------------------------------------
  subroutine calc_mvorb_alpha(mesh, namespace, space, functl, dens, gdens, ispin, rcell_volume, &
    cam_alpha, cam_beta, cam_omega)
    class(mesh_t),         intent(in) :: mesh
    type(namespace_t),     intent(in) :: namespace
    class(space_t),        intent(in) :: space
    type(xc_functional_t), intent(inout) :: functl(:)
    real(real64),          intent(in) :: dens(:,:)
    real(real64),          intent(in) :: gdens(:,:,:)
    integer,               intent(in) :: ispin
    real(real64),          intent(in) :: rcell_volume
    real(real64),          intent(inout) :: cam_alpha, cam_beta, cam_omega
 
    real(real64), allocatable :: gnon(:)
    real(real64) :: tb09_c, alpha
    real(real64) :: gn(space%dim), n
    integer :: ii
    real(real64) :: parameters(3)
 
    push_sub(calc_mvorb_alpha)
 
    safe_allocate(gnon(1:mesh%np))
 
    do ii = 1, mesh%np
      if (ispin == unpolarized) then
        n = dens(ii, 1)
        gn(1:space%dim) = gdens(ii, 1:space%dim, 1)
      else
        n = dens(ii, 1) + dens(ii, 2)
        gn(1:space%dim) = gdens(ii, 1:space%dim, 1) + gdens(ii, 1:space%dim, 2)
      end if
 
      if (n <= 1e-7_real64) then
        gnon(ii) = m_zero
      else
        gnon(ii) = sqrt(sum((gn(1:space%dim)/n)**2))
        gnon(ii) = sqrt(gnon(ii))
      end if
    end do
 
    tb09_c =  dmf_integrate(mesh, gnon)/rcell_volume
 
    safe_deallocate_a(gnon)
 
    select case (functl(func_c)%id)
    case (xc_hyb_gga_xc_mvorb_hse06)
      alpha = 0.121983_real64+0.130711_real64*tb09_c**4
 
      if (alpha > 1) then
        write(message(1), '(a,f6.3,a)') 'MVORB mixing parameter bigger than one (' , alpha ,').'
        call messages_warning(1, namespace=namespace)
        alpha = 0.25_real64
      end if
 
 
      parameters(1) = alpha
      parameters(2) = cam_omega
      parameters(3) = cam_omega
      call xc_f03_func_set_ext_params(functl(func_c)%conf, parameters)
      !The name is confusing. Here alpha is the beta of hybrids in functionals,
      !but is called alpha in the original paper.
      cam_beta = alpha
 
    case (xc_hyb_gga_xc_mvorb_pbeh)
      alpha = -1.00778_real64+1.10507_real64*tb09_c
 
      if (alpha > 1) then
        write(message(1), '(a,f6.3,a)') 'MVORB mixing parameter bigger than one (' , alpha ,').'
        call messages_warning(1, namespace=namespace)
        alpha = 0.25_real64
      end if
      if (alpha < 0) then
        write(message(1), '(a,f6.3,a)') 'MVORB mixing parameter smaller than zero (' , alpha ,').'
        call messages_warning(1, namespace=namespace)
        alpha = 0.25_real64
      end if
 
      parameters(1) = alpha
      call xc_f03_func_set_ext_params(functl(func_c)%conf, parameters)
      cam_alpha = alpha
    case default
      call messages_not_implemented("MVORB density-based mixing for functionals other than PBE0 and HSE06", namespace=namespace)
    end select
 
    pop_sub(calc_mvorb_alpha)
  end subroutine calc_mvorb_alpha
 
 
end module xc_interaction_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: