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