infrared.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"
 
program infrared
  use batch_oct_m
  use command_line_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  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 :: iunit, ierr, ii, jj, iter, read_iter
  real(real64), allocatable :: dipole(:,:)
  complex(real64), allocatable :: ftdipole(:,:)
  type(space_t)     :: space
  real(real64) :: ww, irtotal, dt
  integer :: ifreq, idir, istart, iend, time_steps, energy_steps, ntiter
  type(spectrum_t)  :: spectrum
 
  ! Initialize stuff
  call init_octopus_globals(serial_dummy_comm)
 
  call getopt_init(ierr)
 
  call getopt_end()
 
  call parser_init()
 
  call messages_init()
 
  call io_init()
  call profiling_init(global_namespace)
 
  call unit_system_init(global_namespace)
 
  call spectrum_init(spectrum, global_namespace, &
    default_max_energy  = units_to_atomic(unit_invcm, 10000.0_real64))
 
  space = space_t(global_namespace)
 
  call read_dipole(dipole)
 
  safe_allocate(ftdipole(1:energy_steps, 1:space%dim))
  call fourier(dipole, ftdipole)
 
  !and print the spectrum
  iunit = io_open('td.general/infrared', global_namespace, action='write')
 
100 FORMAT(100('#'))
 
  write(unit = iunit, iostat = ierr, fmt = 100)
  write(unit = iunit, iostat = ierr, fmt = '(8a)')  '# HEADER'
  write(unit = iunit, iostat = ierr, fmt = '(a25,a1,a1,a10)') &
    '# all absorptions are in ',units_out%length%abbrev,'*',units_out%time%abbrev
  write(unit = iunit, iostat = ierr, fmt = '(a1)') '#'
  write(unit = iunit, iostat = ierr, fmt = '(a19,41x,a17)') '#        Energy    ', 'absorption'
  write(unit = iunit, iostat = ierr, fmt = '(a15,13x,a5,15x,a7,13x,a7,13x,a7)') &
    '#        [1/cm]','total','FT(<x>)','FT(<y>)','FT(<z>)'
  write(unit = iunit, iostat = ierr, fmt = 0100)
 
  do ifreq = 1, energy_steps
    ww = (ifreq-1)*spectrum%energy_step + spectrum%min_energy
    irtotal = norm2(abs(ftdipole(ifreq, 1:3)))
    write(unit = iunit, iostat = ierr, fmt = '(5e20.10)') &
      units_from_atomic(unit_invcm, ww), units_from_atomic(units_out%length*units_out%time, ww*irtotal), &
      (units_from_atomic(units_out%length*units_out%time,  ww*abs(ftdipole(ifreq, idir))), idir = 1, 3)
  end do
  call io_close(iunit)
 
  safe_deallocate_a(dipole)
  safe_deallocate_a(ftdipole)
 
  call profiling_end(global_namespace)
  call io_end()
  call messages_end()
 
  call parser_end()
  call global_end()
 
contains
 
  subroutine read_dipole(dipole)
    real(real64), allocatable,  intent(inout)   :: dipole(:, :)
    character(len=100) :: ioerrmsg
    real(real64) :: charge, time
 
    push_sub(read_dipole)
 
    ! Opens the coordinates files.
    iunit = io_open('td.general/multipoles', global_namespace, action='read', status='old')
    call io_skip_header(iunit)
    call spectrum_count_time_steps(global_namespace, iunit, time_steps, dt)
    time_steps = time_steps + 1
 
    ! Find out the iteration numbers corresponding to the time limits.
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
    istart = max(1, istart)
    energy_steps = spectrum_nenergy_steps(spectrum)
 
    safe_allocate(dipole(1:time_steps, 1:space%dim))
 
    call io_skip_header(iunit)
 
    do iter = 1, time_steps
      read(unit = iunit, iostat = ierr, iomsg=ioerrmsg, fmt = *) read_iter, time, &
        charge, (dipole(iter, idir), idir=1, space%dim)
 
      !dipole moment has unit of charge*length, charge has the same unit in both systems
      do ii = 1, space%dim
        dipole(iter, ii) = units_to_atomic(units_out%length, dipole(iter, ii))
      end do
 
      if (ierr /= 0) then
        exit
      end if
    end do
 
    if (ierr == 0) then
      write (message(1), '(a)') "Read dipole moment from '"// &
        trim(io_workpath('td.general/multipoles', global_namespace))//"'."
      call messages_info(1)
    else
      message(1) = "Error reading file '"//trim(io_workpath('td.general/multipoles', global_namespace))//"'"
      message(2) = ioerrmsg
      call messages_fatal(2)
    end if
 
    call io_close(iunit)
 
    pop_sub(read_dipole)
  end subroutine read_dipole
 
  ! -------------------------------------------------
  subroutine fourier(fi, ftfi)
    implicit none
    real(real64), contiguous, intent(inout)  :: fi(:,:)
    complex(real64), intent(out)    :: ftfi(:,:)
 
    real(real64) :: av
    integer :: ifreq
    type(batch_t)     :: dipoleb, ftrealb, ftimagb
    real(real64), allocatable :: ftreal(:,:), ftimag(:,:)
 
    push_sub(fourier)
 
    call batch_init(dipoleb, 1, 1, space%dim, fi)
 
    !apply an envelope
    call spectrum_signal_damp(spectrum_damp_sin, spectrum%damp_factor, istart, iend, spectrum%start_time, dt, dipoleb)
 
    !remove the DC component
    do idir = 1, space%dim
      av = sum(fi(istart:iend, idir)) / (iend - istart + 1)
      !$omp parallel do
      do jj = istart, iend
        fi(jj, idir) = fi(jj, idir) - av
      end do
    end do
 
    write (message(1), '(a)') "Taking the Fourier transform."
    call messages_info(1)
 
    safe_allocate(ftreal(1:energy_steps, 1:space%dim))
    safe_allocate(ftimag(1:energy_steps, 1:space%dim))
    call batch_init(ftrealb, 1, 1, space%dim, ftreal)
    call batch_init(ftimagb, 1, 1, space%dim, ftimag)
 
    !now calculate the FT
    call spectrum_fourier_transform(spectrum_fourier, spectrum_transform_cos, m_zero, &
      istart, iend, spectrum%start_time, dt, dipoleb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftrealb)
 
    call spectrum_fourier_transform(spectrum_fourier, spectrum_transform_sin, m_zero, &
      istart, iend, spectrum%start_time, dt, dipoleb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftimagb)
 
    do idir = 1, space%dim
      !$omp parallel do
      do ifreq = 1, energy_steps
        ftfi(ifreq, idir) = cmplx(ftreal(ifreq, idir), ftimag(ifreq, idir), real64)
      end do
    end do
 
    call dipoleb%end()
    call ftrealb%end()
    call ftimagb%end()
    safe_deallocate_a(ftreal)
    safe_deallocate_a(ftimag)
 
    write (message(1), '(a)') "Done."
    call messages_info(1)
 
    pop_sub(fourier)
  end subroutine fourier
 
end program infrared
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: