Commit b6a0e13f authored by Micael Oliveira's avatar Micael Oliveira
Browse files

Merge branch 'td_pcm' into develop

parents 97765780 2c1715bd
Pipeline #57086802 failed with stage
in 0 seconds
......@@ -62,7 +62,10 @@ module epot_oct_m
epot_generate, &
epot_local_potential, &
epot_precalc_local_potential, &
epot_global_force
epot_global_force, &
epot_have_lasers, &
epot_have_kick, &
epot_have_external_potentials
integer, public, parameter :: &
CLASSICAL_NONE = 0, & !< no classical charges
......@@ -918,6 +921,51 @@ contains
POP_SUB(epot_global_force)
end subroutine epot_global_force
! ---------------------------------------------------------
logical function epot_have_lasers(ep)
type(epot_t), intent(in) :: ep
PUSH_SUB(epot_have_lasers)
epot_have_lasers = .false.
if( ep%no_lasers /= 0 ) epot_have_lasers = .true.
POP_SUB(epot_have_lasers)
end function epot_have_lasers
! ---------------------------------------------------------
logical function epot_have_kick(ep)
type(epot_t), intent(in) :: ep
PUSH_SUB(epot_have_kick)
epot_have_kick = .false.
if( ep%kick%delta_strength /= M_ZERO ) epot_have_kick = .true.
POP_SUB(epot_have_kick)
end function epot_have_kick
! ---------------------------------------------------------
logical function epot_have_external_potentials(ep)
type(epot_t), intent(in) :: ep
PUSH_SUB(epot_have_external_potentials)
epot_have_external_potentials = .false.
if( associated(ep%v_static) .or. associated(ep%E_field) .or. epot_have_lasers(ep) ) epot_have_external_potentials = .true.
POP_SUB(epot_have_external_potentials)
end function epot_have_external_potentials
end module epot_oct_m
!! Local Variables:
......
......@@ -414,11 +414,9 @@ contains
!%End
call parse_variable('HamiltonianApplyPacked', .true., hm%apply_packed)
external_potentials_present = associated(hm%ep%v_static) .or. &
associated(hm%ep%E_field) .or. &
associated(hm%ep%lasers)
external_potentials_present = epot_have_external_potentials(hm%ep)
kick_present = hm%ep%kick%delta_strength /= M_ZERO
kick_present = epot_have_kick(hm%ep)
call pcm_init(hm%pcm, geo, gr, st%qtot, st%val_charge, external_potentials_present, kick_present ) !< initializes PCM
if(hm%pcm%run_pcm .and. hm%theory_level /= KOHN_SHAM_DFT) call messages_not_implemented("PCM for TheoryLevel /= DFT")
......@@ -1327,10 +1325,18 @@ contains
forall (ip = 1:mesh%np) this%hm_base%potential(ip, ispin) = this%vhxc(ip, ispin) + this%ep%vpsl(ip)
!> Adds PCM contributions
if (this%pcm%run_pcm) then
forall (ip = 1:mesh%np)
this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
this%pcm%v_e_rs(ip) + this%pcm%v_n_rs(ip)
end forall
if (this%pcm%solute) then
forall (ip = 1:mesh%np)
this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
this%pcm%v_e_rs(ip) + this%pcm%v_n_rs(ip)
end forall
end if
if (this%pcm%localf) then
forall (ip = 1:mesh%np)
this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
this%pcm%v_ext_rs(ip)
end forall
end if
end if
if(this%bc%abtype == IMAGINARY_ABSORBING) then
......
......@@ -758,18 +758,19 @@ contains
call kick_function_get(mesh, kick, kick_function_interpolate, to_interpolate = .true.)
SAFE_ALLOCATE(kick_function_real(1:mesh%np_part))
kick_function_real = DREAL(kick_function_interpolate)
if ( pcm%kick_like .or. pcm%which_eps == 'drl' ) then
! computing kick-like polarization due to kick or initialize polarization due to kick for the Drude-Lorentz model
if ( pcm%kick_like ) then
! computing kick-like polarization due to kick
call pcm_calc_pot_rs(pcm, mesh, kick = kick%delta_strength * kick_function_real, kick_time = .true.)
else if ( .not.pcm%kick_like .and. pcm%which_eps == 'deb' ) then
! computing the kick-like part of polarization due to kick for Debye dielectric model
pcm%kick_like = .true.
call pcm_calc_pot_rs(pcm, mesh, kick = kick%delta_strength * kick_function_real, kick_time = .true.)
pcm%kick_like = .false.
else if ( .not.pcm%kick_like .and. pcm%which_eps == 'drl' ) then
POP_SUB(kick_pcm_function_get)
return
end if
if( pcm%kick_like .or. pcm%which_eps == 'deb' ) then
kick_pcm_function = pcm%v_kick_rs / kick%delta_strength
end if
kick_pcm_function = pcm%v_kick_rs / kick%delta_strength
end if
POP_SUB(kick_pcm_function_get)
......
This diff is collapsed.
......@@ -456,17 +456,11 @@ module pcm_eom_oct_m
if( .not.allocated(matqq) ) then
SAFE_ALLOCATE(matqq(nts_act,nts_act))
endif
else if( which_eom == 'external' ) then
else if( which_eom == 'external' .or. which_eom == 'justkick' ) then
SAFE_ALLOCATE(matq0_lf(nts_act,nts_act)) !< not used yet
SAFE_ALLOCATE(matqd_lf(nts_act,nts_act))
SAFE_ALLOCATE(matqv_lf(nts_act,nts_act))
if( (.not.allocated(matqq)) .and. which_eps == 'deb' ) then
SAFE_ALLOCATE(matqq(nts_act,nts_act))
endif
else if( which_eom == 'justkick' ) then
SAFE_ALLOCATE(matqv_lf(nts_act,nts_act))
SAFE_ALLOCATE(matqd_lf(nts_act,nts_act))
if( (.not.allocated(matqq)) .and. which_eps == 'deb' ) then
if( .not.allocated(matqq) ) then
SAFE_ALLOCATE(matqq(nts_act,nts_act))
endif
endif
......@@ -544,6 +538,9 @@ module pcm_eom_oct_m
if( allocated(matqd) ) then
SAFE_DEALLOCATE_A(matqd)
endif
if( allocated(matq0_lf) ) then
SAFE_DEALLOCATE_A(matq0_lf)
endif
if( allocated(matqd_lf) ) then
SAFE_DEALLOCATE_A(matqd_lf)
endif
......@@ -664,7 +661,7 @@ module pcm_eom_oct_m
enddo
if( which_eom == 'electron' ) then
matq0=-matmul(scr1,scr4) !< from Eq.(14) and (18) for eps_0 in Ref.1
else if( which_eom == 'external' ) then
else if( which_eom == 'external' .or. which_eom == 'justkick' ) then
matq0_lf=-matmul(scr1,scr4) !< local field analogous !< not used yet
endif
do i=1,nts_act
......
......@@ -162,6 +162,24 @@ contains
call hamiltonian_init(hm, sys%gr, sys%geo, sys%st, sys%ks%theory_level, &
sys%ks%xc_family, family_is_mgga_with_exc(sys%ks%xc, sys%st%d%nspin))
if(hm%pcm%run_pcm .and. calc_mode_id /= CM_GS .and. calc_mode_id /= CM_TD ) then
call messages_not_implemented("PCM for CalculationMode /= gs or td")
else if(hm%pcm%run_pcm .and. calc_mode_id == CM_TD ) then
call messages_experimental("PCM for CalculationMode = td")
else if ( hm%pcm%epsilon_infty /= hm%pcm%epsilon_0 .and. hm%pcm%noneq .and. calc_mode_id == CM_GS ) then
call messages_write('Non-equilbrium PCM is not active in a time-independent run.')
call messages_new_line()
call messages_write('You set epsilon_infty /= epsilon_0, but epsilon_infty is not relevant for CalculationMode = gs.')
call messages_new_line()
call messages_write('By definition, the ground state is in equilibrium with the solvent.')
call messages_new_line()
call messages_write('Therefore, the only relevant dielectric constant is the static one.')
call messages_new_line()
call messages_write('Nevertheless, the dynamical PCM response matrix is evaluated for benchamarking purposes.')
call messages_new_line()
call messages_warning()
end if
if (hm%pcm%run_pcm) then
if ( (sys%mc%par_strategy /= P_STRATEGY_SERIAL).and.(sys%mc%par_strategy /= P_STRATEGY_STATES) ) then
......
......@@ -1203,6 +1203,8 @@ contains
logical :: kick_time
logical :: laser_present, kick_present
PUSH_SUB(v_ks_hartree)
ASSERT(associated(ks%hartree_solver))
......@@ -1230,13 +1232,14 @@ contains
if (hm%pcm%solute) &
call pcm_calc_pot_rs(hm%pcm, ks%gr%mesh, v_h = pot, time_present = ks%calc%time_present)
dt=hm%current_time/hm%pcm%iter
!> Local field effects due to the applied electrostatic potential representing the laser and the kick (if they were).
!! For the laser, the latter is only valid in the long-wavelength limit.
!! Static potentials are included in subroutine hamiltonian_epot_generate (module hamiltonian).
!! The sign convention for typical potentials and kick are different...
if( hm%pcm%localf .and. ks%calc%time_present ) then
if ( associated(hm%ep%lasers) .and. hm%ep%kick%delta_strength /= M_ZERO ) then !< external potential and kick
laser_present = epot_have_lasers( hm%ep )
kick_present = epot_have_kick( hm%ep )
if ( laser_present .and. kick_present ) then !< external potential and kick
SAFE_ALLOCATE(potx(1:ks%gr%mesh%np_part))
SAFE_ALLOCATE(kick(1:ks%gr%mesh%np_part))
SAFE_ALLOCATE(kick_real(1:ks%gr%mesh%np_part))
......@@ -1246,18 +1249,18 @@ contains
call laser_potential(hm%ep%lasers(ii), ks%gr%mesh, potx, ks%calc%time)
end do
kick_real = M_ZERO
kick_time =((hm%pcm%iter-1)*dt <= hm%ep%kick%time) .and. (hm%pcm%iter*dt > hm%ep%kick%time)
if ( hm%pcm%iter > 1 .and. kick_time ) then
kick_time = ((hm%pcm%iter-1)*hm%pcm%dt <= hm%ep%kick%time) .and. (hm%pcm%iter*hm%pcm%dt > hm%ep%kick%time)
if ( kick_time ) then
call kick_function_get(ks%gr%mesh, hm%ep%kick, kick, to_interpolate = .true.)
kick = hm%ep%kick%delta_strength * kick
kick_real = DREAL(kick)
end if
call pcm_calc_pot_rs(hm%pcm, ks%gr%mesh, v_ext = potx, kick = kick_real, time_present = ks%calc%time_present, &
call pcm_calc_pot_rs(hm%pcm, ks%gr%mesh, v_ext = potx, kick = -kick_real, time_present = ks%calc%time_present, &
kick_time = kick_time )
SAFE_DEALLOCATE_A(potx)
SAFE_DEALLOCATE_A(kick)
SAFE_DEALLOCATE_A(kick_real)
else if ( associated(hm%ep%lasers) .and. hm%ep%kick%delta_strength == M_ZERO ) then !< just external potential
else if ( laser_present .and. .not.kick_present ) then !< just external potential
SAFE_ALLOCATE(potx(1:ks%gr%mesh%np_part))
potx = M_ZERO
do ii = 1, hm%ep%no_lasers
......@@ -1265,18 +1268,18 @@ contains
end do
call pcm_calc_pot_rs(hm%pcm, ks%gr%mesh, v_ext = potx, time_present = ks%calc%time_present)
SAFE_DEALLOCATE_A(potx)
else if ( (.not.associated(hm%ep%lasers)) .and. hm%ep%kick%delta_strength /= M_ZERO ) then !< just kick
else if ( .not.laser_present .and. kick_present ) then !< just kick
SAFE_ALLOCATE(kick(1:ks%gr%mesh%np_part))
SAFE_ALLOCATE(kick_real(1:ks%gr%mesh%np_part))
kick = M_ZERO
kick_real = M_ZERO
kick_time =((hm%pcm%iter-1)*dt <= hm%ep%kick%time) .and. (hm%pcm%iter*dt > hm%ep%kick%time)
if ( hm%pcm%iter > 1 .and. kick_time ) then
kick_time =((hm%pcm%iter-1)*hm%pcm%dt <= hm%ep%kick%time) .and. (hm%pcm%iter*hm%pcm%dt > hm%ep%kick%time)
if ( kick_time ) then
call kick_function_get(ks%gr%mesh, hm%ep%kick, kick, to_interpolate = .true.)
kick = hm%ep%kick%delta_strength * kick
kick_real = DREAL(kick)
end if
call pcm_calc_pot_rs(hm%pcm, ks%gr%mesh, kick = kick_real, time_present = ks%calc%time_present, kick_time = kick_time )
call pcm_calc_pot_rs(hm%pcm, ks%gr%mesh, kick = -kick_real, time_present = ks%calc%time_present, kick_time = kick_time )
SAFE_DEALLOCATE_A(kick)
SAFE_DEALLOCATE_A(kick_real)
end if
......
......@@ -31,6 +31,7 @@ module spectrum_oct_m
use messages_oct_m
use minimizer_oct_m
use parser_oct_m
use pcm_oct_m
use profiling_oct_m
use string_oct_m
use types_oct_m
......@@ -601,6 +602,8 @@ contains
type(unit_system_t) :: file_units, ref_file_units
type(batch_t) :: dipoleb, sigmab
type(pcm_min_t) :: pcm
PUSH_SUB(spectrum_cross_section)
! This function gives us back the unit connected to the "multipoles" file, the header information,
......@@ -646,6 +649,13 @@ contains
call spectrum_read_dipole(ref_file, ref_dipole)
end if
! parsing and re-printing to output useful PCM data
call pcm_min_input_parsing_for_spectrum(pcm)
! adding the dipole generated by the PCM polarization charges due solute
if(pcm%localf) &
call spectrum_add_pcm_dipole(dipole, time_steps, dt, nspin)
! Now subtract the initial dipole.
if(present(ref_file)) then
dipole = dipole - ref_dipole
......@@ -662,16 +672,34 @@ contains
no_e = int(spectrum%max_energy / spectrum%energy_step)
SAFE_ALLOCATE(sigma(0:no_e, 1:3, 1:nspin))
if ( pcm%localf ) then
! for numerical reasons, we cannot add the difference (d(t)-d(t0)) of PCM dipoles here -- although it would look cleaner
! in the PCM local field case \sigma(\omega) \propto \Im{\alpha(\omega)\epsilon(\omega)} not just \propto \Im{\alpha(\omega)}
! since the dielectric function is complex as well, we need to compute both the real and imaginary part of the polarizability
call spectrum_times_pcm_epsilon(pcm, dipole, sigma, nspin, spectrum, istart, iend, kick%time, dt, no_e)
write(out_file,'(a57)') "Cross-section spectrum contains full local field effects."
else
call batch_init(dipoleb, 3, 1, nspin, dipole)
call batch_init(sigmab, 3, 1, nspin, sigma)
call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, kick%time, dt, dipoleb)
call spectrum_fourier_transform(spectrum%method, spectrum%transform, spectrum%noise, &
istart + 1, iend + 1, kick%time, dt, dipoleb, 1, no_e + 1, spectrum%energy_step, sigmab)
call batch_end(dipoleb)
call batch_end(sigmab)
end if
if( pcm%run_pcm ) &
call spectrum_over_pcm_refraction_index(pcm, sigma, nspin, spectrum, no_e)
call batch_init(dipoleb, 3, 1, nspin, dipole)
call batch_init(sigmab, 3, 1, nspin, sigma)
call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, kick%time, dt, dipoleb)
call spectrum_fourier_transform(spectrum%method, spectrum%transform, spectrum%noise, &
istart + 1, iend + 1, kick%time, dt, dipoleb, 1, no_e + 1, spectrum%energy_step, sigmab)
call batch_end(dipoleb)
call batch_end(sigmab)
SAFE_DEALLOCATE_A(dipole)
if(present(ref_file)) then
......@@ -787,6 +815,219 @@ contains
POP_SUB(spectrum_read_dipole)
end subroutine spectrum_read_dipole
! ---------------------------------------------------------
subroutine spectrum_add_pcm_dipole(dipole, time_steps, dt, nspin)
FLOAT, intent(inout) :: dipole(0:, :, :)
integer, intent(in) :: time_steps
FLOAT, intent(in) :: dt
integer, intent(in) :: nspin
type(pcm_t) :: pcm
FLOAT :: dipole_pcm(1:3)
integer :: ia, it, ii
! unit io variables
integer :: asc_unit_test
integer :: cavity_unit
integer :: asc_vs_t_unit, asc_vs_t_unit_check
integer :: dipole_vs_t_unit_check, dipole_vs_t_unit_check1
integer :: iocheck
integer :: aux_int
FLOAT :: aux_float, aux_float1, aux_vec(1:3)
character(len=23) :: asc_vs_t_unit_format
character(len=16) :: asc_vs_t_unit_format_tail
PUSH_SUB(spectrum_add_pcm_dipole)
! reading PCM cavity from standard output file in two steps
! first step - counting tesserae
asc_unit_test = io_open(PCM_DIR//'ASC_e.dat', action='read')
pcm%n_tesserae = 0
iocheck = 1
do while( iocheck >= 0 )
read(asc_unit_test,*,IOSTAT=iocheck) aux_vec(1:3), aux_float, aux_int
if( iocheck >= 0 ) pcm%n_tesserae = pcm%n_tesserae + 1
end do
call io_close(asc_unit_test)
! intermezzo - allocating PCM tessellation and polarization charges arrays
SAFE_ALLOCATE(pcm%tess(1:pcm%n_tesserae))
SAFE_ALLOCATE(pcm%q_e(1:pcm%n_tesserae))
SAFE_ALLOCATE(pcm%q_e_in(1:pcm%n_tesserae)) ! with auxiliary role
! second step - reading of PCM tessellation arrays from standard output file
! writing the cavity to debug-purpose file
asc_unit_test = io_open(PCM_DIR//'ASC_e.dat', action='read')
cavity_unit = io_open(PCM_DIR//'cavity_check.xyz', action='write')
write(cavity_unit,'(I3)') pcm%n_tesserae
write(cavity_unit,*)
do ia = 1, pcm%n_tesserae
read(asc_unit_test,*) pcm%tess(ia)%point(1:3), aux_float, aux_int
write(cavity_unit,'(A1,3(1X,F14.8))') 'H', pcm%tess(ia)%point(1:3)
end do
call io_close(asc_unit_test)
call io_close(cavity_unit)
write (asc_vs_t_unit_format_tail,'(I5,A11)') pcm%n_tesserae,'(1X,F14.8))'
write (asc_vs_t_unit_format,'(A)') '(F14.8,'//trim(adjustl(asc_vs_t_unit_format_tail))
! Now, summary: * read the time-dependent PCM polarization charges due to solute electrons, pcm%q_e
! * compute the real-time dipole generated by pcm%q_e, dipole_pcm
! * add it to the real-time molecular dipole
! * write the total dipole and its PCM contribution to debug-purpose files
! N.B. we assume nuclei fixed in time
! opening time-dependent PCM charges standard and debug-purpose file
asc_vs_t_unit = io_open(PCM_DIR//'ASC_e_vs_t.dat', action='read', form='formatted')
asc_vs_t_unit_check = io_open(PCM_DIR//'ASC_e_vs_t_check.dat', action='write', form='formatted')
! opening time-dependent PCM and total dipole debug-purpose files
dipole_vs_t_unit_check = io_open(PCM_DIR//'dipole_e_vs_t_check.dat', action='write', form='formatted')
dipole_vs_t_unit_check1 = io_open(PCM_DIR//'dipole_e_vs_t_check1.dat', action='write', form='formatted')
! reading PCM charges for the zeroth-iteration - not used - pcm%q_e_in is only auxiliary here
read(asc_vs_t_unit,trim(adjustl(asc_vs_t_unit_format))) aux_float1, ( pcm%q_e_in(ia) , ia=1,pcm%n_tesserae )
do it = 1, time_steps
! reading real-time PCM charges due to electrons per timestep
read(asc_vs_t_unit,trim(adjustl(asc_vs_t_unit_format))) aux_float, ( pcm%q_e(ia) , ia=1,pcm%n_tesserae )
! computing real-time PCM dipole per timestep
call pcm_dipole(dipole_pcm(1:3), -pcm%q_e(1:pcm%n_tesserae), pcm%tess, pcm%n_tesserae)
! adding PCM dipole to the molecular dipole per timestep
dipole(it, 1, 1:nspin) = dipole(it, 1, 1:nspin) + dipole_pcm(1)
dipole(it, 2, 1:nspin) = dipole(it, 2, 1:nspin) + dipole_pcm(2)
dipole(it, 3, 1:nspin) = dipole(it, 3, 1:nspin) + dipole_pcm(3)
! since we always have a kick for the optical spectrum in Octopus
! the first-iteration dipole should be equal to the zeroth-iteration one
! in any case, made them equal by hand
if ( it == 1 ) then
dipole(0, 1, 1:nspin) = dipole(1, 1, 1:nspin)
dipole(0, 2, 1:nspin) = dipole(1, 2, 1:nspin)
dipole(0, 3, 1:nspin) = dipole(1, 3, 1:nspin)
end if
! writing real-time PCM charges and dipole, and the total dipole for debug purposes
write(asc_vs_t_unit_check,trim(adjustl(asc_vs_t_unit_format))) aux_float, ( pcm%q_e(ia) , ia=1,pcm%n_tesserae )
write(dipole_vs_t_unit_check,'(F14.8,3(1X,F14.8))') aux_float, dipole_pcm
write(dipole_vs_t_unit_check1,'(F14.8,3(1X,F14.8))') aux_float, dipole(it,:,1)
end do
! closing PCM and debug files
call io_close(asc_vs_t_unit)
call io_close(asc_vs_t_unit_check)
call io_close(dipole_vs_t_unit_check)
call io_close(dipole_vs_t_unit_check1)
! deallocating PCM arrays
SAFE_DEALLOCATE_A(pcm%tess)
SAFE_DEALLOCATE_A(pcm%q_e)
SAFE_DEALLOCATE_A(pcm%q_e_in)
POP_SUB(spectrum_add_pcm_dipole)
end subroutine spectrum_add_pcm_dipole
! ---------------------------------------------------------
subroutine spectrum_times_pcm_epsilon(pcm, dipole, sigma, nspin, spectrum, istart, iend, kick_time, dt, no_e)
type(pcm_min_t) , intent(in) :: pcm
FLOAT, allocatable, intent(inout) :: sigma(:, :, :)
FLOAT, allocatable, intent(in) :: dipole(:, :, :)
integer, intent(in) :: nspin
type(spectrum_t), intent(in) :: spectrum !< check intent in
FLOAT, intent(in) :: kick_time
integer, intent(in) :: istart, iend
FLOAT, intent(in) :: dt
integer, intent(in) :: no_e
FLOAT, allocatable :: sigmap(:, :, :)
type(batch_t) :: dipoleb, sigmab
integer :: ie
CMPLX, allocatable :: eps(:)
PUSH_SUB(spectrum_times_pcm_epsilon)
! imaginary part of the polarizability
call batch_init(dipoleb, 3, 1, nspin, dipole)
call batch_init(sigmab, 3, 1, nspin, sigma)
call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, kick_time, dt, dipoleb)
call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_SIN, spectrum%noise, &
istart + 1, iend + 1, kick_time, dt, dipoleb, 1, no_e + 1, spectrum%energy_step, sigmab)
call batch_end(dipoleb)
call batch_end(sigmab)
! real part of the polarizability
SAFE_ALLOCATE(sigmap(0:no_e, 1:3, 1:nspin))
call batch_init(dipoleb, 3, 1, nspin, dipole)
call batch_init(sigmab, 3, 1, nspin, sigmap)
call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, kick_time, dt, dipoleb)
call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_COS, spectrum%noise, &
istart + 1, iend + 1, kick_time, dt, dipoleb, 1, no_e + 1, spectrum%energy_step, sigmab)
call batch_end(dipoleb)
call batch_end(sigmab)
SAFE_ALLOCATE(eps(0:no_e))
! multiplying by the dielectric function and taking the imaginary part of the product
do ie = 0, no_e
call pcm_eps(pcm, eps(ie), ie*spectrum%energy_step)
sigma(ie, 1:3, 1:nspin) = sigma(ie, 1:3, 1:nspin) * REAL(eps(ie), REAL_PRECISION) + sigmap(ie, 1:3, 1:nspin) *AIMAG(eps(ie))
end do
SAFE_DEALLOCATE_A(sigmap)
SAFE_DEALLOCATE_A(eps)
POP_SUB(spectrum_times_pcm_epsilon)
end subroutine spectrum_times_pcm_epsilon
! ---------------------------------------------------------
subroutine spectrum_over_pcm_refraction_index(pcm, sigma, nspin, spectrum, no_e)
type(pcm_min_t) , intent(in) :: pcm
FLOAT, allocatable, intent(inout) :: sigma(:, :, :)
integer, intent(in) :: nspin
type(spectrum_t), intent(in) :: spectrum !< check intent in
integer, intent(in) :: no_e
integer :: ie
CMPLX, allocatable :: eps(:)
PUSH_SUB(spectrum_over_pcm_refraction_index)
SAFE_ALLOCATE(eps(0:no_e))
! dividing by the refraction index - n(\omega)=\sqrt{\frac{|\epsilon(\omega)|+\Re[\epsilon(\omega)]}{2}}
do ie = 0, no_e
call pcm_eps(pcm, eps(ie), ie*spectrum%energy_step)
sigma(ie, 1:3, 1:nspin) = sigma(ie, 1:3, 1:nspin) / sqrt( CNST(0.5) * ( ABS(eps(ie)) + REAL(eps(ie), REAL_PRECISION) ) )
end do
SAFE_DEALLOCATE_A(eps)
POP_SUB(spectrum_over_pcm_refraction_index)
end subroutine spectrum_over_pcm_refraction_index
! ---------------------------------------------------------
subroutine spectrum_dipole_power(in_file, out_file, spectrum)
......
......@@ -470,7 +470,7 @@ contains