LCOV - code coverage report
Current view: top level - src/emd - rt_propagation_output.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 363 369 98.4 %
Date: 2024-04-20 06:29:22 Functions: 11 11 100.0 %

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

Generated by: LCOV version 1.15