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

Generated by: LCOV version 1.15