xc_photons.F90

Go to the documentation of this file.

!! Copyright (C) 2022 I.-T Lu
!!
!! 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_photons_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use epot_oct_m
  use global_oct_m
  use grid_oct_m
  use kpoints_oct_m
  use lalg_basic_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use photon_mode_oct_m
  use poisson_oct_m
  use restart_oct_m
  use space_oct_m
  use states_elec_oct_m
  use parser_oct_m
  use profiling_oct_m
  use solvers_oct_m
 
  implicit none
 
  private
  public :: xc_photons_t
 
  !
  type xc_photons_t
    private
    integer                      :: method = 0
    real(real64),   allocatable, public :: vpx(:)
    real(real64),         public :: ex
    type(photon_mode_t)          :: pt
    real(real64)                 :: pxlda_kappa
    real(real64)                 :: eta_c
    integer                      :: energy_method = 0
    logical                      :: lcorrelations = .false.    
    logical                      :: llamb_re_mass = .false.    
    logical                      :: llamb_freespace =.false.   
    real(real64)                 :: lamb_omega
 
    logical,              public :: lpfmf = .false. 
    real(real64),   allocatable, public :: mf_vector_potential(:)
    real(real64),   allocatable         :: jp_proj_eo(:,:)
 
  contains
    procedure :: init => xc_photons_init
    procedure :: end => xc_photons_end
    procedure :: wants_to_renormalize_mass => xc_photons_wants_to_renormalize_mass
    procedure :: get_renormalized_mass => xc_photons_get_renormalized_emass
    procedure :: mf_dump => xc_photons_mf_dump
    procedure :: mf_load => xc_photons_mf_load
    procedure :: v_ks => xc_photons_v_ks
    procedure :: add_mean_field => xc_photons_add_mean_field
  end type xc_photons_t
 
  ! the PhotonXCXCMethod
  integer, private, parameter ::  &
    XC_PHOTONS_NONE    = 0,       &
    xc_photons_lda     = 1,       &
    xc_photons_wfs     = 2
 
