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, &
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 :: &
202 MAXWELL_TOTAL_FIELD = 1, &
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
239
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
305 type(restart_t) :: restart_gs
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 (integral 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 (integral 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_gs%init(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
662 call restart_gs%end()
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_world%is_root()) 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_world%is_root()
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 (writ%out(out_multipoles)%mpi_grp%is_root()) then
825 do ist = st%st_start, st%st_end
826 write(filename, '(a,i4.4)') 'td.general/multipoles-ist', ist
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_world%is_root()) 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 (writ%out(iout)%mpi_grp%is_root()) then
904 if (writ%out(iout)%resolve_states) then
905 do ist = writ%out(iout)%hand_start, writ%out(iout)%hand_end
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_world%is_root()) 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
986
987 if (writ%out(out_proj)%write .and. mod(iter, writ%compute_interval) == 0) then
988 if (mpi_world%is_root()) 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
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_world%is_root()) call write_iter_set(writ%out(out_n_ex)%handle, iter)
1073 if (recalculate_gs) then ! Load recomputed GS states
1074 call restart_gs%init(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')
1077 call restart_gs%end()
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 (writ%out(iout)%mpi_grp%is_root()) then
1121 if (writ%out(iout)%resolve_states) then
1122 do ii = writ%out(iout)%hand_start, writ%out(iout)%hand_end
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_world%is_root()) 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_world%is_root()) then ! only first node outputs
1201
1202 if (iter == 0) then
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_world%is_root()) 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)
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_world%is_root()) 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)
1386
1387 if (mpi_world%is_root()) then ! Only first node outputs
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
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().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()) 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_world%is_root().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_world%is_root()) 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_world%is_root()) 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_world%is_root()) 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
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)
1845 push_sub(td_write_acc)
1846
1847 if (iter == 0 .and. mpi_world%is_root()) 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_world%is_root()) 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_world%is_root()) 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_world%is_root()) 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
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_world%is_root()) 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
2023
2024 call write_iter_start(out_laser)
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_world%is_root()) 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
2145 if (.not. mpi_world%is_root()) then
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
2228
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_world%is_root()) 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_world%is_root()) 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_world%is_root()) 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_world%is_root()) then
2382 call td_write_print_header_end(out_proj)
2383 end if
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_world%is_root()) 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_world%is_root()) 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_world%is_root()) 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 gs_nst = 1
2522 do ik = 1, st%nik
2523 do ist = 1, gs_st%nst
2524 if (gs_st%occ(ist, ik) > m_min_occ .and. ist > gs_nst) gs_nst = ist
2525 end do
2526 end do
2527
2528 safe_allocate(projections(1:gs_nst, 1:st%nst))
2529
2530 safe_allocate(nex_kpt(1:st%nik))
2531 nex_kpt = m_zero
2532 do ik = st%d%kpt%start, st%d%kpt%end
2533 ikpt = st%d%get_kpoint_index(ik)
2534 call zstates_elec_calc_projections(st, gs_st, namespace, mesh, ik, projections, gs_nst)
2535 do ist = 1, gs_nst
2536 weight = st%kweights(ik) * gs_st%occ(ist, ik)/ st%smear%el_per_state
2537 do uist = st%st_start, st%st_end
2538 nex_kpt(ikpt) = nex_kpt(ikpt) - weight * st%occ(uist, ik) * abs(projections(ist, uist))**2
2539 end do
2540 end do
2541 nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*st%kweights(ik)
2542 end do
2543
2544 if (st%parallel_in_states .or. st%d%kpt%parallel) then
2545 call comm_allreduce(st%st_kpt_mpi_grp, nex_kpt)
2546 end if
2547
2548 nex = sum(nex_kpt)
2549
2550 if (mpi_world%is_root()) then
2551 call write_iter_start(out_nex)
2552 call write_iter_double(out_nex, nex, 1)
2553 call write_iter_nl(out_nex)
2554
2555 ! now write down the k-resolved part
2556 write(dir, '(a,a,i7.7)') trim(outp%iter_dir),"td.", iter ! name of directory
2557
2558 call io_function_output_global_bz(outp%how(option__output__current_kpt) &
2559 + outp%how(option__output__density_kpt), dir, "n_excited_el_kpt", namespace, &
2560 kpoints, nex_kpt, unit_one, err)
2561 end if
2562
2563 safe_deallocate_a(projections)
2564 safe_deallocate_a(nex_kpt)
2565
2566 pop_sub(td_write_n_ex)
2567 end subroutine td_write_n_ex
2568
2569 ! ---------------------------------------------------------
2574 ! ---------------------------------------------------------
2575 subroutine calc_projections(mesh, st, gs_st, projections)
2576 class(mesh_t), intent(in) :: mesh
2577 type(states_elec_t), intent(inout) :: st
2578 type(states_elec_t), intent(in) :: gs_st
2579 complex(real64), intent(inout) :: projections(1:st%nst, gs_st%st_start:gs_st%nst, 1:st%nik)
2580
2581 integer :: uist, ist, ik
2582 complex(real64), allocatable :: psi(:, :), gspsi(:, :)
2583 push_sub(calc_projections)
2584
2585 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
2586 safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
2587
2588 projections(:,:,:) = m_zero
2589
2590 do ik = st%d%kpt%start, st%d%kpt%end
2591 do ist = st%st_start, st%st_end
2592 call states_elec_get_state(st, mesh, ist, ik, psi)
2593 do uist = gs_st%st_start, gs_st%nst
2594 call states_elec_get_state(gs_st, mesh, uist, ik, gspsi)
2595 projections(ist, uist, ik) = zmf_dotp(mesh, st%d%dim, psi, gspsi, reduce = .false.)
2596 end do
2597 end do
2598 end do
2599
2600 safe_deallocate_a(psi)
2601 safe_deallocate_a(gspsi)
2602
2603 call comm_allreduce(st%dom_st_kpt_mpi_grp, projections)
2604
2605 pop_sub(calc_projections)
2606 end subroutine calc_projections
2607
2608
2609 subroutine td_write_proj_kp(mesh, kpoints, st, gs_st, namespace, iter)
2610 class(mesh_t), intent(in) :: mesh
2611 type(kpoints_t), intent(in) :: kpoints
2612 type(states_elec_t), intent(in) :: st
2613 type(states_elec_t), intent(inout) :: gs_st
2614 type(namespace_t), intent(in) :: namespace
2615 integer, intent(in) :: iter
2616
2617 complex(real64), allocatable :: proj(:,:), psi(:,:,:), gs_psi(:,:,:), temp_state(:,:)
2618 character(len=80) :: filename1, filename2
2619 integer :: ik,ist, jst, file, idim, nk_proj
2620
2621 push_sub(td_write_proj_kp)
2622
2623 write(filename1,'(I10)') iter
2624 filename1 = 'td.general/projections_iter_'//trim(adjustl(filename1))
2625 file = 9845623
2626
2627 safe_allocate(proj(1:gs_st%nst, 1:gs_st%nst))
2628 safe_allocate(psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2629 safe_allocate(gs_psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
2630 safe_allocate(temp_state(1:mesh%np,1:gs_st%d%dim))
2631
2632 ! Project only k-points that have a zero weight.
2633 ! Why? It is unlikely that one is interested in the projections
2634 ! of the Monkhorst-Pack kpoints, but instead we assume that
2635 ! the user has specified a k-path with zero weights
2636 nk_proj = kpoints%nik_skip
2637
2638 do ik = kpoints%reduced%npoints-nk_proj+1, kpoints%reduced%npoints
2639 ! reset arrays
2640 psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)= m_zero
2641 gs_psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)= m_zero
2642 ! open file for writing
2643 if (mpi_world%is_root()) then
2644 write(filename2,'(I10)') ik
2645 filename2 = trim(adjustl(filename1))//'_ik_'//trim(adjustl(filename2))
2646 file = io_open(filename2, namespace, action='write')
2647 end if
2648 ! get all states at ik that are locally stored (ground state and td-states)
2649 do ist=gs_st%st_start,gs_st%st_end
2650 if (state_kpt_is_local(gs_st, ist, ik)) then
2651 call states_elec_get_state(st, mesh, ist, ik,temp_state)
2652 do idim = 1,gs_st%d%dim
2653 psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2654 end do
2655 call states_elec_get_state(gs_st, mesh, ist, ik, temp_state)
2656 do idim = 1,gs_st%d%dim
2657 gs_psi(ist,idim,1:mesh%np) = temp_state(1:mesh%np,idim)
2658 end do
2659 end if
2660 end do
2661 ! collect states at ik from all processes in one array
2662 call comm_allreduce(mpi_world, psi)
2663 call comm_allreduce(mpi_world, gs_psi)
2664
2665 ! compute the overlaps as a matrix product
2666 assert(mesh%np_global*gs_st%d%dim < huge(0_int32))
2667 proj(1:gs_st%nst, 1:gs_st%nst) = m_zero
2668 call zgemm('n', &
2669 'c', &
2670 gs_st%nst, &
2671 gs_st%nst, &
2672 i8_to_i4(mesh%np_global*gs_st%d%dim), &
2673 cmplx(mesh%volume_element, m_zero, real64) , &
2674 psi(1, 1, 1), &
2675 ubound(psi, dim = 1), &
2676 gs_psi(1, 1, 1), &
2677 ubound(gs_psi, dim = 1), &
2678 m_z0, &
2679 proj(1, 1), &
2680 ubound(proj, dim = 1))
2681
2682 ! write to file
2683 if (mpi_world%is_root()) then
2684 do ist = 1, gs_st%nst
2685 do jst = 1, gs_st%nst
2686 write(file,'(I3,1x,I3,1x,e13.6,1x,e13.6,2x)') ist, jst, proj(ist,jst)
2687 end do
2688 end do
2689 call io_close(file)
2690 end if
2691
2692 end do! ik
2693
2694 safe_deallocate_a(proj)
2695 safe_deallocate_a(psi)
2696 safe_deallocate_a(gs_psi)
2697 safe_deallocate_a(temp_state)
2698
2699 pop_sub(td_write_proj_kp)
2700 end subroutine td_write_proj_kp
2701
2702 !---------------------------------------
2703 subroutine td_write_floquet(namespace, space, hm, ext_partners, gr, st, iter)
2704 type(namespace_t), intent(in) :: namespace
2705 class(space_t), intent(in) :: space
2706 type(hamiltonian_elec_t), intent(inout) :: hm
2707 type(partner_list_t), intent(in) :: ext_partners
2708 type(grid_t), intent(in) :: gr
2709 type(states_elec_t), intent(inout) :: st
2710 integer, intent(in) :: iter
2711
2712 complex(real64), allocatable :: hmss(:,:), psi(:,:,:), hpsi(:,:,:), temp_state1(:,:)
2713 complex(real64), allocatable :: HFloquet(:,:,:), HFloq_eff(:,:), temp(:,:)
2714 real(real64), allocatable :: eigenval(:), bands(:,:)
2715 character(len=80) :: filename
2716 integer :: it, nT, ik, ist, in, im, file, idim, nik, ik_count
2717 integer :: Forder, Fdim, m0, n0, n1, nst, ii, jj, lim_nst
2718 logical :: downfolding
2719 type(states_elec_t) :: hm_st
2720
2721 real(real64) :: dt, Tcycle, omega
2722
2723 push_sub(td_write_floquet)
2724
2725 downfolding = .false.
2726
2727 ! this does not depend on propagation, so we do it only once
2728 if (.not. iter == 0) then
2729 pop_sub(td_write_floquet)
2730 return
2731 end if
2732
2733 nst = st%nst
2734
2735 !for now no domain distributionallowed
2736 assert(gr%np == gr%np_global)
2737
2738 ! this is used to initialize the hpsi (more effiecient ways?)
2739 call states_elec_copy(hm_st, st)
2740
2741 !%Variable TDFloquetFrequency
2742 !%Type float
2743 !%Default 0
2744 !%Section Time-Dependent::TD Output
2745 !%Description
2746 !% Frequency for the Floquet analysis, this should be the carrier frequency or integer multiples of it.
2747 !% Other options will work, but likely be nonsense.
2748 !%
2749 !%End
2750 call parse_variable(namespace, 'TDFloquetFrequency', m_zero, omega, units_inp%energy)
2751 call messages_print_var_value('Frequency used for Floquet analysis', omega, namespace=namespace)
2752 if (abs(omega) <= m_epsilon) then
2753 message(1) = "Please give a non-zero value for TDFloquetFrequency"
2754 call messages_fatal(1, namespace=namespace)
2755 end if
2756
2757 ! get time of one cycle
2758 tcycle = m_two * m_pi / omega
2759
2760 !%Variable TDFloquetSample
2761 !%Type integer
2762 !%Default 20
2763 !%Section Time-Dependent::TD Output
2764 !%Description
2765 !% Number of points on which one Floquet cycle is sampled in the time-integral of the Floquet analysis.
2766 !%
2767 !%End
2768 call parse_variable(namespace, 'TDFloquetSample',20 ,nt)
2769 call messages_print_var_value('Number of Floquet time-sampling points', nt, namespace=namespace)
2770 dt = tcycle/real(nt, real64)
2771
2772 !%Variable TDFloquetDimension
2773 !%Type integer
2774 !%Default -1
2775 !%Section Time-Dependent::TD Output
2776 !%Description
2777 !% Order of Floquet Hamiltonian. If negative number is given, downfolding is performed.
2778 !%End
2779 call parse_variable(namespace, 'TDFloquetDimension',-1,forder)
2780 if (forder .ge. 0) then
2781 call messages_print_var_value('Order of multiphoton Floquet-Hamiltonian', forder, namespace=namespace)
2782 !Dimension of multiphoton Floquet-Hamiltonian
2783 fdim = 2 * forder + 1
2784 else
2785 message(1) = 'Floquet-Hamiltonian is downfolded'
2786 call messages_info(1, namespace=namespace)
2787 downfolding = .true.
2788 forder = 1
2789 fdim = 3
2790 end if
2791
2792 dt = tcycle/real(nt, real64)
2793
2794 ! we are only interested for k-point with zero weight
2795 nik = hm%kpoints%nik_skip
2796
2797 safe_allocate(hmss(1:nst,1:nst))
2798 safe_allocate( psi(1:nst,1:st%d%dim,1:gr%np))
2799 safe_allocate(hpsi(1:nst,1:st%d%dim,1:gr%np))
2800 safe_allocate(temp_state1(1:gr%np,1:st%d%dim))
2801
2802 ! multiphoton Floquet Hamiltonian, layout:
2803 ! (H_{-n,-m} ... H_{-n,0} ... H_{-n,m})
2804 ! ( . . . . . )
2805 ! H = (H_{0,-m} ... H_{0,0} ... H_{0,m} )
2806 ! ( . . . . . )
2807 ! (H_{n,-m} ... H_{n,0} ... H_{n,m} )
2808 safe_allocate(hfloquet(1:nik,1:nst*fdim, 1:nst*fdim))
2809 hfloquet(1:nik,1:nst*fdim, 1:nst*fdim) = m_zero
2810
2811 ! perform time-integral over one cycle
2812 do it = 1, nt
2813 ! get non-interacting Hamiltonian at time (offset by one cycle to allow for ramp)
2814 call hm%update(gr, namespace, space, ext_partners, time=tcycle+it*dt)
2815 ! get hpsi
2816 call zhamiltonian_elec_apply_all(hm, namespace, gr, st, hm_st)
2817
2818 ! project Hamiltonian into grounstates for zero weight k-points
2819 ik_count = 0
2820
2821 do ik = hm%kpoints%reduced%npoints-nik+1, hm%kpoints%reduced%npoints
2822 ik_count = ik_count + 1
2823
2824 psi(1:nst, 1:st%d%dim, 1:gr%np)= m_zero
2825 hpsi(1:nst, 1:st%d%dim, 1:gr%np)= m_zero
2826
2827 do ist = st%st_start, st%st_end
2828 if (state_kpt_is_local(st, ist, ik)) then
2829 call states_elec_get_state(st, gr, ist, ik,temp_state1)
2830 do idim = 1, st%d%dim
2831 psi(ist, idim, 1:gr%np) = temp_state1(1:gr%np, idim)
2832 end do
2833 call states_elec_get_state(hm_st, gr, ist, ik,temp_state1)
2834 do idim = 1, st%d%dim
2835 hpsi(ist, idim, 1:gr%np) = temp_state1(1:gr%np, idim)
2836 end do
2837 end if
2838 end do
2839 call comm_allreduce(mpi_world, psi)
2840 call comm_allreduce(mpi_world, hpsi)
2841 assert(gr%np_global*st%d%dim < huge(0_int32))
2842 hmss(1:nst,1:nst) = m_zero
2843 call zgemm( 'n', &
2844 'c', &
2845 nst, &
2846 nst, &
2847 i8_to_i4(gr%np_global*st%d%dim), &
2848 cmplx(gr%volume_element, m_zero, real64) , &
2849 hpsi(1, 1, 1), &
2850 ubound(hpsi, dim = 1), &
2851 psi(1, 1, 1), &
2852 ubound(psi, dim = 1), &
2853 m_z0, &
2854 hmss(1, 1), &
2855 ubound(hmss, dim = 1))
2856
2857 hmss(1:nst,1:nst) = conjg(hmss(1:nst,1:nst))
2858
2859 ! accumulate the Floqeut integrals
2860 do in = -forder, forder
2861 do im = -forder, forder
2862 ii = (in+forder) * nst
2863 jj = (im+forder) * nst
2864 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) = &
2865 hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) + hmss(1:nst, 1:nst) * exp(-(in-im)*m_zi*omega*it*dt)
2866 ! diagonal term
2867 if (in == im) then
2868 do ist = 1, nst
2869 hfloquet(ik_count, ii+ist, ii+ist) = hfloquet(ik_count, ii+ist, ii+ist) + in*omega
2870 end do
2871 end if
2872 end do
2873 end do
2874 end do !ik
2875
2876 end do ! it
2877
2878 hfloquet(:,:,:) = m_one/nt*hfloquet(:,:,:)
2879
2880 ! diagonalize Floquet Hamiltonian
2881 if (downfolding) then
2882 ! here perform downfolding
2883 safe_allocate(hfloq_eff(1:nst,1:nst))
2884 safe_allocate(eigenval(1:nst))
2885 safe_allocate(bands(1:nik,1:nst))
2886
2887 hfloq_eff(1:nst,1:nst) = m_zero
2888 do ik = 1, nik
2889 ! the HFloquet blocks are copied directly out of the super matrix
2890 m0 = nst ! the m=0 start position
2891 n0 = nst ! the n=0 start postion
2892 n1 = 2*nst ! the n=+1 start postion
2893 hfloq_eff(1:nst, 1:nst) = hfloquet(ik, n0+1:n0+nst, m0+1:m0+nst) + &
2894 m_one/omega*(matmul(hfloquet(ik, 1:nst, m0+1:m0+nst), hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst))- &
2895 matmul(hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst), hfloquet(ik, 1:nst, m0+1:m0+nst)))
2896
2897 call lalg_eigensolve(nst, hfloq_eff, eigenval)
2898 bands(ik,1:nst) = eigenval(1:nst)
2899 end do
2900 safe_deallocate_a(hfloq_eff)
2901 else
2902 ! the full Floquet
2903 safe_allocate(eigenval(1:nst*fdim))
2904 safe_allocate(bands(1:nik,1:nst*fdim))
2905 safe_allocate(temp(1:nst*fdim, 1:nst*fdim))
2906
2907 do ik = 1, nik
2908 temp(1:nst*fdim, 1:nst*fdim) = hfloquet(ik, 1:nst*fdim, 1:nst*fdim)
2909 call lalg_eigensolve(nst*fdim, temp, eigenval)
2910 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2911 end do
2912 end if
2913
2914 !write bandstructure to file
2915 if (downfolding) then
2916 lim_nst = nst
2917 filename = "downfolded_floquet_bands"
2918 else
2919 lim_nst = nst*fdim
2920 filename = "floquet_bands"
2921 end if
2922 ! write bands (full or downfolded)
2923 if (mpi_world%is_root()) then
2924 file = 987254
2925 file = io_open(filename, namespace, action = 'write')
2926 do ik = 1,nik
2927 do ist = 1,lim_nst
2928 write(file,'(e13.6, 1x)',advance='no') bands(ik,ist)
2929 end do
2930 write(file,'(1x)')
2931 end do
2932 call io_close(file)
2933 end if
2934
2935 if (.not. downfolding) then
2936 ! for the full Floquet case compute also the trivially shifted
2937 ! Floquet bands for reference (i.e. setting H_{nm}=0 for n!=m)
2938 bands(1:nik, 1:nst*fdim) = m_zero
2939 do ik = 1, nik
2940 temp(1:nst*fdim,1:nst*fdim) = m_zero
2941 do jj = 0, fdim - 1
2942 ii = jj * nst
2943 temp(ii+1:ii+nst, ii+1:ii+nst) = hfloquet(ik, ii+1:ii+nst, ii+1:ii+nst)
2944 end do
2945 call lalg_eigensolve(nst*fdim, temp, eigenval)
2946 bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
2947 end do
2948
2949 if (mpi_world%is_root()) then
2950 filename = 'trivial_floquet_bands'
2951 file = io_open(filename, namespace, action = 'write')
2952 do ik=1,nik
2953 do ist = 1,lim_nst
2954 write(file,'(e13.6, 1x)', advance='no') bands(ik,ist)
2955 end do
2956 write(file,'(1x)')
2957 end do
2958 call io_close(file)
2959 end if
2960 end if
2961
2962 ! reset time in Hamiltonian
2963 call hm%update(gr, namespace, space, ext_partners, time=m_zero)
2964
2965 safe_deallocate_a(hmss)
2966 safe_deallocate_a(psi)
2967 safe_deallocate_a(hpsi)
2968 safe_deallocate_a(temp_state1)
2969 safe_deallocate_a(hfloquet)
2970 safe_deallocate_a(eigenval)
2971 safe_deallocate_a(bands)
2972 safe_deallocate_a(temp)
2973
2974 pop_sub(td_write_floquet)
2975
2976 end subroutine td_write_floquet
2977
2978 ! ---------------------------------------------------------
2979 subroutine td_write_total_current(out_total_current, space, mesh, st, iter)
2980 type(c_ptr), intent(inout) :: out_total_current
2981 class(space_t), intent(in) :: space
2982 class(mesh_t), intent(in) :: mesh
2983 type(states_elec_t), intent(in) :: st
2984 integer, intent(in) :: iter
2985
2986 integer :: idir, ispin
2987 character(len=50) :: aux
2988 real(real64) :: total_current(space%dim), abs_current(space%dim)
2989
2990 push_sub(td_write_total_current)
2991
2992 if (mpi_world%is_root() .and. iter == 0) then
2993 call td_write_print_header_init(out_total_current)
2994
2995 ! first line: column names
2996 call write_iter_header_start(out_total_current)
2997
2998 do idir = 1, space%dim
2999 write(aux, '(a2,a1,a1)') 'I(', index2axis(idir), ')'
3000 call write_iter_header(out_total_current, aux)
3001 end do
3002
3003 do idir = 1, space%dim
3004 write(aux, '(a10,a1,a1)') 'IntAbs(j)(', index2axis(idir), ')'
3005 call write_iter_header(out_total_current, aux)
3006 end do
3007
3008 do ispin = 1, st%d%nspin
3009 do idir = 1, space%dim
3010 write(aux, '(a4,i1,a1,a1,a1)') 'I-sp', ispin, '(', index2axis(idir), ')'
3011 call write_iter_header(out_total_current, aux)
3012 end do
3013 end do
3014
3015 call write_iter_nl(out_total_current)
3016
3017 call td_write_print_header_end(out_total_current)
3018 end if
3019
3020 assert(allocated(st%current))
3021
3022 if (mpi_world%is_root()) then
3023 call write_iter_start(out_total_current)
3024 end if
3025
3026 total_current = 0.0_real64
3027 do idir = 1, space%dim
3028 do ispin = 1, st%d%spin_channels
3029 total_current(idir) = total_current(idir) + dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.)
3030 end do
3031 total_current(idir) = units_from_atomic(units_out%length/units_out%time, total_current(idir))
3032 end do
3033 call mesh%allreduce(total_current, dim = space%dim)
3034
3035 abs_current = 0.0_real64
3036 do idir = 1, space%dim
3037 do ispin = 1, st%d%spin_channels
3038 abs_current(idir) = abs_current(idir) + dmf_integrate(mesh, abs(st%current(:, idir, ispin)), reduce = .false.)
3039 end do
3040 abs_current(idir) = units_from_atomic(units_out%length/units_out%time, abs_current(idir))
3041 end do
3042 call mesh%allreduce(abs_current, dim = space%dim)
3043
3044 if (mpi_world%is_root()) then
3045 call write_iter_double(out_total_current, total_current, space%dim)
3046 call write_iter_double(out_total_current, abs_current, space%dim)
3047 end if
3048
3049 do ispin = 1, st%d%nspin
3050 total_current = units_from_atomic(units_out%length/units_out%time, dmf_integrate(mesh, space%dim, st%current(:, :, ispin)))
3051
3052 if (mpi_world%is_root()) then
3053 call write_iter_double(out_total_current, total_current, space%dim)
3054 end if
3055 end do
3056
3057 if (mpi_world%is_root()) then
3058 call write_iter_nl(out_total_current)
3059 end if
3060
3061 pop_sub(td_write_total_current)
3062 end subroutine td_write_total_current
3063
3064 ! ---------------------------------------------------------
3065
3066 subroutine td_write_total_heat_current(write_obj, space, hm, gr, st, iter)
3067 type(c_ptr), intent(inout) :: write_obj
3068 class(space_t), intent(in) :: space
3069 type(hamiltonian_elec_t), intent(inout) :: hm
3070 type(grid_t), intent(in) :: gr
3071 type(states_elec_t), intent(in) :: st
3072 integer, intent(in) :: iter
3073
3074 integer :: idir, ispin
3075 character(len=50) :: aux
3076 real(real64), allocatable :: heat_current(:, :, :)
3077 real(real64) :: total_current(space%dim)
3078
3079 push_sub(td_write_total_current)
3080
3081 if (mpi_world%is_root() .and. iter == 0) then
3082 call td_write_print_header_init(write_obj)
3083
3084 ! first line: column names
3085 call write_iter_header_start(write_obj)
3086
3087 do idir = 1, space%dim
3088 write(aux, '(a2,i1,a1)') 'Jh(', idir, ')'
3089 call write_iter_header(write_obj, aux)
3090 end do
3091
3092 call write_iter_nl(write_obj)
3093
3094 call td_write_print_header_end(write_obj)
3095 end if
3096
3097 safe_allocate(heat_current(1:gr%np, 1:space%dim, 1:st%d%nspin))
3098
3099 call current_heat_calculate(space, gr%der, hm, st, heat_current)
3100
3101 if (mpi_world%is_root()) call write_iter_start(write_obj)
3102
3103 total_current = 0.0_real64
3104 do idir = 1, space%dim
3105 do ispin = 1, st%d%spin_channels
3106 total_current(idir) = total_current(idir) + dmf_integrate(gr, heat_current(:, idir, ispin))
3107 end do
3108 total_current(idir) = units_from_atomic(units_out%energy*units_out%length/units_out%time, total_current(idir))
3109 end do
3110
3111 safe_deallocate_a(heat_current)
3112
3113 if (mpi_world%is_root()) call write_iter_double(write_obj, total_current, space%dim)
3114
3115 if (mpi_world%is_root()) call write_iter_nl(write_obj)
3116
3117 pop_sub(td_write_total_current)
3118 end subroutine td_write_total_heat_current
3119
3120
3121 ! ---------------------------------------------------------
3122 subroutine td_write_partial_charges(out_partial_charges, mesh, st, ions, iter)
3123 type(c_ptr), intent(inout) :: out_partial_charges
3124 class(mesh_t), intent(in) :: mesh
3125 type(states_elec_t), intent(in) :: st
3126 type(ions_t), intent(in) :: ions
3127 integer, intent(in) :: iter
3128
3129 integer :: idir
3130 character(len=50) :: aux
3131 real(real64), allocatable :: hirshfeld_charges(:)
3132
3133 push_sub(td_write_partial_charges)
3134
3135 safe_allocate(hirshfeld_charges(1:ions%natoms))
3136
3137 call partial_charges_calculate(mesh, st, ions, hirshfeld_charges)
3138
3139 if (mpi_world%is_root()) then
3140
3141 if (iter == 0) then
3142
3143 call td_write_print_header_init(out_partial_charges)
3144
3145 ! first line: column names
3146 call write_iter_header_start(out_partial_charges)
3147
3148 do idir = 1, ions%natoms
3149 write(aux, '(a13,i3,a1)') 'hirshfeld(atom=', idir, ')'
3150 call write_iter_header(out_partial_charges, aux)
3151 end do
3152
3153 call write_iter_nl(out_partial_charges)
3154
3155 call td_write_print_header_end(out_partial_charges)
3156 end if
3157
3158 call write_iter_start(out_partial_charges)
3159
3160 call write_iter_double(out_partial_charges, hirshfeld_charges, ions%natoms)
3162 call write_iter_nl(out_partial_charges)
3163 end if
3164
3165 safe_deallocate_a(hirshfeld_charges)
3166
3168 end subroutine td_write_partial_charges
3169
3170 ! ---------------------------------------------------------
3171 subroutine td_write_q(out_q, space, ks, iter)
3172 type(c_ptr), intent(inout) :: out_q
3173 class(space_t), intent(in) :: space
3174 type(v_ks_t), intent(in) :: ks
3175 integer, intent(in) :: iter
3176
3177 integer :: ii
3178 character(len=50) :: aux
3179
3180 push_sub(td_write_q)
3181
3182 if (mpi_world%is_root()) then
3183 if (iter == 0) then
3184 call td_write_print_header_init(out_q)
3185 call write_iter_header_start(out_q)
3186 do ii = 1, ks%pt%nmodes
3187 write(aux, '(a1,i3,a3)') 'q', ii, '(t)'
3188 call write_iter_header(out_q, aux)
3189 end do
3190 do ii = 1, ks%pt%nmodes
3191 write(aux, '(a1,i3,a3)') 'p', ii, '(t)'
3192 call write_iter_header(out_q, aux)
3193 end do
3194 do ii = 1, space%dim
3195 write(aux, '(a3,i3,a3)') 'f_pt', ii, '(t)'
3196 call write_iter_header(out_q, aux)
3197 end do
3198 call write_iter_nl(out_q)
3199 call td_write_print_header_end(out_q)
3200 end if
3201
3202 call write_iter_start(out_q)
3203 call write_iter_double(out_q, ks%pt_mx%pt_q, ks%pt%nmodes)
3204 call write_iter_double(out_q, ks%pt_mx%pt_p, ks%pt%nmodes)
3205 call write_iter_double(out_q, ks%pt_mx%fmf, space%dim)
3206 call write_iter_nl(out_q)
3207 end if
3208
3209 pop_sub(td_write_q)
3210 end subroutine td_write_q
3211
3212
3213 ! ---------------------------------------------------------
3214 subroutine td_write_mxll_field(out_mxll, space, hm, dt, iter)
3215 type(c_ptr), intent(inout) :: out_mxll
3216 class(space_t), intent(in) :: space
3217 type(hamiltonian_elec_t),intent(in) :: hm
3218 real(real64), intent(in) :: dt
3219 integer, intent(in) :: iter
3220
3221 integer :: idir
3222 real(real64) :: field(space%dim)
3223 character(len=80) :: aux
3224 character(len=1) :: field_char
3225
3226 push_sub(td_write_mxll_field)
3227
3228 if (.not. mpi_world%is_root()) then
3229 pop_sub(td_write_mxll_field)
3230 return ! only first node outputs
3231 end if
3232
3233 if (iter == 0) then
3234 call td_write_print_header_init(out_mxll)
3235
3236 write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, dt), &
3237 " [", trim(units_abbrev(units_out%time)), "]"
3238 call write_iter_string(out_mxll, aux)
3239 call write_iter_nl(out_mxll)
3240
3241 call write_iter_header_start(out_mxll)
3242 select case (hm%mxll%coupling_mode)
3243 case (length_gauge_dipole, multipolar_expansion)
3244 if (hm%mxll%add_electric_dip) field_char = 'E'
3245 if (hm%mxll%add_magnetic_dip) field_char = 'B'
3246 do idir = 1, space%dim
3247 write(aux, '(a,i1,a)') field_char // '(', idir, ')'
3248 call write_iter_header(out_mxll, aux)
3249 end do
3250 case (velocity_gauge_dipole)
3251 do idir = 1, space%dim
3252 write(aux, '(a,i1,a)') 'A(', idir, ')'
3253 call write_iter_header(out_mxll, aux)
3254 end do
3255 end select
3256 call write_iter_nl(out_mxll)
3257
3258 call write_iter_string(out_mxll, '#[Iter n.]')
3259 call write_iter_header(out_mxll, '[' // trim(units_abbrev(units_out%time)) // ']')
3260
3261 ! Note that we do not print out units of E, A or B, but rather units of e*E, e*A or e*B
3262 ! (force, energy and magnetic flux density, respectively).
3263 select case (hm%mxll%coupling_mode)
3264 case (length_gauge_dipole, multipolar_expansion)
3265 if (hm%mxll%add_electric_dip) aux = '[' // trim(units_abbrev(units_out%force)) // ']'
3266 if (hm%mxll%add_magnetic_dip) aux = '[' // trim(units_abbrev(unit_one/units_out%length**2)) // ']'
3267 do idir = 1, space%dim
3268 call write_iter_header(out_mxll, aux)
3269 end do
3270 case (velocity_gauge_dipole)
3271 aux = '[' // trim(units_abbrev(units_out%energy)) // ']'
3272 do idir = 1, space%dim
3273 call write_iter_header(out_mxll, aux)
3274 end do
3275 end select
3276 call write_iter_nl(out_mxll)
3277 call td_write_print_header_end(out_mxll)
3278 end if
3279
3280 call write_iter_start(out_mxll)
3281
3282 field = m_zero
3283 select case (hm%mxll%coupling_mode)
3284 case (length_gauge_dipole, multipolar_expansion)
3285 if (hm%mxll%add_electric_dip) field = units_from_atomic(units_out%force, hm%mxll%e_field_dip)
3286 if (hm%mxll%add_magnetic_dip) field = units_from_atomic(unit_one/units_out%length**2, hm%mxll%b_field_dip)
3287 call write_iter_double(out_mxll, field, space%dim)
3288 case (velocity_gauge_dipole)
3289 field = units_from_atomic(units_out%energy, hm%mxll%vec_pot_dip)
3290 call write_iter_double(out_mxll, field, space%dim)
3291 end select
3292 call write_iter_nl(out_mxll)
3293
3294 pop_sub(td_write_mxll_field)
3295 end subroutine td_write_mxll_field
3296
3297
3298 ! ---------------------------------------------------------
3299 subroutine td_write_effective_u(out_coords, lda_u, iter)
3300 type(c_ptr), intent(inout) :: out_coords
3301 type(lda_u_t), intent(in) :: lda_u
3302 integer, intent(in) :: iter
3303
3304 integer :: ios, inn
3305 character(len=50) :: aux
3306
3307 if (.not. mpi_world%is_root()) return ! only first node outputs
3308
3310
3311 if (iter == 0) then
3312 call td_write_print_header_init(out_coords)
3313
3314 ! first line: column names
3315 call write_iter_header_start(out_coords)
3316
3317 do ios = 1, lda_u%norbsets
3318 write(aux, '(a2,i3,a1)') 'Ueff(', ios, ')'
3319 call write_iter_header(out_coords, aux)
3320 end do
3321
3322 do ios = 1, lda_u%norbsets
3323 write(aux, '(a2,i3,a1)') 'U(', ios, ')'
3324 call write_iter_header(out_coords, aux)
3325 end do
3326
3327 do ios = 1, lda_u%norbsets
3328 write(aux, '(a2,i3,a1)') 'J(', ios, ')'
3329 call write_iter_header(out_coords, aux)
3330 end do
3331
3332 if (lda_u%intersite) then
3333 do ios = 1, lda_u%norbsets
3334 do inn = 1, lda_u%orbsets(ios)%nneighbors
3335 write(aux, '(a2,i3,a1,i3,a1)') 'V(', ios,'-', inn, ')'
3336 call write_iter_header(out_coords, aux)
3337 end do
3338 end do
3339 end if
3340
3341
3342 call write_iter_nl(out_coords)
3343
3344 ! second line: units
3345 call write_iter_string(out_coords, '#[Iter n.]')
3346 call write_iter_header(out_coords, '[' // trim(units_abbrev(units_out%time)) // ']')
3347 call write_iter_string(out_coords, &
3348 'Effective U ' // trim(units_abbrev(units_out%energy)) // &
3349 ', U in '// trim(units_abbrev(units_out%energy)) // &
3350 ', J in ' // trim(units_abbrev(units_out%energy)))
3351 call write_iter_nl(out_coords)
3352
3353 call td_write_print_header_end(out_coords)
3354 end if
3355
3356 call write_iter_start(out_coords)
3357
3358 do ios = 1, lda_u%norbsets
3359 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3360 lda_u%orbsets(ios)%Ueff), 1)
3361 end do
3362
3363 do ios = 1, lda_u%norbsets
3364 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3365 lda_u%orbsets(ios)%Ubar), 1)
3366 end do
3367
3368 do ios = 1, lda_u%norbsets
3369 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3370 lda_u%orbsets(ios)%Jbar), 1)
3371 end do
3372
3373 if (lda_u%intersite) then
3374 do ios = 1, lda_u%norbsets
3375 do inn = 1, lda_u%orbsets(ios)%nneighbors
3376 call write_iter_double(out_coords, units_from_atomic(units_out%energy, &
3377 lda_u%orbsets(ios)%V_ij(inn,0)), 1)
3378 end do
3379 end do
3380 end if
3381
3382 call write_iter_nl(out_coords)
3383
3384 pop_sub(td_write_effective_u)
3385 end subroutine td_write_effective_u
3386
3388 subroutine td_write_norm_ks_orbitals(file_handle, grid, kpoints, st, iter)
3389 type(c_ptr), intent(inout) :: file_handle
3390 type(grid_t), intent(in) :: grid
3391 type(kpoints_t), intent(in) :: kpoints
3392 type(states_elec_t), intent(in) :: st
3393 integer, intent(in) :: iter
3395 integer :: ik_ispin, ist
3396 character(len=7) :: nkpt_str, nst_str
3397 character(len=7) :: ik_str, ist_str
3398 real(real64), allocatable :: norm_ks(:, :)
3399 real(real64) :: n_electrons
3400
3402
3403 safe_allocate(norm_ks(1:st%nst, 1:st%nik))
3404 call states_elec_calc_norms(grid, kpoints, st, norm_ks)
3405
3406 if (mpi_world%is_root()) then
3407 ! Header
3408 if (iter == 0) then
3409 call td_write_print_header_init(file_handle)
3410
3411 ! Line 1
3412 write(nkpt_str, '(I7)') st%nik
3413 write(nst_str, '(I7)') st%nst
3414 call write_iter_string(file_handle, '# Dimensions. (nstates, nkpt * nspin):')
3415 call write_iter_string(file_handle, trim(adjustl(nst_str)) // ' ' // trim(adjustl(nkpt_str)))
3416 call write_iter_nl(file_handle)
3417
3418 ! Line 2
3419 call write_iter_string(file_handle, '# Norm ordering: (istate, ikpoint_spin)')
3420 call write_iter_nl(file_handle)
3421
3422 ! Line 3
3423 call write_iter_header_start(file_handle)
3424 call write_iter_header(file_handle, 'N_electrons')
3425 do ik_ispin = 1, st%nik
3426 do ist = 1, st%nst
3427 write(ik_str, '(I7)') ik_ispin
3428 write(ist_str, '(I7)') ist
3429 call write_iter_header(file_handle, &
3430 'Norm (' // trim(ist_str) // ',' // trim(ik_str) // ')')
3431 enddo
3432 enddo
3433 call write_iter_nl(file_handle)
3434 call td_write_print_header_end(file_handle)
3435 end if
3436
3437 n_electrons = sum(st%occ * norm_ks**2)
3438
3439 ! Output norms for time step `iter`
3440 call write_iter_start(file_handle)
3441 call write_iter_double(file_handle, n_electrons, 1)
3442 do ik_ispin = 1, st%nik
3443 call write_iter_double(file_handle, norm_ks(:, ik_ispin), size(norm_ks, 1))
3444 enddo
3445 call write_iter_nl(file_handle)
3446
3447 end if
3448
3449 safe_deallocate_a(norm_ks)
3450
3452
3453 end subroutine td_write_norm_ks_orbitals
3454
3455 ! ---------------------------------------------------------
3457 subroutine td_write_cell_parameters(file_handle, ions, iter)
3458 type(c_ptr), intent(inout) :: file_handle
3459 type(ions_t), intent(in) :: ions
3460 integer, intent(in) :: iter
3461
3462 integer :: idir
3463 real(real64) :: tmp(3)
3464
3465 if (.not. mpi_world%is_root()) return ! only first task outputs
3466
3467 push_sub(td_write_cell_parameters)
3468
3469 assert(ions%space%dim == 3)
3470
3471 if (iter == 0) then
3472 call td_write_print_header_init(file_handle)
3473
3474 ! first line: column names
3475 call write_iter_header_start(file_handle)
3476
3477 call write_iter_string(file_handle, '# Iter, a, b, c, volume, alpha, beta, gamma, ' &
3478 // 'a_x, a_y, a_z, b_x, b_y, b_z, c_x, c_y, c_z')
3479
3480 ! second line: units
3481 call write_iter_string(file_handle, '#[Iter n.]')
3482 call write_iter_header(file_handle, '[' // trim(units_abbrev(units_out%time)) // ']')
3483 call write_iter_string(file_handle, &
3484 'Lengths in ' // trim(units_abbrev(units_out%length)) // &
3485 ', Volume in ' // trim(units_abbrev(units_out%length**3)) // &
3486 ', Angles in degree, Lattice vectors in '// trim(units_abbrev(units_out%length)))
3487 call write_iter_nl(file_handle)
3488
3489 call td_write_print_header_end(file_handle)
3490 end if
3491
3492 call write_iter_start(file_handle)
3493
3494 ! Length of the lattice vectors
3495 do idir = 1, 3
3496 tmp(idir) = units_from_atomic(units_out%length, norm2(ions%latt%rlattice(1:3, idir)))
3497 end do
3498 call write_iter_double(file_handle, tmp, 3)
3499
3500 ! Cell volume
3501 tmp(1) = units_from_atomic(units_out%length**3, ions%latt%rcell_volume)
3502 call write_iter_double(file_handle, tmp(1), 1)
3503
3504 ! Lattice angles
3505 call write_iter_double(file_handle, ions%latt%alpha, 1)
3506 call write_iter_double(file_handle, ions%latt%beta, 1)
3507 call write_iter_double(file_handle, ions%latt%gamma, 1)
3508
3509 !Lattice vectors
3510 do idir = 1, 3
3511 tmp(1:3) = units_from_atomic(units_out%length, ions%latt%rlattice(:, idir))
3512 call write_iter_double(file_handle, tmp, 3)
3513 end do
3514 call write_iter_nl(file_handle)
3515
3517 end subroutine td_write_cell_parameters
3518
3519
3520 ! ---------------------------------------------------------
3521 subroutine td_write_mxll_init(writ, namespace, iter, dt)
3522 type(td_write_t), intent(out) :: writ
3523 type(namespace_t), intent(in) :: namespace
3524 integer, intent(in) :: iter
3525 real(real64), intent(in) :: dt
3526
3527 integer :: default, flags, iout, first
3528
3529 push_sub(td_write_mxll_init)
3530
3531 !%Variable MaxwellTDOutput
3532 !%Type flag
3533 !%Default maxwell_energy
3534 !%Section Maxwell::Output
3535 !%Description
3536 !% Defines what should be output during the time-dependent
3537 !% Maxwell simulation. Many of the options can increase the computational
3538 !% cost of the simulation, so only use the ones that you need. In
3539 !% most cases the default value is enough, as it is adapted to the
3540 !% details of the TD run.
3541 !% WARNING: the calculation of the longitudinal or transverse E and B fields
3542 !% can be very expensive, so please consider using the MaxwellOutput block
3543 !% to calculate and output these quantities at certain timesteps.
3544 !%Option maxwell_total_e_field 1
3545 !% Output of the total (longitudinal plus transverse) electric field at
3546 !% the points specified in the MaxwellFieldsCoordinate block
3547 !%Option maxwell_total_b_field 8
3548 !% Output of the total (longitudinal plus transverse) magnetic field at
3549 !% the points specified in the MaxwellFieldsCoordinate block
3550 !%Option maxwell_longitudinal_e_field 64
3551 !% Output of the longitudinal electric field at the points
3552 !% specified in the MaxwellFieldsCoordinate block (can slow down the run)
3553 !%Option maxwell_longitudinal_b_field 512
3554 !% Output of the longitudinal magnetic field at the points
3555 !% specified in the MaxwellFieldsCoordinate block (can slow down the run)
3556 !%Option maxwell_transverse_e_field 4096
3557 !% Output of the transverse electric field at the points
3558 !% specified in the MaxwellFieldsCoordinate block (can slow down the run)
3559 !%Option maxwell_transverse_b_field 32768
3560 !% Output of the transverse magnetic field at the points
3561 !% specified in the MaxwellFieldsCoordinate block (can slow down the run)
3562 !%Option maxwell_energy 262144
3563 !% Output of the electromagnetic field energy into the folder <tt>td.general/maxwell</tt>.
3564 !% WARNING: the transverse and longitudinal energies will be correct only if you request
3565 !% the longitudinal and transverse E or B fields as output. Otherwise they will be set to
3566 !% zero.
3567 !%Option e_field_surface_x 524288
3568 !% Output of the E field sliced along the plane x=0 for each field component
3569 !%Option e_field_surface_y 1048576
3570 !% Output of the E field sliced along the plane y=0 for each field component
3571 !%Option e_field_surface_z 2097152
3572 !% Output of the E field sliced along the plane z=0 for each field component
3573 !%Option b_field_surface_x 4194304
3574 !% Output of the B field sliced along the plane x=0 for each field component
3575 !%Option b_field_surface_y 8388608
3576 !% Output of the B field sliced along the plane y=0 for each field component
3577 !%Option b_field_surface_z 16777216
3578 !% Output of the B field sliced along the plane z=0 for each field component
3579 !%End
3580
3581 default = 2**(out_maxwell_energy - 1)
3582 call parse_variable(namespace, 'MaxwellTDOutput', default, flags)
3583
3584 if (.not. varinfo_valid_option('MaxwellTDOutput', flags, is_flag = .true.)) then
3585 call messages_input_error(namespace, 'MaxwellTDOutput')
3586 end if
3587
3589 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3590 if (writ%out(iout)%write) then
3591 writ%out(iout + 1)%write = .true.
3592 writ%out(iout + 2)%write = .true.
3593 end if
3594 end do
3595
3597 writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
3598 end do
3599
3600 if (iter == 0) then
3601 first = 0
3602 else
3603 first = iter + 1
3604 end if
3605
3606 call io_mkdir('td.general', namespace)
3607
3608 ! total E field
3609 if (writ%out(out_maxwell_total_e_field)%write) then
3610 call write_iter_init(writ%out(out_maxwell_total_e_field)%handle, first, &
3611 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_x", namespace)))
3612 call write_iter_init(writ%out(out_maxwell_total_e_field + 1)%handle, first, &
3613 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_y", namespace)))
3614 call write_iter_init(writ%out(out_maxwell_total_e_field + 2)%handle, first, &
3615 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_z", namespace)))
3616 end if
3617
3618 ! total B field
3619 if (writ%out(out_maxwell_total_b_field)%write) then
3620 call write_iter_init(writ%out(out_maxwell_total_b_field)%handle, first, &
3621 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_x", namespace)))
3622 call write_iter_init(writ%out(out_maxwell_total_b_field + 1)%handle, first, &
3623 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_y", namespace)))
3624 call write_iter_init(writ%out(out_maxwell_total_b_field + 2)%handle, first, &
3625 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_z", namespace)))
3626 end if
3627
3628 ! longitudinal E field
3629 if (writ%out(out_maxwell_long_e_field)%write) then
3630 call write_iter_init(writ%out(out_maxwell_long_e_field)%handle, first, &
3631 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_x", namespace)))
3632 call write_iter_init(writ%out(out_maxwell_long_e_field + 1)%handle, first, &
3633 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_y", namespace)))
3634 call write_iter_init(writ%out(out_maxwell_long_e_field + 2)%handle, first, &
3635 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_z", namespace)))
3636 end if
3637
3638 ! longitudinal B field
3639 if (writ%out(out_maxwell_long_b_field)%write) then
3640 call write_iter_init(writ%out(out_maxwell_long_b_field)%handle, first, &
3641 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_x", namespace)))
3642 call write_iter_init(writ%out(out_maxwell_long_b_field + 1)%handle, first, &
3643 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_y", namespace)))
3644 call write_iter_init(writ%out(out_maxwell_long_b_field + 2)%handle, first, &
3645 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_z", namespace)))
3646 end if
3647
3648 ! transverse E field
3649 if (writ%out(out_maxwell_trans_e_field)%write) then
3650 call write_iter_init(writ%out(out_maxwell_trans_e_field)%handle, first, &
3651 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_x", namespace)))
3652 call write_iter_init(writ%out(out_maxwell_trans_e_field + 1)%handle, first, &
3653 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_y", namespace)))
3654 call write_iter_init(writ%out(out_maxwell_trans_e_field + 2)%handle, first, &
3655 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_z", namespace)))
3656 end if
3657
3658 ! transverse B field
3659 if (writ%out(out_maxwell_trans_b_field)%write) then
3660 call write_iter_init(writ%out(out_maxwell_trans_b_field)%handle, first, &
3661 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_x", namespace)))
3662 call write_iter_init(writ%out(out_maxwell_trans_b_field + 1)%handle, first, &
3663 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_y", namespace)))
3664 call write_iter_init(writ%out(out_maxwell_trans_b_field + 2)%handle, first, &
3665 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_z", namespace)))
3666 end if
3667
3668 ! Maxwell energy
3669 if (writ%out(out_maxwell_energy)%write) then
3670 call write_iter_init(writ%out(out_maxwell_energy)%handle, first, &
3671 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/maxwell_energy", namespace)))
3672 end if
3673
3674 if (writ%out(out_e_field_surface_x)%write) then
3675 call write_iter_init(writ%out(out_e_field_surface_x)%handle, first, &
3676 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-x", namespace)))
3677 end if
3678
3679 if (writ%out(out_e_field_surface_y)%write) then
3680 call write_iter_init(writ%out(out_e_field_surface_y)%handle, first, &
3681 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-y", namespace)))
3682 end if
3683
3684 if (writ%out(out_e_field_surface_z)%write) then
3685 call write_iter_init(writ%out(out_e_field_surface_z)%handle, first, &
3686 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-z", namespace)))
3687 end if
3688
3689 if (writ%out(out_b_field_surface_x)%write) then
3690 call write_iter_init(writ%out(out_b_field_surface_x)%handle, first, &
3691 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-x", namespace)))
3692 end if
3693
3694 if (writ%out(out_b_field_surface_y)%write) then
3695 call write_iter_init(writ%out(out_b_field_surface_y)%handle, first, &
3696 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-y", namespace)))
3697 end if
3698
3699 if (writ%out(out_b_field_surface_z)%write) then
3700 call write_iter_init(writ%out(out_b_field_surface_z)%handle, first, &
3701 units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-z", namespace)))
3702 end if
3703
3704 pop_sub(td_write_mxll_init)
3705 end subroutine td_write_mxll_init
3706
3707
3708 ! ---------------------------------------------------------
3709 subroutine td_write_mxll_end(writ)
3710 type(td_write_t), intent(inout) :: writ
3711
3712 integer :: iout
3713
3714 push_sub(td_write_mxll_end)
3715
3716 do iout = 1, out_maxwell_max
3717 if (writ%out(iout)%write) call write_iter_end(writ%out(iout)%handle)
3718 end do
3719
3720 pop_sub(td_write_mxll_end)
3721 end subroutine td_write_mxll_end
3722
3723
3724 ! ---------------------------------------------------------
3725 subroutine td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
3726 type(td_write_t), intent(inout) :: writ
3727 class(space_t), intent(in) :: space
3728 type(grid_t), intent(inout) :: gr
3729 type(states_mxll_t), intent(inout) :: st
3730 type(hamiltonian_mxll_t), intent(inout) :: hm
3731 type(helmholtz_decomposition_t), intent(inout) :: helmholtz
3732 real(real64), intent(in) :: dt
3733 integer, intent(in) :: iter
3734 type(namespace_t), intent(in) :: namespace
3735
3736
3737 push_sub(td_write_mxll_iter)
3738
3739 call profiling_in("TD_WRITE_ITER_MAXWELL")
3740
3741 if (writ%out(out_maxwell_trans_e_field)%write .or. writ%out(out_maxwell_trans_b_field)%write) then
3742 call helmholtz%get_trans_field(namespace, st%rs_state_trans, total_field=st%rs_state)
3743 call get_rs_state_at_point(st%selected_points_rs_state_trans(:,:), st%rs_state_trans, &
3744 st%selected_points_coordinate(:,:), st, gr)
3745 !TODO: calculate transverse energy
3746 else
3747 hm%energy%energy_trans = m_zero
3748 end if
3749
3750 if (writ%out(out_maxwell_long_e_field)%write .or. writ%out(out_maxwell_long_b_field)%write) then
3751 call helmholtz%get_long_field(namespace, st%rs_state_long, total_field=st%rs_state)
3752 call get_rs_state_at_point(st%selected_points_rs_state_long(:,:), st%rs_state_long, &
3753 st%selected_points_coordinate(:,:), st, gr)
3754 !TODO: calculate longitudinal energy
3755 else
3756 hm%energy%energy_long = m_zero
3757 end if
3758
3759 ! total E field
3760 if (writ%out(out_maxwell_total_e_field)%write) then
3761 call td_write_fields(writ%out(out_maxwell_total_e_field)%handle, space, st, iter, dt, &
3763 call td_write_fields(writ%out(out_maxwell_total_e_field + 1)%handle, space, st, iter, dt, &
3765 call td_write_fields(writ%out(out_maxwell_total_e_field + 2)%handle, space, st, iter, dt, &
3767 end if
3768
3769 ! total B field
3770 if (writ%out(out_maxwell_total_b_field)%write) then
3771 call td_write_fields(writ%out(out_maxwell_total_b_field)%handle, space, st, iter, dt, &
3773 call td_write_fields(writ%out(out_maxwell_total_b_field + 1)%handle, space, st, iter, dt, &
3775 call td_write_fields(writ%out(out_maxwell_total_b_field + 2)%handle, space, st, iter, dt, &
3777 end if
3778
3779 ! Longitudinal E field
3780 if (writ%out(out_maxwell_long_e_field)%write) then
3781 call td_write_fields(writ%out(out_maxwell_long_e_field)%handle, space, st, iter, dt, &
3783 call td_write_fields(writ%out(out_maxwell_long_e_field + 1)%handle, space, st, iter, dt, &
3785 call td_write_fields(writ%out(out_maxwell_long_e_field + 2)%handle, space, st, iter, dt, &
3787 end if
3788
3789 ! Longitudinal B field
3790 if (writ%out(out_maxwell_long_b_field)%write) then
3791 call td_write_fields(writ%out(out_maxwell_long_b_field)%handle, space, st, iter, dt, &
3793 call td_write_fields(writ%out(out_maxwell_long_b_field + 1)%handle, space, st, iter, dt, &
3795 call td_write_fields(writ%out(out_maxwell_long_b_field + 2)%handle, space, st, iter, dt, &
3797 end if
3798
3799 ! Transverse E field
3800 if (writ%out(out_maxwell_trans_e_field)%write) then
3801 call td_write_fields(writ%out(out_maxwell_trans_e_field)%handle, space, st, iter, dt, &
3803 call td_write_fields(writ%out(out_maxwell_trans_e_field + 1)%handle, space, st, iter, dt, &
3805 call td_write_fields(writ%out(out_maxwell_trans_e_field + 2)%handle, space, st, iter, dt, &
3807 end if
3808
3809 ! Transverse B field
3810 if (writ%out(out_maxwell_trans_b_field)%write) then
3811 call td_write_fields(writ%out(out_maxwell_trans_b_field)%handle, space, st, iter, dt, &
3813 call td_write_fields(writ%out(out_maxwell_trans_b_field + 1)%handle, space, st, iter, dt, &
3815 call td_write_fields(writ%out(out_maxwell_trans_b_field + 2)%handle, space, st, iter, dt, &
3817 end if
3818
3819 ! Maxwell energy
3820 if (writ%out(out_maxwell_energy)%write) then
3821 call td_write_maxwell_energy(writ%out(out_maxwell_energy)%handle, hm, iter)
3822 end if
3823
3824 if (writ%out(out_e_field_surface_x)%write) then
3825 call td_write_electric_field_box_surface(writ%out(out_e_field_surface_x)%handle, st, 1, iter)
3826 end if
3827
3828 if (writ%out(out_e_field_surface_y)%write) then
3829 call td_write_electric_field_box_surface(writ%out(out_e_field_surface_y)%handle, st, 2, iter)
3830 end if
3831
3832 if (writ%out(out_e_field_surface_z)%write) then
3833 call td_write_electric_field_box_surface(writ%out(out_e_field_surface_z)%handle, st, 3, iter)
3834 end if
3835
3836 if (writ%out(out_b_field_surface_x)%write) then
3837 call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_x)%handle, st, 1, iter)
3838 end if
3839
3840 if (writ%out(out_b_field_surface_y)%write) then
3841 call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_y)%handle, st, 2, iter)
3842 end if
3843
3844 if (writ%out(out_b_field_surface_z)%write) then
3845 call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_z)%handle, st, 3, iter)
3846 end if
3847
3848 call profiling_out("TD_WRITE_ITER_MAXWELL")
3849
3850 pop_sub(td_write_mxll_iter)
3851 end subroutine td_write_mxll_iter
3852
3853
3854 ! ---------------------------------------------------------
3855 subroutine td_write_maxwell_energy(out_maxwell_energy, hm, iter)
3856 type(c_ptr), intent(inout) :: out_maxwell_energy
3857 type(hamiltonian_mxll_t), intent(in) :: hm
3858 integer, intent(in) :: iter
3859
3860 integer :: ii
3861
3862 integer :: n_columns
3863
3864 if (.not. mpi_world%is_root()) return ! only first node outputs
3865
3866 push_sub(td_write_maxwell_energy)
3867
3868 n_columns = 8
3869
3870 if (iter == 0) then
3871 call td_write_print_header_init(out_maxwell_energy)
3872
3873 ! first line -> column names
3874 call write_iter_header_start(out_maxwell_energy)
3875 call write_iter_header(out_maxwell_energy, 'Total')
3876 call write_iter_header(out_maxwell_energy, 'E**2')
3877 call write_iter_header(out_maxwell_energy, 'B**2')
3878 call write_iter_header(out_maxwell_energy, 'Total+Boundaries')
3879 call write_iter_header(out_maxwell_energy, 'Boundaries')
3880 call write_iter_header(out_maxwell_energy, 'Transversal')
3881 call write_iter_header(out_maxwell_energy, 'Longitudinal')
3882 call write_iter_header(out_maxwell_energy, 'Incident Waves')
3883
3884 call write_iter_nl(out_maxwell_energy)
3885
3886 ! units
3887
3888 call write_iter_string(out_maxwell_energy, '#[Iter n.]')
3889 call write_iter_header(out_maxwell_energy, '[' // trim(units_abbrev(units_out%time)) // ']')
3890
3891 do ii = 1, n_columns
3892 call write_iter_header(out_maxwell_energy, '[' // trim(units_abbrev(units_out%energy)) // ']')
3893 end do
3894 call write_iter_nl(out_maxwell_energy)
3895
3896 call td_write_print_header_end(out_maxwell_energy)
3897 end if
3898
3899 call write_iter_start(out_maxwell_energy)
3900 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy), 1)
3901 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%e_energy), 1)
3902 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%b_energy), 1)
3903 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, &
3904 hm%energy%energy+hm%energy%boundaries), 1)
3905 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%boundaries), 1)
3906 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_trans), 1)
3907 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_long), 1)
3908 call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_plane_waves), 1)
3909 call write_iter_nl(out_maxwell_energy)
3910
3912 end subroutine td_write_maxwell_energy
3913
3914
3915 ! ---------------------------------------------------------
3916 subroutine td_write_electric_field_box_surface(out_field_surf, st, dim, iter)
3917 type(c_ptr), intent(inout) :: out_field_surf
3918 type(states_mxll_t), intent(in) :: st
3919 integer, intent(in) :: dim
3920 integer, intent(in) :: iter
3921
3922 integer :: ii
3923
3924 integer :: n_columns
3925
3926 if (.not. mpi_world%is_root()) return ! only first node outputs
3927
3929
3930 n_columns = 12
3931
3932 if (iter == 0) then
3933 call td_write_print_header_init(out_field_surf)
3934
3935 ! first line -> column names
3936 call write_iter_header_start(out_field_surf)
3937 call write_iter_header(out_field_surf, '- x direction')
3938 call write_iter_header(out_field_surf, '+ x direction')
3939 call write_iter_header(out_field_surf, '- y direction')
3940 call write_iter_header(out_field_surf, '+ y direction')
3941 call write_iter_header(out_field_surf, '- z direction')
3942 call write_iter_header(out_field_surf, '+ z direction')
3943 call write_iter_header(out_field_surf, '- x dir. p. w.')
3944 call write_iter_header(out_field_surf, '+ x dir. p. w.')
3945 call write_iter_header(out_field_surf, '- y dir. p. w.')
3946 call write_iter_header(out_field_surf, '+ y dir. p. w.')
3947 call write_iter_header(out_field_surf, '- z dir. p. w.')
3948 call write_iter_header(out_field_surf, '+ z dir. p. w.')
3949
3950 call write_iter_nl(out_field_surf)
3951
3952 ! units
3953 call write_iter_string(out_field_surf, '#[Iter n.]')
3954 call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%time)) // ']')
3955
3956 do ii = 1, n_columns
3957 call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%energy/units_out%length)) // ']')
3958 end do
3959 call write_iter_nl(out_field_surf)
3960
3961 call td_write_print_header_end(out_field_surf)
3962 end if
3963
3964 call write_iter_start(out_field_surf)
3965 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3966 st%electric_field_box_surface(1,1,dim)), 1)
3967 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3968 st%electric_field_box_surface(2,1,dim)), 1)
3969 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3970 st%electric_field_box_surface(1,2,dim)), 1)
3971 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3972 st%electric_field_box_surface(2,2,dim)), 1)
3973 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3974 st%electric_field_box_surface(1,3,dim)), 1)
3975 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3976 st%electric_field_box_surface(2,3,dim)), 1)
3977 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3978 st%electric_field_box_surface_plane_waves(1,1,dim)), 1)
3979 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3980 st%electric_field_box_surface_plane_waves(2,1,dim)), 1)
3981 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3982 st%electric_field_box_surface_plane_waves(1,2,dim)), 1)
3983 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3984 st%electric_field_box_surface_plane_waves(2,2,dim)), 1)
3985 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3986 st%electric_field_box_surface_plane_waves(1,3,dim)), 1)
3987 call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
3988 st%electric_field_box_surface_plane_waves(2,3,dim)), 1)
3989 call write_iter_nl(out_field_surf)
3990
3993
3994
3995 ! ---------------------------------------------------------
3996 subroutine td_write_magnetic_field_box_surface(out_field_surf, st, dim, iter)
3997 type(c_ptr), intent(inout) :: out_field_surf
3998 type(states_mxll_t), intent(in) :: st
3999 integer, intent(in) :: dim
4000 integer, intent(in) :: iter
4001
4002 integer :: ii
4003
4004 integer :: n_columns
4005
4006 if (.not. mpi_world%is_root()) return ! only first node outputs
4007
4009
4010 n_columns = 12
4012 if (iter == 0) then
4013 call td_write_print_header_init(out_field_surf)
4014
4015 ! first line -> column names
4016 call write_iter_header_start(out_field_surf)
4017 call write_iter_header(out_field_surf, '- x direction')
4018 call write_iter_header(out_field_surf, '+ x direction')
4019 call write_iter_header(out_field_surf, '- y direction')
4020 call write_iter_header(out_field_surf, '+ y direction')
4021 call write_iter_header(out_field_surf, '- z direction')
4022 call write_iter_header(out_field_surf, '+ z direction')
4023 call write_iter_header(out_field_surf, '- x dir. p. w.')
4024 call write_iter_header(out_field_surf, '+ x dir. p. w.')
4025 call write_iter_header(out_field_surf, '- y dir. p. w.')
4026 call write_iter_header(out_field_surf, '+ y dir. p. w.')
4027 call write_iter_header(out_field_surf, '- z dir. p. w.')
4028 call write_iter_header(out_field_surf, '+ z dir. p. w.')
4029
4030 call write_iter_nl(out_field_surf)
4031
4032 ! units
4033 call write_iter_string(out_field_surf, '#[Iter n.]')
4034 call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%time)) // ']')
4035
4036 do ii = 1, n_columns
4037 call write_iter_header(out_field_surf, '[' // trim(units_abbrev(unit_one/units_out%length**2)) // ']')
4038 end do
4039 call write_iter_nl(out_field_surf)
4040
4041 call td_write_print_header_end(out_field_surf)
4042 end if
4043
4044 call write_iter_start(out_field_surf)
4045 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4046 st%magnetic_field_box_surface(1,1,dim)), 1)
4047 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4048 st%magnetic_field_box_surface(2,1,dim)), 1)
4049 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4050 st%magnetic_field_box_surface(1,2,dim)), 1)
4051 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4052 st%magnetic_field_box_surface(2,2,dim)), 1)
4053 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4054 st%magnetic_field_box_surface(1,3,dim)), 1)
4055 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4056 st%magnetic_field_box_surface(2,3,dim)), 1)
4057 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4058 st%magnetic_field_box_surface_plane_waves(1,1,dim)), 1)
4059 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4060 st%magnetic_field_box_surface_plane_waves(2,1,dim)), 1)
4061 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4062 st%magnetic_field_box_surface_plane_waves(1,2,dim)), 1)
4063 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4064 st%magnetic_field_box_surface_plane_waves(2,2,dim)), 1)
4065 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4066 st%magnetic_field_box_surface_plane_waves(1,3,dim)), 1)
4067 call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
4068 st%magnetic_field_box_surface_plane_waves(2,3,dim)), 1)
4069 call write_iter_nl(out_field_surf)
4070
4073
4074 ! ---------------------------------------------------------
4075 subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
4076 type(c_ptr), intent(inout) :: out_fields
4077 class(space_t), intent(in) :: space
4078 type(states_mxll_t), intent(in) :: st
4079 integer, intent(in) :: iter
4080 real(real64), intent(in) :: dt
4081 integer, intent(in) :: e_or_b_field
4082 integer, intent(in) :: field_type
4083 integer, intent(in) :: idir
4084
4085
4086 integer :: id
4087 real(real64) :: field(space%dim), selected_field
4088 character(len=80) :: aux
4089
4090 if (.not. mpi_world%is_root()) return ! only first node outputs
4092 push_sub(td_write_fields)
4093
4094 if (iter == 0) then
4095 call td_write_print_header_init(out_fields)
4096
4097 ! first line
4098 write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, dt), &
4099 " [", trim(units_abbrev(units_out%time)), "]"
4100 call write_iter_string(out_fields, aux)
4101 call write_iter_nl(out_fields)
4102
4103 call write_iter_header_start(out_fields)
4104
4105 do id = 1, st%selected_points_number
4106 select case (e_or_b_field)
4107 case (maxwell_e_field)
4108 write(aux, '(a,i1,a)') 'E(', id, ')'
4109 case (maxwell_b_field)
4110 write(aux, '(a,i1,a)') 'B(', id, ')'
4111 end select
4112 call write_iter_header(out_fields, aux)
4113 end do
4114
4115 call write_iter_nl(out_fields)
4116 call write_iter_string(out_fields, '#[Iter n.]')
4117 call write_iter_header(out_fields, ' [' // trim(units_abbrev(units_out%time)) // ']')
4118
4119 ! Note that we do not print out units of E, B, or A, but rather units of e*E, e*B, e*A.
4120 ! (force, force, and energy, respectively). The reason is that the units of E, B or A
4121 ! are ugly.
4122 aux = ' [' // trim(units_abbrev(units_out%force)) // ']'
4123 do id = 1, st%selected_points_number
4124 call write_iter_header(out_fields, aux)
4125 end do
4126 call write_iter_nl(out_fields)
4127 call td_write_print_header_end(out_fields)
4128 end if
4129
4130 call write_iter_start(out_fields)
4131
4132 do id = 1, st%selected_points_number
4133 select case (e_or_b_field)
4134 case (maxwell_e_field)
4135 ! Output of electric field at selected point
4136 select case (field_type)
4137 case (maxwell_total_field)
4138 call get_electric_field_vector(st%selected_points_rs_state(:,id), field(1:st%dim))
4139 case (maxwell_long_field)
4140 call get_electric_field_vector(st%selected_points_rs_state_long(:,id), field(1:st%dim))
4141 case (maxwell_trans_field)
4142 call get_electric_field_vector(st%selected_points_rs_state_trans(:,id), field(1:st%dim))
4143 end select
4144 selected_field = units_from_atomic(units_out%energy/units_out%length, field(idir))
4145 case (maxwell_b_field)
4146 ! Output of magnetic field at selected point
4147 select case (field_type)
4148 case (maxwell_total_field)
4149 call get_magnetic_field_vector(st%selected_points_rs_state(:,id), st%rs_sign, field(1:st%dim))
4150 case (maxwell_long_field)
4151 call get_magnetic_field_vector(st%selected_points_rs_state_long(:,id), st%rs_sign, field(1:st%dim))
4152 case (maxwell_trans_field)
4153 call get_magnetic_field_vector(st%selected_points_rs_state_trans(:,id), st%rs_sign, field(1:st%dim))
4154 end select
4155 selected_field = units_from_atomic(unit_one/units_out%length**2, field(idir))
4156 end select
4157 call write_iter_double(out_fields, selected_field, 1)
4158 end do
4159
4160 call write_iter_nl(out_fields)
4161
4162 pop_sub(td_write_fields)
4163 end subroutine td_write_fields
4164
4165
4166 !----------------------------------------------------------
4167 subroutine td_write_mxll_free_data(namespace, space, gr, st, hm, helmholtz, outp, iter, time)
4168 type(namespace_t), intent(in) :: namespace
4169 class(space_t), intent(in) :: space
4170 type(grid_t), intent(inout) :: gr
4171 type(states_mxll_t), intent(inout) :: st
4172 type(hamiltonian_mxll_t), intent(inout) :: hm
4173 type(helmholtz_decomposition_t), intent(inout) :: helmholtz
4174 type(output_t), intent(in) :: outp
4175 integer, intent(in) :: iter
4176 real(real64), intent(in) :: time
4177
4178 character(len=256) :: filename
4179
4180 push_sub(td_write_maxwell_free_data)
4181 call profiling_in("TD_WRITE_MAXWELL_DATA")
4182
4183 ! now write down the rest
4184 write(filename, '(a,a,i7.7)') trim(outp%iter_dir),"td.", iter ! name of directory
4185
4186 call output_mxll(outp, namespace, space, gr, st, hm, helmholtz, time, filename)
4187
4188 call profiling_out("TD_WRITE_MAXWELL_DATA")
4189 pop_sub(td_write_maxwell_free_data)
4190 end subroutine td_write_mxll_free_data
4191
4192end module td_write_oct_m
4193
4194!! Local Variables:
4195!! mode: f90
4196!! coding: utf-8
4197!! 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:173
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
Sets the iteration number to the C object.
Definition: write_iter.F90:170
Writes to the corresponding file and adds one to the iteration. Must be called after write_iter_init(...
Definition: write_iter.F90:163
double exp(double __x) __attribute__((__nothrow__
This module contains interfaces for BLAS routines You should not use these routines directly....
Definition: blas.F90:120
This module implements a calculator for the density and defines related functions.
Definition: density.F90:122
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
Definition: density.F90:612
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:193
real(real64), parameter, public m_zero
Definition: global.F90:191
complex(real64), parameter, public m_z0
Definition: global.F90:201
real(real64), parameter, public lmm_r_single_atom
Default local magnetic moments sphere radius for an isolated system.
Definition: global.F90:221
complex(real64), parameter, public m_zi
Definition: global.F90:205
integer, parameter, public max_output_types
Definition: global.F90:151
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public m_one
Definition: global.F90:192
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
Definition: global.F90:214
This module implements the underlying real-space grid.
Definition: grid.F90:119
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:116
subroutine, public io_close(iunit, grp)
Definition: io.F90:466
character(len=max_path_len) function, public io_workpath(path, namespace)
construct path name from given name and namespace
Definition: io.F90:317
subroutine, public io_rm(fname, namespace)
Definition: io.F90:391
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:360
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:401
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:171
integer, parameter, public qkickmode_none
Definition: kick.F90:171
integer, parameter, public qkickmode_sin
Definition: kick.F90:171
subroutine, public kick_function_get(space, mesh, kick, kick_function, iq, to_interpolate)
Definition: kick.F90:994
subroutine, public kick_write(kick, iunit, out)
Definition: kick.F90:892
integer, parameter, public qkickmode_bessel
Definition: kick.F90:171
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:1152
subroutine, public lasers_set_nondipole_parameters(this, ndfield, nd_integration_time)
Set parameters for nondipole SFA calculation.
Definition: lasers.F90:754
logical function, public lasers_with_nondipole_field(lasers)
Check if a nondipole SFA correction should be computed for the given laser.
Definition: lasers.F90:741
integer, parameter, public e_field_electric
Definition: lasers.F90:179
integer, parameter, public e_field_vector_potential
Definition: lasers.F90:179
integer, parameter, public e_field_scalar_potential
Definition: lasers.F90:179
integer pure elemental function, public laser_kind(laser)
Definition: lasers.F90:717
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:1117
integer, parameter, public e_field_magnetic
Definition: lasers.F90:179
integer, parameter, public dft_u_none
Definition: lda_u.F90:203
integer, parameter, public dft_u_acbn0
Definition: lda_u.F90:203
subroutine, public magnetic_local_moments(mesh, st, ions, boundaries, rho, rr, lmm)
Definition: magnetic.F90:261
subroutine, public magnetic_moment(mesh, st, rho, mm)
Definition: magnetic.F90:185
subroutine, public magnetic_total_magnetization(mesh, st, qq, trans_mag)
Definition: magnetic.F90:328
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
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:378
This module defines functions over batches of mesh functions.
Definition: mesh_batch.F90:118
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:120
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1023
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1063
subroutine, public modelmb_sym_all_states(space, mesh, st)
This module contains some common usage patterns of MPI routines.
Definition: mpi_lib.F90:117
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
integer, parameter, public velocity_gauge_dipole
this module contains the low-level part of the output system
Definition: output_low.F90:117
subroutine, public output_modelmb(outp, namespace, space, dir, gr, ions, iter, st)
this module contains the output system
Definition: output.F90:117
subroutine, public output_all(outp, namespace, space, dir, gr, ions, iter, st, hm, ks)
Definition: output.F90:507
subroutine, public output_scalar_pot(outp, namespace, space, dir, mesh, ions, ext_partners, time)
Definition: output.F90:1681
subroutine, public parse_block_string(blk, l, c, res, convert_to_c)
Definition: parser.F90:810
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
integer, parameter, public restart_gs
Definition: restart.F90:156
integer, parameter, public restart_proj
Definition: restart.F90:156
integer, parameter, public restart_type_load
Definition: restart.F90:183
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:167
subroutine, public td_calc_ionch(mesh, st, ch, Nch)
Multiple ionization probabilities calculated form the KS orbital densities C. Ullrich,...
Definition: td_calc.F90:321
subroutine, public td_calc_tvel(gr, st, space, kpoints, vel)
Electronic velocity (to calculate harmonic spectrum...) It is calculated as:
Definition: td_calc.F90:289
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:4092
integer, parameter, public out_total_current
Definition: td_write.F90:203
integer, parameter maxwell_b_field
Definition: td_write.F90:301
integer, parameter, public out_maxwell_max
Definition: td_write.F90:293
subroutine td_write_proj(out_proj, space, mesh, ions, st, gs_st, kick, iter)
Definition: td_write.F90:2385
integer, parameter, public out_q
Definition: td_write.F90:203
integer, parameter, public out_mxll_field
Definition: td_write.F90:203
subroutine calc_projections(mesh, st, gs_st, projections)
This subroutine calculates:
Definition: td_write.F90:2671
integer, parameter out_b_field_surface_y
Definition: td_write.F90:278
subroutine td_write_ionch(out_ionch, mesh, st, iter)
Definition: td_write.F90:2292
integer, parameter, public out_tot_m
Definition: td_write.F90:203
integer, parameter, public out_norm_ks
Definition: td_write.F90:203
integer, parameter out_maxwell_trans_b_field
Definition: td_write.F90:278
integer, parameter, public out_cell_parameters
Definition: td_write.F90:203
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:1574
integer, parameter, public out_proj
Definition: td_write.F90:203
integer, parameter, public out_partial_charges
Definition: td_write.F90:203
integer, parameter, public out_separate_coords
Definition: td_write.F90:203
subroutine td_write_total_current(out_total_current, space, mesh, st, iter)
Definition: td_write.F90:3075
subroutine, public td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, dt)
Definition: td_write.F90:1239
subroutine td_write_effective_u(out_coords, lda_u, iter)
Definition: td_write.F90:3395
integer, parameter maxwell_trans_field
Definition: td_write.F90:296
subroutine td_write_acc(out_acc, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
Definition: td_write.F90:1925
subroutine, public td_write_mxll_init(writ, namespace, iter, dt)
Definition: td_write.F90:3617
subroutine, public td_write_mxll_end(writ)
Definition: td_write.F90:3805
subroutine td_write_mxll_field(out_mxll, space, hm, dt, iter)
Definition: td_write.F90:3310
integer, parameter out_b_field_surface_x
Definition: td_write.F90:278
integer, parameter out_maxwell_long_e_field
Definition: td_write.F90:278
integer, parameter, public out_kp_proj
Definition: td_write.F90:203
integer, parameter, public out_magnets
Definition: td_write.F90:203
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:1527
subroutine td_write_electric_field_box_surface(out_field_surf, st, dim, iter)
Definition: td_write.F90:4012
subroutine td_write_floquet(namespace, space, hm, ext_partners, gr, st, iter)
Definition: td_write.F90:2799
integer, parameter out_e_field_surface_y
Definition: td_write.F90:278
integer, parameter, public out_angular
Definition: td_write.F90:203
subroutine td_write_populations(out_populations, namespace, space, mesh, st, writ, dt, iter)
Definition: td_write.F90:1845
integer, parameter, public out_max
Definition: td_write.F90:203
integer, parameter out_maxwell_long_b_field
Definition: td_write.F90:278
integer, parameter, public out_energy
Definition: td_write.F90:203
integer, parameter, public out_spin
Definition: td_write.F90:203
subroutine td_write_ftchd(out_ftchd, space, mesh, st, kick, iter)
Definition: td_write.F90:1701
subroutine td_write_partial_charges(out_partial_charges, mesh, st, ions, iter)
Definition: td_write.F90:3218
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:273
subroutine td_write_laser(out_laser, space, lasers, dt, iter)
Definition: td_write.F90:2026
integer, parameter out_maxwell_total_b_field
Definition: td_write.F90:278
integer, parameter, public out_ftchd
Definition: td_write.F90:203
integer, parameter, public out_separate_velocity
Definition: td_write.F90:203
subroutine td_write_tot_mag(out_magnets, mesh, st, kick, iter)
Definition: td_write.F90:1389
subroutine td_write_q(out_q, space, ks, iter)
Definition: td_write.F90:3267
integer, parameter, public out_floquet
Definition: td_write.F90:203
subroutine, public td_write_mxll_free_data(namespace, space, gr, st, hm, helmholtz, outp, iter, time)
Definition: td_write.F90:4263
integer, parameter, public out_acc
Definition: td_write.F90:203
integer, parameter, public out_ion_ch
Definition: td_write.F90:203
integer, parameter maxwell_long_field
Definition: td_write.F90:296
integer, parameter, public out_n_ex
Definition: td_write.F90:203
integer, parameter out_b_field_surface_z
Definition: td_write.F90:278
subroutine td_write_proj_kp(mesh, kpoints, st, gs_st, namespace, iter)
Definition: td_write.F90:2705
integer, parameter, public out_temperature
Definition: td_write.F90:203
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:3484
subroutine, public td_write_data(writ)
Definition: td_write.F90:1205
subroutine td_write_total_heat_current(write_obj, space, hm, gr, st, iter)
Definition: td_write.F90:3162
integer, parameter out_e_field_surface_z
Definition: td_write.F90:278
integer, parameter maxwell_total_field
Definition: td_write.F90:296
integer, parameter, public out_coords
Definition: td_write.F90:203
integer, parameter out_maxwell_total_e_field
Definition: td_write.F90:278
integer, parameter, public out_laser
Definition: td_write.F90:203
integer, parameter, public out_eigs
Definition: td_write.F90:203
subroutine td_write_vel(out_vel, gr, st, space, kpoints, iter)
Definition: td_write.F90:1977
integer, parameter, public out_total_heat_current
Definition: td_write.F90:203
integer, parameter out_e_field_surface_x
Definition: td_write.F90:278
subroutine, public td_write_kick(outp, namespace, space, mesh, kick, ions, iter)
Definition: td_write.F90:337
subroutine, public td_write_end(writ)
Definition: td_write.F90:989
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:1456
subroutine td_write_eigs(out_eigs, st, iter)
Definition: td_write.F90:2231
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:1036
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:2561
subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
Definition: td_write.F90:4171
subroutine td_write_spin(out_spin, mesh, st, iter)
Definition: td_write.F90:1280
integer, parameter, public out_vel
Definition: td_write.F90:203
integer, parameter, public out_gauge_field
Definition: td_write.F90:203
subroutine td_write_temperature(out_temperature, ions, iter)
Definition: td_write.F90:1808
integer, parameter maxwell_e_field
Definition: td_write.F90:301
integer, parameter, public out_populations
Definition: td_write.F90:203
subroutine, public td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
Definition: td_write.F90:3821
subroutine td_write_cell_parameters(file_handle, ions, iter)
Write the cell parameters as a function of time.
Definition: td_write.F90:3553
subroutine td_write_local_magnetic_moments(out_magnets, gr, st, ions, lmm_r, iter)
Definition: td_write.F90:1332
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:375
integer, parameter out_maxwell_energy
Definition: td_write.F90:278
subroutine td_write_energy(out_energy, hm, iter, ke)
Definition: td_write.F90:2148
subroutine td_write_maxwell_energy(out_maxwell_energy, hm, iter)
Definition: td_write.F90:3951
integer, parameter, public out_separate_forces
Definition: td_write.F90:203
integer, parameter out_maxwell_trans_e_field
Definition: td_write.F90:278
type(type_t), public type_cmplx
Definition: types.F90:136
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:225
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:120
character pure function, public index2axis(idir)
Definition: utils.F90:205
subroutine, public v_ks_calculate_current(this, calc_cur)
Definition: v_ks.F90:1494
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:749
Explicit interfaces to C functions, defined in write_iter_low.cc.
Definition: write_iter.F90:116
subroutine, public write_iter_header(out, string)
Definition: write_iter.F90:249
subroutine, public write_iter_string(out, string)
Definition: write_iter.F90:265
subroutine, public write_iter_init(out, iter, factor, file)
Definition: write_iter.F90:229
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:171
Describes mesh distribution to nodes.
Definition: mesh.F90:187
This is defined even when running serial.
Definition: mpi.F90:144
Stores all communicators and groups.
Definition: multicomm.F90:208
output handler class
Definition: output_low.F90:166
The states_elec_t class contains all electronic wave functions.
Time-dependent Write Properties.
Definition: td_write.F90:319
int true(void)
subroutine dipole_matrix_elements(dir)
Definition: td_write.F90:2505