dielectric_function.F90

Go to the documentation of this file.

!! Copyright (C) 2008 X. Andrade
!!
!! 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 dielectric_function_main_oct_m
  use batch_oct_m
  use command_line_oct_m
  use global_oct_m
  use io_oct_m
  use lalg_adv_oct_m
  use lattice_vectors_oct_m
  use math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use spectrum_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public :: dielectric_function_main
 
contains
 
  subroutine dielectric_function_main() bind(c)
    integer :: ierr
 
    ! Initialize stuff
    call init_octopus_globals(serial_dummy_comm)
 
    call getopt_init(ierr)
    if (ierr == 0) call getopt_dielectric_function()
    call getopt_end()
 
    call parser_init()
 
    call messages_init()
 
    call io_init()
    call profiling_init(global_namespace)
 
    call unit_system_init(global_namespace)
 
    call dielectric_function_compute()
 
    call profiling_end(global_namespace)
    call io_end()
    call messages_end()
 
    call parser_end()
    call global_end()
 
  contains
 
    subroutine dielectric_function_compute()
 
      integer :: in_file, out_file, ref_file, ii, jj, kk
      integer :: time_steps, time_steps_ref, energy_steps, istart, iend, ntiter
      real(real64)   :: dt, dt_ref, tt, ww, norm_Aext
      real(real64), allocatable :: vecpot(:, :), Aext(:), Eind_real(:, :), Eind_imag(:, :)
      complex(real64), allocatable :: dielectric(:), chi(:), invdielectric(:)
      real(real64), allocatable :: vecpot_ref(:, :)
      type(spectrum_t)  :: spectrum
      type(block_t)     :: blk
      type(space_t)     :: space
      type(lattice_vectors_t) :: latt
      type(batch_t)     :: vecpotb, Eind_realb, Eind_imagb
      character(len=120) :: header
      real(real64) :: start_time
      character(len=MAX_PATH_LEN) :: ref_filename
      complex(real64),allocatable :: Eind(:)
 
      call spectrum_init(spectrum, global_namespace)
 
      space = space_t(global_namespace)
      latt = lattice_vectors_t(global_namespace, space)
 
      safe_allocate(aext(1:space%dim))
      safe_allocate(eind(1:space%dim))
 
      if (parse_block(global_namespace, 'GaugeVectorField', blk) == 0) then
 
        do ii = 1, space%dim
          call parse_block_float(blk, 0, ii - 1, aext(ii))
        end do
 
        call parse_block_end(blk)
 
      else
 
        message(1) = "Cannot find the GaugeVectorField in the input file"
        call messages_fatal(1)
 
      end if
 
      if(space%dim > 1) then
        message(1) = "This program assumes that the gauge field is along a  "
        message(2) = "direction for which the dielectric tensor has no off-diagonal terms."
        message(3) = "If this is not the case the dielectric function and the"
        message(4) = "susceptibility will be wrong."
        call messages_warning(4)
      end if
 
      start_time = spectrum%start_time
      call parse_variable(global_namespace, 'GaugeFieldDelay', start_time, spectrum%start_time)
 
      in_file = io_open('td.general/gauge_field', global_namespace, action='read', status='old')
 
      ! Get the number of iterations and the time-step
      call io_skip_header(in_file)
      call spectrum_count_time_steps(global_namespace, in_file, time_steps, dt)
 
      ! Fix the correct time at which the kick was applied.
      ! This guaranties that the first time has exactly a zero vector potential
      ! If we do not do this, we can get some artifacts for transient absorption.
      spectrum%start_time = ceiling(spectrum%start_time/dt)*dt
 
      if (parse_is_defined(global_namespace, 'TransientAbsorptionReference')) then
        !%Variable TransientAbsorptionReference
        !%Type string
        !%Default "."
        !%Section Utilities::oct-propagation_spectrum
        !%Description
        !% In case of delayed kick, the calculation of the transient absorption requires
        !% to substract a reference calculation, containing the gauge-field without the kick
        !% This reference must be computed using GaugeFieldPropagate=yes and to have
        !% TDOutput = gauge_field.
        !% This variables defined the directory in which the reference gauge_field field is,
        !% relative to the current folder
        !%End
 
        call parse_variable(global_namespace, 'TransientAbsorptionReference', '.', ref_filename)
        ref_file = io_open(trim(ref_filename)//'/gauge_field', global_namespace, action='read', status='old')
        call io_skip_header(ref_file)
        call spectrum_count_time_steps(global_namespace, ref_file, time_steps_ref, dt_ref)
        if (time_steps_ref < time_steps) then
          message(1) = "The reference calculation does not contain enought time steps"
          call messages_fatal(1)
        end if
 
        if (.not. is_close(dt_ref, dt)) then
          message(1) = "The time step of the reference calculation is different from the current calculation"
          call messages_fatal(1)
        end if
 
      end if
 
      ! Add one as zero is included in the time signal
      time_steps = time_steps + 1
 
      safe_allocate(vecpot(1:time_steps, 1:space%dim*3))
 
      call io_skip_header(in_file)
 
      do ii = 1, time_steps
        read(in_file, *) jj, tt, (vecpot(ii, kk), kk = 1, space%dim*3)
      end do
 
      call io_close(in_file)
 
      !We remove the reference
      if (parse_is_defined(global_namespace, 'TransientAbsorptionReference')) then
        time_steps_ref = time_steps_ref + 1
        safe_allocate(vecpot_ref(1:time_steps_ref, 1:space%dim*3))
        call io_skip_header(ref_file)
        do ii = 1, time_steps_ref
          read(ref_file, *) jj, tt, (vecpot_ref(ii, kk), kk = 1, space%dim*3)
        end do
        call io_close(ref_file)
        do ii = 1, time_steps
          do kk = 1, space%dim*3
            vecpot(ii, kk) = vecpot(ii, kk) - vecpot_ref(ii, kk)
          end do
        end do
      end if
 
      write(message(1), '(a, i7, a)') "Info: Read ", time_steps, " steps from file '"// &
        trim(io_workpath('td.general/gauge_field', global_namespace))//"'"
      call messages_info(1)
 
 
      ! Find out the iteration numbers corresponding to the time limits.
      ! Max time correspond to time_steps-1, as we start from 0.
      call spectrum_fix_time_limits(spectrum, time_steps-1, dt, istart, iend, ntiter)
 
      ! We need to fix istart because the damping starts are (istart-1)*dt, which needs to be spectrum%start_time
      istart = istart + 1
 
      energy_steps = spectrum_nenergy_steps(spectrum)
 
      norm_aext = norm2(aext(1:space%dim))
      if (norm_aext <= tiny(1.0_real64)) then
        message(1) = "GaugeVectorField norm is zero (or underflow-level); cannot normalize dielectric response."
        call messages_fatal(1)
      end if
 
      safe_allocate(eind_real(1:energy_steps, 1:space%dim))
      safe_allocate(eind_imag(1:energy_steps, 1:space%dim))
 
      ! We select dA/dt = E
      call batch_init(vecpotb, 1, 1, space%dim, vecpot(:, space%dim+1:space%dim*2))
      call batch_init(eind_realb, 1, 1, space%dim, eind_real)
      call batch_init(eind_imagb, 1, 1, space%dim, eind_imag)
 
      call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart, iend, spectrum%start_time, dt, vecpotb)
 
      call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
        istart, iend, spectrum%start_time, dt, vecpotb, spectrum%min_energy, &
        spectrum%max_energy, spectrum%energy_step, eind_realb)
 
      call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
        istart, iend, spectrum%start_time, dt, vecpotb, spectrum%min_energy, &
        spectrum%max_energy, spectrum%energy_step, eind_imagb)
 
 
      call vecpotb%end()
      call eind_realb%end()
      call eind_imagb%end()
 
      safe_allocate(invdielectric(1:energy_steps))
      safe_allocate(dielectric(1:energy_steps))
      safe_allocate(chi(1:energy_steps))
 
      do kk = 1, energy_steps
        ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
 
        ! We compute 1/\epsilon(\omega), see Eq. (4) in Bertsch et al. PRB 62, 7998
        ! More precisely, we have \epsilon^{-1}_{\alpha\beta} = \frac{E_{tot,\alpha}}{E_{ext,\beta}}
 
        ! Here we can only assume that the Gauge field is along an axis without off-diagonal terms
        ! Else, we need several calulations
        eind = cmplx(eind_real(kk, 1:space%dim), eind_imag(kk, 1:space%dim), real64)
        invdielectric(kk) = dot_product(aext, aext + eind) /norm_aext**2
 
        dielectric(kk) = m_one / invdielectric(kk)
 
        chi(kk) = (dielectric(kk) - m_one)*latt%rcell_volume/(m_four*m_pi)
      end do
 
      out_file = io_open('td.general/inverse_dielectric_function', global_namespace, action='write')
      write(header, '(7a15)') '#        energy', 'Re', 'Im'
 
      write(out_file,'(a)') trim(header)
      do kk = 1, energy_steps
        ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
        write(out_file, '(e15.6)', advance='no') ww
        write(out_file, '(2e15.6)', advance='no') real(invdielectric(kk), real64), aimag(invdielectric(kk))
        write(out_file, '()')
      end do
      call io_close(out_file)
 
      out_file = io_open('td.general/dielectric_function', global_namespace, action='write')
      write(out_file,'(a)') trim(header)
      do kk = 1, energy_steps
        ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
        write(out_file, '(e15.6)', advance='no') ww
        write(out_file, '(2e15.6)', advance='no') real(dielectric(kk), real64), aimag(dielectric(kk))
        write(out_file, '()')
      end do
      call io_close(out_file)
 
      out_file = io_open('td.general/chi', global_namespace, action='write')
      write(out_file,'(a)') trim(header)
      do kk = 1, energy_steps
        ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
        write(out_file, '(e15.6)', advance='no') ww
        write(out_file, '(2e15.6)', advance='no') real(chi(kk), real64), aimag(chi(kk))
        write(out_file, '()')
      end do
      call io_close(out_file)
 
      safe_deallocate_a(dielectric)
      safe_deallocate_a(invdielectric)
      safe_deallocate_a(chi)
      safe_deallocate_a(vecpot)
      safe_deallocate_a(vecpot_ref)
      safe_deallocate_a(aext)
      safe_deallocate_a(eind)
      safe_deallocate_a(eind_real)
      safe_deallocate_a(eind_imag)
 
    end subroutine dielectric_function_compute
  end subroutine dielectric_function_main
end module dielectric_function_main_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: