tdtdm.F90

Go to the documentation of this file.

!! Copyright (C) 2019-2021 N. Tancogne-Dejean
!!
!! 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 tdtdm_main_oct_m
  use batch_oct_m
  use calc_mode_par_oct_m
  use comm_oct_m
  use debug_oct_m
  use electrons_oct_m
  use electron_space_oct_m
  use fft_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use io_function_oct_m
  use kpoints_oct_m
  use lalg_basic_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use restart_oct_m
  use spectrum_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use states_elec_restart_oct_m
  use symmetries_oct_m
  use symm_op_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public :: tdtdm_main
 
contains
 
  subroutine tdtdm_main() bind(c)
    integer :: in_file, ii, jj, kk, ierr, ip_h, irow, ifreq, nrow, it
    integer :: ik, ist, uist, istep, ikpoint, irep, out_file, iop, idim
    integer :: time_steps, energy_steps, istart, iend, ntiter, Nreplica, Ntrans
    real(real64)   :: dt, tt, weight, kpoint(3), kpoint_sym(3), kred(3), kred_sym(3)
    real(real64)   :: xx_h_sym(3)
    integer :: irep_h, ip_h_sym, rankmin
    real(real64)   :: start_time, dmin
    real(real64), allocatable :: Et(:), ftreal(:, :, :), ftimag(:, :, :), tmp(:), omega(:)
    complex(real64), allocatable :: Xiak(:,:,:), Yiak(:,:,:)
    real(real64), allocatable :: proj_r(:,:,:,:), proj_i(:,:,:,:)
    real(real64), allocatable :: proj_r_corr(:,:), proj_i_corr(:,:), centers(:,:)
    complex(real64), allocatable :: tdm(:,:), tdm_1D(:,:,:,:)
    complex(real64), allocatable, target :: psi(:,:), upsi(:,:)
    complex(real64), allocatable ::  phase(:,:,:), ftcmplx(:,:)
    complex(real64), contiguous, pointer :: psi_sym(:,:), upsi_sym(:,:)
    type(spectrum_t)  :: spectrum
    type(electrons_t), pointer :: sys
    type(batch_t)     :: projb_r, projb_i, ftrealb, ftimagb
    character(len=MAX_PATH_LEN) :: fname
    type(states_elec_t), pointer :: st
    type(states_elec_t) :: gs_st
    type(restart_t) :: restart
    type(unit_t) :: fn_unit
    integer :: kpt_start, kpt_end, supercell(3), nomega, ncols
    type(block_t) :: blk
    real(real64) :: pos_h(3), norm
 
    ! Initializion
    call global_init()
    call parser_init()
 
    call messages_init()
    call io_init()
 
    call profiling_init(global_namespace)
 
    call messages_experimental("oct-tdtdm utility")
    call fft_all_init(global_namespace)
    call unit_system_init(global_namespace)
 
    call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
    sys => electrons_t(global_namespace, mpi_world)
 
    call spectrum_init(spectrum, global_namespace)
 
    st => sys%st
 
    if(sys%st%d%ispin == spinors) then
      call messages_not_implemented('oct-tdtdm with spinors')
    end if
 
    if(st%parallel_in_states) then
      call messages_not_implemented("oct-tdtdm with states parallelization")
    end if
 
    if(sys%gr%parallel_in_domains) then
      call messages_not_implemented("oct-tdtdm with domain parallelization")
    end if
 
    !%Variable TDTDMFrequencies
    !%Type block
    !%Section Utilities::oct-tdtdm
    !%Description
    !% This block defines for which frequencies the analysis is performed.
    !%
    !% Each row of the block indicates a frequency.
    !%End
    if (parse_block(global_namespace, 'TDTDMFrequencies', blk) == 0) then
 
      nrow = parse_block_n(blk)
      nomega = nrow
 
      safe_allocate(omega(1:nrow))
      !read frequencies
      do irow = 0, nrow-1
        call parse_block_float(blk, irow, 0, omega(irow+1))
      end do
 
      call parse_block_end(blk)
    else
      message(1) = "oct-tdtdm: TDTDMFrequencies must be defined."
      call messages_fatal(1)
    end if
 
    ! We check that the resonant and antiresonant transitions are contained in the
    ! energy range of the Fourier transforms
    if(any(omega > spectrum%max_energy)) then
      message(1) = "One requested frequecy is larger than PropagationSpectrumMaxEnergy."
      message(2) = "Please increase the value of PropagationSpectrumMaxEnergy."
      call messages_fatal(2)
    end if
    if(any(omega > -spectrum%min_energy)) then
      message(1) = "One requested frequency is larger than -PropagationSpectrumMinEnergy."
      message(2) = "Please decrease the value of PropagationSpectrumMinEnergy."
      call messages_fatal(2)
    end if
 
 
    call states_elec_copy(gs_st, st, exclude_wfns = .true., exclude_eigenval = .true.)
 
    safe_deallocate_a(gs_st%node)
 
    call restart%init(global_namespace, restart_proj, restart_type_load, sys%mc, ierr, mesh=sys%gr)
    if(ierr == 0) call states_elec_look(restart, ii, jj, gs_st%nst, ierr)
    if(ierr /= 0) then
      message(1) = "oct-tdtdm: Unable to read states information."
      call messages_fatal(1)
    end if
 
    ! allocate memory
    safe_allocate(gs_st%occ(1:gs_st%nst, 1:gs_st%nik))
    safe_allocate(gs_st%eigenval(1:gs_st%nst, 1:gs_st%nik))
 
    ! We want all the task to have all the states
    ! States can be distibuted for the states we propagate.
    safe_allocate(gs_st%node(1:gs_st%nst))
    gs_st%node(:)  = 0
    call kpoints_distribute(gs_st, sys%mc)
    call states_elec_distribute_nodes(gs_st, global_namespace, sys%mc)
 
    kpt_start = gs_st%d%kpt%start
    kpt_end = gs_st%d%kpt%end
 
    kpoint = m_zero
 
    gs_st%eigenval = huge(gs_st%eigenval)
    gs_st%occ      = m_zero
    if(gs_st%d%ispin == spinors) then
      safe_deallocate_a(gs_st%spin)
      safe_allocate(gs_st%spin(1:3, 1:gs_st%nst, 1:gs_st%nik))
    end if
 
    call states_elec_allocate_wfns(gs_st, sys%gr, type_cmplx)
    call states_elec_load(restart, global_namespace, sys%space, gs_st, sys%gr, sys%kpoints, fixed_occ=.true., ierr=ierr)
    if(ierr /= 0 .and. ierr /= (gs_st%st_end-gs_st%st_start+1)*(kpt_end-kpt_start+1)*gs_st%d%dim) then
      message(1) = "oct-tdtdm: Unable to read wavefunctions for TDOutput."
      call messages_fatal(1)
    end if
    call restart%end()
 
 
    in_file = io_open('td.general/projections', action='read', status='old')
    call io_skip_header(in_file)
    call spectrum_count_time_steps(global_namespace, in_file, time_steps, dt)
    dt = units_to_atomic(units_out%time, dt)
 
 
    safe_allocate(tmp(1:st%nst*gs_st%nst*st%nik*2))
    safe_allocate(proj_r(1:time_steps, 1:gs_st%nst, 1:st%nst, 1:st%nik))
    safe_allocate(proj_i(1:time_steps, 1:gs_st%nst, 1:st%nst, 1:st%nik))
 
 
    call io_skip_header(in_file)
 
    do ii = 1, time_steps
      read(in_file, *) jj, tt, (tmp(kk), kk = 1, st%nst*gs_st%nst*st%nik*2)
      do ik = 1, st%nik
        do ist = 1, st%nst
          do uist = 1, gs_st%nst
            jj  = (ik-1)*st%nst*gs_st%nst + (ist-1)*gs_st%nst + uist
            proj_r(ii, uist, ist, ik) = tmp((jj-1)*2+1)
            ! Here we add a minus sign, as we want to get <\phi_0 | \psi(t)>
            ! and td_occup computes the complex conjugaute of this
            proj_i(ii, uist, ist, ik) = -tmp((jj-1)*2+2)
          end do
        end do
      end do
    end do
    safe_deallocate_a(tmp)
 
    call io_close(in_file)
 
    write(message(1), '(a, i7, a)') "oct-tdtdm: Read ", time_steps, " steps from file '"// &
      trim(io_workpath('td.general/projections'))//"'"
    call messages_info(1)
 
    start_time = spectrum%start_time
 
    ! Phase correction of the projections before doing the Fourier transforms
    ! See Eq. (5) of Williams et al., JCTC 17, 1795 (2021)
    ! We need to multiply C_ik(t)e^{-ie_kt} (the projection of \phi_i(t) on \phi_k^GS)
    ! by e^{ie_it}, which is obtained by the cc of the projection of \phi_i(t) on \phi_i^GS
    ! Here we only care about optical transitions (so TD occupied to GS unocc)
    safe_allocate(proj_r_corr(1:time_steps, 1:gs_st%nst*st%nst*(kpt_end-kpt_start+1)))
    safe_allocate(proj_i_corr(1:time_steps, 1:gs_st%nst*st%nst*(kpt_end-kpt_start+1)))
    proj_r_corr = m_zero
    proj_i_corr = m_zero
    do ik = kpt_start, kpt_end
      do ist = 1, st%nst
        do uist = ist+1, gs_st%nst
          jj = (ik-kpt_start)*st%nst*gs_st%nst+(ist-1)*gs_st%nst+uist
          do ii = 1, time_steps
            norm = hypot(proj_r(ii, ist, ist, ik),proj_i(ii, ist, ist, ik))
            if (norm > m_epsilon) then
              proj_r_corr(ii, jj) = (proj_r(ii, uist, ist, ik) * proj_r(ii, ist, ist, ik) &
                + proj_i(ii, uist, ist, ik) * proj_i(ii, ist, ist, ik))/norm
              proj_i_corr(ii, jj) =(-proj_r(ii, uist, ist, ik) * proj_i(ii, ist, ist, ik) &
                + proj_i(ii, uist, ist, ik) * proj_r(ii, ist, ist, ik))/norm
            else
              proj_r_corr(ii, jj) = m_zero
              proj_i_corr(ii, jj) = m_zero
            end if
          end do
        end do
      end do
    end do
 
    safe_deallocate_a(proj_r)
    safe_deallocate_a(proj_i)
 
    ! 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(ftreal(1:energy_steps, 1:st%nst*gs_st%nst*(kpt_end-kpt_start+1), 1:2))
    safe_allocate(ftimag(1:energy_steps, 1:st%nst*gs_st%nst*(kpt_end-kpt_start+1), 1:2))
 
    call batch_init(projb_r, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), proj_r_corr)
    call batch_init(projb_i, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), proj_i_corr)
    call batch_init(ftrealb, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), ftreal(:,:,1))
    call batch_init(ftimagb, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), ftimag(:,:,1))
 
    write(message(1), '(a)') "oct-tdtdm: Fourier transforming real part of the projections"
    call messages_info(1)
 
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
      istart, iend, spectrum%start_time, dt, projb_r, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftrealb)
 
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
      istart, iend, spectrum%start_time, dt, projb_r, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftimagb)
 
    call ftrealb%end()
    call ftimagb%end()
 
    safe_allocate(ftcmplx(1:energy_steps, 1:st%nst*gs_st%nst*(kpt_end-kpt_start+1)))
    do ii = 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1)
      ftcmplx(1:energy_steps,ii) =  cmplx(ftreal(1:energy_steps,ii,1), ftimag(1:energy_steps,ii,1), real64)
    end do
 
    write(message(1), '(a)') "oct-tdtdm: Fourier transforming imaginary part of the projections"
    call messages_info(1)
 
    call batch_init(ftrealb, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), ftreal(:,:,2))
    call batch_init(ftimagb, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), ftimag(:,:,2))
 
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
      istart, iend, spectrum%start_time, dt, projb_i, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftrealb)
 
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
      istart, iend, spectrum%start_time, dt, projb_i, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftimagb)
 
    call projb_i%end()
    call projb_r%end()
    call ftrealb%end()
    call ftimagb%end()
    safe_deallocate_a(proj_r_corr)
    safe_deallocate_a(proj_i_corr)
 
    do ii = 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1)
      ftcmplx(1:energy_steps,ii) = ftcmplx(1:energy_steps,ii) + m_zi*ftreal(1:energy_steps,ii,2) - ftimag(1:energy_steps,ii,2)
    end do
 
    safe_deallocate_a(ftreal)
    safe_deallocate_a(ftimag)
 
    write(message(1), '(a)') "oct-tdtdm: Constructing the two-particle wavefunctions."
    call messages_info(1)
 
    !%Variable SupercellDimensions
    !%Type block
    !%Default KPointsGrid
    !%Section Utilities::oct-tdtdm
    !%Description
    !% This block allows to specify the size of the supercell used to plot excitonic wavefunctions.
    !% If not specified, the code uses the number of k-points for defining the size of the supercell.
    !%End
    if (parse_is_defined(sys%namespace, 'SupercellDimensions')) then
      if (parse_block(sys%namespace, 'SupercellDimensions', blk) == 0) then
        ncols = parse_block_cols(blk, 0)
        if (ncols /= sys%space%dim) then
          write(message(1),'(a,i3,a,i3)') 'SupercellDimensions has ', ncols, ' columns but must have ', sys%space%dim
          call messages_fatal(1, namespace=sys%namespace)
        end if
        do ii = 1, sys%space%dim
          call parse_block_integer(blk, 0, ii - 1, supercell(ii))
        end do
 
        call parse_block_end(blk)
      end if
    else
      supercell(1:sys%space%dim) = sys%kpoints%nik_axis(1:sys%space%dim)
    end if
 
    nreplica = product(supercell(1:sys%space%dim))
 
    ! The center of each replica of the unit cell
    safe_allocate(centers(1:sys%space%dim, 1:nreplica))
    irep = 1
    do ii = 0, supercell(1)-1
      do jj = 0, supercell(2)-1
        do kk = 0, supercell(3)-1
          centers(1, irep) = -floor((supercell(1)-1)/m_two)+ii
          centers(2, irep) = -floor((supercell(2)-1)/m_two)+jj
          centers(3, irep) = -floor((supercell(3)-1)/m_two)+kk
          centers(:, irep) = matmul(sys%ions%latt%rlattice, centers(:, irep))
          irep = irep + 1
        end do
      end do
    end do
 
    ! The phase for each center
    irep = 0
    do ik = kpt_start, kpt_end
      ikpoint = gs_st%d%get_kpoint_index(ik)
      irep = max(irep, kpoints_get_num_symmetry_ops(sys%kpoints, ikpoint))
    end do
    safe_allocate(phase(kpt_start:kpt_end, 1:irep, 1:nreplica))
    do irep = 1, nreplica
      do ik = kpt_start, kpt_end
        ikpoint = gs_st%d%get_kpoint_index(ik)
        kpoint(1:sys%space%dim) = sys%kpoints%get_point(ikpoint)
        do ii = 1, kpoints_get_num_symmetry_ops(sys%kpoints, ikpoint)
          iop = kpoints_get_symmetry_ops(sys%kpoints, ikpoint, ii)
 
          if (sys%kpoints%use_symmetries) then !We apply the symmetry
            call kpoints_to_reduced(sys%kpoints%latt, kpoint, kred)
            call symmetries_apply_kpoint_red(sys%kpoints%symm, iop, kred, kred_sym)
            call kpoints_to_absolute(sys%kpoints%latt, kred_sym, kpoint_sym)
          else
            kpoint_sym = kpoint
          end if
          phase(ik, ii, irep) = exp(-m_zi*sum(kpoint_sym(1:sys%space%dim)*centers(:, irep)))
        end do
      end do
    end do
 
    ! Position of the hole, here assumed to be on top of the first atom
    ! To be obtained from the input file
    if(sys%space%dim > 1) then
      call tdtdm_get_hole_position(pos_h, ip_h)
    end if
 
    ntrans = 0
    ! Here we assume that there is a clear gap, so the information at Gamma is enough
    do ist = 1, gs_st%nst
      if(abs(gs_st%occ(ist, 1)) < m_min_occ) cycle
 
      do uist = 1, gs_st%nst
        if(abs(gs_st%occ(uist, 1)) > m_min_occ) cycle
        weight = gs_st%kweights(1) * (gs_st%occ(ist, 1)-gs_st%occ(uist, 1))
        if(abs(weight) < m_min_occ) cycle
        ntrans = ntrans + 1
      end do
    end do
    if(ntrans == 0) then
      write(message(1), '(a)') "oct-tdtdm: No transition found."
      write(message(2), '(a)') "Please check that unoccupied states are included in the ground state calculation."
      call messages_fatal(2)
    end if
 
    safe_allocate(xiak(1:st%nst, 1:gs_st%nst, 1:st%nik))
    safe_allocate(yiak(1:st%nst, 1:gs_st%nst, 1:st%nik))
    safe_allocate(et(1:ntrans*st%nik))
    safe_allocate(psi(1:sys%gr%np, 1:gs_st%d%dim))
    safe_allocate(upsi(1:sys%gr%np, 1:gs_st%d%dim))
 
    if(sys%kpoints%use_symmetries) then
      safe_allocate(psi_sym(1:sys%gr%np, 1:st%d%dim))
      safe_allocate(upsi_sym(1:sys%gr%np, 1:st%d%dim))
    end if
 
    select case(sys%space%dim)
    case(2,3)
      safe_allocate(tdm(1:sys%gr%np, 1:nreplica))
    case(1)
      safe_allocate(tdm_1d(1:sys%gr%np, 1:sys%gr%np, 1:nreplica, 1:nreplica))
    end select
 
    do ifreq = 1, nomega
 
      write(message(1), '(a, f6.4, a)') "oct-tdtdm: Constructing the two-particle wavefunction at ", omega(ifreq), " Ha."
      call messages_info(1)
 
      select case(sys%space%dim)
      case(2,3)
        tdm = m_z0
      case(1)
        tdm_1d = m_z0
      end select
 
      et = m_zero
      xiak = m_z0
      yiak = m_z0
 
      ! Local transition index
      it = (kpt_start-1)*ntrans + 1
 
      do ik = kpt_start, kpt_end
        ikpoint = st%d%get_kpoint_index(ik)
 
        do ist = 1, st%nst
          if(abs(gs_st%occ(ist, ik)) < m_min_occ) cycle
 
          call states_elec_get_state(gs_st, sys%gr, ist, ik, psi)
          if (sys%hm%phase%is_allocated()) then
            call sys%hm%phase%apply_to_single(psi, sys%gr%np, gs_st%d%dim, ik, .false.)
          end if
 
          do uist = 1, gs_st%nst
            if(abs(gs_st%occ(uist, ik)) > m_min_occ) cycle
 
            ! For a given requested frequency, we get the corresponding values of Xia and Yia
            ! One correspond to the +\Omega frequency, the other one to the -\Omega frequency
            ! For Xiak, we use the fact that TF[f*](\Omega) = (TF[f](-\Omega))^*
            jj = (ik-kpt_start)*st%nst*gs_st%nst+(ist-1)*gs_st%nst+uist
            istep = int((+omega(ifreq)-spectrum%min_energy)/spectrum%energy_step)
            xiak(ist, uist, ik) = conjg(ftcmplx(istep, jj))
            istep = int((+omega(ifreq)-spectrum%min_energy)/spectrum%energy_step)
            yiak(ist, uist, ik) = ftcmplx(istep, jj)
 
 
            weight = gs_st%kweights(ik) * (gs_st%occ(ist, ik)-gs_st%occ(uist, ik)) &
              / kpoints_get_num_symmetry_ops(sys%kpoints, ikpoint)
            if(abs(weight) < m_epsilon) cycle
 
            call states_elec_get_state(gs_st, sys%gr, uist, ik, upsi)
            if(sys%hm%phase%is_allocated()) then
              call sys%hm%phase%apply_to_single(upsi, sys%gr%np, st%d%dim, ik, .false.)
            end if
 
            do ii = 1, kpoints_get_num_symmetry_ops(sys%kpoints, ikpoint)
              iop = kpoints_get_symmetry_ops(sys%kpoints, ikpoint, ii)
 
              if(sys%kpoints%use_symmetries) then
                do idim = 1, st%d%dim
                  call zgrid_symmetrize_single(sys%gr, iop, psi(:,idim), psi_sym(:,idim))
                  call zgrid_symmetrize_single(sys%gr, iop, upsi(:,idim), upsi_sym(:,idim))
                end do
 
                ! We need to get the position of the hole after applying the symmetry operation too
                xx_h_sym = symm_op_apply_cart(sys%kpoints%symm%ops(iop), pos_h)
                xx_h_sym = sys%ions%latt%fold_into_cell(xx_h_sym)
                ! At the moment, we ignore rankmin
                assert(.not.sys%gr%parallel_in_domains)
                ip_h_sym = mesh_nearest_point(sys%gr, xx_h_sym, dmin, rankmin)
              else
                psi_sym => psi
                upsi_sym => upsi
                ip_h_sym = ip_h
              end if
 
              ! We now compute the single mode TDTDM
              ! See Eq. (5) of Williams et al., JCTC 17, 1795 (2021)
              ! We take here the complex conjugate of the 2-body wavefunction
              select case(sys%space%dim)
              case(2,3)
                do irep = 1, nreplica
                  call lalg_axpy(sys%gr%np, phase(ik, ii, irep) * weight &
                    * conjg(xiak(ist,uist,ik))*conjg(psi_sym(ip_h_sym,1)), upsi_sym(:, 1), tdm(:,irep))
                  call lalg_axpy(sys%gr%np, phase(ik, ii, irep) * weight &
                    * yiak(ist,uist,ik)*conjg(upsi_sym(ip_h_sym,1)), psi_sym(:, 1), tdm(:,irep))
                end do
              case(1)
                ! In the 1D case, we contruct the full TDTDM of r_e, r_h
                do irep_h = 1, nreplica
                  do irep = 1, nreplica
                    do ip_h = 1, sys%gr%np
                      call lalg_axpy(sys%gr%np, phase(ik, ii, irep) * conjg(phase(ik, ii, irep_h)) &
                        * weight * conjg(xiak(ist,uist,ik)) * conjg(psi_sym(ip_h,1)), &
                        upsi_sym(:, 1), tdm_1d(:, ip_h, irep, irep_h))
                      call lalg_axpy(sys%gr%np, phase(ik, ii, irep) * conjg(phase(ik, ii, irep_h)) &
                        * weight * conjg(yiak(ist,uist,ik)) * conjg(upsi_sym(ip_h,1)), &
                        psi_sym(:, 1), tdm_1d(:, ip_h, irep, irep_h))
                    end do
                  end do
                end do
              end select
 
            end do ! ii
 
            et(it) = gs_st%eigenval(uist, ik) - gs_st%eigenval(ist, ik)
            it = it + 1
          end do
        end do
      end do
 
      if(gs_st%d%kpt%parallel) then
        if(sys%space%dim > 1) then
          call comm_allreduce(gs_st%d%kpt%mpi_grp, tdm)
        else
          call comm_allreduce(gs_st%d%kpt%mpi_grp, tdm_1d)
        end if
        call comm_allreduce(gs_st%d%kpt%mpi_grp, et)
        call comm_allreduce(gs_st%d%kpt%mpi_grp, xiak)
        call comm_allreduce(gs_st%d%kpt%mpi_grp, yiak)
      end if
 
      call tdtdm_output_density()
 
      call tdtdm_excitonic_weight()
 
    end do ! ifreq
 
    safe_deallocate_a(et)
    safe_deallocate_a(xiak)
    safe_deallocate_a(yiak)
    safe_deallocate_a(tdm)
    safe_deallocate_a(tdm_1d)
 
    safe_deallocate_a(psi)
    safe_deallocate_a(upsi)
    if(sys%kpoints%use_symmetries) then
      safe_deallocate_p(psi_sym)
      safe_deallocate_p(upsi_sym)
    end if
    safe_deallocate_a(ftcmplx)
    safe_deallocate_a(centers)
    safe_deallocate_a(phase)
    safe_deallocate_a(omega)
 
    safe_deallocate_p(sys)
    call states_elec_end(gs_st)
    call fft_all_end()
    call io_end()
    call profiling_end(global_namespace)
    call messages_end()
    call parser_end()
    call global_end()
 
  contains
 
    ! -----------------------------------------------------------------
    ! Determines the position of the hole, either from the input or using the
    ! first atom in the cell.
    ! This returns the index of the point in the mesh closest to the position.
    subroutine tdtdm_get_hole_position(xx_h, ip_h)
      real(real64),   intent(out) :: xx_h(1:sys%space%dim)
      integer, intent(out) :: ip_h
 
      real(real64) :: dmin
      integer :: idir, rankmin
 
      push_sub(tdtdm_get_hole_position)
 
      !%Variable TDTDMHoleCoordinates
      !%Type float
      !%Section Utilities::oct-tdtdm
      !%Description
      !% The position of the hole used to compute the TDTDM,
      !% in Cartesian coordinates.
      !% Note that the code will use the closest grid point.
      !%
      !% The coordinates of the hole are specified in the following way
      !% <tt>%TDTDMHoleCoordinates
      !% <br>&nbsp;&nbsp;hole_x | hole_y | hole_z
      !% <br>%</tt>
      !%
      !% If TDTDMHoleCoordinates or TDTDMHoleReducedCoordinates are not specified,
      !% the code will use the coordinate of the first atom in the cell.
      !%End
 
      if(parse_block(global_namespace, 'TDTDMHoleCoordinates', blk) == 0) then
        if(parse_block_cols(blk,0) < sys%space%dim) then
          call messages_input_error(global_namespace, 'TDTDMHoleCoordinates')
        end if
        do idir = 1, sys%space%dim
          call parse_block_float(blk, 0, idir - 1, xx_h(idir), units_inp%length)
        end do
        call parse_block_end(blk)
      else
        !%Variable TDTDMHoleReducedCoordinates
        !%Type float
        !%Section Utilities::oct-tdtdm
        !%Description
        !% Same as TDTDMHoleCoordinates, except that coordinates are given in reduced coordinates
        !%End
 
        if(parse_block(global_namespace, 'TDTDMHoleReducedCoordinates', blk) == 0) then
          if(parse_block_cols(blk,0) < sys%space%dim) then
            call messages_input_error(global_namespace, 'TDTDMHoleReducedCoordinates')
          end if
          do idir = 1, sys%space%dim
            call parse_block_float(blk, 0, idir - 1, xx_h(idir))
          end do
          call parse_block_end(blk)
          xx_h = sys%ions%latt%red_to_cart(xx_h)
        else
          xx_h(1:sys%space%dim) = sys%ions%pos(1:sys%space%dim, 1)
        end if
      end if
 
      ! We bring back the hole into the cell
      xx_h = sys%ions%latt%fold_into_cell(xx_h)
 
      ! At the moment, we ignore rankmin
      assert(.not.sys%gr%parallel_in_domains)
      ip_h = mesh_nearest_point(sys%gr, xx_h, dmin, rankmin)
      write(message(1), '(a, 3(1x,f7.4,a))') "oct-tdtdm: Requesting the hole at (", xx_h(1), &
        ",", xx_h(2), ",", xx_h(3), ")."
      call mesh_r(sys%gr, ip_h, dmin, coords=xx_h)
      write(message(2), '(a, 3(1x,f7.4,a))') "oct-tdtdm: Setting the hole at (", xx_h(1), &
        ",", xx_h(2), ",", xx_h(3), ")."
 
      call messages_info(2)
 
      pop_sub(tdtdm_get_hole_position)
    end subroutine tdtdm_get_hole_position
 
    subroutine tdtdm_output_density()
      real(real64), allocatable :: den(:,:), den_1d(:,:,:,:)
      real(real64) :: norm, xx(3), xx_h(3)
      integer :: iunit
 
      push_sub(tdtdm_output_density)
 
      ! We compute the TDM density
      select case(sys%space%dim)
      case(2,3)
        safe_allocate(den(1:sys%gr%np, 1:nreplica))
        do irep = 1, nreplica
          do ii = 1, sys%gr%np
            den(ii, irep) = real(tdm(ii, irep)*conjg(tdm(ii, irep)), real64)
          end do
        end do
 
        ! Here we renormalize to avoid too small numbers in the outputs
        norm = maxval(den)
        call lalg_scal(sys%gr%np, nreplica, m_one/norm, den)
 
      case(1)
        safe_allocate(den_1d(1:sys%gr%np, 1:sys%gr%np, 1:nreplica, 1:nreplica))
        do irep_h = 1, nreplica
          do irep = 1, nreplica
            do ip_h = 1, sys%gr%np
              do ii = 1, sys%gr%np
                tdm_1d(ii, ip_h, irep, irep_h) = conjg(tdm_1d(ii, ip_h, irep, irep_h))
                den_1d(ii, ip_h, irep, irep_h) = real(tdm_1d(ii, ip_h, irep, irep_h)*conjg(tdm_1d(ii,ip_h, irep, irep_h)), real64)
              end do
            end do
          end do
        end do
      end select
 
      fn_unit = units_out%length**(-sys%space%dim)
 
      select case(sys%space%dim)
      case(2,3)
        write(fname, '(a, f0.4)') 'tdm_density-0', omega(ifreq)
        call io_function_output_supercell(io_function_fill_how("XCrySDen"), "td.general", fname, &
          sys%gr, sys%space, sys%ions%latt, den, centers, supercell, fn_unit, &
          ierr, global_namespace, pos=sys%ions%pos, atoms=sys%ions%atom, grp = st%dom_st_kpt_mpi_grp, extra_atom=pos_h)
 
        call io_function_output_supercell(io_function_fill_how("PlaneZ"), "td.general", fname, &
          sys%gr, sys%space, sys%ions%latt, den, centers, supercell, fn_unit, &
          ierr, global_namespace, grp = st%dom_st_kpt_mpi_grp)
 
        safe_deallocate_a(den)
 
      case(1)
 
        call tdtdm_get_hole_position(pos_h, ip_h)
        irep_h = floor(supercell(1)/m_two)
 
        write(fname, '(a, f0.4)') 'tdm_density-0', omega(ifreq)
        call io_function_output_supercell(io_function_fill_how("AxisX"), "td.general", fname, &
          sys%gr, sys%space, sys%ions%latt, &
          den_1d(:,ip_h,:,irep_h), centers, supercell, fn_unit, ierr, global_namespace, &
          grp = st%dom_st_kpt_mpi_grp)
 
        write(fname, '(a, f0.4)') 'tdm_wfn-0', omega(ifreq)
        call io_function_output_supercell(io_function_fill_how("AxisX"), "td.general", fname, &
          sys%gr, sys%space, sys%ions%latt, &
          tdm_1d(:,ip_h,:,irep_h), centers, supercell, fn_unit, ierr, global_namespace, &
          grp = st%dom_st_kpt_mpi_grp)
 
        assert(.not.sys%gr%parallel_in_domains)
        if (mpi_world%is_root()) then
          write(fname, '(a, f0.4)') 'td.general/tdm_density-0', omega(ifreq)
          iunit = io_open(fname, action='write')
          write(iunit, '(a)', iostat=ierr) '# r_e    r_h    Re(\Psi(r_e,r_h)) Im(\Psi(r_e,r_h)) |\Psi(r_e,r_h)|^2'
 
          do irep_h = 1, nreplica
            do ip_h = 1, sys%gr%np
              xx_h = units_from_atomic(units_out%length, mesh_x_global(sys%gr, i4_to_i8(ip_h)) &
                + centers(1:sys%space%dim, irep_h))
 
              do irep = 1, nreplica
                do ii = 1, sys%gr%np
                  xx = units_from_atomic(units_out%length, mesh_x_global(sys%gr, i4_to_i8(ii)) &
                    + centers(1:sys%space%dim, irep))
                  write(iunit, '(5es23.14E3)', iostat=ierr) xx(1), xx_h(1), &
                    real(units_from_atomic(fn_unit, tdm_1D(ii, ip_h, irep, irep_h)), real64)  ,&
                    aimag(units_from_atomic(fn_unit, tdm_1D(ii, ip_h, irep, irep_h))), &
                    units_from_atomic(fn_unit, den_1D(ii, ip_h, irep, irep_h))
                end do
              end do
            end do
          end do
        end if
 
        safe_deallocate_a(den_1d)
      end select
 
 
      pop_sub(tdtdm_output_density)
    end subroutine tdtdm_output_density
 
    subroutine tdtdm_excitonic_weight()
      real(real64), allocatable :: weight(:,:)
 
      if (.not. mpi_world%is_root()) return
 
      push_sub(tdtdm_excitonic_weight)
 
      safe_allocate(weight(1:st%nik, 1:gs_st%nst))
      weight = m_zero
 
      do ik = 1, st%nik
        do ist = 1, st%nst
          if(abs(gs_st%occ(ist, ik)) < m_min_occ) cycle
 
          do uist = ist+1, gs_st%nst
            if(abs(gs_st%occ(uist, ik)) > m_min_occ) cycle
 
            weight(ik, ist)  = weight(ik, ist) + abs(xiak(ist, uist, ik))**2
            weight(ik, uist) = weight(ik,uist) + abs(yiak(ist, uist, ik))**2
          end do
        end do
      end do
 
      write(fname, '(a, f0.4)') 'td.general/tdm_weights-0', omega(ifreq)
      out_file = io_open(fname, action='write')
      write(out_file, '(a)') '# ik - kx - ky - kz - sum weights - eigenval and weights(ist,ik) '
      do ik = 1, st%nik
        ikpoint = st%d%get_kpoint_index(ik)
        kpoint(1:sys%space%dim) = sys%kpoints%reduced%point1BZ(1:sys%space%dim,ikpoint)
        write(out_file, '(i4,4e15.6)', advance='no') ik, kpoint(1:3), sum(weight(ik, 1:gs_st%nst))
        do uist = 1, gs_st%nst-1
          write(out_file, '(2e15.6)', advance='no') gs_st%eigenval(uist, ik), weight(ik, uist)
        end do
        write(out_file, '(e15.6)') weight(ik, uist)
      end do
      call io_close(out_file)
 
      safe_deallocate_a(weight)
 
      pop_sub(tdtdm_excitonic_weight)
    end subroutine tdtdm_excitonic_weight
  end subroutine tdtdm_main
end module tdtdm_main_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: