48 use,
intrinsic :: iso_fortran_env
108 integer,
parameter,
public :: &
109 OUT_MULTIPOLES = 1, &
146 character(len=100) :: td_file_name(OUT_MAX) = [character(100) :: &
156 "magnetic_moments", &
171 "total_heat_current", &
172 "total_magnetization", &
174 "maxwell_dipole_field", &
175 "norm_wavefunctions", &
178 integer,
parameter :: &
179 OUT_DFTU_EFFECTIVE_U = 1, &
183 integer,
parameter :: &
184 OUT_MAXWELL_TOTAL_E_FIELD = 1, &
198 integer,
parameter,
public :: &
201 integer,
parameter :: &
202 MAXWELL_TOTAL_FIELD = 1, &
206 integer,
parameter :: &
212 type(c_ptr) :: handle
213 type(c_ptr),
allocatable :: mult_handles(:)
215 integer :: hand_start
217 logical ::
write = .false.
218 logical :: resolve_states = .false.
230 real(real64) :: lmm_r
233 integer :: n_excited_states
235 integer :: compute_interval
245 class(
mesh_t),
intent(in) :: mesh
246 type(
kick_t),
intent(in) :: kick
247 type(
ions_t),
intent(in) :: ions
248 integer,
intent(in) :: iter
250 complex(real64),
allocatable :: kick_function(:)
251 character(len=256) :: filename
256 write(filename,
'(a,i7.7)')
"td.", iter
257 if (outp%what(option__output__delta_perturbation))
then
258 safe_allocate(kick_function(1:mesh%np))
260 call zio_function_output(outp%how(option__output__delta_perturbation), filename,
"kick_function", namespace, &
261 space, mesh, kick_function(:),
units_out%energy, err, pos=ions%pos, atoms=ions%atom)
262 safe_deallocate_a(kick_function)
278 subroutine td_write_init(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, ks, ions_move, &
279 with_gauge_field, kick, iter, max_iter, dt, mc)
283 type(
output_t),
intent(inout) :: outp
284 type(
grid_t),
intent(in) :: gr
287 type(
ions_t),
intent(in) :: ions
289 type(
v_ks_t),
intent(inout) :: ks
290 logical,
intent(in) :: ions_move
291 logical,
intent(in) :: with_gauge_field
292 type(
kick_t),
intent(in) :: kick
293 integer,
intent(in) :: iter
294 integer,
intent(in) :: max_iter
295 real(real64),
intent(in) :: dt
299 integer :: ierr, first, ii, ist, jj, flags, iout, default, ifile
304 character(len=MAX_PATH_LEN) :: filename
306 logical :: resolve_states
307 logical,
allocatable :: skip(:)
450 output_options = .false.
451 output_options(out_multipoles) = .
true.
466 writ%out(iout)%write = output_options(iout)
481 if (space%is_periodic() .and. writ%out(
out_angular)%write)
then
496 if (gr%np /= gr%np_global)
then
497 message(1) =
"TDOutput option td_kpoint_occup and td_floquet do not work with domain parallelization"
504 'Forces for systems periodic in 1D are not currently implemented and options that output the forces are not allowed.')
508 if (writ%out(
out_kp_proj)%write .and. hm%kpoints%nik_skip == 0)
then
509 message(1) =
"TDOutput option td_kpoint_occup only work with zero-weight k-points at the moment."
523 call parse_variable(namespace,
'TDOutputResolveStates', .false., resolve_states)
524 if (.not. writ%out(out_multipoles)%write)
then
525 write(
message(1),
'(a)')
"TDOutputResolveStates works only for TDOutput = multipoles."
538 if (writ%lmax < 0)
then
539 write(
message(1),
'(a,i6,a)')
"Input: '", writ%lmax,
"' is not a valid TDMultipoleLmax."
540 message(2) =
'(Must be TDMultipoleLmax >= 0 )'
546 if ((writ%out(
out_acc)%write) .and. ions_move)
then
547 message(1) =
'If harmonic spectrum is to be calculated, atoms should not be allowed to move.'
551 if ((writ%out(
out_q)%write) .and. .not.(ks%has_photons))
then
552 message(1) =
'If q(t) is to be calculated, you need to allow for photon modes.'
557 .or. hm%mxll%add_electric_dip))
then
558 message(1) =
'If the dipolar field has to be written, MaxwellCouplingMode has to be'
559 message(2) =
'"lenght_gauge_dipole" or "velocity_gauge_dipole" and at least one Maxwell system'
560 message(3) =
'must be present.'
564 rmin = ions%min_distance()
574 message(1) =
"Option TDOutput = populations is not implemented for parallel in states."
586 safe_deallocate_a(writ%gs_st%node)
594 writ%gs_st%st_end = writ%gs_st%nst
596 message(1) =
"Unable to read states information."
600 writ%gs_st%st_start = 1
616 call parse_variable(namespace,
'TDProjStateStart', 1, writ%gs_st%st_start)
618 if (st%parallel_in_states .and. writ%out(
out_proj)%write .and. writ%gs_st%st_start > 1)
then
624 writ%gs_st%lnst = writ%gs_st%st_end - writ%gs_st%st_start + 1
628 writ%gs_st%parallel_in_states = .false.
631 safe_allocate(writ%gs_st%occ(1:writ%gs_st%nst, 1:writ%gs_st%nik))
632 safe_allocate(writ%gs_st%eigenval(1:writ%gs_st%nst, 1:writ%gs_st%nik))
636 safe_allocate(writ%gs_st%node(1:writ%gs_st%nst))
637 writ%gs_st%node(:) = 0
639 writ%gs_st%eigenval = huge(writ%gs_st%eigenval)
641 if (writ%gs_st%d%ispin ==
spinors)
then
642 safe_deallocate_a(writ%gs_st%spin)
643 safe_allocate(writ%gs_st%spin(1:3, 1:writ%gs_st%nst, 1:writ%gs_st%nik))
646 safe_allocate(skip(1:writ%gs_st%nst))
648 skip(1:writ%gs_st%st_start-1) = .
true.
652 safe_deallocate_a(skip)
657 if (ierr /= 0 .and. ierr /= (writ%gs_st%st_end-writ%gs_st%st_start+1)*writ%gs_st%nik &
658 *writ%gs_st%d%dim*writ%gs_st%mpi_grp%size)
then
659 message(1) =
"Unable to read wavefunctions for TDOutput."
704 if (
parse_block(namespace,
'TDExcitedStatesToProject', blk) == 0)
then
706 safe_allocate(writ%excited_st(1:writ%n_excited_states))
707 do ist = 1, writ%n_excited_states
712 writ%n_excited_states = 0
726 call parse_variable(namespace,
'TDOutputComputeInterval', 50, writ%compute_interval)
727 if (writ%compute_interval < 0)
then
728 message(1) =
"TDOutputComputeInterval must be >= 0."
744 call io_mkdir(
'td.general', namespace)
755 do ifile = 1, out_max
759 if (writ%out(ifile)%write)
then
763 trim(
io_workpath(
"td.general/"//trim(td_file_name(ifile)), namespace)))
771 if (writ%out(out_multipoles)%write .and. .not. resolve_states)
then
774 trim(
io_workpath(
"td.general/multipoles", namespace)))
778 select case (kick%qkick_mode)
780 write(filename,
'(a)')
'td.general/ftchd.sin'
782 write(filename,
'(a)')
'td.general/ftchd.cos'
784 write(filename,
'(a, SP, I0.3, a, I0.3)')
'td.general/ftchd.l', kick%qbessel_l,
'_m', kick%qbessel_m
786 write(filename,
'(a)')
'td.general/ftchd'
798 call io_rm(
"td.general/laser", namespace=namespace)
815 if (writ%out(out_multipoles)%write .and. resolve_states)
then
817 writ%out(out_multipoles)%hand_start = st%st_start
818 writ%out(out_multipoles)%hand_end = st%st_end
819 writ%out(out_multipoles)%resolve_states = .
true.
820 writ%out(out_multipoles)%mpi_grp = gr%mpi_grp
822 safe_allocate(writ%out(out_multipoles)%mult_handles(st%st_start:st%st_end))
824 if (writ%out(out_multipoles)%mpi_grp%is_root())
then
825 do ist = st%st_start, st%st_end
826 write(filename,
'(a,i4.4)')
'td.general/multipoles-ist', ist
839 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
840 calc_eigenval=.false., time = iter*dt, calc_energy = .false.)
843 if (writ%out(
out_n_ex)%write .and. writ%compute_interval > 0)
then
844 call io_mkdir(outp%iter_dir, namespace)
847 if (all(outp%how == 0) .and. writ%out(
out_n_ex)%write)
then
869 if (hm%lda_u_level ==
dft_u_acbn0) default = default + 2**(out_dftu_effective_u - 1)
877 writ%out_dftu(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
881 if (writ%out_dftu(out_dftu_effective_u)%write)
then
884 trim(
io_workpath(
"td.general/effectiveU", namespace)))
902 if (writ%out(iout)%write)
then
903 if (writ%out(iout)%mpi_grp%is_root())
then
904 if (writ%out(iout)%resolve_states)
then
905 do ist = writ%out(iout)%hand_start, writ%out(iout)%hand_end
908 safe_deallocate_a(writ%out(iout)%mult_handles)
918 if (writ%out_dftu(iout)%write)
call write_iter_end(writ%out_dftu(iout)%handle)
924 do ist = 1, writ%n_excited_states
927 writ%n_excited_states = 0
940 subroutine td_write_iter(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, kick, ks, dt, iter, mc, recalculate_gs)
943 class(
space_t),
intent(in) :: space
945 type(
grid_t),
intent(in) :: gr
948 type(
ions_t),
intent(inout) :: ions
950 type(
kick_t),
intent(in) :: kick
951 type(
v_ks_t),
intent(in) :: ks
952 real(real64),
intent(in) :: dt
953 integer,
intent(in) :: iter
955 logical,
intent(in) :: recalculate_gs
963 if (writ%out(out_multipoles)%write)
then
964 call td_write_multipole(writ%out(out_multipoles), space, gr, ions, st, writ%lmax, kick, iter)
987 if (writ%out(
out_proj)%write .and. mod(iter, writ%compute_interval) == 0)
then
996 if (writ%out(
out_kp_proj)%write .and. mod(iter, writ%compute_interval) == 0)
then
1002 ions%pos, ions%vel, ions%tot_force, iter)
1007 ions%pos, ions%vel, ions%tot_force, iter, 1)
1012 ions%pos, ions%vel, ions%tot_force, iter, 2)
1017 ions%pos, ions%vel, ions%tot_force, iter, 3)
1028 if (writ%out(
out_acc)%write)
then
1029 call td_write_acc(writ%out(
out_acc)%handle, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
1032 if (writ%out(
out_vel)%write)
then
1045 if(
associated(gfield))
then
1071 if (writ%out(
out_n_ex)%write .and. mod(iter, writ%compute_interval) == 0)
then
1073 if (recalculate_gs)
then
1076 ierr, label =
': Houston states for TDOutput')
1080 call td_write_n_ex(writ%out(
out_n_ex)%handle, outp, namespace, gr, hm%kpoints, st, writ%gs_st, iter)
1092 if (writ%out_dftu(out_dftu_effective_u)%write)
then
1096 if (writ%out(
out_q)%write .and. ks%has_photons)
then
1100 if (writ%out(
out_mxll_field)%write .and. hm%mxll%calc_field_dip)
then
1117 do iout = 1, out_max
1119 if (writ%out(iout)%write)
then
1120 if (writ%out(iout)%mpi_grp%is_root())
then
1121 if (writ%out(iout)%resolve_states)
then
1122 do ii = writ%out(iout)%hand_start, writ%out(iout)%hand_end
1134 if (writ%out_dftu(iout)%write)
call write_iter_flush(writ%out_dftu(iout)%handle)
1143 subroutine td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, dt)
1146 type(
grid_t),
intent(in) :: gr
1149 type(
v_ks_t),
intent(inout) :: ks
1151 type(
ions_t),
intent(in) :: ions
1153 integer,
intent(in) :: iter
1154 real(real64),
optional,
intent(in) :: dt
1156 character(len=256) :: filename
1162 if (st%modelmbparticles%nparticle > 0)
then
1167 write(filename,
'(a,a,i7.7)') trim(outp%iter_dir),
"td.", iter
1169 call output_all(outp, namespace, space, filename, gr, ions, iter, st, hm, ks)
1171 call output_modelmb(outp, namespace, space, filename, gr, ions, iter, st)
1173 if (
present(dt))
then
1174 call output_scalar_pot(outp, namespace, space, filename, gr, ions, ext_partners, iter*dt)
1176 if (iter == 0)
call output_scalar_pot(outp, namespace, space, filename, gr, ions, ext_partners)
1185 type(c_ptr),
intent(inout) ::
out_spin
1186 class(
mesh_t),
intent(in) :: mesh
1188 integer,
intent(in) :: iter
1190 character(len=130) :: aux
1191 real(real64) :: spin(3)
1207 if (st%d%ispin ==
spinors)
then
1208 write(aux,
'(a2,18x)')
'Sx'
1210 write(aux,
'(a2,18x)')
'Sy'
1213 write(aux,
'(a2,18x)')
'Sz'
1221 select case (st%d%ispin)
1240 type(
ions_t),
intent(in) :: ions
1241 real(real64),
intent(in) :: lmm_r
1242 integer,
intent(in) :: iter
1245 character(len=50) :: aux
1246 real(real64),
allocatable :: lmm(:,:)
1251 safe_allocate(lmm(1:3, 1:ions%natoms))
1261 do ia = 1, ions%natoms
1262 if (st%d%ispin ==
spinors)
then
1263 write(aux,
'(a2,i2.2,16x)')
'mx', ia
1265 write(aux,
'(a2,i2.2,16x)')
'my', ia
1268 write(aux,
'(a2,i2.2,16x)')
'mz', ia
1277 do ia = 1, ions%natoms
1278 select case (st%d%ispin)
1286 safe_deallocate_a(lmm)
1294 type(c_ptr),
intent(inout) :: out_magnets
1295 class(
mesh_t),
intent(in) :: mesh
1297 type(
kick_t),
intent(in) :: kick
1298 integer,
intent(in) :: iter
1300 complex(real64),
allocatable :: tm(:,:)
1305 safe_allocate(tm(1:6,1:kick%nqvec))
1307 do iq = 1, kick%nqvec
1337 do iq = 1, kick%nqvec
1346 safe_deallocate_a(tm)
1360 subroutine td_write_angular(out_angular, namespace, space, gr, ions, hm, st, kick, iter)
1361 type(c_ptr),
intent(inout) :: out_angular
1363 class(
space_t),
intent(in) :: space
1364 type(
grid_t),
intent(in) :: gr
1365 type(
ions_t),
intent(inout) :: ions
1368 type(
kick_t),
intent(in) :: kick
1369 integer,
intent(in) :: iter
1372 character(len=130) :: aux
1373 real(real64) :: angular(3)
1380 call angular_momentum%setup_dir(idir)
1383 m_two*real(angular_momentum%zstates_elec_expectation_value(namespace, space, gr, hm, st), real64)
1385 safe_deallocate_p(angular_momentum)
1392 write(aux,
'(a15,i2)')
'# nspin ', st%d%nspin
1400 write(aux,
'(a4,18x)')
'<Lx>'
1402 write(aux,
'(a4,18x)')
'<Ly>'
1404 write(aux,
'(a4,18x)')
'<Lz>'
1433 class(
space_t),
intent(in) :: space
1434 type(
grid_t),
intent(in) :: gr
1435 type(
ions_t),
intent(in) :: ions
1437 integer,
intent(in) :: lmax
1438 type(
kick_t),
intent(in) :: kick
1439 integer,
intent(in) :: iter
1442 real(real64),
allocatable :: rho(:,:)
1446 if (out_multip%resolve_states)
then
1447 safe_allocate(rho(1:gr%np_part, 1:st%d%nspin))
1450 do ist = st%st_start, st%st_end
1452 call td_write_multipole_r(out_multip%mult_handles(ist), space, gr, ions, st, lmax, kick, rho, iter, &
1453 mpi_grp = out_multip%mpi_grp)
1456 safe_deallocate_a(rho)
1459 if (
allocated(st%frozen_rho))
then
1460 safe_allocate(rho(1:gr%np, 1:st%d%nspin))
1461 call lalg_copy(gr%np, st%d%nspin, st%rho, rho)
1466 safe_deallocate_a(rho)
1478 subroutine td_write_multipole_r(out_multip, space, mesh, ions, st, lmax, kick, rho, iter, mpi_grp)
1479 type(c_ptr),
intent(inout) :: out_multip
1480 class(
space_t),
intent(in) :: space
1481 class(
mesh_t),
intent(in) :: mesh
1482 type(
ions_t),
intent(in) :: ions
1484 integer,
intent(in) :: lmax
1485 type(
kick_t),
intent(in) :: kick
1486 real(real64),
intent(in) :: rho(:,:)
1487 integer,
intent(in) :: iter
1488 type(
mpi_grp_t),
optional,
intent(in) :: mpi_grp
1491 integer :: is, idir, ll, mm, add_lm
1492 character(len=120) :: aux
1493 real(real64) :: ionic_dipole(ions%space%dim)
1494 real(real64),
allocatable :: multipole(:,:)
1500 assert(.not. (lmax > 1 .and. space%dim > 3))
1503 if (
present(mpi_grp)) mpi_grp_ = mpi_grp
1505 if (mpi_grp_%is_root().and.iter == 0)
then
1508 write(aux,
'(a15,i2)')
'# nspin ', st%d%nspin
1512 write(aux,
'(a15,i2)')
'# lmax ', lmax
1520 do is = 1, st%d%nspin
1521 write(aux,
'(a18,i1,a1)')
'Electronic charge(', is,
')'
1524 do idir = 1, space%dim
1525 write(aux,
'(4a1,i1,a1)')
'<',
index2axis(idir),
'>',
'(', is,
')'
1531 write(aux,
'(a2,i2,a4,i2,a2,i1,a1)')
'l=', ll,
', m=', mm,
' (', is,
')'
1542 do is = 1, st%d%nspin
1545 do idir = 1, space%dim
1561 if (space%dim > 3 .and. lmax == 1)
then
1563 safe_allocate(multipole(1:space%dim+1, 1:st%d%nspin))
1565 safe_allocate(multipole(1:(lmax + 1)**2, 1:st%d%nspin))
1569 do is = 1, st%d%nspin
1574 ionic_dipole = ions%dipole()
1575 do is = 1, st%d%nspin
1576 multipole(2:space%dim+1, is) = -ionic_dipole(1:space%dim)/st%d%nspin - multipole(2:space%dim+1, is)
1580 if (mpi_grp_%is_root())
then
1582 do is = 1, st%d%nspin
1585 do idir = 1, space%dim
1589 add_lm = space%dim + 2
1600 safe_deallocate_a(multipole)
1605 subroutine td_write_ftchd(out_ftchd, space, mesh, st, kick, iter)
1606 type(c_ptr),
intent(inout) :: out_ftchd
1607 class(
space_t),
intent(in) :: space
1608 class(
mesh_t),
intent(in) :: mesh
1610 type(
kick_t),
intent(in) :: kick
1611 integer,
intent(in) :: iter
1613 integer :: is, ip, idir
1614 character(len=120) :: aux, aux2
1615 real(real64) :: ftchd_bessel
1616 complex(real64) :: ftchd
1618 real(real64),
allocatable :: integrand_bessel(:)
1619 complex(real64),
allocatable :: integrand(:)
1623 if (
mpi_world%is_root().and.iter == 0)
then
1626 write(aux,
'(a15, i2)')
'# qkickmode ', kick%qkick_mode
1631 write(aux,
'(a15, i0.3, 1x, i0.3)')
'# ll, mm ', kick%qbessel_l, kick%qbessel_m
1637 write(aux,
'(a15, f9.6)')
'# qlength ', kick%qlength
1639 write(aux,
'(a15)')
'# qvector '
1640 do idir = 1, space%dim
1641 write(aux2,
'(f9.5)') kick%qvector(idir,1)
1642 aux = trim(aux) // trim(aux2)
1648 write(aux,
'(a15,f18.12)')
'# kick strength', kick%delta_strength
1654 write(aux,
'(a17)')
'int(j_l*Y_lm*rho)'
1656 write(aux,
'(a12)')
'Real, Imag'
1673 safe_allocate(integrand(1:mesh%np))
1675 do is = 1, st%d%nspin
1677 integrand(ip) = integrand(ip) + st%rho(ip, is) *
exp(-
m_zi*sum(mesh%x(ip, 1:space%dim)*kick%qvector(1:space%dim, 1)))
1681 safe_deallocate_a(integrand)
1684 safe_allocate(integrand_bessel(1:mesh%np))
1685 integrand_bessel =
m_zero
1686 do is = 1, st%d%nspin
1688 call ylmr_real(mesh%x(ip, 1:3), kick%qbessel_l, kick%qbessel_m, ylm)
1689 integrand_bessel(ip) = integrand_bessel(ip) + st%rho(ip, is) * &
1690 loct_sph_bessel(kick%qbessel_l, kick%qlength*norm2(mesh%x(ip, :)))*ylm
1694 safe_deallocate_a(integrand_bessel)
1713 type(c_ptr),
intent(inout) :: out_temperature
1714 type(
ions_t),
intent(in) :: ions
1715 integer,
intent(in) :: iter
1750 type(c_ptr),
intent(inout) :: out_populations
1752 class(
space_t),
intent(in) :: space
1753 class(
mesh_t),
intent(in) :: mesh
1756 real(real64),
intent(in) :: dt
1757 integer,
intent(in) :: iter
1760 character(len=6) :: excited_name
1761 complex(real64) :: gsp
1762 complex(real64),
allocatable :: excited_state_p(:)
1763 complex(real64),
allocatable :: dotprodmatrix(:, :, :)
1768 safe_allocate(dotprodmatrix(1:writ%gs_st%nst, 1:st%nst, 1:st%nik))
1773 assert(.not. space%is_periodic())
1778 if (writ%n_excited_states > 0)
then
1779 safe_allocate(excited_state_p(1:writ%n_excited_states))
1780 do ist = 1, writ%n_excited_states
1781 excited_state_p(ist) =
zstates_elec_mpdotp(namespace, mesh, writ%excited_st(ist), st, dotprodmatrix)
1793 do ist = 1, writ%n_excited_states
1794 write(excited_name,
'(a2,i3,a1)')
'P(', ist,
')'
1813 do ist = 1, writ%n_excited_states
1820 if (writ%n_excited_states > 0)
then
1821 safe_deallocate_a(excited_state_p)
1823 safe_deallocate_a(dotprodmatrix)
1829 subroutine td_write_acc(out_acc, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
1830 type(c_ptr),
intent(inout) :: out_acc
1832 class(
space_t),
intent(in) :: space
1833 type(
grid_t),
intent(in) :: gr
1834 type(
ions_t),
intent(inout) :: ions
1838 real(real64),
intent(in) :: dt
1839 integer,
intent(in) :: iter
1842 character(len=7) :: aux
1843 real(real64) :: acc(space%dim)
1847 if (iter == 0 .and.
mpi_world%is_root())
then
1852 do idim = 1, space%dim
1853 write(aux,
'(a4,i1,a1)')
'Acc(', idim,
')'
1861 do idim = 1, space%dim
1868 call td_calc_tacc(namespace, space, gr, ions, ext_partners, st, hm, acc, dt*iter)
1881 subroutine td_write_vel(out_vel, gr, st, space, kpoints, iter)
1882 type(c_ptr),
intent(inout) :: out_vel
1883 type(
grid_t),
intent(in) :: gr
1885 type(
space_t),
intent(in) :: space
1887 integer,
intent(in) :: iter
1890 character(len=7) :: aux
1891 real(real64) :: vel(space%dim)
1895 if (iter == 0 .and.
mpi_world%is_root())
then
1900 do idim = 1, space%dim
1901 write(aux,
'(a4,i1,a1)')
'Vel(', idim,
')'
1909 do idim = 1, space%dim
1931 type(c_ptr),
intent(inout) :: out_laser
1932 class(
space_t),
intent(in) :: space
1933 type(
lasers_t),
intent(inout) :: lasers
1934 real(real64),
intent(in) :: dt
1935 integer,
intent(in) :: iter
1938 real(real64) :: field(space%dim)
1939 real(real64) :: ndfield(space%dim)
1940 character(len=80) :: aux
1956 do il = 1, lasers%no_lasers
1959 do idir = 1, space%dim
1960 write(aux,
'(a,i1,a)')
'E(', idir,
')'
1964 do idir = 1, space%dim
1965 write(aux,
'(a,i1,a)')
'B(', idir,
')'
1969 do idir = 1, space%dim
1970 write(aux,
'(a,i1,a)')
'A(', idir,
')'
1974 write(aux,
'(a,i1,a)')
'e(t)'
1980 do idir = 1, space%dim
1981 write(aux,
'(a,i1,a)')
'A^M(', idir,
')'
1993 do il = 1, lasers%no_lasers
1997 do idir = 1, space%dim
2002 do idir = 1, space%dim
2013 do idir = 1, space%dim
2026 do il = 1, lasers%no_lasers
2028 call laser_field(lasers%lasers(il), field(1:space%dim), iter*dt)
2053 type(c_ptr),
intent(inout) :: out_energy
2055 integer,
intent(in) :: iter
2056 real(real64),
intent(in) :: ke
2060 integer :: n_columns
2083 if (hm%pcm%run_pcm)
then
2085 n_columns = n_columns + 1
2090 n_columns = n_columns + 1
2100 do ii = 1, n_columns
2123 hm%energy%int_nn_pcm + hm%energy%int_ne_pcm), 1)
2136 type(c_ptr),
intent(inout) :: out_eigs
2138 integer,
intent(in) :: iter
2141 character(len=68) :: buf
2154 write(buf,
'(a15,i2)')
'# nst ', st%nst
2158 write(buf,
'(a15,i2)')
'# nspin ', st%d%nspin
2164 do is = 1, st%d%kpt%nglobal
2166 write(buf,
'(a,i4)')
'Eigenvalue ',ii
2175 do is = 1, st%d%kpt%nglobal
2185 do is = 1, st%d%kpt%nglobal
2197 type(c_ptr),
intent(inout) :: out_ionch
2198 class(
mesh_t),
intent(in) :: mesh
2200 integer,
intent(in) :: iter
2202 integer :: ii, ist, Nch, ik, idim
2203 character(len=68) :: buf
2204 real(real64),
allocatable :: ch(:), occ(:)
2205 real(real64),
allocatable :: occbuf(:)
2210 nch = st%nst * st%d%kpt%nglobal * st%d%dim
2211 safe_allocate(ch(0: nch))
2212 safe_allocate(occ(0: nch))
2218 do idim = 1, st%d%dim
2219 if (st%st_start <= ist .and. ist <= st%st_end .and. &
2220 st%d%kpt%start <= ik .and. ik <= st%d%kpt%end)
then
2221 occ(ii) = st%occ(ist, ik)
2229 if (st%parallel_in_states)
then
2230 safe_allocate(occbuf(0: nch))
2232 call st%mpi_grp%allreduce(occ(0), occbuf(0), nch+1, mpi_double_precision, mpi_sum)
2234 safe_deallocate_a(occbuf)
2242 safe_deallocate_a(ch)
2255 if (occ(ii)>
m_zero .or. ii == 0)
then
2256 write(buf,
'(a,f4.1,a)')
'Pion(',occ(ii)*ii,
'+, t)'
2266 if (occ(ii)>
m_zero .or. ii == 0)
then
2276 if (occ(ii)>
m_zero .or. ii == 0)
then
2282 safe_deallocate_a(ch)
2283 safe_deallocate_a(occ)
2289 subroutine td_write_proj(out_proj, space, mesh, ions, st, gs_st, kick, iter)
2290 type(c_ptr),
intent(inout) :: out_proj
2292 class(
mesh_t),
intent(in) :: mesh
2293 type(
ions_t),
intent(in) :: ions
2296 type(
kick_t),
intent(in) :: kick
2297 integer,
intent(in) :: iter
2299 complex(real64),
allocatable :: projections(:,:,:)
2300 character(len=80) :: aux
2301 integer :: ik, ist, uist, idir
2309 write(aux,
'(a15,i2)')
'# nspin ', st%d%nspin
2318 write(aux,
'(a,i8)')
"# nik ", st%nik
2322 write(aux,
'(a,2i8)')
"# st ", gs_st%st_start, st%nst
2326 write(aux,
'(a,2i8)')
"# ust ", gs_st%st_start, gs_st%st_end
2332 do ist = gs_st%st_start, st%nst
2340 do ist = gs_st%st_start, st%nst
2341 do uist = gs_st%st_start, gs_st%st_end
2342 write(aux,
'(i4,a,i4)') ist,
' -> ', uist
2353 if (.not. space%is_periodic())
then
2355 safe_allocate(projections(1:st%nst, gs_st%st_start:gs_st%st_end, 1:st%nik))
2356 do idir = 1, space%dim
2362 write(aux,
'(a,i1,a)')
"<i|x_", idir,
"|a>"
2366 do ist = gs_st%st_start, st%st_end
2367 do uist = gs_st%st_start, gs_st%st_end
2377 safe_deallocate_a(projections)
2387 safe_allocate(projections(1:st%nst, gs_st%st_start:gs_st%st_end, 1:st%nik))
2388 projections(:,:,:) =
m_z0
2394 do ist = gs_st%st_start, st%nst
2395 do uist = gs_st%st_start, gs_st%st_end
2404 safe_deallocate_a(projections)
2410 integer,
intent(in) :: dir
2412 integer :: uist, ist, ik, idim
2413 real(real64) :: n_dip(space%dim)
2414 complex(real64),
allocatable :: xpsi(:,:)
2415 complex(real64),
allocatable :: psi(:, :), gspsi(:, :)
2419 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
2420 safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
2421 safe_allocate(xpsi(1:mesh%np, 1:st%d%dim))
2423 do ik = st%d%kpt%start, st%d%kpt%end
2424 do ist = st%st_start, st%st_end
2426 do uist = gs_st%st_start, gs_st%st_end
2429 do idim = 1, st%d%dim
2430 xpsi(1:mesh%np, idim) = mesh%x(1:mesh%np, dir)*gspsi(1:mesh%np, idim)
2432 projections(ist, uist, ik) = -
zmf_dotp(mesh, st%d%dim, psi, xpsi, reduce = .false.)
2438 safe_deallocate_a(xpsi)
2439 safe_deallocate_a(gspsi)
2440 safe_deallocate_a(psi)
2445 n_dip = ions%dipole()
2447 do ist = gs_st%st_start, st%nst
2448 do uist = gs_st%st_start, gs_st%st_end
2449 projections(ist, uist, ik) = projections(ist, uist, ik) - n_dip(dir)
2465 subroutine td_write_n_ex(out_nex, outp, namespace, mesh, kpoints, st, gs_st, iter)
2466 type(c_ptr),
intent(inout) :: out_nex
2469 class(
mesh_t),
intent(in) :: mesh
2473 integer,
intent(in) :: iter
2475 complex(real64),
allocatable :: projections(:,:)
2476 character(len=80) :: aux, dir
2477 integer :: ik, ikpt, ist, uist, err
2478 real(real64) :: Nex, weight
2480 real(real64),
allocatable :: Nex_kpt(:)
2489 write(aux,
'(a15,i2)')
'# nspin ', st%d%nspin
2496 write(aux,
'(a,i8)')
"# nik ", st%nik
2500 write(aux,
'(a,2i8)')
"# st ", gs_st%st_start, st%nst
2504 write(aux,
'(a,2i8)')
"# ust ", gs_st%st_start, gs_st%st_end
2523 do ist = 1, gs_st%nst
2524 if (gs_st%occ(ist, ik) >
m_min_occ .and. ist > gs_nst) gs_nst = ist
2528 safe_allocate(projections(1:gs_nst, 1:st%nst))
2530 safe_allocate(nex_kpt(1:st%nik))
2532 do ik = st%d%kpt%start, st%d%kpt%end
2533 ikpt = st%d%get_kpoint_index(ik)
2536 weight = st%kweights(ik) * gs_st%occ(ist, ik)/ st%smear%el_per_state
2537 do uist = st%st_start, st%st_end
2538 nex_kpt(ikpt) = nex_kpt(ikpt) - weight * st%occ(uist, ik) * abs(projections(ist, uist))**2
2541 nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*st%kweights(ik)
2544 if (st%parallel_in_states .or. st%d%kpt%parallel)
then
2556 write(dir,
'(a,a,i7.7)') trim(outp%iter_dir),
"td.", iter
2559 + outp%how(option__output__density_kpt), dir,
"n_excited_el_kpt", namespace, &
2563 safe_deallocate_a(projections)
2564 safe_deallocate_a(nex_kpt)
2576 class(
mesh_t),
intent(in) :: mesh
2579 complex(real64),
intent(inout) :: projections(1:st%nst, gs_st%st_start:gs_st%nst, 1:st%nik)
2581 integer :: uist, ist, ik
2582 complex(real64),
allocatable :: psi(:, :), gspsi(:, :)
2585 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
2586 safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
2588 projections(:,:,:) =
m_zero
2590 do ik = st%d%kpt%start, st%d%kpt%end
2591 do ist = st%st_start, st%st_end
2593 do uist = gs_st%st_start, gs_st%nst
2595 projections(ist, uist, ik) =
zmf_dotp(mesh, st%d%dim, psi, gspsi, reduce = .false.)
2600 safe_deallocate_a(psi)
2601 safe_deallocate_a(gspsi)
2610 class(
mesh_t),
intent(in) :: mesh
2615 integer,
intent(in) :: iter
2617 complex(real64),
allocatable :: proj(:,:), psi(:,:,:), gs_psi(:,:,:), temp_state(:,:)
2618 character(len=80) :: filename1, filename2
2619 integer :: ik,ist, jst, file, idim, nk_proj
2623 write(filename1,
'(I10)') iter
2624 filename1 =
'td.general/projections_iter_'//trim(adjustl(filename1))
2627 safe_allocate(proj(1:gs_st%nst, 1:gs_st%nst))
2628 safe_allocate(psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2629 safe_allocate(gs_psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2630 safe_allocate(temp_state(1:mesh%np,1:gs_st%d%dim))
2636 nk_proj = kpoints%nik_skip
2638 do ik = kpoints%reduced%npoints-nk_proj+1, kpoints%reduced%npoints
2640 psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)=
m_zero
2641 gs_psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)=
m_zero
2644 write(filename2,
'(I10)') ik
2645 filename2 = trim(adjustl(filename1))//
'_ik_'//trim(adjustl(filename2))
2646 file =
io_open(filename2, namespace, action=
'write')
2649 do ist=gs_st%st_start,gs_st%st_end
2652 do idim = 1,gs_st%d%dim
2653 psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2656 do idim = 1,gs_st%d%dim
2657 gs_psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2666 assert(mesh%np_global*gs_st%d%dim < huge(0_int32))
2667 proj(1:gs_st%nst, 1:gs_st%nst) =
m_zero
2672 i8_to_i4(mesh%np_global*gs_st%d%dim), &
2673 cmplx(mesh%volume_element,
m_zero, real64) , &
2675 ubound(psi, dim = 1), &
2677 ubound(gs_psi, dim = 1), &
2680 ubound(proj, dim = 1))
2684 do ist = 1, gs_st%nst
2685 do jst = 1, gs_st%nst
2686 write(file,
'(I3,1x,I3,1x,e13.6,1x,e13.6,2x)') ist, jst, proj(ist,jst)
2694 safe_deallocate_a(proj)
2695 safe_deallocate_a(psi)
2696 safe_deallocate_a(gs_psi)
2697 safe_deallocate_a(temp_state)
2703 subroutine td_write_floquet(namespace, space, hm, ext_partners, gr, st, iter)
2704 type(namespace_t),
intent(in) :: namespace
2705 class(space_t),
intent(in) :: space
2706 type(hamiltonian_elec_t),
intent(inout) :: hm
2707 type(partner_list_t),
intent(in) :: ext_partners
2708 type(grid_t),
intent(in) :: gr
2709 type(states_elec_t),
intent(inout) :: st
2710 integer,
intent(in) :: iter
2712 complex(real64),
allocatable :: hmss(:,:), psi(:,:,:), hpsi(:,:,:), temp_state1(:,:)
2713 complex(real64),
allocatable :: HFloquet(:,:,:), HFloq_eff(:,:), temp(:,:)
2714 real(real64),
allocatable :: eigenval(:), bands(:,:)
2715 character(len=80) :: filename
2716 integer :: it, nT, ik, ist, in, im, file, idim, nik, ik_count
2717 integer :: Forder, Fdim, m0, n0, n1, nst, ii, jj, lim_nst
2718 logical :: downfolding
2719 type(states_elec_t) :: hm_st
2721 real(real64) :: dt, Tcycle, omega
2725 downfolding = .false.
2728 if (.not. iter == 0)
then
2736 assert(gr%np == gr%np_global)
2739 call states_elec_copy(hm_st, st)
2750 call parse_variable(namespace,
'TDFloquetFrequency', m_zero, omega, units_inp%energy)
2751 call messages_print_var_value(
'Frequency used for Floquet analysis', omega, namespace=namespace)
2752 if (abs(omega) <= m_epsilon)
then
2753 message(1) =
"Please give a non-zero value for TDFloquetFrequency"
2754 call messages_fatal(1, namespace=namespace)
2758 tcycle = m_two * m_pi / omega
2768 call parse_variable(namespace,
'TDFloquetSample',20 ,nt)
2769 call messages_print_var_value(
'Number of Floquet time-sampling points', nt, namespace=namespace)
2770 dt = tcycle/real(nt, real64)
2779 call parse_variable(namespace,
'TDFloquetDimension',-1,forder)
2780 if (forder .ge. 0)
then
2781 call messages_print_var_value(
'Order of multiphoton Floquet-Hamiltonian', forder, namespace=namespace)
2783 fdim = 2 * forder + 1
2785 message(1) =
'Floquet-Hamiltonian is downfolded'
2786 call messages_info(1, namespace=namespace)
2787 downfolding = .
true.
2792 dt = tcycle/real(nt, real64)
2795 nik = hm%kpoints%nik_skip
2797 safe_allocate(hmss(1:nst,1:nst))
2798 safe_allocate( psi(1:nst,1:st%d%dim,1:gr%np))
2799 safe_allocate(hpsi(1:nst,1:st%d%dim,1:gr%np))
2800 safe_allocate(temp_state1(1:gr%np,1:st%d%dim))
2808 safe_allocate(hfloquet(1:nik,1:nst*fdim, 1:nst*fdim))
2809 hfloquet(1:nik,1:nst*fdim, 1:nst*fdim) = m_zero
2814 call hm%update(gr, namespace, space, ext_partners, time=tcycle+it*dt)
2816 call zhamiltonian_elec_apply_all(hm, namespace, gr, st, hm_st)
2821 do ik = hm%kpoints%reduced%npoints-nik+1, hm%kpoints%reduced%npoints
2822 ik_count = ik_count + 1
2824 psi(1:nst, 1:st%d%dim, 1:gr%np)= m_zero
2825 hpsi(1:nst, 1:st%d%dim, 1:gr%np)= m_zero
2827 do ist = st%st_start, st%st_end
2828 if (state_kpt_is_local(st, ist, ik))
then
2829 call states_elec_get_state(st, gr, ist, ik,temp_state1)
2830 do idim = 1, st%d%dim
2831 psi(ist, idim, 1:gr%np) = temp_state1(1:gr%np, idim)
2833 call states_elec_get_state(hm_st, gr, ist, ik,temp_state1)
2834 do idim = 1, st%d%dim
2835 hpsi(ist, idim, 1:gr%np) = temp_state1(1:gr%np, idim)
2839 call comm_allreduce(mpi_world, psi)
2840 call comm_allreduce(mpi_world, hpsi)
2841 assert(gr%np_global*st%d%dim < huge(0_int32))
2842 hmss(1:nst,1:nst) = m_zero
2847 i8_to_i4(gr%np_global*st%d%dim), &
2848 cmplx(gr%volume_element, m_zero, real64) , &
2850 ubound(hpsi, dim = 1), &
2852 ubound(psi, dim = 1), &
2855 ubound(hmss, dim = 1))
2857 hmss(1:nst,1:nst) = conjg(hmss(1:nst,1:nst))
2860 do in = -forder, forder
2861 do im = -forder, forder
2862 ii = (in+forder) * nst
2863 jj = (im+forder) * nst
2864 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) = &
2865 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) + hmss(1:nst, 1:nst) *
exp(-(in-im)*m_zi*omega*it*dt)
2869 hfloquet(ik_count, ii+ist, ii+ist) = hfloquet(ik_count, ii+ist, ii+ist) + in*omega
2878 hfloquet(:,:,:) = m_one/nt*hfloquet(:,:,:)
2881 if (downfolding)
then
2883 safe_allocate(hfloq_eff(1:nst,1:nst))
2884 safe_allocate(eigenval(1:nst))
2885 safe_allocate(bands(1:nik,1:nst))
2887 hfloq_eff(1:nst,1:nst) = m_zero
2893 hfloq_eff(1:nst, 1:nst) = hfloquet(ik, n0+1:n0+nst, m0+1:m0+nst) + &
2894 m_one/omega*(matmul(hfloquet(ik, 1:nst, m0+1:m0+nst), hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst))- &
2895 matmul(hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst), hfloquet(ik, 1:nst, m0+1:m0+nst)))
2897 call lalg_eigensolve(nst, hfloq_eff, eigenval)
2898 bands(ik,1:nst) = eigenval(1:nst)
2900 safe_deallocate_a(hfloq_eff)
2903 safe_allocate(eigenval(1:nst*fdim))
2904 safe_allocate(bands(1:nik,1:nst*fdim))
2905 safe_allocate(temp(1:nst*fdim, 1:nst*fdim))
2908 temp(1:nst*fdim, 1:nst*fdim) = hfloquet(ik, 1:nst*fdim, 1:nst*fdim)
2909 call lalg_eigensolve(nst*fdim, temp, eigenval)
2910 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2915 if (downfolding)
then
2917 filename =
"downfolded_floquet_bands"
2920 filename =
"floquet_bands"
2923 if (mpi_world%is_root())
then
2925 file = io_open(filename, namespace, action =
'write')
2928 write(file,
'(e13.6, 1x)',advance=
'no') bands(ik,ist)
2935 if (.not. downfolding)
then
2938 bands(1:nik, 1:nst*fdim) = m_zero
2940 temp(1:nst*fdim,1:nst*fdim) = m_zero
2943 temp(ii+1:ii+nst, ii+1:ii+nst) = hfloquet(ik, ii+1:ii+nst, ii+1:ii+nst)
2945 call lalg_eigensolve(nst*fdim, temp, eigenval)
2946 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2949 if (mpi_world%is_root())
then
2950 filename =
'trivial_floquet_bands'
2951 file = io_open(filename, namespace, action =
'write')
2954 write(file,
'(e13.6, 1x)', advance=
'no') bands(ik,ist)
2963 call hm%update(gr, namespace, space, ext_partners, time=m_zero)
2965 safe_deallocate_a(hmss)
2966 safe_deallocate_a(psi)
2967 safe_deallocate_a(hpsi)
2968 safe_deallocate_a(temp_state1)
2969 safe_deallocate_a(hfloquet)
2970 safe_deallocate_a(eigenval)
2971 safe_deallocate_a(bands)
2972 safe_deallocate_a(temp)
2980 type(c_ptr),
intent(inout) :: out_total_current
2981 class(space_t),
intent(in) :: space
2982 class(mesh_t),
intent(in) :: mesh
2983 type(states_elec_t),
intent(in) :: st
2984 integer,
intent(in) :: iter
2986 integer :: idir, ispin
2987 character(len=50) :: aux
2988 real(real64) :: total_current(space%dim), abs_current(space%dim)
2992 if (mpi_world%is_root() .and. iter == 0)
then
2993 call td_write_print_header_init(out_total_current)
2996 call write_iter_header_start(out_total_current)
2998 do idir = 1, space%dim
2999 write(aux,
'(a2,a1,a1)')
'I(', index2axis(idir),
')'
3000 call write_iter_header(out_total_current, aux)
3003 do idir = 1, space%dim
3004 write(aux,
'(a10,a1,a1)')
'IntAbs(j)(', index2axis(idir),
')'
3005 call write_iter_header(out_total_current, aux)
3008 do ispin = 1, st%d%nspin
3009 do idir = 1, space%dim
3010 write(aux,
'(a4,i1,a1,a1,a1)')
'I-sp', ispin,
'(', index2axis(idir),
')'
3011 call write_iter_header(out_total_current, aux)
3015 call write_iter_nl(out_total_current)
3017 call td_write_print_header_end(out_total_current)
3020 assert(
allocated(st%current))
3022 if (mpi_world%is_root())
then
3023 call write_iter_start(out_total_current)
3026 total_current = 0.0_real64
3027 do idir = 1, space%dim
3028 do ispin = 1, st%d%spin_channels
3029 total_current(idir) = total_current(idir) + dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.)
3031 total_current(idir) = units_from_atomic(units_out%length/units_out%time, total_current(idir))
3033 call mesh%allreduce(total_current, dim = space%dim)
3035 abs_current = 0.0_real64
3036 do idir = 1, space%dim
3037 do ispin = 1, st%d%spin_channels
3038 abs_current(idir) = abs_current(idir) + dmf_integrate(mesh, abs(st%current(:, idir, ispin)), reduce = .false.)
3040 abs_current(idir) = units_from_atomic(units_out%length/units_out%time, abs_current(idir))
3042 call mesh%allreduce(abs_current, dim = space%dim)
3044 if (mpi_world%is_root())
then
3045 call write_iter_double(out_total_current, total_current, space%dim)
3046 call write_iter_double(out_total_current, abs_current, space%dim)
3049 do ispin = 1, st%d%nspin
3050 total_current = units_from_atomic(units_out%length/units_out%time, dmf_integrate(mesh, space%dim, st%current(:, :, ispin)))
3052 if (mpi_world%is_root())
then
3053 call write_iter_double(out_total_current, total_current, space%dim)
3057 if (mpi_world%is_root())
then
3058 call write_iter_nl(out_total_current)
3067 type(c_ptr),
intent(inout) :: write_obj
3068 class(space_t),
intent(in) :: space
3069 type(hamiltonian_elec_t),
intent(inout) :: hm
3070 type(grid_t),
intent(in) :: gr
3071 type(states_elec_t),
intent(in) :: st
3072 integer,
intent(in) :: iter
3074 integer :: idir, ispin
3075 character(len=50) :: aux
3076 real(real64),
allocatable :: heat_current(:, :, :)
3077 real(real64) :: total_current(space%dim)
3081 if (mpi_world%is_root() .and. iter == 0)
then
3082 call td_write_print_header_init(write_obj)
3085 call write_iter_header_start(write_obj)
3087 do idir = 1, space%dim
3088 write(aux,
'(a2,i1,a1)')
'Jh(', idir,
')'
3089 call write_iter_header(write_obj, aux)
3092 call write_iter_nl(write_obj)
3094 call td_write_print_header_end(write_obj)
3097 safe_allocate(heat_current(1:gr%np, 1:space%dim, 1:st%d%nspin))
3099 call current_heat_calculate(space, gr%der, hm, st, heat_current)
3101 if (mpi_world%is_root())
call write_iter_start(write_obj)
3103 total_current = 0.0_real64
3104 do idir = 1, space%dim
3105 do ispin = 1, st%d%spin_channels
3106 total_current(idir) = total_current(idir) + dmf_integrate(gr, heat_current(:, idir, ispin))
3108 total_current(idir) = units_from_atomic(units_out%energy*units_out%length/units_out%time, total_current(idir))
3111 safe_deallocate_a(heat_current)
3113 if (mpi_world%is_root())
call write_iter_double(write_obj, total_current, space%dim)
3115 if (mpi_world%is_root())
call write_iter_nl(write_obj)
3123 type(c_ptr),
intent(inout) :: out_partial_charges
3124 class(mesh_t),
intent(in) :: mesh
3125 type(states_elec_t),
intent(in) :: st
3126 type(ions_t),
intent(in) :: ions
3127 integer,
intent(in) :: iter
3130 character(len=50) :: aux
3131 real(real64),
allocatable :: hirshfeld_charges(:)
3135 safe_allocate(hirshfeld_charges(1:ions%natoms))
3137 call partial_charges_calculate(mesh, st, ions, hirshfeld_charges)
3139 if (mpi_world%is_root())
then
3143 call td_write_print_header_init(out_partial_charges)
3146 call write_iter_header_start(out_partial_charges)
3148 do idir = 1, ions%natoms
3149 write(aux,
'(a13,i3,a1)')
'hirshfeld(atom=', idir,
')'
3150 call write_iter_header(out_partial_charges, aux)
3153 call write_iter_nl(out_partial_charges)
3155 call td_write_print_header_end(out_partial_charges)
3158 call write_iter_start(out_partial_charges)
3160 call write_iter_double(out_partial_charges, hirshfeld_charges, ions%natoms)
3162 call write_iter_nl(out_partial_charges)
3165 safe_deallocate_a(hirshfeld_charges)
3171 subroutine td_write_q(out_q, space, ks, iter)
3172 type(c_ptr),
intent(inout) :: out_q
3173 class(space_t),
intent(in) :: space
3174 type(v_ks_t),
intent(in) :: ks
3175 integer,
intent(in) :: iter
3178 character(len=50) :: aux
3182 if (mpi_world%is_root())
then
3184 call td_write_print_header_init(out_q)
3185 call write_iter_header_start(out_q)
3186 do ii = 1, ks%pt%nmodes
3187 write(aux,
'(a1,i3,a3)')
'q', ii,
'(t)'
3188 call write_iter_header(out_q, aux)
3190 do ii = 1, ks%pt%nmodes
3191 write(aux,
'(a1,i3,a3)')
'p', ii,
'(t)'
3192 call write_iter_header(out_q, aux)
3194 do ii = 1, space%dim
3195 write(aux,
'(a3,i3,a3)')
'f_pt', ii,
'(t)'
3196 call write_iter_header(out_q, aux)
3198 call write_iter_nl(out_q)
3199 call td_write_print_header_end(out_q)
3202 call write_iter_start(out_q)
3203 call write_iter_double(out_q, ks%pt_mx%pt_q, ks%pt%nmodes)
3204 call write_iter_double(out_q, ks%pt_mx%pt_p, ks%pt%nmodes)
3205 call write_iter_double(out_q, ks%pt_mx%fmf, space%dim)
3206 call write_iter_nl(out_q)
3215 type(c_ptr),
intent(inout) :: out_mxll
3216 class(space_t),
intent(in) :: space
3217 type(hamiltonian_elec_t),
intent(in) :: hm
3218 real(real64),
intent(in) :: dt
3219 integer,
intent(in) :: iter
3222 real(real64) :: field(space%dim)
3223 character(len=80) :: aux
3224 character(len=1) :: field_char
3228 if (.not. mpi_world%is_root())
then
3234 call td_write_print_header_init(out_mxll)
3236 write(aux,
'(a7,e20.12,3a)')
'# dt = ', units_from_atomic(units_out%time, dt), &
3237 " [", trim(units_abbrev(units_out%time)),
"]"
3238 call write_iter_string(out_mxll, aux)
3239 call write_iter_nl(out_mxll)
3241 call write_iter_header_start(out_mxll)
3242 select case (hm%mxll%coupling_mode)
3243 case (length_gauge_dipole, multipolar_expansion)
3244 if (hm%mxll%add_electric_dip) field_char =
'E'
3245 if (hm%mxll%add_magnetic_dip) field_char =
'B'
3246 do idir = 1, space%dim
3247 write(aux,
'(a,i1,a)') field_char //
'(', idir,
')'
3248 call write_iter_header(out_mxll, aux)
3250 case (velocity_gauge_dipole)
3251 do idir = 1, space%dim
3252 write(aux,
'(a,i1,a)')
'A(', idir,
')'
3253 call write_iter_header(out_mxll, aux)
3256 call write_iter_nl(out_mxll)
3258 call write_iter_string(out_mxll,
'#[Iter n.]')
3259 call write_iter_header(out_mxll,
'[' // trim(units_abbrev(units_out%time)) //
']')
3263 select case (hm%mxll%coupling_mode)
3264 case (length_gauge_dipole, multipolar_expansion)
3265 if (hm%mxll%add_electric_dip) aux =
'[' // trim(units_abbrev(units_out%force)) //
']'
3266 if (hm%mxll%add_magnetic_dip) aux =
'[' // trim(units_abbrev(unit_one/units_out%length**2)) //
']'
3267 do idir = 1, space%dim
3268 call write_iter_header(out_mxll, aux)
3270 case (velocity_gauge_dipole)
3271 aux =
'[' // trim(units_abbrev(units_out%energy)) //
']'
3272 do idir = 1, space%dim
3273 call write_iter_header(out_mxll, aux)
3276 call write_iter_nl(out_mxll)
3277 call td_write_print_header_end(out_mxll)
3280 call write_iter_start(out_mxll)
3283 select case (hm%mxll%coupling_mode)
3284 case (length_gauge_dipole, multipolar_expansion)
3285 if (hm%mxll%add_electric_dip) field = units_from_atomic(units_out%force, hm%mxll%e_field_dip)
3286 if (hm%mxll%add_magnetic_dip) field = units_from_atomic(unit_one/units_out%length**2, hm%mxll%b_field_dip)
3287 call write_iter_double(out_mxll, field, space%dim)
3288 case (velocity_gauge_dipole)
3289 field = units_from_atomic(units_out%energy, hm%mxll%vec_pot_dip)
3290 call write_iter_double(out_mxll, field, space%dim)
3292 call write_iter_nl(out_mxll)
3300 type(c_ptr),
intent(inout) :: out_coords
3301 type(lda_u_t),
intent(in) :: lda_u
3302 integer,
intent(in) :: iter
3305 character(len=50) :: aux
3307 if (.not. mpi_world%is_root())
return
3312 call td_write_print_header_init(out_coords)
3315 call write_iter_header_start(out_coords)
3317 do ios = 1, lda_u%norbsets
3318 write(aux,
'(a2,i3,a1)')
'Ueff(', ios,
')'
3319 call write_iter_header(out_coords, aux)
3322 do ios = 1, lda_u%norbsets
3323 write(aux,
'(a2,i3,a1)')
'U(', ios,
')'
3324 call write_iter_header(out_coords, aux)
3327 do ios = 1, lda_u%norbsets
3328 write(aux,
'(a2,i3,a1)')
'J(', ios,
')'
3329 call write_iter_header(out_coords, aux)
3332 if (lda_u%intersite)
then
3333 do ios = 1, lda_u%norbsets
3334 do inn = 1, lda_u%orbsets(ios)%nneighbors
3335 write(aux,
'(a2,i3,a1,i3,a1)')
'V(', ios,
'-', inn,
')'
3336 call write_iter_header(out_coords, aux)
3342 call write_iter_nl(out_coords)
3345 call write_iter_string(out_coords,
'#[Iter n.]')
3346 call write_iter_header(out_coords,
'[' // trim(units_abbrev(units_out%time)) //
']')
3347 call write_iter_string(out_coords, &
3348 'Effective U ' // trim(units_abbrev(units_out%energy)) // &
3349 ', U in '// trim(units_abbrev(units_out%energy)) // &
3350 ', J in ' // trim(units_abbrev(units_out%energy)))
3351 call write_iter_nl(out_coords)
3353 call td_write_print_header_end(out_coords)
3356 call write_iter_start(out_coords)
3358 do ios = 1, lda_u%norbsets
3359 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3360 lda_u%orbsets(ios)%Ueff), 1)
3363 do ios = 1, lda_u%norbsets
3364 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3365 lda_u%orbsets(ios)%Ubar), 1)
3368 do ios = 1, lda_u%norbsets
3369 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3370 lda_u%orbsets(ios)%Jbar), 1)
3373 if (lda_u%intersite)
then
3374 do ios = 1, lda_u%norbsets
3375 do inn = 1, lda_u%orbsets(ios)%nneighbors
3376 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3377 lda_u%orbsets(ios)%V_ij(inn,0)), 1)
3382 call write_iter_nl(out_coords)
3389 type(c_ptr),
intent(inout) :: file_handle
3390 type(grid_t),
intent(in) :: grid
3391 type(kpoints_t),
intent(in) :: kpoints
3392 type(states_elec_t),
intent(in) :: st
3393 integer,
intent(in) :: iter
3395 integer :: ik_ispin, ist
3396 character(len=7) :: nkpt_str, nst_str
3397 character(len=7) :: ik_str, ist_str
3398 real(real64),
allocatable :: norm_ks(:, :)
3399 real(real64) :: n_electrons
3403 safe_allocate(norm_ks(1:st%nst, 1:st%nik))
3404 call states_elec_calc_norms(grid, kpoints, st, norm_ks)
3406 if (mpi_world%is_root())
then
3409 call td_write_print_header_init(file_handle)
3412 write(nkpt_str,
'(I7)') st%nik
3413 write(nst_str,
'(I7)') st%nst
3414 call write_iter_string(file_handle,
'# Dimensions. (nstates, nkpt * nspin):')
3415 call write_iter_string(file_handle, trim(adjustl(nst_str)) //
' ' // trim(adjustl(nkpt_str)))
3416 call write_iter_nl(file_handle)
3419 call write_iter_string(file_handle,
'# Norm ordering: (istate, ikpoint_spin)')
3420 call write_iter_nl(file_handle)
3423 call write_iter_header_start(file_handle)
3424 call write_iter_header(file_handle,
'N_electrons')
3425 do ik_ispin = 1, st%nik
3427 write(ik_str,
'(I7)') ik_ispin
3428 write(ist_str,
'(I7)') ist
3429 call write_iter_header(file_handle, &
3430 'Norm (' // trim(ist_str) //
',' // trim(ik_str) //
')')
3433 call write_iter_nl(file_handle)
3434 call td_write_print_header_end(file_handle)
3437 n_electrons = sum(st%occ * norm_ks**2)
3440 call write_iter_start(file_handle)
3441 call write_iter_double(file_handle, n_electrons, 1)
3442 do ik_ispin = 1, st%nik
3443 call write_iter_double(file_handle, norm_ks(:, ik_ispin),
size(norm_ks, 1))
3445 call write_iter_nl(file_handle)
3449 safe_deallocate_a(norm_ks)
3458 type(c_ptr),
intent(inout) :: file_handle
3459 type(ions_t),
intent(in) :: ions
3460 integer,
intent(in) :: iter
3463 real(real64) :: tmp(3)
3465 if (.not. mpi_world%is_root())
return
3469 assert(ions%space%dim == 3)
3472 call td_write_print_header_init(file_handle)
3475 call write_iter_header_start(file_handle)
3477 call write_iter_string(file_handle,
'# Iter, a, b, c, volume, alpha, beta, gamma, ' &
3478 //
'a_x, a_y, a_z, b_x, b_y, b_z, c_x, c_y, c_z')
3481 call write_iter_string(file_handle,
'#[Iter n.]')
3482 call write_iter_header(file_handle,
'[' // trim(units_abbrev(units_out%time)) //
']')
3483 call write_iter_string(file_handle, &
3484 'Lengths in ' // trim(units_abbrev(units_out%length)) // &
3485 ', Volume in ' // trim(units_abbrev(units_out%length**3)) // &
3486 ', Angles in degree, Lattice vectors in '// trim(units_abbrev(units_out%length)))
3487 call write_iter_nl(file_handle)
3489 call td_write_print_header_end(file_handle)
3492 call write_iter_start(file_handle)
3496 tmp(idir) = units_from_atomic(units_out%length, norm2(ions%latt%rlattice(1:3, idir)))
3498 call write_iter_double(file_handle, tmp, 3)
3501 tmp(1) = units_from_atomic(units_out%length**3, ions%latt%rcell_volume)
3502 call write_iter_double(file_handle, tmp(1), 1)
3505 call write_iter_double(file_handle, ions%latt%alpha, 1)
3506 call write_iter_double(file_handle, ions%latt%beta, 1)
3507 call write_iter_double(file_handle, ions%latt%gamma, 1)
3511 tmp(1:3) = units_from_atomic(units_out%length, ions%latt%rlattice(:, idir))
3512 call write_iter_double(file_handle, tmp, 3)
3514 call write_iter_nl(file_handle)
3523 type(namespace_t),
intent(in) :: namespace
3524 integer,
intent(in) :: iter
3525 real(real64),
intent(in) :: dt
3527 integer :: default, flags, iout, first
3582 call parse_variable(namespace,
'MaxwellTDOutput', default, flags)
3584 if (.not. varinfo_valid_option(
'MaxwellTDOutput', flags, is_flag = .
true.))
then
3585 call messages_input_error(namespace,
'MaxwellTDOutput')
3589 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3590 if (writ%out(iout)%write)
then
3591 writ%out(iout + 1)%write = .
true.
3592 writ%out(iout + 2)%write = .
true.
3597 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3606 call io_mkdir(
'td.general', namespace)
3611 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_x", namespace)))
3613 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_y", namespace)))
3615 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_z", namespace)))
3621 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_x", namespace)))
3623 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_y", namespace)))
3625 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_z", namespace)))
3631 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_x", namespace)))
3633 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_y", namespace)))
3635 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_z", namespace)))
3641 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_x", namespace)))
3643 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_y", namespace)))
3645 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_z", namespace)))
3651 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_x", namespace)))
3653 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_y", namespace)))
3655 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_z", namespace)))
3661 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_x", namespace)))
3663 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_y", namespace)))
3665 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_z", namespace)))
3671 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/maxwell_energy", namespace)))
3676 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-x", namespace)))
3681 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-y", namespace)))
3686 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-z", namespace)))
3691 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-x", namespace)))
3696 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-y", namespace)))
3701 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-z", namespace)))
3717 if (writ%out(iout)%write)
call write_iter_end(writ%out(iout)%handle)
3725 subroutine td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
3727 class(space_t),
intent(in) :: space
3728 type(grid_t),
intent(inout) :: gr
3729 type(states_mxll_t),
intent(inout) :: st
3730 type(hamiltonian_mxll_t),
intent(inout) :: hm
3731 type(helmholtz_decomposition_t),
intent(inout) :: helmholtz
3732 real(real64),
intent(in) :: dt
3733 integer,
intent(in) :: iter
3734 type(namespace_t),
intent(in) :: namespace
3739 call profiling_in(
"TD_WRITE_ITER_MAXWELL")
3742 call helmholtz%get_trans_field(namespace, st%rs_state_trans, total_field=st%rs_state)
3743 call get_rs_state_at_point(st%selected_points_rs_state_trans(:,:), st%rs_state_trans, &
3744 st%selected_points_coordinate(:,:), st, gr)
3747 hm%energy%energy_trans = m_zero
3751 call helmholtz%get_long_field(namespace, st%rs_state_long, total_field=st%rs_state)
3752 call get_rs_state_at_point(st%selected_points_rs_state_long(:,:), st%rs_state_long, &
3753 st%selected_points_coordinate(:,:), st, gr)
3756 hm%energy%energy_long = m_zero
3848 call profiling_out(
"TD_WRITE_ITER_MAXWELL")
3857 type(hamiltonian_mxll_t),
intent(in) :: hm
3858 integer,
intent(in) :: iter
3862 integer :: n_columns
3864 if (.not. mpi_world%is_root())
return
3889 call write_iter_header(
out_maxwell_energy,
'[' // trim(units_abbrev(units_out%time)) //
']')
3891 do ii = 1, n_columns
3892 call write_iter_header(
out_maxwell_energy,
'[' // trim(units_abbrev(units_out%energy)) //
']')
3900 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy), 1)
3901 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%e_energy), 1)
3902 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%b_energy), 1)
3904 hm%energy%energy+hm%energy%boundaries), 1)
3905 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%boundaries), 1)
3906 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_trans), 1)
3907 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_long), 1)
3908 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_plane_waves), 1)
3917 type(c_ptr),
intent(inout) :: out_field_surf
3918 type(states_mxll_t),
intent(in) :: st
3919 integer,
intent(in) :: dim
3920 integer,
intent(in) :: iter
3924 integer :: n_columns
3926 if (.not. mpi_world%is_root())
return
3933 call td_write_print_header_init(out_field_surf)
3936 call write_iter_header_start(out_field_surf)
3937 call write_iter_header(out_field_surf,
'- x direction')
3938 call write_iter_header(out_field_surf,
'+ x direction')
3939 call write_iter_header(out_field_surf,
'- y direction')
3940 call write_iter_header(out_field_surf,
'+ y direction')
3941 call write_iter_header(out_field_surf,
'- z direction')
3942 call write_iter_header(out_field_surf,
'+ z direction')
3943 call write_iter_header(out_field_surf,
'- x dir. p. w.')
3944 call write_iter_header(out_field_surf,
'+ x dir. p. w.')
3945 call write_iter_header(out_field_surf,
'- y dir. p. w.')
3946 call write_iter_header(out_field_surf,
'+ y dir. p. w.')
3947 call write_iter_header(out_field_surf,
'- z dir. p. w.')
3948 call write_iter_header(out_field_surf,
'+ z dir. p. w.')
3950 call write_iter_nl(out_field_surf)
3953 call write_iter_string(out_field_surf,
'#[Iter n.]')
3954 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%time)) //
']')
3956 do ii = 1, n_columns
3957 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%energy/units_out%length)) //
']')
3959 call write_iter_nl(out_field_surf)
3961 call td_write_print_header_end(out_field_surf)
3964 call write_iter_start(out_field_surf)
3965 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3966 st%electric_field_box_surface(1,1,dim)), 1)
3967 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3968 st%electric_field_box_surface(2,1,dim)), 1)
3969 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3970 st%electric_field_box_surface(1,2,dim)), 1)
3971 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3972 st%electric_field_box_surface(2,2,dim)), 1)
3973 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3974 st%electric_field_box_surface(1,3,dim)), 1)
3975 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3976 st%electric_field_box_surface(2,3,dim)), 1)
3977 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3978 st%electric_field_box_surface_plane_waves(1,1,dim)), 1)
3979 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3980 st%electric_field_box_surface_plane_waves(2,1,dim)), 1)
3981 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3982 st%electric_field_box_surface_plane_waves(1,2,dim)), 1)
3983 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3984 st%electric_field_box_surface_plane_waves(2,2,dim)), 1)
3985 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3986 st%electric_field_box_surface_plane_waves(1,3,dim)), 1)
3987 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3988 st%electric_field_box_surface_plane_waves(2,3,dim)), 1)
3989 call write_iter_nl(out_field_surf)
3997 type(c_ptr),
intent(inout) :: out_field_surf
3998 type(states_mxll_t),
intent(in) :: st
3999 integer,
intent(in) :: dim
4000 integer,
intent(in) :: iter
4004 integer :: n_columns
4006 if (.not. mpi_world%is_root())
return
4013 call td_write_print_header_init(out_field_surf)
4016 call write_iter_header_start(out_field_surf)
4017 call write_iter_header(out_field_surf,
'- x direction')
4018 call write_iter_header(out_field_surf,
'+ x direction')
4019 call write_iter_header(out_field_surf,
'- y direction')
4020 call write_iter_header(out_field_surf,
'+ y direction')
4021 call write_iter_header(out_field_surf,
'- z direction')
4022 call write_iter_header(out_field_surf,
'+ z direction')
4023 call write_iter_header(out_field_surf,
'- x dir. p. w.')
4024 call write_iter_header(out_field_surf,
'+ x dir. p. w.')
4025 call write_iter_header(out_field_surf,
'- y dir. p. w.')
4026 call write_iter_header(out_field_surf,
'+ y dir. p. w.')
4027 call write_iter_header(out_field_surf,
'- z dir. p. w.')
4028 call write_iter_header(out_field_surf,
'+ z dir. p. w.')
4030 call write_iter_nl(out_field_surf)
4033 call write_iter_string(out_field_surf,
'#[Iter n.]')
4034 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%time)) //
']')
4036 do ii = 1, n_columns
4037 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(unit_one/units_out%length**2)) //
']')
4039 call write_iter_nl(out_field_surf)
4041 call td_write_print_header_end(out_field_surf)
4044 call write_iter_start(out_field_surf)
4045 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4046 st%magnetic_field_box_surface(1,1,dim)), 1)
4047 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4048 st%magnetic_field_box_surface(2,1,dim)), 1)
4049 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4050 st%magnetic_field_box_surface(1,2,dim)), 1)
4051 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4052 st%magnetic_field_box_surface(2,2,dim)), 1)
4053 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4054 st%magnetic_field_box_surface(1,3,dim)), 1)
4055 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4056 st%magnetic_field_box_surface(2,3,dim)), 1)
4057 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4058 st%magnetic_field_box_surface_plane_waves(1,1,dim)), 1)
4059 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4060 st%magnetic_field_box_surface_plane_waves(2,1,dim)), 1)
4061 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4062 st%magnetic_field_box_surface_plane_waves(1,2,dim)), 1)
4063 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4064 st%magnetic_field_box_surface_plane_waves(2,2,dim)), 1)
4065 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4066 st%magnetic_field_box_surface_plane_waves(1,3,dim)), 1)
4067 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4068 st%magnetic_field_box_surface_plane_waves(2,3,dim)), 1)
4069 call write_iter_nl(out_field_surf)
4075 subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
4076 type(c_ptr),
intent(inout) :: out_fields
4077 class(space_t),
intent(in) :: space
4078 type(states_mxll_t),
intent(in) :: st
4079 integer,
intent(in) :: iter
4080 real(real64),
intent(in) :: dt
4081 integer,
intent(in) :: e_or_b_field
4082 integer,
intent(in) :: field_type
4083 integer,
intent(in) :: idir
4087 real(real64) :: field(space%dim), selected_field
4088 character(len=80) :: aux
4090 if (.not. mpi_world%is_root())
return
4095 call td_write_print_header_init(out_fields)
4098 write(aux,
'(a7,e20.12,3a)')
'# dt = ', units_from_atomic(units_out%time, dt), &
4099 " [", trim(units_abbrev(units_out%time)),
"]"
4100 call write_iter_string(out_fields, aux)
4101 call write_iter_nl(out_fields)
4103 call write_iter_header_start(out_fields)
4105 do id = 1, st%selected_points_number
4106 select case (e_or_b_field)
4108 write(aux,
'(a,i1,a)')
'E(', id,
')'
4110 write(aux,
'(a,i1,a)')
'B(', id,
')'
4112 call write_iter_header(out_fields, aux)
4115 call write_iter_nl(out_fields)
4116 call write_iter_string(out_fields,
'#[Iter n.]')
4117 call write_iter_header(out_fields,
' [' // trim(units_abbrev(units_out%time)) //
']')
4122 aux =
' [' // trim(units_abbrev(units_out%force)) //
']'
4123 do id = 1, st%selected_points_number
4124 call write_iter_header(out_fields, aux)
4126 call write_iter_nl(out_fields)
4127 call td_write_print_header_end(out_fields)
4130 call write_iter_start(out_fields)
4132 do id = 1, st%selected_points_number
4133 select case (e_or_b_field)
4136 select case (field_type)
4138 call get_electric_field_vector(st%selected_points_rs_state(:,id), field(1:st%dim))
4140 call get_electric_field_vector(st%selected_points_rs_state_long(:,id), field(1:st%dim))
4142 call get_electric_field_vector(st%selected_points_rs_state_trans(:,id), field(1:st%dim))
4144 selected_field = units_from_atomic(units_out%energy/units_out%length, field(idir))
4147 select case (field_type)
4149 call get_magnetic_field_vector(st%selected_points_rs_state(:,id), st%rs_sign, field(1:st%dim))
4151 call get_magnetic_field_vector(st%selected_points_rs_state_long(:,id), st%rs_sign, field(1:st%dim))
4153 call get_magnetic_field_vector(st%selected_points_rs_state_trans(:,id), st%rs_sign, field(1:st%dim))
4155 selected_field = units_from_atomic(unit_one/units_out%length**2, field(idir))
4157 call write_iter_double(out_fields, selected_field, 1)
4160 call write_iter_nl(out_fields)
4168 type(namespace_t),
intent(in) :: namespace
4169 class(space_t),
intent(in) :: space
4170 type(grid_t),
intent(inout) :: gr
4171 type(states_mxll_t),
intent(inout) :: st
4172 type(hamiltonian_mxll_t),
intent(inout) :: hm
4173 type(helmholtz_decomposition_t),
intent(inout) :: helmholtz
4174 type(output_t),
intent(in) :: outp
4175 integer,
intent(in) :: iter
4176 real(real64),
intent(in) :: time
4178 character(len=256) :: filename
4180 push_sub(td_write_maxwell_free_data)
4181 call profiling_in(
"TD_WRITE_MAXWELL_DATA")
4184 write(filename,
'(a,a,i7.7)') trim(outp%iter_dir),
"td.", iter
4186 call output_mxll(outp, namespace, space, gr, st, hm, helmholtz, time, filename)
4188 call profiling_out(
"TD_WRITE_MAXWELL_DATA")
4189 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)
construct path name from given name and 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.
subroutine, public dmf_multipoles(mesh, ff, lmax, multipole, mask)
This routine calculates the multipoles of a function ff.
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.
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)
subroutine, public parse_block_string(blk, l, c, res, convert_to_c)
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
integer, parameter, public restart_proj
integer, parameter, public restart_type_load
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
integer, parameter, public out_cell_parameters
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)
subroutine td_write_floquet(namespace, space, hm, ext_partners, gr, st, 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
subroutine, public td_write_mxll_free_data(namespace, space, gr, st, hm, helmholtz, outp, iter, time)
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 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)
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_cell_parameters(file_handle, ions, iter)
Write the cell parameters as a function of time.
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.
subroutine, public write_iter_header(out, string)
subroutine, public write_iter_string(out, string)
subroutine, public write_iter_init(out, iter, factor, file)
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)