contains
 
  ! ---------------------------------------------------------
  !
  subroutine xc_photons_init(xc_photons, namespace, xc_photon, space, gr, st)
    class(xc_photons_t),     intent(out)   :: xc_photons
    type(namespace_t),       intent(in)    :: namespace
    integer,                 intent(in)    :: xc_photon
    class(space_t),          intent(in)    :: space
    type(grid_t),            intent(in)    :: gr
    type(states_elec_t),     intent(in)    :: st
 
    push_sub(xc_photons_init)
 
    xc_photons%lpfmf = .false.
 
    call messages_experimental("XCPhotonFunctional /= none")
 
    call photon_mode_init(xc_photons%pt, namespace, space%dim, .true.)
    call photon_mode_set_n_electrons(xc_photons%pt, st%qtot)
 
    select case(xc_photon)
    case(option__xcphotonfunctional__photon_x_lda)
      xc_photons%method    = xc_photons_lda
      xc_photons%lcorrelations = .false.
    case(option__xcphotonfunctional__photon_xc_lda)
      xc_photons%method    = xc_photons_lda
      xc_photons%lcorrelations = .true.
    case(option__xcphotonfunctional__photon_x_wfn)
      xc_photons%method    = xc_photons_wfs
      xc_photons%lcorrelations = .false.
    case(option__xcphotonfunctional__photon_xc_wfn)
      xc_photons%method    = xc_photons_wfs
      xc_photons%lcorrelations = .true.
    case default
      xc_photons%method    = xc_photons_none
      return
    end select
 
 
    if (xc_photons%method == xc_photons_lda) then
 
      !%Variable PhotonXCLDAKappa
      !%Type float
      !%Default 1.0
      !%Section Hamiltonian::XC
      !%Description
      !% the scaling factor for px-LDA potential
      !%End
      call parse_variable(namespace, 'PhotonXCLDAKappa', m_one, xc_photons%pxlda_kappa)
 
    end if
 
    !%Variable PhotonXCEnergyMethod
    !%Type integer
    !%Default 1
    !%Section Hamiltonian::XC
    !%Description
    !% There are different ways to calculate the energy,
    !%Option virial    1
    !% (modified) virial approach</br>
    !% <math>
    !%   (E_{\rm{px}}^{\rm{virial}} = \frac{1}{2}\int d\mathbf{r}\ \mathbf{r}\cdot[
    !%   -\rho(\mathbf{r})\nabla v_{\rm{px}}(\mathbf{r})])
    !%  </math></br>
    !%Option expectation_value 2
    !% expectation value w.tr.t. the wave functions (valid only for 1 electron)</br>
    !% <math>
    !%   E_{\rm{px}}[\rho] = -\sum_{\alpha=1}^{M_{p}}\frac{\tilde{\lambda}_{\alpha}^{2}}{2\tilde{\omega}_{\alpha}^{2}}
    !%   \langle (\tilde{\mathbf{{\varepsilon}}}_{\alpha}\cdot\hat{\mathbf{J}}_{\rm{p}})\Phi[\rho]
    !%   | (\tilde{\mathbf{{\varepsilon}}}_{\alpha}\cdot\hat{\mathbf{J}}_{\rm{p}})\Phi[\rho] \rangle
    !% </math></br>
    !%  This option only works for the wave function based electron-photon functionals
    !%Option LDA             3
    !% energy from electron density</br>
    !% <math>
    !%   E_{\rm pxLDA}[\rho] = \frac{-2\pi^{2}}{(d+2)({2V_{d}})^{\frac{2}{d}}}
    !%   \sum_{\alpha=1}^{M_{p}}\frac{\tilde{\lambda}_{\alpha}^{2}}{\tilde{\omega}_{\alpha}^{2}}
    !%   \int d\mathbf{r}\ \rho^{\frac{2+d}{d}}(\mathbf{r})
    !% </math></br>
    !% This option only works with LDA electron-photon functionals.
    !%End
 
    call parse_variable(namespace, 'PhotonXCEnergyMethod', 1, xc_photons%energy_method)
 
    if( xc_photons%method == xc_photons_wfs .and. xc_photons%energy_method == option__photonxcenergymethod__lda ) then
      message(1) = "Calculating the electron-photon energy from the LDA expression"
      message(2) = "is not implemented for wave function based electron-photon functionals"
      call messages_fatal(2, namespace=namespace)
    end if
 
 
    if (xc_photons%lcorrelations) then
 
      !%Variable PhotonXCEtaC
      !%Type float
      !%Default 1.0
      !%Section Hamiltonian::XC
      !%Description
      !% The scaling factor for the px potential to reduce the weak coupling perturbation regime
      !%End
 
      if (parse_is_defined(namespace, 'PhotonXCEtaC')) then
        call parse_variable(namespace, 'PhotonXCEtaC', m_one, xc_photons%eta_c)
      else
        message(1) = "Defining PhotonXCEtaC is required for photon functionals containing correlation."
        call messages_fatal(1, namespace=namespace)
      end if
 
    else
 
      xc_photons%eta_c = m_one
 
    end if
 
    ! This variable will keep vpx across iterations
    safe_allocate(xc_photons%vpx(1:gr%np_part))
 
    xc_photons%vpx = m_zero
 
    !%Variable PhotonXCLambShift
    !%Type logical
    !%Default .false.
    !%Section Hamiltonian::XC
    !%Description
    !% to deal with the photon free exchange potential for continuum mode in free space
    !%End
 
    call parse_variable(namespace, 'PhotonXCLambShift', .false., xc_photons%llamb_freespace)
    call messages_experimental("PhotonXCLambShift", namespace=namespace)
 
    if (xc_photons%llamb_freespace) then
 
      !%Variable PhotonXCLambShiftOmegaCutoff
      !%Type float
      !%Default 0.0
      !%Section Hamiltonian::XC
      !%Description
      !% the cutoff frequency (Ha) for Lamb shift
      !%End
 
      call parse_variable(namespace, 'PhotonXCLambShiftOmegaCutoff', m_zero, xc_photons%lamb_omega)
 
      !%Variable PhotonXCLambShiftRenormalizeMass
      !%Type logical
      !%Default .false.
      !%Section Hamiltonian::XC
      !%Description
      !% to deal with the photon free exchange potential for continuum mode in free space
      !%End
 
      call parse_variable(namespace, 'PhotonXCLambShiftRenormalizeMass', .false., xc_photons%llamb_re_mass)
 
    end if
 
    ! compute the dressed photon modes
    call photon_mode_dressed(xc_photons%pt)
 
 
    pop_sub(xc_photons_init)
 
  end subroutine xc_photons_init
 
  ! ---------------------------------------------------------
  subroutine xc_photons_end(this)
    class(xc_photons_t), intent(inout) :: this
 
    push_sub(xc_photons_end)
 
    call photon_mode_end(this%pt)
    safe_deallocate_a(this%vpx)
 
    if (allocated(this%mf_vector_potential)) then
      safe_deallocate_a(this%mf_vector_potential)
    end if
    if (allocated(this%jp_proj_eo)) then
      safe_deallocate_a(this%jp_proj_eo)
    end if
 
    pop_sub(xc_photons_end)
  end subroutine xc_photons_end
 
  !
  subroutine xc_photons_v_ks(xc_photons, namespace, total_density, gr, space, psolver, ep, st)
    class(xc_photons_t),              intent(inout) :: xc_photons
    type(namespace_t),                intent(in)    :: namespace
    real(real64), pointer, contiguous,       intent(in)    :: total_density(:)
    class(grid_t),                    intent(in)    :: gr
    type(space_t),                    intent(in)    :: space
    type(poisson_t),                  intent(in)    :: psolver
    type(epot_t),                     intent(in)    :: ep
    type(states_elec_t),              intent(inout) :: st
 
    integer :: ia
 
    push_sub(xc_photons_v_ks)
 
    xc_photons%lpfmf = xc_photons%method > 0
 
    xc_photons%vpx = m_zero
    xc_photons%ex = m_zero
 
    if ( .not. allocated(xc_photons%mf_vector_potential) ) then
      safe_allocate(xc_photons%mf_vector_potential(1:space%dim))
      xc_photons%mf_vector_potential = m_zero
    end if
    if ( .not. allocated(xc_photons%jp_proj_eo)) then
      safe_allocate(xc_photons%jp_proj_eo(1:xc_photons%pt%nmodes,1:2))
      xc_photons%jp_proj_eo = m_zero
    end if
 
 
    select case(xc_photons%method)
    case(xc_photons_lda)  ! LDA approximation for px potential
      call photon_free_vpx_lda(namespace, xc_photons, total_density, gr, space, psolver)
    case(xc_photons_wfs)  ! wave function approxmation for px potential
      call photon_free_vpx_wfc(namespace, xc_photons, total_density, gr, space, st)
    case(xc_photons_none) ! no photon-exchange potential
      call messages_write('Photon-free px potential is not computed', new_line = .true.)
      call messages_info()
    case default
      assert(.false.)
    end select
 
 
    if (.not. xc_photons%llamb_freespace) then
      ! add the constant energy shift
      do ia = 1, xc_photons%pt%nmodes
        xc_photons%ex = xc_photons%ex + 0.5_real64 * (xc_photons%pt%dressed_omega(ia)-xc_photons%pt%omega(ia))
      end do
    end if
 
    pop_sub(xc_photons_v_ks)
  end subroutine xc_photons_v_ks
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  !
  ! The Lamb shift code is experimental and untested
  subroutine photon_free_vpx_lda(namespace, xc_photons, total_density, gr, space, psolver)
    type(namespace_t),                intent(in)    :: namespace
    type(xc_photons_t),               intent(inout) :: xc_photons
    real(real64), pointer, contiguous,       intent(in)    :: total_density(:)
    class(grid_t),                    intent(in)    :: gr
    type(space_t),                    intent(in)    :: space
    type(poisson_t),                  intent(in)    :: psolver
 
    integer :: ia, ip, iter
    real(real64) :: unit_volume, r, res,  presum, prefact
    real(real64) :: xx(space%dim), prefactor_lamb
    real(real64), allocatable :: prefactor(:)
    real(real64), allocatable :: rho_aux(:)
    real(real64), allocatable :: grad_rho_aux(:,:)
    real(real64), allocatable :: px_source(:)
    real(real64), allocatable :: tmp1(:)
    real(real64), allocatable :: tmp2(:,:)
    real(real64), allocatable :: tmp3(:)
    real(real64), allocatable :: grad_vpx(:,:)
    real(real64), allocatable :: epsgrad_epsgrad_rho_aux(:)
    real(real64), allocatable :: epx_force_module(:)
 
    real(real64), parameter :: threshold = 1e-7_real64
 
    push_sub(photon_free_vpx_lda)
 
    if (xc_photons%energy_method == 2 .and. xc_photons%pt%n_electrons >1) then
      call messages_not_implemented("expectation value for energy for pxLDA more than 1 electron", namespace=namespace)
    end if
 
    xc_photons%vpx = m_zero
    xc_photons%ex = m_zero
 
    safe_allocate(prefactor(1:xc_photons%pt%nmodes))
    prefactor = m_zero
    ! here we will use only one spin channel
    safe_allocate(rho_aux(1:gr%np_part))
    safe_allocate(grad_rho_aux(1:gr%np, 1:xc_photons%pt%dim))
    safe_allocate(px_source(1:gr%np_part))
    safe_allocate(tmp1(1:gr%np_part))
    safe_allocate(tmp2(1:gr%np, 1:xc_photons%pt%dim))
    safe_allocate(tmp3(1:gr%np_part))
    safe_allocate(grad_vpx(1:gr%np, 1:xc_photons%pt%dim))
    grad_vpx = m_zero
    safe_allocate(epx_force_module(1:gr%np_part))
    epx_force_module = m_zero
 
    select case(xc_photons%pt%dim)
    case(1)
      unit_volume = m_two
    case(2)
      unit_volume = m_pi
    case(3)
      unit_volume = m_four*m_pi/m_three
    case default
      call messages_not_implemented("LDA px more than 3 dimension", namespace=namespace)
    end select
 
    !$omp parallel do
    do ip=1, gr%np
      rho_aux(ip) = ( abs(total_density(ip))/(m_two*unit_volume) )**(m_two/(xc_photons%pt%dim*m_one))
    end do
    !$omp end parallel do
 
 
    ! compute the electron-photon exchange potential
 
    if (xc_photons%llamb_freespace) then
 
      ! Note: The Lamb shift part is currently experimental and untested!
      ! compute the electron-photon exchange potential
 
      prefactor_lamb = -(8.0_real64*m_pi*m_third) * xc_photons%lamb_omega / (p_c**3)
 
      !$OMP parallel do
      do ip=1,gr%np
        xc_photons%vpx(ip) = prefactor_lamb*rho_aux(ip)
      end do
      !$OMP end parallel do
 
    else
 
      do ia = 1, xc_photons%pt%nmodes
        prefactor(ia) = -m_two*(m_pi * xc_photons%pt%dressed_lambda(ia) / xc_photons%pt%dressed_omega(ia))**2
      end do
 
      select case(xc_photons%pt%dim)
      case(1)
        ! solve the pxLDA potential using the analytical form in 1D
        px_source = m_zero
        do ia = 1, xc_photons%pt%nmodes
          !$OMP parallel do
          do ip=1,gr%np_part
            px_source(ip) = px_source(ip) +  prefactor(ia)
          end do
          !$OMP end parallel do
        end do
 
        !$OMP parallel do
        do ip=1,gr%np
          xc_photons%vpx(ip) = px_source(ip)*rho_aux(ip)
        end do
        !$OMP end parallel do
      case(2)
        call get_px_source(px_source)
        ! for 2D we solve the Poisson equation using the conjugate gradient method
        mesh_aux => gr%der%mesh
        iter = 1000
        call dconjugate_gradients(gr%np, xc_photons%vpx(:), px_source, laplacian_op, dmf_dotp_aux, iter, res, threshold)
        write(message(1),'(a,i6,a)')  "Info: CG converged with ", iter, " iterations."
        write(message(2),'(a,e14.6)')    "Info: The residue is ", res
        call messages_info(2, namespace=namespace)
 
      case(3)
        ! for 3D we use thepoisson solver including the prefactor (-4*pi)
        ! therefore, we need to remove the factor in advance
        call get_px_source(px_source)
 
        call lalg_scal(gr%np, m_one/(-m_four*m_pi), px_source)
 
        ! solve the Poisson equation
        call dpoisson_solve(psolver, namespace, xc_photons%vpx(:), px_source)
      case default
        assert(.false.)
      end select
 
    end if
 
    ! scaling the potential
    call lalg_scal(gr%np, (xc_photons%eta_c * xc_photons%pxlda_kappa), xc_photons%vpx)
 
    ! compute electron-photon energy
 
    select case (xc_photons%energy_method)
    case(1) ! compute the epx energy from the virial relation
 
      do ia = 1, xc_photons%pt%nmodes
 
        ! compute the electron-photon force
        !$omp parallel do
        do ip = 1, gr%np
          epx_force_module(ip) = -prefactor(ia)*m_two*abs(total_density(ip))*rho_aux(ip)/(xc_photons%pt%dim*m_one+m_two)
        end do
        !$omp end parallel do
 
        call dderivatives_grad(gr%der, epx_force_module(:), tmp2)
        call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, tmp2, xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia), m_zero, tmp1)
 
        !$omp parallel do private(r, xx)
        do ip = 1, gr%np
          call mesh_r(gr, ip, r, coords=xx)
          tmp3(ip) = tmp1(ip)*dot_product(xx(1:xc_photons%pt%dim), xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia))
        end do
        !$omp end parallel do
 
        xc_photons%ex = xc_photons%ex + m_half*dmf_integrate(gr, tmp3)
      end do
 
      xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * xc_photons%ex
 
    case(2) ! compute the energy as expetation value wrt to the wave functions
 
      rho_aux(1:gr%np) = sqrt(abs( total_density(1:gr%np)))
 
      safe_allocate(epsgrad_epsgrad_rho_aux(1:gr%np))
 
      call dderivatives_grad(gr%der, rho_aux(1:gr%np_part), grad_rho_aux)
 
      do ia = 1, xc_photons%pt%nmodes
        prefact = (xc_photons%pt%dressed_lambda(ia)**2) / (m_two*xc_photons%pt%dressed_omega(ia)**2)
 
        call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
        call dderivatives_grad(gr%der, tmp1, tmp2)
        call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, tmp2, xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia), m_zero, tmp1)
 
        !$OMP parallel do
        do ip=1, gr%np
          epsgrad_epsgrad_rho_aux(ip) = epsgrad_epsgrad_rho_aux(ip) + prefact*tmp1(ip)
        end do
        !$OMP end parallel do
 
      end do
 
      xc_photons%ex = xc_photons%eta_c * dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux(1:gr%np))
 
      safe_deallocate_a(epsgrad_epsgrad_rho_aux)
 
    case(3) ! integrate the aux electron density over the volume
 
      !$OMP parallel do
      do ip=1,gr%np
        tmp1(ip) = abs( total_density(ip))**((m_one*xc_photons%pt%dim+m_two)/(xc_photons%pt%dim*m_one))
      end do
      !$OMP end parallel do
 
      ! sum over the prefactors
      presum = m_zero
      do ia = 1, xc_photons%pt%nmodes
        presum = presum + prefactor(ia)
      end do
      presum = presum * (m_one/(m_two*unit_volume))**(m_two/(xc_photons%pt%dim*m_one)) / (xc_photons%pt%dim*m_one+m_two)
      presum = presum * xc_photons%eta_c * xc_photons%pxlda_kappa
 
      xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * presum * dmf_integrate(gr, tmp1)
 
    end select
 
    safe_deallocate_a(prefactor)
    safe_deallocate_a(rho_aux)
    safe_deallocate_a(grad_rho_aux)
    safe_deallocate_a(px_source)
    safe_deallocate_a(tmp1)
    safe_deallocate_a(tmp2)
    safe_deallocate_a(tmp3)
    safe_deallocate_a(grad_vpx)
    safe_deallocate_a(epx_force_module)
 
    pop_sub(photon_free_vpx_lda)
 
  contains
    ! ---------------------------------------------------------
    subroutine laplacian_op(x, hx)
      real(real64), contiguous, intent(in)    :: x(:)
      real(real64), contiguous, intent(out)   :: Hx(:)
 
      real(real64), allocatable :: tmpx(:)
 
      safe_allocate(tmpx(1:gr%np_part))
      call lalg_copy(gr%np, x, tmpx)
      call dderivatives_lapl(gr%der, tmpx, hx)
      safe_deallocate_a(tmpx)
 
    end subroutine laplacian_op
 
    subroutine get_px_source(px_source)
      real(real64), contiguous, intent(out) :: px_source(:)
 
      call dderivatives_grad(gr%der, rho_aux(1:gr%np_part), grad_rho_aux)
 
      px_source = m_zero
      do ia = 1, xc_photons%pt%nmodes
        call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
        call dderivatives_grad(gr%der, tmp1, tmp2)
        call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, tmp2, xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia), m_zero, tmp1)
        call lalg_axpy(gr%np, prefactor(ia), tmp1, px_source)
      end do
 
    end subroutine get_px_source
 
  end subroutine photon_free_vpx_lda
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  !
  !
  ! The Lamb shift code is experimental and untested
 
  subroutine photon_free_vpx_wfc(namespace, xc_photons, total_density, gr, space, st)
    type(namespace_t),                intent(in)    :: namespace
    type(xc_photons_t),               intent(inout) :: xc_photons
    real(real64), pointer, contiguous,       intent(in)    :: total_density(:)
    class(grid_t),                    intent(in)    :: gr
    type(space_t),                    intent(in)    :: space
    type(states_elec_t),              intent(in)    :: st
 
    integer :: ia, ip
    real(real64) :: prefactor_lamb
    real(real64) :: xx(space%dim), r
    real(real64), allocatable :: prefactor(:)
    real(real64), allocatable :: rho_aux(:)
    real(real64), allocatable :: grad_rho_aux(:,:)
    real(real64), allocatable :: grad_vpx(:,:)
    real(real64), allocatable :: epsgrad_epsgrad_rho_aux(:)
    real(real64), allocatable :: tmp1(:)
    real(real64), allocatable :: tmp2(:,:)
    real(real64)              ::  shift
 
    push_sub(photon_free_vpx_wfc)
 
    if (st%d%nspin >1) then
      call messages_not_implemented("PhotonXCXCMethod = wavefunction for polarized and spinor cases", namespace=namespace)
    end if
 
    if (xc_photons%pt%n_electrons >1) then
      call messages_not_implemented("PhotonXCXCMethod = wavefunction for more than 1 electron", namespace=namespace)
    end if
 
    xc_photons%vpx = m_zero
    xc_photons%ex = m_zero
 
    safe_allocate(prefactor(1:xc_photons%pt%nmodes))
    prefactor = m_zero
    safe_allocate(rho_aux(1:gr%np_part))
    rho_aux = m_zero
    safe_allocate(grad_rho_aux(1:gr%np, 1:xc_photons%pt%dim))
    grad_rho_aux = m_zero
    safe_allocate(epsgrad_epsgrad_rho_aux(1:gr%np))
    safe_allocate(grad_vpx(1:gr%np, 1:xc_photons%pt%dim))
    safe_allocate(tmp1(1:gr%np_part))
    safe_allocate(tmp2(1:gr%np_part, 1:xc_photons%pt%dim))
 
 
    rho_aux(1:gr%np) = sqrt(abs( total_density(1:gr%np)))
    !$OMP parallel do
    do ip = 1, gr%np
      rho_aux(ip) = safe_tol(rho_aux(ip),1e-18_real64)
    end do
    !$OMP end parallel do
 
    if (xc_photons%llamb_freespace) then
 
      ! experimental Lamb shift mode
 
      prefactor_lamb = ( m_two/(m_three*m_pi) ) * xc_photons%lamb_omega / p_c**3
      call dderivatives_lapl(gr%der, rho_aux(1:gr%np_part), epsgrad_epsgrad_rho_aux)
      call lalg_scal(gr%np, prefactor_lamb, epsgrad_epsgrad_rho_aux)
 
    else
 
      do ia = 1, xc_photons%pt%nmodes
        prefactor(ia) = (xc_photons%pt%dressed_lambda(ia)**2) / (m_two*xc_photons%pt%dressed_omega(ia)**2)
      end do
 
      call dderivatives_grad(gr%der, rho_aux, grad_rho_aux)
 
      epsgrad_epsgrad_rho_aux = m_zero
      do ia = 1, xc_photons%pt%nmodes
        tmp1 = m_zero
        call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
        call dderivatives_grad(gr%der, tmp1, tmp2)
        call lalg_gemv(gr%np, xc_photons%pt%dim, prefactor(ia), tmp2, xc_photons%pt%dressed_pol(:, ia), &
          m_one, epsgrad_epsgrad_rho_aux)
      end do
 
    end if
 
    !$OMP parallel do
    do ip = 1, gr%np
      xc_photons%vpx(ip) = epsgrad_epsgrad_rho_aux(ip)/rho_aux(ip)
    end do
    !$OMP end parallel do
 
    if(st%eigenval(1,1) < m_huge .and. .not. space%is_periodic()) then
      shift = m_two * st%eigenval(1,1) * prefactor(1)
      !$OMP parallel do
      do ip=1, gr%np
        xc_photons%vpx(ip) = xc_photons%vpx(ip) + shift
      end do
      !$OMP end parallel do
    end if
 
    call lalg_scal(gr%np, xc_photons%eta_c, xc_photons%vpx(:))
 
    select case (xc_photons%energy_method)
    case(1) ! virial
 
      call dderivatives_grad(gr%der, xc_photons%vpx(:), grad_vpx)
      !$OMP parallel do private(r, xx)
      do ip = 1, gr%np
        call mesh_r(gr, ip, r, coords=xx)
        tmp1(ip) = - total_density(ip)*dot_product(xx(1:xc_photons%pt%dim), grad_vpx(ip,1:xc_photons%pt%dim))
      end do
      !$OMP end parallel do
 
      xc_photons%ex = m_half*dmf_integrate(gr, tmp1)
 
    case(2) ! expectation_value
      xc_photons%ex = xc_photons%eta_c * dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux)
 
    case default
      assert(.false.)
    end select
 
    safe_deallocate_a(prefactor)
    safe_deallocate_a(rho_aux)
    safe_deallocate_a(grad_rho_aux)
    safe_deallocate_a(grad_vpx)
    safe_deallocate_a(epsgrad_epsgrad_rho_aux)
    safe_deallocate_a(tmp1)
    safe_deallocate_a(tmp2)
 
    pop_sub(photon_free_vpx_wfc)
 
  end subroutine photon_free_vpx_wfc
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  !
  subroutine xc_photons_add_mean_field(xc_photons, gr, space, kpoints, st, time, dt)
    class(xc_photons_t),              intent(inout) :: xc_photons
    class(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
 
 
    integer :: ia, idir, ispin
    real(real64)   :: pol_dot_jp
    real(real64), allocatable :: current(:,:,:)
    real(real64), allocatable :: jp(:)
 
    push_sub(xc_photons_add_mean_field)
 
    ! compute the dressed photon-free vector potential
    xc_photons%mf_vector_potential = m_zero
 
    safe_allocate(current(1:gr%np_part, 1:space%dim, 1:st%d%nspin))
    current = m_zero
    ! here we use the paramagnetic current; note that the physical current here
    ! only contains the paramagnetic current.
    call states_elec_calc_quantities(gr, st, kpoints, .false., paramagnetic_current = current)
 
    ! compute the paramagnetic current
    safe_allocate(jp(1:space%dim))
    do idir = 1, space%dim
      jp(idir) = m_zero
      do ispin = 1, st%d%spin_channels
        jp(idir) = jp(idir) + dmf_integrate(gr%der%mesh, current(:,idir,ispin))
      end do
    end do
 
    ! update the 'projected' current
    do ia = 1, xc_photons%pt%nmodes
      pol_dot_jp = dot_product(xc_photons%pt%dressed_pol(1:space%dim, ia),jp(1:space%dim))
      xc_photons%jp_proj_eo(ia,1) = xc_photons%jp_proj_eo(ia,1) + &
        cos(xc_photons%pt%dressed_omega(ia)*( time-dt))*pol_dot_jp*dt
      xc_photons%jp_proj_eo(ia,2) = xc_photons%jp_proj_eo(ia,2) + &
        sin(xc_photons%pt%dressed_omega(ia)*( time-dt))*pol_dot_jp*dt
    end do
 
    do ia = 1, xc_photons%pt%nmodes
      xc_photons%mf_vector_potential(1:xc_photons%pt%dim) = xc_photons%mf_vector_potential(1:xc_photons%pt%dim) &
        + (-p_c*(xc_photons%pt%dressed_lambda(ia)**2) / xc_photons%pt%dressed_omega(ia)) &
        * (xc_photons%jp_proj_eo(ia,1)*sin(xc_photons%pt%dressed_omega(ia)* time) &
        - xc_photons%jp_proj_eo(ia,2)*cos(xc_photons%pt%dressed_omega(ia)* time)) &
        * xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia)
    end do
 
    safe_deallocate_a(current)
    safe_deallocate_a(jp)
 
    pop_sub(xc_photons_add_mean_field)
 
  end subroutine xc_photons_add_mean_field
  ! ---------------------------------------------------------
 
 
  !
  logical pure function xc_photons_wants_to_renormalize_mass(xc_photons) result (renorm)
    class(xc_photons_t), intent(in) :: xc_photons
 
    renorm = (xc_photons%method>0) .and. xc_photons%llamb_freespace .and. xc_photons%llamb_re_mass
 
  end function xc_photons_wants_to_renormalize_mass
 
  real(real64) pure function xc_photons_get_renormalized_emass(xc_photons) result(mass)
    class(xc_photons_t), intent(in) :: xc_photons
 
    mass = m_one - (m_four*xc_photons%lamb_omega) / (3.0*m_pi * p_c**3)
 
  end function xc_photons_get_renormalized_emass
 
  !
  subroutine xc_photons_mf_dump(xc_photons, restart, ierr)
    class(xc_photons_t),      intent(in)  :: xc_photons
    type(restart_t),          intent(in)  :: restart
    integer,                  intent(out) :: ierr
 
    character(len=80), allocatable :: lines(:)
    integer :: iunit, err, jj, nmodes
    integer :: pt_dim
 
    push_sub(photon_free_mf_dump)
    nmodes = xc_photons%pt%nmodes
    pt_dim = xc_photons%pt%dim
 
    safe_allocate(lines(1:nmodes+pt_dim))
 
    ierr = 0
 
    iunit = restart_open(restart, 'photon_free_mf')
 
    do jj = 1, pt_dim
      write(lines(jj), '(2x, es19.12)') xc_photons%mf_vector_potential(jj)
    end do
 
    do jj = 1, nmodes
      write(lines(jj+pt_dim), '(a10,1x,I8,a1,2x,2(es19.12,2x))') 'Mode ', jj, ":", xc_photons%jp_proj_eo(jj,1:2)
    end do
 
    call restart_write(restart, iunit, lines, nmodes+pt_dim, err)
    if (err /= 0) ierr = ierr + 1
    call restart_close(restart, iunit)
 
    safe_deallocate_a(lines)
 
    pop_sub(photon_free_mf_dump)
  end subroutine xc_photons_mf_dump
 
