casida_spectrum.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2012-2013 D. Strubbe
!!
!! 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 casida_spectrum
  use command_line_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use ions_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 unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use varinfo_oct_m
 
  implicit none
 
  type casida_spectrum_t
    real(real64) :: br, energy_step, min_energy, max_energy
    type(space_t) :: space
    integer :: ispin
  end type casida_spectrum_t
 
  integer :: ierr, idir, jdir, iatom
  type(space_t) :: space
  type(casida_spectrum_t) :: cs
  real(real64), allocatable :: rotation(:,:), rot2(:,:), identity(:,:), coord(:)
  type(block_t) :: blk
  type(ions_t),     pointer :: ions
 
  ! Initialize stuff
  call init_octopus_globals(serial_dummy_comm)
 
  call getopt_init(ierr)
  if (ierr == 0) call getopt_casida_spectrum()
  call getopt_end()
 
  call parser_init()
  call messages_init()
  call io_init()
  call profiling_init(global_namespace)
  call unit_system_init(global_namespace)
  space = space_t(global_namespace)
 
  allocate(rotation(space%dim, space%dim))
  allocate(rot2(space%dim, space%dim))
  allocate(identity(space%dim, space%dim))
  allocate(coord(space%dim))
 
  ! Reads the spin components. This is read here, as well as in states_init.
  call parse_variable(global_namespace, 'SpinComponents', 1, cs%ispin)
  if (.not. varinfo_valid_option('SpinComponents', cs%ispin)) call messages_input_error(global_namespace, 'SpinComponents')
  cs%ispin = min(2, cs%ispin)
 
  !%Variable CasidaSpectrumBroadening
  !%Type float
  !%Default 0.005 Ha
  !%Section Utilities::oct-casida_spectrum
  !%Description
  !% Width of the Lorentzian used to broaden the excitations.
  !%End
  call parse_variable(global_namespace, 'CasidaSpectrumBroadening', 0.005_real64, cs%br, units_inp%energy)
 
  call messages_print_var_value('CasidaSpectrumBroadening', cs%br, unit = units_out%energy)
 
  !%Variable CasidaSpectrumEnergyStep
  !%Type float
  !%Default 0.001 Ha
  !%Section Utilities::oct-casida_spectrum
  !%Description
  !% Sampling rate for the spectrum.
  !%End
  call parse_variable(global_namespace, 'CasidaSpectrumEnergyStep', 0.001_real64, cs%energy_step, units_inp%energy)
 
  call messages_print_var_value('CasidaSpectrumEnergyStep', cs%energy_step, unit = units_out%energy)
 
  !%Variable CasidaSpectrumMinEnergy
  !%Type float
  !%Default 0.0
  !%Section Utilities::oct-casida_spectrum
  !%Description
  !% The broadening is done for energies greater than <tt>CasidaSpectrumMinEnergy</tt>.
  !%End
  call parse_variable(global_namespace, 'CasidaSpectrumMinEnergy', m_zero, cs%min_energy, units_inp%energy)
 
  call messages_print_var_value('CasidaSpectrumMinEnergy', cs%min_energy, unit = units_out%energy)
 
  !%Variable CasidaSpectrumMaxEnergy
  !%Type float
  !%Default 1.0 Ha
  !%Section Utilities::oct-casida_spectrum
  !%Description
  !% The broadening is done for energies smaller than <tt>CasidaSpectrumMaxEnergy</tt>.
  !%End
  call parse_variable(global_namespace, 'CasidaSpectrumMaxEnergy', m_one, cs%max_energy, units_inp%energy)
 
  call messages_print_var_value('CasidaSpectrumMaxEnergy', cs%max_energy, unit = units_out%energy)
 
  identity = m_zero
  do idir = 1, space%dim
    identity(idir, idir) = m_one
  end do
 
  !%Variable CasidaSpectrumRotationMatrix
  !%Type block
  !%Default identity
  !%Section Utilities::oct-casida_spectrum
  !%Description
  !% Supply a rotation matrix to apply to the transition dipoles in generating the spectrum. The rotated atomic structure
  !% will also be output. Size of matrix must be <tt>Dimensions</tt>.
  !%End
 
  if (parse_block(global_namespace, 'CasidaSpectrumRotationMatrix', blk) == 0) then
    do idir = 1, space%dim
      do jdir = 1, space%dim
        call parse_block_float(blk, idir - 1,  jdir - 1, rotation(idir, jdir))
      end do
    end do
    call parse_block_end(blk)
 
    message(1) = "Info: Applying rotation matrix"
    call messages_info(1)
    call output_tensor(rotation, cs%space%dim, unit_one, write_average = .false., namespace=global_namespace)
 
    ! allowing inversions is fine
    rot2(:,:) = abs(matmul(transpose(rotation), rotation))
    if (any(abs(rot2(:,:) - identity(:,:)) > 1e-6_real64)) then
      write(message(1),'(a,es13.6)') "Rotation matrix is not orthogonal. max discrepancy in product = ", &
        maxval(abs(rot2(:,:) - identity(:,:)))
      call messages_warning(1)
    end if
 
    ! apply rotation to geometry
    ions => ions_t(global_namespace)
    do iatom = 1, ions%natoms
      coord = ions%pos(:, iatom)
      ions%pos(:, iatom) = matmul(rotation, coord)
    end do
    call ions%write_xyz(trim(casida_dir)//'rotated')
    safe_deallocate_p(ions)
  else
    rotation(:,:) = identity(:,:)
  end if
 
  call calc_broad(cs, casida_dir, 'eps_diff', .true.)
  call calc_broad(cs, casida_dir, 'petersilka', .false.)
  call calc_broad(cs, casida_dir, 'tamm_dancoff', .false.)
  call calc_broad(cs, casida_dir, 'variational', .false.)
  call calc_broad(cs, casida_dir, 'casida', .false.)
 
  deallocate(rotation)
  deallocate(rot2)
  deallocate(identity)
  deallocate(coord)
 
  call profiling_end(global_namespace)
  call io_end()
  call messages_end()
 
  call parser_end()
  call global_end()
 
contains
 
  ! ---------------------------------------------------------
  subroutine calc_broad(cs, dir, fname, extracols)
    type(casida_spectrum_t), intent(in) :: cs
    character(len=*),        intent(in) :: dir
    character(len=*),        intent(in) :: fname
    logical,                 intent(in) :: extracols
 
    real(real64), allocatable :: spectrum(:,:)
    real(real64) :: omega, energy, re_tm(space%dim), im_tm(space%dim), ff(space%dim+1), tm_sq(space%dim)
    integer :: istep, nsteps, iunit, trash(3), idir, ncols, ios
    character(len=256) :: string
    logical :: is_complex
 
    push_sub(calc_broad)
 
    nsteps = int((cs%max_energy - cs%min_energy) / cs%energy_step)
    safe_allocate(spectrum(1:space%dim+1, 1:nsteps))
    spectrum = m_zero
 
    iunit = io_open(trim(dir)// fname, global_namespace, action='read', status='old', die = .false.)
 
    if (iunit < 0) then
      message(1) = 'Cannot open file "'//trim(dir)//trim(fname)//'".'
      message(2) = 'The '//trim(fname)//' spectrum was not generated.'
      call messages_warning(2)
      return
    end if
 
    ! For Casida, CV(2), Tamm-Dancoff, Petersilka: first column is the index of the excitation
    ! For eps_diff: first two columns are occ and unocc states, then spin if spin-polarized
    ncols = 1
    if (extracols) then
      ncols = ncols + cs%ispin
    end if
 
    read(iunit, *) ! skip header
 
    read(iunit, '(a)') string
    ! complex has this many columns; real lacks the imaginary columns (im_tm)
    read(string, *, iostat = ios) trash(1:ncols), energy, (re_tm(idir), im_tm(idir), idir = 1, space%dim), ff(space%dim+1)
    is_complex = (ios == 0)
    ! put it back to read the data in the loop
    backspace(iunit)
 
    do
      if (is_complex) then
        read(iunit, *, iostat = ios) trash(1:ncols), energy, (re_tm(idir), im_tm(idir), idir = 1, space%dim), ff(space%dim+1)
        im_tm = units_to_atomic(units_out%length, im_tm)
      else
        read(iunit, *, iostat = ios) trash(1:ncols), energy, (re_tm(idir),              idir = 1, space%dim), ff(space%dim+1)
      end if
      re_tm = units_to_atomic(units_out%length, re_tm)
      ! ff, the last column, is a dimensionless number
 
      if (ios < 0) then
        exit ! end of file
      else if (ios > 0) then
        message(1) = "Error parsing file " // trim(fname)
        call messages_fatal(1)
      end if
 
      energy = units_to_atomic(units_out%energy, energy)
      ! transition matrix elements by themselves are dependent on gauge in degenerate subspaces
      ! make into oscillator strengths, as in casida_inc.F90 X(oscillator_strengths), and like the last column
      tm_sq = (matmul(rotation, re_tm))**2
      if (is_complex) then
        tm_sq = tm_sq + (matmul(rotation, im_tm))**2
      end if
      ff(1:space%dim) = m_two * energy * tm_sq
 
      do istep = 1, nsteps
        omega = cs%min_energy + real(istep-1, real64) *cs%energy_step
        spectrum(:, istep) = spectrum(:, istep) + ff*cs%br/((omega-energy)**2 + cs%br**2)/m_pi ! Lorentzian
      end do
    end do
    call io_close(iunit)
 
    ! print spectra
    iunit = io_open(trim(dir)//"/spectrum."//fname, global_namespace, action='write')
 
    write(iunit, '(a2,a12)', advance = 'no') '# ', 'E [' // trim(units_abbrev(units_out%energy)) // ']'
    do idir = 1, space%dim
      write(iunit, '(a14)', advance = 'no') '<' // index2axis(idir) // index2axis(idir) // '>'
    end do
    write(iunit, '(a14)') '<f>'
 
    do istep = 1, nsteps
      write(iunit, '(99es14.6)') units_from_atomic(units_out%energy, cs%min_energy + real(istep - 1, real64)  &
        *cs%energy_step), (units_from_atomic(unit_one/units_out%energy, spectrum(idir, istep)), idir = 1, space%dim+1)
    end do
 
    call io_close(iunit)
 
    safe_deallocate_a(spectrum)
    pop_sub(calc_broad)
  end subroutine calc_broad
 
end program casida_spectrum
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: