xc_sic.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2022 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 xc_sic_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use lalg_basic_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use multicomm_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 varinfo_oct_m
  use xc_oct_m
  use xc_f03_lib_m
  use xc_oep_oct_m
  use xc_kernel_oct_m
  use xc_vxc_oct_m
 
  implicit none
 
  private
  public ::                     &
    xc_sic_t,                   &
    xc_sic_init,                &
    xc_sic_calc_adsic,          &
    xc_sic_write_info,          &
    xc_sic_add_fxc_adsic,       &
    xc_sic_end
 
  integer, parameter, public :: &
    SIC_NONE   = 1,     &  !< no self-interaction correction
    sic_pz_oep = 2,     &  
    sic_amaldi = 3,     &  
    sic_adsic  = 4         
 
  type xc_sic_t
    private
    integer,        public :: level = sic_none  
    real(real64),   public :: amaldi_factor
    type(xc_oep_t), public :: oep
  end type xc_sic_t
 
contains
 
  ! ---------------------------------------------------------
  !
  subroutine xc_sic_init(sic, namespace, gr, st, mc, space)
    type(xc_sic_t),      intent(out)   :: sic
    type(namespace_t),   intent(in)    :: namespace
    type(grid_t),        intent(inout) :: gr
    type(states_elec_t), intent(in)    :: st
    type(multicomm_t),   intent(in)    :: mc
    class(space_t),      intent(in)    :: space
 
 
    push_sub(xc_sic_init)
 
    !%Variable SICCorrection
    !%Type integer
    !%Default sic_none
    !%Section Hamiltonian::XC
    !%Description
    !% This variable controls which form of self-interaction correction to use. Note that
    !% this correction will be applied to the functional chosen by <tt>XCFunctional</tt>.
    !%Option sic_none 1
    !% No self-interaction correction.
    !%Option sic_pz 2
    !% Perdew-Zunger SIC, handled by the OEP technique.
    !% J. P. Perdew and Alex Zunger, Phys. Rev. B 23, 5048 (1981)
    !% Extension to the spinor case follows Tancogne-Dejean et al., J. Chem. Phys. 159, 224110 (2023)
    !%
    !% Note that the current implement uses canonical orbitals and not minimizing orbitals.
    !% Please check <tt>SCDMforPZSIC</tt> for using SCDM-based Wannier orbitals instead of canonical orbitals.
    !%Option sic_amaldi 3
    !% Amaldi correction term. Not implemeneted for spinors.
    !% E. Fermi and E. Amaldi, Mem. Reale Accad. Italia 6, 119 (1934)
    !%Option sic_adsic 4
    !% Average-density SIC.
    !% C. Legrand <i>et al.</i>, <i>J. Phys. B</i> <b>35</b>, 1115 (2002).
    !% Extension to the spinor case follows Tancogne-Dejean et al., J. Chem. Phys. 159, 224110 (2023)
    !%End
    call parse_variable(namespace, 'SICCorrection', sic_none, sic%level)
    if (.not. varinfo_valid_option('SICCorrection', sic%level)) call messages_input_error(namespace, 'SICCorrection')
 
    ! check whether we should introduce the Amaldi SIC correction
    sic%amaldi_factor = m_one
    if (sic%level == sic_amaldi) then
      sic%amaldi_factor = (st%qtot - m_one)/st%qtot
      if(st%d%ispin == spinors) then
        call messages_not_implemented("Amaldi SIC with non-collinear spins")
      end if
    end if
 
    if(sic%level == sic_pz_oep) then
      call xc_oep_init(sic%oep, namespace, gr, st, mc, space, oep_type = oep_type_sic)
 
      if(st%nik > st%d%spin_channels) then
        call messages_not_implemented("PZ-SIC with k-points")
      end if
    end if
 
    if (allocated(st%rho_core)) then
      call messages_not_implemented('SIC with nonlinear core corrections')
    end if
 
    if (allocated(st%frozen_rho) .and. (sic%level == sic_pz_oep .or. sic%level == sic_amaldi)) then
      call messages_not_implemented('PZ-SIC with frozen orbitals')
    end if
 
    if (space%is_periodic() .and. sic%level /= sic_none) then
      call messages_not_implemented("SIC corrections in periodic systems")
    end if
 
    pop_sub(xc_sic_init)
  end subroutine xc_sic_init
 
  ! ---------------------------------------------------------
  subroutine xc_sic_end(sic)
    type(xc_sic_t),  intent(inout) :: sic
 
    if (sic%level == sic_none) return
 
    push_sub(xc_sic_end)
 
    if(sic%level == sic_pz_oep) call xc_oep_end(sic%oep)
 
    pop_sub(xc_sic_end)
  end subroutine xc_sic_end
 
 
  ! ---------------------------------------------------------
  subroutine xc_sic_write_info(sic, iunit, namespace)
    type(xc_sic_t),              intent(in) :: sic
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    if (sic%level == sic_none) return
 
    push_sub(xc_sic_write_info)
 
    call messages_print_var_option('SICCorrection', sic%level, iunit=iunit, namespace=namespace)
 
    pop_sub(xc_sic_write_info)
  end subroutine xc_sic_write_info
 
  ! ---------------------------------------------------------
  subroutine xc_sic_calc_adsic(sic, namespace, space, gr, st, hm, xc, density, vxc, ex, ec)
    type(xc_sic_t),              intent(in) :: sic
    type(namespace_t),           intent(in) :: namespace
    class(space_t),              intent(in) :: space
    type(grid_t),                intent(in) :: gr
    type(states_elec_t),         intent(in) :: st
    type(hamiltonian_elec_t),    intent(in) :: hm
    type(xc_t),               intent(inout) :: xc
    real(real64), contiguous,           intent(in) :: density(:,:)
    real(real64), contiguous,        intent(inout) :: vxc(:,:)
    real(real64), optional,          intent(inout) :: ex, ec
 
    integer            :: ispin, ist, ik, ip
    real(real64), allocatable :: vxc_sic(:,:),  vh_sic(:), rho(:, :)
    real(real64) :: ex_sic, ec_sic, qsp(2)
    real(real64) :: dtot, dpol, vpol, wpol
    real(real64)  :: nup
    real(real64), parameter :: adsic_dpol_reg = 1.0e-10_real64
 
    push_sub(xc_sic_calc_adsic)
 
    assert(sic%level == sic_adsic)
    assert(present(ex) .eqv. present(ec))
 
    if (st%d%ispin == spinors .and. .not. in_family(hm%xc%family, [xc_family_lda, xc_family_gga])) then
      write(message(1),'(a)') 'ADSIC with non-collinear spin is currently only possible'
      write(message(2),'(a)') 'with LDA and GGA functionals.'
      call messages_fatal(2, namespace=namespace)
    end if
 
    if (xc_is_not_size_consistent(xc, namespace)) then
      call messages_not_implemented('ADSIC with size inconsistent functionals',  namespace=namespace)
    end if
 
    ! We compute here the number of electrons per spin channel
    qsp = m_zero
    if( .not. allocated(st%frozen_rho)) then
      select case (st%d%ispin)
      case (unpolarized, spin_polarized)
        do ist = 1, st%nst
          do ik = 1, st%nik
            ispin = st%d%get_spin_index(ik)
            qsp(ispin) = qsp(ispin) + st%occ(ist, ik) * st%kweights(ik)
          end do
        end do
      end select
    else
      ! In the case of the frozen density, we can only get the charge from the integral
      ! of the total density, including valence and frozen density
      qsp(1:st%d%spin_channels) = dmf_integrate(gr, st%d%spin_channels, density)
    end if
 
    safe_allocate(vxc_sic(1:gr%np, 1:2))
    safe_allocate(vh_sic(1:gr%np))
    safe_allocate(rho(1:gr%np, 1:2))
    ! We first compute the average xc self-interction error and we substract it
    select case (st%d%ispin)
    case (unpolarized, spin_polarized)
      do ispin = 1, st%d%spin_channels
        if (abs(qsp(ispin)) <= m_min_occ) cycle
 
        rho = m_zero
        vxc_sic = m_zero
 
        rho(:, ispin) = density(:, ispin) / qsp(ispin)
        if(present(ex)) then
          ex_sic = m_zero
          ec_sic = m_zero
          ! This needs always to be called for the spin-polarized case
          ! force_host is needed to ensure the correct density is used for the libxc call on GPU (see comment in xc_vxc_inc.F90)
          call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
            rho, spin_polarized, hm%ions%latt%rcell_volume, vxc_sic, ex = ex_sic, ec = ec_sic, force_host=.true.)
          ex = ex - ex_sic * qsp(ispin)
          ec = ec - ec_sic * qsp(ispin)
        else
          ! This needs always to be called for the spin-polarized case
          ! force_host is needed to ensure the correct density is used for the libxc call on GPU (see comment in xc_vxc_inc.F90)
          call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
            rho, spin_polarized, hm%ions%latt%rcell_volume, vxc_sic, force_host=.true.)
        end if
 
        call lalg_axpy(gr%np, -m_one, vxc_sic(:, ispin), vxc(:, ispin))
 
        ! We now substract the averaged Hartree self-interaction error
        ! See Eq. 15 in [Pietezak and Vieira, Theoretical Chemistry Accounts (2021) 140:130]
        vh_sic = m_zero
        call dpoisson_solve(hm%psolver, namespace, vh_sic, rho(:, ispin), all_nodes=.false.)
        call lalg_axpy(gr%np, -m_one, vh_sic, vxc(:, ispin))
 
        ! Compute the corresponding energy contribution
        if(present(ex)) then
          ex = ex - m_half*dmf_dotp(gr, rho(:,ispin), vh_sic) * qsp(ispin)
        end if
 
      end do
 
    case (spinors)
      ! Here we only treat the case of LDA/GGA. We rotate the average density in the local frame
      ! And we then compute the SIC correction from it
      ! This cannot excerce any xc torque, by construction
      assert(in_family(hm%xc%family, [xc_family_lda, xc_family_gga]))
 
      do ispin = 1, 2
        rho = m_zero
        vxc_sic = m_zero
        ! Averaged density in the local frame
        do ip = 1, gr%np
          dtot = density(ip, 1) + density(ip, 2)
          dpol = sqrt((density(ip, 1) - density(ip, 2))**2 + &
            m_four*(density(ip, 3)**2 + density(ip, 4)**2))
          if(ispin == 1) then
            rho(ip, 1) = max(m_half*(dtot + dpol), m_zero)
          else
            rho(ip, 2) = max(m_half*(dtot - dpol), m_zero)
          end if
        end do
        nup = dmf_integrate(gr, rho(:,ispin))
        if (nup <= m_min_occ) cycle
        call lalg_scal(gr%np, m_one/nup, rho(:,ispin))
 
        ! This needs always to be called for the spin-polarized case
        if(present(ex) .and. present(ec)) then
          ex_sic = m_zero
          ec_sic = m_zero
          call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
            rho, spin_polarized, hm%ions%latt%rcell_volume, vxc_sic, ex = ex_sic, ec = ec_sic, force_host=.true.)
          ex = ex - ex_sic * nup
          ec = ec - ec_sic * nup
        else
          call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
            rho, spin_polarized, hm%ions%latt%rcell_volume, vxc_sic, force_host=.true.)
        end if
 
        ! Select only the potential correspond to this spin channel
        if(ispin == 2) then
          vxc_sic(:, 1) = m_zero
        else
          vxc_sic(:, 2) = m_zero
        end if
 
        vh_sic = m_zero
        call dpoisson_solve(hm%psolver, namespace, vh_sic, rho(:, ispin), all_nodes=.false.)
        call lalg_axpy(gr%np, m_one, vh_sic, vxc_sic(:, ispin))
        ! Compute the corresponding energy contribution
        if(present(ex)) then
          ex = ex - m_half*dmf_dotp(gr, rho(:,ispin), vh_sic) * nup
        end if
 
        do ip = 1, gr%np
          dpol = sqrt((density(ip, 1) - density(ip, 2))**2 + &
            m_four*(density(ip, 3)**2 + density(ip, 4)**2))
          ! See lda_process in xc_vxc_inc.F90
          if (dpol > xc_tiny*(density(ip, 1)+density(ip, 2))) then
            ! Unlike for the LDA potential, the two local-frame SIC potentials do not become
            ! equal as |m| -> 0 (they are normalized by different electron numbers), so the
            ! magnetic part of the correction must be damped where the local frame direction
            ! m/|m| is not well defined.
            wpol = dpol**2/(dpol**2 + (adsic_dpol_reg*(density(ip, 1)+density(ip, 2)))**2)
            vpol = wpol*(vxc_sic(ip, 1) - vxc_sic(ip, 2))*(density(ip, 1) - density(ip, 2))/(safe_tol(dpol, xc_tiny))
 
            vxc(ip, 1) = vxc(ip, 1) - m_half*(vxc_sic(ip, 1) + vxc_sic(ip, 2) + vpol)
            vxc(ip, 2) = vxc(ip, 2) - m_half*(vxc_sic(ip, 1) + vxc_sic(ip, 2) - vpol)
            vxc(ip, 3) = vxc(ip, 3) - wpol*(vxc_sic(ip, 1) - vxc_sic(ip, 2))*density(ip, 3)/(safe_tol(dpol, xc_tiny))
            vxc(ip, 4) = vxc(ip, 4) - wpol*(vxc_sic(ip, 1) - vxc_sic(ip, 2))*density(ip, 4)/(safe_tol(dpol, xc_tiny))
          else
            vxc(ip, 1) = vxc(ip, 1) - m_half*(vxc_sic(ip, 1) + vxc_sic(ip, 2))
            vxc(ip, 2) = vxc(ip, 2) - m_half*(vxc_sic(ip, 1) + vxc_sic(ip, 2))
          end if
        end do
      end do
 
 
    end select
 
    safe_deallocate_a(vxc_sic)
    safe_deallocate_a(vh_sic)
    safe_deallocate_a(rho)
 
    pop_sub(xc_sic_calc_adsic)
  end subroutine xc_sic_calc_adsic
 
  ! ---------------------------------------------------------
  subroutine xc_sic_add_fxc_adsic(namespace, xc, st, gr, rho, fxc, fxc_grad, fxc_grad_spin, triplet)
    type(namespace_t),        intent(in)    :: namespace
    type(xc_t),               intent(in)    :: xc
    type(states_elec_t),      intent(in)    :: st
    type(grid_t),             intent(in)    :: gr
    real(real64),             intent(in)    :: rho(:,:)
    real(real64), contiguous, intent(inout) :: fxc(:,:,:)
    real(real64), contiguous, intent(inout) :: fxc_grad(:,:,:,:,:)
    real(real64), contiguous, intent(inout) :: fxc_grad_spin(:,:,:,:)
    logical,                  intent(in)    :: triplet
 
    real(real64), allocatable :: rho_averaged(:, :)
    real(real64), allocatable :: fxc_sic(:,:,:)
    real(real64), allocatable :: fxc_grad_sic(:,:,:,:,:)
    real(real64), allocatable :: fxc_grad_spin_sic(:,:,:,:)
    real(real64) :: qtot
    integer :: ispin
 
    push_sub(xc_sic_add_fxc_adsic)
 
    !Check spin and triplets
    assert(st%d%ispin /= spinors)
    assert(.not. allocated(st%frozen_rho))
 
    if ( bitand(xc%kernel_family, xc_family_lda) == 0 .and. bitand(xc%kernel_family, xc_family_gga) == 0 ) then
      message(1) = "fxc calculation with ADSIC not implemented beyond LDAs and GGAs."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (xc_is_not_size_consistent(xc, namespace)) then
      call messages_not_implemented('ADSIC with size inconsistent functionals',  namespace=namespace)
    end if
 
    if (st%d%ispin == spinors) then
      call messages_not_implemented('ADSIC fxc with non-collinear spin')
    end if
 
 
    ! This needs always to be called for the spin-polarized case
    ! SAFE_ALLOCATE(rho_averaged(1:gr%np, 1:2))
    safe_allocate(rho_averaged(1:gr%np_part, 1:2))
    safe_allocate(fxc_sic(1:gr%np, 1:2, 1:2))
    safe_allocate(fxc_grad_sic(1:gr%np, 1:gr%der%dim, 1:gr%der%dim, 1:2, 1:2))
    safe_allocate(fxc_grad_spin_sic(1:gr%np, 1:gr%der%dim, 1:2, 1:2))
 
    do ispin = 1, st%d%nspin
      rho_averaged = m_zero
      qtot = dmf_integrate(gr, rho(:, ispin))
      if (abs(qtot) <= m_min_occ) cycle
 
      call lalg_copy(gr%np, rho(:, ispin), rho_averaged(:, ispin))
      call lalg_scal(gr%np, m_one/qtot, rho_averaged(:, ispin))
 
      if(triplet) then
        call lalg_scal(gr%np, m_half, rho_averaged(:, 1))
        call lalg_copy(gr%np, rho_averaged(:, 1), rho_averaged(:, 2))
      endif
 
      fxc_sic = m_zero
      fxc_grad_sic = m_zero
      fxc_grad_spin_sic = m_zero
      call xc_get_fxc(xc, gr, namespace, rho_averaged, spin_polarized, fxc_sic, &
        fxc_grad_sic, fxc_grad_spin_sic)
 
      ! For Casida triplet, fxc = 1/2 (fxc_up[n_up] - fxc_up [n_down])
      if (triplet) then
        call lalg_axpy(gr%np, -m_half/qtot, fxc_sic(:, 1, 1), fxc(:, 1, 1))
        call lalg_axpy(gr%np, +m_half/qtot, fxc_sic(:, 1, 2), fxc(:, 1, 1))
        call lalg_axpy(gr%np, gr%der%dim, gr%der%dim, -m_half/qtot, fxc_grad_sic(:, :, :, 1, 1), fxc_grad(:, :, :, 1, 1))
        call lalg_axpy(gr%np, gr%der%dim, gr%der%dim, +m_half/qtot, fxc_grad_sic(:, :, :, 1, 2), fxc_grad(:, :, :, 1, 1))
      else
        call lalg_axpy(gr%np, -m_one/qtot, fxc_sic(:, ispin, ispin), fxc(:, ispin, ispin))
        call lalg_axpy(gr%np, gr%der%dim, gr%der%dim, -m_one/qtot, fxc_grad_sic(:, :, :, ispin, ispin), &
          fxc_grad(:, :, :, ispin, ispin))
        if (st%d%nspin > 1) then
          call lalg_axpy(gr%np, gr%der%dim, -m_one/qtot, fxc_grad_spin_sic(:, :, ispin, ispin), &
            fxc_grad_spin(:, :, ispin, ispin))
        end if
      end if
 
    end do
 
    safe_deallocate_a(rho_averaged)
    safe_deallocate_a(fxc_sic)
    safe_deallocate_a(fxc_grad_sic)
    safe_deallocate_a(fxc_grad_spin_sic)
 
    pop_sub(xc_sic_add_fxc_adsic)
  end subroutine xc_sic_add_fxc_adsic
 
end module xc_sic_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: