LCOV - code coverage report
Current view: top level - src/emd - rt_propagation_output.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 331 338 97.9 %
Date: 2023-03-30 11:55:16 Functions: 10 10 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routine for the real time propagation output.
      10             : !> \author Florian Schiffmann (02.09)
      11             : ! **************************************************************************************************
      12             : 
      13             : MODULE rt_propagation_output
      14             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      15             :    USE cp_control_types,                ONLY: dft_control_type
      16             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      17             :                                               dbcsr_allocate_matrix_set,&
      18             :                                               dbcsr_deallocate_matrix_set
      19             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add
      20             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      21             :                                               cp_fm_struct_double,&
      22             :                                               cp_fm_struct_release,&
      23             :                                               cp_fm_struct_type
      24             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      25             :                                               cp_fm_get_info,&
      26             :                                               cp_fm_release,&
      27             :                                               cp_fm_set_all,&
      28             :                                               cp_fm_type
      29             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      30             :                                               cp_logger_get_default_io_unit,&
      31             :                                               cp_logger_get_default_unit_nr,&
      32             :                                               cp_logger_type
      33             :    USE cp_output_handling,              ONLY: cp_iter_string,&
      34             :                                               cp_p_file,&
      35             :                                               cp_print_key_finished_output,&
      36             :                                               cp_print_key_should_output,&
      37             :                                               cp_print_key_unit_nr
      38             :    USE cp_realspace_grid_cube,          ONLY: cp_pw_to_cube
      39             :    USE dbcsr_api,                       ONLY: &
      40             :         dbcsr_add, dbcsr_binary_write, dbcsr_checksum, dbcsr_copy, dbcsr_create, &
      41             :         dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_filter, dbcsr_get_info, &
      42             :         dbcsr_get_occupation, dbcsr_init_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      43             :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_scale, &
      44             :         dbcsr_set, dbcsr_type
      45             :    USE efield_utils,                    ONLY: make_field
      46             :    USE input_constants,                 ONLY: ehrenfest,&
      47             :                                               real_time_propagation
      48             :    USE input_section_types,             ONLY: section_get_ivals,&
      49             :                                               section_vals_get_subs_vals,&
      50             :                                               section_vals_type
      51             :    USE kahan_sum,                       ONLY: accurate_sum
      52             :    USE kinds,                           ONLY: default_path_length,&
      53             :                                               dp
      54             :    USE machine,                         ONLY: m_flush
      55             :    USE message_passing,                 ONLY: mp_comm_type
      56             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      57             :    USE particle_list_types,             ONLY: particle_list_type
      58             :    USE particle_types,                  ONLY: particle_type
      59             :    USE pw_env_types,                    ONLY: pw_env_get,&
      60             :                                               pw_env_type
      61             :    USE pw_methods,                      ONLY: pw_zero
      62             :    USE pw_pool_types,                   ONLY: pw_pool_create_pw,&
      63             :                                               pw_pool_give_back_pw,&
      64             :                                               pw_pool_type
      65             :    USE pw_types,                        ONLY: COMPLEXDATA1D,&
      66             :                                               REALDATA3D,&
      67             :                                               REALSPACE,&
      68             :                                               RECIPROCALSPACE,&
      69             :                                               pw_type
      70             :    USE qs_environment_types,            ONLY: get_qs_env,&
      71             :                                               qs_environment_type
      72             :    USE qs_kind_types,                   ONLY: get_qs_kind_set,&
      73             :                                               qs_kind_type
      74             :    USE qs_linres_current,               ONLY: calculate_jrho_resp
      75             :    USE qs_linres_types,                 ONLY: current_env_type
      76             :    USE qs_mo_io,                        ONLY: write_rt_mos_to_restart
      77             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      78             :                                               qs_rho_type
      79             :    USE qs_scf_post_gpw,                 ONLY: qs_scf_post_moments,&
      80             :                                               write_mo_dependent_results,&
      81             :                                               write_mo_free_results
      82             :    USE qs_scf_post_tb,                  ONLY: scf_post_calculation_tb
      83             :    USE qs_scf_types,                    ONLY: qs_scf_env_type
      84             :    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
      85             :                                               qs_subsys_type
      86             :    USE rt_propagation_types,            ONLY: get_rtp,&
      87             :                                               rt_prop_type
      88             :    USE rt_propagation_utils,            ONLY: calculate_P_imaginary,&
      89             :                                               write_rtp_mo_cubes,&
      90             :                                               write_rtp_mos_to_output_unit
      91             : #include "../base/base_uses.f90"
      92             : 
      93             :    IMPLICIT NONE
      94             : 
      95             :    PRIVATE
      96             : 
      97             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_output'
      98             : 
      99             :    PUBLIC :: rt_prop_output, &
     100             :              rt_convergence, &
     101             :              rt_convergence_density, &
     102             :              report_density_occupation
     103             : 
     104             : CONTAINS
     105             : 
     106             : ! **************************************************************************************************
     107             : !> \brief ...
     108             : !> \param qs_env ...
     109             : !> \param run_type ...
     110             : !> \param delta_iter ...
     111             : !> \param used_time ...
     112             : ! **************************************************************************************************
     113        4412 :    SUBROUTINE rt_prop_output(qs_env, run_type, delta_iter, used_time)
     114             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     115             :       INTEGER, INTENT(in)                                :: run_type
     116             :       REAL(dp), INTENT(in), OPTIONAL                     :: delta_iter, used_time
     117             : 
     118             :       INTEGER                                            :: n_electrons, nspin, output_unit, spin
     119             :       REAL(dp)                                           :: orthonormality, tot_rho_r
     120        2206 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: qs_tot_rho_r
     121        2206 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     122        2206 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     123             :       TYPE(cp_logger_type), POINTER                      :: logger
     124        2206 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, P_im, rho_new
     125             :       TYPE(dft_control_type), POINTER                    :: dft_control
     126        2206 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     127        2206 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     128             :       TYPE(qs_rho_type), POINTER                         :: rho
     129             :       TYPE(rt_prop_type), POINTER                        :: rtp
     130             :       TYPE(section_vals_type), POINTER                   :: dft_section, input, rtp_section
     131             : 
     132        2206 :       NULLIFY (logger, dft_control)
     133             : 
     134        4412 :       logger => cp_get_default_logger()
     135             :       CALL get_qs_env(qs_env, &
     136             :                       rtp=rtp, &
     137             :                       matrix_s=matrix_s, &
     138             :                       input=input, &
     139             :                       rho=rho, &
     140             :                       particle_set=particle_set, &
     141             :                       atomic_kind_set=atomic_kind_set, &
     142             :                       qs_kind_set=qs_kind_set, &
     143        2206 :                       dft_control=dft_control)
     144             : 
     145        2206 :       rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
     146             : 
     147        2206 :       CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
     148        2206 :       n_electrons = n_electrons - dft_control%charge
     149             : 
     150        2206 :       CALL qs_rho_get(rho_struct=rho, tot_rho_r=qs_tot_rho_r)
     151             : 
     152        2206 :       tot_rho_r = accurate_sum(qs_tot_rho_r)
     153             : 
     154             :       output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
     155        2206 :                                          extension=".scfLog")
     156             : 
     157        2206 :       IF (output_unit > 0) THEN
     158             :          WRITE (output_unit, FMT="(/,(T3,A,T40,I5))") &
     159        1103 :             "Information at iteration step:", rtp%iter
     160             :          WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
     161        1103 :             "Total electronic density (r-space): ", &
     162        1103 :             tot_rho_r, &
     163             :             tot_rho_r + &
     164        2206 :             REAL(n_electrons, dp)
     165             :          WRITE (UNIT=output_unit, FMT="((T3,A,T59,F22.14))") &
     166        1103 :             "Total energy:", rtp%energy_new
     167        1103 :          IF (run_type == ehrenfest) &
     168             :             WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
     169         624 :             "Energy difference to previous iteration step:", rtp%energy_new - rtp%energy_old
     170        1103 :          IF (run_type == real_time_propagation) &
     171             :             WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
     172         479 :             "Energy difference to initial state:", rtp%energy_new - rtp%energy_old
     173        1103 :          IF (PRESENT(delta_iter)) &
     174             :             WRITE (UNIT=output_unit, FMT="((T3,A,T61,E20.6))") &
     175        1103 :             "Convergence:", delta_iter
     176        1103 :          IF (rtp%converged) THEN
     177         277 :             IF (run_type == real_time_propagation) &
     178             :                WRITE (UNIT=output_unit, FMT="((T3,A,T61,F12.2))") &
     179         151 :                "Time needed for propagation:", used_time
     180             :             WRITE (UNIT=output_unit, FMT="(/,(T3,A,3X,F16.14))") &
     181         277 :                "CONVERGENCE REACHED", rtp%energy_new - rtp%energy_old
     182             :          END IF
     183             :       END IF
     184             : 
     185        2206 :       IF (rtp%converged) THEN
     186         554 :          IF (.NOT. rtp%linear_scaling) THEN
     187         372 :             CALL get_rtp(rtp=rtp, mos_new=mos_new)
     188             :             CALL rt_calculate_orthonormality(orthonormality, &
     189         372 :                                              mos_new, matrix_s(1)%matrix)
     190         372 :             IF (output_unit > 0) &
     191             :                WRITE (output_unit, FMT="(/,(T3,A,T60,F20.10))") &
     192         186 :                "Max deviation from orthonormalization:", orthonormality
     193             :          END IF
     194             :       END IF
     195             : 
     196        2206 :       IF (output_unit > 0) &
     197        1103 :          CALL m_flush(output_unit)
     198             :       CALL cp_print_key_finished_output(output_unit, logger, rtp_section, &
     199        2206 :                                         "PRINT%PROGRAM_RUN_INFO")
     200             : 
     201        2206 :       IF (rtp%converged) THEN
     202         554 :          dft_section => section_vals_get_subs_vals(input, "DFT")
     203         554 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     204             :                                               dft_section, "REAL_TIME_PROPAGATION%PRINT%FIELD"), cp_p_file)) &
     205           4 :             CALL print_field_applied(qs_env, dft_section)
     206         554 :          CALL make_moment(qs_env)
     207         554 :          IF (.NOT. dft_control%qs_control%dftb) THEN
     208         434 :             CALL write_available_results(qs_env=qs_env, rtp=rtp)
     209             :          END IF
     210         554 :          IF (rtp%linear_scaling) THEN
     211         182 :             CALL get_rtp(rtp=rtp, rho_new=rho_new)
     212         182 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     213             :                                                  dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART"), cp_p_file)) THEN
     214          96 :                CALL write_rt_p_to_restart(rho_new, .FALSE.)
     215             :             END IF
     216         182 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     217             :                                                  dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY"), cp_p_file)) THEN
     218           2 :                CALL write_rt_p_to_restart(rho_new, .TRUE.)
     219             :             END IF
     220         182 :             IF (.NOT. dft_control%qs_control%dftb) THEN
     221             :                !Not sure if these things could also work with dftb or not
     222         182 :                IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     223             :                                                     dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
     224          64 :                   DO spin = 1, SIZE(rho_new)/2
     225          64 :                      CALL rt_current(qs_env, rho_new(2*spin)%matrix, dft_section, spin, SIZE(rho_new)/2)
     226             :                   END DO
     227             :                END IF
     228             :             END IF
     229             :          ELSE
     230         372 :             CALL get_rtp(rtp=rtp, mos_new=mos_new)
     231         372 :             IF (.NOT. dft_control%qs_control%dftb) THEN
     232         252 :                IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     233             :                                                     dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
     234           8 :                   NULLIFY (P_im)
     235           8 :                   nspin = SIZE(mos_new)/2
     236           8 :                   CALL dbcsr_allocate_matrix_set(P_im, nspin)
     237          16 :                   DO spin = 1, nspin
     238           8 :                      CALL dbcsr_init_p(P_im(spin)%matrix)
     239          16 :                      CALL dbcsr_create(P_im(spin)%matrix, template=matrix_s(1)%matrix, matrix_type="N")
     240             :                   END DO
     241           8 :                   CALL calculate_P_imaginary(qs_env, rtp, P_im)
     242          16 :                   DO spin = 1, nspin
     243          16 :                      CALL rt_current(qs_env, P_im(spin)%matrix, dft_section, spin, nspin)
     244             :                   END DO
     245           8 :                   CALL dbcsr_deallocate_matrix_set(P_im)
     246             :                END IF
     247             :             END IF
     248             :             CALL write_rt_mos_to_restart(qs_env%mos, mos_new, particle_set, &
     249         372 :                                          dft_section, qs_kind_set)
     250             :          END IF
     251             :       END IF
     252             : 
     253        2206 :       rtp%energy_old = rtp%energy_new
     254             : 
     255        2206 :       IF (.NOT. rtp%converged .AND. rtp%iter >= dft_control%rtp_control%max_iter) &
     256             :          CALL cp_abort(__LOCATION__, "EMD did not converge, either increase MAX_ITER "// &
     257           0 :                        "or use a smaller TIMESTEP")
     258             : 
     259        2206 :    END SUBROUTINE rt_prop_output
     260             : 
     261             : ! **************************************************************************************************
     262             : !> \brief computes the effective orthonormality of a set of mos given an s-matrix
     263             : !>        orthonormality is the max deviation from unity of the C^T S C
     264             : !> \param orthonormality ...
     265             : !> \param mos_new ...
     266             : !> \param matrix_s ...
     267             : !> \author Florian Schiffmann (02.09)
     268             : ! **************************************************************************************************
     269         372 :    SUBROUTINE rt_calculate_orthonormality(orthonormality, mos_new, matrix_s)
     270             :       REAL(KIND=dp), INTENT(out)                         :: orthonormality
     271             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     272             :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
     273             : 
     274             :       CHARACTER(len=*), PARAMETER :: routineN = 'rt_calculate_orthonormality'
     275             : 
     276             :       INTEGER                                            :: handle, i, im, ispin, j, k, n, &
     277             :                                                             ncol_local, nrow_local, nspin, re
     278         372 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     279             :       REAL(KIND=dp)                                      :: alpha, max_alpha, max_beta
     280             :       TYPE(cp_fm_struct_type), POINTER                   :: tmp_fm_struct
     281             :       TYPE(cp_fm_type)                                   :: overlap_re, svec_im, svec_re
     282             : 
     283         372 :       NULLIFY (tmp_fm_struct)
     284             : 
     285         372 :       CALL timeset(routineN, handle)
     286             : 
     287         372 :       nspin = SIZE(mos_new)/2
     288         372 :       max_alpha = 0.0_dp
     289         372 :       max_beta = 0.0_dp
     290         826 :       DO ispin = 1, nspin
     291         454 :          re = ispin*2 - 1
     292         454 :          im = ispin*2
     293             :          ! get S*C
     294         454 :          CALL cp_fm_create(svec_re, mos_new(im)%matrix_struct)
     295         454 :          CALL cp_fm_create(svec_im, mos_new(im)%matrix_struct)
     296             :          CALL cp_fm_get_info(mos_new(im), &
     297         454 :                              nrow_global=n, ncol_global=k)
     298             :          CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(re), &
     299         454 :                                       svec_re, k)
     300             :          CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(im), &
     301         454 :                                       svec_im, k)
     302             : 
     303             :          ! get C^T (S*C)
     304             :          CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
     305             :                                   para_env=mos_new(re)%matrix_struct%para_env, &
     306         454 :                                   context=mos_new(re)%matrix_struct%context)
     307         454 :          CALL cp_fm_create(overlap_re, tmp_fm_struct)
     308             : 
     309         454 :          CALL cp_fm_struct_release(tmp_fm_struct)
     310             : 
     311             :          CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(re), &
     312         454 :                             svec_re, 0.0_dp, overlap_re)
     313             :          CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(im), &
     314         454 :                             svec_im, 1.0_dp, overlap_re)
     315             : 
     316         454 :          CALL cp_fm_release(svec_re)
     317         454 :          CALL cp_fm_release(svec_im)
     318             : 
     319             :          CALL cp_fm_get_info(overlap_re, nrow_local=nrow_local, ncol_local=ncol_local, &
     320         454 :                              row_indices=row_indices, col_indices=col_indices)
     321        1196 :          DO i = 1, nrow_local
     322        4198 :             DO j = 1, ncol_local
     323        3002 :                alpha = overlap_re%local_data(i, j)
     324        3002 :                IF (row_indices(i) .EQ. col_indices(j)) alpha = alpha - 1.0_dp
     325        3744 :                max_alpha = MAX(max_alpha, ABS(alpha))
     326             :             END DO
     327             :          END DO
     328        1734 :          CALL cp_fm_release(overlap_re)
     329             :       END DO
     330         372 :       CALL mos_new(1)%matrix_struct%para_env%max(max_alpha)
     331         372 :       CALL mos_new(1)%matrix_struct%para_env%max(max_beta)
     332         372 :       orthonormality = max_alpha
     333             : 
     334         372 :       CALL timestop(handle)
     335             : 
     336         372 :    END SUBROUTINE rt_calculate_orthonormality
     337             : 
     338             : ! **************************************************************************************************
     339             : !> \brief computes the convergence criterion for RTP and EMD
     340             : !> \param rtp ...
     341             : !> \param matrix_s Overlap matrix without the derivatives
     342             : !> \param delta_mos ...
     343             : !> \param delta_eps ...
     344             : !> \author Florian Schiffmann (02.09)
     345             : ! **************************************************************************************************
     346             : 
     347        1390 :    SUBROUTINE rt_convergence(rtp, matrix_s, delta_mos, delta_eps)
     348             :       TYPE(rt_prop_type), POINTER                        :: rtp
     349             :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     350             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: delta_mos
     351             :       REAL(dp), INTENT(out)                              :: delta_eps
     352             : 
     353             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'rt_convergence'
     354             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     355             : 
     356             :       INTEGER                                            :: handle, i, icol, im, ispin, j, lcol, &
     357             :                                                             lrow, nao, newdim, nmo, nspin, re
     358             :       LOGICAL                                            :: double_col, double_row
     359             :       REAL(KIND=dp)                                      :: alpha, max_alpha
     360             :       TYPE(cp_fm_struct_type), POINTER                   :: newstruct, newstruct1, tmp_fm_struct
     361             :       TYPE(cp_fm_type)                                   :: work, work1, work2
     362        1390 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     363             : 
     364        1390 :       NULLIFY (tmp_fm_struct)
     365             : 
     366        1390 :       CALL timeset(routineN, handle)
     367             : 
     368        1390 :       CALL get_rtp(rtp=rtp, mos_new=mos_new)
     369             : 
     370        1390 :       nspin = SIZE(delta_mos)/2
     371        1390 :       max_alpha = 0.0_dp
     372             : 
     373        4814 :       DO i = 1, SIZE(mos_new)
     374        4814 :          CALL cp_fm_scale_and_add(-one, delta_mos(i), one, mos_new(i))
     375             :       END DO
     376             : 
     377        3102 :       DO ispin = 1, nspin
     378        1712 :          re = ispin*2 - 1
     379        1712 :          im = ispin*2
     380             : 
     381        1712 :          double_col = .TRUE.
     382        1712 :          double_row = .FALSE.
     383             :          CALL cp_fm_struct_double(newstruct, &
     384             :                                   delta_mos(re)%matrix_struct, &
     385             :                                   delta_mos(re)%matrix_struct%context, &
     386             :                                   double_col, &
     387        1712 :                                   double_row)
     388             : 
     389        1712 :          CALL cp_fm_create(work, matrix_struct=newstruct)
     390        1712 :          CALL cp_fm_create(work1, matrix_struct=newstruct)
     391             : 
     392             :          CALL cp_fm_get_info(delta_mos(re), ncol_local=lcol, ncol_global=nmo, &
     393        1712 :                              nrow_global=nao)
     394        1712 :          CALL cp_fm_get_info(work, ncol_global=newdim)
     395             : 
     396        1712 :          CALL cp_fm_set_all(work, zero, zero)
     397             : 
     398        7080 :          DO icol = 1, lcol
     399       40040 :             work%local_data(:, icol) = delta_mos(re)%local_data(:, icol)
     400       41752 :             work%local_data(:, icol + lcol) = delta_mos(im)%local_data(:, icol)
     401             :          END DO
     402             : 
     403        1712 :          CALL cp_dbcsr_sm_fm_multiply(matrix_s, work, work1, ncol=newdim)
     404             : 
     405        1712 :          CALL cp_fm_release(work)
     406             : 
     407             :          CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nmo, &
     408             :                                   para_env=delta_mos(re)%matrix_struct%para_env, &
     409        1712 :                                   context=delta_mos(re)%matrix_struct%context)
     410             :          CALL cp_fm_struct_double(newstruct1, &
     411             :                                   tmp_fm_struct, &
     412             :                                   delta_mos(re)%matrix_struct%context, &
     413             :                                   double_col, &
     414        1712 :                                   double_row)
     415             : 
     416        1712 :          CALL cp_fm_create(work, matrix_struct=newstruct1)
     417        1712 :          CALL cp_fm_create(work2, matrix_struct=newstruct1)
     418             : 
     419             :          CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(re), &
     420        1712 :                             work1, zero, work)
     421             : 
     422             :          CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(im), &
     423        1712 :                             work1, zero, work2)
     424             : 
     425        1712 :          CALL cp_fm_get_info(work, nrow_local=lrow)
     426        4396 :          DO i = 1, lrow
     427       15172 :             DO j = 1, lcol
     428             :                alpha = SQRT((work%local_data(i, j) + work2%local_data(i, j + lcol))**2 + &
     429       10776 :                             (work%local_data(i, j + lcol) - work2%local_data(i, j))**2)
     430       13460 :                max_alpha = MAX(max_alpha, ABS(alpha))
     431             :             END DO
     432             :          END DO
     433             : 
     434        1712 :          CALL cp_fm_release(work)
     435        1712 :          CALL cp_fm_release(work1)
     436        1712 :          CALL cp_fm_release(work2)
     437        1712 :          CALL cp_fm_struct_release(tmp_fm_struct)
     438        1712 :          CALL cp_fm_struct_release(newstruct)
     439        6526 :          CALL cp_fm_struct_release(newstruct1)
     440             : 
     441             :       END DO
     442             : 
     443        1390 :       CALL delta_mos(1)%matrix_struct%para_env%max(max_alpha)
     444        1390 :       delta_eps = SQRT(max_alpha)
     445             : 
     446        1390 :       CALL timestop(handle)
     447             : 
     448        1390 :    END SUBROUTINE rt_convergence
     449             : 
     450             : ! **************************************************************************************************
     451             : !> \brief computes the convergence criterion for RTP and EMD based on the density matrix
     452             : !> \param rtp ...
     453             : !> \param delta_P ...
     454             : !> \param delta_eps ...
     455             : !> \author Samuel Andermatt (02.14)
     456             : ! **************************************************************************************************
     457             : 
     458        1632 :    SUBROUTINE rt_convergence_density(rtp, delta_P, delta_eps)
     459             : 
     460             :       TYPE(rt_prop_type), POINTER                        :: rtp
     461             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: delta_P
     462             :       REAL(dp), INTENT(out)                              :: delta_eps
     463             : 
     464             :       CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence_density'
     465             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     466             : 
     467             :       INTEGER                                            :: col_atom, group_handle, handle, i, &
     468             :                                                             ispin, row_atom
     469             :       REAL(dp)                                           :: alpha, max_alpha
     470         816 :       REAL(dp), DIMENSION(:), POINTER                    :: block_values
     471             :       TYPE(dbcsr_iterator_type)                          :: iter
     472         816 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_new
     473             :       TYPE(dbcsr_type), POINTER                          :: tmp
     474             :       TYPE(mp_comm_type)                                 :: group
     475             : 
     476         816 :       CALL timeset(routineN, handle)
     477             : 
     478         816 :       CALL get_rtp(rtp=rtp, rho_new=rho_new)
     479             : 
     480        2940 :       DO i = 1, SIZE(rho_new)
     481        2940 :          CALL dbcsr_add(delta_P(i)%matrix, rho_new(i)%matrix, one, -one)
     482             :       END DO
     483             :       !get the maximum value of delta_P
     484        2940 :       DO i = 1, SIZE(delta_P)
     485             :          !square all entries of both matrices
     486        2124 :          CALL dbcsr_iterator_start(iter, delta_P(i)%matrix)
     487       12356 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
     488       10232 :             CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
     489      712140 :             block_values = block_values*block_values
     490             :          END DO
     491        2940 :          CALL dbcsr_iterator_stop(iter)
     492             :       END DO
     493             :       NULLIFY (tmp)
     494         816 :       ALLOCATE (tmp)
     495         816 :       CALL dbcsr_create(tmp, template=delta_P(1)%matrix, matrix_type="N")
     496        1878 :       DO ispin = 1, SIZE(delta_P)/2
     497        1062 :          CALL dbcsr_desymmetrize(delta_P(2*ispin - 1)%matrix, tmp)
     498        1878 :          CALL dbcsr_add(delta_P(2*ispin)%matrix, tmp, one, one)
     499             :       END DO
     500             :       !the absolute values are now in the even entries of delta_P
     501         816 :       max_alpha = zero
     502        1878 :       DO ispin = 1, SIZE(delta_P)/2
     503        1062 :          CALL dbcsr_iterator_start(iter, delta_P(2*ispin)%matrix)
     504        6275 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
     505        5213 :             CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
     506      365095 :             alpha = MAXVAL(block_values)
     507        5213 :             IF (alpha > max_alpha) max_alpha = alpha
     508             :          END DO
     509        1878 :          CALL dbcsr_iterator_stop(iter)
     510             :       END DO
     511         816 :       CALL dbcsr_get_info(delta_P(1)%matrix, group=group_handle)
     512         816 :       CALL group%set_handle(group_handle)
     513         816 :       CALL group%max(max_alpha)
     514         816 :       delta_eps = SQRT(max_alpha)
     515         816 :       CALL dbcsr_deallocate_matrix(tmp)
     516         816 :       CALL timestop(handle)
     517             : 
     518         816 :    END SUBROUTINE rt_convergence_density
     519             : 
     520             : ! **************************************************************************************************
     521             : !> \brief interface to qs_moments. Does only work for nonperiodic dipole
     522             : !> \param qs_env ...
     523             : !> \author Florian Schiffmann (02.09)
     524             : ! **************************************************************************************************
     525             : 
     526         554 :    SUBROUTINE make_moment(qs_env)
     527             : 
     528             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     529             : 
     530             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'make_moment'
     531             : 
     532             :       INTEGER                                            :: handle, output_unit
     533             :       TYPE(cp_logger_type), POINTER                      :: logger
     534             :       TYPE(dft_control_type), POINTER                    :: dft_control
     535             : 
     536         554 :       CALL timeset(routineN, handle)
     537             : 
     538         554 :       NULLIFY (dft_control)
     539             : 
     540         554 :       logger => cp_get_default_logger()
     541         554 :       output_unit = cp_logger_get_default_io_unit(logger)
     542         554 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     543         554 :       IF (dft_control%qs_control%dftb) THEN
     544         120 :          CALL scf_post_calculation_tb(qs_env, "DFTB", .FALSE.)
     545         434 :       ELSE IF (dft_control%qs_control%xtb) THEN
     546          60 :          CALL scf_post_calculation_tb(qs_env, "xTB", .FALSE.)
     547             :       ELSE
     548         374 :          CALL qs_scf_post_moments(qs_env%input, logger, qs_env, output_unit)
     549             :       END IF
     550         554 :       CALL timestop(handle)
     551             : 
     552         554 :    END SUBROUTINE make_moment
     553             : 
     554             : ! **************************************************************************************************
     555             : !> \brief Reports the sparsity pattern of the complex density matrix
     556             : !> \param filter_eps ...
     557             : !> \param rho ...
     558             : !> \author Samuel Andermatt (09.14)
     559             : ! **************************************************************************************************
     560             : 
     561         182 :    SUBROUTINE report_density_occupation(filter_eps, rho)
     562             : 
     563             :       REAL(KIND=dp)                                      :: filter_eps
     564             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho
     565             : 
     566             :       CHARACTER(len=*), PARAMETER :: routineN = 'report_density_occupation'
     567             : 
     568             :       INTEGER                                            :: handle, i, im, ispin, re, unit_nr
     569             :       REAL(KIND=dp)                                      :: eps, occ
     570             :       TYPE(cp_logger_type), POINTER                      :: logger
     571         182 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: tmp
     572             : 
     573         182 :       CALL timeset(routineN, handle)
     574             : 
     575         182 :       logger => cp_get_default_logger()
     576         182 :       unit_nr = cp_logger_get_default_io_unit(logger)
     577         182 :       NULLIFY (tmp)
     578         182 :       CALL dbcsr_allocate_matrix_set(tmp, SIZE(rho))
     579         686 :       DO i = 1, SIZE(rho)
     580         504 :          CALL dbcsr_init_p(tmp(i)%matrix)
     581         504 :          CALL dbcsr_create(tmp(i)%matrix, template=rho(i)%matrix)
     582         686 :          CALL dbcsr_copy(tmp(i)%matrix, rho(i)%matrix)
     583             :       END DO
     584         434 :       DO ispin = 1, SIZE(rho)/2
     585         252 :          re = 2*ispin - 1
     586         252 :          im = 2*ispin
     587         252 :          eps = MAX(filter_eps, 1.0E-11_dp)
     588        2338 :          DO WHILE (eps < 1.1_dp)
     589        2086 :             CALL dbcsr_filter(tmp(re)%matrix, eps)
     590        2086 :             occ = dbcsr_get_occupation(tmp(re)%matrix)
     591        3129 :             IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
     592        2086 :                ispin, " eps ", eps, " real: ", occ
     593        2086 :             eps = eps*10
     594             :          END DO
     595         252 :          eps = MAX(filter_eps, 1.0E-11_dp)
     596        2520 :          DO WHILE (eps < 1.1_dp)
     597        2086 :             CALL dbcsr_filter(tmp(im)%matrix, eps)
     598        2086 :             occ = dbcsr_get_occupation(tmp(im)%matrix)
     599        3129 :             IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
     600        2086 :                ispin, " eps ", eps, " imag: ", occ
     601        2086 :             eps = eps*10.0_dp
     602             :          END DO
     603             :       END DO
     604         182 :       CALL dbcsr_deallocate_matrix_set(tmp)
     605         182 :       CALL timestop(handle)
     606             : 
     607         182 :    END SUBROUTINE report_density_occupation
     608             : 
     609             : ! **************************************************************************************************
     610             : !> \brief Writes the density matrix and the atomic positions to a restart file
     611             : !> \param rho_new ...
     612             : !> \param history ...
     613             : !> \author Samuel Andermatt (09.14)
     614             : ! **************************************************************************************************
     615             : 
     616          98 :    SUBROUTINE write_rt_p_to_restart(rho_new, history)
     617             : 
     618             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_new
     619             :       LOGICAL                                            :: history
     620             : 
     621             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_p_to_restart'
     622             : 
     623             :       CHARACTER(LEN=default_path_length)                 :: file_name, project_name
     624             :       INTEGER                                            :: handle, im, ispin, re, unit_nr
     625             :       REAL(KIND=dp)                                      :: cs_pos
     626             :       TYPE(cp_logger_type), POINTER                      :: logger
     627             : 
     628          98 :       CALL timeset(routineN, handle)
     629          98 :       logger => cp_get_default_logger()
     630          98 :       IF (logger%para_env%is_source()) THEN
     631          49 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     632             :       ELSE
     633             :          unit_nr = -1
     634             :       END IF
     635             : 
     636          98 :       project_name = logger%iter_info%project_name
     637         234 :       DO ispin = 1, SIZE(rho_new)/2
     638         136 :          re = 2*ispin - 1
     639         136 :          im = 2*ispin
     640         136 :          IF (history) THEN
     641             :             WRITE (file_name, '(A,I0,A)') &
     642           2 :                TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
     643             :          ELSE
     644         134 :             WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
     645             :          END IF
     646         136 :          cs_pos = dbcsr_checksum(rho_new(re)%matrix, pos=.TRUE.)
     647         136 :          IF (unit_nr > 0) THEN
     648          68 :             WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
     649             :          END IF
     650         136 :          CALL dbcsr_binary_write(rho_new(re)%matrix, file_name)
     651         136 :          IF (history) THEN
     652             :             WRITE (file_name, '(A,I0,A)') &
     653           2 :                TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
     654             :          ELSE
     655         134 :             WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
     656             :          END IF
     657         136 :          cs_pos = dbcsr_checksum(rho_new(im)%matrix, pos=.TRUE.)
     658         136 :          IF (unit_nr > 0) THEN
     659          68 :             WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
     660             :          END IF
     661         234 :          CALL dbcsr_binary_write(rho_new(im)%matrix, file_name)
     662             :       END DO
     663             : 
     664          98 :       CALL timestop(handle)
     665             : 
     666          98 :    END SUBROUTINE write_rt_p_to_restart
     667             : 
     668             : ! **************************************************************************************************
     669             : !> \brief Collocation of the current and printing of it in a cube file
     670             : !> \param qs_env ...
     671             : !> \param P_im ...
     672             : !> \param dft_section ...
     673             : !> \param spin ...
     674             : !> \param nspin ...
     675             : !> \author Samuel Andermatt (06.15)
     676             : ! **************************************************************************************************
     677          44 :    SUBROUTINE rt_current(qs_env, P_im, dft_section, spin, nspin)
     678             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     679             :       TYPE(dbcsr_type), POINTER                          :: P_im
     680             :       TYPE(section_vals_type), POINTER                   :: dft_section
     681             :       INTEGER                                            :: spin, nspin
     682             : 
     683             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'rt_current'
     684             : 
     685             :       CHARACTER(len=1)                                   :: char_spin
     686             :       CHARACTER(len=14)                                  :: ext
     687             :       CHARACTER(len=2)                                   :: sdir
     688             :       INTEGER                                            :: dir, handle, print_unit
     689          44 :       INTEGER, DIMENSION(:), POINTER                     :: stride(:)
     690             :       LOGICAL                                            :: mpi_io
     691             :       TYPE(cp_logger_type), POINTER                      :: logger
     692             :       TYPE(current_env_type)                             :: current_env
     693             :       TYPE(dbcsr_type), POINTER                          :: tmp, zero
     694             :       TYPE(particle_list_type), POINTER                  :: particles
     695             :       TYPE(pw_env_type), POINTER                         :: pw_env
     696             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     697             :       TYPE(pw_type)                                      :: gs, rs
     698             :       TYPE(qs_subsys_type), POINTER                      :: subsys
     699             : 
     700          44 :       CALL timeset(routineN, handle)
     701             : 
     702          44 :       logger => cp_get_default_logger()
     703          44 :       CALL get_qs_env(qs_env=qs_env, subsys=subsys, pw_env=pw_env)
     704          44 :       CALL qs_subsys_get(subsys, particles=particles)
     705          44 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     706             : 
     707          44 :       NULLIFY (zero, tmp)
     708          44 :       ALLOCATE (zero, tmp)
     709          44 :       CALL dbcsr_create(zero, template=P_im)
     710          44 :       CALL dbcsr_copy(zero, P_im)
     711          44 :       CALL dbcsr_set(zero, 0.0_dp)
     712          44 :       CALL dbcsr_create(tmp, template=P_im)
     713          44 :       CALL dbcsr_copy(tmp, P_im)
     714          44 :       IF (nspin == 1) THEN
     715          28 :          CALL dbcsr_scale(tmp, 0.5_dp)
     716             :       END IF
     717          44 :       current_env%gauge = -1
     718          44 :       current_env%gauge_init = .FALSE.
     719          44 :       CALL pw_pool_create_pw(auxbas_pw_pool, rs, use_data=REALDATA3D, in_space=REALSPACE)
     720          44 :       CALL pw_pool_create_pw(auxbas_pw_pool, gs, use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
     721             : 
     722             :       NULLIFY (stride)
     723          44 :       ALLOCATE (stride(3))
     724             : 
     725         176 :       DO dir = 1, 3
     726             : 
     727         132 :          CALL pw_zero(rs)
     728         132 :          CALL pw_zero(gs)
     729             : 
     730         132 :          CALL calculate_jrho_resp(zero, tmp, zero, zero, dir, dir, rs, gs, qs_env, current_env, retain_rsgrid=.TRUE.)
     731             : 
     732         528 :          stride = section_get_ivals(dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT%STRIDE")
     733             : 
     734         132 :          IF (dir == 1) THEN
     735          44 :             sdir = "-x"
     736          88 :          ELSEIF (dir == 2) THEN
     737          44 :             sdir = "-y"
     738             :          ELSE
     739          44 :             sdir = "-z"
     740             :          END IF
     741         132 :          WRITE (char_spin, "(I1)") spin
     742             : 
     743         132 :          ext = "-SPIN-"//char_spin//sdir//".cube"
     744         132 :          mpi_io = .TRUE.
     745             :          print_unit = cp_print_key_unit_nr(logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
     746             :                                            extension=ext, file_status="REPLACE", file_action="WRITE", &
     747         132 :                                            log_filename=.FALSE., mpi_io=mpi_io)
     748             : 
     749             :          CALL cp_pw_to_cube(rs, print_unit, "EMD current", particles=particles, stride=stride, &
     750         132 :                             mpi_io=mpi_io)
     751             : 
     752             :          CALL cp_print_key_finished_output(print_unit, logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
     753         176 :                                            mpi_io=mpi_io)
     754             : 
     755             :       END DO
     756             : 
     757          44 :       CALL pw_pool_give_back_pw(auxbas_pw_pool, rs)
     758          44 :       CALL pw_pool_give_back_pw(auxbas_pw_pool, gs)
     759             : 
     760          44 :       CALL dbcsr_deallocate_matrix(zero)
     761          44 :       CALL dbcsr_deallocate_matrix(tmp)
     762             : 
     763          44 :       DEALLOCATE (stride)
     764             : 
     765          44 :       CALL timestop(handle)
     766             : 
     767        3212 :    END SUBROUTINE rt_current
     768             : 
     769             : ! **************************************************************************************************
     770             : !> \brief Interface routine to trigger writing of results available from normal
     771             : !>        SCF. Can write MO-dependent and MO free results (needed for call from
     772             : !>        the linear scaling code)
     773             : !>        Update: trigger also some of prints for time-dependent runs
     774             : !> \param qs_env ...
     775             : !> \param rtp ...
     776             : !> \par History
     777             : !>      2022-11 Update [Guillaume Le Breton]
     778             : ! **************************************************************************************************
     779         434 :    SUBROUTINE write_available_results(qs_env, rtp)
     780             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     781             :       TYPE(rt_prop_type), POINTER                        :: rtp
     782             : 
     783             :       CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results'
     784             : 
     785             :       INTEGER                                            :: handle
     786             :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     787             : 
     788         434 :       CALL timeset(routineN, handle)
     789             : 
     790         434 :       CALL get_qs_env(qs_env, scf_env=scf_env)
     791         434 :       IF (rtp%linear_scaling) THEN
     792         182 :          CALL write_mo_free_results(qs_env)
     793             :       ELSE
     794         252 :          CALL write_mo_free_results(qs_env)
     795         252 :          CALL write_mo_dependent_results(qs_env, scf_env)
     796             :          ! Time-dependent MO print
     797         252 :          CALL write_rtp_mos_to_output_unit(qs_env, rtp)
     798         252 :          CALL write_rtp_mo_cubes(qs_env, rtp)
     799             :       END IF
     800             : 
     801         434 :       CALL timestop(handle)
     802             : 
     803         434 :    END SUBROUTINE write_available_results
     804             : 
     805             : ! **************************************************************************************************
     806             : !> \brief Print the field applied to the system. Either the electric
     807             : !>        field or the vector potential depending on the gauge used
     808             : !> \param qs_env ...
     809             : !> \param dft_section ...
     810             : !> \par History
     811             : !>      2023-01  Created [Guillaume Le Breton]
     812             : ! **************************************************************************************************
     813           4 :    SUBROUTINE print_field_applied(qs_env, dft_section)
     814             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     815             :       TYPE(section_vals_type), POINTER                   :: dft_section
     816             : 
     817             :       CHARACTER(LEN=3), DIMENSION(3)                     :: rlab
     818             :       CHARACTER(LEN=default_path_length)                 :: filename
     819             :       INTEGER                                            :: i, output_unit, unit_nr
     820             :       REAL(kind=dp)                                      :: field(3), to_write(3)
     821             :       TYPE(cp_logger_type), POINTER                      :: logger
     822             :       TYPE(dft_control_type), POINTER                    :: dft_control
     823             : 
     824           4 :       NULLIFY (dft_control)
     825             : 
     826           4 :       logger => cp_get_default_logger()
     827           4 :       output_unit = cp_logger_get_default_io_unit(logger)
     828             : 
     829           4 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     830             : 
     831             :       unit_nr = cp_print_key_unit_nr(logger, dft_section, &
     832           4 :                                      "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat")
     833             : 
     834           4 :       IF (output_unit > 0) THEN
     835           2 :          IF (unit_nr /= output_unit) THEN
     836           2 :             INQUIRE (UNIT=unit_nr, NAME=filename)
     837             :             WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
     838           2 :                "FIELD", "The field applied is written to the file:", &
     839           4 :                TRIM(filename)
     840             :             WRITE (UNIT=unit_nr, FMT="(/,(T2,A,T40,I6))") &
     841           2 :                "Real time propagation step:", qs_env%sim_step
     842             :          ELSE
     843           0 :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED"
     844             :          END IF
     845             : 
     846           8 :          rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
     847             : 
     848           2 :          IF (dft_control%apply_efield_field) THEN
     849           0 :             WRITE (unit_nr, "(T3,A)") "Electric Field (LG) in atomic units:"
     850           0 :             CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
     851           0 :             to_write = field
     852           2 :          ELSE IF (dft_control%apply_vector_potential) THEN
     853           0 :             WRITE (unit_nr, "(T3,A)") "Vector potential (VG) in atomic units:"
     854           0 :             to_write = dft_control%rtp_control%vec_pot
     855             :          ELSE
     856           2 :             WRITE (unit_nr, "(T3,A)") "No electric field applied"
     857           2 :             to_write = 0._dp
     858             :          END IF
     859             : 
     860             :          WRITE (unit_nr, "(T5,3(A,A,E16.8,1X))") &
     861           8 :             (TRIM(rlab(i)), "=", to_write(i), i=1, 3)
     862             :       END IF
     863           4 :    END SUBROUTINE print_field_applied
     864             : 
     865             : END MODULE rt_propagation_output

Generated by: LCOV version 1.15