xc_oep_photon.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
!!
!! 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_photon_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 linear_response_oct_m
  use linear_solver_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use multigrid_oct_m
  use namespace_oct_m
  use parser_oct_m
  use photon_mode_oct_m
  use poisson_oct_m
  use profiling_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 scf_tol_oct_m
  use varinfo_oct_m
  use xc_oct_m
  use xc_oep_oct_m
  use xc_f03_lib_m
  use xc_functional_oct_m
 
  implicit none
 
  private
  public ::                     &
    xc_oep_photon_t,            &
    xc_oep_photon_init,         &
    xc_oep_photon_end,          &
    dxc_oep_photon_calc,        &
    zxc_oep_photon_calc
 
  type, extends(xc_oep_t) :: xc_oep_photon_t
    private
    logical                      :: coc_translation
    type(photon_mode_t),  public :: pt
    type(lr_t)                   :: photon_lr
    logical,              public :: rm_ee_interaction
  end type xc_oep_photon_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine xc_oep_photon_init(oep, namespace, family, gr, st, mc, space)
    type(xc_oep_photon_t),  intent(inout) :: oep
    type(namespace_t),      intent(in)    :: namespace
    integer,                intent(in)    :: family
    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_oep_photon_init)
 
    if(bitand(family, xc_family_oep) == 0) then
      oep%level = oep_level_none
      pop_sub(xc_oep_photon_init)
      return
    end if
 
    oep%type = oep_type_photons
 
    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
      call messages_experimental("EnablePhotons = yes")
      call photon_mode_init(oep%pt, namespace, space%dim)
      call photon_mode_set_n_electrons(oep%pt, st%qtot)
      if (oep%pt%nmodes > 1) then
        call messages_not_implemented('Photon OEP for more than one photon mode.')
      end if
      if (oep%level == oep_level_full .and. st%d%nspin /= unpolarized) then
        call messages_not_implemented('Spin-polarized calculations with photon OEP.')
      end if
 
      if (states_are_complex(st)) then
        call messages_not_implemented('Photon OEP with complex wavefunctions.')
      end if
 
      if (st%nik > st%d%ispin .and. oep%level == oep_level_full) then
        call messages_not_implemented("Full OEP for periodic systems", namespace=namespace)
      end if
      if (st%nik > st%d%ispin .and. st%d%ispin==spinors) then
        call messages_not_implemented("OEP for periodic systems with spinors", namespace=namespace)
      end if
 
      safe_allocate(oep%pt%correlator(1:gr%np, 1:oep%pt%nmodes))
      oep%pt%correlator = m_zero
 
      call photon_mode_compute_dipoles(oep%pt, gr)
    end if
 
    if(oep%level == oep_level_full) then
 
      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 1.0
      !%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==3) .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
 
    !%Variable KLIPhotonCOC
    !%Type logical
    !%Default .false.
    !%Section Hamiltonian::XC
    !%Description
    !% Activate the center of charge translation of the electric dipole operator which should avoid the dependence of the photon KLI on an permanent dipole.
    !%End
 
    ! when performing full OEP, we need to solve a linear equation
    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)
    call lr_init(oep%photon_lr)
    call parse_variable(namespace, 'KLIPhotonCOC', .false., oep%coc_translation)
 
    !%Variable OEPRemoveElectron
    !%Type logical
    !%Default .false.
    !%Section Hamiltonian::XC
    !%Description
    !% Remove electron-electron interaction in OEP-Photon calculations
    !%End
 
    call parse_variable(namespace, 'OEPRemoveElectron', .false., oep%rm_ee_interaction)
 
    ! 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%kpt%parallel) then
      call messages_not_implemented("OEP parallel in spin/k-points", namespace=namespace)
    end if
 
    pop_sub(xc_oep_photon_init)
  end subroutine xc_oep_photon_init
 
 
  ! ---------------------------------------------------------
  subroutine xc_oep_photon_end(oep)
    type(xc_oep_photon_t), intent(inout) :: oep
 
    push_sub(xc_oep_photon_end)
 
    call xc_oep_end(oep)
 
    if(oep%level /= oep_level_none) then
      call lr_dealloc(oep%photon_lr)
      call photon_mode_end(oep%pt)
    end if
 
    pop_sub(xc_oep_photon_end)
  end subroutine xc_oep_photon_end
 
#include "xc_oep_qed_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "xc_kli_photon_inc.F90"
#include "xc_oep_photon_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "xc_kli_photon_inc.F90"
#include "xc_oep_photon_inc.F90"
 
end module xc_oep_photon_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: