photon_mode.F90

Go to the documentation of this file.

!! Copyright (C) 2017 Johannes Flick
!! Copyright (C) 2019 F. Buchholz, M. Oliveira
!! Copyright (C) 2020 S. Ohlmann
!!
!! 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 photon_mode_oct_m
  use comm_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use lalg_adv_oct_m
  use, intrinsic :: iso_fortran_env
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                       &
    photon_mode_t,                &
    photon_mode_init,             &
    photon_mode_end,              &
    photon_mode_write_info,       &
    photon_mode_add_poisson_terms,&
    photon_mode_compute_dipoles,  &
    photon_mode_set_n_electrons,  &
    photon_mode_dressed
 
  type photon_mode_t
    ! All components are public by default
    integer               :: nmodes
    integer               :: dim
    real(real64), allocatable    :: omega(:)
    real(real64), allocatable    :: lambda(:)
    real(real64), allocatable    :: pol(:,:)
    real(real64), allocatable    :: pol_dipole(:,:)
    real(real64)          :: ex
    real(real64), allocatable    :: number(:)
    real(real64), allocatable    :: correlator(:,:)
    real(real64)          :: n_electrons
    real(real64), pointer        :: pt_coord_q0(:) => null()   
    real(real64), pointer        :: pt_momen_p0(:) => null()   
    real(real64)          :: mu
    ! dressed photons for photon-free approach
    logical               :: use_photon_exchange
    real(real64), allocatable    :: dressed_omega(:)
    real(real64), allocatable    :: dressed_lambda(:)
    real(real64), allocatable    :: dressed_pol(:,:)
    real(real64), allocatable    :: rotated_omega(:)
 
    logical               :: has_q0_p0 = .false.
  end type photon_mode_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine photon_mode_init(this, namespace, dim, photon_free)
    type(photon_mode_t),  intent(out) :: this
    type(namespace_t),    intent(in)  :: namespace
    integer,              intent(in)  :: dim
    logical, optional,    intent(in)  :: photon_free
 
    type(block_t)         :: blk
    integer               :: ii, idir, iunit, ncols
    logical               :: file_exists
    character(256)        :: filename
 
    push_sub(photon_mode_init)
 
    this%nmodes = 0
 
    this%dim = dim
    this%has_q0_p0 = .false.
 
    !%Variable PhotonmodesFilename
    !%Type string
    !%Default "photonmodes"
    !%Section Linear Response::Casida
    !%Description
    !% Filename for photon modes in text format
    !%  - first line contains 2 integers: number of photon modes and number of
    !%    columns
    !%  - each further line contains the given number of floats for one photon
    !%    mode
    !%End
    call parse_variable(namespace, 'PhotonmodesFilename', 'photonmodes', filename)
    inquire(file=trim(filename), exist=file_exists)
    if (file_exists) then
      if (mpi_grp_is_root(mpi_world)) then
        message(1) = 'Opening '//trim(filename)
        call messages_info(1, namespace=namespace)
        ! open file on root
        iunit = io_open(trim(filename), namespace, action='read', form='formatted')
 
        ! get dimensions from first line
        read(iunit, *) this%nmodes, ncols
 
        write(message(1), '(3a,i7,a,i3,a)') 'Reading file ', trim(filename), ' with ', &
          this%nmodes, ' photon modes and ', ncols, ' columns.'
        call messages_info(1, namespace=namespace)
 
        safe_allocate(this%omega(1:this%nmodes))
        safe_allocate(this%lambda(1:this%nmodes))
        safe_allocate(this%pol(1:this%nmodes,1:3))
 
        ! now read in all modes
        do ii = 1, this%nmodes
          if (ncols == 5) then
            read(iunit, *) this%omega(ii), this%lambda(ii), &
              this%pol(ii,1), this%pol(ii,2), this%pol(ii,3)
          else if (ncols == 7) then
            this%has_q0_p0 = .true.
            if (.not. associated(this%pt_coord_q0)) then
              safe_allocate(this%pt_coord_q0(1:this%nmodes))
            end if
            if (.not. associated(this%pt_momen_p0)) then
              safe_allocate(this%pt_momen_p0(1:this%nmodes))
            end if
            read(iunit, *) this%omega(ii), this%lambda(ii), &
              this%pol(ii,1), this%pol(ii,2), this%pol(ii,3), &
              this%pt_coord_q0(ii), this%pt_momen_p0(ii)
          else
            ! error if not 5 columns
            message(1) = 'Error: unexpected number of columns in file:'
            message(2) = filename
            call messages_fatal(2, namespace=namespace)
          end if
 
          ! Normalize polarization vector
          this%pol(ii,:) = this%pol(ii,:)/norm2(this%pol(ii,:))
        end do
        call io_close(iunit)
      end if
      ! broadcast first array dimensions, then allocate and broadcast arrays
      call mpi_world%bcast(this%nmodes, 1, mpi_integer, 0)
      call mpi_world%bcast(ncols, 1, mpi_integer, 0)
      if (.not. mpi_grp_is_root(mpi_world)) then
        safe_allocate(this%omega(1:this%nmodes))
        safe_allocate(this%lambda(1:this%nmodes))
        safe_allocate(this%pol(1:this%nmodes,1:3))
        if (ncols == 7) then
          safe_allocate(this%pt_coord_q0(1:this%nmodes))
          safe_allocate(this%pt_momen_p0(1:this%nmodes))
        end if
      end if
      call mpi_world%bcast(this%omega(1), this%nmodes, mpi_double_precision, 0)
      call mpi_world%bcast(this%lambda(1), this%nmodes, mpi_double_precision, 0)
      call mpi_world%bcast(this%pol(1,1), this%nmodes*3, mpi_double_precision, 0)
      if (ncols == 7) then
        call mpi_world%bcast(this%pt_coord_q0(1), this%nmodes, mpi_double_precision, 0)
        call mpi_world%bcast(this%pt_momen_p0(1), this%nmodes, mpi_double_precision, 0)
      end if
    else
      if (.not. parse_is_defined(namespace, 'PhotonModes')) then
        call messages_write('You need to specify the correct external photon modes file')
        call messages_write('or define the PhotonModes variable!')
        call messages_fatal(namespace=namespace)
      end if
    end if
 
    !%Variable PhotonModes
    !%Type block
    !%Section Hamiltonian::XC
    !%Description
    !% Each line of the block should specify one photon mode. The syntax is the following:
    !%
    !% %PhotonModes
    !%  omega1 | lambda1| PolX1 | PolY1 | PolZ1
    !%  ...
    !% %
    !%
    !% The first column is the mode frequency, in units of energy.
    !% The second column is the coupling strength, in units of energy.
    !% The remaining columns specify the polarization direction of the mode.
    !% If the polarization vector should be normalized to one. If that is not the case
    !% the code will normalize it.
    !%End
 
    if (.not. file_exists) then
      if (parse_block(namespace, 'PhotonModes', blk) == 0) then
 
        this%nmodes = parse_block_n(blk)
        safe_allocate(this%omega(1:this%nmodes))
        safe_allocate(this%lambda(1:this%nmodes))
        safe_allocate(this%pol(1:this%nmodes, 1:this%dim))
 
        this%pol = m_zero
 
        do ii = 1, this%nmodes
          ncols = parse_block_cols(blk, ii-1)
 
          ! Read line
          call parse_block_float(blk, ii-1, 0, this%omega(ii), units_inp%energy)  ! frequency
          call parse_block_float(blk, ii-1, 1, this%lambda(ii), units_inp%energy) ! coupling strength
          do idir = 1, this%dim
            call parse_block_float(blk, ii-1, idir + 1, this%pol(ii, idir)) ! polarization vector components
          end do
 
          if (ncols > this%dim + 2) then
            this%has_q0_p0 = .true.
            if (.not. associated(this%pt_coord_q0)) then
              safe_allocate(this%pt_coord_q0(1:this%nmodes))
            end if
            if (.not. associated(this%pt_momen_p0)) then
              safe_allocate(this%pt_momen_p0(1:this%nmodes))
            end if
            call parse_block_float(blk, ii-1, this%dim+2, this%pt_coord_q0(ii))   !row, column
            call parse_block_float(blk, ii-1, this%dim+3, this%pt_momen_p0(ii))   !row, column
          end if
          ! Sanity check
          if ((ncols /= this%dim + 2) .and. (ncols /= this%dim + 4)) then
            call messages_input_error(namespace, 'PhotonModes', 'Incorrect number of columns')
          end if
 
          ! Normalize polarization vector
          this%pol(ii,:) = this%pol(ii,:)/norm2(this%pol(ii,:))
 
        end do
        call parse_block_end(blk)
      else
        call messages_write('You need to specify the photon modes!')
        call messages_fatal(namespace=namespace)
      end if
    end if
 
    !%Variable TDPhotonicTimeScale
    !%Type float
    !%Default 1.0
    !%Section Time-Dependent::Propagation
    !%Description
    !% This variable defines the factor between the timescale of photonic
    !% and electronic movement.
    !% for more details see the documentation of TDIonicTimeScale
    !% If you also use TDIonicTimeScale, we advise to set
    !% TDPhotonicTimeScale = TDIonicTimeScale, in the case the
    !% photon frequency is in a vibrational energy range.
    !% Important: The electronic time step will be the value of
    !% <tt>TDTimeStep</tt> divided by this variable, so if you have determined an
    !% optimal electronic time step (that we can call <i>dte</i>), it is
    !% recommended that you define your time step as:
    !%
    !% <tt>TDTimeStep</tt> = <i>dte</i> * <tt>TDPhotonicTimeScale</tt>
    !%
    !% so you will always use the optimal electronic time step
    !% (<a href=http://arxiv.org/abs/0710.3321>more details</a>).
    !%End
 
    call parse_variable(namespace, 'TDPhotonicTimeScale', 1.0_real64, this%mu)
 
    this%ex = m_zero
    safe_allocate(this%number(1:this%nmodes))
    this%number = m_zero
 
    ! photon-free approach, which requires the dressed photon modes
    this%use_photon_exchange = .false.
 
    this%use_photon_exchange = optional_default(photon_free, .false.)
    if (this%use_photon_exchange) then
      safe_allocate(this%dressed_omega(1:this%nmodes))
      safe_allocate(this%dressed_lambda(1:this%nmodes))
      safe_allocate(this%dressed_pol(1:this%dim, 1:this%nmodes))
      safe_allocate(this%rotated_omega(1:this%nmodes))
    end if
 
    pop_sub(photon_mode_init)
  end subroutine photon_mode_init
 
  ! ---------------------------------------------------------
  subroutine photon_mode_end(this)
    type(photon_mode_t), intent(inout) :: this
 
    push_sub(photon_mode_end)
 
    safe_deallocate_a(this%correlator)
 
    safe_deallocate_a(this%omega)
    safe_deallocate_a(this%lambda)
    safe_deallocate_a(this%number)
 
    safe_deallocate_a(this%pol)
    safe_deallocate_a(this%pol_dipole)
 
    if (this%has_q0_p0) then
      safe_deallocate_p(this%pt_coord_q0)
      safe_deallocate_p(this%pt_momen_p0)
    end if
 
    if (this%use_photon_exchange) then
      safe_deallocate_a(this%dressed_omega)
      safe_deallocate_a(this%dressed_lambda)
      safe_deallocate_a(this%dressed_pol)
      safe_deallocate_a(this%rotated_omega)
    end if
 
    pop_sub(photon_mode_end)
  end subroutine photon_mode_end
 
  !-----------------------------------------------------------------
  subroutine photon_mode_write_info(this, iunit, namespace)
    type(photon_mode_t),           intent(in) :: this
    integer,             optional, intent(in) :: iunit
    type(namespace_t),   optional, intent(in) :: namespace
 
    integer :: im, idir
 
    push_sub(photon_mode_write_info)
 
    call messages_print_with_emphasis(msg="Photon Modes", iunit=iunit, namespace=namespace)
    write(iunit, '(6x,a,10x,a,3x)', advance='no') 'Omega', 'Lambda'
    write(iunit, '(1x,a,i1,a)') ('Pol.(', idir, ')', idir = 1, this%dim)
    do im = 1, this%nmodes
      write(iunit, '(1x,f14.12)', advance='no') this%omega(im)
      write(iunit, '(1x,f14.12)', advance='no') this%lambda(im)
      write(iunit, '(2X,f5.3,1X)') (this%pol(im, idir), idir = 1, this%dim)
    end do
    call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
 
    pop_sub(photon_mode_write_info)
  end subroutine photon_mode_write_info
 
  !-----------------------------------------------------------------
  subroutine photon_mode_set_n_electrons(this, qtot)
    type(photon_mode_t), intent(inout)  :: this
    real(real64),        intent(in)     :: qtot
 
    push_sub(photon_mode_set_n_electrons)
 
    this%n_electrons = qtot
 
    pop_sub(photon_mode_set_n_electrons)
  end subroutine photon_mode_set_n_electrons
 
  !-----------------------------------------------------------------
  subroutine photon_mode_compute_dipoles(this, mesh)
    type(photon_mode_t), intent(inout)  :: this
    class(mesh_t),       intent(in)     :: mesh
 
    integer :: i, ip
 
    push_sub(photon_mode_compute_dipoles)
 
    safe_allocate(this%pol_dipole(1:mesh%np, 1:this%nmodes))
 
    do i = 1, this%nmodes
      ! Calculate polarization dipole
      do ip = 1, mesh%np
        this%pol_dipole(ip, i) = dot_product(this%pol(i, 1:this%dim), mesh%x(ip, 1:this%dim))
      end do
    end do
 
    pop_sub(photon_mode_compute_dipoles)
  end subroutine photon_mode_compute_dipoles
 
  !-----------------------------------------------------------------
  subroutine photon_mode_add_poisson_terms(this, mesh, rho, pot)
    type(photon_mode_t), intent(in)    :: this
    type(mesh_t),        intent(in)    :: mesh
    real(real64),        intent(in)    :: rho(:)
    real(real64),        intent(inout) :: pot(:)
 
    integer :: ip
    real(real64) :: lx, ld, dipole(mesh%box%dim)
 
    push_sub(photon_mode_add_poisson_terms)
 
    ! Currently this only works with one photon mode
    assert(this%nmodes == 1)
    assert(this%n_electrons > m_epsilon)
 
    call dmf_dipole(mesh, rho, dipole)
    ld = dot_product(this%pol(1, 1:this%dim), dipole(1:this%dim))*this%lambda(1)
 
    do ip = 1, mesh%np
      lx = this%pol_dipole(ip, 1)*this%lambda(1)
      pot(ip) = pot(ip) - this%omega(1)/sqrt(this%n_electrons)*(mesh%x(ip, this%dim + 1)*ld + lx*dipole(this%dim + 1)) + lx*ld
    end do
 
    pop_sub(photon_mode_add_poisson_terms)
  end subroutine photon_mode_add_poisson_terms
 
  !-----------------------------------------------------------------
  subroutine photon_mode_dressed(this)
    type(photon_mode_t), intent(inout)    :: this
 
    integer :: ia, iap
    real(real64) :: eps_epsp, wmat_off
    real(real64), allocatable :: wmat(:,:)
    real(real64), allocatable :: dressed_omega_sq(:)
    real(real64), allocatable :: lambda_pol(:,:)
    real(real64), allocatable :: dressed_lambda_pol(:,:)
 
    push_sub(photon_mode_dressed)
 
    safe_allocate(wmat(1:this%nmodes,1:this%nmodes))
    safe_allocate(dressed_omega_sq(1:this%nmodes))
    safe_allocate(lambda_pol(1:this%nmodes,1:this%dim))
    safe_allocate(dressed_lambda_pol(1:this%dim, 1:this%nmodes))
 
    ! compute \lambda**2 * polarization array
    do ia = 1, this%nmodes
      lambda_pol(ia, 1:this%dim) = this%lambda(ia) * this%pol(ia, 1:this%dim)
    end do
 
    ! compute the photon mode matrix, the W matrix
    do ia = 1, this%nmodes
      do iap = ia, this%nmodes
        eps_epsp = dot_product(this%pol(ia, 1:this%dim), this%pol(iap, 1:this%dim))
        ! TODO: n_electrons should be the total number of electrons including the core levels
        wmat_off = this%n_electrons * this%lambda(ia) * this%lambda(iap) * eps_epsp
        wmat(ia,iap) = wmat_off
        wmat(iap,ia) = wmat_off
      end do
      wmat(ia,ia) = this%omega(ia)**2 + wmat(ia,ia)
    end do
 
    ! diagonalize the W matrix
    call lalg_eigensolve(this%nmodes, wmat, dressed_omega_sq)
 
    do ia = 1, this%nmodes
      this%dressed_omega(ia) = sqrt(dressed_omega_sq(ia))
    end do
 
    ! rotate the ( coupling times polarization )
    call dgemm('T','N', this%dim, this%nmodes, this%nmodes,&
      m_one, lambda_pol, this%nmodes, &
      wmat,this%nmodes, &
      m_zero,dressed_lambda_pol,this%dim)
 
    ! compute the rotated coupling and polarization direction
    do ia = 1, this%nmodes
      this%dressed_lambda(ia) = norm2(dressed_lambda_pol(1:this%dim, ia))
      !
      if (this%dressed_lambda(ia) .lt. m_epsilon) THEN
        ! to avoid dividing by zero
        this%dressed_pol(1:this%dim, ia) = m_zero
      else
        this%dressed_pol(1:this%dim, ia) = dressed_lambda_pol(1:this%dim, ia)/this%dressed_lambda(ia)
      end if
    end do
 
    ! rotate the bare frequency
    this%rotated_omega = m_zero
    ! here we define the rotated bare frequency from the dressed frequency and coupling
    do ia = 1, this%nmodes
      this%rotated_omega(ia) = sqrt(this%dressed_omega(ia)**2 - this%n_electrons * this%dressed_lambda(ia)**2)
    end do
 
    safe_deallocate_a(wmat)
    safe_deallocate_a(dressed_omega_sq)
    safe_deallocate_a(lambda_pol)
    safe_deallocate_a(dressed_lambda_pol)
 
    pop_sub(photon_mode_dressed)
 
  end subroutine photon_mode_dressed
 
end module photon_mode_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: