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 :: &
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
241 subroutine td_write_kick(outp, namespace, space, mesh, kick, ions, iter)
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))
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
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
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)
1238 type(
grid_t),
intent(in) :: gr
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
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)
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(:)
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)
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)
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)
2290 type(c_ptr),
intent(inout) :: out_proj
2291 class(
space_t),
intent(in) :: space
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
2521 if (.not. st%parallel_in_states)
then
2522 do ik = st%d%kpt%start, st%d%kpt%end
2523 do ist = 1, gs_st%nst
2524 if (gs_st%occ(ist, ik) >
m_min_occ) gs_nst = ist
2531 safe_allocate(projections(1:gs_nst, 1:st%nst))
2533 safe_allocate(nex_kpt(1:st%nik))
2535 do ik = st%d%kpt%start, st%d%kpt%end
2536 ikpt = st%d%get_kpoint_index(ik)
2539 weight = st%kweights(ik) * gs_st%occ(ist, ik)/ st%smear%el_per_state
2540 do uist = st%st_start, st%st_end
2541 nex_kpt(ikpt) = nex_kpt(ikpt) - weight * st%occ(uist, ik) * abs(projections(ist, uist))**2
2545 nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*
m_half*st%kweights(ik)
2547 nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*st%kweights(ik)
2551 if (st%parallel_in_states .or. st%d%kpt%parallel)
then
2563 write(dir,
'(a,a,i7.7)') trim(outp%iter_dir),
"td.", iter
2566 + outp%how(option__output__density_kpt), dir,
"n_excited_el_kpt", namespace, &
2570 safe_deallocate_a(projections)
2571 safe_deallocate_a(nex_kpt)
2583 class(
mesh_t),
intent(in) :: mesh
2586 complex(real64),
intent(inout) :: projections(1:st%nst, gs_st%st_start:gs_st%nst, 1:st%nik)
2588 integer :: uist, ist, ik
2589 complex(real64),
allocatable :: psi(:, :), gspsi(:, :)
2592 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
2593 safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
2595 projections(:,:,:) =
m_zero
2597 do ik = st%d%kpt%start, st%d%kpt%end
2598 do ist = st%st_start, st%st_end
2600 do uist = gs_st%st_start, gs_st%nst
2602 projections(ist, uist, ik) =
zmf_dotp(mesh, st%d%dim, psi, gspsi, reduce = .false.)
2607 safe_deallocate_a(psi)
2608 safe_deallocate_a(gspsi)
2617 class(
mesh_t),
intent(in) :: mesh
2622 integer,
intent(in) :: iter
2624 complex(real64),
allocatable :: proj(:,:), psi(:,:,:), gs_psi(:,:,:), temp_state(:,:)
2625 character(len=80) :: filename1, filename2
2626 integer :: ik,ist, jst, file, idim, nk_proj
2630 write(filename1,
'(I10)') iter
2631 filename1 =
'td.general/projections_iter_'//trim(adjustl(filename1))
2634 safe_allocate(proj(1:gs_st%nst, 1:gs_st%nst))
2635 safe_allocate(psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2636 safe_allocate(gs_psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2637 safe_allocate(temp_state(1:mesh%np,1:gs_st%d%dim))
2643 nk_proj = kpoints%nik_skip
2645 do ik = kpoints%reduced%npoints-nk_proj+1, kpoints%reduced%npoints
2647 psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)=
m_zero
2648 gs_psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)=
m_zero
2651 write(filename2,
'(I10)') ik
2652 filename2 = trim(adjustl(filename1))//
'_ik_'//trim(adjustl(filename2))
2653 file =
io_open(filename2, namespace, action=
'write')
2656 do ist=gs_st%st_start,gs_st%st_end
2659 do idim = 1,gs_st%d%dim
2660 psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2663 do idim = 1,gs_st%d%dim
2664 gs_psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2673 assert(mesh%np_global*gs_st%d%dim < huge(0_int32))
2674 proj(1:gs_st%nst, 1:gs_st%nst) =
m_zero
2679 i8_to_i4(mesh%np_global*gs_st%d%dim), &
2680 cmplx(mesh%volume_element,
m_zero, real64) , &
2682 ubound(psi, dim = 1), &
2684 ubound(gs_psi, dim = 1), &
2687 ubound(proj, dim = 1))
2691 do ist = 1, gs_st%nst
2692 do jst = 1, gs_st%nst
2693 write(file,
'(I3,1x,I3,1x,e13.6,1x,e13.6,2x)') ist, jst, proj(ist,jst)
2701 safe_deallocate_a(proj)
2702 safe_deallocate_a(psi)
2703 safe_deallocate_a(gs_psi)
2704 safe_deallocate_a(temp_state)
2710 subroutine td_write_floquet(namespace, space, hm, ext_partners, mesh, st, iter)
2711 type(namespace_t),
intent(in) :: namespace
2712 class(space_t),
intent(in) :: space
2713 type(hamiltonian_elec_t),
intent(inout) :: hm
2714 type(partner_list_t),
intent(in) :: ext_partners
2715 class(mesh_t),
intent(in) :: mesh
2716 type(states_elec_t),
intent(inout) :: st
2717 integer,
intent(in) :: iter
2719 complex(real64),
allocatable :: hmss(:,:), psi(:,:,:), hpsi(:,:,:), temp_state1(:,:)
2720 complex(real64),
allocatable :: HFloquet(:,:,:), HFloq_eff(:,:), temp(:,:)
2721 real(real64),
allocatable :: eigenval(:), bands(:,:)
2722 character(len=80) :: filename
2723 integer :: it, nT, ik, ist, in, im, file, idim, nik, ik_count
2724 integer :: Forder, Fdim, m0, n0, n1, nst, ii, jj, lim_nst
2725 logical :: downfolding
2726 type(states_elec_t) :: hm_st
2728 real(real64) :: dt, Tcycle, omega
2732 downfolding = .false.
2735 if (.not. iter == 0)
then
2743 assert(mesh%np == mesh%np_global)
2746 call states_elec_copy(hm_st, st)
2757 call parse_variable(namespace,
'TDFloquetFrequency', m_zero, omega, units_inp%energy)
2758 call messages_print_var_value(
'Frequency used for Floquet analysis', omega, namespace=namespace)
2759 if (abs(omega) <= m_epsilon)
then
2760 message(1) =
"Please give a non-zero value for TDFloquetFrequency"
2761 call messages_fatal(1, namespace=namespace)
2765 tcycle = m_two * m_pi / omega
2775 call parse_variable(namespace,
'TDFloquetSample',20 ,nt)
2776 call messages_print_var_value(
'Number of Floquet time-sampling points', nt, namespace=namespace)
2777 dt = tcycle/real(nt, real64)
2786 call parse_variable(namespace,
'TDFloquetDimension',-1,forder)
2787 if (forder .ge. 0)
then
2788 call messages_print_var_value(
'Order of multiphoton Floquet-Hamiltonian', forder, namespace=namespace)
2790 fdim = 2 * forder + 1
2792 message(1) =
'Floquet-Hamiltonian is downfolded'
2793 call messages_info(1, namespace=namespace)
2794 downfolding = .
true.
2799 dt = tcycle/real(nt, real64)
2802 nik = hm%kpoints%nik_skip
2804 safe_allocate(hmss(1:nst,1:nst))
2805 safe_allocate( psi(1:nst,1:st%d%dim,1:mesh%np))
2806 safe_allocate(hpsi(1:nst,1:st%d%dim,1:mesh%np))
2807 safe_allocate(temp_state1(1:mesh%np,1:st%d%dim))
2815 safe_allocate(hfloquet(1:nik,1:nst*fdim, 1:nst*fdim))
2816 hfloquet(1:nik,1:nst*fdim, 1:nst*fdim) = m_zero
2821 call hm%update(mesh, namespace, space, ext_partners, time=tcycle+it*dt)
2823 call zhamiltonian_elec_apply_all(hm, namespace, mesh, st, hm_st)
2828 do ik = hm%kpoints%reduced%npoints-nik+1, hm%kpoints%reduced%npoints
2829 ik_count = ik_count + 1
2831 psi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
2832 hpsi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
2834 do ist = st%st_start, st%st_end
2835 if (state_kpt_is_local(st, ist, ik))
then
2836 call states_elec_get_state(st, mesh, ist, ik,temp_state1)
2837 do idim = 1, st%d%dim
2838 psi(ist, idim, 1:mesh%np) = temp_state1(1:mesh%np,idim)
2840 call states_elec_get_state(hm_st, mesh, ist, ik,temp_state1)
2841 do idim = 1, st%d%dim
2842 hpsi(ist, idim,1:mesh%np) = temp_state1(1:mesh%np,idim)
2846 call comm_allreduce(mpi_world, psi)
2847 call comm_allreduce(mpi_world, hpsi)
2848 assert(mesh%np_global*st%d%dim < huge(0_int32))
2849 hmss(1:nst,1:nst) = m_zero
2854 i8_to_i4(mesh%np_global*st%d%dim), &
2855 cmplx(mesh%volume_element, m_zero, real64) , &
2857 ubound(hpsi, dim = 1), &
2859 ubound(psi, dim = 1), &
2862 ubound(hmss, dim = 1))
2864 hmss(1:nst,1:nst) = conjg(hmss(1:nst,1:nst))
2867 do in = -forder, forder
2868 do im = -forder, forder
2869 ii = (in+forder) * nst
2870 jj = (im+forder) * nst
2871 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) = &
2872 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) + hmss(1:nst, 1:nst) *
exp(-(in-im)*m_zi*omega*it*dt)
2876 hfloquet(ik_count, ii+ist, ii+ist) = hfloquet(ik_count, ii+ist, ii+ist) + in*omega
2885 hfloquet(:,:,:) = m_one/nt*hfloquet(:,:,:)
2888 if (downfolding)
then
2890 safe_allocate(hfloq_eff(1:nst,1:nst))
2891 safe_allocate(eigenval(1:nst))
2892 safe_allocate(bands(1:nik,1:nst))
2894 hfloq_eff(1:nst,1:nst) = m_zero
2900 hfloq_eff(1:nst, 1:nst) = hfloquet(ik, n0+1:n0+nst, m0+1:m0+nst) + &
2901 m_one/omega*(matmul(hfloquet(ik, 1:nst, m0+1:m0+nst), hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst))- &
2902 matmul(hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst), hfloquet(ik, 1:nst, m0+1:m0+nst)))
2904 call lalg_eigensolve(nst, hfloq_eff, eigenval)
2905 bands(ik,1:nst) = eigenval(1:nst)
2907 safe_deallocate_a(hfloq_eff)
2910 safe_allocate(eigenval(1:nst*fdim))
2911 safe_allocate(bands(1:nik,1:nst*fdim))
2912 safe_allocate(temp(1:nst*fdim, 1:nst*fdim))
2915 temp(1:nst*fdim, 1:nst*fdim) = hfloquet(ik, 1:nst*fdim, 1:nst*fdim)
2916 call lalg_eigensolve(nst*fdim, temp, eigenval)
2917 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2922 if (downfolding)
then
2924 filename =
"downfolded_floquet_bands"
2927 filename =
"floquet_bands"
2930 if (mpi_world%rank == 0)
then
2932 file = io_open(filename, namespace, action =
'write')
2935 write(file,
'(e13.6, 1x)',advance=
'no') bands(ik,ist)
2942 if (.not. downfolding)
then
2945 bands(1:nik, 1:nst*fdim) = m_zero
2947 temp(1:nst*fdim,1:nst*fdim) = m_zero
2950 temp(ii+1:ii+nst, ii+1:ii+nst) = hfloquet(ik, ii+1:ii+nst, ii+1:ii+nst)
2952 call lalg_eigensolve(nst*fdim, temp, eigenval)
2953 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2956 if (mpi_world%rank == 0)
then
2957 filename =
'trivial_floquet_bands'
2958 file = io_open(filename, namespace, action =
'write')
2961 write(file,
'(e13.6, 1x)', advance=
'no') bands(ik,ist)
2970 call hm%update(mesh, namespace, space, ext_partners, time=m_zero)
2972 safe_deallocate_a(hmss)
2973 safe_deallocate_a(psi)
2974 safe_deallocate_a(hpsi)
2975 safe_deallocate_a(temp_state1)
2976 safe_deallocate_a(hfloquet)
2977 safe_deallocate_a(eigenval)
2978 safe_deallocate_a(bands)
2979 safe_deallocate_a(temp)
2987 type(c_ptr),
intent(inout) :: out_total_current
2988 class(space_t),
intent(in) :: space
2989 class(mesh_t),
intent(in) :: mesh
2990 type(states_elec_t),
intent(in) :: st
2991 integer,
intent(in) :: iter
2993 integer :: idir, ispin
2994 character(len=50) :: aux
2995 real(real64) :: total_current(space%dim), abs_current(space%dim)
2999 if (mpi_grp_is_root(mpi_world) .and. iter == 0)
then
3000 call td_write_print_header_init(out_total_current)
3003 call write_iter_header_start(out_total_current)
3005 do idir = 1, space%dim
3006 write(aux,
'(a2,a1,a1)')
'I(', index2axis(idir),
')'
3007 call write_iter_header(out_total_current, aux)
3010 do idir = 1, space%dim
3011 write(aux,
'(a10,a1,a1)')
'IntAbs(j)(', index2axis(idir),
')'
3012 call write_iter_header(out_total_current, aux)
3015 do ispin = 1, st%d%nspin
3016 do idir = 1, space%dim
3017 write(aux,
'(a4,i1,a1,a1,a1)')
'I-sp', ispin,
'(', index2axis(idir),
')'
3018 call write_iter_header(out_total_current, aux)
3022 call write_iter_nl(out_total_current)
3024 call td_write_print_header_end(out_total_current)
3027 assert(
allocated(st%current))
3029 if (mpi_grp_is_root(mpi_world))
then
3030 call write_iter_start(out_total_current)
3033 total_current = 0.0_real64
3034 do idir = 1, space%dim
3035 do ispin = 1, st%d%spin_channels
3036 total_current(idir) = total_current(idir) + dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.)
3038 total_current(idir) = units_from_atomic(units_out%length/units_out%time, total_current(idir))
3040 call mesh%allreduce(total_current, dim = space%dim)
3042 abs_current = 0.0_real64
3043 do idir = 1, space%dim
3044 do ispin = 1, st%d%spin_channels
3045 abs_current(idir) = abs_current(idir) + dmf_integrate(mesh, abs(st%current(:, idir, ispin)), reduce = .false.)
3047 abs_current(idir) = units_from_atomic(units_out%length/units_out%time, abs_current(idir))
3049 call mesh%allreduce(abs_current, dim = space%dim)
3051 if (mpi_grp_is_root(mpi_world))
then
3052 call write_iter_double(out_total_current, total_current, space%dim)
3053 call write_iter_double(out_total_current, abs_current, space%dim)
3056 do ispin = 1, st%d%nspin
3057 total_current = units_from_atomic(units_out%length/units_out%time, dmf_integrate(mesh, space%dim, st%current(:, :, ispin)))
3059 if (mpi_grp_is_root(mpi_world))
then
3060 call write_iter_double(out_total_current, total_current, space%dim)
3064 if (mpi_grp_is_root(mpi_world))
then
3065 call write_iter_nl(out_total_current)
3074 type(c_ptr),
intent(inout) :: write_obj
3075 class(space_t),
intent(in) :: space
3076 type(hamiltonian_elec_t),
intent(inout) :: hm
3077 type(grid_t),
intent(in) :: gr
3078 type(states_elec_t),
intent(in) :: st
3079 integer,
intent(in) :: iter
3081 integer :: idir, ispin
3082 character(len=50) :: aux
3083 real(real64),
allocatable :: heat_current(:, :, :)
3084 real(real64) :: total_current(space%dim)
3088 if (mpi_grp_is_root(mpi_world) .and. iter == 0)
then
3089 call td_write_print_header_init(write_obj)
3092 call write_iter_header_start(write_obj)
3094 do idir = 1, space%dim
3095 write(aux,
'(a2,i1,a1)')
'Jh(', idir,
')'
3096 call write_iter_header(write_obj, aux)
3099 call write_iter_nl(write_obj)
3101 call td_write_print_header_end(write_obj)
3104 safe_allocate(heat_current(1:gr%np, 1:space%dim, 1:st%d%nspin))
3106 call current_heat_calculate(space, gr%der, hm, st, heat_current)
3108 if (mpi_grp_is_root(mpi_world))
call write_iter_start(write_obj)
3110 total_current = 0.0_real64
3111 do idir = 1, space%dim
3112 do ispin = 1, st%d%spin_channels
3113 total_current(idir) = total_current(idir) + dmf_integrate(gr, heat_current(:, idir, ispin))
3115 total_current(idir) = units_from_atomic(units_out%energy*units_out%length/units_out%time, total_current(idir))
3118 safe_deallocate_a(heat_current)
3120 if (mpi_grp_is_root(mpi_world))
call write_iter_double(write_obj, total_current, space%dim)
3122 if (mpi_grp_is_root(mpi_world))
call write_iter_nl(write_obj)
3130 type(c_ptr),
intent(inout) :: out_partial_charges
3131 class(mesh_t),
intent(in) :: mesh
3132 type(states_elec_t),
intent(in) :: st
3133 type(ions_t),
intent(in) :: ions
3134 integer,
intent(in) :: iter
3137 character(len=50) :: aux
3138 real(real64),
allocatable :: hirshfeld_charges(:)
3142 safe_allocate(hirshfeld_charges(1:ions%natoms))
3144 call partial_charges_calculate(mesh, st, ions, hirshfeld_charges)
3146 if (mpi_grp_is_root(mpi_world))
then
3150 call td_write_print_header_init(out_partial_charges)
3153 call write_iter_header_start(out_partial_charges)
3155 do idir = 1, ions%natoms
3156 write(aux,
'(a13,i3,a1)')
'hirshfeld(atom=', idir,
')'
3157 call write_iter_header(out_partial_charges, aux)
3160 call write_iter_nl(out_partial_charges)
3162 call td_write_print_header_end(out_partial_charges)
3165 call write_iter_start(out_partial_charges)
3167 call write_iter_double(out_partial_charges, hirshfeld_charges, ions%natoms)
3169 call write_iter_nl(out_partial_charges)
3172 safe_deallocate_a(hirshfeld_charges)
3178 subroutine td_write_q(out_q, space, ks, iter)
3179 type(c_ptr),
intent(inout) :: out_q
3180 class(space_t),
intent(in) :: space
3181 type(v_ks_t),
intent(in) :: ks
3182 integer,
intent(in) :: iter
3185 character(len=50) :: aux
3189 if (mpi_grp_is_root(mpi_world))
then
3191 call td_write_print_header_init(out_q)
3192 call write_iter_header_start(out_q)
3193 do ii = 1, ks%pt%nmodes
3194 write(aux,
'(a1,i3,a3)')
'q', ii,
'(t)'
3195 call write_iter_header(out_q, aux)
3197 do ii = 1, ks%pt%nmodes
3198 write(aux,
'(a1,i3,a3)')
'p', ii,
'(t)'
3199 call write_iter_header(out_q, aux)
3201 do ii = 1, space%dim
3202 write(aux,
'(a3,i3,a3)')
'f_pt', ii,
'(t)'
3203 call write_iter_header(out_q, aux)
3205 call write_iter_nl(out_q)
3206 call td_write_print_header_end(out_q)
3209 call write_iter_start(out_q)
3210 call write_iter_double(out_q, ks%pt_mx%pt_q, ks%pt%nmodes)
3211 call write_iter_double(out_q, ks%pt_mx%pt_p, ks%pt%nmodes)
3212 call write_iter_double(out_q, ks%pt_mx%fmf, space%dim)
3213 call write_iter_nl(out_q)
3222 type(c_ptr),
intent(inout) :: out_mxll
3223 class(space_t),
intent(in) :: space
3224 type(hamiltonian_elec_t),
intent(in) :: hm
3225 real(real64),
intent(in) :: dt
3226 integer,
intent(in) :: iter
3229 real(real64) :: field(space%dim)
3230 character(len=80) :: aux
3231 character(len=1) :: field_char
3235 if (.not. mpi_grp_is_root(mpi_world))
then
3241 call td_write_print_header_init(out_mxll)
3243 write(aux,
'(a7,e20.12,3a)')
'# dt = ', units_from_atomic(units_out%time, dt), &
3244 " [", trim(units_abbrev(units_out%time)),
"]"
3245 call write_iter_string(out_mxll, aux)
3246 call write_iter_nl(out_mxll)
3248 call write_iter_header_start(out_mxll)
3249 select case (hm%mxll%coupling_mode)
3250 case (length_gauge_dipole, multipolar_expansion)
3251 if (hm%mxll%add_electric_dip) field_char =
'E'
3252 if (hm%mxll%add_magnetic_dip) field_char =
'B'
3253 do idir = 1, space%dim
3254 write(aux,
'(a,i1,a)') field_char //
'(', idir,
')'
3255 call write_iter_header(out_mxll, aux)
3257 case (velocity_gauge_dipole)
3258 do idir = 1, space%dim
3259 write(aux,
'(a,i1,a)')
'A(', idir,
')'
3260 call write_iter_header(out_mxll, aux)
3263 call write_iter_nl(out_mxll)
3265 call write_iter_string(out_mxll,
'#[Iter n.]')
3266 call write_iter_header(out_mxll,
'[' // trim(units_abbrev(units_out%time)) //
']')
3270 select case (hm%mxll%coupling_mode)
3271 case (length_gauge_dipole, multipolar_expansion)
3272 if (hm%mxll%add_electric_dip) aux =
'[' // trim(units_abbrev(units_out%force)) //
']'
3273 if (hm%mxll%add_magnetic_dip) aux =
'[' // trim(units_abbrev(unit_one/units_out%length**2)) //
']'
3274 do idir = 1, space%dim
3275 call write_iter_header(out_mxll, aux)
3277 case (velocity_gauge_dipole)
3278 aux =
'[' // trim(units_abbrev(units_out%energy)) //
']'
3279 do idir = 1, space%dim
3280 call write_iter_header(out_mxll, aux)
3283 call write_iter_nl(out_mxll)
3284 call td_write_print_header_end(out_mxll)
3287 call write_iter_start(out_mxll)
3290 select case (hm%mxll%coupling_mode)
3291 case (length_gauge_dipole, multipolar_expansion)
3292 if (hm%mxll%add_electric_dip) field = units_from_atomic(units_out%force, hm%mxll%e_field_dip)
3293 if (hm%mxll%add_magnetic_dip) field = units_from_atomic(unit_one/units_out%length**2, hm%mxll%b_field_dip)
3294 call write_iter_double(out_mxll, field, space%dim)
3295 case (velocity_gauge_dipole)
3296 field = units_from_atomic(units_out%energy, hm%mxll%vec_pot_dip)
3297 call write_iter_double(out_mxll, field, space%dim)
3299 call write_iter_nl(out_mxll)
3307 type(c_ptr),
intent(inout) :: out_coords
3308 type(lda_u_t),
intent(in) :: lda_u
3309 integer,
intent(in) :: iter
3312 character(len=50) :: aux
3314 if (.not. mpi_grp_is_root(mpi_world))
return
3319 call td_write_print_header_init(out_coords)
3322 call write_iter_header_start(out_coords)
3324 do ios = 1, lda_u%norbsets
3325 write(aux,
'(a2,i3,a1)')
'Ueff(', ios,
')'
3326 call write_iter_header(out_coords, aux)
3329 do ios = 1, lda_u%norbsets
3330 write(aux,
'(a2,i3,a1)')
'U(', ios,
')'
3331 call write_iter_header(out_coords, aux)
3334 do ios = 1, lda_u%norbsets
3335 write(aux,
'(a2,i3,a1)')
'J(', ios,
')'
3336 call write_iter_header(out_coords, aux)
3339 if (lda_u%intersite)
then
3340 do ios = 1, lda_u%norbsets
3341 do inn = 1, lda_u%orbsets(ios)%nneighbors
3342 write(aux,
'(a2,i3,a1,i3,a1)')
'V(', ios,
'-', inn,
')'
3343 call write_iter_header(out_coords, aux)
3349 call write_iter_nl(out_coords)
3352 call write_iter_string(out_coords,
'#[Iter n.]')
3353 call write_iter_header(out_coords,
'[' // trim(units_abbrev(units_out%time)) //
']')
3354 call write_iter_string(out_coords, &
3355 'Effective U ' // trim(units_abbrev(units_out%energy)) // &
3356 ', U in '// trim(units_abbrev(units_out%energy)) // &
3357 ', J in ' // trim(units_abbrev(units_out%energy)))
3358 call write_iter_nl(out_coords)
3360 call td_write_print_header_end(out_coords)
3363 call write_iter_start(out_coords)
3365 do ios = 1, lda_u%norbsets
3366 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3367 lda_u%orbsets(ios)%Ueff), 1)
3370 do ios = 1, lda_u%norbsets
3371 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3372 lda_u%orbsets(ios)%Ubar), 1)
3375 do ios = 1, lda_u%norbsets
3376 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3377 lda_u%orbsets(ios)%Jbar), 1)
3380 if (lda_u%intersite)
then
3381 do ios = 1, lda_u%norbsets
3382 do inn = 1, lda_u%orbsets(ios)%nneighbors
3383 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3384 lda_u%orbsets(ios)%V_ij(inn,0)), 1)
3389 call write_iter_nl(out_coords)
3396 type(c_ptr),
intent(inout) :: file_handle
3397 type(grid_t),
intent(in) :: grid
3398 type(kpoints_t),
intent(in) :: kpoints
3399 type(states_elec_t),
intent(in) :: st
3400 integer,
intent(in) :: iter
3402 integer :: ik_ispin, ist
3403 character(len=7) :: nkpt_str, nst_str
3404 character(len=7) :: ik_str, ist_str
3405 real(real64),
allocatable :: norm_ks(:, :)
3406 real(real64) :: n_electrons
3410 safe_allocate(norm_ks(1:st%nst, 1:st%nik))
3411 call states_elec_calc_norms(grid, kpoints, st, norm_ks)
3413 if (mpi_grp_is_root(mpi_world))
then
3416 call td_write_print_header_init(file_handle)
3419 write(nkpt_str,
'(I7)') st%nik
3420 write(nst_str,
'(I7)') st%nst
3421 call write_iter_string(file_handle,
'# Dimensions. (nstates, nkpt * nspin):')
3422 call write_iter_string(file_handle, trim(adjustl(nst_str)) //
' ' // trim(adjustl(nkpt_str)))
3423 call write_iter_nl(file_handle)
3426 call write_iter_string(file_handle,
'# Norm ordering: (istate, ikpoint_spin)')
3427 call write_iter_nl(file_handle)
3430 call write_iter_header_start(file_handle)
3431 call write_iter_header(file_handle,
'N_electrons')
3432 do ik_ispin = 1, st%nik
3434 write(ik_str,
'(I7)') ik_ispin
3435 write(ist_str,
'(I7)') ist
3436 call write_iter_header(file_handle, &
3437 'Norm (' // trim(ist_str) //
',' // trim(ik_str) //
')')
3440 call write_iter_nl(file_handle)
3441 call td_write_print_header_end(file_handle)
3444 n_electrons = sum(st%occ * norm_ks**2)
3447 call write_iter_start(file_handle)
3448 call write_iter_double(file_handle, n_electrons, 1)
3449 do ik_ispin = 1, st%nik
3450 call write_iter_double(file_handle, norm_ks(:, ik_ispin),
size(norm_ks, 1))
3452 call write_iter_nl(file_handle)
3456 safe_deallocate_a(norm_ks)
3465 type(c_ptr),
intent(inout) :: file_handle
3466 type(ions_t),
intent(in) :: ions
3467 integer,
intent(in) :: iter
3470 real(real64) :: tmp(3)
3472 if (.not. mpi_grp_is_root(mpi_world))
return
3476 assert(ions%space%dim == 3)
3479 call td_write_print_header_init(file_handle)
3482 call write_iter_header_start(file_handle)
3484 call write_iter_string(file_handle,
'# Iter, a, b, c, volume, alpha, beta, gamma, ' &
3485 //
'a_x, a_y, a_z, b_x, b_y, b_z, c_x, c_y, c_z')
3488 call write_iter_string(file_handle,
'#[Iter n.]')
3489 call write_iter_header(file_handle,
'[' // trim(units_abbrev(units_out%time)) //
']')
3490 call write_iter_string(file_handle, &
3491 'Lengths in ' // trim(units_abbrev(units_out%length)) // &
3492 ', Volume in ' // trim(units_abbrev(units_out%length**3)) // &
3493 ', Angles in degree, Lattice vectors in '// trim(units_abbrev(units_out%length)))
3494 call write_iter_nl(file_handle)
3496 call td_write_print_header_end(file_handle)
3499 call write_iter_start(file_handle)
3503 tmp(idir) = units_from_atomic(units_out%length, norm2(ions%latt%rlattice(1:3, idir)))
3505 call write_iter_double(file_handle, tmp, 3)
3508 tmp(1) = units_from_atomic(units_out%length**3, ions%latt%rcell_volume)
3509 call write_iter_double(file_handle, tmp(1), 1)
3512 call write_iter_double(file_handle, ions%latt%alpha, 1)
3513 call write_iter_double(file_handle, ions%latt%beta, 1)
3514 call write_iter_double(file_handle, ions%latt%gamma, 1)
3518 tmp(1:3) = units_from_atomic(units_out%length, ions%latt%rlattice(:, idir))
3519 call write_iter_double(file_handle, tmp, 3)
3521 call write_iter_nl(file_handle)
3530 type(namespace_t),
intent(in) :: namespace
3531 integer,
intent(in) :: iter
3532 real(real64),
intent(in) :: dt
3534 integer :: default, flags, iout, first
3589 call parse_variable(namespace,
'MaxwellTDOutput', default, flags)
3591 if (.not. varinfo_valid_option(
'MaxwellTDOutput', flags, is_flag = .
true.))
then
3592 call messages_input_error(namespace,
'MaxwellTDOutput')
3596 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3597 if (writ%out(iout)%write)
then
3598 writ%out(iout + 1)%write = .
true.
3599 writ%out(iout + 2)%write = .
true.
3604 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3613 call io_mkdir(
'td.general', namespace)
3618 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_x", namespace)))
3620 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_y", namespace)))
3622 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_e_field_z", namespace)))
3628 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_x", namespace)))
3630 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_y", namespace)))
3632 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/total_b_field_z", namespace)))
3638 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_x", namespace)))
3640 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_y", namespace)))
3642 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_e_field_z", namespace)))
3648 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_x", namespace)))
3650 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_y", namespace)))
3652 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/longitudinal_b_field_z", namespace)))
3658 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_x", namespace)))
3660 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_y", namespace)))
3662 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_e_field_z", namespace)))
3668 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_x", namespace)))
3670 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_y", namespace)))
3672 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/transverse_b_field_z", namespace)))
3678 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/maxwell_energy", namespace)))
3683 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-x", namespace)))
3688 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-y", namespace)))
3693 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/electric_field_surface-z", namespace)))
3698 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-x", namespace)))
3703 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-y", namespace)))
3708 units_from_atomic(units_out%time, dt), trim(io_workpath(
"td.general/magnetic_field_surface-z", namespace)))
3724 if (writ%out(iout)%write)
call write_iter_end(writ%out(iout)%handle)
3732 subroutine td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
3734 class(space_t),
intent(in) :: space
3735 type(grid_t),
intent(inout) :: gr
3736 type(states_mxll_t),
intent(inout) :: st
3737 type(hamiltonian_mxll_t),
intent(inout) :: hm
3738 type(helmholtz_decomposition_t),
intent(inout) :: helmholtz
3739 real(real64),
intent(in) :: dt
3740 integer,
intent(in) :: iter
3741 type(namespace_t),
intent(in) :: namespace
3746 call profiling_in(
"TD_WRITE_ITER_MAXWELL")
3749 call helmholtz%get_trans_field(namespace, st%rs_state_trans, total_field=st%rs_state)
3750 call get_rs_state_at_point(st%selected_points_rs_state_trans(:,:), st%rs_state_trans, &
3751 st%selected_points_coordinate(:,:), st, gr)
3754 hm%energy%energy_trans = m_zero
3758 call helmholtz%get_long_field(namespace, st%rs_state_long, total_field=st%rs_state)
3759 call get_rs_state_at_point(st%selected_points_rs_state_long(:,:), st%rs_state_long, &
3760 st%selected_points_coordinate(:,:), st, gr)
3763 hm%energy%energy_long = m_zero
3855 call profiling_out(
"TD_WRITE_ITER_MAXWELL")
3864 type(hamiltonian_mxll_t),
intent(in) :: hm
3865 integer,
intent(in) :: iter
3869 integer :: n_columns
3871 if (.not. mpi_grp_is_root(mpi_world))
return
3896 call write_iter_header(
out_maxwell_energy,
'[' // trim(units_abbrev(units_out%time)) //
']')
3898 do ii = 1, n_columns
3899 call write_iter_header(
out_maxwell_energy,
'[' // trim(units_abbrev(units_out%energy)) //
']')
3907 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy), 1)
3908 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%e_energy), 1)
3909 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%b_energy), 1)
3911 hm%energy%energy+hm%energy%boundaries), 1)
3912 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%boundaries), 1)
3913 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_trans), 1)
3914 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_long), 1)
3915 call write_iter_double(
out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_plane_waves), 1)
3924 type(c_ptr),
intent(inout) :: out_field_surf
3925 type(states_mxll_t),
intent(in) :: st
3926 integer,
intent(in) :: dim
3927 integer,
intent(in) :: iter
3931 integer :: n_columns
3933 if (.not. mpi_grp_is_root(mpi_world))
return
3940 call td_write_print_header_init(out_field_surf)
3943 call write_iter_header_start(out_field_surf)
3944 call write_iter_header(out_field_surf,
'- x direction')
3945 call write_iter_header(out_field_surf,
'+ x direction')
3946 call write_iter_header(out_field_surf,
'- y direction')
3947 call write_iter_header(out_field_surf,
'+ y direction')
3948 call write_iter_header(out_field_surf,
'- z direction')
3949 call write_iter_header(out_field_surf,
'+ z direction')
3950 call write_iter_header(out_field_surf,
'- x dir. p. w.')
3951 call write_iter_header(out_field_surf,
'+ x dir. p. w.')
3952 call write_iter_header(out_field_surf,
'- y dir. p. w.')
3953 call write_iter_header(out_field_surf,
'+ y dir. p. w.')
3954 call write_iter_header(out_field_surf,
'- z dir. p. w.')
3955 call write_iter_header(out_field_surf,
'+ z dir. p. w.')
3957 call write_iter_nl(out_field_surf)
3960 call write_iter_string(out_field_surf,
'#[Iter n.]')
3961 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%time)) //
']')
3963 do ii = 1, n_columns
3964 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%energy/units_out%length)) //
']')
3966 call write_iter_nl(out_field_surf)
3968 call td_write_print_header_end(out_field_surf)
3971 call write_iter_start(out_field_surf)
3972 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3973 st%electric_field_box_surface(1,1,dim)), 1)
3974 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3975 st%electric_field_box_surface(2,1,dim)), 1)
3976 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3977 st%electric_field_box_surface(1,2,dim)), 1)
3978 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3979 st%electric_field_box_surface(2,2,dim)), 1)
3980 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3981 st%electric_field_box_surface(1,3,dim)), 1)
3982 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3983 st%electric_field_box_surface(2,3,dim)), 1)
3984 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3985 st%electric_field_box_surface_plane_waves(1,1,dim)), 1)
3986 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3987 st%electric_field_box_surface_plane_waves(2,1,dim)), 1)
3988 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3989 st%electric_field_box_surface_plane_waves(1,2,dim)), 1)
3990 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3991 st%electric_field_box_surface_plane_waves(2,2,dim)), 1)
3992 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3993 st%electric_field_box_surface_plane_waves(1,3,dim)), 1)
3994 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3995 st%electric_field_box_surface_plane_waves(2,3,dim)), 1)
3996 call write_iter_nl(out_field_surf)
4004 type(c_ptr),
intent(inout) :: out_field_surf
4005 type(states_mxll_t),
intent(in) :: st
4006 integer,
intent(in) :: dim
4007 integer,
intent(in) :: iter
4011 integer :: n_columns
4013 if (.not. mpi_grp_is_root(mpi_world))
return
4020 call td_write_print_header_init(out_field_surf)
4023 call write_iter_header_start(out_field_surf)
4024 call write_iter_header(out_field_surf,
'- x direction')
4025 call write_iter_header(out_field_surf,
'+ x direction')
4026 call write_iter_header(out_field_surf,
'- y direction')
4027 call write_iter_header(out_field_surf,
'+ y direction')
4028 call write_iter_header(out_field_surf,
'- z direction')
4029 call write_iter_header(out_field_surf,
'+ z direction')
4030 call write_iter_header(out_field_surf,
'- x dir. p. w.')
4031 call write_iter_header(out_field_surf,
'+ x dir. p. w.')
4032 call write_iter_header(out_field_surf,
'- y dir. p. w.')
4033 call write_iter_header(out_field_surf,
'+ y dir. p. w.')
4034 call write_iter_header(out_field_surf,
'- z dir. p. w.')
4035 call write_iter_header(out_field_surf,
'+ z dir. p. w.')
4037 call write_iter_nl(out_field_surf)
4040 call write_iter_string(out_field_surf,
'#[Iter n.]')
4041 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(units_out%time)) //
']')
4043 do ii = 1, n_columns
4044 call write_iter_header(out_field_surf,
'[' // trim(units_abbrev(unit_one/units_out%length**2)) //
']')
4046 call write_iter_nl(out_field_surf)
4048 call td_write_print_header_end(out_field_surf)
4051 call write_iter_start(out_field_surf)
4052 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4053 st%magnetic_field_box_surface(1,1,dim)), 1)
4054 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4055 st%magnetic_field_box_surface(2,1,dim)), 1)
4056 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4057 st%magnetic_field_box_surface(1,2,dim)), 1)
4058 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4059 st%magnetic_field_box_surface(2,2,dim)), 1)
4060 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4061 st%magnetic_field_box_surface(1,3,dim)), 1)
4062 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4063 st%magnetic_field_box_surface(2,3,dim)), 1)
4064 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4065 st%magnetic_field_box_surface_plane_waves(1,1,dim)), 1)
4066 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4067 st%magnetic_field_box_surface_plane_waves(2,1,dim)), 1)
4068 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4069 st%magnetic_field_box_surface_plane_waves(1,2,dim)), 1)
4070 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4071 st%magnetic_field_box_surface_plane_waves(2,2,dim)), 1)
4072 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4073 st%magnetic_field_box_surface_plane_waves(1,3,dim)), 1)
4074 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4075 st%magnetic_field_box_surface_plane_waves(2,3,dim)), 1)
4076 call write_iter_nl(out_field_surf)
4082 subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
4083 type(c_ptr),
intent(inout) :: out_fields
4084 class(space_t),
intent(in) :: space
4085 type(states_mxll_t),
intent(in) :: st
4086 integer,
intent(in) :: iter
4087 real(real64),
intent(in) :: dt
4088 integer,
intent(in) :: e_or_b_field
4089 integer,
intent(in) :: field_type
4090 integer,
intent(in) :: idir
4094 real(real64) :: field(space%dim), selected_field
4095 character(len=80) :: aux
4097 if (.not. mpi_grp_is_root(mpi_world))
return
4102 call td_write_print_header_init(out_fields)
4105 write(aux,
'(a7,e20.12,3a)')
'# dt = ', units_from_atomic(units_out%time, dt), &
4106 " [", trim(units_abbrev(units_out%time)),
"]"
4107 call write_iter_string(out_fields, aux)
4108 call write_iter_nl(out_fields)
4110 call write_iter_header_start(out_fields)
4112 do id = 1, st%selected_points_number
4113 select case (e_or_b_field)
4115 write(aux,
'(a,i1,a)')
'E(', id,
')'
4117 write(aux,
'(a,i1,a)')
'B(', id,
')'
4119 call write_iter_header(out_fields, aux)
4122 call write_iter_nl(out_fields)
4123 call write_iter_string(out_fields,
'#[Iter n.]')
4124 call write_iter_header(out_fields,
' [' // trim(units_abbrev(units_out%time)) //
']')
4129 aux =
' [' // trim(units_abbrev(units_out%force)) //
']'
4130 do id = 1, st%selected_points_number
4131 call write_iter_header(out_fields, aux)
4133 call write_iter_nl(out_fields)
4134 call td_write_print_header_end(out_fields)
4137 call write_iter_start(out_fields)
4139 do id = 1, st%selected_points_number
4140 select case (e_or_b_field)
4143 select case (field_type)
4145 call get_electric_field_vector(st%selected_points_rs_state(:,id), field(1:st%dim))
4147 call get_electric_field_vector(st%selected_points_rs_state_long(:,id), field(1:st%dim))
4149 call get_electric_field_vector(st%selected_points_rs_state_trans(:,id), field(1:st%dim))
4151 selected_field = units_from_atomic(units_out%energy/units_out%length, field(idir))
4154 select case (field_type)
4156 call get_magnetic_field_vector(st%selected_points_rs_state(:,id), st%rs_sign, field(1:st%dim))
4158 call get_magnetic_field_vector(st%selected_points_rs_state_long(:,id), st%rs_sign, field(1:st%dim))
4160 call get_magnetic_field_vector(st%selected_points_rs_state_trans(:,id), st%rs_sign, field(1:st%dim))
4162 selected_field = units_from_atomic(unit_one/units_out%length**2, field(idir))
4164 call write_iter_double(out_fields, selected_field, 1)
4167 call write_iter_nl(out_fields)
4176 type(namespace_t),
intent(in) :: namespace
4177 class(space_t),
intent(in) :: space
4178 type(grid_t),
intent(inout) :: gr
4179 type(states_mxll_t),
intent(inout) :: st
4180 type(hamiltonian_mxll_t),
intent(inout) :: hm
4181 type(helmholtz_decomposition_t),
intent(inout) :: helmholtz
4182 type(output_t),
intent(in) :: outp
4183 integer,
intent(in) :: iter
4184 real(real64),
intent(in) :: time
4186 character(len=256) :: filename
4188 push_sub(td_write_maxwell_free_data)
4189 call profiling_in(
"TD_WRITE_MAXWELL_DATA")
4192 write(filename,
'(a,a,i7.7)') trim(outp%iter_dir),
"td.", iter
4194 call output_mxll(outp, namespace, space, gr, st, hm, helmholtz, time, filename)
4196 call profiling_out(
"TD_WRITE_MAXWELL_DATA")
4197 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.
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.
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
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)
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_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.
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)