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"
 
program dielectric_function
  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
 
  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', die=.false.)
    if (in_file < 0) then
      message(1) = "Cannot open file '"//trim(io_workpath('td.general/gauge_field', global_namespace))//"'"
      call messages_fatal(1)
    end if
 
    ! 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', die=.false.)
      if (ref_file < 0) then
        message(1) = "Cannot open reference file '"// &
          trim(io_workpath(trim(ref_filename)//'/gauge_field', global_namespace))//"'"
        call messages_fatal(1)
      end if
      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))
 
    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 program dielectric_function
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: