oscillator_strength.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"
 
module oscillator_strength_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use kick_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use messages_oct_m
  use minimizer_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use spectrum_oct_m
  use string_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  integer, parameter ::   &
    SINE_TRANSFORM   = 1, &
    cosine_transform = 2
 
  type local_operator_t
    integer :: n_multipoles
    integer, allocatable :: l(:), m(:)
    real(real64),   allocatable :: weight(:)
  end type local_operator_t
 
  integer             :: observable(2)
  type(unit_system_t) :: units
  real(real64), allocatable  :: ot(:)
  type(kick_t)        :: kick
  integer             :: time_steps
  real(real64)        :: total_time
  integer             :: mode
  real(real64)        :: dt
 
contains
 
  ! ---------------------------------------------------------
  subroutine ft(omega, power)
    real(real64), intent(in)   :: omega
    real(real64), intent(out)  :: power
 
    integer :: j
    real(real64) :: x
    power = m_zero
 
    push_sub(ft)
 
    select case (mode)
 
    case (sine_transform)
 
      do j = 0, time_steps
        x = sin(omega*j*dt)
        power = power + x*ot(j)
      end do
      power = power*dt / (dt*time_steps)
 
    case (cosine_transform)
 
      do j = 0, time_steps
        x = cos(omega*j*dt)
        power = power + x*ot(j)
      end do
      power = power*dt / (dt*time_steps)
 
    end select
 
    pop_sub(ft)
  end subroutine ft
 
 
  ! ---------------------------------------------------------
  subroutine ft2(omega, power)
    real(real64), intent(in)   :: omega
    real(real64), intent(out)  :: power
 
    integer :: j
    real(real64) :: x
    power = m_zero
 
    push_sub(ft2)
 
    select case (mode)
 
    case (sine_transform)
 
      do j = 0, time_steps
        x = sin(omega*j*dt)
        power = power + x*ot(j)
      end do
      ! The function should be *minus* the sine Fourier transform, since this
      ! is the function to be minimized.
      power = - (power*dt / (dt*time_steps))**2
 
    case (cosine_transform)
 
      do j = 0, time_steps
        x = cos(omega*j*dt)
        power = power + x*ot(j)
      end do
      power = - (power*dt / (dt*time_steps))**2
 
    end select
 
    pop_sub(ft2)
  end subroutine ft2
 
 
  ! ---------------------------------------------------------
  subroutine local_operator_copy(o, i)
    type(local_operator_t), intent(inout) :: o
    type(local_operator_t), intent(inout) :: i
 
    integer :: j
 
    push_sub(local_operator_copy)
 
    o%n_multipoles = i%n_multipoles
    safe_allocate(     o%l(1:o%n_multipoles))
    safe_allocate(     o%m(1:o%n_multipoles))
    safe_allocate(o%weight(1:o%n_multipoles))
 
    do j = 1, o%n_multipoles
      o%l(j) = i%l(j)
      o%m(j) = i%m(j)
      o%weight(j) = i%weight(j)
    end do
 
    pop_sub(local_operator_copy)
  end subroutine local_operator_copy
 
  ! ---------------------------------------------------------
  subroutine read_resonances_file(order, ffile, namespace, search_interval, final_time, nfrequencies)
    integer, intent(inout)       :: order
    character(len=*), intent(in) :: ffile
    type(namespace_t), intent(in):: namespace
    real(real64),   intent(inout)       :: search_interval
    real(real64),   intent(in)          :: final_time
    integer, intent(in)          :: nfrequencies
 
    real(real64) :: dummy, leftbound, rightbound, w, power, dw
    integer :: iunit, nresonances, ios, i, j, k, npairs, nspin, order_in_file, nw_subtracted, ierr
    logical :: file_exists
    real(real64), allocatable :: wij(:), omega(:), c0i(:)
 
    push_sub(read_resonances_file)
 
    if (order /= 2) then
      write(message(1),'(a)') 'The run mode #3 is only compatible with the analysis of the'
      write(message(2),'(a)') 'second-order response.'
      call messages_fatal(2)
    end if
 
    ! First, let us check that the file "ot" exists.
    inquire(file="ot", exist  = file_exists)
    if (.not. file_exists) then
      write(message(1),'(a)') "Could not find 'ot' file."
      call messages_fatal(1)
    end if
 
    ! Now, we should find out which units the file "ot" has.
    call unit_system_from_file(units, "ot", namespace, ierr)
    if (ierr /= 0) then
      write(message(1),'(a)') "Could not retrieve units in the 'ot' file."
      call messages_fatal(1)
    end if
 
    mode = cosine_transform
 
    iunit = io_open(trim(ffile), namespace, action='read', status='old', die=.false.)
    if (iunit == 0) then
      write(message(1),'(a)') 'Could not open '//trim(ffile)//' file.'
      call messages_fatal(1)
    end if
 
    call io_skip_header(iunit)
    ! Count the number of resonances
    nresonances = 0
    do
      read(iunit, *, iostat = ios) dummy, dummy
      if (ios /= 0) exit
      nresonances = nresonances + 1
    end do
 
    npairs = (nresonances*(nresonances-1))/2
 
    safe_allocate(omega(1:nresonances))
    safe_allocate(  c0i(1:nresonances))
    safe_allocate(wij(1:npairs))
 
    call io_skip_header(iunit)
    do i = 1, nresonances
      read(iunit, *) omega(i), c0i(i)
    end do
 
    k = 1
    do i = 1, nresonances
      do j = i + 1, nresonances
        wij(k) = omega(j) - omega(i)
        k = k + 1
      end do
    end do
 
    if (search_interval > m_zero) then
      search_interval = units_to_atomic(units%energy, search_interval)
    else
      search_interval = m_half
    end if
 
    call read_ot(namespace, nspin, order_in_file, nw_subtracted)
 
    if (order_in_file /= order) then
      write(message(1), '(a)') 'The ot file should contain the second-order response in this run mode.'
      call messages_fatal(1)
    end if
 
    if (final_time > m_zero) then
      total_time = units_to_atomic(units%time, final_time)
      if (total_time > dt*time_steps) then
        total_time = dt*time_steps
        write(message(1), '(a)')        'The requested total time to process is larger than the time available in the input file.'
        write(message(2), '(a,f8.4,a)') 'The time has been adjusted to ', &
          units_from_atomic(units%time, total_time), units_abbrev(units%time)
        call messages_warning(2, namespace=namespace)
      end if
      time_steps = int(total_time / dt)
      total_time = time_steps * dt
    else
      total_time = dt*time_steps
    end if
 
    dw = (rightbound-leftbound) / (nfrequencies - 1)
 
    ! First, subtract zero resonance...
    w = m_zero
    call resonance_second_order(w, power, nw_subtracted, leftbound, rightbound, m_zero, m_zero)
    call modify_ot(time_steps, dt, order, ot, w, power)
    nw_subtracted = nw_subtracted + 1
 
    ! Then, get all the others...
    k = 1
    do i = 1, nresonances
      do j = i + 1, nresonances
        leftbound = wij(k) - search_interval
        rightbound = wij(k) + search_interval
        call find_resonance(wij(k), leftbound, rightbound, nfrequencies)
        call resonance_second_order(wij(k), power, nw_subtracted, leftbound, rightbound, c0i(i), c0i(j))
        call modify_ot(time_steps, dt, order, ot, wij(k), power)
        nw_subtracted = nw_subtracted + 1
        k = k + 1
      end do
    end do
 
    safe_deallocate_a(wij)
    safe_deallocate_a(c0i)
    safe_deallocate_a(omega)
    call io_close(iunit)
 
    pop_sub(read_resonances_file)
  end subroutine read_resonances_file
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine analyze_signal(order, omega, search_interval, final_time, nresonances, nfrequencies, &
    damping, namespace)
    integer, intent(inout) :: order
    real(real64),   intent(inout) :: omega
    real(real64),   intent(inout) :: search_interval
    real(real64),   intent(inout) :: final_time
    integer, intent(inout) :: nresonances
    integer, intent(inout) :: nfrequencies
    real(real64),   intent(in)    :: damping
    type(namespace_t), intent(in) :: namespace
 
    real(real64) :: leftbound, rightbound, dw, power
    real(real64), allocatable :: w(:), c0I2(:)
    integer :: nspin, i, ierr, order_in_file, nw_subtracted
    logical :: file_exists
 
    push_sub(analyze_signal)
 
    ! First, let us check that the file "ot" exists.
    inquire(file="ot", exist  = file_exists)
    if (.not. file_exists) then
      write(message(1),'(a)') "Could not find 'ot' file."
      call messages_fatal(1)
    end if
 
    ! Now, we should find out which units the file "ot" has.
    call unit_system_from_file(units, "ot", namespace, ierr)
    if (ierr /= 0) then
      write(message(1),'(a)') "Could not retrieve units in the 'ot' file."
      call messages_fatal(1)
    end if
 
    if (omega > m_zero) then
      omega = units_to_atomic(units%energy, omega)
    else
      omega = m_half
    end if
 
    if (search_interval > m_zero) then
      search_interval = units_to_atomic(units%energy, search_interval)
    else
      search_interval = m_half
    end if
 
    leftbound = omega - search_interval
    rightbound = omega + search_interval
 
    call read_ot(namespace, nspin, order_in_file, nw_subtracted)
 
    if (order_in_file /= order) then
      write(message(1), '(a)') 'Internal error in analyze_signal'
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (mod(order, 2) == 1) then
      mode = sine_transform
    else
      mode = cosine_transform
    end if
 
    if (final_time > m_zero) then
      total_time = units_to_atomic(units%time, final_time)
      if (total_time > dt*time_steps) then
        total_time = dt*time_steps
        write(message(1), '(a)')        'The requested total time to process is larger than the time available in the input file.'
        write(message(2), '(a,f8.4,a)') 'The time has been adjusted to ', &
          units_from_atomic(units%time, total_time), units_abbrev(units%time)
        call messages_warning(2, namespace=namespace)
      end if
      time_steps = int(total_time / dt)
      total_time = time_steps * dt
    else
      total_time = dt*time_steps
    end if
 
    dw = (rightbound-leftbound) / (nfrequencies - 1)
 
    safe_allocate(   w(1:nresonances))
    safe_allocate(c0i2(1:nresonances))
    w = m_zero
    c0i2 = m_zero
 
    i = 1
    do
      if (nw_subtracted >= nresonances) exit
 
      if (mode == cosine_transform .and. nw_subtracted == 0) then
        omega = m_zero
      else
        call find_resonance(omega, leftbound, rightbound, nfrequencies)
      end if
 
      select case (order)
      case (1)
        call resonance_first_order(omega, power, nw_subtracted, dw, leftbound, rightbound)
      case (2)
        call resonance_second_order(omega, power, nw_subtracted, leftbound, rightbound, m_zero, m_zero)
      end select
 
      w(i) = omega
      c0i2(i) = power
 
      call modify_ot(time_steps, dt, order, ot, omega, power)
 
      nw_subtracted = nw_subtracted + 1
      i = i + 1
    end do
 
    select case (order)
    case (1)
      call write_polarizability(namespace, nfrequencies, nresonances, dw, w, c0i2, damping)
    end select
 
    safe_deallocate_a(ot)
    safe_deallocate_a(w)
    safe_deallocate_a(c0i2)
 
    pop_sub(analyze_signal)
  end subroutine analyze_signal
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine write_polarizability(namespace, nfrequencies, nresonances, dw, w, c0I2, gamma)
    type(namespace_t), intent(in) :: namespace
    integer, intent(in) :: nfrequencies, nresonances
    real(real64), intent(in)   :: dw
    real(real64), intent(in)   :: w(nresonances), c0I2(nresonances)
    real(real64), intent(in)   :: gamma
 
    integer :: iunit, i, j
    real(real64) :: e
    complex(real64) :: pol
 
    push_sub(write_polarizability)
 
    iunit = io_open('polarizability', namespace, action = 'write', status='replace', die=.false.)
    write(iunit, '(a)') '# Polarizability file. Generated using the SOS formula with the following data:'
    write(iunit, '(a)') '#'
 
    do i = 1, nresonances
      write(iunit, '(a1,3e20.12)') '#', w(i), sqrt(abs(c0i2(i))), c0i2(i)
    end do
 
    write(iunit, '(a10,f12.6)') '# Gamma = ', gamma
    write(iunit, '(a)')         '#'
 
    do i = 1, nfrequencies
      e = (i - 1) * dw
      pol = m_z0
      do j = 1, nresonances
        pol = pol + c0i2(j) * (m_one / (w(j) - e - m_zi * gamma) + m_one/(w(j) + e + m_zi * gamma))
      end do
      write(iunit, '(3e20.12)') e, pol
    end do
 
    call io_close(iunit)
 
    pop_sub(write_polarizability)
  end subroutine write_polarizability
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  ! \todo This subroutine should be simplified.
  subroutine find_resonance(omega, leftbound, rightbound, nfrequencies)
    real(real64), intent(inout) :: omega
    real(real64), intent(in)    :: leftbound, rightbound
    integer, intent(in)  :: nfrequencies
 
    integer :: i, ierr
    real(real64) :: dw, w, aw, min_aw, min_w, omega_orig, min_w1, min_w2
    real(real64), allocatable :: warray(:), tarray(:)
 
    push_sub(find_resonance)
 
    safe_allocate(warray(1:nfrequencies))
    safe_allocate(tarray(1:nfrequencies))
 
    warray = m_zero
    tarray = m_zero
    dw = (rightbound-leftbound) / (nfrequencies - 1)
    do i = 1, nfrequencies
      w = leftbound + (i-1)*dw
      warray(i) = w
      call ft2(w, aw)
      tarray(i) = aw
    end do
 
    min_w = omega
    min_aw = m_zero
    do i = 1, nfrequencies
      w = leftbound + (i-1)*dw
      if (tarray(i) < min_aw) then
        min_aw = tarray(i)
        min_w = w
      end if
    end do
 
    omega_orig = omega
    omega = min_w
    min_w1 = min_w - 2*dw
    min_w2 = min_w + 2*dw
    call loct_1dminimize(min_w1, min_w2, omega, ft2, ierr)
 
    if (ierr /= 0) then
      write(message(1),'(a)') 'Could not find a maximum.'
      write(message(2),'(a)')
      write(message(3), '(a,f12.8,a,f12.8,a)') '   Search interval = [', &
        units_from_atomic(units%energy, leftbound), ',', units_from_atomic(units%energy, rightbound), ']'
      write(message(4), '(a,f12.4,a)')         '   Search discretization = ', &
        units_from_atomic(units%energy, dw), ' '//trim(units_abbrev(units%energy))
      call messages_fatal(4)
    end if
 
    safe_deallocate_a(warray)
    safe_deallocate_a(tarray)
 
    pop_sub(find_resonance)
  end subroutine find_resonance
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine resonance_first_order(omega, power, nw_subtracted, dw, leftbound, rightbound)
    real(real64), intent(in)               :: omega
    real(real64), intent(out)              :: power
    integer, intent(in)             :: nw_subtracted
    real(real64), intent(in)               :: dw, leftbound, rightbound
 
    push_sub(resonance_first_order)
 
    call ft(omega, power)
 
    select case (mode)
    case (sine_transform)
      power = power / (m_one - sin(m_two*omega*total_time)/(m_two*omega*total_time))
    case (cosine_transform)
      call messages_not_implemented("resonance first order cosine transform")
    end select
 
    write(message(1), '(a)')                 '******************************************************************'
    write(message(2), '(a,i3)')              'Resonance #', nw_subtracted + 1
    write(message(3), '(a,f12.8,a,f12.8,a)') 'omega    = ', units_from_atomic(units_out%energy, omega), &
      ' '//trim(units_abbrev(units_out%energy))//' = ',  omega, ' Ha'
    write(message(4), '(a,f12.8,a,f12.8,a)') 'C(omega) = ', units_from_atomic(units_out%length**2, power), &
      ' '//trim(units_abbrev(units_out%length**2))//' =', power, ' b^2'
    write(message(5), '(a,f12.8,a,f12.8,a)') '<0|P|I>  = ', units_from_atomic(units_out%length, sqrt(abs(power))), &
      ' '//trim(units_abbrev(units_out%length))//' = ', sqrt(abs(power)),' b'
    write(message(6), '(a,f12.8)')           'f[O->I]  = ', m_two*omega*power
    write(message(7), '(a)')
    write(message(8), '(a,f12.8,a,f12.8,a)') '   Search interval = [', units_from_atomic(units_out%energy, leftbound), ',', &
      units_from_atomic(units_out%energy, rightbound), ']'
    write(message(9), '(a,f12.4,a)')         '   Search discretization = ', units_from_atomic(units_out%energy, dw), &
      ' '//trim(units_abbrev(units_out%energy))
    write(message(10), '(a)')                '******************************************************************'
    write(message(11), '(a)')
    call messages_info(11)
 
    pop_sub(resonance_first_order)
 
  end subroutine resonance_first_order
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine resonance_second_order(omega, power, nw_subtracted, leftbound, rightbound, c01, c02)
    real(real64), intent(in)               :: omega
    real(real64), intent(out)              :: power
    integer, intent(in)             :: nw_subtracted
    real(real64), intent(in)               :: leftbound, rightbound
    real(real64), intent(in)               :: c01, c02
 
    push_sub(resonance_second_order)
 
    call ft(omega, power)
    select case (mode)
    case (sine_transform)
      power = power / (m_one - sin(m_two*omega*total_time)/(m_two*omega*total_time))
    case (cosine_transform)
      ! WARNING: there is some difference between the omega=0 case and the rest.
      if (abs(omega) > m_epsilon) then
        power = power / (m_one + sin(m_two*omega*total_time)/(m_two*omega*total_time))
      else
        power = power / m_two
      end if
    end select
 
    write(message(1), '(a)')                 '******************************************************************'
    write(message(2), '(a,i3)')              'Resonance #', nw_subtracted + 1
    write(message(3), '(a,f12.8,a,f12.8,a)') 'omega    = ', units_from_atomic(units_out%energy, omega), &
      ' '//trim(units_abbrev(units_out%energy))//' = ', omega, ' Ha'
    write(message(4), '(a,f12.8,a,f12.8,a)') 'C(omega) = ', units_from_atomic(units_out%length**3, power), &
      ' '//trim(units_abbrev(units_out%length**3))//' = ', power, ' b^3'
    call messages_info(4)
 
    if (abs(c01*c02) > m_epsilon) then
      write(message(1), '(a,f12.8)')         '    C(omega)/(C0i*C0j) = ', power / (c01 * c02)
      call messages_info(1)
    end if
 
    write(message(1), '(a)')
    write(message(2), '(a,f12.8,a,f12.8,a)') '   Search interval = [', units_from_atomic(units_out%energy, leftbound), ',', &
      units_from_atomic(units_out%energy, rightbound), ']'
    write(message(3), '(a)')                 '******************************************************************'
    write(message(4), '(a)')
    call messages_info(4)
 
    pop_sub(resonance_second_order)
  end subroutine resonance_second_order
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine generate_signal(namespace, order, observable)
    type(namespace_t), intent(in) :: namespace
    integer, intent(in) :: order
    integer, intent(in) :: observable(2)
 
    logical :: file_exists
    integer :: i, j, nspin, time_steps, lmax, nfiles, k, add_lm, l, m, max_add_lm
    integer, allocatable :: iunit(:)
    real(real64) :: dt, lambda, dump, o0
    type(unit_t) :: mp_unit
    real(real64), allocatable :: q(:), mu(:), qq(:, :), c(:)
    character(len=20) :: filename
    type(kick_t) :: kick
    type(unit_system_t) :: units
    real(real64), allocatable :: multipole(:, :, :), ot(:), dipole(:, :)
    type(local_operator_t) :: kick_operator
    type(local_operator_t) :: obs
 
    push_sub(generate_signal)
 
    ! Find out how many files do we have
    nfiles = 0
    do
      write(filename,'(a11,i1)') 'multipoles.', nfiles+1
      inquire(file=trim(filename), exist  = file_exists)
      if (.not. file_exists) exit
      nfiles = nfiles + 1
    end do
 
    ! WARNING: Check that order is smaller or equal to nfiles
    if (nfiles == 0) then
      write(message(1),'(a)') 'No multipoles.x file was found'
      call messages_fatal(1, namespace=namespace)
    end if
    if (order > nfiles) then
      write(message(1),'(a)') 'The order that you ask for is higher than the number'
      write(message(2),'(a)') 'of multipoles.x file that you supply.'
      call messages_fatal(2, namespace=namespace)
    end if
 
    ! Open the files.
    safe_allocate(iunit(1:nfiles))
    do j = 1, nfiles
      write(filename,'(a11,i1)') 'multipoles.', j
      iunit(j) = io_open(trim(filename), namespace, action='read', status='old', die=.false.)
    end do
 
    safe_allocate( q(1:nfiles))
    safe_allocate(mu(1:nfiles))
    safe_allocate(qq(1:nfiles, 1:nfiles))
    safe_allocate( c(1:nfiles))
 
    c        = m_zero
    c(order) = m_one
 
    ! Get the basic info from the first file
    call spectrum_mult_info(namespace, iunit(1), nspin, kick, time_steps, dt, units, lmax=lmax)
 
    ! Sets the kick operator...
    if (kick%n_multipoles > 0) then
      kick_operator%n_multipoles = kick%n_multipoles
      safe_allocate(     kick_operator%l(1:kick_operator%n_multipoles))
      safe_allocate(     kick_operator%m(1:kick_operator%n_multipoles))
      safe_allocate(kick_operator%weight(1:kick_operator%n_multipoles))
      do i = 1, kick_operator%n_multipoles
        kick_operator%l(i) = kick%l(i)
        kick_operator%m(i) = kick%m(i)
        kick_operator%weight(i) = kick%weight(i)
      end do
    else
      kick_operator%n_multipoles = 3
      safe_allocate(     kick_operator%l(1:kick_operator%n_multipoles))
      safe_allocate(     kick_operator%m(1:kick_operator%n_multipoles))
      safe_allocate(kick_operator%weight(1:kick_operator%n_multipoles))
      kick_operator%l(1:3) = 1
      kick_operator%m(1) = -1
      kick_operator%m(2) =  0
      kick_operator%m(3) =  1
      ! WARNING: not sure if m = -1 => y, and m = 1 => x. What is the convention?
      kick_operator%weight(1) = -sqrt((m_four*m_pi)/m_three) * kick%pol(2, kick%pol_dir)
      kick_operator%weight(2) =  sqrt((m_four*m_pi)/m_three) * kick%pol(3, kick%pol_dir)
      kick_operator%weight(3) = -sqrt((m_four*m_pi)/m_three) * kick%pol(1, kick%pol_dir)
    end if
 
    ! Sets the observation operator
    select case (observable(1))
    case (-1)
      ! This means that the "observation operator" should be equal
      ! to the "perturbation operator", i.e., the kick.
      call local_operator_copy(obs, kick_operator)
    case (0)
      ! This means that the observable is the dipole operator; observable(2) determines
      ! if it is x, y or z.
      obs%n_multipoles = 1
      safe_allocate(obs%l(1:1))
      safe_allocate(obs%m(1:1))
      safe_allocate(obs%weight(1:1))
      obs%l(1) = 1
      select case (observable(2))
      case (1)
        obs%m(1) = -1
        obs%weight(1) = -sqrt((m_four*m_pi)/m_three)
      case (2)
        obs%m(1) =  1
        obs%weight(1) = sqrt((m_four*m_pi)/m_three)
      case (3)
        obs%m(1) =  0
        obs%weight(1) = -sqrt((m_four*m_pi)/m_three)
      end select
    case default
      ! This means that the observation operator is (l,m) = (observable(1), observable(2))
      obs%n_multipoles = 1
      safe_allocate(obs%l(1:1))
      safe_allocate(obs%m(1:1))
      safe_allocate(obs%weight(1:1))
      obs%weight(1) = m_one
      obs%l(1) = observable(1)
      obs%m(1) = observable(2)
    end select
 
    lambda = abs(kick%delta_strength)
    q(1) = kick%delta_strength / lambda
 
    do j = 2, nfiles
      call spectrum_mult_info(namespace, iunit(j), nspin, kick, time_steps, dt, units, lmax=lmax)
      q(j) = kick%delta_strength / lambda
    end do
 
    do i = 1, nfiles
      do j = 1, nfiles
        qq(i, j) = q(j)**i
      end do
    end do
 
    call lalg_inverse(nfiles, qq, 'dir')
 
    mu = matmul(qq, c)
 
    if (kick%n_multipoles > 0) then
      lmax = maxval(kick%l(1:obs%n_multipoles))
      max_add_lm = (lmax+1)**2-1
      safe_allocate(multipole(1:max_add_lm, 0:time_steps, 1:nspin))
      ! The units have nothing to do with the perturbing kick??
      mp_unit = units%length**kick%l(1)
    else
      max_add_lm = 3
      safe_allocate(multipole(1:3, 0:time_steps, 1:nspin))
      mp_unit = units%length
    end if
    safe_allocate(ot(0:time_steps))
    multipole = m_zero
    ot = m_zero
    o0 = m_zero
 
    safe_allocate(dipole(1:3, 1:nspin))
 
    do j = 1, nfiles
      call io_skip_header(iunit(j))
 
      do i = 0, time_steps
        select case (nspin)
        case (1)
          read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm)
        case (2)
          read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm), &
            dump, (multipole(add_lm, i, 2), add_lm = 1, max_add_lm)
        case (4)
          read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm), &
            dump, (multipole(add_lm, i, 2), add_lm = 1, max_add_lm), &
            dump, (multipole(add_lm, i, 3), add_lm = 1, max_add_lm), &
            dump, (multipole(add_lm, i, 4), add_lm = 1, max_add_lm)
        end select
        multipole(1:max_add_lm, i, :) = units_to_atomic(mp_unit, multipole(1:max_add_lm, i, :))
 
        ! The dipole is treated differently in the multipoles file: first of all,
        ! the program should have printed *minus* the dipole operator.
        dipole(1:3, 1:nspin) = - multipole(1:3, i, 1:nspin)
        ! And then it contains the "cartesian" dipole, opposed to the spherical dipole:
        multipole(1, i, 1:nspin) = -sqrt(m_three/(m_four*m_pi)) * dipole(2, 1:nspin)
        multipole(2, i, 1:nspin) =  sqrt(m_three/(m_four*m_pi)) * dipole(3, 1:nspin)
        multipole(3, i, 1:nspin) = -sqrt(m_three/(m_four*m_pi)) * dipole(1, 1:nspin)
 
        dump = m_zero
        do k = 1, obs%n_multipoles
          add_lm = 1
          l = 1
          lcycle: do
            do m = -l, l
              if (l == obs%l(k) .and. m == obs%m(k)) exit lcycle
              add_lm = add_lm + 1
            end do
            l = l + 1
          end do lcycle
          ! Warning: it should not add the nspin components?
          dump = dump + obs%weight(k) * sum(multipole(add_lm, i, 1:nspin))
        end do
 
        if (i == 0) o0 = dump
        ot(i) = ot(i) + mu(j)*(dump - o0)
      end do
 
    end do
 
    ot = ot / lambda**order
 
    call write_ot(namespace, nspin, time_steps, dt, kick, order, ot(0:time_steps), observable)
 
    ! Close files and exit.
    do j = 1, nfiles
      call io_close(iunit(j))
    end do
    safe_deallocate_a(iunit)
    safe_deallocate_a(q)
    safe_deallocate_a(mu)
    safe_deallocate_a(qq)
    safe_deallocate_a(c)
    safe_deallocate_a(ot)
    safe_deallocate_a(multipole)
    safe_deallocate_a(dipole)
 
    pop_sub(generate_signal)
  end subroutine generate_signal
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine modify_ot(time_steps, dt, order, ot, omega, power)
    integer,             intent(in)    :: time_steps
    real(real64),        intent(in)    :: dt
    integer,             intent(in)    :: order
    real(real64),        intent(inout) :: ot(0:time_steps)
    real(real64),        intent(in)    :: omega
    real(real64),        intent(in)    :: power
 
    integer :: i
 
    push_sub(modify_ot)
 
    select case (mod(order, 2))
    case (1)
      do i = 0, time_steps
        ot(i) = ot(i) - m_two * power * sin(omega*i*dt)
      end do
    case (0)
      if (abs(omega) <= m_epsilon) then
        do i = 0, time_steps
          ot(i) = ot(i) - power * cos(omega*i*dt)
        end do
      else
        do i = 0, time_steps
          ot(i) = ot(i) - m_two * power * cos(omega*i*dt)
        end do
      end if
    end select
 
    pop_sub(modify_ot)
  end subroutine modify_ot
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine write_ot(namespace, nspin, time_steps, dt, kick, order, ot, observable)
    type(namespace_t),   intent(in) :: namespace
    integer,             intent(in) :: nspin, time_steps
    real(real64),        intent(in) :: dt
    type(kick_t),        intent(in) :: kick
    integer,             intent(in) :: order
    real(real64),        intent(in) :: ot(0:time_steps)
    integer,             intent(in) :: observable(2)
 
    integer :: iunit, i
    character(len=20) :: header_string
    type(unit_t) :: ot_unit
 
    push_sub(write_ot)
 
    iunit = io_open('ot', namespace, action='write', status='replace')
 
    write(iunit, '(a15,i2)')      '# nspin        ', nspin
    write(iunit, '(a15,i2)')      '# Order        ', order
    write(iunit, '(a28,i2)')      '# Frequencies subtracted = ', 0
    select case (observable(1))
    case (-1)
      write(iunit,'(a)') '# Observable operator = kick operator'
      if (kick%n_multipoles > 0) then
        ot_unit = units_out%length**kick%l(1)
      else
        ot_unit = units_out%length
      end if
    case (0)
      select case (observable(2))
      case (1)
        write(iunit,'(a)') '# O = x'
      case (2)
        write(iunit,'(a)') '# O = y'
      case (3)
        write(iunit,'(a)') '# O = z'
      end select
      ot_unit = units_out%length
    case default
      ot_unit = units_out%length**observable(1)
      write(iunit, '(a12,i1,a1,i2,a1)') '# (l, m) = (', observable(1),',',observable(2),')'
    end select
    call kick_write(kick, iunit)
 
    ! Units
    write(iunit,'(a1)', advance = 'no') '#'
    write(iunit,'(a20)', advance = 'no') str_center('t', 19)
    write(iunit,'(a20)', advance = 'yes') str_center('<O>(t)', 20)
    write(iunit,'(a1)', advance = 'no') '#'
    write(header_string, '(a)') '['//trim(units_abbrev(units_out%time))//']'
    write(iunit,'(a20)', advance = 'no')  str_center(trim(header_string), 19)
    write(header_string, '(a)') '['//trim(units_abbrev(units_out%length))//']'
    write(iunit,'(a20)', advance = 'yes')  str_center(trim(header_string), 20)
 
    do i = 0, time_steps
      write(iunit, '(2e20.8)') units_from_atomic(units_out%time, i*dt), units_from_atomic(ot_unit, ot(i))
    end do
 
    call io_close(iunit)
    pop_sub(write_ot)
  end subroutine write_ot
 
 
  ! ---------------------------------------------------------
  subroutine read_ot(namespace, nspin, order, nw_subtracted)
    type(namespace_t), intent(in) :: namespace
    integer, intent(out) :: nspin
    integer, intent(out) :: order
    integer, intent(out) :: nw_subtracted
 
    integer :: iunit, i, ierr
    character(len=100) :: line
    character(len=12)  :: dummychar
    real(real64) :: dummy, t1, t2
    type(unit_t) :: ot_unit
 
    push_sub(read_ot)
 
    iunit = io_open('ot', namespace, action='read', status='old')
    if (iunit == 0) then
      write(message(1),'(a)') 'A file called ot should be present and was not found.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    read(iunit, '(15x,i2)') nspin
    read(iunit, '(15x,i2)') order
    read(iunit, '(28x,i2)') nw_subtracted
    read(iunit, '(a)')      line
 
    i = index(line, 'Observable')
    if (index(line, 'Observable') /= 0) then
      observable(1) = -1
    elseif (index(line, '# O =') /= 0) then
      observable(1) = 0
      if (index(line,'x') /= 0) then
        observable(2) = 1
      elseif (index(line,'y') /= 0) then
        observable(2) = 2
      elseif (index(line,'z') /= 0) then
        observable(2) = 3
      end if
    elseif (index(line, '# (l, m) = ') /= 0) then
      read(line,'(a12,i1,a1,i2,a1)') dummychar(1:12), observable(1), dummychar(1:1), observable(2), dummychar(1:1)
    else
      write(message(1),'(a)') 'Problem reading "ot" file: could not figure out the shape'
      write(message(2),'(a)') 'of the observation operator.'
      call messages_fatal(2, namespace=namespace)
    end if
 
    call kick_read(kick, iunit, namespace)
    read(iunit, '(a)')  line
    read(iunit, '(a)')  line
    call io_skip_header(iunit)
 
    ! Figure out about the units of the file
    call unit_system_from_file(units, "ot", namespace, ierr)
    if (ierr /= 0) then
      write(message(1), '(a)') 'Could not figure out the units in file "ot".'
      call messages_fatal(1, namespace=namespace)
    end if
 
    select case (observable(1))
    case (-1)
      if (kick%n_multipoles > 0) then
        ot_unit = units_out%length**kick%l(1)
      else
        ot_unit = units_out%length
      end if
    case (0)
      ot_unit = units_out%length
    case default
      ot_unit = units_out%length**observable(1)
    end select
 
    ! count number of time_steps
    time_steps = 0
    do
      read(iunit, *, end=100) dummy
      time_steps = time_steps + 1
      if (time_steps == 1) t1 = dummy
      if (time_steps == 2) t2 = dummy
    end do
100 continue
    dt = units_to_atomic(units%time, (t2 - t1)) ! units_out is OK
 
    call io_skip_header(iunit)
 
    safe_allocate(ot(0:time_steps))
 
    do i = 0, time_steps-1
      read(iunit, *) dummy, ot(i)
      ot(i) = units_to_atomic(ot_unit, ot(i))
    end do
 
    pop_sub(read_ot)
  end subroutine read_ot
 
 
  ! ---------------------------------------------------------
  subroutine print_omega_file(namespace, omega, search_interval, final_time, nfrequencies)
    type(namespace_t), intent(in) :: namespace
    real(real64),   intent(inout) :: omega
    real(real64),   intent(inout) :: search_interval
    real(real64),   intent(inout) :: final_time
    integer, intent(inout) :: nfrequencies
 
    integer :: iunit, i, ierr, nspin, order, nw_subtracted
    logical :: file_exists
    character(len=20) :: header_string
    real(real64), allocatable :: warray(:), tarray(:)
    real(real64) :: leftbound, rightbound, dw, w, aw
 
    push_sub(print_omega_file)
 
    ! First, let us check that the file "ot" exists.
    inquire(file="ot", exist  = file_exists)
    if (.not. file_exists) then
      write(message(1),'(a)') "Could not find 'ot' file."
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! Now, we should find out which units the file "ot" has.
    call unit_system_from_file(units, "ot", namespace, ierr)
    if (ierr /= 0) then
      write(message(1),'(a)') "Could not retrieve units in the 'ot' file."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (omega > m_zero) then
      omega = units_to_atomic(units%energy, omega)
    else
      omega = m_half
    end if
 
    if (search_interval > m_zero) then
      search_interval = units_to_atomic(units%energy, search_interval)
    else
      search_interval = m_half
    end if
 
    leftbound = omega - search_interval
    rightbound = omega + search_interval
 
    safe_allocate(warray(1:nfrequencies))
    safe_allocate(tarray(1:nfrequencies))
 
    call read_ot(namespace, nspin, order, nw_subtracted)
 
    if (mod(order, 2) == 1) then
      mode = sine_transform
    else
      mode = cosine_transform
    end if
 
    if (final_time > m_zero) then
      total_time = units_to_atomic(units%time, final_time)
      if (total_time > dt*time_steps) then
        total_time = dt*time_steps
        write(message(1), '(a)')        'The requested total time to process is larger than the time available in the input file.'
        write(message(2), '(a,f8.4,a)') 'The time has been adjusted to ', &
          units_from_atomic(units_out%time, total_time), trim(units_abbrev(units_out%time))
        call messages_warning(2, namespace=namespace)
      end if
      time_steps = int(total_time / dt)
      total_time = time_steps * dt
    else
      total_time = dt*time_steps
    end if
 
    warray = m_zero
    tarray = m_zero
    dw = (rightbound-leftbound) / (nfrequencies - 1)
    do i = 1, nfrequencies
      w = leftbound + (i-1)*dw
      warray(i) = w
      call ft(w, aw)
      tarray(i) = aw
    end do
 
    iunit = io_open('omega', namespace, action='write', status='replace')
    write(iunit, '(a15,i2)')      '# nspin        ', nspin
    call kick_write(kick, iunit)
    write(iunit, '(a)') '#%'
    write(iunit, '(a1,a20)', advance = 'no') '#', str_center("omega", 20)
    write(header_string,'(a)') 'F(omega)'
    write(iunit, '(a20)', advance = 'yes') str_center(trim(header_string), 20)
    write(iunit, '(a1,a20)', advance = 'no') '#', str_center('['//trim(units_abbrev(units_out%energy)) // ']', 20)
    ! Here we should print the units of the transform.
    write(iunit, '(a)', advance = 'yes')
 
    do i = 1, nfrequencies
      write(iunit,'(2e20.8)') units_from_atomic(units_out%energy, warray(i)), &
        tarray(i)
    end do
 
    safe_deallocate_a(warray)
    safe_deallocate_a(tarray)
 
    pop_sub(print_omega_file)
  end subroutine print_omega_file
  ! ---------------------------------------------------------
 
end module oscillator_strength_oct_m
 
! ---------------------------------------------------------
! ---------------------------------------------------------
! ---------------------------------------------------------
program oscillator_strength
  use command_line_oct_m
  use global_oct_m
  use io_oct_m
  use messages_oct_m
  use oscillator_strength_oct_m
 
  implicit none
 
  integer :: run_mode, order, nfrequencies, ierr, nresonances
  real(real64) :: omega, search_interval, final_time, damping
  integer, parameter ::                    &
    ANALYZE_NTHORDER_SIGNAL           = 1, &
    generate_nthorder_signal          = 2, &
    read_resonances_from_file         = 3, &
    generate_omega_file               = 4
  character(len=100) :: ffile
 
  ! Reads the information passed through the command line options (if available).
  call getopt_init(ierr)
  if (ierr /= 0) then
    message(1) = "Your Fortran compiler doesn't support command-line arguments;"
    message(2) = "the oct-oscillator-strength command is not available."
    call messages_fatal(2)
  end if
 
  ! Set the default values
  run_mode        = analyze_nthorder_signal
  omega           = - m_one
  search_interval = - m_one
  order           = 1
  nfrequencies    = 1000
  final_time      = - m_one
  nresonances     = 1
  observable(1)   = -1
  observable(2)   = 0
  ffile           = ""
  damping         = 0.1_real64/27.2114_real64 ! This is the usual damping factor used in the literature.
 
  ! Get the parameters from the command line.
  call getopt_oscillator_strength(run_mode, omega, search_interval,             &
    order, nresonances, nfrequencies, final_time, &
    observable(1), observable(2), damping, ffile)
  call getopt_end()
 
  ! Initialize stuff
  call init_octopus_globals(serial_dummy_comm)
  call parser_init()
  call io_init(defaults = .true.)
 
  select case (run_mode)
  case (generate_nthorder_signal)
    call generate_signal(global_namespace, order, observable)
  case (analyze_nthorder_signal)
    call analyze_signal(order, omega, search_interval, final_time, nresonances, nfrequencies, damping, &
      global_namespace)
  case (read_resonances_from_file)
    call read_resonances_file(order, ffile, global_namespace, search_interval, final_time, nfrequencies)
  case (generate_omega_file)
    call print_omega_file(global_namespace, omega, search_interval, final_time, nfrequencies)
  case default
  end select
 
  call io_end()
  call parser_end()
  call global_end()
 
end program oscillator_strength
! ---------------------------------------------------------
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: