28 use,
intrinsic :: iso_fortran_env
43 integer :: iunit, ierr, ii, jj, iter, read_iter, ntime, nvel, ivel
44 integer :: istart, iend, energy_steps, out_file
45 real(real64),
allocatable :: time(:), velocities(:, :)
46 real(real64),
allocatable :: total_current(:, :), ftcurr(:, :, :), curr(:, :)
47 real(real64),
allocatable :: heat_current(:,:), ftheatcurr(:,:,:), heatcurr(:,:)
48 complex(real64),
allocatable :: invdielectric(:, :)
49 type(ions_t),
pointer :: ions
50 type(space_t) :: space
51 type(spectrum_t) :: spectrum
53 type(batch_t) :: currb, ftcurrb, heatcurrb, ftheatcurrb
54 real(real64) :: ww, curtime, deltat, integral(1:2), v0
55 character(len=MAX_PATH_LEN) :: ref_filename
56 integer :: ref_file, time_steps_ref, kk
57 real(real64),
allocatable :: velcm(:), vel0(:), current(:), current_ref(:, :)
58 real(real64) :: dt_ref, tt, start_time
59 integer :: ifreq, max_freq
60 integer :: skip, smearing
61 logical :: from_forces
62 character(len=120) :: header
63 real(real64) :: excess_charge, qtot, javerage
113 if (spectrum%end_time <
m_zero) spectrum%end_time = huge(spectrum%end_time)
119 qtot = -(ions%val_charge() + excess_charge)
121 if (from_forces)
then
123 call messages_write(
'Info: Reading coordinates from td.general/coordinates')
133 nvel = ions%natoms*ions%space%dim
135 safe_allocate(time(1:ntime))
136 safe_allocate(velocities(1:nvel, 1:ntime))
147 read(unit = iunit, iostat = ierr, fmt = *) read_iter, curtime, &
148 ((ions%pos(jj, ii), jj = 1, 3), ii = 1, ions%natoms), &
149 ((ions%vel(jj, ii), jj = 1, 3), ii = 1, ions%natoms)
153 if (ierr /= 0 .or. curtime >= spectrum%end_time)
then
159 assert(iter == read_iter + 1)
161 if (curtime >= spectrum%start_time .and. mod(iter, skip) == 0)
then
163 time(ntime) = curtime
165 do ii = 1, ions%natoms
166 do jj = 1, ions%space%dim
185 call messages_write(
'Info: Reading total current from td.general/total_current')
190 if (iunit /= -1)
then
199 safe_allocate(total_current(1:space%dim, 1:ntime))
200 safe_allocate(heat_current(1:space%dim, 1:ntime))
201 safe_allocate(time(1:ntime))
204 read(iunit, *) read_iter, time(iter), (total_current(ii, iter), ii = 1, space%dim)
215 ref_file =
io_open(trim(ref_filename)//
'/total_current', action=
'read', status=
'old')
218 if (time_steps_ref < ntime)
then
219 message(1) =
"The reference calculation does not contain enought time steps"
223 if (.not.
is_close(dt_ref, time(2)-time(1)))
then
224 message(1) =
"The time step of the reference calculation is different from the current calculation"
229 time_steps_ref = time_steps_ref + 1
230 safe_allocate(current_ref(1:space%dim, 1:time_steps_ref))
232 do ii = 1, time_steps_ref
233 read(ref_file, *) jj, tt, (current_ref(kk, ii), kk = 1, space%dim)
238 total_current(kk, iter) = total_current(kk, iter) - current_ref(kk, iter)
241 safe_deallocate_a(current_ref)
243 start_time = spectrum%start_time
254 javerage = javerage + total_current(kk, iter)
256 javerage = javerage/ntime
258 total_current(kk, iter) = total_current(kk, iter) - javerage
265 call messages_write(
"Cannot find the 'td.general/total_current' file.")
266 call messages_write(
"Conductivity will only be calculated from the forces")
270 safe_allocate(total_current(1:space%dim, 1:ntime))
271 safe_allocate(heat_current(1:space%dim, 1:ntime))
272 safe_allocate(time(1:ntime))
274 total_current(1:space%dim, 1:ntime) =
m_zero
280 if (iunit /= -1)
then
291 read(iunit, *) read_iter, time(iter), (heat_current(ii, iter), ii = 1, space%dim)
298 call messages_write(
"Cannot find the 'td.general/heat_current' file.")
299 call messages_write(
"Thermal conductivity will only be calculated from the forces")
302 heat_current(1:space%dim, 1:ntime) =
m_zero
308 safe_allocate(curr(ntime, 1:space%dim))
309 safe_allocate(heatcurr(ntime, 1:space%dim))
310 safe_allocate(vel0(1:space%dim))
311 safe_allocate(velcm(1:space%dim))
312 safe_allocate(current(1:space%dim))
319 if (from_forces)
then
323 do ii = 1, ions%natoms
325 vel0(jj) = vel0(jj) + velocities(ivel, iter)/real(ions%natoms, real64)
335 do ii = 1, ions%natoms
337 velcm(jj) = velcm(jj) + velocities(ivel, iter)/real(ions%natoms, real64)
338 current(jj) = current(jj) + &
339 ions%mass(ii)/ions%latt%rcell_volume*(velocities(ivel, iter) - vel0(jj))
344 integral(1) = integral(1) + deltat/vel0(1)*(vel0(1)*qtot/ions%latt%rcell_volume + current(1))
345 integral(2) = integral(2) + &
346 deltat/vel0(1)*(vel0(1)*qtot/ions%latt%rcell_volume - total_current(1, iter)/ions%latt%rcell_volume)
348 curr(iter, 1:space%dim) = vel0(1:space%dim)*qtot/ions%latt%rcell_volume + current(1:space%dim)
351 curr(iter, 1:space%dim) = total_current(1:space%dim, iter)/ions%latt%rcell_volume
352 heatcurr(iter,1:space%dim) = heat_current(1:space%dim, iter)/ions%latt%rcell_volume
355 if (from_forces)
write(iunit,*) iter, iter*deltat, curr(iter, 1:space%dim)
359 safe_deallocate_a(velcm)
360 safe_deallocate_a(current)
363 if (from_forces)
call io_close(iunit)
367 istart = max(1, istart)
370 safe_allocate(ftcurr(1:energy_steps, 1:space%dim, 1:2))
373 call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart, iend, spectrum%start_time, deltat, currb)
375 call batch_init(ftcurrb, 1, 1, space%dim, ftcurr(:, :, 1))
377 istart, iend, spectrum%start_time, deltat, currb, spectrum%min_energy, spectrum%max_energy, &
378 spectrum%energy_step, ftcurrb)
381 call batch_init(ftcurrb, 1, 1, space%dim, ftcurr(:, :, 2))
383 istart, iend, spectrum%start_time, deltat, currb, spectrum%min_energy, spectrum%max_energy, &
384 spectrum%energy_step, ftcurrb)
394 write(unit = iunit, iostat = ierr, fmt =
'(a)') &
395 '###########################################################################################################################'
396 write(unit = iunit, iostat = ierr, fmt =
'(8a)')
'# HEADER'
397 write(unit = iunit, iostat = ierr, fmt =
'(a,a,a)') &
399 write(unit = iunit, iostat = ierr, fmt =
'(a)') &
400 '###########################################################################################################################'
402 v0 = norm2(vel0(1:space%dim))
403 if (.not. from_forces .or. v0 < epsilon(v0)) v0 =
m_one
412 vel0(1:space%dim) = vel0(1:space%dim) /
p_c
413 v0 = norm2(vel0(1:space%dim))
417 do ifreq = 1, energy_steps
418 ww = spectrum%energy_step*(ifreq - 1) + spectrum%min_energy
420 transpose(ftcurr(ifreq, 1:space%dim, 1:2)/v0)
429 safe_allocate(invdielectric(1:space%dim, 1:energy_steps))
430 do ifreq = 1, energy_steps
431 ww = max((ifreq-1)*spectrum%energy_step + spectrum%min_energy,
m_epsilon)
433 invdielectric(1:space%dim, ifreq) = (vel0(1:space%dim) +
m_four *
m_pi * &
434 cmplx(ftcurr(ifreq, 1:space%dim, 2),-ftcurr(ifreq, 1:space%dim, 1), real64) / ww) / v0
439 select case (space%dim)
441 write(header,
'(3a15)')
'# energy',
'Re x',
'Im x'
443 write(header,
'(5a15)')
'# energy',
'Re x',
'Im x',
'Re y',
'Im y'
445 write(header,
'(7a15)')
'# energy',
'Re x',
'Im x',
'Re y',
'Im y',
'Re z',
'Im z'
447 write(out_file,
'(a)') trim(header)
448 do ifreq = 1, energy_steps
449 ww = (ifreq-1)*spectrum%energy_step + spectrum%min_energy
450 select case (space%dim)
452 write(out_file,
'(3e15.6)') ww, &
453 real(invdielectric(1, ifreq), real64), aimag(invdielectric(1, ifreq))
455 write(out_file,
'(5e15.6)') ww, &
456 real(invdielectric(1, ifreq), real64), aimag(invdielectric(1, ifreq)), &
457 real(invdielectric(2, ifreq), real64), aimag(invdielectric(2, ifreq))
459 write(out_file,
'(7e15.6)') ww, &
460 real(invdielectric(1, ifreq), real64), aimag(invdielectric(1, ifreq)), &
461 real(invdielectric(2, ifreq), real64), aimag(invdielectric(2, ifreq)), &
462 real(invdielectric(3, ifreq), real64), aimag(invdielectric(3, ifreq))
467 safe_deallocate_a(ftcurr)
468 safe_deallocate_a(invdielectric)
472 safe_allocate(ftheatcurr(1:energy_steps, 1:space%dim, 1:2))
476 call batch_init(heatcurrb, 1, 1, space%dim, heatcurr)
480 call batch_init(ftheatcurrb, 1, 1, space%dim, ftheatcurr(:, :, 1))
482 1, ntime,
m_zero, deltat, heatcurrb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftheatcurrb)
483 call ftheatcurrb%end()
485 call batch_init(ftheatcurrb, 1, 1, space%dim, ftheatcurr(:, :, 2))
487 1, ntime,
m_zero, deltat, heatcurrb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftheatcurrb)
488 call ftheatcurrb%end()
497 write(unit = iunit, iostat = ierr, fmt =
'(a)') &
498 '###########################################################################################################################'
499 write(unit = iunit, iostat = ierr, fmt =
'(8a)')
'# HEADER'
500 write(unit = iunit, iostat = ierr, fmt =
'(a,a,a)') &
502 write(unit = iunit, iostat = ierr, fmt =
'(a)') &
503 '###########################################################################################################################'
505 v0 = norm2(vel0(1:space%dim))
506 if (.not. from_forces .or. v0 < epsilon(v0)) v0 =
m_one
508 do ifreq = 1, energy_steps
509 ww = spectrum%energy_step*(ifreq - 1) + spectrum%min_energy
511 transpose(ftheatcurr(ifreq, 1:space%dim, 1:2)/v0)
517 safe_deallocate_a(ftheatcurr)
518 safe_deallocate_a(vel0)
522 safe_deallocate_p(ions)
524 safe_deallocate_a(total_current)
525 safe_deallocate_a(heat_current)
526 safe_deallocate_a(time)
initialize a batch with existing memory
This module implements batches of mesh functions.
subroutine, public getopt_init(ierr)
Initializes the getopt machinery. Must be called before attempting to parse the options....
subroutine, public getopt_end
subroutine, public global_end()
Finalise parser varinfo file, and MPI.
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_pi
some mathematical constants
type(mpi_comm), parameter, public serial_dummy_comm
Alias MPI_COMM_UNDEFINED for the specific use case of initialising Octopus utilities with no MPI supp...
subroutine, public init_octopus_globals(comm)
Initialise Octopus-specific global constants and files. This routine performs no initialisation calls...
real(real64), parameter, public m_epsilon
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
real(real64), parameter, public m_one
This module implements the underlying real-space grid.
subroutine, public io_init(defaults)
If the argument defaults is present and set to true, then the routine will not try to read anything f...
subroutine, public io_close(iunit, grp)
subroutine, public io_skip_header(iunit)
subroutine, public io_end()
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
This module is intended to contain "only mathematical" functions and procedures.
subroutine, public messages_end()
subroutine, public messages_init(output_dir)
subroutine, public messages_warning(no_lines, all_nodes, namespace)
subroutine, public messages_obsolete_variable(namespace, name, rep)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_experimental(name, namespace)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
type(namespace_t), public global_namespace
logical function, public parse_is_defined(namespace, name)
subroutine, public parser_init()
Initialise the Octopus parser.
subroutine, public parser_end()
End the Octopus parser.
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public profiling_end(namespace)
subroutine, public profiling_init(namespace)
Create profiling subdirectory.
integer, parameter, public smear_semiconductor
subroutine, public spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
subroutine, public spectrum_fourier_transform(method, transform, noise, time_start, time_end, t0, time_step, time_function, energy_start, energy_end, energy_step, energy_function)
Computes the sine, cosine, (or "exponential") Fourier transform of the real function given in the tim...
subroutine, public spectrum_init(spectrum, namespace, default_energy_step, default_max_energy)
subroutine, public spectrum_signal_damp(damp_type, damp_factor, time_start, time_end, t0, time_step, time_function)
integer, parameter, public spectrum_transform_cos
integer, parameter, public spectrum_transform_sin
subroutine, public spectrum_count_time_steps(namespace, iunit, time_steps, dt)
pure integer function, public spectrum_nenergy_steps(spectrum)
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
character(len=20) pure function, public units_abbrev(this)
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
subroutine, public unit_system_init(namespace)