v_ks_photon.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 v_ks_photon_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 ions_oct_m
  use kpoints_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use messages_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use photon_mode_oct_m
  use photon_mode_mf_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use xc_oct_m
  use xc_f03_lib_m
  use xc_functional_oct_m
  use xc_interaction_oct_m
  use xc_oep_oct_m
  use xc_oep_photon_oct_m
  use xc_photons_oct_m
 
  implicit none
 
  private
  public :: v_ks_photon_t
 
  type v_ks_photon_t
    private
    integer                      :: xc_photon = 0
    type(xc_photons_t)           :: xc_photons
    type(xc_oep_photon_t)        :: oep_photon
    logical                      :: has_photons = .false.
    logical                      :: include_hartree = .true.
    type(photon_mode_t), pointer :: pt => null()
    type(mf_t)                   :: pt_mx
  contains
    procedure :: init             => v_ks_photon_init
    procedure :: init_xc          => v_ks_photon_init_xc
    procedure :: setup            => v_ks_photon_setup
    procedure :: end              => v_ks_photon_end
    procedure :: mf_calc          => v_ks_photon_mf_calc
    procedure :: add_px           => v_ks_photon_add_px
    procedure :: oep_calc         => v_ks_photon_oep_calc
    procedure :: add_mf_potential => v_ks_photon_add_mf_potential
    procedure :: add_mean_field   => v_ks_photon_add_mean_field
    procedure :: restart_dump     => v_ks_photon_restart_dump
    procedure :: restart_load     => v_ks_photon_restart_load
    procedure :: write_info       => v_ks_photon_write_info
    procedure :: active           => v_ks_photon_active
    procedure :: set_active       => v_ks_photon_set_active
    procedure :: functional       => v_ks_photon_functional
    procedure :: includes_hartree => v_ks_photon_include_hartree
    procedure :: nmodes           => v_ks_photon_nmodes
    procedure :: oep_full         => v_ks_photon_oep_full
    procedure :: correlator       => v_ks_photon_correlator
    procedure :: get_pt_mx        => v_ks_photon_get_pt_mx
    procedure :: get_xc_photons   => v_ks_photon_get_xc_photons
  end type v_ks_photon_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_init(this, namespace)
    class(v_ks_photon_t), intent(inout) :: this
    type(namespace_t),   intent(in)    :: namespace
 
    push_sub(v_ks_photon_init)
 
    this%oep_photon%level = oep_level_none
 
    !%Variable XCPhotonFunctional
    !%Type integer
    !%Default 0
    !%Section Hamiltonian::XC
    !%Description
    !% Defines the exchange and correlation functionals to be used for the QEDFT
    !% description of the electron-photon system.
    !%Option none      0
    !% No functional is used
    !%Option photon_x_lda      10
    !% Exchange-only local density approcimation
    !%Option photon_xc_lda     11
    !% Exchange-correlation local density approcimation
    !%Option photon_x_wfn      20
    !% Exchange-only based on wave functions
    !%Option photon_xc_wfn     21
    !% Exchange-correlation based on wave functions
    !%End
    call parse_variable(namespace, 'XCPhotonFunctional', option__xcphotonfunctional__none, this%xc_photon)
 
    !%Variable XCPhotonIncludeHartree
    !%Type logical
    !%Default yes
    !%Section Hamiltonian::XC
    !%Description
    !% Use the Hartree potential and energy in calculations
    !%End
    call parse_variable(namespace, 'XCPhotonIncludeHartree', .true., this%include_hartree)
 
    if (.not. this%include_hartree) then
      call messages_write('turn off hartree potential and energy')
      call messages_warning(namespace=namespace)
    end if
 
    pop_sub(v_ks_photon_init)
  end subroutine v_ks_photon_init
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_init_xc(this, namespace, space, gr, st)
    class(v_ks_photon_t), intent(inout) :: this
    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
 
    push_sub(v_ks_photon_init_xc)
 
    if (this%xc_photon /= 0) then
      ! initialize the photon free variables
      call this%xc_photons%init(namespace, this%xc_photon, space, gr, st)
      this%oep_photon%level = oep_level_none
    end if
 
    pop_sub(v_ks_photon_init_xc)
  end subroutine v_ks_photon_init_xc
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_setup(this, namespace, modes, gr, st, ions, mc, space, xc, xc_family, oep_level)
    class(v_ks_photon_t),         intent(inout) :: this
    type(namespace_t),           intent(in)    :: namespace
    type(photon_mode_t), target, intent(inout) :: modes
    type(grid_t),                intent(inout) :: gr
    type(states_elec_t),         intent(in)    :: st
    type(ions_t),                intent(in)    :: ions
    type(multicomm_t),           intent(in)    :: mc
    class(space_t),              intent(in)    :: space
    type(xc_t),                  intent(in)    :: xc
    integer,                     intent(in)    :: xc_family
    integer,                     intent(in)    :: oep_level
 
    push_sub(v_ks_photon_setup)
 
    this%pt => modes
    ! Temporary creation that should go in the system initialization later
    call photon_mode_set_n_electrons(modes, st%qtot)
    write(message(1), '(a,i5,a)') 'Happy to have ', modes%nmodes, ' photon modes with us.'
    call messages_info(1)
    call mf_init(this%pt_mx, gr, st, ions, this%pt)
    ! OEP for photons
    if (bitand(xc_family, xc_family_oep) /= 0 .and. xc%functional(func_x, 1)%id == xc_oep_x) then
      this%oep_photon%level = oep_level
      call xc_oep_photon_init(this%oep_photon, namespace, xc_family, gr, st, mc, space)
    else
      this%oep_photon%level = oep_level_none
    end if
 
    pop_sub(v_ks_photon_setup)
  end subroutine v_ks_photon_setup
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_end(this)
    class(v_ks_photon_t), intent(inout) :: this
 
    push_sub(v_ks_photon_end)
 
    call xc_oep_photon_end(this%oep_photon)
    if (this%has_photons) then
      call mf_end(this%pt_mx)
    end if
    if (this%xc_photon /= 0) then
      call this%xc_photons%end()
    end if
 
    pop_sub(v_ks_photon_end)
  end subroutine v_ks_photon_end
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_mf_calc(this, gr, st, ions, time)
    class(v_ks_photon_t), intent(inout) :: this
    type(grid_t),        intent(inout) :: gr
    type(states_elec_t), intent(inout) :: st
    type(ions_t),        intent(in)    :: ions
    real(real64),        intent(in)    :: time
 
    push_sub(v_ks_photon_mf_calc)
 
    call mf_calc(this%pt_mx, gr, st, ions, this%pt, time)
 
    pop_sub(v_ks_photon_mf_calc)
  end subroutine v_ks_photon_mf_calc
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_add_px(this, namespace, total_density, gr, space, psolver, st, spin_channels, vxc, e_photon_exchange)
    class(v_ks_photon_t), intent(inout) :: this
    type(namespace_t),   intent(in)    :: namespace
    real(real64), pointer, contiguous, intent(in) :: total_density(:)
    type(grid_t),        intent(in)    :: gr
    class(space_t),      intent(in)    :: space
    type(poisson_t),     intent(in)    :: psolver
    type(states_elec_t), intent(inout) :: st
    integer,             intent(in)    :: spin_channels
    real(real64),        intent(inout) :: vxc(:, :)
    real(real64),        intent(out)   :: e_photon_exchange
 
    integer :: ispin
 
    push_sub(v_ks_photon_add_px)
 
    call this%xc_photons%v_ks(namespace, total_density, gr, space, psolver, st)
 
    ! add the photon-free px potential into the xc potential
    do ispin = 1, spin_channels
      call lalg_axpy(gr%np, m_one, this%xc_photons%vpx(1:gr%np), vxc(1:gr%np, ispin))
    end do
 
    ! photon-exchange energy
    e_photon_exchange = this%xc_photons%ex
 
    pop_sub(v_ks_photon_add_px)
  end subroutine v_ks_photon_add_px
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_oep_calc(this, namespace, xc, gr, hm, st, space, vxc)
    class(v_ks_photon_t),      intent(inout) :: this
    type(namespace_t),        intent(in)    :: namespace
    type(xc_t),               intent(inout) :: xc
    type(grid_t),             intent(in)    :: gr
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(states_elec_t),      intent(inout) :: st
    class(space_t),           intent(in)    :: space
    real(real64), contiguous, intent(inout) :: vxc(:, :)
 
    push_sub(v_ks_photon_oep_calc)
 
    if (this%oep_photon%level /= oep_level_none) then
      if (states_are_real(st)) then
        call dxc_oep_photon_calc(this%oep_photon, namespace, xc, gr, &
          hm, st, space, hm%energy%exchange, hm%energy%correlation, vxc = vxc)
      else
        call zxc_oep_photon_calc(this%oep_photon, namespace, xc, gr, &
          hm, st, space, hm%energy%exchange, hm%energy%correlation, vxc = vxc)
      end if
      hm%energy%photon_exchange = this%oep_photon%pt%ex
    end if
 
    pop_sub(v_ks_photon_oep_calc)
  end subroutine v_ks_photon_oep_calc
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_add_mf_potential(this, gr, vhxc, ispin, photon_forces)
    class(v_ks_photon_t), intent(in)    :: this
    type(grid_t),        intent(in)    :: gr
    real(real64),        intent(inout) :: vhxc(:, :)
    integer,             intent(in)    :: ispin
    real(real64),        intent(out)   :: photon_forces(:)
 
    push_sub(v_ks_photon_add_mf_potential)
 
    if (associated(this%pt_mx%vmf)) then
      call lalg_axpy(gr%np, m_one, this%pt_mx%vmf, vhxc(:, 1))
      if (ispin > unpolarized) then
        call lalg_axpy(gr%np, m_one, this%pt_mx%vmf, vhxc(:, 2))
      end if
    end if
    photon_forces(:) = this%pt_mx%fmf(1:size(photon_forces))
 
    pop_sub(v_ks_photon_add_mf_potential)
  end subroutine v_ks_photon_add_mf_potential
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_add_mean_field(this, gr, space, kpoints, st, time, dt)
    class(v_ks_photon_t), intent(inout) :: this
    type(grid_t),        intent(in)    :: gr
    class(space_t),      intent(in)    :: space
    type(kpoints_t),     intent(in)    :: kpoints
    type(states_elec_t), intent(inout) :: st
    real(real64),        intent(in)    :: time
    real(real64),        intent(in)    :: dt
 
    push_sub(v_ks_photon_add_mean_field)
 
    if (this%xc_photon /= 0) then
      call this%xc_photons%add_mean_field(gr, space, kpoints, st, time, dt)
    end if
 
    pop_sub(v_ks_photon_add_mean_field)
  end subroutine v_ks_photon_add_mean_field
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_restart_dump(this, restart, gr, dt, ierr)
    class(v_ks_photon_t), intent(in)    :: this
    type(restart_t),     intent(in)    :: restart
    type(grid_t),        intent(in)    :: gr
    real(real64),        intent(in)    :: dt
    integer,             intent(inout) :: ierr
 
    integer :: err
 
    push_sub(v_ks_photon_restart_dump)
 
    if (this%has_photons) then
      call mf_photons_dump(restart, this%pt_mx, gr, dt, this%pt, err)
      if (err /= 0) ierr = ierr + 16
    end if
    if (this%xc_photon /= 0) then
      ! photon-free mean field
      call this%xc_photons%mf_dump(restart, err)
      if (err /= 0) ierr = ierr + 32
    end if
 
    pop_sub(v_ks_photon_restart_dump)
  end subroutine v_ks_photon_restart_dump
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_restart_load(this, restart, gr, space, ierr)
    class(v_ks_photon_t), intent(inout) :: this
    type(restart_t),     intent(in)    :: restart
    type(grid_t),        intent(in)    :: gr
    class(space_t),      intent(in)    :: space
    integer,             intent(inout) :: ierr
 
    integer :: err
 
    push_sub(v_ks_photon_restart_load)
 
    if (this%has_photons) then
      call mf_photons_load(restart, this%pt_mx, gr, err)
      if (err /= 0) ierr = ierr + 16
    end if
    if (this%xc_photon /= 0) then
      call this%xc_photons%mf_load(restart, space, err)
      if (err /= 0) ierr = ierr + 32
    end if
 
    pop_sub(v_ks_photon_restart_load)
  end subroutine v_ks_photon_restart_load
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_write_info(this, iunit)
    class(v_ks_photon_t), intent(in) :: this
    integer,             intent(in) :: iunit
 
    push_sub(v_ks_photon_write_info)
 
    if (this%oep_photon%level == oep_level_full .or. this%oep_photon%level == oep_level_kli) then
      write(iunit, '(a)') 'Photon observables:'
      write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon number = ', this%oep_photon%pt%number(1)
      write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon ex. = ', this%oep_photon%pt%ex
      write(iunit, '(1x)')
    end if
 
    pop_sub(v_ks_photon_write_info)
  end subroutine v_ks_photon_write_info
 
  ! ---------------------------------------------------------
  logical pure function v_ks_photon_active(this)
    class(v_ks_photon_t), intent(in) :: this
    v_ks_photon_active = this%has_photons
  end function v_ks_photon_active
 
  ! ---------------------------------------------------------
  subroutine v_ks_photon_set_active(this, val)
    class(v_ks_photon_t), intent(inout) :: this
    logical,             intent(in)    :: val
    this%has_photons = val
  end subroutine v_ks_photon_set_active
 
  ! ---------------------------------------------------------
  integer pure function v_ks_photon_functional(this)
    class(v_ks_photon_t), intent(in) :: this
    v_ks_photon_functional = this%xc_photon
  end function v_ks_photon_functional
 
  ! ---------------------------------------------------------
  logical pure function v_ks_photon_include_hartree(this)
    class(v_ks_photon_t), intent(in) :: this
    v_ks_photon_include_hartree = this%include_hartree
  end function v_ks_photon_include_hartree
 
  ! ---------------------------------------------------------
  integer pure function v_ks_photon_nmodes(this)
    class(v_ks_photon_t), intent(in) :: this
    v_ks_photon_nmodes = this%pt%nmodes
  end function v_ks_photon_nmodes
 
  ! ---------------------------------------------------------
  logical pure function v_ks_photon_oep_full(this)
    class(v_ks_photon_t), intent(in) :: this
    v_ks_photon_oep_full = (this%oep_photon%level == oep_level_full)
  end function v_ks_photon_oep_full
 
  ! ---------------------------------------------------------
  function v_ks_photon_correlator(this) result(corr)
    class(v_ks_photon_t), target, intent(in) :: this
    real(real64), pointer :: corr(:)
    corr => this%oep_photon%pt%correlator(:, 1)
  end function v_ks_photon_correlator
 
  ! ---------------------------------------------------------
  function v_ks_photon_get_pt_mx(this) result(pt_mx)
    class(v_ks_photon_t), target, intent(in) :: this
    type(mf_t), pointer :: pt_mx
    pt_mx => this%pt_mx
  end function v_ks_photon_get_pt_mx
 
  ! ---------------------------------------------------------
  function v_ks_photon_get_xc_photons(this) result(xcp)
    class(v_ks_photon_t), target, intent(in) :: this
    type(xc_photons_t), pointer :: xcp
    xcp => this%xc_photons
  end function v_ks_photon_get_xc_photons
 
end module v_ks_photon_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: