Octopus
td_write.F90
Go to the documentation of this file.
1! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module td_write_oct_m
22 use blas_oct_m
23 use comm_oct_m
25 use debug_oct_m
33 use global_oct_m
34 use grid_oct_m
35 use output_oct_m
44 use io_oct_m
46 use ions_oct_m
47 use kick_oct_m
48 use, intrinsic :: iso_fortran_env
50 use lasers_oct_m
53 use lda_u_oct_m
56 use math_oct_m
59 use mesh_oct_m
62 use mpi_oct_m
67 use parser_oct_m
73 use space_oct_m
83 use types_oct_m
84 use unit_oct_m
86 use utils_oct_m
88 use v_ks_oct_m
90
91 implicit none
92
93 private
94 public :: &
95 td_write_t, &
106
108 integer, parameter, public :: &
109 OUT_MULTIPOLES = 1, &
110 out_angular = 2, &
111 out_spin = 3, &
112 out_populations = 4, &
113 out_coords = 5, &
114 out_acc = 6, &
115 out_laser = 7, &
116 out_energy = 8, &
117 out_proj = 9, &
118 out_magnets = 10, &
119 out_gauge_field = 11, &
120 out_temperature = 12, &
121 out_ftchd = 13, &
122 out_vel = 14, &
123 out_eigs = 15, &
124 out_ion_ch = 16, &
125 out_total_current = 17, &
126 out_partial_charges = 18, &
127 out_kp_proj = 19, &
128 out_floquet = 20, &
129 out_n_ex = 21, &
130 out_separate_coords = 22, &
132 out_separate_forces = 24, &
134 out_tot_m = 26, &
135 out_q = 27, &
136 out_mxll_field = 28, &
137 out_norm_ks = 29, &
138 ! This constant should always equal the last integer above
139 ! as it defines the total numer of output options
140 out_max = 29
141
142
144 ! Entries must be consistent with the ordering above
145 character(len=100) :: td_file_name(OUT_MAX) = [character(100) :: &
146 "NULL", & ! Do not include OUT_MULTIPOLES (extension is a function of state index)
147 "angular", &
148 "spin", &
149 "populations", &
150 "coordinates", &
151 "acceleration", &
152 "NULL", & ! Do not include OUT_LASER
153 "energy", &
154 "projections", &
155 "magnetic_moments", &
156 "gauge_field", &
157 "temperature", &
158 "NULL", & ! Do not include OUT_FTCHD
159 "velocity", &
160 "eigenvalues", &
161 "ion_ch", &
162 "total_current", &
163 "partial_charges", &
164 "projections", &
165 "floquet_bands", &
166 "n_ex", &
167 "onlyCoordinates", &
168 "onlyVelocities", &
169 "onlyForces", &
170 "total_heat_current", &
171 "total_magnetization", &
172 "photons_q", &
173 "maxwell_dipole_field", &
174 "norm_wavefunctions"]
175
176 integer, parameter :: &
177 OUT_DFTU_EFFECTIVE_U = 1, &
178 out_dftu_max = 1
179
180
181 integer, parameter :: &
182 OUT_MAXWELL_TOTAL_E_FIELD = 1, &
188 out_maxwell_energy = 19, &
195
196 integer, parameter, public :: &
197 OUT_MAXWELL_MAX = 25
198
199 integer, parameter :: &
200 MAXWELL_TOTAL_FIELD = 1, &
203
204 integer, parameter :: &
205 maxwell_e_field = 1, &
207
209 type td_write_prop_t
210 type(c_ptr) :: handle
211 type(c_ptr), allocatable :: mult_handles(:)
212 type(mpi_grp_t) :: mpi_grp
213 integer :: hand_start
214 integer :: hand_end
215 logical :: write = .false.
216 logical :: resolve_states = .false.
217 end type td_write_prop_t
218
222 type td_write_t
223 private
224 type(td_write_prop_t), public :: out(out_max)
225 type(td_write_prop_t) :: out_dftu(out_dftu_max)
226
227 integer :: lmax
228 real(real64) :: lmm_r
229 type(states_elec_t) :: gs_st
230 ! calculate the projections(s) onto it.
231 integer :: n_excited_states
232 type(excited_states_t), allocatable :: excited_st(:)
233 integer :: compute_interval
234 end type td_write_t
235
236contains
237
239 subroutine td_write_kick(outp, namespace, space, mesh, kick, ions, iter)
240 type(output_t), intent(in) :: outp
241 type(namespace_t), intent(in) :: namespace
242 type(electron_space_t), intent(in) :: space
243 class(mesh_t), intent(in) :: mesh
244 type(kick_t), intent(in) :: kick
245 type(ions_t), intent(in) :: ions
246 integer, intent(in) :: iter
247
248 complex(real64), allocatable :: kick_function(:)
249 character(len=256) :: filename
250 integer :: err
251
252 push_sub(td_write_kick)
253
254 write(filename, '(a,i7.7)') "td.", iter ! name of directory
255 if (outp%what(option__output__delta_perturbation)) then
256 safe_allocate(kick_function(1:mesh%np))
257 call kick_function_get(space, mesh, kick, kick_function, 1)
258 call zio_function_output(outp%how(option__output__delta_perturbation), filename, "kick_function", namespace, &
259 space, mesh, kick_function(:), units_out%energy, err, pos=ions%pos, atoms=ions%atom)
260 safe_deallocate_a(kick_function)
261 end if
262
263 pop_sub(td_write_kick)
264 end subroutine td_write_kick
265
266
276 subroutine td_write_init(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, ks, ions_move, &
277 with_gauge_field, kick, iter, max_iter, dt, mc)
278 type(td_write_t), target, intent(out) :: writ
279 type(namespace_t), intent(in) :: namespace
280 type(electron_space_t), intent(in) :: space
281 type(output_t), intent(inout) :: outp
282 type(grid_t), intent(in) :: gr
283 type(states_elec_t), intent(inout) :: st
284 type(hamiltonian_elec_t), intent(inout) :: hm
285 type(ions_t), intent(in) :: ions
286 type(partner_list_t), intent(in) :: ext_partners
287 type(v_ks_t), intent(inout) :: ks
288 logical, intent(in) :: ions_move
289 logical, intent(in) :: with_gauge_field
290 type(kick_t), intent(in) :: kick
291 integer, intent(in) :: iter
292 integer, intent(in) :: max_iter
293 real(real64), intent(in) :: dt
294 type(multicomm_t), intent(in) :: mc
295
296 real(real64) :: rmin
297 integer :: ierr, first, ii, ist, jj, flags, iout, default, ifile
298 logical :: output_options(max_output_types)
299 integer :: output_interval(0:max_output_types)
300 integer(int64) :: how(0:max_output_types)
301 type(block_t) :: blk
302 character(len=MAX_PATH_LEN) :: filename
304 logical :: resolve_states
305 logical, allocatable :: skip(:)
309 !%Variable TDOutput
310 !%Type block
311 !%Default multipoles + energy (+ others depending on other options)
312 !%Section Time-Dependent::TD Output
313 !%Description
314 !% Defines what should be output during the time-dependent
315 !% simulation. Many of the options can increase the computational
316 !% cost of the simulation, so only use the ones that you need. In
317 !% most cases the default value is enough, as it is adapted to the
318 !% details of the TD run. If the ions are allowed to be moved, additionally
319 !% the geometry and the temperature are output. If a laser is
320 !% included it will output by default.
321 !%
322 !% Note: the output files generated by this option are updated
323 !% every <tt>RestartWriteInterval</tt> steps.
324 !%
325 !% Example:
326 !% <br><br><tt>%TDOutput
327 !% <br>&nbsp;&nbsp;multipoles
328 !% <br>&nbsp;&nbsp;energy
329 !% <br>%<br></tt>
330 !%
331 !%Option multipoles 1
332 !% Outputs the (electric) multipole moments of the density to the file <tt>td.general/multipoles</tt>.
333 !% This is required to, <i>e.g.</i>, calculate optical absorption spectra of finite systems. The
334 !% maximum value of <math>l</math> can be set with the variable <tt>TDMultipoleLmax</tt>.
335 !%Option angular 2
336 !% Outputs the orbital angular momentum of the system to <tt>td.general/angular</tt>.
337 !% This can be used to calculate circular dicrhoism. Not implemented for periodic systems.
338 !% In case of a nonlocal pseudopotential, the gauge used can be specified using
339 !% the variable <tt>MagneticGaugeCorrection</tt>.
340 !%Option spin 3
341 !% (Experimental) Outputs the expectation value of the spin, which can be used to calculate magnetic
342 !% circular dichroism.
343 !%Option populations 4
344 !% (Experimental) Outputs the projection of the time-dependent
345 !% Kohn-Sham Slater determinant onto the ground state (or
346 !% approximations to the excited states) to the file
347 !% <tt>td.general/populations</tt>. Note that the calculation of
348 !% populations is expensive in memory and computer time, so it
349 !% should only be used if it is really needed. See <tt>TDExcitedStatesToProject</tt>.
350 !%Option geometry 5
351 !% If set (and if the atoms are allowed to move), outputs the coordinates, velocities,
352 !% and forces of the atoms to the the file <tt>td.general/coordinates</tt>. On by default if <tt>MoveIons = yes</tt>.
353 !%Option dipole_acceleration 6
354 !% When set, outputs the acceleration of the electronic dipole, calculated from the Ehrenfest theorem,
355 !% in the file <tt>td.general/acceleration</tt>. This file can then be
356 !% processed by the utility <tt>oct-harmonic-spectrum</tt> in order to obtain the harmonic spectrum.
357 !%Option laser 7
358 !% If set, outputs the laser field to the file <tt>td.general/laser</tt>.
359 !% On by default if <tt>TDExternalFields</tt> is set.
360 !%Option energy 8
361 !% If set, <tt>octopus</tt> outputs the different components of the energy
362 !% to the file <tt>td.general/energy</tt>. Will be zero except for every <tt>TDEnergyUpdateIter</tt> iterations.
363 !%Option td_occup 9
364 !% (Experimental) If set, outputs the projections of the
365 !% time-dependent Kohn-Sham wavefunctions onto the static
366 !% (zero-time) wavefunctions to the file
367 !% <tt>td.general/projections.XXX</tt>. Only use this option if
368 !% you really need it, as it might be computationally expensive. See <tt>TDProjStateStart</tt>.
369 !% The output interval of this quantity is controled by the variable <tt>TDOutputComputeInterval</tt>
370 !% In case of states parallelization, all the ground-state states are stored by each task.
371 !%Option local_mag_moments 10
372 !% If set, outputs the local magnetic moments, integrated in sphere centered around each atom.
373 !% The radius of the sphere can be set with <tt>LocalMagneticMomentsSphereRadius</tt>.
374 !%Option gauge_field 11
375 !% If set, outputs the vector gauge field corresponding to a spatially uniform (but time-dependent)
376 !% external electrical potential. This is only useful in a time-dependent periodic run.
377 !% On by default if <tt>GaugeVectorField</tt> is set.
378 !%Option temperature 12
379 !% If set, the ionic temperature at each step is printed. On by default if <tt>MoveIons = yes</tt>.
380 !%Option ftchd 13
381 !% Write Fourier transform of the electron density to the file <tt>ftchd.X</tt>,
382 !% where X depends on the kick (e.g. with sin-shaped perturbation X=sin).
383 !% This is needed for calculating the dynamic structure factor.
384 !% In the case that the kick mode is qbessel, the written quantity is integral over
385 !% density, multiplied by spherical Bessel function times real spherical harmonic.
386 !% On by default if <tt>TDMomentumTransfer</tt> is set.
387 !%Option dipole_velocity 14
388 !% When set, outputs the dipole velocity, calculated from the Ehrenfest theorem,
389 !% in the file <tt>td.general/velocity</tt>. This file can then be
390 !% processed by the utility <tt>oct-harmonic-spectrum</tt> in order to obtain the harmonic spectrum.
391 !%Option eigenvalues 15
392 !% Write the KS eigenvalues.
393 !%Option ionization_channels 16
394 !% Write the multiple-ionization channels using the KS orbital densities as proposed in
395 !% C. Ullrich, Journal of Molecular Structure: THEOCHEM 501, 315 (2000).
396 !%Option total_current 17
397 !% Output the total current (average of the current density over the cell).
398 !%Option partial_charges 18
399 !% Bader and Hirshfeld partial charges. The output file is called 'td.general/partial_charges'.
400 !%Option td_kpoint_occup 19
401 !% Project propagated Kohn-Sham states to the states at t=0 given in the directory
402 !% restart_proj (see %RestartOptions). This is an alternative to the option
403 !% td_occup, with a formating more suitable for k-points and works only in
404 !% k- and/or state parallelization
405 !%Option td_floquet 20
406 !% Compute non-interacting Floquet bandstructure according to further options:
407 !% TDFloquetFrequency, TDFloquetSample, TDFloquetDimension.
408 !% This is done only once per td-run at t=0.
409 !% works only in k- and/or state parallelization
410 !%Option n_excited_el 21
411 !% Output the number of excited electrons, based on the projections
412 !% of the time evolved wave-functions on the ground-state wave-functions.
413 !% The output interval of this quantity is controled by the variable <tt>TDOutputComputeInterval</tt>
414 !% Note that if one sets RecalculateGSDuringEvolution=yes,
415 !% the code will recompute the GS states
416 !% and use them for the computing the number of excited electrons.
417 !% This is useful if ions are moving or if one wants to get the number of excited electrons in the basis
418 !% of the instantaneous eigenvalues of the Hamiltonian (Houston states).
419 !%Option coordinates_sep 22
420 !% Writes geometries in a separate file.
421 !%Option velocities_sep 23
422 !% Writes velocities in a separate file.
423 !%Option forces_sep 24
424 !% Writes forces in a separate file.
425 !%Option total_heat_current 25
426 !% Output the total heat current (average of the heat current density over the cell).
427 !%Option total_magnetization 26
428 !% Writes the total magnetization, where the total magnetization is calculated at the momentum
429 !% defined by <tt>TDMomentumTransfer</tt>.
430 !% This is used to extract the magnon frequency in case of a magnon kick.
431 !%Option photons_q 27
432 !% Writes photons_q in a separate file.
433 !%Option maxwell_field 28
434 !% Writes total electric field (if coupling is in length_geuge) or vector potential
435 !% (if coupling is in velocity_gauge) coming from the interaction with Maxwell systems
436 !% (only if Maxwell-TDDFT coupling is in dipole).
437 !%Option norm_ks_orbitals 29
438 !% Writes the norm of each Kohn-Sham orbital.
439 !% The data is ordered per row as:
440 !% Iteration time (state 1 kpoint 1) (state2 kpoint1) ... (state-Nstates kpoint1) (state1 kpoint2) ... (state-Nstates kpoint-nkpt)
441 !% noting that the kpoint index will also include the spin index for spin-polarised calculations.
442 !%End
443
444 !defaults
445 output_options = .false.
446 output_options(out_multipoles) = .true.
447 output_options(out_energy) = .true.
448 if (ions_move) then
449 output_options(out_coords) = .true.
450 output_options(out_temperature) = .true.
451 end if
452 if (with_gauge_field) output_options(out_gauge_field) = .true.
453 if (list_has_lasers(ext_partners)) output_options(out_laser) = .true.
454 if (kick%qkick_mode /= qkickmode_none) output_options(out_ftchd) = .true.
455
456 call io_function_read_what_how_when(namespace, space, output_options, how, output_interval, &
457 'TDOutput')
458
459 ! Define which files to write
460 do iout = 1, out_max
461 writ%out(iout)%write = output_options(iout)
462 end do
463
464 ! experimental stuff
465 if (writ%out(out_spin)%write) call messages_experimental('TDOutput = spin', namespace=namespace)
466 if (writ%out(out_populations)%write) call messages_experimental('TDOutput = populations', namespace=namespace)
467 if (writ%out(out_proj)%write) call messages_experimental('TDOutput = td_occup', namespace=namespace)
468 if (writ%out(out_ion_ch)%write) call messages_experimental('TDOutput = ionization_channels', namespace=namespace)
469 if (writ%out(out_partial_charges)%write) call messages_experimental('TDOutput = partial_charges', namespace=namespace)
470 if (writ%out(out_kp_proj)%write) call messages_experimental('TDOutput = td_kpoint_occup', namespace=namespace)
471 if (writ%out(out_floquet)%write) call messages_experimental('TDOutput = td_floquet', namespace=namespace)
472 if (writ%out(out_n_ex)%write) call messages_experimental('TDOutput = n_excited_el', namespace=namespace)
473 if (writ%out(out_q)%write) call messages_experimental('TDOutput = photons_q', namespace=namespace)
474 if (writ%out(out_mxll_field)%write) call messages_experimental('TDOutput = maxwell_field', namespace=namespace)
475
476 if (space%is_periodic() .and. writ%out(out_angular)%write) then
477 call messages_not_implemented("TDOutput angular for periodic systems", namespace=namespace)
478 end if
479
480 !See comment in zstates_elec_mpdotp
481 if (space%is_periodic() .and. writ%out(out_populations)%write) then
482 call messages_not_implemented("TDOutput populations for periodic systems", namespace=namespace)
483 end if
484
485 if (writ%out(out_kp_proj)%write.or.writ%out(out_floquet)%write) then
486 ! make sure this is not domain distributed
487 if (gr%np /= gr%np_global) then
488 message(1) = "TDOutput option td_kpoint_occup and td_floquet do not work with domain parallelization"
489 call messages_fatal(1, namespace=namespace)
490 end if
491 end if
492
493 if ((writ%out(out_separate_forces)%write .or. writ%out(out_coords)%write) .and. space%periodic_dim == 1) then
494 call messages_input_error(namespace, 'TDOutput', &
495 'Forces for systems periodic in 1D are not currently implemented and options that output the forces are not allowed.')
496 end if
497
498 ! NTD: The implementation of the option should be really redone properly
499 if (writ%out(out_kp_proj)%write .and. hm%kpoints%nik_skip == 0) then
500 message(1) = "TDOutput option td_kpoint_occup only work with zero-weight k-points at the moment."
501 call messages_fatal(1, namespace=namespace)
502 end if
503
504 !%Variable TDOutputResolveStates
505 !%Type logical
506 !%Default No
507 !%Section Time-Dependent::TD Output
508 !%Description
509 !% Defines whether the output should be resolved by states.
510 !%
511 !% So far only TDOutput = multipoles is supported.
512 !%
513 !%End
514 call parse_variable(namespace, 'TDOutputResolveStates', .false., resolve_states)
515 if (.not. writ%out(out_multipoles)%write) then
516 write(message(1),'(a)') "TDOutputResolveStates works only for TDOutput = multipoles."
517 call messages_warning(1, namespace=namespace)
518 end if
519
520 !%Variable TDMultipoleLmax
521 !%Type integer
522 !%Default 1
523 !%Section Time-Dependent::TD Output
524 !%Description
525 !% Maximum electric multipole of the density output to the file <tt>td.general/multipoles</tt>
526 !% during a time-dependent simulation. Must be non-negative.
527 !%End
528 call parse_variable(namespace, 'TDMultipoleLmax', 1, writ%lmax)
529 if (writ%lmax < 0) then
530 write(message(1), '(a,i6,a)') "Input: '", writ%lmax, "' is not a valid TDMultipoleLmax."
531 message(2) = '(Must be TDMultipoleLmax >= 0 )'
532 call messages_fatal(2, namespace=namespace)
533 end if
534 call messages_obsolete_variable(namespace, 'TDDipoleLmax', 'TDMultipoleLmax')
535
536 ! Compatibility test
537 if ((writ%out(out_acc)%write) .and. ions_move) then
538 message(1) = 'If harmonic spectrum is to be calculated, atoms should not be allowed to move.'
539 call messages_fatal(1, namespace=namespace)
540 end if
541
542 if ((writ%out(out_q)%write) .and. .not.(ks%has_photons)) then
543 message(1) = 'If q(t) is to be calculated, you need to allow for photon modes.'
544 call messages_fatal(1, namespace=namespace)
545 end if
546
547 if ((writ%out(out_mxll_field)%write) .and. .not. (hm%mxll%coupling_mode == velocity_gauge_dipole &
548 .or. hm%mxll%add_electric_dip)) then
549 message(1) = 'If the dipolar field has to be written, MaxwellCouplingMode has to be'
550 message(2) = '"lenght_gauge_dipole" or "velocity_gauge_dipole" and at least one Maxwell system'
551 message(3) = 'must be present.'
552 call messages_fatal(3, namespace=namespace)
553 end if
554
555 rmin = ions%min_distance()
556
557 ! This variable is documented in scf/scf.F90
558 call parse_variable(namespace, 'LocalMagneticMomentsSphereRadius', min(m_half*rmin, lmm_r_single_atom), writ%lmm_r, &
559 unit=units_inp%length)
560
561 if (writ%out(out_proj)%write .or. writ%out(out_populations)%write &
562 .or.writ%out(out_kp_proj)%write .or. writ%out(out_n_ex)%write) then
563
564 if (st%parallel_in_states .and. writ%out(out_populations)%write) then
565 message(1) = "Option TDOutput = populations is not implemented for parallel in states."
566 call messages_fatal(1, namespace=namespace)
567 end if
568
569 if (writ%out(out_proj)%write .or. writ%out(out_populations)%write) then
570 call states_elec_copy(writ%gs_st, st, exclude_wfns = .true., exclude_eigenval = .true.)
571 else
572 ! we want the same layout of gs_st as st
573 call states_elec_copy(writ%gs_st, st, special=.true.)
574 end if
575
576 ! clean up all the stuff we have to reallocate
577 safe_deallocate_a(writ%gs_st%node)
578
579 call restart_init(restart_gs, namespace, restart_proj, restart_type_load, mc, ierr, mesh=gr)
580
581 if (writ%out(out_proj)%write .or. writ%out(out_populations)%write) then
582 if (ierr == 0) then
583 call states_elec_look(restart_gs, ii, jj, writ%gs_st%nst, ierr)
584 end if
585 writ%gs_st%st_end = writ%gs_st%nst
586 if (ierr /= 0) then
587 message(1) = "Unable to read states information."
588 call messages_fatal(1, namespace=namespace)
589 end if
590
591 writ%gs_st%st_start = 1
592 ! do this only when not calculating populations, since all states are needed then
593 if (.not. writ%out(out_populations)%write) then
594 ! We will store the ground-state Kohn-Sham system for all processors.
595 !%Variable TDProjStateStart
596 !%Type integer
597 !%Default 1
598 !%Section Time-Dependent::TD Output
599 !%Description
600 !% To be used with <tt>TDOutput = td_occup</tt>. Not available if <tt>TDOutput = populations</tt>.
601 !% Only output projections to states above <tt>TDProjStateStart</tt>. Usually one is only interested
602 !% in particle-hole projections around the HOMO, so there is no need to calculate (and store)
603 !% the projections of all TD states onto all static states. This sets a lower limit. The upper limit
604 !% is set by the number of states in the propagation and the number of unoccupied states
605 !% available.
606 !%End
607 call parse_variable(namespace, 'TDProjStateStart', 1, writ%gs_st%st_start)
608
609 if (st%parallel_in_states .and. writ%out(out_proj)%write .and. writ%gs_st%st_start > 1) then
610 call messages_not_implemented("TDOutput = td_occup for parallel in states with TDProjStateStart > 1", &
611 namespace=namespace)
612 end if
613 end if
614
615 writ%gs_st%lnst = writ%gs_st%st_end - writ%gs_st%st_start + 1
616
617 call states_elec_deallocate_wfns(writ%gs_st)
618
619 writ%gs_st%parallel_in_states = .false.
620
621 ! allocate memory
622 safe_allocate(writ%gs_st%occ(1:writ%gs_st%nst, 1:writ%gs_st%nik))
623 safe_allocate(writ%gs_st%eigenval(1:writ%gs_st%nst, 1:writ%gs_st%nik))
624
625 !We want all the task to have all the states
626 !States can be distibuted for the states we propagate.
627 safe_allocate(writ%gs_st%node(1:writ%gs_st%nst))
628 writ%gs_st%node(:) = 0
629
630 writ%gs_st%eigenval = huge(writ%gs_st%eigenval)
631 writ%gs_st%occ = m_zero
632 if (writ%gs_st%d%ispin == spinors) then
633 safe_deallocate_a(writ%gs_st%spin)
634 safe_allocate(writ%gs_st%spin(1:3, 1:writ%gs_st%nst, 1:writ%gs_st%nik))
635 end if
636
637 safe_allocate(skip(1:writ%gs_st%nst))
638 skip = .false.
639 skip(1:writ%gs_st%st_start-1) = .true.
640
641 call states_elec_allocate_wfns(writ%gs_st, gr, type_cmplx, packed=.true., skip=skip)
642
643 safe_deallocate_a(skip)
644 end if
645
646 call states_elec_load(restart_gs, namespace, space, writ%gs_st, gr, hm%kpoints, ierr, label = ': gs for TDOutput')
647
648 if (ierr /= 0 .and. ierr /= (writ%gs_st%st_end-writ%gs_st%st_start+1)*writ%gs_st%nik &
649 *writ%gs_st%d%dim*writ%gs_st%mpi_grp%size) then
650 message(1) = "Unable to read wavefunctions for TDOutput."
651 call messages_fatal(1, namespace=namespace)
652 end if
654 end if
655
656 ! Build the excited states...
657 if (writ%out(out_populations)%write) then
658 !%Variable TDExcitedStatesToProject
659 !%Type block
660 !%Section Time-Dependent::TD Output
661 !%Description
662 !% <b>[WARNING: This is a *very* experimental feature]</b>
663 !% To be used with <tt>TDOutput = populations</tt>.
664 !% The population of the excited states
665 !% (as defined by <Phi_I|Phi(t)> where |Phi(t)> is the many-body time-dependent state at
666 !% time <i>t</i>, and |Phi_I> is the excited state of interest) can be approximated -- it is not clear
667 !% how well -- by substituting for those real many-body states the time-dependent Kohn-Sham
668 !% determinant and some modification of the Kohn-Sham ground-state determinant (<i>e.g.</i>,
669 !% a simple HOMO-LUMO substitution, or the Casida ansatz for excited states in linear-response
670 !% theory. If you set <tt>TDOutput</tt> to contain <tt>populations</tt>, you may ask for these approximated
671 !% populations for a number of excited states, which will be described in the files specified
672 !% in this block: each line should be the name of a file that contains one excited state.
673 !%
674 !% This file structure is the one written by the casida run mode, in the files in the directory <tt>*_excitations</tt>.
675 !% The file describes the "promotions" from occupied
676 !% to unoccupied levels that change the initial Slater determinant
677 !% structure specified in ground_state. These promotions are a set
678 !% of electron-hole pairs. Each line in the file, after an optional header, has four
679 !% columns:
680 !%
681 !% <i>i a <math>\sigma</math> weight</i>
682 !%
683 !% where <i>i</i> should be an occupied state, <i>a</i> an unoccupied one, and <math>\sigma</math>
684 !% the spin of the corresponding orbital. This pair is then associated with a
685 !% creation-annihilation pair <math>a^{\dagger}_{a,\sigma} a_{i,\sigma}</math>, so that the
686 !% excited state will be a linear combination in the form:
687 !%
688 !% <math>\left|{\rm ExcitedState}\right> =
689 !% \sum weight(i,a,\sigma) a^{\dagger}_{a,\sigma} a_{i,\sigma} \left|{\rm GroundState}\right></math>
690 !%
691 !% where <i>weight</i> is the number in the fourth column.
692 !% These weights should be normalized to one; otherwise the routine
693 !% will normalize them, and write a warning.
694 !%End
695 if (parse_block(namespace, 'TDExcitedStatesToProject', blk) == 0) then
696 writ%n_excited_states = parse_block_n(blk)
697 safe_allocate(writ%excited_st(1:writ%n_excited_states))
698 do ist = 1, writ%n_excited_states
699 call parse_block_string(blk, ist-1, 0, filename)
700 call excited_states_init(writ%excited_st(ist), writ%gs_st, trim(filename), namespace)
701 end do
702 else
703 writ%n_excited_states = 0
704 end if
705 end if
706
707 !%Variable TDOutputComputeInterval
708 !%Type integer
709 !%Default 50
710 !%Section Time-Dependent::TD Output
711 !%Description
712 !% The TD output requested are computed
713 !% when the iteration number is a multiple of the <tt>TDOutputComputeInterval</tt> variable.
714 !% Must be >= 0. If it is 0, then no output is written.
715 !% Implemented only for projections and number of excited electrons for the moment.
716 !%End
717 call parse_variable(namespace, 'TDOutputComputeInterval', 50, writ%compute_interval)
718 if (writ%compute_interval < 0) then
719 message(1) = "TDOutputComputeInterval must be >= 0."
720 call messages_fatal(1, namespace=namespace)
721 end if
722
723 ! ----------------------------------------------
724 ! Initialize write file handlers
725 ! ----------------------------------------------
726
727 ! Step
728 if (iter == 0) then
729 first = 0
730 else
731 first = iter + 1
732 end if
733
734 ! Root dir for TD files
735 call io_mkdir('td.general', namespace)
736
737 ! ----------------------------------------------
738 ! Initialize write that are MPI agnostic
739 ! ----------------------------------------------
740
741 ! Default
742 writ%out(:)%mpi_grp = mpi_world
743
744 if (mpi_grp_is_root(mpi_world)) then
745
746 do ifile = 1, out_max
747 ! Exceptions that are handled later
748 if (any([out_multipoles, out_laser, out_ftchd] == ifile)) cycle
749
750 if (writ%out(ifile)%write) then
751 call write_iter_init(writ%out(ifile)%handle, &
752 first, &
753 units_from_atomic(units_out%time, dt), &
754 trim(io_workpath("td.general/"//trim(td_file_name(ifile)), namespace)))
755 end if
756
757 enddo
758
759 ! Exceptions
760
761 ! Multipoles, when MPI is not state-parallel
762 if (writ%out(out_multipoles)%write .and. .not. resolve_states) then
763 call write_iter_init(writ%out(out_multipoles)%handle, &
764 first, units_from_atomic(units_out%time, dt), &
765 trim(io_workpath("td.general/multipoles", namespace)))
766 end if
767
768 if (writ%out(out_ftchd)%write) then
769 select case (kick%qkick_mode)
770 case (qkickmode_sin)
771 write(filename, '(a)') 'td.general/ftchd.sin'
772 case (qkickmode_cos)
773 write(filename, '(a)') 'td.general/ftchd.cos'
774 case (qkickmode_bessel)
775 write(filename, '(a, SP, I0.3, a, I0.3)') 'td.general/ftchd.l', kick%qbessel_l, '_m', kick%qbessel_m
776 case default
777 write(filename, '(a)') 'td.general/ftchd'
778 end select
779 call write_iter_init(writ%out(out_ftchd)%handle, &
780 first, units_from_atomic(units_out%time, dt), trim(io_workpath(filename, namespace)))
781 end if
782
783 if (writ%out(out_laser)%write) then
784 if(associated(list_get_lasers(ext_partners))) then
785 ! The laser file is written for the full propagation in one go, so that
786 ! the user can check that the laser is correct and as intended before letting
787 ! the code run for a possibly large period of time. This is done even after
788 ! a restart, so that it takes into account any changes to max_iter.
789 call io_rm("td.general/laser", namespace=namespace)
790 call write_iter_init(writ%out(out_laser)%handle, 0, &
791 units_from_atomic(units_out%time, dt), &
792 trim(io_workpath("td.general/laser", namespace)))
793 do ii = 0, max_iter
794 call td_write_laser(writ%out(out_laser)%handle, space, list_get_lasers(ext_partners), dt, ii)
795 end do
796 call write_iter_end(writ%out(out_laser)%handle)
797 end if
798 end if
799
800 end if ! mpi_grp_is_root(mpi_world)
801
802 ! ----------------------------------------------
803 ! Initialize write that are MPI group-specific
804 ! ----------------------------------------------
805
806 if (writ%out(out_multipoles)%write .and. resolve_states) then
807 !resolve state contribution on multipoles
808 writ%out(out_multipoles)%hand_start = st%st_start
809 writ%out(out_multipoles)%hand_end = st%st_end
810 writ%out(out_multipoles)%resolve_states = .true.
811 writ%out(out_multipoles)%mpi_grp = gr%mpi_grp
812
813 safe_allocate(writ%out(out_multipoles)%mult_handles(st%st_start:st%st_end))
814
815 if (mpi_grp_is_root(writ%out(out_multipoles)%mpi_grp)) then
816 do ist = st%st_start, st%st_end
817 write(filename, '(a,i4.4)') 'td.general/multipoles-ist', ist
818 call write_iter_init(writ%out(out_multipoles)%mult_handles(ist), &
819 first, units_from_atomic(units_out%time, dt), &
820 trim(io_workpath(filename, namespace)))
821 end do
822 end if
823 end if
824
825 ! Misc operations
826
827 if (writ%out(out_total_current)%write .or. writ%out(out_total_heat_current)%write) then
828 !TODO: we should only compute the current here, not v_ks
829 call v_ks_calculate_current(ks, .true.)
830 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
831 calc_eigenval=.false., time = iter*dt, calc_energy = .false.)
832 end if
833
834 if (writ%out(out_n_ex)%write .and. writ%compute_interval > 0) then
835 call io_mkdir(outp%iter_dir, namespace)
836 end if
837
838 if (all(outp%how == 0) .and. writ%out(out_n_ex)%write) then
839 call io_function_read_what_how_when(namespace, space, outp%what, outp%how, outp%output_interval)
840 end if
841
842 ! --------------------------------------------------------
843 ! All steps of the above routine, but specific to DFT+U
844 ! --------------------------------------------------------
845
846 !%Variable TDOutputDFTU
847 !%Type flag
848 !%Default none
849 !%Section Time-Dependent::TD Output
850 !%Description
851 !% Defines what should be output during the time-dependent
852 !% simulation, related to DFT+U.
853 !%
854 !% Note: the output files generated by this option are updated
855 !% every <tt>RestartWriteInterval</tt> steps.
856 !%Option effective_u 1
857 !% Writes the effective U for each orbital set as a function of time.
858 !%End
859 default = 0
860 if (hm%lda_u_level == dft_u_acbn0) default = default + 2**(out_dftu_effective_u - 1)
861 call parse_variable(namespace, 'TDOutputDFTU', default, flags)
862
863 if (.not. varinfo_valid_option('TDOutputDFTU', flags, is_flag = .true.)) then
864 call messages_input_error(namespace, 'TDOutputDFTU')
865 end if
866
867 do iout = 1, out_dftu_max
868 writ%out_dftu(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
869 end do
870
871 if (mpi_grp_is_root(mpi_world)) then
872 if (writ%out_dftu(out_dftu_effective_u)%write) then
873 call write_iter_init(writ%out_dftu(out_dftu_effective_u)%handle, &
874 first, units_from_atomic(units_out%time, dt), &
875 trim(io_workpath("td.general/effectiveU", namespace)))
876 end if
877 end if
878
879 pop_sub(td_write_init)
880 end subroutine td_write_init
881
882
883 ! ---------------------------------------------------------
884 subroutine td_write_end(writ)
885 type(td_write_t), intent(inout) :: writ
886
887 integer :: ist, iout
888
889 push_sub(td_write_end)
890
891 do iout = 1, out_max
892 if (iout == out_laser) cycle
893 if (writ%out(iout)%write) then
894 if (mpi_grp_is_root(writ%out(iout)%mpi_grp)) then
895 if (writ%out(iout)%resolve_states) then
896 do ist = writ%out(iout)%hand_start, writ%out(iout)%hand_end
897 call write_iter_end(writ%out(iout)%mult_handles(ist))
898 end do
899 safe_deallocate_a(writ%out(iout)%mult_handles)
900 else
901 call write_iter_end(writ%out(iout)%handle)
902 end if
903 end if
904 end if
905 end do
906
907 if (mpi_grp_is_root(mpi_world)) then
908 do iout = 1, out_dftu_max
909 if (writ%out_dftu(iout)%write) call write_iter_end(writ%out_dftu(iout)%handle)
910 end do
911
912 end if
913
914 if (writ%out(out_populations)%write) then
915 do ist = 1, writ%n_excited_states
916 call excited_states_kill(writ%excited_st(ist))
917 end do
918 writ%n_excited_states = 0
919 end if
920
921 if (writ%out(out_proj)%write .or. writ%out(out_populations)%write &
922 .or. writ%out(out_n_ex)%write .or. writ%out(out_kp_proj)%write) then
923 call states_elec_end(writ%gs_st)
924 end if
925
926 pop_sub(td_write_end)
927 end subroutine td_write_end
928
929
930 ! ---------------------------------------------------------
931 subroutine td_write_iter(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, kick, ks, dt, iter, mc, recalculate_gs)
932 type(td_write_t), intent(inout) :: writ
933 type(namespace_t), intent(in) :: namespace
934 class(space_t), intent(in) :: space
935 type(output_t), intent(in) :: outp
936 type(grid_t), intent(in) :: gr
937 type(states_elec_t), intent(inout) :: st
938 type(hamiltonian_elec_t), intent(inout) :: hm
939 type(ions_t), intent(inout) :: ions
940 type(partner_list_t), intent(in) :: ext_partners
941 type(kick_t), intent(in) :: kick
942 type(v_ks_t), intent(in) :: ks
943 real(real64), intent(in) :: dt
944 integer, intent(in) :: iter
945 type(multicomm_t), intent(in) :: mc
946 logical, intent(in) :: recalculate_gs
947
948 type(gauge_field_t), pointer :: gfield
949 integer :: ierr
950 type(restart_t) :: restart_gs
951
952 push_sub_with_profile(td_write_iter)
953
954 if (writ%out(out_multipoles)%write) then
955 call td_write_multipole(writ%out(out_multipoles), space, gr, ions, st, writ%lmax, kick, iter)
956 end if
957
958 if (writ%out(out_ftchd)%write) then
959 call td_write_ftchd(writ%out(out_ftchd)%handle, space, gr, st, kick, iter)
960 end if
961
962 if (writ%out(out_angular)%write) then
963 call td_write_angular(writ%out(out_angular)%handle, namespace, space, gr, ions, hm, st, kick, iter)
964 end if
965
966 if (writ%out(out_spin)%write) then
967 call td_write_spin(writ%out(out_spin)%handle, gr, st, iter)
968 end if
969
970 if (writ%out(out_magnets)%write) then
971 call td_write_local_magnetic_moments(writ%out(out_magnets)%handle, gr, st, ions, writ%lmm_r, iter)
972 end if
973
974 if (writ%out(out_tot_m)%write) then
975 call td_write_tot_mag(writ%out(out_tot_m)%handle, gr, st, kick, iter)
976 end if
978 if (writ%out(out_proj)%write .and. mod(iter, writ%compute_interval) == 0) then
979 if (mpi_grp_is_root(mpi_world)) call write_iter_set(writ%out(out_proj)%handle, iter)
980 call td_write_proj(writ%out(out_proj)%handle, space, gr, ions, st, writ%gs_st, kick, iter)
981 end if
982
983 if (writ%out(out_floquet)%write) then
984 call td_write_floquet(namespace, space, hm, ext_partners, gr, st, iter)
985 end if
986
987 if (writ%out(out_kp_proj)%write .and. mod(iter, writ%compute_interval) == 0) then
988 call td_write_proj_kp(gr, hm%kpoints, st, writ%gs_st, namespace, iter)
989 end if
990
991 if (writ%out(out_coords)%write) then
992 call td_write_coordinates(writ%out(out_coords)%handle, ions%natoms, ions%space, &
993 ions%pos, ions%vel, ions%tot_force, iter)
994 end if
995
996 if (writ%out(out_separate_coords)%write) then
997 call td_write_sep_coordinates(writ%out(out_separate_coords)%handle, ions%natoms, ions%space, &
998 ions%pos, ions%vel, ions%tot_force, iter, 1)
999 end if
1000
1001 if (writ%out(out_separate_velocity)%write) then
1002 call td_write_sep_coordinates(writ%out(out_separate_velocity)%handle, ions%natoms, ions%space, &
1003 ions%pos, ions%vel, ions%tot_force, iter, 2)
1004 end if
1005
1006 if (writ%out(out_separate_forces)%write) then
1007 call td_write_sep_coordinates(writ%out(out_separate_forces)%handle, ions%natoms, ions%space, &
1008 ions%pos, ions%vel, ions%tot_force, iter, 3)
1009 end if
1010
1011 if (writ%out(out_temperature)%write) then
1012 call td_write_temperature(writ%out(out_temperature)%handle, ions, iter)
1013 end if
1014
1015 if (writ%out(out_populations)%write) then
1016 call td_write_populations(writ%out(out_populations)%handle, namespace, space, gr, st, writ, dt, iter)
1017 end if
1018
1019 if (writ%out(out_acc)%write) then
1020 call td_write_acc(writ%out(out_acc)%handle, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
1021 end if
1022
1023 if (writ%out(out_vel)%write) then
1024 call td_write_vel(writ%out(out_vel)%handle, gr, st, space, hm%kpoints, iter)
1025 end if
1026
1027 ! td_write_laser no longer called here, because the whole laser is printed
1028 ! out at the beginning.
1029
1030 if (writ%out(out_energy)%write) then
1031 call td_write_energy(writ%out(out_energy)%handle, hm, iter, ions%kinetic_energy)
1032 end if
1033
1034 if (writ%out(out_gauge_field)%write) then
1035 gfield => list_get_gauge_field(ext_partners)
1036 if(associated(gfield)) then
1037 call gauge_field_output_write(gfield, writ%out(out_gauge_field)%handle, iter)
1038 end if
1039 end if
1040
1041 if (writ%out(out_eigs)%write) then
1042 call td_write_eigs(writ%out(out_eigs)%handle, st, iter)
1043 end if
1044
1045 if (writ%out(out_ion_ch)%write) then
1046 call td_write_ionch(writ%out(out_ion_ch)%handle, gr, st, iter)
1047 end if
1048
1049 if (writ%out(out_total_current)%write) then
1050 call td_write_total_current(writ%out(out_total_current)%handle, space, gr, st, iter)
1051 end if
1052
1053 if (writ%out(out_total_heat_current)%write) then
1054 call td_write_total_heat_current(writ%out(out_total_heat_current)%handle, space, hm, gr, st, iter)
1055 end if
1056
1057 if (writ%out(out_partial_charges)%write) then
1058 call td_write_partial_charges(writ%out(out_partial_charges)%handle, gr, st, &
1059 ions, iter)
1060 end if
1061
1062 if (writ%out(out_n_ex)%write .and. mod(iter, writ%compute_interval) == 0) then
1063 if (mpi_grp_is_root(mpi_world)) call write_iter_set(writ%out(out_n_ex)%handle, iter)
1064 if (recalculate_gs) then ! Load recomputed GS states
1065 call restart_init(restart_gs, namespace, restart_proj, restart_type_load, mc, ierr, mesh=gr)
1066 call states_elec_load(restart_gs, namespace, space, writ%gs_st, gr, hm%kpoints, &
1067 ierr, label = ': Houston states for TDOutput')
1069 end if
1070
1071 call td_write_n_ex(writ%out(out_n_ex)%handle, outp, namespace, gr, hm%kpoints, st, writ%gs_st, iter)
1072 end if
1073
1074 if (writ%out(out_norm_ks)%write) then
1075 call td_write_norm_ks_orbitals(writ%out(out_norm_ks)%handle, gr, hm%kpoints, st, iter)
1076 end if
1077
1078 !DFT+U outputs
1079 if (writ%out_dftu(out_dftu_effective_u)%write) then
1080 call td_write_effective_u(writ%out_dftu(out_dftu_effective_u)%handle, hm%lda_u, iter)
1081 end if
1082
1083 if (writ%out(out_q)%write .and. ks%has_photons) then
1084 call td_write_q(writ%out(out_q)%handle, space, ks, iter)
1085 end if
1086
1087 if (writ%out(out_mxll_field)%write .and. hm%mxll%calc_field_dip) then
1088 call td_write_mxll_field(writ%out(out_mxll_field)%handle, space, hm, dt, iter)
1089 end if
1090
1091 pop_sub_with_profile(td_write_iter)
1092 end subroutine td_write_iter
1093
1094
1095 ! ---------------------------------------------------------
1096 subroutine td_write_data(writ)
1097 type(td_write_t), intent(inout) :: writ
1098
1099 integer :: iout, ii
1100
1101 push_sub(td_write_data)
1102 call profiling_in("TD_WRITE_DATA")
1103
1104 do iout = 1, out_max
1105 if (iout == out_laser) cycle
1106 if (writ%out(iout)%write) then
1107 if (mpi_grp_is_root(writ%out(iout)%mpi_grp)) then
1108 if (writ%out(iout)%resolve_states) then
1109 do ii = writ%out(iout)%hand_start, writ%out(iout)%hand_end
1110 call write_iter_flush(writ%out(iout)%mult_handles(ii))
1111 end do
1112 else
1113 call write_iter_flush(writ%out(iout)%handle)
1114 end if
1115 end if
1116 end if
1117 end do
1118
1119 if (mpi_grp_is_root(mpi_world)) then
1120 do iout = 1, out_dftu_max
1121 if (writ%out_dftu(iout)%write) call write_iter_flush(writ%out_dftu(iout)%handle)
1122 end do
1123 end if
1124
1125 call profiling_out("TD_WRITE_DATA")
1126 pop_sub(td_write_data)
1127 end subroutine td_write_data
1128
1129 ! ---------------------------------------------------------
1130 subroutine td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, dt)
1131 type(namespace_t), intent(in) :: namespace
1132 type(electron_space_t), intent(in) :: space
1133 type(grid_t), intent(in) :: gr
1134 type(states_elec_t), intent(inout) :: st
1135 type(hamiltonian_elec_t), intent(inout) :: hm
1136 type(v_ks_t), intent(inout) :: ks
1137 type(output_t), intent(in) :: outp
1138 type(ions_t), intent(in) :: ions
1139 type(partner_list_t), intent(in) :: ext_partners
1140 integer, intent(in) :: iter
1141 real(real64), optional, intent(in) :: dt
1142
1143 character(len=256) :: filename
1144
1145 push_sub(td_write_output)
1146 call profiling_in("TD_WRITE_OUTPUT")
1147
1148 ! TODO this now overwrites wf inside st. If this is not wanted need to add an optional overwrite=no flag
1149 if (st%modelmbparticles%nparticle > 0) then
1150 call modelmb_sym_all_states(space, gr, st)
1151 end if
1152
1153 ! now write down the rest
1154 write(filename, '(a,a,i7.7)') trim(outp%iter_dir),"td.", iter ! name of directory
1155
1156 call output_all(outp, namespace, space, filename, gr, ions, iter, st, hm, ks)
1157
1158 call output_modelmb(outp, namespace, space, filename, gr, ions, iter, st)
1159
1160 if (present(dt)) then
1161 call output_scalar_pot(outp, namespace, space, filename, gr, ions, ext_partners, iter*dt)
1162 else
1163 if (iter == 0) call output_scalar_pot(outp, namespace, space, filename, gr, ions, ext_partners)
1164 end if
1165
1166 call profiling_out("TD_WRITE_OUTPUT")
1167 pop_sub(td_write_output)
1168 end subroutine td_write_output
1169
1170 ! ---------------------------------------------------------
1171 subroutine td_write_spin(out_spin, mesh, st, iter)
1172 type(c_ptr), intent(inout) :: out_spin
1173 class(mesh_t), intent(in) :: mesh
1174 type(states_elec_t), intent(in) :: st
1175 integer, intent(in) :: iter
1176
1177 character(len=130) :: aux
1178 real(real64) :: spin(3)
1179
1180 push_sub(td_write_spin)
1181
1182 ! The expectation value of the spin operator is half the total magnetic moment
1183 ! This has to be calculated by all nodes
1184 call magnetic_moment(mesh, st, st%rho, spin)
1185 spin = m_half*spin
1186
1187 if (mpi_grp_is_root(mpi_world)) then ! only first node outputs
1188
1189 if (iter == 0) then
1191
1192 !second line -> columns name
1194 if (st%d%ispin == spinors) then
1195 write(aux, '(a2,18x)') 'Sx'
1196 call write_iter_header(out_spin, aux)
1197 write(aux, '(a2,18x)') 'Sy'
1198 call write_iter_header(out_spin, aux)
1199 end if
1200 write(aux, '(a2,18x)') 'Sz'
1201 call write_iter_header(out_spin, aux)
1203
1205 end if
1206
1208 select case (st%d%ispin)
1209 case (spin_polarized)
1210 call write_iter_double(out_spin, spin(3), 1)
1211 case (spinors)
1212 call write_iter_double(out_spin, spin(1:3), 3)
1213 end select
1215
1216 end if
1217
1218 pop_sub(td_write_spin)
1219 end subroutine td_write_spin
1220
1221
1222 ! ---------------------------------------------------------
1223 subroutine td_write_local_magnetic_moments(out_magnets, gr, st, ions, lmm_r, iter)
1224 type(c_ptr), intent(inout) :: out_magnets
1225 type(grid_t), intent(in) :: gr
1226 type(states_elec_t), intent(in) :: st
1227 type(ions_t), intent(in) :: ions
1228 real(real64), intent(in) :: lmm_r
1229 integer, intent(in) :: iter
1230
1231 integer :: ia
1232 character(len=50) :: aux
1233 real(real64), allocatable :: lmm(:,:)
1234
1236
1237 !get the atoms` magnetization. This has to be calculated by all nodes
1238 safe_allocate(lmm(1:3, 1:ions%natoms))
1239 call magnetic_local_moments(gr, st, ions, gr%der%boundaries, st%rho, lmm_r, lmm)
1240
1241 if (mpi_grp_is_root(mpi_world)) then ! only first node outputs
1242
1243 if (iter == 0) then
1245
1246 !second line -> columns name
1248 do ia = 1, ions%natoms
1249 if (st%d%ispin == spinors) then
1250 write(aux, '(a2,i2.2,16x)') 'mx', ia
1252 write(aux, '(a2,i2.2,16x)') 'my', ia
1254 end if
1255 write(aux, '(a2,i2.2,16x)') 'mz', ia
1257 end do
1259
1261 end if
1262
1264 do ia = 1, ions%natoms
1265 select case (st%d%ispin)
1266 case (spin_polarized)
1267 call write_iter_double(out_magnets, lmm(3, ia), 1)
1268 case (spinors)
1269 call write_iter_double(out_magnets, lmm(1:3, ia), 3)
1270 end select
1271 end do
1273 safe_deallocate_a(lmm)
1274 end if
1275
1277 end subroutine td_write_local_magnetic_moments
1278
1279 ! ---------------------------------------------------------
1280 subroutine td_write_tot_mag(out_magnets, mesh, st, kick, iter)
1281 type(c_ptr), intent(inout) :: out_magnets
1282 class(mesh_t), intent(in) :: mesh
1283 type(states_elec_t), intent(in) :: st
1284 type(kick_t), intent(in) :: kick
1285 integer, intent(in) :: iter
1286
1287 complex(real64), allocatable :: tm(:,:)
1288 integer :: ii, iq
1289
1290 push_sub(td_write_tot_mag)
1291
1292 safe_allocate(tm(1:6,1:kick%nqvec))
1293
1294 do iq = 1, kick%nqvec
1295 call magnetic_total_magnetization(mesh, st, kick%qvector(:,iq), tm(1:6,iq))
1296 end do
1297
1298 if (mpi_grp_is_root(mpi_world)) then ! only first node outputs
1299
1300 if (iter == 0) then
1301 call td_write_print_header_init(out_magnets)
1302 call kick_write(kick, out = out_magnets)
1303
1304 !second line -> columns name
1305 call write_iter_header_start(out_magnets)
1306 call write_iter_header(out_magnets, 'Re[m_x(q)]')
1307 call write_iter_header(out_magnets, 'Im[m_x(q)]')
1308 call write_iter_header(out_magnets, 'Re[m_y(q)]')
1309 call write_iter_header(out_magnets, 'Im[m_y(q)]')
1310 call write_iter_header(out_magnets, 'Re[m_z(q)]')
1311 call write_iter_header(out_magnets, 'Im[m_z(q)]')
1312 call write_iter_header(out_magnets, 'Re[m_x(-q)]')
1313 call write_iter_header(out_magnets, 'Im[m_x(-q)]')
1314 call write_iter_header(out_magnets, 'Re[m_y(-q)]')
1315 call write_iter_header(out_magnets, 'Im[m_y(-q)]')
1316 call write_iter_header(out_magnets, 'Re[m_z(-q)]')
1317 call write_iter_header(out_magnets, 'Im[m_z(-q)]')
1318 call write_iter_nl(out_magnets)
1319
1320 call td_write_print_header_end(out_magnets)
1321 end if
1322
1323 call write_iter_start(out_magnets)
1324 do iq = 1, kick%nqvec
1325 do ii = 1, 6
1326 call write_iter_double(out_magnets, real(tm(ii, iq), real64), 1)
1327 call write_iter_double(out_magnets, aimag(tm(ii, iq)), 1)
1328 end do
1329 end do
1330 call write_iter_nl(out_magnets)
1331 end if
1332
1333 safe_deallocate_a(tm)
1334
1335 pop_sub(td_write_tot_mag)
1336 end subroutine td_write_tot_mag
1337
1338
1339 ! ---------------------------------------------------------
1347 subroutine td_write_angular(out_angular, namespace, space, gr, ions, hm, st, kick, iter)
1348 type(c_ptr), intent(inout) :: out_angular
1349 type(namespace_t), intent(in) :: namespace
1350 class(space_t), intent(in) :: space
1351 type(grid_t), intent(in) :: gr
1352 type(ions_t), intent(inout) :: ions
1353 type(hamiltonian_elec_t), intent(inout) :: hm
1354 type(states_elec_t), intent(inout) :: st
1355 type(kick_t), intent(in) :: kick
1356 integer, intent(in) :: iter
1357
1358 integer :: idir
1359 character(len=130) :: aux
1360 real(real64) :: angular(3)
1361 class(perturbation_magnetic_t), pointer :: angular_momentum
1362
1363 push_sub(td_write_angular)
1364
1365 angular_momentum => perturbation_magnetic_t(namespace, ions)
1366 do idir = 1, 3
1367 call angular_momentum%setup_dir(idir)
1368 !we have to multiply by 2, because the perturbation returns L/2
1369 angular(idir) = &
1370 m_two*real(angular_momentum%zstates_elec_expectation_value(namespace, space, gr, hm, st), real64)
1371 end do
1372 safe_deallocate_p(angular_momentum)
1374 if (mpi_grp_is_root(mpi_world)) then ! Only first node outputs
1375
1376 if (iter == 0) then
1377 call td_write_print_header_init(out_angular)
1378
1379 write(aux, '(a15,i2)') '# nspin ', st%d%nspin
1380 call write_iter_string(out_angular, aux)
1381 call write_iter_nl(out_angular)
1382
1383 call kick_write(kick, out = out_angular)
1384
1385 !second line -> columns name
1386 call write_iter_header_start(out_angular)
1387 write(aux, '(a4,18x)') '<Lx>'
1388 call write_iter_header(out_angular, aux)
1389 write(aux, '(a4,18x)') '<Ly>'
1390 call write_iter_header(out_angular, aux)
1391 write(aux, '(a4,18x)') '<Lz>'
1392 call write_iter_header(out_angular, aux)
1393 call write_iter_nl(out_angular)
1394
1395 !third line -> should hold the units.
1396 call write_iter_string(out_angular, '#[Iter n.]')
1397 call write_iter_header(out_angular, '[' // trim(units_abbrev(units_out%time)) // ']')
1398 call write_iter_header(out_angular, '[hbar]')
1399 call write_iter_header(out_angular, '[hbar]')
1400 call write_iter_header(out_angular, '[hbar]')
1401 call write_iter_nl(out_angular)
1402
1403 call td_write_print_header_end(out_angular)
1404 end if
1405
1406 call write_iter_start(out_angular)
1407 call write_iter_double(out_angular, angular(1:3), 3)
1408 call write_iter_nl(out_angular)
1409
1410 end if
1411
1412 pop_sub(td_write_angular)
1413 end subroutine td_write_angular
1414
1415
1418 subroutine td_write_multipole(out_multip, space, gr, ions, st, lmax, kick, iter)
1419 type(td_write_prop_t), intent(inout) :: out_multip
1420 class(space_t), intent(in) :: space
1421 type(grid_t), intent(in) :: gr
1422 type(ions_t), intent(in) :: ions
1423 type(states_elec_t), intent(in) :: st
1424 integer, intent(in) :: lmax
1425 type(kick_t), intent(in) :: kick
1426 integer, intent(in) :: iter
1427
1428 integer :: ist
1429 real(real64), allocatable :: rho(:,:)
1430
1431 push_sub(td_write_multipole)
1432
1433 if (out_multip%resolve_states) then
1434 safe_allocate(rho(1:gr%np_part, 1:st%d%nspin))
1435 rho(:,:) = m_zero
1436
1437 do ist = st%st_start, st%st_end
1438 call density_calc(st, gr, rho, istin = ist)
1439 call td_write_multipole_r(out_multip%mult_handles(ist), space, gr, ions, st, lmax, kick, rho, iter, &
1440 mpi_grp = out_multip%mpi_grp)
1441 end do
1442
1443 safe_deallocate_a(rho)
1444
1445 else
1446 if (allocated(st%frozen_rho)) then
1447 safe_allocate(rho(1:gr%np, 1:st%d%nspin))
1448 call lalg_copy(gr%np, st%d%nspin, st%rho, rho)
1449 call lalg_axpy(gr%np, st%d%nspin, m_one, st%frozen_rho, rho)
1450
1451 call td_write_multipole_r(out_multip%handle, space, gr, ions, st, lmax, kick, rho, iter)
1452
1453 safe_deallocate_a(rho)
1454 else
1455 call td_write_multipole_r(out_multip%handle, space, gr, ions, st, lmax, kick, st%rho, iter)
1456 end if
1457
1458 end if
1459
1460 pop_sub(td_write_multipole)
1461 end subroutine td_write_multipole
1462
1463
1465 subroutine td_write_multipole_r(out_multip, space, mesh, ions, st, lmax, kick, rho, iter, mpi_grp)
1466 type(c_ptr), intent(inout) :: out_multip
1467 class(space_t), intent(in) :: space
1468 class(mesh_t), intent(in) :: mesh
1469 type(ions_t), intent(in) :: ions
1470 type(states_elec_t), intent(in) :: st
1471 integer, intent(in) :: lmax
1472 type(kick_t), intent(in) :: kick
1473 real(real64), intent(in) :: rho(:,:)
1474 integer, intent(in) :: iter
1475 type(mpi_grp_t), optional, intent(in) :: mpi_grp
1476
1477
1478 integer :: is, idir, ll, mm, add_lm
1479 character(len=120) :: aux
1480 real(real64) :: ionic_dipole(ions%space%dim)
1481 real(real64), allocatable :: multipole(:,:)
1482 type(mpi_grp_t) :: mpi_grp_
1483
1484 push_sub(td_write_multipole_r)
1485
1486 ! We cannot output multipoles beyond the dipole for higher dimensions
1487 assert(.not. (lmax > 1 .and. space%dim > 3))
1488
1489 mpi_grp_ = mpi_world
1490 if (present(mpi_grp)) mpi_grp_ = mpi_grp
1491
1492 if (mpi_grp_is_root(mpi_grp_).and.iter == 0) then
1493 call td_write_print_header_init(out_multip)
1494
1495 write(aux, '(a15,i2)') '# nspin ', st%d%nspin
1496 call write_iter_string(out_multip, aux)
1497 call write_iter_nl(out_multip)
1498
1499 write(aux, '(a15,i2)') '# lmax ', lmax
1500 call write_iter_string(out_multip, aux)
1501 call write_iter_nl(out_multip)
1502
1503 call kick_write(kick, out = out_multip)
1504
1505 call write_iter_header_start(out_multip)
1506
1507 do is = 1, st%d%nspin
1508 write(aux,'(a18,i1,a1)') 'Electronic charge(', is,')'
1509 call write_iter_header(out_multip, aux)
1510 if (lmax > 0) then
1511 do idir = 1, space%dim
1512 write(aux, '(4a1,i1,a1)') '<', index2axis(idir), '>', '(', is,')'
1513 call write_iter_header(out_multip, aux)
1514 end do
1515 end if
1516 do ll = 2, lmax
1517 do mm = -ll, ll
1518 write(aux, '(a2,i2,a4,i2,a2,i1,a1)') 'l=', ll, ', m=', mm, ' (', is,')'
1519 call write_iter_header(out_multip, aux)
1520 end do
1521 end do
1522 end do
1523 call write_iter_nl(out_multip)
1524
1525 ! units
1526 call write_iter_string(out_multip, '#[Iter n.]')
1527 call write_iter_header(out_multip, '[' // trim(units_abbrev(units_out%time)) // ']')
1528
1529 do is = 1, st%d%nspin
1530 call write_iter_header(out_multip, 'Electrons')
1531 if (lmax > 0) then
1532 do idir = 1, space%dim
1533 call write_iter_header(out_multip, '[' // trim(units_abbrev(units_out%length)) // ']')
1534 end do
1535 end if
1536 do ll = 2, lmax
1537 do mm = -ll, ll
1538 write(aux, '(a,a2,i1)') trim(units_abbrev(units_out%length)), "**", ll
1539 call write_iter_header(out_multip, '[' // trim(aux) // ']')
1540 end do
1541 end do
1542 end do
1543 call write_iter_nl(out_multip)
1544
1545 call td_write_print_header_end(out_multip)
1546 end if
1547
1548 if (space%dim > 3 .and. lmax == 1) then
1549 ! For higher dimensions we can only have up to the dipole
1550 safe_allocate(multipole(1:space%dim+1, 1:st%d%nspin))
1551 else
1552 safe_allocate(multipole(1:(lmax + 1)**2, 1:st%d%nspin))
1553 end if
1554 multipole(:,:) = m_zero
1555
1556 do is = 1, st%d%nspin
1557 call dmf_multipoles(mesh, rho(:,is), lmax, multipole(:,is))
1558 end do
1559
1560 if (lmax > 0) then
1561 ionic_dipole = ions%dipole()
1562 do is = 1, st%d%nspin
1563 multipole(2:space%dim+1, is) = -ionic_dipole(1:space%dim)/st%d%nspin - multipole(2:space%dim+1, is)
1564 end do
1565 end if
1566
1567 if (mpi_grp_is_root(mpi_grp_)) then
1568 call write_iter_start(out_multip)
1569 do is = 1, st%d%nspin
1570 call write_iter_double(out_multip, units_from_atomic(units_out%length**0, multipole(1, is)), 1)
1571 if (lmax > 0) then
1572 do idir = 1, space%dim
1573 call write_iter_double(out_multip, units_from_atomic(units_out%length, multipole(1+idir, is)), 1)
1574 end do
1575 end if
1576 add_lm = space%dim + 2
1577 do ll = 2, lmax
1578 do mm = -ll, ll
1579 call write_iter_double(out_multip, units_from_atomic(units_out%length**ll, multipole(add_lm, is)), 1)
1580 add_lm = add_lm + 1
1581 end do
1582 end do
1583 end do
1584 call write_iter_nl(out_multip)
1585 end if
1586
1587 safe_deallocate_a(multipole)
1588 pop_sub(td_write_multipole_r)
1589 end subroutine td_write_multipole_r
1590
1591 ! ---------------------------------------------------------
1592 subroutine td_write_ftchd(out_ftchd, space, mesh, st, kick, iter)
1593 type(c_ptr), intent(inout) :: out_ftchd
1594 class(space_t), intent(in) :: space
1595 class(mesh_t), intent(in) :: mesh
1596 type(states_elec_t), intent(in) :: st
1597 type(kick_t), intent(in) :: kick
1598 integer, intent(in) :: iter
1599
1600 integer :: is, ip, idir
1601 character(len=120) :: aux, aux2
1602 real(real64) :: ftchd_bessel
1603 complex(real64) :: ftchd
1604 real(real64) :: ylm
1605 real(real64), allocatable :: integrand_bessel(:)
1606 complex(real64), allocatable :: integrand(:)
1607
1608 push_sub(td_write_ftchd)
1609
1610 if (mpi_grp_is_root(mpi_world).and.iter == 0) then
1611 call td_write_print_header_init(out_ftchd)
1612
1613 write(aux,'(a15, i2)') '# qkickmode ', kick%qkick_mode
1614 call write_iter_string(out_ftchd, aux)
1615 call write_iter_nl(out_ftchd)
1616
1617 if (kick%qkick_mode == qkickmode_bessel) then
1618 write(aux,'(a15, i0.3, 1x, i0.3)') '# ll, mm ', kick%qbessel_l, kick%qbessel_m
1619 call write_iter_string(out_ftchd, aux)
1620 call write_iter_nl(out_ftchd)
1621 end if
1622
1623 if (kick%qkick_mode == qkickmode_bessel) then
1624 write(aux, '(a15, f9.6)') '# qlength ', kick%qlength
1625 else ! sin or cos
1626 write(aux, '(a15)') '# qvector '
1627 do idir = 1, space%dim
1628 write(aux2, '(f9.5)') kick%qvector(idir,1)
1629 aux = trim(aux) // trim(aux2)
1630 end do
1631 end if
1632 call write_iter_string(out_ftchd, aux)
1633 call write_iter_nl(out_ftchd)
1634
1635 write(aux, '(a15,f18.12)') '# kick strength', kick%delta_strength
1636 call write_iter_string(out_ftchd, aux)
1637 call write_iter_nl(out_ftchd)
1638
1639 call write_iter_header_start(out_ftchd)
1640 if (kick%qkick_mode == qkickmode_bessel) then
1641 write(aux,'(a17)') 'int(j_l*Y_lm*rho)'
1642 else
1643 write(aux,'(a12)') 'Real, Imag'
1644 end if
1645 call write_iter_header(out_ftchd, aux)
1646 call write_iter_nl(out_ftchd)
1647
1648 ! units
1649 call write_iter_string(out_ftchd, '#[Iter n.]')
1650 call write_iter_header(out_ftchd, '[' // trim(units_abbrev(units_out%time)) // ']')
1651 call write_iter_nl(out_ftchd)
1652 call td_write_print_header_end(out_ftchd)
1653
1654 end if
1655
1656 ftchd = m_zero
1657
1658 ! If kick mode is exp, sin, or cos, apply the normal Fourier transform
1659 if (kick%qkick_mode /= qkickmode_bessel) then
1660 safe_allocate(integrand(1:mesh%np))
1661 integrand = m_zero
1662 do is = 1, st%d%nspin
1663 do ip = 1, mesh%np
1664 integrand(ip) = integrand(ip) + st%rho(ip, is) * exp(-m_zi*sum(mesh%x(ip, 1:space%dim)*kick%qvector(1:space%dim, 1)))
1665 end do
1666 end do
1667 ftchd = zmf_integrate(mesh, integrand)
1668 safe_deallocate_a(integrand)
1669 else
1670 ftchd_bessel = m_zero
1671 safe_allocate(integrand_bessel(1:mesh%np))
1672 integrand_bessel = m_zero
1673 do is = 1, st%d%nspin
1674 do ip = 1, mesh%np
1675 call ylmr_real(mesh%x(ip, 1:3), kick%qbessel_l, kick%qbessel_m, ylm)
1676 integrand_bessel(ip) = integrand_bessel(ip) + st%rho(ip, is) * &
1677 loct_sph_bessel(kick%qbessel_l, kick%qlength*norm2(mesh%x(ip, :)))*ylm
1678 end do
1679 end do
1680 ftchd_bessel = dmf_integrate(mesh, integrand_bessel)
1681 safe_deallocate_a(integrand_bessel)
1682 end if
1683
1684 if (mpi_grp_is_root(mpi_world)) then
1685 call write_iter_start(out_ftchd)
1686 if (kick%qkick_mode == qkickmode_bessel) then
1687 call write_iter_double(out_ftchd, ftchd_bessel, 1)
1688 else ! exp, sin, cos
1689 call write_iter_double(out_ftchd, real(ftchd), 1)
1690 call write_iter_double(out_ftchd, aimag(ftchd), 1)
1691 end if
1692 call write_iter_nl(out_ftchd)
1693 end if
1694
1695 pop_sub(td_write_ftchd)
1696 end subroutine td_write_ftchd
1697
1698 ! ---------------------------------------------------------
1699 subroutine td_write_temperature(out_temperature, ions, iter)
1700 type(c_ptr), intent(inout) :: out_temperature
1701 type(ions_t), intent(in) :: ions
1702 integer, intent(in) :: iter
1703
1704 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
1705
1706 push_sub(td_write_temperature)
1707
1708 if (iter == 0) then
1709 call td_write_print_header_init(out_temperature)
1710
1711 ! first line: column names
1712 call write_iter_header_start(out_temperature)
1713 call write_iter_header(out_temperature, 'Temperature')
1714 call write_iter_nl(out_temperature)
1715
1716 ! second line: units
1717 call write_iter_string(out_temperature, '#[Iter n.]')
1718 call write_iter_header(out_temperature, '[' // trim(units_abbrev(units_out%time)) // ']')
1719 call write_iter_string(out_temperature, ' [K]')
1720 call write_iter_nl(out_temperature)
1721
1722 call td_write_print_header_end(out_temperature)
1723 end if
1724
1725 call write_iter_start(out_temperature)
1726
1728
1729 call write_iter_nl(out_temperature)
1730
1731 pop_sub(td_write_temperature)
1732 end subroutine td_write_temperature
1733
1734
1735 ! ---------------------------------------------------------
1736 subroutine td_write_populations(out_populations, namespace, space, mesh, st, writ, dt, iter)
1737 type(c_ptr), intent(inout) :: out_populations
1738 type(namespace_t), intent(in) :: namespace
1739 class(space_t), intent(in) :: space
1740 class(mesh_t), intent(in) :: mesh
1741 type(states_elec_t), intent(inout) :: st
1742 type(td_write_t), intent(in) :: writ
1743 real(real64), intent(in) :: dt
1744 integer, intent(in) :: iter
1745
1746 integer :: ist
1747 character(len=6) :: excited_name
1748 complex(real64) :: gsp
1749 complex(real64), allocatable :: excited_state_p(:)
1750 complex(real64), allocatable :: dotprodmatrix(:, :, :)
1751
1752
1753 push_sub(td_write_populations)
1754
1755 safe_allocate(dotprodmatrix(1:writ%gs_st%nst, 1:st%nst, 1:st%nik))
1756 call zstates_elec_matrix(writ%gs_st, st, mesh, dotprodmatrix)
1757
1758
1759 !See comment in zstates_elec_mpdotp
1760 assert(.not. space%is_periodic())
1761
1762 ! all processors calculate the projection
1763 gsp = zstates_elec_mpdotp(namespace, mesh, writ%gs_st, st, dotprodmatrix)
1764
1765 if (writ%n_excited_states > 0) then
1766 safe_allocate(excited_state_p(1:writ%n_excited_states))
1767 do ist = 1, writ%n_excited_states
1768 excited_state_p(ist) = zstates_elec_mpdotp(namespace, mesh, writ%excited_st(ist), st, dotprodmatrix)
1769 end do
1770 end if
1771
1772 if (mpi_grp_is_root(mpi_world)) then
1773 if (iter == 0) then
1774 call td_write_print_header_init(out_populations)
1775
1776 ! first line -> column names
1777 call write_iter_header_start(out_populations)
1778 call write_iter_header(out_populations, 'Re<Phi_gs|Phi(t)>')
1779 call write_iter_header(out_populations, 'Im<Phi_gs|Phi(t)>')
1780 do ist = 1, writ%n_excited_states
1781 write(excited_name,'(a2,i3,a1)') 'P(', ist,')'
1782 call write_iter_header(out_populations, 'Re<'//excited_name//'|Phi(t)>')
1783 call write_iter_header(out_populations, 'Im<'//excited_name//'|Phi(t)>')
1784 end do
1785 call write_iter_nl(out_populations)
1786
1787 ! second line -> units
1788 call write_iter_string(out_populations, '#[Iter n.]')
1789 call write_iter_header(out_populations, '[' // trim(units_abbrev(units_out%time)) // ']')
1790 call write_iter_nl(out_populations)
1791
1792 call td_write_print_header_end(out_populations)
1793 end if
1794
1795 ! cannot call write_iter_start, for the step is not 1
1796 call write_iter_int(out_populations, iter, 1)
1797 call write_iter_double(out_populations, units_from_atomic(units_out%time, iter*dt), 1)
1798 call write_iter_double(out_populations, real(gsp), 1)
1799 call write_iter_double(out_populations, aimag(gsp), 1)
1800 do ist = 1, writ%n_excited_states
1801 call write_iter_double(out_populations, real(excited_state_p(ist)), 1)
1802 call write_iter_double(out_populations, aimag(excited_state_p(ist)), 1)
1803 end do
1804 call write_iter_nl(out_populations)
1805 end if
1806
1807 if (writ%n_excited_states > 0) then
1808 safe_deallocate_a(excited_state_p)
1809 end if
1810 safe_deallocate_a(dotprodmatrix)
1811 pop_sub(td_write_populations)
1812 end subroutine td_write_populations
1813
1814
1815 ! ---------------------------------------------------------
1816 subroutine td_write_acc(out_acc, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
1817 type(c_ptr), intent(inout) :: out_acc
1818 type(namespace_t), intent(in) :: namespace
1819 class(space_t), intent(in) :: space
1820 type(grid_t), intent(in) :: gr
1821 type(ions_t), intent(inout) :: ions
1822 type(states_elec_t), intent(inout) :: st
1823 type(hamiltonian_elec_t), intent(inout) :: hm
1824 type(partner_list_t), intent(in) :: ext_partners
1825 real(real64), intent(in) :: dt
1826 integer, intent(in) :: iter
1827
1828 integer :: idim
1829 character(len=7) :: aux
1830 real(real64) :: acc(space%dim)
1831
1832 push_sub(td_write_acc)
1833
1834 if (iter == 0 .and. mpi_grp_is_root(mpi_world)) then
1835 call td_write_print_header_init(out_acc)
1836
1837 ! first line -> column names
1838 call write_iter_header_start(out_acc)
1839 do idim = 1, space%dim
1840 write(aux, '(a4,i1,a1)') 'Acc(', idim, ')'
1841 call write_iter_header(out_acc, aux)
1842 end do
1843 call write_iter_nl(out_acc)
1844
1845 ! second line: units
1846 call write_iter_string(out_acc, '#[Iter n.]')
1847 call write_iter_header(out_acc, '[' // trim(units_abbrev(units_out%time)) // ']')
1848 do idim = 1, space%dim
1849 call write_iter_header(out_acc, '[' // trim(units_abbrev(units_out%acceleration)) // ']')
1850 end do
1851 call write_iter_nl(out_acc)
1852 call td_write_print_header_end(out_acc)
1853 end if
1854
1855 call td_calc_tacc(namespace, space, gr, ions, ext_partners, st, hm, acc, dt*iter)
1856
1857 if (mpi_grp_is_root(mpi_world)) then
1858 call write_iter_start(out_acc)
1859 acc = units_from_atomic(units_out%acceleration, acc)
1860 call write_iter_double(out_acc, acc, space%dim)
1861 call write_iter_nl(out_acc)
1862 end if
1863
1864 pop_sub(td_write_acc)
1865 end subroutine td_write_acc
1866
1867 ! ---------------------------------------------------------
1868 subroutine td_write_vel(out_vel, gr, st, space, kpoints, iter)
1869 type(c_ptr), intent(inout) :: out_vel
1870 type(grid_t), intent(in) :: gr
1871 type(states_elec_t), intent(in) :: st
1872 type(space_t), intent(in) :: space
1873 type(kpoints_t), intent(in) :: kpoints
1874 integer, intent(in) :: iter
1875
1876 integer :: idim
1877 character(len=7) :: aux
1878 real(real64) :: vel(space%dim)
1879
1880 push_sub(td_write_vel)
1881
1882 if (iter == 0 .and. mpi_grp_is_root(mpi_world)) then
1883 call td_write_print_header_init(out_vel)
1884
1885 ! first line -> column names
1886 call write_iter_header_start(out_vel)
1887 do idim = 1, space%dim
1888 write(aux, '(a4,i1,a1)') 'Vel(', idim, ')'
1889 call write_iter_header(out_vel, aux)
1890 end do
1891 call write_iter_nl(out_vel)
1892
1893 ! second line: units
1894 call write_iter_string(out_vel, '#[Iter n.]')
1895 call write_iter_header(out_vel, '[' // trim(units_abbrev(units_out%time)) // ']')
1896 do idim = 1, space%dim
1897 call write_iter_header(out_vel, '[' // trim(units_abbrev(units_out%velocity)) // ']')
1898 end do
1899 call write_iter_nl(out_vel)
1900 call td_write_print_header_end(out_vel)
1901 end if
1902
1903 call td_calc_tvel(gr, st, space, kpoints, vel)
1904
1905 if (mpi_grp_is_root(mpi_world)) then
1906 call write_iter_start(out_vel)
1907 vel = units_from_atomic(units_out%velocity, vel)
1908 call write_iter_double(out_vel, vel, space%dim)
1909 call write_iter_nl(out_vel)
1910 end if
1911
1912 pop_sub(td_write_vel)
1913 end subroutine td_write_vel
1914
1915
1916 ! ---------------------------------------------------------
1917 subroutine td_write_laser(out_laser, space, lasers, dt, iter)
1918 type(c_ptr), intent(inout) :: out_laser
1919 class(space_t), intent(in) :: space
1920 type(lasers_t), intent(inout) :: lasers
1921 real(real64), intent(in) :: dt
1922 integer, intent(in) :: iter
1923
1924 integer :: il, idir
1925 real(real64) :: field(space%dim)
1926 real(real64) :: ndfield(space%dim)
1927 character(len=80) :: aux
1928
1929 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
1930
1931 ! no PUSH SUB, called too often
1932
1933 if (iter == 0) then
1934 call td_write_print_header_init(out_laser)
1935
1936 ! first line
1937 write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, dt), &
1938 " [", trim(units_abbrev(units_out%time)), "]"
1939 call write_iter_string(out_laser, aux)
1940 call write_iter_nl(out_laser)
1941
1942 call write_iter_header_start(out_laser)
1943 do il = 1, lasers%no_lasers
1944 select case (laser_kind(lasers%lasers(il)))
1945 case (e_field_electric)
1946 do idir = 1, space%dim
1947 write(aux, '(a,i1,a)') 'E(', idir, ')'
1948 call write_iter_header(out_laser, aux)
1949 end do
1950 case (e_field_magnetic)
1951 do idir = 1, space%dim
1952 write(aux, '(a,i1,a)') 'B(', idir, ')'
1953 call write_iter_header(out_laser, aux)
1954 end do
1956 do idir = 1, space%dim
1957 write(aux, '(a,i1,a)') 'A(', idir, ')'
1958 call write_iter_header(out_laser, aux)
1959 end do
1961 write(aux, '(a,i1,a)') 'e(t)'
1962 call write_iter_header(out_laser, aux)
1963 end select
1964 end do
1965
1966 if (lasers_with_nondipole_field(lasers)) then
1967 do idir = 1, space%dim
1968 write(aux, '(a,i1,a)') 'A^M(', idir, ')'
1969 call write_iter_header(out_laser, aux)
1970 end do
1971 end if
1972 call write_iter_nl(out_laser)
1973
1974 call write_iter_string(out_laser, '#[Iter n.]')
1975 call write_iter_header(out_laser, '[' // trim(units_abbrev(units_out%time)) // ']')
1976
1977 ! Note that we do not print out units of E, B, or A, but rather units of e*E, e*B, e*A.
1978 ! (force, force, and energy, respectively). The reason is that the units of E, B or A
1979 ! are ugly.
1980 do il = 1, lasers%no_lasers
1981 select case (laser_kind(lasers%lasers(il)))
1983 aux = '[' // trim(units_abbrev(units_out%force)) // ']'
1984 do idir = 1, space%dim
1985 call write_iter_header(out_laser, aux)
1986 end do
1988 aux = '[' // trim(units_abbrev(units_out%energy)) // ']'
1989 do idir = 1, space%dim
1990 call write_iter_header(out_laser, aux)
1991 end do
1993 aux = '[adim]'
1994 call write_iter_header(out_laser, aux)
1995 end select
1996 end do
1997
1998 if (lasers_with_nondipole_field(lasers)) then
1999 aux = '[' // trim(units_abbrev(units_out%energy)) // ']'
2000 do idir = 1, space%dim
2001 call write_iter_header(out_laser, aux)
2002 end do
2003 end if
2004
2005
2006 call write_iter_nl(out_laser)
2007
2008 call td_write_print_header_end(out_laser)
2009 end if
2011 call write_iter_start(out_laser)
2012
2013 do il = 1, lasers%no_lasers
2014 field = m_zero
2015 call laser_field(lasers%lasers(il), field(1:space%dim), iter*dt)
2016 select case (laser_kind(lasers%lasers(il)))
2018 field = units_from_atomic(units_out%force, field)
2019 call write_iter_double(out_laser, field, space%dim)
2021 field = units_from_atomic(units_out%energy, field)
2022 call write_iter_double(out_laser, field, space%dim)
2024 call write_iter_double(out_laser, field(1), 1)
2025 end select
2026 end do
2027
2028 if (lasers_with_nondipole_field(lasers)) then
2029 call lasers_nondipole_laser_field_step(lasers, ndfield, iter*dt)
2030 call lasers_set_nondipole_parameters(lasers, ndfield, iter*dt)
2031 call write_iter_double(out_laser, ndfield, space%dim)
2032 end if
2033 call write_iter_nl(out_laser)
2034
2035 end subroutine td_write_laser
2036
2037
2038 ! ---------------------------------------------------------
2039 subroutine td_write_energy(out_energy, hm, iter, ke)
2040 type(c_ptr), intent(inout) :: out_energy
2041 type(hamiltonian_elec_t), intent(in) :: hm
2042 integer, intent(in) :: iter
2043 real(real64), intent(in) :: ke
2044
2045 integer :: ii
2046
2047 integer :: n_columns
2048
2049 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
2050
2051 push_sub(td_write_energy)
2052
2053 n_columns = 9
2054
2055 if (iter == 0) then
2056 call td_write_print_header_init(out_energy)
2057
2058 ! first line -> column names
2059 call write_iter_header_start(out_energy)
2060 call write_iter_header(out_energy, 'Total')
2061 call write_iter_header(out_energy, 'Kinetic (ions)')
2062 call write_iter_header(out_energy, 'Ion-Ion')
2063 call write_iter_header(out_energy, 'Electronic')
2064 call write_iter_header(out_energy, 'Eigenvalues')
2065 call write_iter_header(out_energy, 'Hartree')
2066 call write_iter_header(out_energy, 'Int[n v_xc]')
2067 call write_iter_header(out_energy, 'Exchange')
2068 call write_iter_header(out_energy, 'Correlation')
2069
2070 if (hm%pcm%run_pcm) then
2071 call write_iter_header(out_energy, 'E_M-solvent')
2072 n_columns = n_columns + 1
2073 end if
2074
2075 if (hm%lda_u_level /= dft_u_none) then
2076 call write_iter_header(out_energy, 'Hubbard')
2077 n_columns = n_columns + 1
2078 end if
2079
2080 call write_iter_nl(out_energy)
2081
2082 ! units
2083
2084 call write_iter_string(out_energy, '#[Iter n.]')
2085 call write_iter_header(out_energy, '[' // trim(units_abbrev(units_out%time)) // ']')
2086
2087 do ii = 1, n_columns
2088 call write_iter_header(out_energy, '[' // trim(units_abbrev(units_out%energy)) // ']')
2089 end do
2090 call write_iter_nl(out_energy)
2091
2092
2093 call td_write_print_header_end(out_energy)
2094 end if
2095
2096 call write_iter_start(out_energy)
2097 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%total+ke), 1)
2098 call write_iter_double(out_energy, units_from_atomic(units_out%energy, ke), 1)
2099 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%ep%eii), 1)
2100 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%total-hm%ep%eii), 1)
2101 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%eigenvalues), 1)
2102 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%hartree), 1)
2103 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%intnvxc), 1)
2104 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%exchange), 1)
2105 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%correlation), 1)
2106
2107 !adding the molecule-solvent electrostatic interaction
2108 if (hm%pcm%run_pcm) call write_iter_double(out_energy, &
2109 units_from_atomic(units_out%energy, hm%energy%int_ee_pcm + hm%energy%int_en_pcm + &
2110 hm%energy%int_nn_pcm + hm%energy%int_ne_pcm), 1)
2111
2112 if (hm%lda_u_level /= dft_u_none) then
2113 call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%dft_u), 1)
2114 end if
2115
2116 call write_iter_nl(out_energy)
2117
2118 pop_sub(td_write_energy)
2119 end subroutine td_write_energy
2120
2121 ! ---------------------------------------------------------
2122 subroutine td_write_eigs(out_eigs, st, iter)
2123 type(c_ptr), intent(inout) :: out_eigs
2124 type(states_elec_t), intent(in) :: st
2125 integer, intent(in) :: iter
2126
2127 integer :: ii, is
2128 character(len=68) :: buf
2129
2130 push_sub(td_write_eigs)
2131
2133 pop_sub(td_write_eigs)
2134 return ! only first node outputs
2135 end if
2136
2137
2138 if (iter == 0) then
2139 call td_write_print_header_init(out_eigs)
2140
2141 write(buf, '(a15,i2)') '# nst ', st%nst
2142 call write_iter_string(out_eigs, buf)
2143 call write_iter_nl(out_eigs)
2144
2145 write(buf, '(a15,i2)') '# nspin ', st%d%nspin
2146 call write_iter_string(out_eigs, buf)
2147 call write_iter_nl(out_eigs)
2148
2149 ! first line -> column names
2150 call write_iter_header_start(out_eigs)
2151 do is = 1, st%d%kpt%nglobal
2152 do ii = 1, st%nst
2153 write(buf, '(a,i4)') 'Eigenvalue ',ii
2154 call write_iter_header(out_eigs, buf)
2155 end do
2156 end do
2157 call write_iter_nl(out_eigs)
2158
2159 ! second line: units
2160 call write_iter_string(out_eigs, '#[Iter n.]')
2161 call write_iter_header(out_eigs, '[' // trim(units_abbrev(units_out%time)) // ']')
2162 do is = 1, st%d%kpt%nglobal
2163 do ii = 1, st%nst
2164 call write_iter_header(out_eigs, '[' // trim(units_abbrev(units_out%energy)) // ']')
2165 end do
2166 end do
2167 call write_iter_nl(out_eigs)
2168 call td_write_print_header_end(out_eigs)
2169 end if
2170
2171 call write_iter_start(out_eigs)
2172 do is = 1, st%d%kpt%nglobal
2173 do ii =1 , st%nst
2174 call write_iter_double(out_eigs, units_from_atomic(units_out%energy, st%eigenval(ii,is)), 1)
2175 end do
2176 end do
2177 call write_iter_nl(out_eigs)
2178
2179 pop_sub(td_write_eigs)
2180 end subroutine td_write_eigs
2181
2182 ! ---------------------------------------------------------
2183 subroutine td_write_ionch(out_ionch, mesh, st, iter)
2184 type(c_ptr), intent(inout) :: out_ionch
2185 class(mesh_t), intent(in) :: mesh
2186 type(states_elec_t), intent(in) :: st
2187 integer, intent(in) :: iter
2188
2189 integer :: ii, ist, Nch, ik, idim
2190 character(len=68) :: buf
2191 real(real64), allocatable :: ch(:), occ(:)
2192 real(real64), allocatable :: occbuf(:)
2193
2194 push_sub(td_write_ionch)
2195
2196
2197 nch = st%nst * st%d%kpt%nglobal * st%d%dim
2198 safe_allocate(ch(0: nch))
2199 safe_allocate(occ(0: nch))
2200
2201 occ(:) = m_zero
2202 ii = 1
2203 do ik = 1, st%nik
2204 do ist = 1, st%nst
2205 do idim = 1, st%d%dim
2206 if (st%st_start <= ist .and. ist <= st%st_end .and. &
2207 st%d%kpt%start <= ik .and. ik <= st%d%kpt%end) then
2208 occ(ii) = st%occ(ist, ik)
2209 end if
2210 ii = ii+1
2211 end do
2212 end do
2213 end do
2214
2216 if (st%parallel_in_states) then
2217 safe_allocate(occbuf(0: nch))
2218 occbuf(:) = m_zero
2219 call st%mpi_grp%allreduce(occ(0), occbuf(0), nch+1, mpi_double_precision, mpi_sum)
2220 occ(:) = occbuf(:)
2221 safe_deallocate_a(occbuf)
2222 end if
2223
2224 !Calculate the channels
2225 call td_calc_ionch(mesh, st, ch, nch)
2226
2227
2228 if (.not. mpi_grp_is_root(mpi_world)) then
2229 safe_deallocate_a(ch)
2230 pop_sub(td_write_ionch)
2231 return ! only first node outputs
2232 end if
2233
2234
2235 if (iter == 0) then
2236 call td_write_print_header_init(out_ionch)
2237
2238 ! first line -> column names
2239 call write_iter_header_start(out_ionch)
2240
2241 do ii = 0, nch
2242 if (occ(ii)>m_zero .or. ii == 0) then
2243 write(buf, '(a,f4.1,a)') 'Pion(',occ(ii)*ii,'+, t)'
2244 call write_iter_header(out_ionch, buf)
2245 end if
2246 end do
2247 call write_iter_nl(out_ionch)
2248
2249 ! second line: units
2250 call write_iter_string(out_ionch, '#[Iter n.]')
2251 call write_iter_header(out_ionch, '[' // trim(units_abbrev(units_out%time)) // ']')
2252 do ii = 0, nch
2253 if (occ(ii)>m_zero .or. ii == 0) then
2254 call write_iter_header(out_ionch, '[' // trim(units_abbrev(unit_one)) // ']')
2255 end if
2256 end do
2257 call write_iter_nl(out_ionch)
2258 call td_write_print_header_end(out_ionch)
2259 end if
2260
2261 call write_iter_start(out_ionch)
2262 do ii =0 , nch
2263 if (occ(ii)>m_zero .or. ii == 0) then
2264 call write_iter_double(out_ionch, units_from_atomic(unit_one, ch(ii)), 1)
2265 end if
2266 end do
2267 call write_iter_nl(out_ionch)
2268
2269 safe_deallocate_a(ch)
2270 safe_deallocate_a(occ)
2271
2272 pop_sub(td_write_ionch)
2273 end subroutine td_write_ionch
2274
2275 ! ---------------------------------------------------------
2276 subroutine td_write_proj(out_proj, space, mesh, ions, st, gs_st, kick, iter)
2277 type(c_ptr), intent(inout) :: out_proj
2278 class(space_t), intent(in) :: space
2279 class(mesh_t), intent(in) :: mesh
2280 type(ions_t), intent(in) :: ions
2281 type(states_elec_t), intent(inout) :: st
2282 type(states_elec_t), intent(in) :: gs_st
2283 type(kick_t), intent(in) :: kick
2284 integer, intent(in) :: iter
2285
2286 complex(real64), allocatable :: projections(:,:,:)
2287 character(len=80) :: aux
2288 integer :: ik, ist, uist, idir
2289
2290 push_sub(td_write_proj)
2291
2292 if (iter == 0) then
2293 if (mpi_grp_is_root(mpi_world)) then
2294 call td_write_print_header_init(out_proj)
2295
2296 write(aux, '(a15,i2)') '# nspin ', st%d%nspin
2297 call write_iter_string(out_proj, aux)
2298 call write_iter_nl(out_proj)
2299
2300 call kick_write(kick, out = out_proj)
2301
2302 call write_iter_string(out_proj, "#%")
2303 call write_iter_nl(out_proj)
2304
2305 write(aux, '(a,i8)') "# nik ", st%nik
2306 call write_iter_string(out_proj, aux)
2307 call write_iter_nl(out_proj)
2308
2309 write(aux, '(a,2i8)') "# st ", gs_st%st_start, st%nst
2310 call write_iter_string(out_proj, aux)
2311 call write_iter_nl(out_proj)
2312
2313 write(aux, '(a,2i8)') "# ust ", gs_st%st_start, gs_st%st_end
2314 call write_iter_string(out_proj, aux)
2315 call write_iter_nl(out_proj)
2316
2317 do ik = 1, st%nik
2318 call write_iter_string(out_proj, "# w(ik)*occ(ist,ik) ")
2319 do ist = gs_st%st_start, st%nst
2320 call write_iter_double(out_proj, st%kweights(ik)*st%occ(ist, ik), 1)
2321 end do
2322 call write_iter_nl(out_proj)
2323 end do
2324
2325 call write_iter_header_start(out_proj)
2326 do ik = 1, st%nik
2327 do ist = gs_st%st_start, st%nst
2328 do uist = gs_st%st_start, gs_st%st_end
2329 write(aux, '(i4,a,i4)') ist, ' -> ', uist
2330 call write_iter_header(out_proj, 'Re {'//trim(aux)//'}')
2331 call write_iter_header(out_proj, 'Im {'//trim(aux)//'}')
2332 end do
2333 end do
2334 end do
2335 call write_iter_nl(out_proj)
2336
2337 end if
2338
2339 !The dipole matrix elements cannot be computed like that for solids
2340 if (.not. space%is_periodic()) then
2341
2342 safe_allocate(projections(1:st%nst, gs_st%st_start:gs_st%st_end, 1:st%nik))
2343 do idir = 1, space%dim
2344 projections = m_z0
2345
2346 call dipole_matrix_elements(idir)
2347
2348 if (mpi_grp_is_root(mpi_world)) then
2349 write(aux, '(a,i1,a)') "<i|x_", idir, "|a>"
2350 call write_iter_string(out_proj, "# ------")
2351 call write_iter_header(out_proj, aux)
2352 do ik = 1, st%nik
2353 do ist = gs_st%st_start, st%st_end
2354 do uist = gs_st%st_start, gs_st%st_end
2355 call write_iter_double(out_proj, real(projections(ist, uist, ik), real64), 1)
2356 call write_iter_double(out_proj, aimag(projections(ist, uist, ik)), 1)
2357 end do
2358 end do
2359 end do
2360 call write_iter_nl(out_proj)
2361
2362 end if
2363 end do
2364 safe_deallocate_a(projections)
2365
2366 end if
2367
2368 if (mpi_grp_is_root(mpi_world)) then
2370 end if
2371
2372 end if
2373
2374 safe_allocate(projections(1:st%nst, gs_st%st_start:gs_st%st_end, 1:st%nik))
2375 projections(:,:,:) = m_z0
2376 call calc_projections(mesh, st, gs_st, projections)
2377
2378 if (mpi_grp_is_root(mpi_world)) then
2379 call write_iter_start(out_proj)
2380 do ik = 1, st%nik
2381 do ist = gs_st%st_start, st%nst
2382 do uist = gs_st%st_start, gs_st%st_end
2383 call write_iter_double(out_proj, real(projections(ist, uist, ik), real64), 1)
2384 call write_iter_double(out_proj, aimag(projections(ist, uist, ik)), 1)
2385 end do
2386 end do
2387 end do
2388 call write_iter_nl(out_proj)
2389 end if
2390
2391 safe_deallocate_a(projections)
2392 pop_sub(td_write_proj)
2393
2394 contains
2395 ! ---------------------------------------------------------
2396 subroutine dipole_matrix_elements(dir)
2397 integer, intent(in) :: dir
2398
2399 integer :: uist, ist, ik, idim
2400 real(real64) :: n_dip(space%dim)
2401 complex(real64), allocatable :: xpsi(:,:)
2402 complex(real64), allocatable :: psi(:, :), gspsi(:, :)
2403
2405
2406 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
2407 safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
2408 safe_allocate(xpsi(1:mesh%np, 1:st%d%dim))
2409
2410 do ik = st%d%kpt%start, st%d%kpt%end
2411 do ist = st%st_start, st%st_end
2412 call states_elec_get_state(st, mesh, ist, ik, psi)
2413 do uist = gs_st%st_start, gs_st%st_end
2414 call states_elec_get_state(gs_st, mesh, uist, ik, gspsi)
2415
2416 do idim = 1, st%d%dim
2417 xpsi(1:mesh%np, idim) = mesh%x(1:mesh%np, dir)*gspsi(1:mesh%np, idim)
2418 end do
2419 projections(ist, uist, ik) = -zmf_dotp(mesh, st%d%dim, psi, xpsi, reduce = .false.)
2420
2421 end do
2422 end do
2423 end do
2424
2425 safe_deallocate_a(xpsi)
2426 safe_deallocate_a(gspsi)
2427 safe_deallocate_a(psi)
2428
2429 call comm_allreduce(st%dom_st_kpt_mpi_grp, projections)
2430
2431 ! n_dip is not defined for more than space%dim
2432 n_dip = ions%dipole()
2433 do ik = 1, st%nik
2434 do ist = gs_st%st_start, st%nst
2435 do uist = gs_st%st_start, gs_st%st_end
2436 projections(ist, uist, ik) = projections(ist, uist, ik) - n_dip(dir)
2437 end do
2438 end do
2439 end do
2440
2441
2443 end subroutine dipole_matrix_elements
2444
2445 end subroutine td_write_proj
2446
2447 ! ---------------------------------------------------------
2451 ! ---------------------------------------------------------
2452 subroutine td_write_n_ex(out_nex, outp, namespace, mesh, kpoints, st, gs_st, iter)
2453 type(c_ptr), intent(inout) :: out_nex
2454 type(output_t), intent(in) :: outp
2455 type(namespace_t), intent(in) :: namespace
2456 class(mesh_t), intent(in) :: mesh
2457 type(kpoints_t), intent(in) :: kpoints
2458 type(states_elec_t), intent(inout) :: st
2459 type(states_elec_t), intent(inout) :: gs_st
2460 integer, intent(in) :: iter
2461
2462 complex(real64), allocatable :: projections(:,:)
2463 character(len=80) :: aux, dir
2464 integer :: ik, ikpt, ist, uist, err
2465 real(real64) :: Nex, weight
2466 integer :: gs_nst
2467 real(real64), allocatable :: Nex_kpt(:)
2468
2469
2470 push_sub(td_write_n_ex)
2471
2472 if (iter == 0) then
2473 if (mpi_grp_is_root(mpi_world)) then
2474 call td_write_print_header_init(out_nex)
2475
2476 write(aux, '(a15,i2)') '# nspin ', st%d%nspin
2477 call write_iter_string(out_nex, aux)
2478 call write_iter_nl(out_nex)
2479
2480 call write_iter_string(out_nex, "#%")
2481 call write_iter_nl(out_nex)
2482
2483 write(aux, '(a,i8)') "# nik ", st%nik
2484 call write_iter_string(out_nex, aux)
2485 call write_iter_nl(out_nex)
2486
2487 write(aux, '(a,2i8)') "# st ", gs_st%st_start, st%nst
2488 call write_iter_string(out_nex, aux)
2489 call write_iter_nl(out_nex)
2490
2491 write(aux, '(a,2i8)') "# ust ", gs_st%st_start, gs_st%st_end
2492 call write_iter_string(out_nex, aux)
2493 call write_iter_nl(out_nex)
2494
2495 call write_iter_header_start(out_nex)
2496 call write_iter_header(out_nex, '# iter t Nex(t)')
2497 call write_iter_nl(out_nex)
2498
2499 end if
2500
2501 if (mpi_grp_is_root(mpi_world)) then
2502 call td_write_print_header_end(out_nex)
2503 end if
2504
2505 end if
2506
2507 !We only need the occupied GS states
2508 if (.not. st%parallel_in_states) then
2509 do ik = st%d%kpt%start, st%d%kpt%end
2510 do ist = 1, gs_st%nst
2511 if (gs_st%occ(ist, ik) > m_min_occ) gs_nst = ist
2512 end do
2513 end do
2514 else
2515 gs_nst = gs_st%nst
2516 end if
2517
2518 safe_allocate(projections(1:gs_nst, 1:st%nst))
2519
2520 safe_allocate(nex_kpt(1:st%nik))
2521 nex_kpt = m_zero
2522 do ik = st%d%kpt%start, st%d%kpt%end
2523 ikpt = st%d%get_kpoint_index(ik)
2524 call zstates_elec_calc_projections(st, gs_st, namespace, mesh, ik, projections, gs_nst)
2525 do ist = 1, gs_nst
2526 weight = st%kweights(ik) * gs_st%occ(ist, ik)/ st%smear%el_per_state
2527 do uist = st%st_start, st%st_end
2528 nex_kpt(ikpt) = nex_kpt(ikpt) - weight * st%occ(uist, ik) * abs(projections(ist, uist))**2
2529 end do
2530 end do
2531 if (st%d%ispin == spin_polarized) then
2532 nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*m_half*st%kweights(ik)
2533 else
2534 nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*st%kweights(ik)
2535 end if
2536 end do
2537
2538 if (st%parallel_in_states .or. st%d%kpt%parallel) then
2539 call comm_allreduce(st%st_kpt_mpi_grp, nex_kpt)
2540 end if
2541
2542 nex = sum(nex_kpt)
2543
2544 if (mpi_grp_is_root(mpi_world)) then
2545 call write_iter_start(out_nex)
2546 call write_iter_double(out_nex, nex, 1)
2547 call write_iter_nl(out_nex)
2548
2549 ! now write down the k-resolved part
2550 write(dir, '(a,a,i7.7)') trim(outp%iter_dir),"td.", iter ! name of directory
2551
2552 call io_function_output_global_bz(outp%how(option__output__current_kpt) &
2553 + outp%how(option__output__density_kpt), dir, "n_excited_el_kpt", namespace, &
2554 kpoints, nex_kpt, unit_one, err)
2555 end if
2556
2557 safe_deallocate_a(projections)
2558 safe_deallocate_a(nex_kpt)
2559
2560 pop_sub(td_write_n_ex)
2561 end subroutine td_write_n_ex
2562
2563 ! ---------------------------------------------------------
2568 ! ---------------------------------------------------------
2569 subroutine calc_projections(mesh, st, gs_st, projections)
2570 class(mesh_t), intent(in) :: mesh
2571 type(states_elec_t), intent(inout) :: st
2572 type(states_elec_t), intent(in) :: gs_st
2573 complex(real64), intent(inout) :: projections(1:st%nst, gs_st%st_start:gs_st%nst, 1:st%nik)
2574
2575 integer :: uist, ist, ik
2576 complex(real64), allocatable :: psi(:, :), gspsi(:, :)
2577 push_sub(calc_projections)
2578
2579 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
2580 safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
2581
2582 projections(:,:,:) = m_zero
2583
2584 do ik = st%d%kpt%start, st%d%kpt%end
2585 do ist = st%st_start, st%st_end
2586 call states_elec_get_state(st, mesh, ist, ik, psi)
2587 do uist = gs_st%st_start, gs_st%nst
2588 call states_elec_get_state(gs_st, mesh, uist, ik, gspsi)
2589 projections(ist, uist, ik) = zmf_dotp(mesh, st%d%dim, psi, gspsi, reduce = .false.)
2590 end do
2591 end do
2592 end do
2593
2594 safe_deallocate_a(psi)
2595 safe_deallocate_a(gspsi)
2596
2597 call comm_allreduce(st%dom_st_kpt_mpi_grp, projections)
2598
2599 pop_sub(calc_projections)
2600 end subroutine calc_projections
2601
2602
2603 subroutine td_write_proj_kp(mesh, kpoints, st, gs_st, namespace, iter)
2604 class(mesh_t), intent(in) :: mesh
2605 type(kpoints_t), intent(in) :: kpoints
2606 type(states_elec_t), intent(in) :: st
2607 type(states_elec_t), intent(inout) :: gs_st
2608 type(namespace_t), intent(in) :: namespace
2609 integer, intent(in) :: iter
2610
2611 complex(real64), allocatable :: proj(:,:), psi(:,:,:), gs_psi(:,:,:), temp_state(:,:)
2612 character(len=80) :: filename1, filename2
2613 integer :: ik,ist, jst, file, idim, nk_proj
2614
2615 push_sub(td_write_proj_kp)
2616
2617 write(filename1,'(I10)') iter
2618 filename1 = 'td.general/projections_iter_'//trim(adjustl(filename1))
2619 file = 9845623
2620
2621 safe_allocate(proj(1:gs_st%nst, 1:gs_st%nst))
2622 safe_allocate(psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2623 safe_allocate(gs_psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2624 safe_allocate(temp_state(1:mesh%np,1:gs_st%d%dim))
2625
2626 ! Project only k-points that have a zero weight.
2627 ! Why? It is unlikely that one is interested in the projections
2628 ! of the Monkhorst-Pack kpoints, but instead we assume that
2629 ! the user has specified a k-path with zero weights
2630 nk_proj = kpoints%nik_skip
2631
2632 do ik = kpoints%reduced%npoints-nk_proj+1, kpoints%reduced%npoints
2633 ! reset arrays
2634 psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)= m_zero
2635 gs_psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)= m_zero
2636 ! open file for writing
2637 if (mpi_world%rank == 0) then
2638 write(filename2,'(I10)') ik
2639 filename2 = trim(adjustl(filename1))//'_ik_'//trim(adjustl(filename2))
2640 file = io_open(filename2, namespace, action='write')
2641 end if
2642 ! get all states at ik that are locally stored (ground state and td-states)
2643 do ist=gs_st%st_start,gs_st%st_end
2644 if (state_kpt_is_local(gs_st, ist, ik)) then
2645 call states_elec_get_state(st, mesh, ist, ik,temp_state)
2646 do idim = 1,gs_st%d%dim
2647 psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2648 end do
2649 call states_elec_get_state(gs_st, mesh, ist, ik, temp_state)
2650 do idim = 1,gs_st%d%dim
2651 gs_psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2652 end do
2653 end if
2654 end do
2655 ! collect states at ik from all processes in one array
2656 call comm_allreduce(mpi_world, psi)
2657 call comm_allreduce(mpi_world, gs_psi)
2658
2659 ! compute the overlaps as a matrix product
2660 assert(mesh%np_global*gs_st%d%dim < huge(0_int32))
2661 proj(1:gs_st%nst, 1:gs_st%nst) = m_zero
2662 call zgemm('n', &
2663 'c', &
2664 gs_st%nst, &
2665 gs_st%nst, &
2666 i8_to_i4(mesh%np_global*gs_st%d%dim), &
2667 cmplx(mesh%volume_element, m_zero, real64) , &
2668 psi(1, 1, 1), &
2669 ubound(psi, dim = 1), &
2670 gs_psi(1, 1, 1), &
2671 ubound(gs_psi, dim = 1), &
2672 m_z0, &
2673 proj(1, 1), &
2674 ubound(proj, dim = 1))
2675
2676 ! write to file
2677 if (mpi_world%rank == 0) then
2678 do ist = 1, gs_st%nst
2679 do jst = 1, gs_st%nst
2680 write(file,'(I3,1x,I3,1x,e13.6,1x,e13.6,2x)') ist, jst, proj(ist,jst)
2681 end do
2682 end do
2683 call io_close(file)
2684 end if
2685
2686 end do! ik
2687
2688 safe_deallocate_a(proj)
2689 safe_deallocate_a(psi)
2690 safe_deallocate_a(gs_psi)
2691 safe_deallocate_a(temp_state)
2692
2693 pop_sub(td_write_proj_kp)
2694 end subroutine td_write_proj_kp
2695
2696 !---------------------------------------
2697 subroutine td_write_floquet(namespace, space, hm, ext_partners, mesh, st, iter)
2698 type(namespace_t), intent(in) :: namespace
2699 class(space_t), intent(in) :: space
2700 type(hamiltonian_elec_t), intent(inout) :: hm
2701 type(partner_list_t), intent(in) :: ext_partners
2702 class(mesh_t), intent(in) :: mesh
2703 type(states_elec_t), intent(inout) :: st
2704 integer, intent(in) :: iter
2705
2706 complex(real64), allocatable :: hmss(:,:), psi(:,:,:), hpsi(:,:,:), temp_state1(:,:)
2707 complex(real64), allocatable :: HFloquet(:,:,:), HFloq_eff(:,:), temp(:,:)
2708 real(real64), allocatable :: eigenval(:), bands(:,:)
2709 character(len=80) :: filename
2710 integer :: it, nT, ik, ist, in, im, file, idim, nik, ik_count
2711 integer :: Forder, Fdim, m0, n0, n1, nst, ii, jj, lim_nst
2712 logical :: downfolding
2713 type(states_elec_t) :: hm_st
2714
2715 real(real64) :: dt, Tcycle, omega
2716
2717 push_sub(td_write_floquet)
2718
2719 downfolding = .false.
2720
2721 ! this does not depend on propagation, so we do it only once
2722 if (.not. iter == 0) then
2723 pop_sub(td_write_floquet)
2724 return
2725 end if
2726
2727 nst = st%nst
2728
2729 !for now no domain distributionallowed
2730 assert(mesh%np == mesh%np_global)
2731
2732 ! this is used to initialize the hpsi (more effiecient ways?)
2733 call states_elec_copy(hm_st, st)
2734
2735 !%Variable TDFloquetFrequency
2736 !%Type float
2737 !%Default 0
2738 !%Section Time-Dependent::TD Output
2739 !%Description
2740 !% Frequency for the Floquet analysis, this should be the carrier frequency or integer multiples of it.
2741 !% Other options will work, but likely be nonsense.
2742 !%
2743 !%End
2744 call parse_variable(namespace, 'TDFloquetFrequency', m_zero, omega, units_inp%energy)
2745 call messages_print_var_value('Frequency used for Floquet analysis', omega, namespace=namespace)
2746 if (abs(omega) <= m_epsilon) then
2747 message(1) = "Please give a non-zero value for TDFloquetFrequency"
2748 call messages_fatal(1, namespace=namespace)
2749 end if
2750
2751 ! get time of one cycle
2752 tcycle = m_two * m_pi / omega
2753
2754 !%Variable TDFloquetSample
2755 !%Type integer
2756 !%Default 20
2757 !%Section Time-Dependent::TD Output
2758 !%Description
2759 !% Number of points on which one Floquet cycle is sampled in the time-integral of the Floquet analysis.
2760 !%
2761 !%End
2762 call parse_variable(namespace, 'TDFloquetSample',20 ,nt)
2763 call messages_print_var_value('Number of Floquet time-sampling points', nt, namespace=namespace)
2764 dt = tcycle/real(nt, real64)
2765
2766 !%Variable TDFloquetDimension
2767 !%Type integer
2768 !%Default -1
2769 !%Section Time-Dependent::TD Output
2770 !%Description
2771 !% Order of Floquet Hamiltonian. If negative number is given, downfolding is performed.
2772 !%End
2773 call parse_variable(namespace, 'TDFloquetDimension',-1,forder)
2774 if (forder .ge. 0) then
2775 call messages_print_var_value('Order of multiphoton Floquet-Hamiltonian', forder, namespace=namespace)
2776 !Dimension of multiphoton Floquet-Hamiltonian
2777 fdim = 2 * forder + 1
2778 else
2779 message(1) = 'Floquet-Hamiltonian is downfolded'
2780 call messages_info(1, namespace=namespace)
2781 downfolding = .true.
2782 forder = 1
2783 fdim = 3
2784 end if
2785
2786 dt = tcycle/real(nt, real64)
2787
2788 ! we are only interested for k-point with zero weight
2789 nik = hm%kpoints%nik_skip
2791 safe_allocate(hmss(1:nst,1:nst))
2792 safe_allocate( psi(1:nst,1:st%d%dim,1:mesh%np))
2793 safe_allocate(hpsi(1:nst,1:st%d%dim,1:mesh%np))
2794 safe_allocate(temp_state1(1:mesh%np,1:st%d%dim))
2795
2796 ! multiphoton Floquet Hamiltonian, layout:
2797 ! (H_{-n,-m} ... H_{-n,0} ... H_{-n,m})
2798 ! ( . . . . . )
2799 ! H = (H_{0,-m} ... H_{0,0} ... H_{0,m} )
2800 ! ( . . . . . )
2801 ! (H_{n,-m} ... H_{n,0} ... H_{n,m} )
2802 safe_allocate(hfloquet(1:nik,1:nst*fdim, 1:nst*fdim))
2803 hfloquet(1:nik,1:nst*fdim, 1:nst*fdim) = m_zero
2804
2805 ! perform time-integral over one cycle
2806 do it = 1, nt
2807 ! get non-interacting Hamiltonian at time (offset by one cycle to allow for ramp)
2808 call hm%update(mesh, namespace, space, ext_partners, time=tcycle+it*dt)
2809 ! get hpsi
2810 call zhamiltonian_elec_apply_all(hm, namespace, mesh, st, hm_st)
2811
2812 ! project Hamiltonian into grounstates for zero weight k-points
2813 ik_count = 0
2814
2815 do ik = hm%kpoints%reduced%npoints-nik+1, hm%kpoints%reduced%npoints
2816 ik_count = ik_count + 1
2817
2818 psi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
2819 hpsi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
2820
2821 do ist = st%st_start, st%st_end
2822 if (state_kpt_is_local(st, ist, ik)) then
2823 call states_elec_get_state(st, mesh, ist, ik,temp_state1)
2824 do idim = 1, st%d%dim
2825 psi(ist, idim, 1:mesh%np) = temp_state1(1:mesh%np,idim)
2826 end do
2827 call states_elec_get_state(hm_st, mesh, ist, ik,temp_state1)
2828 do idim = 1, st%d%dim
2829 hpsi(ist, idim,1:mesh%np) = temp_state1(1:mesh%np,idim)
2830 end do
2831 end if
2832 end do
2833 call comm_allreduce(mpi_world, psi)
2834 call comm_allreduce(mpi_world, hpsi)
2835 assert(mesh%np_global*st%d%dim < huge(0_int32))
2836 hmss(1:nst,1:nst) = m_zero
2837 call zgemm( 'n', &
2838 'c', &
2839 nst, &
2840 nst, &
2841 i8_to_i4(mesh%np_global*st%d%dim), &
2842 cmplx(mesh%volume_element, m_zero, real64) , &
2843 hpsi(1, 1, 1), &
2844 ubound(hpsi, dim = 1), &
2845 psi(1, 1, 1), &
2846 ubound(psi, dim = 1), &
2847 m_z0, &
2848 hmss(1, 1), &
2849 ubound(hmss, dim = 1))
2850
2851 hmss(1:nst,1:nst) = conjg(hmss(1:nst,1:nst))
2852
2853 ! accumulate the Floqeut integrals
2854 do in = -forder, forder
2855 do im = -forder, forder
2856 ii = (in+forder) * nst
2857 jj = (im+forder) * nst
2858 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) = &
2859 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) + hmss(1:nst, 1:nst) * exp(-(in-im)*m_zi*omega*it*dt)
2860 ! diagonal term
2861 if (in == im) then
2862 do ist = 1, nst
2863 hfloquet(ik_count, ii+ist, ii+ist) = hfloquet(ik_count, ii+ist, ii+ist) + in*omega
2864 end do
2865 end if
2866 end do
2867 end do
2868 end do !ik
2869
2870 end do ! it
2871
2872 hfloquet(:,:,:) = m_one/nt*hfloquet(:,:,:)
2873
2874 ! diagonalize Floquet Hamiltonian
2875 if (downfolding) then
2876 ! here perform downfolding
2877 safe_allocate(hfloq_eff(1:nst,1:nst))
2878 safe_allocate(eigenval(1:nst))
2879 safe_allocate(bands(1:nik,1:nst))
2880
2881 hfloq_eff(1:nst,1:nst) = m_zero
2882 do ik = 1, nik
2883 ! the HFloquet blocks are copied directly out of the super matrix
2884 m0 = nst ! the m=0 start position
2885 n0 = nst ! the n=0 start postion
2886 n1 = 2*nst ! the n=+1 start postion
2887 hfloq_eff(1:nst, 1:nst) = hfloquet(ik, n0+1:n0+nst, m0+1:m0+nst) + &
2888 m_one/omega*(matmul(hfloquet(ik, 1:nst, m0+1:m0+nst), hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst))- &
2889 matmul(hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst), hfloquet(ik, 1:nst, m0+1:m0+nst)))
2890
2891 call lalg_eigensolve(nst, hfloq_eff, eigenval)
2892 bands(ik,1:nst) = eigenval(1:nst)
2893 end do
2894 safe_deallocate_a(hfloq_eff)
2895 else
2896 ! the full Floquet
2897 safe_allocate(eigenval(1:nst*fdim))
2898 safe_allocate(bands(1:nik,1:nst*fdim))
2899 safe_allocate(temp(1:nst*fdim, 1:nst*fdim))
2900
2901 do ik = 1, nik
2902 temp(1:nst*fdim, 1:nst*fdim) = hfloquet(ik, 1:nst*fdim, 1:nst*fdim)
2903 call lalg_eigensolve(nst*fdim, temp, eigenval)
2904 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2905 end do
2906 end if
2907
2908 !write bandstructure to file
2909 if (downfolding) then
2910 lim_nst = nst
2911 filename = "downfolded_floquet_bands"
2912 else
2913 lim_nst = nst*fdim
2914 filename = "floquet_bands"
2915 end if
2916 ! write bands (full or downfolded)
2917 if (mpi_world%rank == 0) then
2918 file = 987254
2919 file = io_open(filename, namespace, action = 'write')
2920 do ik = 1,nik
2921 do ist = 1,lim_nst
2922 write(file,'(e13.6, 1x)',advance='no') bands(ik,ist)
2923 end do
2924 write(file,'(1x)')
2925 end do
2926 call io_close(file)
2927 end if
2928
2929 if (.not. downfolding) then
2930 ! for the full Floquet case compute also the trivially shifted
2931 ! Floquet bands for reference (i.e. setting H_{nm}=0 for n!=m)
2932 bands(1:nik, 1:nst*fdim) = m_zero
2933 do ik = 1, nik
2934 temp(1:nst*fdim,1:nst*fdim) = m_zero
2935 do jj = 0, fdim - 1
2936 ii = jj * nst
2937 temp(ii+1:ii+nst, ii+1:ii+nst) = hfloquet(ik, ii+1:ii+nst, ii+1:ii+nst)
2938 end do
2939 call lalg_eigensolve(nst*fdim, temp, eigenval)
2940 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2941 end do
2942
2943 if (mpi_world%rank == 0) then
2944 filename = 'trivial_floquet_bands'
2945 file = io_open(filename, namespace, action = 'write')
2946 do ik=1,nik
2947 do ist = 1,lim_nst
2948 write(file,'(e13.6, 1x)', advance='no') bands(ik,ist)
2949 end do
2950 write(file,'(1x)')
2951 end do
2952 call io_close(file)
2953 end if
2954 end if
2955
2956 ! reset time in Hamiltonian
2957 call hm%update(mesh, namespace, space, ext_partners, time=m_zero)
2958
2959 safe_deallocate_a(hmss)
2960 safe_deallocate_a(psi)
2961 safe_deallocate_a(hpsi)
2962 safe_deallocate_a(temp_state1)
2963 safe_deallocate_a(hfloquet)
2964 safe_deallocate_a(eigenval)
2965 safe_deallocate_a(bands)
2966 safe_deallocate_a(temp)
2967
2968 pop_sub(td_write_floquet)
2969
2970 end subroutine td_write_floquet
2971
2972 ! ---------------------------------------------------------
2973 subroutine td_write_total_current(out_total_current, space, mesh, st, iter)
2974 type(c_ptr), intent(inout) :: out_total_current
2975 class(space_t), intent(in) :: space
2976 class(mesh_t), intent(in) :: mesh
2977 type(states_elec_t), intent(in) :: st
2978 integer, intent(in) :: iter
2979
2980 integer :: idir, ispin
2981 character(len=50) :: aux
2982 real(real64) :: total_current(space%dim), abs_current(space%dim)
2983
2984 push_sub(td_write_total_current)
2985
2986 if (mpi_grp_is_root(mpi_world) .and. iter == 0) then
2987 call td_write_print_header_init(out_total_current)
2988
2989 ! first line: column names
2990 call write_iter_header_start(out_total_current)
2991
2992 do idir = 1, space%dim
2993 write(aux, '(a2,a1,a1)') 'I(', index2axis(idir), ')'
2994 call write_iter_header(out_total_current, aux)
2995 end do
2996
2997 do idir = 1, space%dim
2998 write(aux, '(a10,a1,a1)') 'IntAbs(j)(', index2axis(idir), ')'
2999 call write_iter_header(out_total_current, aux)
3000 end do
3001
3002 do ispin = 1, st%d%nspin
3003 do idir = 1, space%dim
3004 write(aux, '(a4,i1,a1,a1,a1)') 'I-sp', ispin, '(', index2axis(idir), ')'
3005 call write_iter_header(out_total_current, aux)
3006 end do
3007 end do
3008
3009 call write_iter_nl(out_total_current)
3010
3011 call td_write_print_header_end(out_total_current)
3012 end if
3013
3014 assert(allocated(st%current))
3015
3016 if (mpi_grp_is_root(mpi_world)) then
3017 call write_iter_start(out_total_current)
3018 end if
3019
3020 total_current = 0.0_real64
3021 do idir = 1, space%dim
3022 do ispin = 1, st%d%spin_channels
3023 total_current(idir) = total_current(idir) + dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.)
3024 end do
3025 total_current(idir) = units_from_atomic(units_out%length/units_out%time, total_current(idir))
3026 end do
3027 if (mesh%parallel_in_domains) then
3028 call mesh%allreduce(total_current, dim = space%dim)
3029 end if
3030
3031 abs_current = 0.0_real64
3032 do idir = 1, space%dim
3033 do ispin = 1, st%d%spin_channels
3034 abs_current(idir) = abs_current(idir) + dmf_integrate(mesh, abs(st%current(:, idir, ispin)), reduce = .false.)
3035 end do
3036 abs_current(idir) = units_from_atomic(units_out%length/units_out%time, abs_current(idir))
3037 end do
3038 if (mesh%parallel_in_domains) then
3039 call mesh%allreduce(abs_current, dim = space%dim)
3040 end if
3041
3042 if (mpi_grp_is_root(mpi_world)) then
3043 call write_iter_double(out_total_current, total_current, space%dim)
3044 call write_iter_double(out_total_current, abs_current, space%dim)
3045 end if
3046
3047 do ispin = 1, st%d%nspin
3048 total_current = 0.0_real64
3049 do idir = 1, space%dim
3050 total_current(idir) = units_from_atomic(units_out%length/units_out%time, &
3051 dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.))
3052 end do
3053 if (mesh%parallel_in_domains) then
3054 call mesh%allreduce(total_current, dim = space%dim)
3055 end if
3056 if (mpi_grp_is_root(mpi_world)) then
3057 call write_iter_double(out_total_current, total_current, space%dim)
3058 end if
3059 end do
3060
3061 if (mpi_grp_is_root(mpi_world)) then
3062 call write_iter_nl(out_total_current)
3063 end if
3064
3065 pop_sub(td_write_total_current)
3067
3068 ! ---------------------------------------------------------
3069
3070 subroutine td_write_total_heat_current(write_obj, space, hm, gr, st, iter)
3071 type(c_ptr), intent(inout) :: write_obj
3072 class(space_t), intent(in) :: space
3073 type(hamiltonian_elec_t), intent(inout) :: hm
3074 type(grid_t), intent(in) :: gr
3075 type(states_elec_t), intent(in) :: st
3076 integer, intent(in) :: iter
3077
3078 integer :: idir, ispin
3079 character(len=50) :: aux
3080 real(real64), allocatable :: heat_current(:, :, :)
3081 real(real64) :: total_current(space%dim)
3082
3083 push_sub(td_write_total_current)
3084
3085 if (mpi_grp_is_root(mpi_world) .and. iter == 0) then
3086 call td_write_print_header_init(write_obj)
3087
3088 ! first line: column names
3089 call write_iter_header_start(write_obj)
3090
3091 do idir = 1, space%dim
3092 write(aux, '(a2,i1,a1)') 'Jh(', idir, ')'
3093 call write_iter_header(write_obj, aux)
3094 end do
3095
3096 call write_iter_nl(write_obj)
3097
3098 call td_write_print_header_end(write_obj)
3099 end if
3100
3101 safe_allocate(heat_current(1:gr%np, 1:space%dim, 1:st%d%nspin))
3102
3103 call current_heat_calculate(space, gr%der, hm, st, heat_current)
3104
3105 if (mpi_grp_is_root(mpi_world)) call write_iter_start(write_obj)
3106
3107 total_current = 0.0_real64
3108 do idir = 1, space%dim
3109 do ispin = 1, st%d%spin_channels
3110 total_current(idir) = total_current(idir) + dmf_integrate(gr, heat_current(:, idir, ispin))
3111 end do
3112 total_current(idir) = units_from_atomic(units_out%energy*units_out%length/units_out%time, total_current(idir))
3113 end do
3114
3115 safe_deallocate_a(heat_current)
3116
3117 if (mpi_grp_is_root(mpi_world)) call write_iter_double(write_obj, total_current, space%dim)
3118
3119 if (mpi_grp_is_root(mpi_world)) call write_iter_nl(write_obj)
3120
3121 pop_sub(td_write_total_current)
3122 end subroutine td_write_total_heat_current
3123
3124
3125 ! ---------------------------------------------------------
3126 subroutine td_write_partial_charges(out_partial_charges, mesh, st, ions, iter)
3127 type(c_ptr), intent(inout) :: out_partial_charges
3128 class(mesh_t), intent(in) :: mesh
3129 type(states_elec_t), intent(in) :: st
3130 type(ions_t), intent(in) :: ions
3131 integer, intent(in) :: iter
3132
3133 integer :: idir
3134 character(len=50) :: aux
3135 real(real64), allocatable :: hirshfeld_charges(:)
3136
3137 push_sub(td_write_partial_charges)
3138
3139 safe_allocate(hirshfeld_charges(1:ions%natoms))
3140
3141 call partial_charges_calculate(mesh, st, ions, hirshfeld_charges)
3142
3143 if (mpi_grp_is_root(mpi_world)) then
3144
3145 if (iter == 0) then
3146
3147 call td_write_print_header_init(out_partial_charges)
3148
3149 ! first line: column names
3150 call write_iter_header_start(out_partial_charges)
3151
3152 do idir = 1, ions%natoms
3153 write(aux, '(a13,i3,a1)') 'hirshfeld(atom=', idir, ')'
3154 call write_iter_header(out_partial_charges, aux)
3155 end do
3156
3157 call write_iter_nl(out_partial_charges)
3158
3159 call td_write_print_header_end(out_partial_charges)
3160 end if
3161
3162 call write_iter_start(out_partial_charges)
3164 call write_iter_double(out_partial_charges, hirshfeld_charges, ions%natoms)
3165
3166 call write_iter_nl(out_partial_charges)
3167 end if
3168
3169 safe_deallocate_a(hirshfeld_charges)
3170
3172 end subroutine td_write_partial_charges
3173
3174 ! ---------------------------------------------------------
3175 subroutine td_write_q(out_q, space, ks, iter)
3176 type(c_ptr), intent(inout) :: out_q
3177 class(space_t), intent(in) :: space
3178 type(v_ks_t), intent(in) :: ks
3179 integer, intent(in) :: iter
3180
3181 integer :: ii
3182 character(len=50) :: aux
3183
3184 push_sub(td_write_q)
3185
3186 if (mpi_grp_is_root(mpi_world)) then
3187 if (iter == 0) then
3188 call td_write_print_header_init(out_q)
3189 call write_iter_header_start(out_q)
3190 do ii = 1, ks%pt%nmodes
3191 write(aux, '(a1,i3,a3)') 'q', ii, '(t)'
3192 call write_iter_header(out_q, aux)
3193 end do
3194 do ii = 1, ks%pt%nmodes
3195 write(aux, '(a1,i3,a3)') 'p', ii, '(t)'
3196 call write_iter_header(out_q, aux)
3197 end do
3198 do ii = 1, space%dim
3199 write(aux, '(a3,i3,a3)') 'f_pt', ii, '(t)'
3200 call write_iter_header(out_q, aux)
3201 end do
3202 call write_iter_nl(out_q)
3203 call td_write_print_header_end(out_q)
3204 end if
3205
3206 call write_iter_start(out_q)
3207 call write_iter_double(out_q, ks%pt_mx%pt_q, ks%pt%nmodes)
3208 call write_iter_double(out_q, ks%pt_mx%pt_p, ks%pt%nmodes)
3209 call write_iter_double(out_q, ks%pt_mx%fmf, space%dim)
3210 call write_iter_nl(out_q)
3211 end if
3212
3213 pop_sub(td_write_q)
3214 end subroutine td_write_q
3215
3216
3217 ! ---------------------------------------------------------
3218 subroutine td_write_mxll_field(out_mxll, space, hm, dt, iter)
3219 type(c_ptr), intent(inout) :: out_mxll
3220 class(space_t), intent(in) :: space
3221 type(hamiltonian_elec_t),intent(in) :: hm
3222 real(real64), intent(in) :: dt
3223 integer, intent(in) :: iter
3224
3225 integer :: idir
3226 real(real64) :: field(space%dim)
3227 character(len=80) :: aux
3228 character(len=1) :: field_char
3229
3230 push_sub(td_write_mxll_field)
3231
3232 if (.not. mpi_grp_is_root(mpi_world)) then
3233 pop_sub(td_write_mxll_field)
3234 return ! only first node outputs
3235 end if
3236
3237 if (iter == 0) then
3238 call td_write_print_header_init(out_mxll)
3239
3240 write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, dt), &
3241 " [", trim(units_abbrev(units_out%time)), "]"
3242 call write_iter_string(out_mxll, aux)
3243 call write_iter_nl(out_mxll)
3244
3245 call write_iter_header_start(out_mxll)
3246 select case (hm%mxll%coupling_mode)
3247 case (length_gauge_dipole, multipolar_expansion)
3248 if (hm%mxll%add_electric_dip) field_char = 'E'
3249 if (hm%mxll%add_magnetic_dip) field_char = 'B'
3250 do idir = 1, space%dim
3251 write(aux, '(a,i1,a)') field_char // '(', idir, ')'
3252 call write_iter_header(out_mxll, aux)
3253 end do
3254 case (velocity_gauge_dipole)
3255 do idir = 1, space%dim
3256 write(aux, '(a,i1,a)') 'A(', idir, ')'
3257 call write_iter_header(out_mxll, aux)
3258 end do
3259 end select
3260 call write_iter_nl(out_mxll)
3261
3262 call write_iter_string(out_mxll, '#[Iter n.]')
3263 call write_iter_header(out_mxll, '[' // trim(units_abbrev(units_out%time)) // ']')
3264
3265 ! Note that we do not print out units of E, A or B, but rather units of e*E, e*A or e*B
3266 ! (force, energy and magnetic flux density, respectively).
3267 select case (hm%mxll%coupling_mode)
3268 case (length_gauge_dipole, multipolar_expansion)
3269 if (hm%mxll%add_electric_dip) aux = '[' // trim(units_abbrev(units_out%force)) // ']'
3270 if (hm%mxll%add_magnetic_dip) aux = '[' // trim(units_abbrev(unit_one/units_out%length**2)) // ']'
3271 do idir = 1, space%dim
3272 call write_iter_header(out_mxll, aux)
3273 end do
3274 case (velocity_gauge_dipole)
3275 aux = '[' // trim(units_abbrev(units_out%energy)) // ']'
3276 do idir = 1, space%dim
3277 call write_iter_header(out_mxll, aux)
3278 end do
3279 end select
3280 call write_iter_nl(out_mxll)
3281 call td_write_print_header_end(out_mxll)
3282 end if
3283
3284 call write_iter_start(out_mxll)
3285
3286 field = m_zero
3287 select case (hm%mxll%coupling_mode)
3288 case (length_gauge_dipole, multipolar_expansion)
3289 if (hm%mxll%add_electric_dip) field = units_from_atomic(units_out%force, hm%mxll%e_field_dip)
3290 if (hm%mxll%add_magnetic_dip) field = units_from_atomic(unit_one/units_out%length**2, hm%mxll%b_field_dip)
3291 call write_iter_double(out_mxll, field, space%dim)
3292 case (velocity_gauge_dipole)
3293 field = units_from_atomic(units_out%energy, hm%mxll%vec_pot_dip)
3294 call write_iter_double(out_mxll, field, space%dim)
3295 end select
3296 call write_iter_nl(out_mxll)
3297
3298 pop_sub(td_write_mxll_field)
3299 end subroutine td_write_mxll_field
3300
3301
3302 ! ---------------------------------------------------------
3303 subroutine td_write_effective_u(out_coords, lda_u, iter)
3304 type(c_ptr), intent(inout) :: out_coords
3305 type(lda_u_t), intent(in) :: lda_u
3306 integer, intent(in) :: iter
3307
3308 integer :: ios, inn
3309 character(len=50) :: aux
3310
3311 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
3312
3313 push_sub(td_write_effective_u)
3314
3315 if (iter == 0) then
3316 call td_write_print_header_init(out_coords)
3317
3318 ! first line: column names
3319 call write_iter_header_start(out_coords)
3320
3321 do ios = 1, lda_u%norbsets
3322 write(aux, '(a2,i3,a1)') 'Ueff(', ios, ')'
3323 call write_iter_header(out_coords, aux)
3324 end do
3325
3326 do ios = 1, lda_u%norbsets
3327 write(aux, '(a2,i3,a1)') 'U(', ios, ')'
3328 call write_iter_header(out_coords, aux)
3329 end do
3330
3331 do ios = 1, lda_u%norbsets
3332 write(aux, '(a2,i3,a1)') 'J(', ios, ')'
3333 call write_iter_header(out_coords, aux)
3334 end do
3335
3336 if (lda_u%intersite) then
3337 do ios = 1, lda_u%norbsets
3338 do inn = 1, lda_u%orbsets(ios)%nneighbors
3339 write(aux, '(a2,i3,a1,i3,a1)') 'V(', ios,'-', inn, ')'
3340 call write_iter_header(out_coords, aux)
3341 end do
3342 end do
3343 end if
3344
3345
3346 call write_iter_nl(out_coords)
3347
3348 ! second line: units
3349 call write_iter_string(out_coords, '#[Iter n.]')
3350 call write_iter_header(out_coords, '[' // trim(units_abbrev(units_out%time)) // ']')
3351 call write_iter_string(out_coords, &
3352 'Effective U ' // trim(units_abbrev(units_out%energy)) // &
3353 ', U in '// trim(units_abbrev(units_out%energy)) // &
3354 ', J in ' // trim(units_abbrev(units_out%energy)))
3355 call write_iter_nl(out_coords)
3356
3357 call td_write_print_header_end(out_coords)
3358 end if
3359
3360 call write_iter_start(out_coords)
3361
3362 do ios = 1, lda_u%norbsets
3363 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3364 lda_u%orbsets(ios)%Ueff), 1)
3365 end do
3366
3367 do ios = 1, lda_u%norbsets
3368 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3369 lda_u%orbsets(ios)%Ubar), 1)
3370 end do
3371
3372 do ios = 1, lda_u%norbsets
3373 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3374 lda_u%orbsets(ios)%Jbar), 1)
3375 end do
3376
3377 if (lda_u%intersite) then
3378 do ios = 1, lda_u%norbsets
3379 do inn = 1, lda_u%orbsets(ios)%nneighbors
3380 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3381 lda_u%orbsets(ios)%V_ij(inn,0)), 1)
3382 end do
3383 end do
3384 end if
3385
3386 call write_iter_nl(out_coords)
3387
3388 pop_sub(td_write_effective_u)
3389 end subroutine td_write_effective_u
3390
3392 subroutine td_write_norm_ks_orbitals(file_handle, grid, kpoints, st, iter)
3393 type(c_ptr), intent(inout) :: file_handle
3394 type(grid_t), intent(in) :: grid
3395 type(kpoints_t), intent(in) :: kpoints
3396 type(states_elec_t), intent(in) :: st
3397 integer, intent(in) :: iter
3398
3399 integer :: ik_ispin, ist
3400 character(len=7) :: nkpt_str, nst_str
3401 character(len=7) :: ik_str, ist_str
3402 real(real64), allocatable :: norm_ks(:, :)
3403 real(real64) :: n_electrons
3404
3406
3407 safe_allocate(norm_ks(1:st%nst, 1:st%nik))
3408 call states_elec_calc_norms(grid, kpoints, st, norm_ks)
3409
3410 if (mpi_grp_is_root(mpi_world)) then
3411 ! Header
3412 if (iter == 0) then
3413 call td_write_print_header_init(file_handle)
3414
3415 ! Line 1
3416 write(nkpt_str, '(I7)') st%nik
3417 write(nst_str, '(I7)') st%nst
3418 call write_iter_string(file_handle, '# Dimensions. (nstates, nkpt * nspin):')
3419 call write_iter_string(file_handle, trim(adjustl(nst_str)) // ' ' // trim(adjustl(nkpt_str)))
3420 call write_iter_nl(file_handle)
3421
3422 ! Line 2
3423 call write_iter_string(file_handle, '# Norm ordering: (istate, ikpoint_spin)')
3424 call write_iter_nl(file_handle)
3425
3426 ! Line 3
3427 call write_iter_header_start(file_handle)
3428 call write_iter_header(file_handle, 'N_electrons')
3429 do ik_ispin = 1, st%nik
3430 do ist = 1, st%nst
3431 write(ik_str, '(I7)') ik_ispin
3432 write(ist_str, '(I7)') ist
3433 call write_iter_header(file_handle, &
3434 'Norm (' // trim(ist_str) // ',' // trim(ik_str) // ')')
3435 enddo
3436 enddo
3437 call write_iter_nl(file_handle)
3438 call td_write_print_header_end(file_handle)
3439 end if
3440
3441 n_electrons = sum(st%occ * norm_ks**2)
3442
3443 ! Output norms for time step `iter`
3444 call write_iter_start(file_handle)
3445 call write_iter_double(file_handle, n_electrons, 1)
3446 do ik_ispin = 1, st%nik
3447 call write_iter_double(file_handle, norm_ks(:, ik_ispin), size(norm_ks, 1))
3448 enddo
3449 call write_iter_nl(file_handle)
3450
3451 end if
3452
3453 safe_deallocate_a(norm_ks)
3454
3456
3457 end subroutine td_write_norm_ks_orbitals
3458
3459 ! ---------------------------------------------------------
3460 subroutine td_write_mxll_init(writ, namespace, iter, dt)
3461 type(td_write_t), intent(out) :: writ
3462 type(namespace_t), intent(in) :: namespace
3463 integer, intent(in) :: iter
3464 real(real64), intent(in) :: dt
3465
3466 integer :: default, flags, iout, first
3467
3468 push_sub(td_write_mxll_init)
3469
3470 !%Variable MaxwellTDOutput
3471 !%Type flag
3472 !%Default maxwell_energy
3473 !%Section Maxwell::Output
3474 !%Description
3475 !% Defines what should be output during the time-dependent
3476 !% Maxwell simulation. Many of the options can increase the computational
3477 !% cost of the simulation, so only use the ones that you need. In
3478 !% most cases the default value is enough, as it is adapted to the
3479 !% details of the TD run.
3480 !% WARNING: the calculation of the longitudinal or transverse E and B fields
3481 !% can be very expensive, so please consider using the MaxwellOutput block
3482 !% to calculate and output these quantities at certain timesteps.
3483 !%Option maxwell_total_e_field 1
3484 !% Output of the total (longitudinal plus transverse) electric field at
3485 !% the points specified in the MaxwellFieldsCoordinate block
3486 !%Option maxwell_total_b_field 8
3487 !% Output of the total (longitudinal plus transverse) magnetic field at
3488 !% the points specified in the MaxwellFieldsCoordinate block
3489 !%Option maxwell_longitudinal_e_field 64
3490 !% Output of the longitudinal electric field at the points
3491 !% specified in the MaxwellFieldsCoordinate block (can slow down the run)
3492 !%Option maxwell_longitudinal_b_field 512
3493 !% Output of the longitudinal magnetic field at the points
3494 !% specified in the MaxwellFieldsCoordinate block (can slow down the run)
3495 !%Option maxwell_transverse_e_field 4096
3496 !% Output of the transverse electric field at the points
3497 !% specified in the MaxwellFieldsCoordinate block (can slow down the run)
3498 !%Option maxwell_transverse_b_field 32768
3499 !% Output of the transverse magnetic field at the points
3500 !% specified in the MaxwellFieldsCoordinate block (can slow down the run)
3501 !%Option maxwell_energy 262144
3502 !% Output of the electromagnetic field energy into the folder <tt>td.general/maxwell</tt>.
3503 !% WARNING: the transverse and longitudinal energies will be correct only if you request
3504 !% the longitudinal and transverse E or B fields as output. Otherwise they will be set to
3505 !% zero.
3506 !%Option e_field_surface_x 524288
3507 !% Output of the E field sliced along the plane x=0 for each field component
3508 !%Option e_field_surface_y 1048576
3509 !% Output of the E field sliced along the plane y=0 for each field component
3510 !%Option e_field_surface_z 2097152
3511 !% Output of the E field sliced along the plane z=0 for each field component
3512 !%Option b_field_surface_x 4194304
3513 !% Output of the B field sliced along the plane x=0 for each field component
3514 !%Option b_field_surface_y 8388608
3515 !% Output of the B field sliced along the plane y=0 for each field component
3516 !%Option b_field_surface_z 16777216
3517 !% Output of the B field sliced along the plane z=0 for each field component
3518 !%End
3519
3520 default = 2**(out_maxwell_energy - 1)
3521 call parse_variable(namespace, 'MaxwellTDOutput', default, flags)
3522
3523 if (.not. varinfo_valid_option('MaxwellTDOutput', flags, is_flag = .true.)) then
3524 call messages_input_error(namespace, 'MaxwellTDOutput')
3525 end if
3526
3528 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3529 if (writ%out(iout)%write) then
3530 writ%out(iout + 1)%write = .true.
3531 writ%out(iout + 2)%write = .true.
3532 end if
3533 end do
3534
3536 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3537 end do
3538
3539 if (iter == 0) then
3540 first = 0
3541 else
3542 first = iter + 1
3543 end if
3544
3545 call io_mkdir('td.general', namespace)
3546
3547 ! total E field
3548 if (writ%out(out_maxwell_total_e_field)%write) then
3549 call write_iter_init(writ%out(out_maxwell_total_e_field)%handle, first, &
3550 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_x", namespace)))
3551 call write_iter_init(writ%out(out_maxwell_total_e_field + 1)%handle, first, &
3552 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_y", namespace)))
3553 call write_iter_init(writ%out(out_maxwell_total_e_field + 2)%handle, first, &
3554 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_z", namespace)))
3555 end if
3556
3557 ! total B field
3558 if (writ%out(out_maxwell_total_b_field)%write) then
3559 call write_iter_init(writ%out(out_maxwell_total_b_field)%handle, first, &
3560 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_x", namespace)))
3561 call write_iter_init(writ%out(out_maxwell_total_b_field + 1)%handle, first, &
3562 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_y", namespace)))
3563 call write_iter_init(writ%out(out_maxwell_total_b_field + 2)%handle, first, &
3564 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_z", namespace)))
3565 end if
3566
3567 ! longitudinal E field
3568 if (writ%out(out_maxwell_long_e_field)%write) then
3569 call write_iter_init(writ%out(out_maxwell_long_e_field)%handle, first, &
3570 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_x", namespace)))
3571 call write_iter_init(writ%out(out_maxwell_long_e_field + 1)%handle, first, &
3572 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_y", namespace)))
3573 call write_iter_init(writ%out(out_maxwell_long_e_field + 2)%handle, first, &
3574 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_z", namespace)))
3575 end if
3576
3577 ! longitudinal B field
3578 if (writ%out(out_maxwell_long_b_field)%write) then
3579 call write_iter_init(writ%out(out_maxwell_long_b_field)%handle, first, &
3580 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_x", namespace)))
3581 call write_iter_init(writ%out(out_maxwell_long_b_field + 1)%handle, first, &
3582 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_y", namespace)))
3583 call write_iter_init(writ%out(out_maxwell_long_b_field + 2)%handle, first, &
3584 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_z", namespace)))
3585 end if
3586
3587 ! transverse E field
3588 if (writ%out(out_maxwell_trans_e_field)%write) then
3589 call write_iter_init(writ%out(out_maxwell_trans_e_field)%handle, first, &
3590 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_x", namespace)))
3591 call write_iter_init(writ%out(out_maxwell_trans_e_field + 1)%handle, first, &
3592 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_y", namespace)))
3593 call write_iter_init(writ%out(out_maxwell_trans_e_field + 2)%handle, first, &
3594 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_z", namespace)))
3595 end if
3596
3597 ! transverse B field
3598 if (writ%out(out_maxwell_trans_b_field)%write) then
3599 call write_iter_init(writ%out(out_maxwell_trans_b_field)%handle, first, &
3600 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_x", namespace)))
3601 call write_iter_init(writ%out(out_maxwell_trans_b_field + 1)%handle, first, &
3602 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_y", namespace)))
3603 call write_iter_init(writ%out(out_maxwell_trans_b_field + 2)%handle, first, &
3604 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_z", namespace)))
3605 end if
3606
3607 ! Maxwell energy
3608 if (writ%out(out_maxwell_energy)%write) then
3609 call write_iter_init(writ%out(out_maxwell_energy)%handle, first, &
3610 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/maxwell_energy", namespace)))
3611 end if
3612
3613 if (writ%out(out_e_field_surface_x)%write) then
3614 call write_iter_init(writ%out(out_e_field_surface_x)%handle, first, &
3615 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-x", namespace)))
3616 end if
3617
3618 if (writ%out(out_e_field_surface_y)%write) then
3619 call write_iter_init(writ%out(out_e_field_surface_y)%handle, first, &
3620 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-y", namespace)))
3621 end if
3622
3623 if (writ%out(out_e_field_surface_z)%write) then
3624 call write_iter_init(writ%out(out_e_field_surface_z)%handle, first, &
3625 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-z", namespace)))
3626 end if
3627
3628 if (writ%out(out_b_field_surface_x)%write) then
3629 call write_iter_init(writ%out(out_b_field_surface_x)%handle, first, &
3630 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-x", namespace)))
3631 end if
3632
3633 if (writ%out(out_b_field_surface_y)%write) then
3634 call write_iter_init(writ%out(out_b_field_surface_y)%handle, first, &
3635 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-y", namespace)))
3636 end if
3637
3638 if (writ%out(out_b_field_surface_z)%write) then
3639 call write_iter_init(writ%out(out_b_field_surface_z)%handle, first, &
3640 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-z", namespace)))
3641 end if
3642
3643 pop_sub(td_write_mxll_init)
3644 end subroutine td_write_mxll_init
3645
3646
3647 ! ---------------------------------------------------------
3648 subroutine td_write_mxll_end(writ)
3649 type(td_write_t), intent(inout) :: writ
3650
3651 integer :: iout
3652
3653 push_sub(td_write_mxll_end)
3654
3655 do iout = 1, out_maxwell_max
3656 if (writ%out(iout)%write) call write_iter_end(writ%out(iout)%handle)
3657 end do
3658
3659 pop_sub(td_write_mxll_end)
3660 end subroutine td_write_mxll_end
3661
3662
3663 ! ---------------------------------------------------------
3664 subroutine td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
3665 type(td_write_t), intent(inout) :: writ
3666 class(space_t), intent(in) :: space
3667 type(grid_t), intent(inout) :: gr
3668 type(states_mxll_t), intent(inout) :: st
3669 type(hamiltonian_mxll_t), intent(inout) :: hm
3670 type(helmholtz_decomposition_t), intent(inout) :: helmholtz
3671 real(real64), intent(in) :: dt
3672 integer, intent(in) :: iter
3673 type(namespace_t), intent(in) :: namespace
3674
3675
3676 push_sub(td_write_mxll_iter)
3677
3678 call profiling_in("TD_WRITE_ITER_MAXWELL")
3679
3680 if (writ%out(out_maxwell_trans_e_field)%write .or. writ%out(out_maxwell_trans_b_field)%write) then
3681 call helmholtz%get_trans_field(namespace, st%rs_state_trans, total_field=st%rs_state)
3682 call get_rs_state_at_point(st%selected_points_rs_state_trans(:,:), st%rs_state_trans, &
3683 st%selected_points_coordinate(:,:), st, gr)
3684 !TODO: calculate transverse energy
3685 else
3686 hm%energy%energy_trans = m_zero
3687 end if
3688
3689 if (writ%out(out_maxwell_long_e_field)%write .or. writ%out(out_maxwell_long_b_field)%write) then
3690 call helmholtz%get_long_field(namespace, st%rs_state_long, total_field=st%rs_state)
3691 call get_rs_state_at_point(st%selected_points_rs_state_long(:,:), st%rs_state_long, &
3692 st%selected_points_coordinate(:,:), st, gr)
3693 !TODO: calculate longitudinal energy
3694 else
3695 hm%energy%energy_long = m_zero
3696 end if
3697
3698 ! total E field
3699 if (writ%out(out_maxwell_total_e_field)%write) then
3700 call td_write_fields(writ%out(out_maxwell_total_e_field)%handle, space, st, iter, dt, &
3702 call td_write_fields(writ%out(out_maxwell_total_e_field + 1)%handle, space, st, iter, dt, &
3704 call td_write_fields(writ%out(out_maxwell_total_e_field + 2)%handle, space, st, iter, dt, &
3706 end if
3707
3708 ! total B field
3709 if (writ%out(out_maxwell_total_b_field)%write) then
3710 call td_write_fields(writ%out(out_maxwell_total_b_field)%handle, space, st, iter, dt, &
3712 call td_write_fields(writ%out(out_maxwell_total_b_field + 1)%handle, space, st, iter, dt, &
3714 call td_write_fields(writ%out(out_maxwell_total_b_field + 2)%handle, space, st, iter, dt, &
3716 end if
3717
3718 ! Longitudinal E field
3719 if (writ%out(out_maxwell_long_e_field)%write) then
3720 call td_write_fields(writ%out(out_maxwell_long_e_field)%handle, space, st, iter, dt, &
3722 call td_write_fields(writ%out(out_maxwell_long_e_field + 1)%handle, space, st, iter, dt, &
3724 call td_write_fields(writ%out(out_maxwell_long_e_field + 2)%handle, space, st, iter, dt, &
3726 end if
3727
3728 ! Longitudinal B field
3729 if (writ%out(out_maxwell_long_b_field)%write) then
3730 call td_write_fields(writ%out(out_maxwell_long_b_field)%handle, space, st, iter, dt, &
3732 call td_write_fields(writ%out(out_maxwell_long_b_field + 1)%handle, space, st, iter, dt, &
3734 call td_write_fields(writ%out(out_maxwell_long_b_field + 2)%handle, space, st, iter, dt, &
3736 end if
3737
3738 ! Transverse E field
3739 if (writ%out(out_maxwell_trans_e_field)%write) then
3740 call td_write_fields(writ%out(out_maxwell_trans_e_field)%handle, space, st, iter, dt, &
3742 call td_write_fields(writ%out(out_maxwell_trans_e_field + 1)%handle, space, st, iter, dt, &
3744 call td_write_fields(writ%out(out_maxwell_trans_e_field + 2)%handle, space, st, iter, dt, &
3746 end if
3747
3748 ! Transverse B field
3749 if (writ%out(out_maxwell_trans_b_field)%write) then
3750 call td_write_fields(writ%out(out_maxwell_trans_b_field)%handle, space, st, iter, dt, &
3752 call td_write_fields(writ%out(out_maxwell_trans_b_field + 1)%handle, space, st, iter, dt, &
3754 call td_write_fields(writ%out(out_maxwell_trans_b_field + 2)%handle, space, st, iter, dt, &
3756 end if
3758 ! Maxwell energy
3759 if (writ%out(out_maxwell_energy)%write) then
3760 call td_write_maxwell_energy(writ%out(out_maxwell_energy)%handle, hm, iter)
3761 end if
3762
3763 if (writ%out(out_e_field_surface_x)%write) then
3764 call td_write_electric_field_box_surface(writ%out(out_e_field_surface_x)%handle, st, 1, iter)
3765 end if
3766
3767 if (writ%out(out_e_field_surface_y)%write) then
3768 call td_write_electric_field_box_surface(writ%out(out_e_field_surface_y)%handle, st, 2, iter)
3769 end if
3770
3771 if (writ%out(out_e_field_surface_z)%write) then
3772 call td_write_electric_field_box_surface(writ%out(out_e_field_surface_z)%handle, st, 3, iter)
3773 end if
3774
3775 if (writ%out(out_b_field_surface_x)%write) then
3776 call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_x)%handle, st, 1, iter)
3777 end if
3778
3779 if (writ%out(out_b_field_surface_y)%write) then
3780 call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_y)%handle, st, 2, iter)
3781 end if
3782
3783 if (writ%out(out_b_field_surface_z)%write) then
3784 call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_z)%handle, st, 3, iter)
3785 end if
3786
3787 call profiling_out("TD_WRITE_ITER_MAXWELL")
3788
3789 pop_sub(td_write_mxll_iter)
3790 end subroutine td_write_mxll_iter
3791
3792
3793 ! ---------------------------------------------------------
3794 subroutine td_write_maxwell_energy(out_maxwell_energy, hm, iter)
3795 type(c_ptr), intent(inout) :: out_maxwell_energy
3796 type(hamiltonian_mxll_t), intent(in) :: hm
3797 integer, intent(in) :: iter
3798
3799 integer :: ii
3800
3801 integer :: n_columns
3802
3803 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
3804
3805 push_sub(td_write_maxwell_energy)
3806
3807 n_columns = 8
3808
3809 if (iter == 0) then
3810 call td_write_print_header_init(out_maxwell_energy)
3811
3812 ! first line -> column names
3813 call write_iter_header_start(out_maxwell_energy)
3814 call write_iter_header(out_maxwell_energy, 'Total')
3815 call write_iter_header(out_maxwell_energy, 'E**2')
3816 call write_iter_header(out_maxwell_energy, 'B**2')
3817 call write_iter_header(out_maxwell_energy, 'Total+Boundaries')
3818 call write_iter_header(out_maxwell_energy, 'Boundaries')
3819 call write_iter_header(out_maxwell_energy, 'Transversal')
3820 call write_iter_header(out_maxwell_energy, 'Longitudinal')
3821 call write_iter_header(out_maxwell_energy, 'Incident Waves')
3822
3823 call write_iter_nl(out_maxwell_energy)
3824
3825 ! units
3826
3827 call write_iter_string(out_maxwell_energy, '#[Iter n.]')
3828 call write_iter_header(out_maxwell_energy, '[' // trim(units_abbrev(units_out%time)) // ']')
3829
3830 do ii = 1, n_columns
3831 call write_iter_header(out_maxwell_energy, '[' // trim(units_abbrev(units_out%energy)) // ']')
3832 end do
3833 call write_iter_nl(out_maxwell_energy)
3834
3835 call td_write_print_header_end(out_maxwell_energy)
3836 end if
3837
3838 call write_iter_start(out_maxwell_energy)
3839 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy), 1)
3840 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%e_energy), 1)
3841 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%b_energy), 1)
3842 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, &
3843 hm%energy%energy+hm%energy%boundaries), 1)
3844 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%boundaries), 1)
3845 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_trans), 1)
3846 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_long), 1)
3847 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_plane_waves), 1)
3848 call write_iter_nl(out_maxwell_energy)
3849
3851 end subroutine td_write_maxwell_energy
3852
3853
3854 ! ---------------------------------------------------------
3855 subroutine td_write_electric_field_box_surface(out_field_surf, st, dim, iter)
3856 type(c_ptr), intent(inout) :: out_field_surf
3857 type(states_mxll_t), intent(in) :: st
3858 integer, intent(in) :: dim
3859 integer, intent(in) :: iter
3860
3861 integer :: ii
3862
3863 integer :: n_columns
3864
3865 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
3866
3868
3869 n_columns = 12
3870
3871 if (iter == 0) then
3872 call td_write_print_header_init(out_field_surf)
3873
3874 ! first line -> column names
3875 call write_iter_header_start(out_field_surf)
3876 call write_iter_header(out_field_surf, '- x direction')
3877 call write_iter_header(out_field_surf, '+ x direction')
3878 call write_iter_header(out_field_surf, '- y direction')
3879 call write_iter_header(out_field_surf, '+ y direction')
3880 call write_iter_header(out_field_surf, '- z direction')
3881 call write_iter_header(out_field_surf, '+ z direction')
3882 call write_iter_header(out_field_surf, '- x dir. p. w.')
3883 call write_iter_header(out_field_surf, '+ x dir. p. w.')
3884 call write_iter_header(out_field_surf, '- y dir. p. w.')
3885 call write_iter_header(out_field_surf, '+ y dir. p. w.')
3886 call write_iter_header(out_field_surf, '- z dir. p. w.')
3887 call write_iter_header(out_field_surf, '+ z dir. p. w.')
3888
3889 call write_iter_nl(out_field_surf)
3890
3891 ! units
3892 call write_iter_string(out_field_surf, '#[Iter n.]')
3893 call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%time)) // ']')
3894
3895 do ii = 1, n_columns
3896 call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%energy/units_out%length)) // ']')
3897 end do
3898 call write_iter_nl(out_field_surf)
3899
3900 call td_write_print_header_end(out_field_surf)
3901 end if
3902
3903 call write_iter_start(out_field_surf)
3904 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3905 st%electric_field_box_surface(1,1,dim)), 1)
3906 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3907 st%electric_field_box_surface(2,1,dim)), 1)
3908 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3909 st%electric_field_box_surface(1,2,dim)), 1)
3910 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3911 st%electric_field_box_surface(2,2,dim)), 1)
3912 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3913 st%electric_field_box_surface(1,3,dim)), 1)
3914 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3915 st%electric_field_box_surface(2,3,dim)), 1)
3916 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3917 st%electric_field_box_surface_plane_waves(1,1,dim)), 1)
3918 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3919 st%electric_field_box_surface_plane_waves(2,1,dim)), 1)
3920 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3921 st%electric_field_box_surface_plane_waves(1,2,dim)), 1)
3922 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3923 st%electric_field_box_surface_plane_waves(2,2,dim)), 1)
3924 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3925 st%electric_field_box_surface_plane_waves(1,3,dim)), 1)
3926 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3927 st%electric_field_box_surface_plane_waves(2,3,dim)), 1)
3928 call write_iter_nl(out_field_surf)
3929
3932
3933
3934 ! ---------------------------------------------------------
3935 subroutine td_write_magnetic_field_box_surface(out_field_surf, st, dim, iter)
3936 type(c_ptr), intent(inout) :: out_field_surf
3937 type(states_mxll_t), intent(in) :: st
3938 integer, intent(in) :: dim
3939 integer, intent(in) :: iter
3940
3941 integer :: ii
3942
3943 integer :: n_columns
3944
3945 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
3946
3949 n_columns = 12
3950
3951 if (iter == 0) then
3952 call td_write_print_header_init(out_field_surf)
3953
3954 ! first line -> column names
3955 call write_iter_header_start(out_field_surf)
3956 call write_iter_header(out_field_surf, '- x direction')
3957 call write_iter_header(out_field_surf, '+ x direction')
3958 call write_iter_header(out_field_surf, '- y direction')
3959 call write_iter_header(out_field_surf, '+ y direction')
3960 call write_iter_header(out_field_surf, '- z direction')
3961 call write_iter_header(out_field_surf, '+ z direction')
3962 call write_iter_header(out_field_surf, '- x dir. p. w.')
3963 call write_iter_header(out_field_surf, '+ x dir. p. w.')
3964 call write_iter_header(out_field_surf, '- y dir. p. w.')
3965 call write_iter_header(out_field_surf, '+ y dir. p. w.')
3966 call write_iter_header(out_field_surf, '- z dir. p. w.')
3967 call write_iter_header(out_field_surf, '+ z dir. p. w.')
3968
3969 call write_iter_nl(out_field_surf)
3970
3971 ! units
3972 call write_iter_string(out_field_surf, '#[Iter n.]')
3973 call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%time)) // ']')
3974
3975 do ii = 1, n_columns
3976 call write_iter_header(out_field_surf, '[' // trim(units_abbrev(unit_one/units_out%length**2)) // ']')
3977 end do
3978 call write_iter_nl(out_field_surf)
3979
3980 call td_write_print_header_end(out_field_surf)
3981 end if
3982
3983 call write_iter_start(out_field_surf)
3984 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3985 st%magnetic_field_box_surface(1,1,dim)), 1)
3986 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3987 st%magnetic_field_box_surface(2,1,dim)), 1)
3988 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3989 st%magnetic_field_box_surface(1,2,dim)), 1)
3990 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3991 st%magnetic_field_box_surface(2,2,dim)), 1)
3992 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3993 st%magnetic_field_box_surface(1,3,dim)), 1)
3994 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3995 st%magnetic_field_box_surface(2,3,dim)), 1)
3996 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3997 st%magnetic_field_box_surface_plane_waves(1,1,dim)), 1)
3998 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
3999 st%magnetic_field_box_surface_plane_waves(2,1,dim)), 1)
4000 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4001 st%magnetic_field_box_surface_plane_waves(1,2,dim)), 1)
4002 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4003 st%magnetic_field_box_surface_plane_waves(2,2,dim)), 1)
4004 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4005 st%magnetic_field_box_surface_plane_waves(1,3,dim)), 1)
4006 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4007 st%magnetic_field_box_surface_plane_waves(2,3,dim)), 1)
4008 call write_iter_nl(out_field_surf)
4009
4012
4013 ! ---------------------------------------------------------
4014 subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
4015 type(c_ptr), intent(inout) :: out_fields
4016 class(space_t), intent(in) :: space
4017 type(states_mxll_t), intent(in) :: st
4018 integer, intent(in) :: iter
4019 real(real64), intent(in) :: dt
4020 integer, intent(in) :: e_or_b_field
4021 integer, intent(in) :: field_type
4022 integer, intent(in) :: idir
4023
4024
4025 integer :: id
4026 real(real64) :: field(space%dim), selected_field
4027 character(len=80) :: aux
4029 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
4030
4031 push_sub(td_write_fields)
4032
4033 if (iter == 0) then
4034 call td_write_print_header_init(out_fields)
4035
4036 ! first line
4037 write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, dt), &
4038 " [", trim(units_abbrev(units_out%time)), "]"
4039 call write_iter_string(out_fields, aux)
4040 call write_iter_nl(out_fields)
4041
4042 call write_iter_header_start(out_fields)
4043
4044 do id = 1, st%selected_points_number
4045 select case (e_or_b_field)
4046 case (maxwell_e_field)
4047 write(aux, '(a,i1,a)') 'E(', id, ')'
4048 case (maxwell_b_field)
4049 write(aux, '(a,i1,a)') 'B(', id, ')'
4050 end select
4051 call write_iter_header(out_fields, aux)
4052 end do
4053
4054 call write_iter_nl(out_fields)
4055 call write_iter_string(out_fields, '#[Iter n.]')
4056 call write_iter_header(out_fields, ' [' // trim(units_abbrev(units_out%time)) // ']')
4057
4058 ! Note that we do not print out units of E, B, or A, but rather units of e*E, e*B, e*A.
4059 ! (force, force, and energy, respectively). The reason is that the units of E, B or A
4060 ! are ugly.
4061 aux = ' [' // trim(units_abbrev(units_out%force)) // ']'
4062 do id = 1, st%selected_points_number
4063 call write_iter_header(out_fields, aux)
4064 end do
4065 call write_iter_nl(out_fields)
4066 call td_write_print_header_end(out_fields)
4067 end if
4068
4069 call write_iter_start(out_fields)
4070
4071 do id = 1, st%selected_points_number
4072 select case (e_or_b_field)
4073 case (maxwell_e_field)
4074 ! Output of electric field at selected point
4075 select case (field_type)
4076 case (maxwell_total_field)
4077 call get_electric_field_vector(st%selected_points_rs_state(:,id), field(1:st%dim))
4078 case (maxwell_long_field)
4079 call get_electric_field_vector(st%selected_points_rs_state_long(:,id), field(1:st%dim))
4080 case (maxwell_trans_field)
4081 call get_electric_field_vector(st%selected_points_rs_state_trans(:,id), field(1:st%dim))
4082 end select
4083 selected_field = units_from_atomic(units_out%energy/units_out%length, field(idir))
4084 case (maxwell_b_field)
4085 ! Output of magnetic field at selected point
4086 select case (field_type)
4087 case (maxwell_total_field)
4088 call get_magnetic_field_vector(st%selected_points_rs_state(:,id), st%rs_sign, field(1:st%dim))
4089 case (maxwell_long_field)
4090 call get_magnetic_field_vector(st%selected_points_rs_state_long(:,id), st%rs_sign, field(1:st%dim))
4091 case (maxwell_trans_field)
4092 call get_magnetic_field_vector(st%selected_points_rs_state_trans(:,id), st%rs_sign, field(1:st%dim))
4093 end select
4094 selected_field = units_from_atomic(unit_one/units_out%length**2, field(idir))
4095 end select
4096 call write_iter_double(out_fields, selected_field, 1)
4097 end do
4098
4099 call write_iter_nl(out_fields)
4100
4101 pop_sub(td_write_fields)
4102 end subroutine td_write_fields
4103
4104
4105 !----------------------------------------------------------
4106 subroutine td_write_mxll_free_data(writ, namespace, space, gr, st, hm, helmholtz, outp, iter, time)
4107 type(td_write_t), intent(inout) :: writ
4108 type(namespace_t), intent(in) :: namespace
4109 class(space_t), intent(in) :: space
4110 type(grid_t), intent(inout) :: gr
4111 type(states_mxll_t), intent(inout) :: st
4112 type(hamiltonian_mxll_t), intent(inout) :: hm
4113 type(helmholtz_decomposition_t), intent(inout) :: helmholtz
4114 type(output_t), intent(in) :: outp
4115 integer, intent(in) :: iter
4116 real(real64), intent(in) :: time
4117
4118 character(len=256) :: filename
4119
4120 push_sub(td_write_maxwell_free_data)
4121 call profiling_in("TD_WRITE_MAXWELL_DATA")
4122
4123 ! now write down the rest
4124 write(filename, '(a,a,i7.7)') trim(outp%iter_dir),"td.", iter ! name of directory
4125
4126 call output_mxll(outp, namespace, space, gr, st, hm, helmholtz, time, filename)
4127
4128 call profiling_out("TD_WRITE_MAXWELL_DATA")
4129 pop_sub(td_write_maxwell_free_data)
4130 end subroutine td_write_mxll_free_data
4131
4132end module td_write_oct_m
4133
4134!! Local Variables:
4135!! mode: f90
4136!! coding: utf-8
4137!! End:
ssize_t ssize_t write(int __fd, const void *__buf, size_t __n) __attribute__((__access__(__read_only__
constant times a vector plus a vector
Definition: lalg_basic.F90:171
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:186
Sets the iteration number to the C object.
Definition: write_iter.F90:174
Writes to the corresponding file and adds one to the iteration. Must be called after write_iter_init(...
Definition: write_iter.F90:167
double exp(double __x) __attribute__((__nothrow__
This module contains interfaces for BLAS routines You should not use these routines directly....
Definition: blas.F90:118
This module implements a calculator for the density and defines related functions.
Definition: density.F90:120
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
Definition: density.F90:609
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
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
complex(real64), parameter, public m_z0
Definition: global.F90:198
real(real64), parameter, public lmm_r_single_atom
Default local magnetic moments sphere radius for an isolated system.
Definition: global.F90:216
complex(real64), parameter, public m_zi
Definition: global.F90:202
integer, parameter, public max_output_types
Definition: global.F90:148
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
Definition: global.F90:211
This module implements the underlying real-space grid.
Definition: grid.F90:117
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)
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
character(len=max_path_len) function, public io_workpath(path, namespace)
Definition: io.F90:270
subroutine, public io_rm(fname, namespace)
Definition: io.F90:342
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:311
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
real(real64) function, public ion_dynamics_temperature(ions)
This function returns the ionic temperature in energy units.
integer, parameter, public qkickmode_cos
Definition: kick.F90:169
integer, parameter, public qkickmode_none
Definition: kick.F90:169
integer, parameter, public qkickmode_sin
Definition: kick.F90:169
subroutine, public kick_function_get(space, mesh, kick, kick_function, iq, to_interpolate)
Definition: kick.F90:993
subroutine, public kick_write(kick, iunit, out)
Definition: kick.F90:891
integer, parameter, public qkickmode_bessel
Definition: kick.F90:169
subroutine, public lasers_nondipole_laser_field_step(this, field, time)
Retrieves the NDSFA vector_potential correction. The nondipole field is obtained for consecutive time...
Definition: lasers.F90:1150
subroutine, public lasers_set_nondipole_parameters(this, ndfield, nd_integration_time)
Set parameters for nondipole SFA calculation.
Definition: lasers.F90:752
logical function, public lasers_with_nondipole_field(lasers)
Check if a nondipole SFA correction should be computed for the given laser.
Definition: lasers.F90:739
integer, parameter, public e_field_electric
Definition: lasers.F90:177
integer, parameter, public e_field_vector_potential
Definition: lasers.F90:177
integer, parameter, public e_field_scalar_potential
Definition: lasers.F90:177
integer pure elemental function, public laser_kind(laser)
Definition: lasers.F90:715
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...
Definition: lasers.F90:1115
integer, parameter, public e_field_magnetic
Definition: lasers.F90:177
integer, parameter, public dft_u_none
Definition: lda_u.F90:201
integer, parameter, public dft_u_acbn0
Definition: lda_u.F90:201
subroutine, public magnetic_local_moments(mesh, st, ions, boundaries, rho, rr, lmm)
Definition: magnetic.F90:265
subroutine, public magnetic_moment(mesh, st, rho, mm)
Definition: magnetic.F90:183
subroutine, public magnetic_total_magnetization(mesh, st, qq, trans_mag)
Definition: magnetic.F90:334
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public ylmr_real(xx, li, mi, ylm)
This is a Numerical Recipes-based subroutine computes real spherical harmonics ylm at position (x,...
Definition: math.F90:373
This module defines functions over batches of mesh functions.
Definition: mesh_batch.F90:116
This module defines various routines, operating on mesh functions.
complex(real64) function, public zmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
subroutine, public dmf_multipoles(mesh, ff, lmax, multipole, mask)
This routine calculates the multipoles of a function ff.
real(real64) function, public dmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1045
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public modelmb_sym_all_states(space, mesh, st)
This module contains some common usage patterns of MPI routines.
Definition: mpi_lib.F90:115
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
Definition: mpi.F90:430
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:266
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
integer, parameter, public velocity_gauge_dipole
this module contains the low-level part of the output system
Definition: output_low.F90:115
subroutine, public output_modelmb(outp, namespace, space, dir, gr, ions, iter, st)
this module contains the output system
Definition: output.F90:115
subroutine, public output_all(outp, namespace, space, dir, gr, ions, iter, st, hm, ks)
Definition: output.F90:493
subroutine, public output_scalar_pot(outp, namespace, space, dir, mesh, ions, ext_partners, time)
Definition: output.F90:2465
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
integer, parameter, public restart_gs
Definition: restart.F90:200
subroutine, public restart_init(restart, namespace, data_type, type, mc, ierr, mesh, dir, exact)
Initializes a restart object.
Definition: restart.F90:516
integer, parameter, public restart_proj
Definition: restart.F90:200
integer, parameter, public restart_type_load
Definition: restart.F90:245
subroutine, public restart_end(restart)
Definition: restart.F90:722
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:
Definition: td_calc.F90:165
subroutine, public td_calc_ionch(mesh, st, ch, Nch)
Multiple ionization probabilities calculated form the KS orbital densities C. Ullrich,...
Definition: td_calc.F90:319
subroutine, public td_calc_tvel(gr, st, space, kpoints, vel)
Electronic velocity (to calculate harmonic spectrum...) It is calculated as:
Definition: td_calc.F90:287
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)
Definition: td_write.F90:4029
integer, parameter, public out_total_current
Definition: td_write.F90:201
integer, parameter maxwell_b_field
Definition: td_write.F90:297
integer, parameter, public out_maxwell_max
Definition: td_write.F90:289
subroutine td_write_proj(out_proj, space, mesh, ions, st, gs_st, kick, iter)
Definition: td_write.F90:2370
integer, parameter, public out_q
Definition: td_write.F90:201
integer, parameter, public out_mxll_field
Definition: td_write.F90:201
subroutine calc_projections(mesh, st, gs_st, projections)
This subroutine calculates:
Definition: td_write.F90:2663
integer, parameter out_b_field_surface_y
Definition: td_write.F90:274
subroutine td_write_ionch(out_ionch, mesh, st, iter)
Definition: td_write.F90:2277
integer, parameter, public out_tot_m
Definition: td_write.F90:201
integer, parameter, public out_norm_ks
Definition: td_write.F90:201
integer, parameter out_maxwell_trans_b_field
Definition: td_write.F90:274
subroutine td_write_multipole_r(out_multip, space, mesh, ions, st, lmax, kick, rho, iter, mpi_grp)
Write multipoles to the corresponding file.
Definition: td_write.F90:1559
integer, parameter, public out_proj
Definition: td_write.F90:201
integer, parameter, public out_partial_charges
Definition: td_write.F90:201
integer, parameter, public out_separate_coords
Definition: td_write.F90:201
subroutine td_write_total_current(out_total_current, space, mesh, st, iter)
Definition: td_write.F90:3067
subroutine, public td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, dt)
Definition: td_write.F90:1224
subroutine td_write_effective_u(out_coords, lda_u, iter)
Definition: td_write.F90:3397
integer, parameter maxwell_trans_field
Definition: td_write.F90:292
subroutine td_write_acc(out_acc, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
Definition: td_write.F90:1910
subroutine, public td_write_mxll_init(writ, namespace, iter, dt)
Definition: td_write.F90:3554
subroutine, public td_write_mxll_end(writ)
Definition: td_write.F90:3742
subroutine td_write_mxll_field(out_mxll, space, hm, dt, iter)
Definition: td_write.F90:3312
integer, parameter out_b_field_surface_x
Definition: td_write.F90:274
integer, parameter out_maxwell_long_e_field
Definition: td_write.F90:274
integer, parameter, public out_kp_proj
Definition: td_write.F90:201
integer, parameter, public out_magnets
Definition: td_write.F90:201
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...
Definition: td_write.F90:1512
subroutine td_write_electric_field_box_surface(out_field_surf, st, dim, iter)
Definition: td_write.F90:3949
integer, parameter out_e_field_surface_y
Definition: td_write.F90:274
integer, parameter, public out_angular
Definition: td_write.F90:201
subroutine td_write_populations(out_populations, namespace, space, mesh, st, writ, dt, iter)
Definition: td_write.F90:1830
integer, parameter, public out_max
Definition: td_write.F90:201
integer, parameter out_maxwell_long_b_field
Definition: td_write.F90:274
integer, parameter, public out_energy
Definition: td_write.F90:201
integer, parameter, public out_spin
Definition: td_write.F90:201
subroutine td_write_ftchd(out_ftchd, space, mesh, st, kick, iter)
Definition: td_write.F90:1686
subroutine td_write_partial_charges(out_partial_charges, mesh, st, ions, iter)
Definition: td_write.F90:3220
integer, parameter out_dftu_max
For the Maxwell fields we increment in steps of 3 to leave room for x, y, and z output.
Definition: td_write.F90:269
subroutine td_write_laser(out_laser, space, lasers, dt, iter)
Definition: td_write.F90:2011
integer, parameter out_maxwell_total_b_field
Definition: td_write.F90:274
integer, parameter, public out_ftchd
Definition: td_write.F90:201
integer, parameter, public out_separate_velocity
Definition: td_write.F90:201
subroutine td_write_tot_mag(out_magnets, mesh, st, kick, iter)
Definition: td_write.F90:1374
subroutine td_write_q(out_q, space, ks, iter)
Definition: td_write.F90:3269
integer, parameter, public out_floquet
Definition: td_write.F90:201
integer, parameter, public out_acc
Definition: td_write.F90:201
integer, parameter, public out_ion_ch
Definition: td_write.F90:201
integer, parameter maxwell_long_field
Definition: td_write.F90:292
integer, parameter, public out_n_ex
Definition: td_write.F90:201
integer, parameter out_b_field_surface_z
Definition: td_write.F90:274
subroutine td_write_proj_kp(mesh, kpoints, st, gs_st, namespace, iter)
Definition: td_write.F90:2697
integer, parameter, public out_temperature
Definition: td_write.F90:201
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.
Definition: td_write.F90:3486
subroutine, public td_write_data(writ)
Definition: td_write.F90:1190
subroutine, public td_write_mxll_free_data(writ, namespace, space, gr, st, hm, helmholtz, outp, iter, time)
Definition: td_write.F90:4200
subroutine td_write_total_heat_current(write_obj, space, hm, gr, st, iter)
Definition: td_write.F90:3164
integer, parameter out_e_field_surface_z
Definition: td_write.F90:274
integer, parameter maxwell_total_field
Definition: td_write.F90:292
integer, parameter, public out_coords
Definition: td_write.F90:201
integer, parameter out_maxwell_total_e_field
Definition: td_write.F90:274
integer, parameter, public out_laser
Definition: td_write.F90:201
integer, parameter, public out_eigs
Definition: td_write.F90:201
subroutine td_write_vel(out_vel, gr, st, space, kpoints, iter)
Definition: td_write.F90:1962
subroutine td_write_floquet(namespace, space, hm, ext_partners, mesh, st, iter)
Definition: td_write.F90:2791
integer, parameter, public out_total_heat_current
Definition: td_write.F90:201
integer, parameter out_e_field_surface_x
Definition: td_write.F90:274
subroutine, public td_write_kick(outp, namespace, space, mesh, kick, ions, iter)
Definition: td_write.F90:333
subroutine, public td_write_end(writ)
Definition: td_write.F90:978
subroutine td_write_angular(out_angular, namespace, space, gr, ions, hm, st, kick, iter)
Computes and outputs the orbital angular momentum defined by.
Definition: td_write.F90:1441
subroutine td_write_eigs(out_eigs, st, iter)
Definition: td_write.F90:2216
subroutine, public td_write_iter(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, kick, ks, dt, iter, mc, recalculate_gs)
Definition: td_write.F90:1025
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...
Definition: td_write.F90:2546
subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
Definition: td_write.F90:4108
subroutine td_write_spin(out_spin, mesh, st, iter)
Definition: td_write.F90:1265
integer, parameter, public out_vel
Definition: td_write.F90:201
integer, parameter, public out_gauge_field
Definition: td_write.F90:201
subroutine td_write_temperature(out_temperature, ions, iter)
Definition: td_write.F90:1793
integer, parameter maxwell_e_field
Definition: td_write.F90:297
integer, parameter, public out_populations
Definition: td_write.F90:201
subroutine, public td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
Definition: td_write.F90:3758
subroutine td_write_local_magnetic_moments(out_magnets, gr, st, ions, lmm_r, iter)
Definition: td_write.F90:1317
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.
Definition: td_write.F90:371
integer, parameter out_maxwell_energy
Definition: td_write.F90:274
subroutine td_write_energy(out_energy, hm, iter, ke)
Definition: td_write.F90:2133
subroutine td_write_maxwell_energy(out_maxwell_energy, hm, iter)
Definition: td_write.F90:3888
integer, parameter, public out_separate_forces
Definition: td_write.F90:201
integer, parameter out_maxwell_trans_e_field
Definition: td_write.F90:274
type(type_t), public type_cmplx
Definition: types.F90:134
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:223
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.
Definition: utils.F90:118
character pure function, public index2axis(idir)
Definition: utils.F90:202
subroutine, public v_ks_calculate_current(this, calc_cur)
Definition: v_ks.F90:1436
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Definition: v_ks.F90:730
Explicit interfaces to C functions, defined in write_iter_low.cc.
Definition: write_iter.F90:114
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:168
Describes mesh distribution to nodes.
Definition: mesh.F90:186
This is defined even when running serial.
Definition: mpi.F90:142
Stores all communicators and groups.
Definition: multicomm.F90:206
output handler class
Definition: output_low.F90:164
The states_elec_t class contains all electronic wave functions.
Time-dependent Write Properties.
Definition: td_write.F90:315
int true(void)
subroutine dipole_matrix_elements(dir)
Definition: td_write.F90:2490