48 use,
intrinsic :: iso_fortran_env
108 integer,
parameter,
public :: &
109 OUT_MULTIPOLES = 1, &
145 character(len=100) :: td_file_name(OUT_MAX) = [character(100) :: &
155 "magnetic_moments", &
170 "total_heat_current", &
171 "total_magnetization", &
173 "maxwell_dipole_field", &
174 "norm_wavefunctions"]
176 integer,
parameter :: &
177 OUT_DFTU_EFFECTIVE_U = 1, &
181 integer,
parameter :: &
182 OUT_MAXWELL_TOTAL_E_FIELD = 1, &
196 integer,
parameter,
public :: &
199 integer,
parameter :: &
200 MAXWELL_TOTAL_FIELD = 1, &
204 integer,
parameter :: &
210 type(c_ptr) :: handle
211 type(c_ptr),
allocatable :: mult_handles(:)
213 integer :: hand_start
215 logical ::
write = .false.
216 logical :: resolve_states = .false.
228 real(real64) :: lmm_r
231 integer :: n_excited_states
233 integer :: compute_interval
239 subroutine td_write_kick(outp, namespace, space, mesh, kick, ions, iter)
243 class(
mesh_t),
intent(in) :: mesh
244 type(
kick_t),
intent(in) :: kick
245 type(
ions_t),
intent(in) :: ions
246 integer,
intent(in) :: iter
248 complex(real64),
allocatable :: kick_function(:)
249 character(len=256) :: filename
254 write(filename,
'(a,i7.7)')
"td.", iter
255 if (outp%what(option__output__delta_perturbation))
then
256 safe_allocate(kick_function(1:mesh%np))
258 call zio_function_output(outp%how(option__output__delta_perturbation), filename,
"kick_function", namespace, &
259 space, mesh, kick_function(:),
units_out%energy, err, pos=ions%pos, atoms=ions%atom)
260 safe_deallocate_a(kick_function)
276 subroutine td_write_init(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, ks, ions_move, &
277 with_gauge_field, kick, iter, max_iter, dt, mc)
281 type(
output_t),
intent(inout) :: outp
282 type(
grid_t),
intent(in) :: gr
285 type(
ions_t),
intent(in) :: ions
287 type(
v_ks_t),
intent(inout) :: ks
288 logical,
intent(in) :: ions_move
289 logical,
intent(in) :: with_gauge_field
290 type(
kick_t),
intent(in) :: kick
291 integer,
intent(in) :: iter
292 integer,
intent(in) :: max_iter
293 real(real64),
intent(in) :: dt
297 integer :: ierr, first, ii, ist, jj, flags, iout, default, ifile
302 character(len=MAX_PATH_LEN) :: filename
304 logical :: resolve_states
305 logical,
allocatable :: skip(:)
445 output_options = .false.
446 output_options(out_multipoles) = .
true.
461 writ%out(iout)%write = output_options(iout)
476 if (space%is_periodic() .and. writ%out(
out_angular)%write)
then
487 if (gr%np /= gr%np_global)
then
488 message(1) =
"TDOutput option td_kpoint_occup and td_floquet do not work with domain parallelization"
495 'Forces for systems periodic in 1D are not currently implemented and options that output the forces are not allowed.')
499 if (writ%out(
out_kp_proj)%write .and. hm%kpoints%nik_skip == 0)
then
500 message(1) =
"TDOutput option td_kpoint_occup only work with zero-weight k-points at the moment."
514 call parse_variable(namespace,
'TDOutputResolveStates', .false., resolve_states)
515 if (.not. writ%out(out_multipoles)%write)
then
516 write(
message(1),
'(a)')
"TDOutputResolveStates works only for TDOutput = multipoles."
529 if (writ%lmax < 0)
then
530 write(
message(1),
'(a,i6,a)')
"Input: '", writ%lmax,
"' is not a valid TDMultipoleLmax."
531 message(2) =
'(Must be TDMultipoleLmax >= 0 )'
537 if ((writ%out(
out_acc)%write) .and. ions_move)
then
538 message(1) =
'If harmonic spectrum is to be calculated, atoms should not be allowed to move.'
542 if ((writ%out(
out_q)%write) .and. .not.(ks%has_photons))
then
543 message(1) =
'If q(t) is to be calculated, you need to allow for photon modes.'
548 .or. hm%mxll%add_electric_dip))
then
549 message(1) =
'If the dipolar field has to be written, MaxwellCouplingMode has to be'
550 message(2) =
'"lenght_gauge_dipole" or "velocity_gauge_dipole" and at least one Maxwell system'
551 message(3) =
'must be present.'
555 rmin = ions%min_distance()
565 message(1) =
"Option TDOutput = populations is not implemented for parallel in states."
577 safe_deallocate_a(writ%gs_st%node)
585 writ%gs_st%st_end = writ%gs_st%nst
587 message(1) =
"Unable to read states information."
591 writ%gs_st%st_start = 1
607 call parse_variable(namespace,
'TDProjStateStart', 1, writ%gs_st%st_start)
609 if (st%parallel_in_states .and. writ%out(
out_proj)%write .and. writ%gs_st%st_start > 1)
then
615 writ%gs_st%lnst = writ%gs_st%st_end - writ%gs_st%st_start + 1
619 writ%gs_st%parallel_in_states = .false.
622 safe_allocate(writ%gs_st%occ(1:writ%gs_st%nst, 1:writ%gs_st%nik))
623 safe_allocate(writ%gs_st%eigenval(1:writ%gs_st%nst, 1:writ%gs_st%nik))
627 safe_allocate(writ%gs_st%node(1:writ%gs_st%nst))
628 writ%gs_st%node(:) = 0
630 writ%gs_st%eigenval = huge(writ%gs_st%eigenval)
632 if (writ%gs_st%d%ispin ==
spinors)
then
633 safe_deallocate_a(writ%gs_st%spin)
634 safe_allocate(writ%gs_st%spin(1:3, 1:writ%gs_st%nst, 1:writ%gs_st%nik))
637 safe_allocate(skip(1:writ%gs_st%nst))
639 skip(1:writ%gs_st%st_start-1) = .
true.
643 safe_deallocate_a(skip)
648 if (ierr /= 0 .and. ierr /= (writ%gs_st%st_end-writ%gs_st%st_start+1)*writ%gs_st%nik &
649 *writ%gs_st%d%dim*writ%gs_st%mpi_grp%size)
then
650 message(1) =
"Unable to read wavefunctions for TDOutput."
695 if (
parse_block(namespace,
'TDExcitedStatesToProject', blk) == 0)
then
697 safe_allocate(writ%excited_st(1:writ%n_excited_states))
698 do ist = 1, writ%n_excited_states
703 writ%n_excited_states = 0
717 call parse_variable(namespace,
'TDOutputComputeInterval', 50, writ%compute_interval)
718 if (writ%compute_interval < 0)
then
719 message(1) =
"TDOutputComputeInterval must be >= 0."
735 call io_mkdir(
'td.general', namespace)
746 do ifile = 1, out_max
750 if (writ%out(ifile)%write)
then
754 trim(
io_workpath(
"td.general/"//trim(td_file_name(ifile)), namespace)))
762 if (writ%out(out_multipoles)%write .and. .not. resolve_states)
then
765 trim(
io_workpath(
"td.general/multipoles", namespace)))
769 select case (kick%qkick_mode)
771 write(filename,
'(a)')
'td.general/ftchd.sin'
773 write(filename,
'(a)')
'td.general/ftchd.cos'
775 write(filename,
'(a, SP, I0.3, a, I0.3)')
'td.general/ftchd.l', kick%qbessel_l,
'_m', kick%qbessel_m
777 write(filename,
'(a)')
'td.general/ftchd'
789 call io_rm(
"td.general/laser", namespace=namespace)
806 if (writ%out(out_multipoles)%write .and. resolve_states)
then
808 writ%out(out_multipoles)%hand_start = st%st_start
809 writ%out(out_multipoles)%hand_end = st%st_end
810 writ%out(out_multipoles)%resolve_states = .
true.
811 writ%out(out_multipoles)%mpi_grp = gr%mpi_grp
813 safe_allocate(writ%out(out_multipoles)%mult_handles(st%st_start:st%st_end))
816 do ist = st%st_start, st%st_end
817 write(filename,
'(a,i4.4)')
'td.general/multipoles-ist', ist
830 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
831 calc_eigenval=.false., time = iter*dt, calc_energy = .false.)
834 if (writ%out(
out_n_ex)%write .and. writ%compute_interval > 0)
then
835 call io_mkdir(outp%iter_dir, namespace)
838 if (all(outp%how == 0) .and. writ%out(
out_n_ex)%write)
then
860 if (hm%lda_u_level ==
dft_u_acbn0) default = default + 2**(out_dftu_effective_u - 1)
868 writ%out_dftu(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
872 if (writ%out_dftu(out_dftu_effective_u)%write)
then
875 trim(
io_workpath(
"td.general/effectiveU", namespace)))
893 if (writ%out(iout)%write)
then
895 if (writ%out(iout)%resolve_states)
then
896 do ist = writ%out(iout)%hand_start, writ%out(iout)%hand_end
899 safe_deallocate_a(writ%out(iout)%mult_handles)
909 if (writ%out_dftu(iout)%write)
call write_iter_end(writ%out_dftu(iout)%handle)
915 do ist = 1, writ%n_excited_states
918 writ%n_excited_states = 0
931 subroutine td_write_iter(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, kick, ks, dt, iter, mc, recalculate_gs)
934 class(
space_t),
intent(in) :: space
936 type(
grid_t),
intent(in) :: gr
939 type(
ions_t),
intent(inout) :: ions
941 type(
kick_t),
intent(in) :: kick
942 type(
v_ks_t),
intent(in) :: ks
943 real(real64),
intent(in) :: dt
944 integer,
intent(in) :: iter
946 logical,
intent(in) :: recalculate_gs
954 if (writ%out(out_multipoles)%write)
then
955 call td_write_multipole(writ%out(out_multipoles), space, gr, ions, st, writ%lmax, kick, iter)
978 if (writ%out(
out_proj)%write .and. mod(iter, writ%compute_interval) == 0)
then
987 if (writ%out(
out_kp_proj)%write .and. mod(iter, writ%compute_interval) == 0)
then
993 ions%pos, ions%vel, ions%tot_force, iter)
998 ions%pos, ions%vel, ions%tot_force, iter, 1)
1003 ions%pos, ions%vel, ions%tot_force, iter, 2)
1008 ions%pos, ions%vel, ions%tot_force, iter, 3)
1019 if (writ%out(
out_acc)%write)
then
1020 call td_write_acc(writ%out(
out_acc)%handle, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
1023 if (writ%out(
out_vel)%write)
then
1036 if(
associated(gfield))
then
1062 if (writ%out(
out_n_ex)%write .and. mod(iter, writ%compute_interval) == 0)
then
1064 if (recalculate_gs)
then
1067 ierr, label =
': Houston states for TDOutput')
1071 call td_write_n_ex(writ%out(
out_n_ex)%handle, outp, namespace, gr, hm%kpoints, st, writ%gs_st, iter)
1079 if (writ%out_dftu(out_dftu_effective_u)%write)
then
1083 if (writ%out(
out_q)%write .and. ks%has_photons)
then
1087 if (writ%out(
out_mxll_field)%write .and. hm%mxll%calc_field_dip)
then
1104 do iout = 1, out_max
1106 if (writ%out(iout)%write)
then
1108 if (writ%out(iout)%resolve_states)
then
1109 do ii = writ%out(iout)%hand_start, writ%out(iout)%hand_end
1121 if (writ%out_dftu(iout)%write)
call write_iter_flush(writ%out_dftu(iout)%handle)
1130 subroutine td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, dt)
1133 type(
grid_t),
intent(in) :: gr
1136 type(
v_ks_t),
intent(inout) :: ks
1138 type(
ions_t),
intent(in) :: ions
1140 integer,
intent(in) :: iter
1141 real(real64),
optional,
intent(in) :: dt
1143 character(len=256) :: filename
1149 if (st%modelmbparticles%nparticle > 0)
then
1154 write(filename,
'(a,a,i7.7)') trim(outp%iter_dir),
"td.", iter
1156 call output_all(outp, namespace, space, filename, gr, ions, iter, st, hm, ks)
1158 call output_modelmb(outp, namespace, space, filename, gr, ions, iter, st)
1160 if (
present(dt))
then
1161 call output_scalar_pot(outp, namespace, space, filename, gr, ions, ext_partners, iter*dt)
1163 if (iter == 0)
call output_scalar_pot(outp, namespace, space, filename, gr, ions, ext_partners)
1172 type(c_ptr),
intent(inout) ::
out_spin
1173 class(
mesh_t),
intent(in) :: mesh
1175 integer,
intent(in) :: iter
1177 character(len=130) :: aux
1178 real(real64) :: spin(3)
1194 if (st%d%ispin ==
spinors)
then
1195 write(aux,
'(a2,18x)')
'Sx'
1197 write(aux,
'(a2,18x)')
'Sy'
1200 write(aux,
'(a2,18x)')
'Sz'
1208 select case (st%d%ispin)
1225 type(
grid_t),
intent(in) :: gr
1227 type(
ions_t),
intent(in) :: ions
1228 real(real64),
intent(in) :: lmm_r
1229 integer,
intent(in) :: iter
1232 character(len=50) :: aux
1233 real(real64),
allocatable :: lmm(:,:)
1238 safe_allocate(lmm(1:3, 1:ions%natoms))
1248 do ia = 1, ions%natoms
1249 if (st%d%ispin ==
spinors)
then
1250 write(aux,
'(a2,i2.2,16x)')
'mx', ia
1252 write(aux,
'(a2,i2.2,16x)')
'my', ia
1255 write(aux,
'(a2,i2.2,16x)')
'mz', ia
1264 do ia = 1, ions%natoms
1265 select case (st%d%ispin)
1273 safe_deallocate_a(lmm)
1281 type(c_ptr),
intent(inout) :: out_magnets
1282 class(
mesh_t),
intent(in) :: mesh
1284 type(
kick_t),
intent(in) :: kick
1285 integer,
intent(in) :: iter
1287 complex(real64),
allocatable :: tm(:,:)
1292 safe_allocate(tm(1:6,1:kick%nqvec))
1294 do iq = 1, kick%nqvec
1324 do iq = 1, kick%nqvec
1333 safe_deallocate_a(tm)
1347 subroutine td_write_angular(out_angular, namespace, space, gr, ions, hm, st, kick, iter)
1348 type(c_ptr),
intent(inout) :: out_angular
1350 class(
space_t),
intent(in) :: space
1351 type(
grid_t),
intent(in) :: gr
1352 type(
ions_t),
intent(inout) :: ions
1355 type(
kick_t),
intent(in) :: kick
1356 integer,
intent(in) :: iter
1359 character(len=130) :: aux
1360 real(real64) :: angular(3)
1367 call angular_momentum%setup_dir(idir)
1370 m_two*real(angular_momentum%zstates_elec_expectation_value(namespace, space, gr, hm, st), real64)
1372 safe_deallocate_p(angular_momentum)
1379 write(aux,
'(a15,i2)')
'# nspin ', st%d%nspin
1387 write(aux,
'(a4,18x)')
'<Lx>'
1389 write(aux,
'(a4,18x)')
'<Ly>'
1391 write(aux,
'(a4,18x)')
'<Lz>'
1420 class(
space_t),
intent(in) :: space
1421 type(
grid_t),
intent(in) :: gr
1422 type(
ions_t),
intent(in) :: ions
1424 integer,
intent(in) :: lmax
1425 type(
kick_t),
intent(in) :: kick
1426 integer,
intent(in) :: iter
1429 real(real64),
allocatable :: rho(:,:)
1433 if (out_multip%resolve_states)
then
1434 safe_allocate(rho(1:gr%np_part, 1:st%d%nspin))
1437 do ist = st%st_start, st%st_end
1439 call td_write_multipole_r(out_multip%mult_handles(ist), space, gr, ions, st, lmax, kick, rho, iter, &
1440 mpi_grp = out_multip%mpi_grp)
1443 safe_deallocate_a(rho)
1446 if (
allocated(st%frozen_rho))
then
1447 safe_allocate(rho(1:gr%np, 1:st%d%nspin))
1448 call lalg_copy(gr%np, st%d%nspin, st%rho, rho)
1453 safe_deallocate_a(rho)
1465 subroutine td_write_multipole_r(out_multip, space, mesh, ions, st, lmax, kick, rho, iter, mpi_grp)
1466 type(c_ptr),
intent(inout) :: out_multip
1467 class(
space_t),
intent(in) :: space
1468 class(
mesh_t),
intent(in) :: mesh
1469 type(
ions_t),
intent(in) :: ions
1471 integer,
intent(in) :: lmax
1472 type(
kick_t),
intent(in) :: kick
1473 real(real64),
intent(in) :: rho(:,:)
1474 integer,
intent(in) :: iter
1475 type(
mpi_grp_t),
optional,
intent(in) :: mpi_grp
1478 integer :: is, idir, ll, mm, add_lm
1479 character(len=120) :: aux
1480 real(real64) :: ionic_dipole(ions%space%dim)
1481 real(real64),
allocatable :: multipole(:,:)
1487 assert(.not. (lmax > 1 .and. space%dim > 3))
1490 if (
present(mpi_grp)) mpi_grp_ = mpi_grp
1495 write(aux,
'(a15,i2)')
'# nspin ', st%d%nspin
1499 write(aux,
'(a15,i2)')
'# lmax ', lmax
1507 do is = 1, st%d%nspin
1508 write(aux,
'(a18,i1,a1)')
'Electronic charge(', is,
')'
1511 do idir = 1, space%dim
1512 write(aux,
'(4a1,i1,a1)')
'<',
index2axis(idir),
'>',
'(', is,
')'
1518 write(aux,
'(a2,i2,a4,i2,a2,i1,a1)')
'l=', ll,
', m=', mm,
' (', is,
')'
1529 do is = 1, st%d%nspin
1532 do idir = 1, space%dim
1548 if (space%dim > 3 .and. lmax == 1)
then
1550 safe_allocate(multipole(1:space%dim+1, 1:st%d%nspin))
1552 safe_allocate(multipole(1:(lmax + 1)**2, 1:st%d%nspin))
1556 do is = 1, st%d%nspin
1561 ionic_dipole = ions%dipole()
1562 do is = 1, st%d%nspin
1563 multipole(2:space%dim+1, is) = -ionic_dipole(1:space%dim)/st%d%nspin - multipole(2:space%dim+1, is)
1569 do is = 1, st%d%nspin
1572 do idir = 1, space%dim
1576 add_lm = space%dim + 2
1587 safe_deallocate_a(multipole)
1592 subroutine td_write_ftchd(out_ftchd, space, mesh, st, kick, iter)
1593 type(c_ptr),
intent(inout) :: out_ftchd
1594 class(
space_t),
intent(in) :: space
1595 class(
mesh_t),
intent(in) :: mesh
1597 type(
kick_t),
intent(in) :: kick
1598 integer,
intent(in) :: iter
1600 integer :: is, ip, idir
1601 character(len=120) :: aux, aux2
1602 real(real64) :: ftchd_bessel
1603 complex(real64) :: ftchd
1605 real(real64),
allocatable :: integrand_bessel(:)
1606 complex(real64),
allocatable :: integrand(:)
1613 write(aux,
'(a15, i2)')
'# qkickmode ', kick%qkick_mode
1618 write(aux,
'(a15, i0.3, 1x, i0.3)')
'# ll, mm ', kick%qbessel_l, kick%qbessel_m
1624 write(aux,
'(a15, f9.6)')
'# qlength ', kick%qlength
1626 write(aux,
'(a15)')
'# qvector '
1627 do idir = 1, space%dim
1628 write(aux2,
'(f9.5)') kick%qvector(idir,1)
1629 aux = trim(aux) // trim(aux2)
1635 write(aux,
'(a15,f18.12)')
'# kick strength', kick%delta_strength
1641 write(aux,
'(a17)')
'int(j_l*Y_lm*rho)'
1643 write(aux,
'(a12)')
'Real, Imag'
1660 safe_allocate(integrand(1:mesh%np))
1662 do is = 1, st%d%nspin
1664 integrand(ip) = integrand(ip) + st%rho(ip, is) *
exp(-
m_zi*sum(mesh%x(ip, 1:space%dim)*kick%qvector(1:space%dim, 1)))
1668 safe_deallocate_a(integrand)
1671 safe_allocate(integrand_bessel(1:mesh%np))
1672 integrand_bessel =
m_zero
1673 do is = 1, st%d%nspin
1675 call ylmr_real(mesh%x(ip, 1:3), kick%qbessel_l, kick%qbessel_m, ylm)
1676 integrand_bessel(ip) = integrand_bessel(ip) + st%rho(ip, is) * &
1677 loct_sph_bessel(kick%qbessel_l, kick%qlength*norm2(mesh%x(ip, :)))*ylm
1681 safe_deallocate_a(integrand_bessel)
1700 type(c_ptr),
intent(inout) :: out_temperature
1701 type(
ions_t),
intent(in) :: ions
1702 integer,
intent(in) :: iter
1737 type(c_ptr),
intent(inout) :: out_populations
1739 class(
space_t),
intent(in) :: space
1740 class(
mesh_t),
intent(in) :: mesh
1743 real(real64),
intent(in) :: dt
1744 integer,
intent(in) :: iter
1747 character(len=6) :: excited_name
1748 complex(real64) :: gsp
1749 complex(real64),
allocatable :: excited_state_p(:)
1750 complex(real64),
allocatable :: dotprodmatrix(:, :, :)
1755 safe_allocate(dotprodmatrix(1:writ%gs_st%nst, 1:st%nst, 1:st%nik))
1760 assert(.not. space%is_periodic())
1765 if (writ%n_excited_states > 0)
then
1766 safe_allocate(excited_state_p(1:writ%n_excited_states))
1767 do ist = 1, writ%n_excited_states
1768 excited_state_p(ist) =
zstates_elec_mpdotp(namespace, mesh, writ%excited_st(ist), st, dotprodmatrix)
1780 do ist = 1, writ%n_excited_states
1781 write(excited_name,
'(a2,i3,a1)')
'P(', ist,
')'
1800 do ist = 1, writ%n_excited_states
1807 if (writ%n_excited_states > 0)
then
1808 safe_deallocate_a(excited_state_p)
1810 safe_deallocate_a(dotprodmatrix)
1816 subroutine td_write_acc(out_acc, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
1817 type(c_ptr),
intent(inout) :: out_acc
1819 class(
space_t),
intent(in) :: space
1820 type(
grid_t),
intent(in) :: gr
1821 type(
ions_t),
intent(inout) :: ions
1825 real(real64),
intent(in) :: dt
1826 integer,
intent(in) :: iter
1829 character(len=7) :: aux
1830 real(real64) :: acc(space%dim)
1839 do idim = 1, space%dim
1840 write(aux,
'(a4,i1,a1)')
'Acc(', idim,
')'
1848 do idim = 1, space%dim
1855 call td_calc_tacc(namespace, space, gr, ions, ext_partners, st, hm, acc, dt*iter)
1868 subroutine td_write_vel(out_vel, gr, st, space, kpoints, iter)
1869 type(c_ptr),
intent(inout) :: out_vel
1870 type(
grid_t),
intent(in) :: gr
1872 type(
space_t),
intent(in) :: space
1874 integer,
intent(in) :: iter
1877 character(len=7) :: aux
1878 real(real64) :: vel(space%dim)
1887 do idim = 1, space%dim
1888 write(aux,
'(a4,i1,a1)')
'Vel(', idim,
')'
1896 do idim = 1, space%dim
1918 type(c_ptr),
intent(inout) :: out_laser
1919 class(
space_t),
intent(in) :: space
1920 type(
lasers_t),
intent(inout) :: lasers
1921 real(real64),
intent(in) :: dt
1922 integer,
intent(in) :: iter
1925 real(real64) :: field(space%dim)
1926 real(real64) :: ndfield(space%dim)
1927 character(len=80) :: aux
1943 do il = 1, lasers%no_lasers
1946 do idir = 1, space%dim
1947 write(aux,
'(a,i1,a)')
'E(', idir,
')'
1951 do idir = 1, space%dim
1952 write(aux,
'(a,i1,a)')
'B(', idir,
')'
1956 do idir = 1, space%dim
1957 write(aux,
'(a,i1,a)')
'A(', idir,
')'
1961 write(aux,
'(a,i1,a)')
'e(t)'
1967 do idir = 1, space%dim
1968 write(aux,
'(a,i1,a)')
'A^M(', idir,
')'
1980 do il = 1, lasers%no_lasers
1984 do idir = 1, space%dim
1989 do idir = 1, space%dim
2000 do idir = 1, space%dim
2013 do il = 1, lasers%no_lasers
2015 call laser_field(lasers%lasers(il), field(1:space%dim), iter*dt)
2040 type(c_ptr),
intent(inout) :: out_energy
2042 integer,
intent(in) :: iter
2043 real(real64),
intent(in) :: ke
2047 integer :: n_columns
2070 if (hm%pcm%run_pcm)
then
2072 n_columns = n_columns + 1
2077 n_columns = n_columns + 1
2087 do ii = 1, n_columns
2110 hm%energy%int_nn_pcm + hm%energy%int_ne_pcm), 1)
2123 type(c_ptr),
intent(inout) :: out_eigs
2125 integer,
intent(in) :: iter
2128 character(len=68) :: buf
2141 write(buf,
'(a15,i2)')
'# nst ', st%nst
2145 write(buf,
'(a15,i2)')
'# nspin ', st%d%nspin
2151 do is = 1, st%d%kpt%nglobal
2153 write(buf,
'(a,i4)')
'Eigenvalue ',ii
2162 do is = 1, st%d%kpt%nglobal
2172 do is = 1, st%d%kpt%nglobal
2184 type(c_ptr),
intent(inout) :: out_ionch
2185 class(
mesh_t),
intent(in) :: mesh
2187 integer,
intent(in) :: iter
2189 integer :: ii, ist, Nch, ik, idim
2190 character(len=68) :: buf
2191 real(real64),
allocatable :: ch(:), occ(:)
2192 real(real64),
allocatable :: occbuf(:)
2197 nch = st%nst * st%d%kpt%nglobal * st%d%dim
2198 safe_allocate(ch(0: nch))
2199 safe_allocate(occ(0: nch))
2205 do idim = 1, st%d%dim
2206 if (st%st_start <= ist .and. ist <= st%st_end .and. &
2207 st%d%kpt%start <= ik .and. ik <= st%d%kpt%end)
then
2208 occ(ii) = st%occ(ist, ik)
2216 if (st%parallel_in_states)
then
2217 safe_allocate(occbuf(0: nch))
2219 call st%mpi_grp%allreduce(occ(0), occbuf(0), nch+1, mpi_double_precision, mpi_sum)
2221 safe_deallocate_a(occbuf)
2229 safe_deallocate_a(ch)
2242 if (occ(ii)>
m_zero .or. ii == 0)
then
2243 write(buf,
'(a,f4.1,a)')
'Pion(',occ(ii)*ii,
'+, t)'
2253 if (occ(ii)>
m_zero .or. ii == 0)
then
2263 if (occ(ii)>
m_zero .or. ii == 0)
then
2269 safe_deallocate_a(ch)
2270 safe_deallocate_a(occ)
2277 type(c_ptr),
intent(inout) :: out_proj
2278 class(
space_t),
intent(in) :: space
2279 class(
mesh_t),
intent(in) :: mesh
2280 type(
ions_t),
intent(in) :: ions
2283 type(
kick_t),
intent(in) :: kick
2284 integer,
intent(in) :: iter
2286 complex(real64),
allocatable :: projections(:,:,:)
2287 character(len=80) :: aux
2288 integer :: ik, ist, uist, idir
2296 write(aux,
'(a15,i2)')
'# nspin ', st%d%nspin
2305 write(aux,
'(a,i8)')
"# nik ", st%nik
2309 write(aux,
'(a,2i8)')
"# st ", gs_st%st_start, st%nst
2313 write(aux,
'(a,2i8)')
"# ust ", gs_st%st_start, gs_st%st_end
2319 do ist = gs_st%st_start, st%nst
2327 do ist = gs_st%st_start, st%nst
2328 do uist = gs_st%st_start, gs_st%st_end
2329 write(aux,
'(i4,a,i4)') ist,
' -> ', uist
2340 if (.not. space%is_periodic())
then
2342 safe_allocate(projections(1:st%nst, gs_st%st_start:gs_st%st_end, 1:st%nik))
2343 do idir = 1, space%dim
2349 write(aux,
'(a,i1,a)')
"<i|x_", idir,
"|a>"
2353 do ist = gs_st%st_start, st%st_end
2354 do uist = gs_st%st_start, gs_st%st_end
2364 safe_deallocate_a(projections)
2374 safe_allocate(projections(1:st%nst, gs_st%st_start:gs_st%st_end, 1:st%nik))
2375 projections(:,:,:) =
m_z0
2381 do ist = gs_st%st_start, st%nst
2382 do uist = gs_st%st_start, gs_st%st_end
2391 safe_deallocate_a(projections)
2397 integer,
intent(in) :: dir
2399 integer :: uist, ist, ik, idim
2400 real(real64) :: n_dip(space%dim)
2401 complex(real64),
allocatable :: xpsi(:,:)
2402 complex(real64),
allocatable :: psi(:, :), gspsi(:, :)
2406 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
2407 safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
2408 safe_allocate(xpsi(1:mesh%np, 1:st%d%dim))
2410 do ik = st%d%kpt%start, st%d%kpt%end
2411 do ist = st%st_start, st%st_end
2413 do uist = gs_st%st_start, gs_st%st_end
2416 do idim = 1, st%d%dim
2417 xpsi(1:mesh%np, idim) = mesh%x(1:mesh%np, dir)*gspsi(1:mesh%np, idim)
2419 projections(ist, uist, ik) = -
zmf_dotp(mesh, st%d%dim, psi, xpsi, reduce = .false.)
2425 safe_deallocate_a(xpsi)
2426 safe_deallocate_a(gspsi)
2427 safe_deallocate_a(psi)
2432 n_dip = ions%dipole()
2434 do ist = gs_st%st_start, st%nst
2435 do uist = gs_st%st_start, gs_st%st_end
2436 projections(ist, uist, ik) = projections(ist, uist, ik) - n_dip(dir)
2452 subroutine td_write_n_ex(out_nex, outp, namespace, mesh, kpoints, st, gs_st, iter)
2453 type(c_ptr),
intent(inout) :: out_nex
2456 class(
mesh_t),
intent(in) :: mesh
2460 integer,
intent(in) :: iter
2462 complex(real64),
allocatable :: projections(:,:)
2463 character(len=80) :: aux, dir
2464 integer :: ik, ikpt, ist, uist, err
2465 real(real64) :: Nex, weight
2467 real(real64),
allocatable :: Nex_kpt(:)
2476 write(aux,
'(a15,i2)')
'# nspin ', st%d%nspin
2483 write(aux,
'(a,i8)')
"# nik ", st%nik
2487 write(aux,
'(a,2i8)')
"# st ", gs_st%st_start, st%nst
2491 write(aux,
'(a,2i8)')
"# ust ", gs_st%st_start, gs_st%st_end
2508 if (.not. st%parallel_in_states)
then
2509 do ik = st%d%kpt%start, st%d%kpt%end
2510 do ist = 1, gs_st%nst
2511 if (gs_st%occ(ist, ik) >
m_min_occ) gs_nst = ist
2518 safe_allocate(projections(1:gs_nst, 1:st%nst))
2520 safe_allocate(nex_kpt(1:st%nik))
2522 do ik = st%d%kpt%start, st%d%kpt%end
2523 ikpt = st%d%get_kpoint_index(ik)
2526 weight = st%kweights(ik) * gs_st%occ(ist, ik)/ st%smear%el_per_state
2527 do uist = st%st_start, st%st_end
2528 nex_kpt(ikpt) = nex_kpt(ikpt) - weight * st%occ(uist, ik) * abs(projections(ist, uist))**2
2532 nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*
m_half*st%kweights(ik)
2534 nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*st%kweights(ik)
2538 if (st%parallel_in_states .or. st%d%kpt%parallel)
then
2550 write(dir,
'(a,a,i7.7)') trim(outp%iter_dir),
"td.", iter
2553 + outp%how(option__output__density_kpt), dir,
"n_excited_el_kpt", namespace, &
2557 safe_deallocate_a(projections)
2558 safe_deallocate_a(nex_kpt)
2570 class(
mesh_t),
intent(in) :: mesh
2573 complex(real64),
intent(inout) :: projections(1:st%nst, gs_st%st_start:gs_st%nst, 1:st%nik)
2575 integer :: uist, ist, ik
2576 complex(real64),
allocatable :: psi(:, :), gspsi(:, :)
2579 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
2580 safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
2582 projections(:,:,:) =
m_zero
2584 do ik = st%d%kpt%start, st%d%kpt%end
2585 do ist = st%st_start, st%st_end
2587 do uist = gs_st%st_start, gs_st%nst
2589 projections(ist, uist, ik) =
zmf_dotp(mesh, st%d%dim, psi, gspsi, reduce = .false.)
2594 safe_deallocate_a(psi)
2595 safe_deallocate_a(gspsi)
2604 class(
mesh_t),
intent(in) :: mesh
2609 integer,
intent(in) :: iter
2611 complex(real64),
allocatable :: proj(:,:), psi(:,:,:), gs_psi(:,:,:), temp_state(:,:)
2612 character(len=80) :: filename1, filename2
2613 integer :: ik,ist, jst, file, idim, nk_proj
2617 write(filename1,
'(I10)') iter
2618 filename1 =
'td.general/projections_iter_'//trim(adjustl(filename1))
2621 safe_allocate(proj(1:gs_st%nst, 1:gs_st%nst))
2622 safe_allocate(psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2623 safe_allocate(gs_psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2624 safe_allocate(temp_state(1:mesh%np,1:gs_st%d%dim))
2630 nk_proj = kpoints%nik_skip
2632 do ik = kpoints%reduced%npoints-nk_proj+1, kpoints%reduced%npoints
2634 psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)=
m_zero
2635 gs_psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)=
m_zero
2638 write(filename2,
'(I10)') ik
2639 filename2 = trim(adjustl(filename1))//
'_ik_'//trim(adjustl(filename2))
2640 file =
io_open(filename2, namespace, action=
'write')
2643 do ist=gs_st%st_start,gs_st%st_end
2646 do idim = 1,gs_st%d%dim
2647 psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2650 do idim = 1,gs_st%d%dim
2651 gs_psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2660 assert(mesh%np_global*gs_st%d%dim < huge(0_int32))
2661 proj(1:gs_st%nst, 1:gs_st%nst) =
m_zero
2666 i8_to_i4(mesh%np_global*gs_st%d%dim), &
2667 cmplx(mesh%volume_element,
m_zero, real64) , &
2669 ubound(psi, dim = 1), &
2671 ubound(gs_psi, dim = 1), &
2674 ubound(proj, dim = 1))
2678 do ist = 1, gs_st%nst
2679 do jst = 1, gs_st%nst
2680 write(file,
'(I3,1x,I3,1x,e13.6,1x,e13.6,2x)') ist, jst, proj(ist,jst)
2688 safe_deallocate_a(proj)
2689 safe_deallocate_a(psi)
2690 safe_deallocate_a(gs_psi)
2691 safe_deallocate_a(temp_state)
2697 subroutine td_write_floquet(namespace, space, hm, ext_partners, mesh, st, iter)
2698 type(namespace_t),
intent(in) :: namespace
2699 class(space_t),
intent(in) :: space
2700 type(hamiltonian_elec_t),
intent(inout) :: hm
2701 type(partner_list_t),
intent(in) :: ext_partners
2702 class(mesh_t),
intent(in) :: mesh
2703 type(states_elec_t),
intent(inout) :: st
2704 integer,
intent(in) :: iter
2706 complex(real64),
allocatable :: hmss(:,:), psi(:,:,:), hpsi(:,:,:), temp_state1(:,:)
2707 complex(real64),
allocatable :: HFloquet(:,:,:), HFloq_eff(:,:), temp(:,:)
2708 real(real64),
allocatable :: eigenval(:), bands(:,:)
2709 character(len=80) :: filename
2710 integer :: it, nT, ik, ist, in, im, file, idim, nik, ik_count
2711 integer :: Forder, Fdim, m0, n0, n1, nst, ii, jj, lim_nst
2712 logical :: downfolding
2713 type(states_elec_t) :: hm_st
2715 real(real64) :: dt, Tcycle, omega
2719 downfolding = .false.
2722 if (.not. iter == 0)
then
2730 assert(mesh%np == mesh%np_global)
2733 call states_elec_copy(hm_st, st)
2744 call parse_variable(namespace,
'TDFloquetFrequency', m_zero, omega, units_inp%energy)
2745 call messages_print_var_value(
'Frequency used for Floquet analysis', omega, namespace=namespace)
2746 if (abs(omega) <= m_epsilon)
then
2747 message(1) =
"Please give a non-zero value for TDFloquetFrequency"
2748 call messages_fatal(1, namespace=namespace)
2752 tcycle = m_two * m_pi / omega
2762 call parse_variable(namespace,
'TDFloquetSample',20 ,nt)
2763 call messages_print_var_value(
'Number of Floquet time-sampling points', nt, namespace=namespace)
2764 dt = tcycle/real(nt, real64)
2773 call parse_variable(namespace,
'TDFloquetDimension',-1,forder)
2774 if (forder .ge. 0)
then
2775 call messages_print_var_value(
'Order of multiphoton Floquet-Hamiltonian', forder, namespace=namespace)
2777 fdim = 2 * forder + 1
2779 message(1) =
'Floquet-Hamiltonian is downfolded'
2780 call messages_info(1, namespace=namespace)
2781 downfolding = .
true.
2786 dt = tcycle/real(nt, real64)
2789 nik = hm%kpoints%nik_skip
2791 safe_allocate(hmss(1:nst,1:nst))
2792 safe_allocate( psi(1:nst,1:st%d%dim,1:mesh%np))
2793 safe_allocate(hpsi(1:nst,1:st%d%dim,1:mesh%np))
2794 safe_allocate(temp_state1(1:mesh%np,1:st%d%dim))
2802 safe_allocate(hfloquet(1:nik,1:nst*fdim, 1:nst*fdim))
2803 hfloquet(1:nik,1:nst*fdim, 1:nst*fdim) = m_zero
2808 call hm%update(mesh, namespace, space, ext_partners, time=tcycle+it*dt)
2810 call zhamiltonian_elec_apply_all(hm, namespace, mesh, st, hm_st)
2815 do ik = hm%kpoints%reduced%npoints-nik+1, hm%kpoints%reduced%npoints
2816 ik_count = ik_count + 1
2818 psi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
2819 hpsi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
2821 do ist = st%st_start, st%st_end
2822 if (state_kpt_is_local(st, ist, ik))
then
2823 call states_elec_get_state(st, mesh, ist, ik,temp_state1)
2824 do idim = 1, st%d%dim
2825 psi(ist, idim, 1:mesh%np) = temp_state1(1:mesh%np,idim)
2827 call states_elec_get_state(hm_st, mesh, ist, ik,temp_state1)
2828 do idim = 1, st%d%dim
2829 hpsi(ist, idim,1:mesh%np) = temp_state1(1:mesh%np,idim)
2833 call comm_allreduce(mpi_world, psi)
2834 call comm_allreduce(mpi_world, hpsi)
2835 assert(mesh%np_global*st%d%dim < huge(0_int32))
2836 hmss(1:nst,1:nst) = m_zero
2841 i8_to_i4(mesh%np_global*st%d%dim), &
2842 cmplx(mesh%volume_element, m_zero, real64) , &
2844 ubound(hpsi, dim = 1), &
2846 ubound(psi, dim = 1), &
2849 ubound(hmss, dim = 1))
2851 hmss(1:nst,1:nst) = conjg(hmss(1:nst,1:nst))
2854 do in = -forder, forder
2855 do im = -forder, forder
2856 ii = (in+forder) * nst
2857 jj = (im+forder) * nst
2858 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) = &
2859 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) + hmss(1:nst, 1:nst) *
exp(-(in-im)*m_zi*omega*it*dt)
2863 hfloquet(ik_count, ii+ist, ii+ist) = hfloquet(ik_count, ii+ist, ii+ist) + in*omega
2872 hfloquet(:,:,:) = m_one/nt*hfloquet(:,:,:)
2875 if (downfolding)
then
2877 safe_allocate(hfloq_eff(1:nst,1:nst))
2878 safe_allocate(eigenval(1:nst))
2879 safe_allocate(bands(1:nik,1:nst))
2881 hfloq_eff(1:nst,1:nst) = m_zero
2887 hfloq_eff(1:nst, 1:nst) = hfloquet(ik, n0+1:n0+nst, m0+1:m0+nst) + &
2888 m_one/omega*(matmul(hfloquet(ik, 1:nst, m0+1:m0+nst), hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst))- &
2889 matmul(hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst), hfloquet(ik, 1:nst, m0+1:m0+nst)))
2891 call lalg_eigensolve(nst, hfloq_eff, eigenval)
2892 bands(ik,1:nst) = eigenval(1:nst)
2894 safe_deallocate_a(hfloq_eff)
2897 safe_allocate(eigenval(1:nst*fdim))
2898 safe_allocate(bands(1:nik,1:nst*fdim))
2899 safe_allocate(temp(1:nst*fdim, 1:nst*fdim))
2902 temp(1:nst*fdim, 1:nst*fdim) = hfloquet(ik, 1:nst*fdim, 1:nst*fdim)
2903 call lalg_eigensolve(nst*fdim, temp, eigenval)
2904 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2909 if (downfolding)
then
2911 filename =
"downfolded_floquet_bands"
2914 filename =
"floquet_bands"
2917 if (mpi_world%rank == 0)
then
2919 file = io_open(filename, namespace, action =
'write')
2922 write(file,
'(e13.6, 1x)',advance=
'no') bands(ik,ist)
2929 if (.not. downfolding)
then
2932 bands(1:nik, 1:nst*fdim) = m_zero
2934 temp(1:nst*fdim,1:nst*fdim) = m_zero
2937 temp(ii+1:ii+nst, ii+1:ii+nst) = hfloquet(ik, ii+1:ii+nst, ii+1:ii+nst)
2939 call lalg_eigensolve(nst*fdim, temp, eigenval)
2940 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2943 if (mpi_world%rank == 0)
then
2944 filename =
'trivial_floquet_bands'
2945 file = io_open(filename, namespace, action =
'write')
2948 write(file,
'(e13.6, 1x)', advance=
'no') bands(ik,ist)
2957 call hm%update(mesh, namespace, space, ext_partners, time=m_zero)
2959 safe_deallocate_a(hmss)
2960 safe_deallocate_a(psi)
2961 safe_deallocate_a(hpsi)
2962 safe_deallocate_a(temp_state1)
2963 safe_deallocate_a(hfloquet)
2964 safe_deallocate_a(eigenval)
2965 safe_deallocate_a(bands)
2966 safe_deallocate_a(temp)
2974 type(c_ptr),
intent(inout) :: out_total_current
2975 class(space_t),
intent(in) :: space
2976 class(mesh_t),
intent(in) :: mesh
2977 type(states_elec_t),
intent(in) :: st
2978 integer,
intent(in) :: iter
2980 integer :: idir, ispin
2981 character(len=50) :: aux
2982 real(real64) :: total_current(space%dim), abs_current(space%dim)
2986 if (mpi_grp_is_root(mpi_world) .and. iter == 0)
then
2987 call td_write_print_header_init(out_total_current)
2990 call write_iter_header_start(out_total_current)
2992 do idir = 1, space%dim
2993 write(aux,
'(a2,a1,a1)')
'I(', index2axis(idir),
')'
2994 call write_iter_header(out_total_current, aux)
2997 do idir = 1, space%dim
2998 write(aux,
'(a10,a1,a1)')
'IntAbs(j)(', index2axis(idir),
')'
2999 call write_iter_header(out_total_current, aux)
3002 do ispin = 1, st%d%nspin
3003 do idir = 1, space%dim
3004 write(aux,
'(a4,i1,a1,a1,a1)')
'I-sp', ispin,
'(', index2axis(idir),
')'
3005 call write_iter_header(out_total_current, aux)
3009 call write_iter_nl(out_total_current)
3011 call td_write_print_header_end(out_total_current)
3014 assert(
allocated(st%current))
3016 if (mpi_grp_is_root(mpi_world))
then
3017 call write_iter_start(out_total_current)
3020 total_current = 0.0_real64
3021 do idir = 1, space%dim
3022 do ispin = 1, st%d%spin_channels
3023 total_current(idir) = total_current(idir) + dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.)
3025 total_current(idir) = units_from_atomic(units_out%length/units_out%time, total_current(idir))
3027 if (mesh%parallel_in_domains)
then
3028 call mesh%allreduce(total_current, dim = space%dim)
3031 abs_current = 0.0_real64
3032 do idir = 1, space%dim
3033 do ispin = 1, st%d%spin_channels
3034 abs_current(idir) = abs_current(idir) + dmf_integrate(mesh, abs(st%current(:, idir, ispin)), reduce = .false.)
3036 abs_current(idir) = units_from_atomic(units_out%length/units_out%time, abs_current(idir))
3038 if (mesh%parallel_in_domains)
then
3039 call mesh%allreduce(abs_current, dim = space%dim)
3042 if (mpi_grp_is_root(mpi_world))
then
3043 call write_iter_double(out_total_current, total_current, space%dim)
3044 call write_iter_double(out_total_current, abs_current, space%dim)
3047 do ispin = 1, st%d%nspin
3048 total_current = 0.0_real64
3049 do idir = 1, space%dim
3050 total_current(idir) = units_from_atomic(units_out%length/units_out%time, &
3051 dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.))
3053 if (mesh%parallel_in_domains)
then
3054 call mesh%allreduce(total_current, dim = space%dim)
3056 if (mpi_grp_is_root(mpi_world))
then
3057 call write_iter_double(out_total_current, total_current, space%dim)
3061 if (mpi_grp_is_root(mpi_world))
then
3062 call write_iter_nl(out_total_current)
3071 type(c_ptr),
intent(inout) :: write_obj
3072 class(space_t),
intent(in) :: space
3073 type(hamiltonian_elec_t),
intent(inout) :: hm
3074 type(grid_t),
intent(in) :: gr
3075 type(states_elec_t),
intent(in) :: st
3076 integer,
intent(in) :: iter
3078 integer :: idir, ispin
3079 character(len=50) :: aux
3080 real(real64),
allocatable :: heat_current(:, :, :)
3081 real(real64) :: total_current(space%dim)
3085 if (mpi_grp_is_root(mpi_world) .and. iter == 0)
then
3086 call td_write_print_header_init(write_obj)
3089 call write_iter_header_start(write_obj)
3091 do idir = 1, space%dim
3092 write(aux,
'(a2,i1,a1)')
'Jh(', idir,
')'
3093 call write_iter_header(write_obj, aux)
3096 call write_iter_nl(write_obj)
3098 call td_write_print_header_end(write_obj)
3101 safe_allocate(heat_current(1:gr%np, 1:space%dim, 1:st%d%nspin))
3103 call current_heat_calculate(space, gr%der, hm, st, heat_current)
3105 if (mpi_grp_is_root(mpi_world))
call write_iter_start(write_obj)
3107 total_current = 0.0_real64
3108 do idir = 1, space%dim
3109 do ispin = 1, st%d%spin_channels
3110 total_current(idir) = total_current(idir) + dmf_integrate(gr, heat_current(:, idir, ispin))
3112 total_current(idir) = units_from_atomic(units_out%energy*units_out%length/units_out%time, total_current(idir))
3115 safe_deallocate_a(heat_current)
3117 if (mpi_grp_is_root(mpi_world))
call write_iter_double(write_obj, total_current, space%dim)
3119 if (mpi_grp_is_root(mpi_world))
call write_iter_nl(write_obj)
3127 type(c_ptr),
intent(inout) :: out_partial_charges
3128 class(mesh_t),
intent(in) :: mesh
3129 type(states_elec_t),
intent(in) :: st
3130 type(ions_t),
intent(in) :: ions
3131 integer,
intent(in) :: iter
3134 character(len=50) :: aux
3135 real(real64),
allocatable :: hirshfeld_charges(:)
3139 safe_allocate(hirshfeld_charges(1:ions%natoms))
3141 call partial_charges_calculate(mesh, st, ions, hirshfeld_charges)
3143 if (mpi_grp_is_root(mpi_world))
then
3147 call td_write_print_header_init(out_partial_charges)
3150 call write_iter_header_start(out_partial_charges)
3152 do idir = 1, ions%natoms
3153 write(aux,
'(a13,i3,a1)')
'hirshfeld(atom=', idir,
')'
3154 call write_iter_header(out_partial_charges, aux)
3157 call write_iter_nl(out_partial_charges)
3159 call td_write_print_header_end(out_partial_charges)
3162 call write_iter_start(out_partial_charges)
3164 call write_iter_double(out_partial_charges, hirshfeld_charges, ions%natoms)
3166 call write_iter_nl(out_partial_charges)
3169 safe_deallocate_a(hirshfeld_charges)
3175 subroutine td_write_q(out_q, space, ks, iter)
3176 type(c_ptr),
intent(inout) :: out_q
3177 class(space_t),
intent(in) :: space
3178 type(v_ks_t),
intent(in) :: ks
3179 integer,
intent(in) :: iter
3182 character(len=50) :: aux
3186 if (mpi_grp_is_root(mpi_world))
then
3188 call td_write_print_header_init(out_q)
3189 call write_iter_header_start(out_q)
3190 do ii = 1, ks%pt%nmodes
3191 write(aux,
'(a1,i3,a3)')
'q', ii,
'(t)'
3192 call write_iter_header(out_q, aux)
3194 do ii = 1, ks%pt%nmodes
3195 write(aux,
'(a1,i3,a3)')
'p', ii,
'(t)'
3196 call write_iter_header(out_q, aux)
3198 do ii = 1, space%dim
3199 write(aux,
'(a3,i3,a3)')
'f_pt', ii,
'(t)'
3200 call write_iter_header(out_q, aux)
3202 call write_iter_nl(out_q)
3203 call td_write_print_header_end(out_q)
3206 call write_iter_start(out_q)
3207 call write_iter_double(out_q, ks%pt_mx%pt_q, ks%pt%nmodes)
3208 call write_iter_double(out_q, ks%pt_mx%pt_p, ks%pt%nmodes)
3209 call write_iter_double(out_q, ks%pt_mx%fmf, space%dim)
3210 call write_iter_nl(out_q)
3219 type(c_ptr),
intent(inout) :: out_mxll
3220 class(space_t),
intent(in) :: space
3221 type(hamiltonian_elec_t),
intent(in) :: hm
3222 real(real64),
intent(in) :: dt
3223 integer,
intent(in) :: iter
3226 real(real64) :: field(space%dim)
3227 character(len=80) :: aux
3228 character(len=1) :: field_char
3232 if (.not. mpi_grp_is_root(mpi_world))
then
3238 call td_write_print_header_init(out_mxll)
3240 write(aux,
'(a7,e20.12,3a)')
'# dt = ', units_from_atomic(units_out%time, dt), &
3241 " [", trim(units_abbrev(units_out%time)),
"]"
3242 call write_iter_string(out_mxll, aux)
3243 call write_iter_nl(out_mxll)
3245 call write_iter_header_start(out_mxll)
3246 select case (hm%mxll%coupling_mode)
3247 case (length_gauge_dipole, multipolar_expansion)
3248 if (hm%mxll%add_electric_dip) field_char =
'E'
3249 if (hm%mxll%add_magnetic_dip) field_char =
'B'
3250 do idir = 1, space%dim
3251 write(aux,
'(a,i1,a)') field_char //
'(', idir,
')'
3252 call write_iter_header(out_mxll, aux)
3254 case (velocity_gauge_dipole)
3255 do idir = 1, space%dim
3256 write(aux,
'(a,i1,a)')
'A(', idir,
')'
3257 call write_iter_header(out_mxll, aux)
3260 call write_iter_nl(out_mxll)
3262 call write_iter_string(out_mxll,
'#[Iter n.]')
3263 call write_iter_header(out_mxll,
'[' // trim(units_abbrev(units_out%time)) //
']')
3267 select case (hm%mxll%coupling_mode)
3268 case (length_gauge_dipole, multipolar_expansion)
3269 if (hm%mxll%add_electric_dip) aux =
'[' // trim(units_abbrev(units_out%force)) //
']'
3270 if (hm%mxll%add_magnetic_dip) aux =
'[' // trim(units_abbrev(unit_one/units_out%length**2)) //
']'
3271 do idir = 1, space%dim
3272 call write_iter_header(out_mxll, aux)
3274 case (velocity_gauge_dipole)
3275 aux =
'[' // trim(units_abbrev(units_out%energy)) //
']'
3276 do idir = 1, space%dim
3277 call write_iter_header(out_mxll, aux)
3280 call write_iter_nl(out_mxll)
3281 call td_write_print_header_end(out_mxll)
3284 call write_iter_start(out_mxll)
3287 select case (hm%mxll%coupling_mode)
3288 case (length_gauge_dipole, multipolar_expansion)
3289 if (hm%mxll%add_electric_dip) field = units_from_atomic(units_out%force, hm%mxll%e_field_dip)
3290 if (hm%mxll%add_magnetic_dip) field = units_from_atomic(unit_one/units_out%length**2, hm%mxll%b_field_dip)
3291 call write_iter_double(out_mxll, field, space%dim)
3292 case (velocity_gauge_dipole)
3293 field = units_from_atomic(units_out%energy, hm%mxll%vec_pot_dip)
3294 call write_iter_double(out_mxll, field, space%dim)
3296 call write_iter_nl(out_mxll)
3304 type(c_ptr),
intent(inout) :: out_coords
3305 type(lda_u_t),
intent(in) :: lda_u
3306 integer,
intent(in) :: iter
3309 character(len=50) :: aux
3311 if (.not. mpi_grp_is_root(mpi_world))
return
3316 call td_write_print_header_init(out_coords)
3319 call write_iter_header_start(out_coords)
3321 do ios = 1, lda_u%norbsets
3322 write(aux,
'(a2,i3,a1)')
'Ueff(', ios,
')'
3323 call write_iter_header(out_coords, aux)
3326 do ios = 1, lda_u%norbsets
3327 write(aux,
'(a2,i3,a1)')
'U(', ios,
')'
3328 call write_iter_header(out_coords, aux)
3331 do ios = 1, lda_u%norbsets
3332 write(aux,
'(a2,i3,a1)')
'J(', ios,
')'
3333 call write_iter_header(out_coords, aux)
3336 if (lda_u%intersite)
then
3337 do ios = 1, lda_u%norbsets
3338 do inn = 1, lda_u%orbsets(ios)%nneighbors
3339 write(aux,
'(a2,i3,a1,i3,a1)')
'V(', ios,
'-', inn,
')'
3340 call write_iter_header(out_coords, aux)
3346 call write_iter_nl(out_coords)
3349 call write_iter_string(out_coords,
'#[Iter n.]')
3350 call write_iter_header(out_coords,
'[' // trim(units_abbrev(units_out%time)) //
']')
3351 call write_iter_string(out_coords, &
3352 'Effective U ' // trim(units_abbrev(units_out%energy)) // &
3353 ', U in '// trim(units_abbrev(units_out%energy)) // &
3354 ', J in ' // trim(units_abbrev(units_out%energy)))
3355 call write_iter_nl(out_coords)
3357 call td_write_print_header_end(out_coords)
3360 call write_iter_start(out_coords)
3362 do ios = 1, lda_u%norbsets
3363 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3364 lda_u%orbsets(ios)%Ueff), 1)
3367 do ios = 1, lda_u%norbsets
3368 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3369 lda_u%orbsets(ios)%Ubar), 1)
3372 do ios = 1, lda_u%norbsets
3373 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3374 lda_u%orbsets(ios)%Jbar), 1)
3377 if (lda_u%intersite)
then
3378 do ios = 1, lda_u%norbsets
3379 do inn = 1, lda_u%orbsets(ios)%nneighbors
3380 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3381 lda_u%orbsets(ios)%V_ij(inn,0)), 1)
3386 call write_iter_nl(out_coords)
3393 type(c_ptr),
intent(inout) :: file_handle
3394 type(grid_t),
intent(in) :: grid
3395 type(kpoints_t),
intent(in) :: kpoints
3396 type(states_elec_t),
intent(in) :: st
3397 integer,
intent(in) :: iter
3399 integer :: ik_ispin, ist
3400 character(len=7) :: nkpt_str, nst_str
3401 character(len=7) :: ik_str, ist_str
3402 real(real64),
allocatable :: norm_ks(:, :)
3403 real(real64) :: n_electrons
3407 safe_allocate(norm_ks(1:st%nst, 1:st%nik))
3408 call states_elec_calc_norms(grid, kpoints, st, norm_ks)
3410 if (mpi_grp_is_root(mpi_world))
then
3413 call td_write_print_header_init(file_handle)
3416 write(nkpt_str,
'(I7)') st%nik
3417 write(nst_str,
'(I7)') st%nst
3418 call write_iter_string(file_handle,
'# Dimensions. (nstates, nkpt * nspin):')
3419 call write_iter_string(file_handle, trim(adjustl(nst_str)) //
' ' // trim(adjustl(nkpt_str)))
3420 call write_iter_nl(file_handle)
3423 call write_iter_string(file_handle,
'# Norm ordering: (istate, ikpoint_spin)')
3424 call write_iter_nl(file_handle)
3427 call write_iter_header_start(file_handle)
3428 call write_iter_header(file_handle,
'N_electrons')
3429 do ik_ispin = 1, st%nik
3431 write(ik_str,
'(I7)') ik_ispin
3432 write(ist_str,
'(I7)') ist
3433 call write_iter_header(file_handle, &
3434 'Norm (' // trim(ist_str) //
',' // trim(ik_str) //
')')
3437 call write_iter_nl(file_handle)
3438 call td_write_print_header_end(file_handle)
3441 n_electrons = sum(st%occ * norm_ks**2)
3444 call write_iter_start(file_handle)
3445 call write_iter_double(file_handle, n_electrons, 1)
3446 do ik_ispin = 1, st%nik
3447 call write_iter_double(file_handle, norm_ks(:, ik_ispin),
size(norm_ks, 1))
3449 call write_iter_nl(file_handle)
3453 safe_deallocate_a(norm_ks)
3462 type(namespace_t),
intent(in) :: namespace
3463 integer,
intent(in) :: iter
3464 real(real64),
intent(in) :: dt
3466 integer :: default, flags, iout, first
3521 call parse_variable(namespace,
'MaxwellTDOutput', default, flags)
3523 if (.not. varinfo_valid_option(
'MaxwellTDOutput', flags, is_flag = .
true.))
then
3524 call messages_input_error(namespace,
'MaxwellTDOutput')
3528 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3529 if (writ%out(iout)%write)
then
3530 writ%out(iout + 1)%write = .
true.
3531 writ%out(iout + 2)%write = .
true.
3536 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3545 call io_mkdir(
'td.general', namespace)
3550 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_x", namespace)))
3552 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_y", namespace)))
3554 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_z", namespace)))
3560 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_x", namespace)))
3562 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_y", namespace)))
3564 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_z", namespace)))
3570 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_x", namespace)))
3572 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_y", namespace)))
3574 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_z", namespace)))
3580 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_x", namespace)))
3582 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_y", namespace)))
3584 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_z", namespace)))
3590 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_x", namespace)))
3592 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_y", namespace)))
3594 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_z", namespace)))
3600 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_x", namespace)))
3602 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_y", namespace)))
3604 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_z", namespace)))
3610 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/maxwell_energy", namespace)))
3615 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-x", namespace)))
3620 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-y", namespace)))
3625 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-z", namespace)))
3630 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-x", namespace)))
3635 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-y", namespace)))
3640 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-z", namespace)))
3656 if (writ%out(iout)%write)
call write_iter_end(writ%out(iout)%handle)
3664 subroutine td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
3666 class(space_t),
intent(in) :: space
3667 type(grid_t),
intent(inout) :: gr
3668 type(states_mxll_t),
intent(inout) :: st
3669 type(hamiltonian_mxll_t),
intent(inout) :: hm
3670 type(helmholtz_decomposition_t),
intent(inout) :: helmholtz
3671 real(real64),
intent(in) :: dt
3672 integer,
intent(in) :: iter
3673 type(namespace_t),
intent(in) :: namespace
3678 call profiling_in(
"TD_WRITE_ITER_MAXWELL")
3681 call helmholtz%get_trans_field(namespace, st%rs_state_trans, total_field=st%rs_state)
3682 call get_rs_state_at_point(st%selected_points_rs_state_trans(:,:), st%rs_state_trans, &
3683 st%selected_points_coordinate(:,:), st, gr)
3686 hm%energy%energy_trans = m_zero
3690 call helmholtz%get_long_field(namespace, st%rs_state_long, total_field=st%rs_state)
3691 call get_rs_state_at_point(st%selected_points_rs_state_long(:,:), st%rs_state_long, &
3692 st%selected_points_coordinate(:,:), st, gr)
3695 hm%energy%energy_long = m_zero
3787 call profiling_out(
"TD_WRITE_ITER_MAXWELL")
3796 type(hamiltonian_mxll_t),
intent(in) :: hm
3797 integer,
intent(in) :: iter
3801 integer :: n_columns
3803 if (.not. mpi_grp_is_root(mpi_world))
return
3828 call write_iter_header(
out_maxwell_energy,
'[' // trim(units_abbrev(units_out%time)) //
']')
3830 do ii = 1, n_columns
3831 call write_iter_header(
out_maxwell_energy,
'[' // trim(units_abbrev(units_out%energy)) //
']')
3839 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy), 1)
3840 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%e_energy), 1)
3841 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%b_energy), 1)
3843 hm%energy%energy+hm%energy%boundaries), 1)
3844 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%boundaries), 1)
3845 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_trans), 1)
3846 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_long), 1)
3847 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_plane_waves), 1)
3856 type(c_ptr),
intent(inout) :: out_field_surf
3857 type(states_mxll_t),
intent(in) :: st
3858 integer,
intent(in) :: dim
3859 integer,
intent(in) :: iter
3863 integer :: n_columns
3865 if (.not. mpi_grp_is_root(mpi_world))
return
3872 call td_write_print_header_init(out_field_surf)
3875 call write_iter_header_start(out_field_surf)
3876 call write_iter_header(out_field_surf,
'- x direction')
3877 call write_iter_header(out_field_surf,
'+ x direction')
3878 call write_iter_header(out_field_surf,
'- y direction')
3879 call write_iter_header(out_field_surf,
'+ y direction')
3880 call write_iter_header(out_field_surf,
'- z direction')
3881 call write_iter_header(out_field_surf,
'+ z direction')
3882 call write_iter_header(out_field_surf,
'- x dir. p. w.')
3883 call write_iter_header(out_field_surf,
'+ x dir. p. w.')
3884 call write_iter_header(out_field_surf,
'- y dir. p. w.')
3885 call write_iter_header(out_field_surf,
'+ y dir. p. w.')
3886 call write_iter_header(out_field_surf,
'- z dir. p. w.')
3887 call write_iter_header(out_field_surf,
'+ z dir. p. w.')
3889 call write_iter_nl(out_field_surf)
3892 call write_iter_string(out_field_surf,
'#[Iter n.]')
3893 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%time)) //
']')
3895 do ii = 1, n_columns
3896 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%energy/units_out%length)) //
']')
3898 call write_iter_nl(out_field_surf)
3900 call td_write_print_header_end(out_field_surf)
3903 call write_iter_start(out_field_surf)
3904 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3905 st%electric_field_box_surface(1,1,dim)), 1)
3906 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3907 st%electric_field_box_surface(2,1,dim)), 1)
3908 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3909 st%electric_field_box_surface(1,2,dim)), 1)
3910 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3911 st%electric_field_box_surface(2,2,dim)), 1)
3912 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3913 st%electric_field_box_surface(1,3,dim)), 1)
3914 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3915 st%electric_field_box_surface(2,3,dim)), 1)
3916 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3917 st%electric_field_box_surface_plane_waves(1,1,dim)), 1)
3918 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3919 st%electric_field_box_surface_plane_waves(2,1,dim)), 1)
3920 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3921 st%electric_field_box_surface_plane_waves(1,2,dim)), 1)
3922 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3923 st%electric_field_box_surface_plane_waves(2,2,dim)), 1)
3924 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3925 st%electric_field_box_surface_plane_waves(1,3,dim)), 1)
3926 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3927 st%electric_field_box_surface_plane_waves(2,3,dim)), 1)
3928 call write_iter_nl(out_field_surf)
3936 type(c_ptr),
intent(inout) :: out_field_surf
3937 type(states_mxll_t),
intent(in) :: st
3938 integer,
intent(in) :: dim
3939 integer,
intent(in) :: iter
3943 integer :: n_columns
3945 if (.not. mpi_grp_is_root(mpi_world))
return
3952 call td_write_print_header_init(out_field_surf)
3955 call write_iter_header_start(out_field_surf)
3956 call write_iter_header(out_field_surf,
'- x direction')
3957 call write_iter_header(out_field_surf,
'+ x direction')
3958 call write_iter_header(out_field_surf,
'- y direction')
3959 call write_iter_header(out_field_surf,
'+ y direction')
3960 call write_iter_header(out_field_surf,
'- z direction')
3961 call write_iter_header(out_field_surf,
'+ z direction')
3962 call write_iter_header(out_field_surf,
'- x dir. p. w.')
3963 call write_iter_header(out_field_surf,
'+ x dir. p. w.')
3964 call write_iter_header(out_field_surf,
'- y dir. p. w.')
3965 call write_iter_header(out_field_surf,
'+ y dir. p. w.')
3966 call write_iter_header(out_field_surf,
'- z dir. p. w.')
3967 call write_iter_header(out_field_surf,
'+ z dir. p. w.')
3969 call write_iter_nl(out_field_surf)
3972 call write_iter_string(out_field_surf,
'#[Iter n.]')
3973 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%time)) //
']')
3975 do ii = 1, n_columns
3976 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(unit_one/units_out%length**2)) //
']')
3978 call write_iter_nl(out_field_surf)
3980 call td_write_print_header_end(out_field_surf)
3983 call write_iter_start(out_field_surf)
3984 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3985 st%magnetic_field_box_surface(1,1,dim)), 1)
3986 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3987 st%magnetic_field_box_surface(2,1,dim)), 1)
3988 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3989 st%magnetic_field_box_surface(1,2,dim)), 1)
3990 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3991 st%magnetic_field_box_surface(2,2,dim)), 1)
3992 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3993 st%magnetic_field_box_surface(1,3,dim)), 1)
3994 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3995 st%magnetic_field_box_surface(2,3,dim)), 1)
3996 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3997 st%magnetic_field_box_surface_plane_waves(1,1,dim)), 1)
3998 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3999 st%magnetic_field_box_surface_plane_waves(2,1,dim)), 1)
4000 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4001 st%magnetic_field_box_surface_plane_waves(1,2,dim)), 1)
4002 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4003 st%magnetic_field_box_surface_plane_waves(2,2,dim)), 1)
4004 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4005 st%magnetic_field_box_surface_plane_waves(1,3,dim)), 1)
4006 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4007 st%magnetic_field_box_surface_plane_waves(2,3,dim)), 1)
4008 call write_iter_nl(out_field_surf)
4014 subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
4015 type(c_ptr),
intent(inout) :: out_fields
4016 class(space_t),
intent(in) :: space
4017 type(states_mxll_t),
intent(in) :: st
4018 integer,
intent(in) :: iter
4019 real(real64),
intent(in) :: dt
4020 integer,
intent(in) :: e_or_b_field
4021 integer,
intent(in) :: field_type
4022 integer,
intent(in) :: idir
4026 real(real64) :: field(space%dim), selected_field
4027 character(len=80) :: aux
4029 if (.not. mpi_grp_is_root(mpi_world))
return
4034 call td_write_print_header_init(out_fields)
4037 write(aux,
'(a7,e20.12,3a)')
'# dt = ', units_from_atomic(units_out%time, dt), &
4038 " [", trim(units_abbrev(units_out%time)),
"]"
4039 call write_iter_string(out_fields, aux)
4040 call write_iter_nl(out_fields)
4042 call write_iter_header_start(out_fields)
4044 do id = 1, st%selected_points_number
4045 select case (e_or_b_field)
4047 write(aux,
'(a,i1,a)')
'E(', id,
')'
4049 write(aux,
'(a,i1,a)')
'B(', id,
')'
4051 call write_iter_header(out_fields, aux)
4054 call write_iter_nl(out_fields)
4055 call write_iter_string(out_fields,
'#[Iter n.]')
4056 call write_iter_header(out_fields,
' [' // trim(units_abbrev(units_out%time)) //
']')
4061 aux =
' [' // trim(units_abbrev(units_out%force)) //
']'
4062 do id = 1, st%selected_points_number
4063 call write_iter_header(out_fields, aux)
4065 call write_iter_nl(out_fields)
4066 call td_write_print_header_end(out_fields)
4069 call write_iter_start(out_fields)
4071 do id = 1, st%selected_points_number
4072 select case (e_or_b_field)
4075 select case (field_type)
4077 call get_electric_field_vector(st%selected_points_rs_state(:,id), field(1:st%dim))
4079 call get_electric_field_vector(st%selected_points_rs_state_long(:,id), field(1:st%dim))
4081 call get_electric_field_vector(st%selected_points_rs_state_trans(:,id), field(1:st%dim))
4083 selected_field = units_from_atomic(units_out%energy/units_out%length, field(idir))
4086 select case (field_type)
4088 call get_magnetic_field_vector(st%selected_points_rs_state(:,id), st%rs_sign, field(1:st%dim))
4090 call get_magnetic_field_vector(st%selected_points_rs_state_long(:,id), st%rs_sign, field(1:st%dim))
4092 call get_magnetic_field_vector(st%selected_points_rs_state_trans(:,id), st%rs_sign, field(1:st%dim))
4094 selected_field = units_from_atomic(unit_one/units_out%length**2, field(idir))
4096 call write_iter_double(out_fields, selected_field, 1)
4099 call write_iter_nl(out_fields)
4108 type(namespace_t),
intent(in) :: namespace
4109 class(space_t),
intent(in) :: space
4110 type(grid_t),
intent(inout) :: gr
4111 type(states_mxll_t),
intent(inout) :: st
4112 type(hamiltonian_mxll_t),
intent(inout) :: hm
4113 type(helmholtz_decomposition_t),
intent(inout) :: helmholtz
4114 type(output_t),
intent(in) :: outp
4115 integer,
intent(in) :: iter
4116 real(real64),
intent(in) :: time
4118 character(len=256) :: filename
4120 push_sub(td_write_maxwell_free_data)
4121 call profiling_in(
"TD_WRITE_MAXWELL_DATA")
4124 write(filename,
'(a,a,i7.7)') trim(outp%iter_dir),
"td.", iter
4126 call output_mxll(outp, namespace, space, gr, st, hm, helmholtz, time, filename)
4128 call profiling_out(
"TD_WRITE_MAXWELL_DATA")
4129 pop_sub(td_write_maxwell_free_data)
ssize_t ssize_t write(int __fd, const void *__buf, size_t __n) __attribute__((__access__(__read_only__
constant times a vector plus a vector
Copies a vector x, to a vector y.
Sets the iteration number to the C object.
Writes to the corresponding file and adds one to the iteration. Must be called after write_iter_init(...
double exp(double __x) __attribute__((__nothrow__
This module contains interfaces for BLAS routines You should not use these routines directly....
This module implements a calculator for the density and defines related functions.
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
integer, parameter, public spinors
integer, parameter, public spin_polarized
subroutine, public excited_states_kill(excited_state)
Kills an excited_state structure.
subroutine, public excited_states_init(excited_state, ground_state, filename, namespace)
Fills in an excited_state structure, by reading a file called "filename". This file describes the "pr...
type(gauge_field_t) function, pointer, public list_get_gauge_field(partners)
logical function, public list_has_lasers(partners)
type(lasers_t) function, pointer, public list_get_lasers(partners)
subroutine, public gauge_field_output_write(this, out_gauge, iter)
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
complex(real64), parameter, public m_z0
real(real64), parameter, public lmm_r_single_atom
Default local magnetic moments sphere radius for an isolated system.
complex(real64), parameter, public m_zi
integer, parameter, public max_output_types
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
This module implements the underlying real-space grid.
The Helmholtz decomposition is intended to contain "only mathematical" functions and procedures to co...
This module defines classes and functions for interaction partners.
subroutine, public zio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
subroutine, public io_function_read_what_how_when(namespace, space, what, how, output_interval, what_tag_in, how_tag_in, output_interval_tag_in, ignore_error)
subroutine, public io_close(iunit, grp)
character(len=max_path_len) function, public io_workpath(path, namespace)
subroutine, public io_rm(fname, namespace)
subroutine, public io_mkdir(fname, namespace, parents)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
real(real64) function, public ion_dynamics_temperature(ions)
This function returns the ionic temperature in energy units.
integer, parameter, public qkickmode_cos
integer, parameter, public qkickmode_none
integer, parameter, public qkickmode_sin
subroutine, public kick_function_get(space, mesh, kick, kick_function, iq, to_interpolate)
subroutine, public kick_write(kick, iunit, out)
integer, parameter, public qkickmode_bessel
subroutine, public lasers_nondipole_laser_field_step(this, field, time)
Retrieves the NDSFA vector_potential correction. The nondipole field is obtained for consecutive time...
subroutine, public lasers_set_nondipole_parameters(this, ndfield, nd_integration_time)
Set parameters for nondipole SFA calculation.
logical function, public lasers_with_nondipole_field(lasers)
Check if a nondipole SFA correction should be computed for the given laser.
integer, parameter, public e_field_electric
integer, parameter, public e_field_vector_potential
integer, parameter, public e_field_scalar_potential
integer pure elemental function, public laser_kind(laser)
subroutine, public laser_field(laser, field, time)
Retrieves the value of either the electric or the magnetic field. If the laser is given by a scalar p...
integer, parameter, public e_field_magnetic
integer, parameter, public dft_u_none
integer, parameter, public dft_u_acbn0
subroutine, public magnetic_local_moments(mesh, st, ions, boundaries, rho, rr, lmm)
subroutine, public magnetic_moment(mesh, st, rho, mm)
subroutine, public magnetic_total_magnetization(mesh, st, qq, trans_mag)
This module is intended to contain "only mathematical" functions and procedures.
subroutine, public ylmr_real(xx, li, mi, ylm)
This is a Numerical Recipes-based subroutine computes real spherical harmonics ylm at position (x,...
This module defines functions over batches of mesh functions.
This module defines various routines, operating on mesh functions.
complex(real64) function, public zmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
subroutine, public dmf_multipoles(mesh, ff, lmax, multipole, mask)
This routine calculates the multipoles of a function ff.
real(real64) function, public dmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_not_implemented(feature, namespace)
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 modelmb_sym_all_states(space, mesh, st)
This module contains some common usage patterns of MPI routines.
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
type(mpi_grp_t), public mpi_world
This module handles the communicators for the various parallelization strategies.
integer, parameter, public velocity_gauge_dipole
this module contains the low-level part of the output system
subroutine, public output_modelmb(outp, namespace, space, dir, gr, ions, iter, st)
this module contains the output system
subroutine, public output_all(outp, namespace, space, dir, gr, ions, iter, st, hm, ks)
subroutine, public output_scalar_pot(outp, namespace, space, dir, mesh, ions, ext_partners, time)
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
integer, parameter, public restart_gs
subroutine, public restart_init(restart, namespace, data_type, type, mc, ierr, mesh, dir, exact)
Initializes a restart object.
integer, parameter, public restart_proj
integer, parameter, public restart_type_load
subroutine, public restart_end(restart)
subroutine, public zstates_elec_matrix(st1, st2, mesh, aa)
subroutine, public zstates_elec_calc_projections(st, gs_st, namespace, mesh, ik, proj, gs_nst)
This routine computes the projection between two set of states.
This module handles spin dimensions of the states and the k-point distribution.
logical function, public state_kpt_is_local(st, ist, ik)
check whether a given state (ist, ik) is on the local node
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine, public states_elec_deallocate_wfns(st)
Deallocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
subroutine, public states_elec_look(restart, nik, dim, nst, ierr)
Reads the 'states' file in the restart directory, and finds out the nik, dim, and nst contained in it...
This module handles reading and writing restart information for the states_elec_t.
subroutine, public states_elec_load(restart, namespace, space, st, mesh, kpoints, ierr, iter, lr, lowest_missing, label, verbose, skip)
returns in ierr: <0 => Fatal error, or nothing read =0 => read all wavefunctions >0 => could only rea...
subroutine, public td_calc_tacc(namespace, space, gr, ions, ext_partners, st, hm, acc, time)
Electronic acceleration (to calculate harmonic spectrum...) It is calculated as:
subroutine, public td_calc_ionch(mesh, st, ch, Nch)
Multiple ionization probabilities calculated form the KS orbital densities C. Ullrich,...
subroutine, public td_calc_tvel(gr, st, space, kpoints, vel)
Electronic velocity (to calculate harmonic spectrum...) It is calculated as:
subroutine, public td_write_coordinates(out_coords, natoms, space, pos, vel, tot_forces, iter)
subroutine, public td_write_sep_coordinates(out_coords, natoms, space, pos, vel, tot_forces, iter, which)
subroutine, public td_write_print_header_init(out)
subroutine, public td_write_print_header_end(out)
subroutine td_write_magnetic_field_box_surface(out_field_surf, st, dim, iter)
integer, parameter, public out_total_current
integer, parameter maxwell_b_field
integer, parameter, public out_maxwell_max
subroutine td_write_proj(out_proj, space, mesh, ions, st, gs_st, kick, iter)
integer, parameter, public out_q
integer, parameter, public out_mxll_field
subroutine calc_projections(mesh, st, gs_st, projections)
This subroutine calculates:
integer, parameter out_b_field_surface_y
subroutine td_write_ionch(out_ionch, mesh, st, iter)
integer, parameter, public out_tot_m
integer, parameter, public out_norm_ks
integer, parameter out_maxwell_trans_b_field
subroutine td_write_multipole_r(out_multip, space, mesh, ions, st, lmax, kick, rho, iter, mpi_grp)
Write multipoles to the corresponding file.
integer, parameter, public out_proj
integer, parameter, public out_partial_charges
integer, parameter, public out_separate_coords
subroutine td_write_total_current(out_total_current, space, mesh, st, iter)
subroutine, public td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, dt)
subroutine td_write_effective_u(out_coords, lda_u, iter)
integer, parameter maxwell_trans_field
subroutine td_write_acc(out_acc, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
subroutine, public td_write_mxll_init(writ, namespace, iter, dt)
subroutine, public td_write_mxll_end(writ)
subroutine td_write_mxll_field(out_mxll, space, hm, dt, iter)
integer, parameter out_b_field_surface_x
integer, parameter out_maxwell_long_e_field
integer, parameter, public out_kp_proj
integer, parameter, public out_magnets
subroutine td_write_multipole(out_multip, space, gr, ions, st, lmax, kick, iter)
Top-level routine that write multipoles to file, or files depending on whether a state-resolved outpu...
subroutine td_write_electric_field_box_surface(out_field_surf, st, dim, iter)
integer, parameter out_e_field_surface_y
integer, parameter, public out_angular
subroutine td_write_populations(out_populations, namespace, space, mesh, st, writ, dt, iter)
integer, parameter, public out_max
integer, parameter out_maxwell_long_b_field
integer, parameter, public out_energy
integer, parameter, public out_spin
subroutine td_write_ftchd(out_ftchd, space, mesh, st, kick, iter)
subroutine td_write_partial_charges(out_partial_charges, mesh, st, ions, iter)
integer, parameter out_dftu_max
For the Maxwell fields we increment in steps of 3 to leave room for x, y, and z output.
subroutine td_write_laser(out_laser, space, lasers, dt, iter)
integer, parameter out_maxwell_total_b_field
integer, parameter, public out_ftchd
integer, parameter, public out_separate_velocity
subroutine td_write_tot_mag(out_magnets, mesh, st, kick, iter)
subroutine td_write_q(out_q, space, ks, iter)
integer, parameter, public out_floquet
integer, parameter, public out_acc
integer, parameter, public out_ion_ch
integer, parameter maxwell_long_field
integer, parameter, public out_n_ex
integer, parameter out_b_field_surface_z
subroutine td_write_proj_kp(mesh, kpoints, st, gs_st, namespace, iter)
integer, parameter, public out_temperature
subroutine td_write_norm_ks_orbitals(file_handle, grid, kpoints, st, iter)
Write the norm of the KS orbitals to file as a function of time step.
subroutine, public td_write_data(writ)
subroutine, public td_write_mxll_free_data(writ, namespace, space, gr, st, hm, helmholtz, outp, iter, time)
subroutine td_write_total_heat_current(write_obj, space, hm, gr, st, iter)
integer, parameter out_e_field_surface_z
integer, parameter maxwell_total_field
integer, parameter, public out_coords
integer, parameter out_maxwell_total_e_field
integer, parameter, public out_laser
integer, parameter, public out_eigs
subroutine td_write_vel(out_vel, gr, st, space, kpoints, iter)
subroutine td_write_floquet(namespace, space, hm, ext_partners, mesh, st, iter)
integer, parameter, public out_total_heat_current
integer, parameter out_e_field_surface_x
subroutine, public td_write_kick(outp, namespace, space, mesh, kick, ions, iter)
subroutine, public td_write_end(writ)
subroutine td_write_angular(out_angular, namespace, space, gr, ions, hm, st, kick, iter)
Computes and outputs the orbital angular momentum defined by.
subroutine td_write_eigs(out_eigs, st, iter)
subroutine, public td_write_iter(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, kick, ks, dt, iter, mc, recalculate_gs)
subroutine td_write_n_ex(out_nex, outp, namespace, mesh, kpoints, st, gs_st, iter)
This routine computes the total number of excited electrons based on projections on the GS orbitals T...
subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
subroutine td_write_spin(out_spin, mesh, st, iter)
integer, parameter, public out_vel
integer, parameter, public out_gauge_field
subroutine td_write_temperature(out_temperature, ions, iter)
integer, parameter maxwell_e_field
integer, parameter, public out_populations
subroutine, public td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
subroutine td_write_local_magnetic_moments(out_magnets, gr, st, ions, lmm_r, iter)
subroutine, public td_write_init(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, ks, ions_move, with_gauge_field, kick, iter, max_iter, dt, mc)
Initialize files to write when prograting in time.
integer, parameter out_maxwell_energy
subroutine td_write_energy(out_energy, hm, iter, ke)
subroutine td_write_maxwell_energy(out_maxwell_energy, hm, iter)
integer, parameter, public out_separate_forces
integer, parameter out_maxwell_trans_e_field
type(type_t), public type_cmplx
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
type(unit_t), public unit_kelvin
For converting energies into temperatures.
type(unit_system_t), public units_inp
the units systems for reading and writing
type(unit_t), public unit_one
some special units required for particular quantities
This module is intended to contain simple general-purpose utility functions and procedures.
character pure function, public index2axis(idir)
subroutine, public v_ks_calculate_current(this, calc_cur)
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Explicit interfaces to C functions, defined in write_iter_low.cc.
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Describes mesh distribution to nodes.
This is defined even when running serial.
Stores all communicators and groups.
The states_elec_t class contains all electronic wave functions.
Time-dependent Write Properties.
subroutine dipole_matrix_elements(dir)