! ---------------------------------------------------------
 
  !
  subroutine xc_photons_mf_load(xc_photons, restart, space, ierr)
    class(xc_photons_t),      intent(inout) :: xc_photons
    type(restart_t),          intent(in)    :: restart
    class(space_t),           intent(in)    :: space
    integer,                  intent(out)   :: ierr
 
    character(len=80), allocatable :: lines(:)
    character(len=7) :: sdummy
    integer :: idummy
    integer :: iunit, err, jj, nmodes
    integer :: pt_dim
 
    push_sub(photon_free_mf_load)
 
    ierr = 0
    nmodes = xc_photons%pt%nmodes
    pt_dim = xc_photons%pt%dim
 
    if (restart_skip(restart)) then
      ierr = -1
      pop_sub(photon_free_mf_load)
      return
    end if
 
    message(1) = "Debug: Reading Photon-Free Photons restart."
    call messages_info(1, namespace=restart%namespace, debug_only=.true.)
 
    if ( .not. allocated(xc_photons%jp_proj_eo)) then
      safe_allocate(xc_photons%jp_proj_eo(1:xc_photons%pt%nmodes, 1:2))
      xc_photons%jp_proj_eo = m_zero
    end if
    if ( .not. allocated(xc_photons%mf_vector_potential)) then
      safe_allocate(xc_photons%mf_vector_potential(1:space%dim))
      xc_photons%mf_vector_potential = m_zero
    end if
 
    safe_allocate(lines(1:nmodes+pt_dim))
    iunit = restart_open(restart, 'photon_free_mf')
    call restart_read(restart, iunit, lines, nmodes+pt_dim, err)
    if (err /= 0) then
      ierr = ierr + 1
    else
      do jj = 1, pt_dim
        read(lines(jj),'(2x, es19.12)') xc_photons%mf_vector_potential(jj)
      end do
 
      do jj = 1, nmodes
        read(lines(jj+pt_dim), '(a10,1x,I8,a1,2x,2(es19.12,2x))') sdummy, idummy, sdummy, xc_photons%jp_proj_eo(jj,1:2)
      end do
    end if
    call restart_close(restart, iunit)
 
    message(1) = "Debug: Reading Photons restart done."
    call messages_info(1, namespace=restart%namespace, debug_only=.true.)
 
    safe_deallocate_a(lines)
 
    pop_sub(photon_free_mf_load)
  end subroutine xc_photons_mf_load
 
end module xc_photons_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: