xc_oep.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2012-2013 M. Gruning, P. Melo, M. Oliveira
!! 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_oep_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use exchange_operator_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_elec_base_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lalg_basic_oct_m
  use lalg_adv_oct_m
  use lattice_vectors_oct_m
  use linear_response_oct_m
  use linear_solver_oct_m
  use math_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use scdm_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_calc_oct_m
  use states_elec_dim_oct_m
  use scdm_oct_m
  use scf_tol_oct_m
  use varinfo_oct_m
  use xc_oct_m
  use xc_f03_lib_m
  use xc_functional_oct_m
 
  implicit none
 
  private
  public ::                     &
    xc_oep_t,                   &
    xc_oep_init,                &
    xc_oep_end,                 &
    xc_oep_write_info,          &
    dxc_oep_calc,               &
    zxc_oep_calc,               &
    xc_oep_analyzeeigen,        &
    xc_oep_spinfactor,          &
    dxc_oep_mix,                &
    zxc_oep_mix
 
  integer, public, parameter ::  &
    OEP_LEVEL_NONE   = 1,           &
    oep_level_kli    = 3,           &
    oep_level_full   = 5
 
  integer, public, parameter ::  &
    OEP_MIXING_SCHEME_CONST = 1, &
    oep_mixing_scheme_bb    = 2, &
    oep_mixing_scheme_dens  = 3
 
  integer, public, parameter ::  &
    OEP_TYPE_EXX     = 1,          &
    oep_type_mgga    = 2,          &
    oep_type_sic     = 3,          &
    oep_type_photons = 4             ! one-photon OEP
 
  type xc_oep_t
    integer,              public :: level
    real(real64)                 :: mixing
    integer                      :: mixing_scheme
    type(lr_t)                   :: lr
    type(linear_solver_t)        :: solver
    type(scf_tol_t)              :: scftol
    integer                      :: eigen_n
    integer, allocatable         :: eigen_type(:), eigen_index(:)
    real(real64)                 :: socc, sfact
    real(real64),   allocatable, public :: vxc(:,:), uxc_bar(:,:,:)
    real(real64),   allocatable         :: dlxc(:, :, :, :)
    complex(real64),   allocatable         :: zlxc(:, :, :, :)
    real(real64),         public :: norm2ss
    real(real64),   allocatable         :: vxc_old(:,:), ss_old(:,:)
    integer                      :: noccst
 
    type(scdm_t)                 :: scdm
    logical                      :: scdm_for_pzsic
 
    integer, public              :: type = -1
  end type xc_oep_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine xc_oep_init(oep, namespace, gr, st, mc, space, oep_type)
    type(xc_oep_t),      intent(inout) :: oep
    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
    integer,             intent(in)    :: oep_type
 
    push_sub(xc_oep_init)
 
    !%Variable OEPLevel
    !%Type integer
    !%Default oep_kli
    !%Section Hamiltonian::XC
    !%Description
    !% At what level shall <tt>Octopus</tt> handle the optimized effective potential (OEP) equation.
    !%Option oep_none 1
    !% Do not solve OEP equation.
    !%Option oep_kli 3
    !% Krieger-Li-Iafrate (KLI) approximation.
    !% Ref: JB Krieger, Y Li, GJ Iafrate, <i>Phys. Lett. A</i> <b>146</b>, 256 (1990).
    !%Option oep_full 5
    !% (Experimental) Full solution of OEP equation using the Sternheimer approach.
    !% The linear solver will be controlled by the variables in section <tt>Linear Response::Solver</tt>,
    !% and the iterations for OEP by <tt>Linear Response::SCF in LR calculations</tt> and variable
    !% <tt>OEPMixing</tt>. Note that default for <tt>LRMaximumIter</tt> is set to 10.
    !% Ref: S. Kuemmel and J. Perdew, <i>Phys. Rev. Lett.</i> <b>90</b>, 043004 (2003).
    !%End
    call messages_obsolete_variable(namespace, 'OEP_Level', 'OEPLevel')
    call parse_variable(namespace, 'OEPLevel', oep_level_kli, oep%level)
    if (.not. varinfo_valid_option('OEPLevel', oep%level)) call messages_input_error(namespace, 'OEPLevel')
 
    if (oep%level == oep_level_none) then
      pop_sub(xc_oep_init)
      return
    end if
 
    oep%type = oep_type
 
    if (oep%level == oep_level_full) then
 
      if (st%d%nspin == spinors) then
        call messages_not_implemented("Full OEP with spinors", namespace=namespace)
      end if
 
      call messages_experimental("Full OEP")
 
      !%Variable OEPMixing
      !%Type float
      !%Default 1.0
      !%Section Hamiltonian::XC
      !%Description
      !% The linear mixing factor used to solve the Sternheimer
      !% equation in the full OEP procedure.
      !%End
      call messages_obsolete_variable(namespace, 'OEP_Mixing', 'OEPMixing')
      call parse_variable(namespace, 'OEPMixing', m_one, oep%mixing)
 
      !%Variable OEPMixingScheme
      !%Type integer
      !%Default oep_mixing_scheme_const
      !%Section Hamiltonian::XC
      !%Description
      !% Different Mixing Schemes are possible
      !%Option oep_mixing_scheme_const 1
      !% Use a constant
      !% Reference: S. Kuemmel and J. Perdew, <i>Phys. Rev. Lett.</i> <b>90</b>, 4, 043004 (2003)
      !%Option oep_mixing_scheme_bb 2
      !% Use the Barzilai-Borwein (BB) Method
      !% Reference: T. W. Hollins, S. J. Clark, K. Refson, and N. I. Gidopoulos,
      !% <i>Phys. Rev. B</i> <b>85</b>, 235126 (2012)
      !%Option oep_mixing_scheme_dens 3
      !% Use the inverse of the electron density
      !% Reference: S. Kuemmel and J. Perdew, <i>Phys. Rev. B</i> <b>68</b>, 035103 (2003)
      !%End
      call parse_variable(namespace, 'OEPMixingScheme', oep_mixing_scheme_const, oep%mixing_scheme)
 
      if (oep%mixing_scheme == oep_mixing_scheme_bb) then
        safe_allocate(oep%vxc_old(1:gr%np,st%d%ispin))
        safe_allocate(oep%ss_old(1:gr%np,st%d%ispin))
        oep%vxc_old = m_zero
        oep%ss_old = m_zero
      end if
 
      oep%norm2ss = m_zero
    end if
 
    ! obtain the spin factors
    call xc_oep_spinfactor(oep, st%d%nspin)
 
    ! This variable will keep vxc across iterations
    if ((st%d%ispin == spinors) .or. oep%level == oep_level_full) then
      safe_allocate(oep%vxc(1:gr%np,st%d%nspin))
    else
      safe_allocate(oep%vxc(1:gr%np,1:min(st%d%nspin, 2)))
    end if
    oep%vxc = m_zero
 
    ! when performing full OEP, we need to solve a linear equation
    if (oep%level == oep_level_full) then
      call scf_tol_init(oep%scftol, namespace, st%qtot, def_maximumiter=10)
      call linear_solver_init(oep%solver, namespace, gr, states_are_real(st), mc, space)
      call lr_init(oep%lr)
    end if
 
    if (oep%type == oep_type_sic) then
      !%Variable SCDMforPZSIC
      !%Type logical
      !%Default .false.
      !%Section Hamiltonian::XC
      !%Description
      !% If set to .true., the code will use the SCDM method to Wannierize the orbitals before
      !% computing the orbital-dependent SIC correction
      !%End
      call parse_variable(namespace, 'SCDMforPZSIC', .false., oep%scdm_for_pzsic)
      if (oep%scdm_for_pzsic) call scdm_init(oep%scdm, namespace)
    end if
 
    ! the linear equation has to be more converged if we are to attain the required precision
    !oep%lr%conv_abs_dens = oep%lr%conv_abs_dens / (oep%mixing)
 
    if (st%d%nspin == spinors) then
      call messages_experimental("OEP with spinors")
    end if
 
    if (st%d%kpt%parallel .and. oep%level == oep_level_full) then
      call messages_experimental("OEP Full parallel in spin/k-points", namespace=namespace)
    end if
 
    if (st%d%kpt%parallel .and. oep%type == oep_type_sic) then
      call messages_not_implemented("PZ-SIC with k-point parallelization")
    end if
 
    pop_sub(xc_oep_init)
  end subroutine xc_oep_init
 
 
  ! ---------------------------------------------------------
  subroutine xc_oep_end(oep)
    class(xc_oep_t), intent(inout) :: oep
 
    push_sub(xc_oep_end)
 
    if (oep%level /= oep_level_none) then
      safe_deallocate_a(oep%vxc)
      if (oep%level == oep_level_full .or. oep%type == oep_type_photons) then
        call lr_dealloc(oep%lr)
        call linear_solver_end(oep%solver)
      end if
      if (oep%level == oep_level_full .and. oep%mixing_scheme == oep_mixing_scheme_bb) then
        safe_deallocate_a(oep%vxc_old)
        safe_deallocate_a(oep%ss_old)
      end if
    end if
 
    pop_sub(xc_oep_end)
  end subroutine xc_oep_end
 
 
  ! ---------------------------------------------------------
  subroutine xc_oep_write_info(oep, iunit, namespace)
    type(xc_oep_t),              intent(in) :: oep
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    if (oep%level == oep_level_none) return
 
    push_sub(xc_oep_write_info)
 
    call messages_print_var_option('OEPLevel', oep%level, iunit=iunit, namespace=namespace)
 
    pop_sub(xc_oep_write_info)
  end subroutine xc_oep_write_info
 
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine xc_oep_spinfactor(oep, nspin)
    class(xc_oep_t), intent(inout) :: oep
    integer,         intent(in)    :: nspin
 
    push_sub(xc_oep_spinfactor)
 
    select case (nspin)
    case (1) ! we need to correct or the spin occupancies
      oep%socc  = m_half
      oep%sfact = m_two
    case (2, 4)
      oep%socc  = m_one
      oep%sfact = m_one
    case default ! cannot handle any other case
      assert(.false.)
    end select
 
    pop_sub(xc_oep_spinfactor)
  end subroutine xc_oep_spinfactor
 
 
  ! ---------------------------------------------------------
  subroutine xc_oep_analyzeeigen(oep, st, is)
    class(xc_oep_t),      intent(inout) :: oep
    type(states_elec_t),  intent(in)    :: st
    integer,              intent(in)    :: is
 
    integer  :: ist
    real(real64) :: max_eigen
 
    push_sub(xc_oep_analyzeeigen)
 
    ! The arrays are reuse for different spins, we reset them
    oep%eigen_type = -1
    oep%eigen_index = -1
    oep%noccst = 0
    oep%eigen_n = 0
 
    ! find out the top occupied state, to correct for the asymptotics
    ! of the potential
    max_eigen = -m_huge
    do ist = 1, st%nst
      if ((st%occ(ist, is) > m_min_occ).and.(st%eigenval(ist, is) > max_eigen)) then
        max_eigen = st%eigenval(ist, is)
      end if
    end do
 
    do ist = 1, st%nst
      ! criterion for degeneracy
      if (abs(st%eigenval(ist, is)-max_eigen) <= 1e-3_real64) then
        oep%eigen_type(ist) = 2
      else
        if (st%occ(ist, is) > m_min_occ) then
          oep%eigen_type(ist) = 1
          oep%eigen_index(oep%eigen_n+1) = ist
          oep%eigen_n = oep%eigen_n + 1
        else
          oep%eigen_type(ist) = 0
        end if
      end if
    end do
 
    ! find how many states are occupied.
    do ist = 1, st%nst
      if (st%occ(ist, is) > m_min_occ) oep%noccst = ist
    end do
 
    ! We check that all states have been assigned
    assert(all(oep%eigen_type >= 0))
 
    pop_sub(xc_oep_analyzeeigen)
  end subroutine xc_oep_analyzeeigen
 
 
#include "xc_kli_pauli_inc.F90"
#include "xc_oep_sic_pauli_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "xc_kli_inc.F90"
#include "xc_oep_sic_inc.F90"
#include "xc_oep_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "xc_kli_inc.F90"
#include "xc_oep_sic_inc.F90"
#include "xc_oep_inc.F90"
 
end module xc_oep_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: