subroutine output_spectrum(neigv, eigv, spec_low, spec_high)
use util_mod, only: get_free_unit
integer, intent(in) :: neigv
real(dp), intent(in) :: eigv(neigv)
! Data for the testsuite to use.
real(dp), intent(out) :: spec_low, spec_high
integer :: i, j, min_vec, temp_unit
complex(dp) :: omega
complex(dp) :: spectral_weight
complex(dp), parameter :: iU = (0.0, 1.0)
if (iProcIndex /= root) return
temp_unit = get_free_unit()
open(temp_unit, file='SPECTRAL_DATA', status='replace')
write(temp_unit, '(1x,a53)') "Eigenvalues and left and right transition amplitudes:"
do i = 1, neigv
write(temp_unit, '(1x,i7,5x,f15.10,5x,f15.10,5x,f15.10)') &
i, eigv(i), trans_amps_left(i), trans_amps_right(i)
end do
close(temp_unit)
write(stdout, '(1x,a5,18X,a15)') "Omega", "Spectral_weight"
! Do we include the ground state in the spectrum or not?
if (tIncludeGroundSpectral) then
min_vec = 1
else
min_vec = 2
end if
omega = min_omega_spectral
if (tIWSpec) omega = omega * iU
do i = 1, nomega_spectral + 1
spectral_weight = 0.0_dp
do j = min_vec, neigv
if (tIWSpec) then
spectral_weight = spectral_weight + (trans_amps_left(j) * &
trans_amps_right(j)) * (1.0 / ((spectral_ground_energy - eigv(j)) &
+ omega) + 1.0 / (-1.0 * (spectral_ground_energy - eigv(j)) &
+ omega))
else
spectral_weight = spectral_weight + &
(trans_amps_left(j) * trans_amps_right(j) * spectral_broadening) / &
(pi * (spectral_broadening**2 + (spectral_ground_energy - eigv(j) + omega)**2))
end if
end do
write(stdout, '(f18.12, 4x, f18.12, 4x, f18.12, 4x, f18.12)') real(real(omega)), &
real(aimag(omega)), real(real(spectral_weight)), real(aimag(spectral_weight))
if (tIWSpec) then
omega = omega + iU * delta_omega_spectral
else
omega = omega + delta_omega_spectral
end if
! Store the values of the spectrum for the highest and lowest
! values of omega for the testsuite to use.
if (i == 1) spec_low = real(spectral_weight, dp)
if (i == nomega_spectral + 1) spec_high = real(spectral_weight, dp)
end do
end subroutine output_spectrum