LCOV - code coverage report
Current view: top level - src/emd - rt_bse.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:2fce0f8) Lines: 366 445 82.2 %
Date: 2024-12-21 06:28:57 Functions: 14 17 82.4 %

          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 Routines for the propagation via RT-BSE method.
      10             : !> \note  The control is handed directly from cp2k_runs
      11             : !> \author Stepan Marek (12.23)
      12             : ! **************************************************************************************************
      13             : 
      14             : MODULE rt_bse
      15             :    USE cp_control_types, ONLY: dft_control_type, &
      16             :                                rtp_control_type
      17             :    USE qs_environment_types, ONLY: get_qs_env, &
      18             :                                    qs_environment_type
      19             :    USE force_env_types, ONLY: force_env_type
      20             :    USE qs_mo_types, ONLY: mo_set_type, &
      21             :                           init_mo_set
      22             :    USE atomic_kind_types, ONLY: atomic_kind_type
      23             :    USE qs_kind_types, ONLY: qs_kind_type
      24             :    USE particle_types, ONLY: particle_type
      25             :    USE rt_propagation_types, ONLY: get_rtp, &
      26             :                                    rt_prop_type
      27             :    USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
      28             :    USE cp_fm_types, ONLY: cp_fm_type, &
      29             :                           cp_fm_p_type, &
      30             :                           cp_fm_to_fm, &
      31             :                           cp_fm_create, &
      32             :                           cp_fm_set_all, &
      33             :                           cp_fm_release, &
      34             :                           cp_fm_get_element, &
      35             :                           cp_fm_set_element, &
      36             :                           cp_fm_get_info, &
      37             :                           cp_fm_write_unformatted, &
      38             :                           cp_fm_read_unformatted, &
      39             :                           cp_fm_write_formatted
      40             :    USE cp_cfm_types, ONLY: cp_cfm_type, &
      41             :                            cp_cfm_p_type, &
      42             :                            cp_fm_to_cfm, &
      43             :                            cp_cfm_to_cfm, &
      44             :                            cp_cfm_to_fm, &
      45             :                            cp_cfm_create, &
      46             :                            cp_cfm_get_info, &
      47             :                            cp_cfm_release
      48             :    USE kinds, ONLY: dp, &
      49             :                     default_path_length
      50             :    USE cp_dbcsr_api, ONLY: dbcsr_p_type, &
      51             :                            dbcsr_type, &
      52             :                            dbcsr_print, &
      53             :                            dbcsr_has_symmetry, &
      54             :                            dbcsr_desymmetrize, &
      55             :                            dbcsr_create, &
      56             :                            dbcsr_release, &
      57             :                            dbcsr_copy, &
      58             :                            dbcsr_scale, &
      59             :                            dbcsr_add, &
      60             :                            dbcsr_set, &
      61             :                            dbcsr_clear, &
      62             :                            dbcsr_setname, &
      63             :                            dbcsr_iterator_type, &
      64             :                            dbcsr_iterator_start, &
      65             :                            dbcsr_iterator_stop, &
      66             :                            dbcsr_iterator_blocks_left, &
      67             :                            dbcsr_iterator_next_block, &
      68             :                            dbcsr_put_block, &
      69             :                            dbcsr_reserve_blocks, &
      70             :                            dbcsr_get_num_blocks, &
      71             :                            dbcsr_get_block_p, &
      72             :                            dbcsr_get_info
      73             :    USE OMP_LIB, ONLY: omp_get_thread_num, &
      74             :                       omp_get_num_threads, &
      75             :                       omp_set_num_threads, &
      76             :                       omp_get_max_threads
      77             :    USE dbt_api, ONLY: dbt_clear, &
      78             :                       dbt_contract, &
      79             :                       dbt_copy_matrix_to_tensor, &
      80             :                       dbt_copy_tensor_to_matrix
      81             :    USE libint_2c_3c, ONLY: libint_potential_type
      82             :    USE mp2_ri_2c, ONLY: RI_2c_integral_mat
      83             :    USE qs_tensors, ONLY: neighbor_list_3c_destroy, &
      84             :                          build_2c_integrals, &
      85             :                          build_2c_neighbor_lists
      86             :    USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type, &
      87             :                                      release_neighbor_list_sets
      88             :    USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm, &
      89             :                                   copy_fm_to_dbcsr, &
      90             :                                   dbcsr_allocate_matrix_set, &
      91             :                                   dbcsr_deallocate_matrix_set, &
      92             :                                   cp_dbcsr_sm_fm_multiply, &
      93             :                                   copy_cfm_to_dbcsr, &
      94             :                                   copy_dbcsr_to_cfm
      95             :    USE cp_fm_basic_linalg, ONLY: cp_fm_scale, &
      96             :                                  cp_fm_invert, &
      97             :                                  cp_fm_trace, &
      98             :                                  cp_fm_transpose, &
      99             :                                  cp_fm_norm, &
     100             :                                  cp_fm_column_scale, &
     101             :                                  cp_fm_scale_and_add
     102             :    USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add, &
     103             :                                   cp_cfm_scale, &
     104             :                                   cp_cfm_transpose, &
     105             :                                   cp_cfm_norm, &
     106             :                                   cp_cfm_trace, &
     107             :                                   cp_cfm_column_scale
     108             :    USE cp_fm_diag, ONLY: cp_fm_geeig
     109             :    USE cp_cfm_diag, ONLY: cp_cfm_geeig
     110             :    USE parallel_gemm_api, ONLY: parallel_gemm
     111             :    USE qs_moments, ONLY: build_local_moment_matrix, &
     112             :                          build_berry_moment_matrix
     113             :    USE moments_utils, ONLY: get_reference_point
     114             :    USE qs_mo_io, ONLY: read_mo_set_from_restart
     115             :    USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
     116             :    USE force_env_methods, ONLY: force_env_calc_energy_force
     117             :    USE qs_energy_init, ONLY: qs_energies_init
     118             :    USE qs_energy_utils, ONLY: qs_energies_properties
     119             :    USE efield_utils, ONLY: make_field
     120             :    USE rt_propagator_init, ONLY: rt_initialize_rho_from_mos
     121             :    USE rt_propagation_methods, ONLY: s_matrices_create
     122             :    USE gw_large_cell_Gamma, ONLY: compute_3c_integrals
     123             :    USE gw_integrals, ONLY: build_3c_integral_block
     124             :    USE matrix_exp, ONLY: taylor_full_complex
     125             :    USE message_passing, ONLY: mp_para_env_type
     126             :    USE input_section_types, ONLY: section_vals_get_subs_vals, &
     127             :                                   section_vals_type
     128             :    USE cell_types, ONLY: cell_type
     129             :    USE input_constants, ONLY: rtp_bse_ham_ks, &
     130             :                               rtp_bse_ham_g0w0, &
     131             :                               use_mom_ref_coac, &
     132             :                               use_mom_ref_zero, &
     133             :                               do_taylor, &
     134             :                               do_bch, &
     135             :                               do_exact, &
     136             :                               use_rt_restart
     137             :    USE rt_bse_types, ONLY: rtbse_env_type, &
     138             :                            create_rtbse_env, &
     139             :                            release_rtbse_env, &
     140             :                            multiply_cfm_fm, &
     141             :                            multiply_fm_cfm
     142             :    USE rt_bse_io, ONLY: output_moments, &
     143             :                         read_moments, &
     144             :                         read_field, &
     145             :                         output_moments_ft, &
     146             :                         output_polarizability, &
     147             :                         output_field, &
     148             :                         output_mos_contravariant, &
     149             :                         read_restart, &
     150             :                         output_restart, &
     151             :                         print_timestep_info, &
     152             :                         print_etrs_info_header, &
     153             :                         print_etrs_info, &
     154             :                         print_rtbse_header_info
     155             :    USE qs_ks_types, ONLY: qs_ks_env_type
     156             :    USE qs_integrate_potential_product, ONLY: integrate_v_rspace
     157             :    USE qs_collocate_density, ONLY: calculate_rho_elec
     158             :    USE cp_files, ONLY: open_file, &
     159             :                        file_exists, &
     160             :                        close_file
     161             :    USE physcon, ONLY: femtoseconds
     162             :    USE mathconstants, ONLY: twopi
     163             : 
     164             : #include "../base/base_uses.f90"
     165             : 
     166             :    IMPLICIT NONE
     167             : 
     168             :    PRIVATE
     169             : 
     170             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse"
     171             : 
     172             :    #:include "rt_bse_macros.fypp"
     173             : 
     174             :    PUBLIC :: run_propagation_bse
     175             : 
     176             :    INTERFACE get_sigma
     177             :       MODULE PROCEDURE get_sigma_dbcsr, get_sigma_real, get_sigma_complex
     178             :    END INTERFACE
     179             : 
     180             : CONTAINS
     181             : 
     182             : ! **************************************************************************************************
     183             : !> \brief Runs the electron-only real time BSE propagation
     184             : !> \param qs_env Quickstep environment data, containing the config from input files
     185             : !> \param force_env Force environment data, entry point of the calculation
     186             : ! **************************************************************************************************
     187          12 :    SUBROUTINE run_propagation_bse(qs_env, force_env)
     188             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     189             :       TYPE(force_env_type), POINTER                      :: force_env
     190             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'run_propagation_bse'
     191             :       TYPE(rtbse_env_type), POINTER                      :: rtbse_env
     192             :       INTEGER                                            :: i, j, k, handle
     193             :       LOGICAL                                            :: converged
     194             :       REAL(kind=dp)                                      :: metric, enum_re, enum_im, &
     195             :                                                             idempotence_dev, a_metric_1, a_metric_2
     196             : 
     197          12 :       CALL timeset(routineN, handle)
     198             : 
     199             :       ! Run the initial SCF calculation / read SCF restart information
     200          12 :       CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., consistent_energies=.FALSE.)
     201             : 
     202             :       ! Allocate all persistant storage and read input that does not need further processing
     203          12 :       CALL create_rtbse_env(rtbse_env, qs_env, force_env)
     204             : 
     205          12 :       IF (rtbse_env%unit_nr > 0) CALL print_rtbse_header_info(rtbse_env)
     206             : 
     207             :       ! Initialize non-trivial values
     208             :       !  - calculates the moment operators
     209             :       !  - populates the initial density matrix
     210             :       !     - reads the restart density if requested
     211             :       !  - reads the field and moment trace from previous runs
     212             :       !  - populates overlap and inverse overlap matrices
     213             :       !  - calculates/populates the G0W0/KS Hamiltonian, respectively
     214             :       !  - calculates the Hartree reference potential
     215             :       !  - calculates the COHSEX reference self-energy
     216             :       !  - prints some info about loaded files into the output
     217          12 :       CALL initialize_rtbse_env(rtbse_env)
     218             : 
     219             :       ! Setup the time based on the starting step
     220             :       ! Assumes identical dt between two runs
     221          12 :       rtbse_env%sim_time = REAL(rtbse_env%sim_start, dp)*rtbse_env%sim_dt
     222             :       ! Output 0 time moments and field
     223          12 :       IF (.NOT. rtbse_env%restart_extracted) THEN
     224           8 :          CALL output_field(rtbse_env)
     225           8 :          CALL output_moments(rtbse_env, rtbse_env%rho)
     226             :       END IF
     227             : 
     228             :       ! Do not apply the delta kick if we are doing a restart calculation
     229          12 :       IF (rtbse_env%dft_control%rtp_control%apply_delta_pulse .AND. (.NOT. rtbse_env%restart_extracted)) THEN
     230           4 :          CALL apply_delta_pulse(rtbse_env)
     231             :       END IF
     232             : 
     233             :       ! ********************** Start the time loop **********************
     234         620 :       DO i = rtbse_env%sim_start, rtbse_env%sim_nsteps
     235             : 
     236             :          ! Update the simulation time
     237         608 :          rtbse_env%sim_time = REAL(i, dp)*rtbse_env%sim_dt
     238         608 :          rtbse_env%sim_step = i
     239             :          ! Start the enforced time reversal method
     240             :          ! This method determines the density matrix at time (t+dt) by guessing the effective Hamiltonian at (t + dt)
     241             :          ! and using the Hamiltonian at time (t), it propagates density from time (t) while ensuring that the density
     242             :          ! at (t + dt/2) is the same for both forward and backwards propagation. Then, density at (t + dt) is
     243             :          ! used to calculate the new Hamiltonian at (t+dt), which is then used to get the new propagator, and so on
     244             :          ! until the density matrix does not change within certain limit
     245             :          ! Pseudocode of the algorithm
     246             :          !      rho_M = exp(-i S^(-1) H[rho(t)] dt/2) rho(t) exp(i H[rho(t)] S^(-1) dt/2)
     247             :          !      rho(t+dt, 0) = rho_M
     248             :          !      for j in 1,max_self_iter
     249             :          !              rho(t+dt,j) = exp(- i S^(-1) H[rho(t+dt,j-1)] dt/2) rho_M exp(i H [rho(t+dt,j-1)] S^(-1) dt/2)
     250             :          !              if ||rho(t+dt,j) - rho(t+dt,j-1)|| < epsilon
     251             :          !                      break
     252             :          ! *************** Determine rho_M ***************
     253             :          ! Update the effective Hamiltonian
     254             :          !  - sets the correct external field, updates Hartree and COHSEX terms to current density
     255         608 :          CALL update_effective_ham(rtbse_env, rtbse_env%rho)
     256             :          ! Convert the effective hamiltonian into the exponential
     257             :          ! TODO : Store the result more explicitly - explicit dependence on input/output matrices
     258         608 :          CALL ham_to_exp(rtbse_env)
     259             :          ! Propagate the density to mid-point
     260         608 :          CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rtbse_env%rho, rtbse_env%rho_M)
     261             :          ! In initial iteration, copy rho_M to rho_new_last - that is our initial guess
     262        1216 :          DO j = 1, rtbse_env%n_spin
     263        1216 :             CALL cp_cfm_to_cfm(rtbse_env%rho_M(j), rtbse_env%rho_new_last(j))
     264             :          END DO
     265             :          ! *********** Start the self-consistent loop ************************
     266             :          ! The guess Hamiltonian uses the guess density and field from the next timestep
     267         608 :          rtbse_env%sim_time = REAL(i + 1, dp)*rtbse_env%sim_dt
     268         608 :          rtbse_env%sim_step = i + 1
     269         608 :          converged = .FALSE.
     270             :          ! Prints the grid for etrs info dumping - only for log level > medium
     271         608 :          CALL print_etrs_info_header(rtbse_env)
     272        1824 :          DO k = 1, rtbse_env%etrs_max_iter
     273        1824 :             CALL update_effective_ham(rtbse_env, rtbse_env%rho_new_last)
     274        1824 :             CALL ham_to_exp(rtbse_env)
     275        1824 :             CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rtbse_env%rho_M, rtbse_env%rho_new)
     276             :             ! *** Self-consistency check ***
     277        1824 :             metric = rho_metric(rtbse_env%rho_new, rtbse_env%rho_new_last, rtbse_env%n_spin)
     278             :             ! ETRS info - only for log level > medium
     279        1824 :             CALL print_etrs_info(rtbse_env, k, metric)
     280        1824 :             IF (metric < rtbse_env%etrs_threshold) THEN
     281             :                converged = .TRUE.
     282             :                EXIT
     283             :             ELSE
     284             :                ! Copy rho_new to rho_new_last
     285        2432 :                DO j = 1, rtbse_env%n_spin
     286             :                   ! Leaving for free convergence
     287        2432 :                   CALL cp_cfm_to_cfm(rtbse_env%rho_new(j), rtbse_env%rho_new_last(j))
     288             :                END DO
     289             :             END IF
     290             :          END DO
     291         608 :          CALL get_electron_number(rtbse_env, rtbse_env%rho_new, enum_re, enum_im)
     292             :          IF (.FALSE.) THEN
     293             :             ! Not all of these are used, but they are all good metrics to check the convergence in problematic cases
     294             :             ! TODO : Allow for conditional warning
     295             :             CALL get_idempotence_deviation(rtbse_env, rtbse_env%rho_new, idempotence_dev)
     296             :             DO j = 1, rtbse_env%n_spin
     297             :                CALL cp_cfm_to_fm(rtbse_env%sigma_SEX(j), rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
     298             :                CALL antiherm_metric(real_fm=rtbse_env%real_workspace(1), imag_fm=rtbse_env%real_workspace(2), &
     299             :                                     workspace=rtbse_env%rho_workspace, metric=a_metric_1)
     300             :                CALL antiherm_metric(real_fm=rtbse_env%hartree_curr(j), &
     301             :                                     workspace=rtbse_env%rho_workspace, metric=a_metric_2)
     302             :             END DO
     303             :          END IF
     304         608 :          CALL print_timestep_info(rtbse_env, i, metric, enum_re, k)
     305         608 :          CPASSERT(converged)
     306        1216 :          DO j = 1, rtbse_env%n_spin
     307        1216 :             CALL cp_cfm_to_cfm(rtbse_env%rho_new(j), rtbse_env%rho(j))
     308             :          END DO
     309             :          ! Print the updated field
     310         608 :          CALL output_field(rtbse_env)
     311             :          ! If needed, print out the density matrix in MO basis
     312         608 :          CALL output_mos_contravariant(rtbse_env, rtbse_env%rho, rtbse_env%rho_section)
     313             :          ! Also handles outputting to memory
     314         608 :          CALL output_moments(rtbse_env, rtbse_env%rho)
     315             :          ! Output restart files, so that the restart starts at the following time index
     316        1228 :          CALL output_restart(rtbse_env, rtbse_env%rho, i + 1)
     317             :       END DO
     318             :       ! ********************** End the time loop **********************
     319             : 
     320             :       ! Carry out the FT
     321          12 :       CALL output_moments_ft(rtbse_env)
     322          12 :       CALL output_polarizability(rtbse_env)
     323             : 
     324             :       ! Deallocate everything
     325          12 :       CALL release_rtbse_env(rtbse_env)
     326             : 
     327          12 :       CALL timestop(handle)
     328          12 :    END SUBROUTINE run_propagation_bse
     329             : 
     330             : ! **************************************************************************************************
     331             : !> \brief Calculates the initial values, based on restart/scf density, and other non-trivial values
     332             : !> \param rtbse_env RT-BSE environment
     333             : !> \param qs_env Quickstep environment (needed for reference to previous calculations)
     334             : !> \author Stepan Marek (09.24)
     335             : ! **************************************************************************************************
     336          12 :    SUBROUTINE initialize_rtbse_env(rtbse_env)
     337             :       TYPE(rtbse_env_type), POINTER                :: rtbse_env
     338             :       CHARACTER(len=*), PARAMETER                  :: routineN = "initialize_rtbse_env"
     339             :       TYPE(post_scf_bandstructure_type), POINTER   :: bs_env
     340          12 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER    :: moments_dbcsr_p
     341          12 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER    :: matrix_s
     342          12 :       REAL(kind=dp), DIMENSION(:), POINTER         :: custom_ref_point, &
     343          12 :                                                       occupations
     344             :       REAL(kind=dp), DIMENSION(3)                  :: rpoint
     345             :       INTEGER                                      :: i, k, handle
     346             : 
     347          12 :       CALL timeset(routineN, handle)
     348             : 
     349             :       ! Get pointers to parameters from qs_env
     350          12 :       CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env, matrix_s=matrix_s)
     351             : 
     352             :       ! ****** START MOMENTS OPERATOR CALCULATION
     353             :       ! Construct moments from dbcsr
     354             :       NULLIFY (moments_dbcsr_p)
     355          48 :       ALLOCATE (moments_dbcsr_p(3))
     356          48 :       DO k = 1, 3
     357             :          ! Make sure the pointer is empty
     358          36 :          NULLIFY (moments_dbcsr_p(k)%matrix)
     359             :          ! Allocate a new matrix that the pointer points to
     360          36 :          ALLOCATE (moments_dbcsr_p(k)%matrix)
     361             :          ! Create the matrix storage - matrix copies the structure of overlap matrix
     362          48 :          CALL dbcsr_copy(moments_dbcsr_p(k)%matrix, matrix_s(1)%matrix)
     363             :       END DO
     364             :       ! Run the moment calculation
     365             :       ! check for presence to prevent memory errors
     366             :       NULLIFY (custom_ref_point)
     367          48 :       ALLOCATE (custom_ref_point(3), source=0.0_dp)
     368          12 :       rpoint(:) = 0.0_dp
     369          12 :       CALL get_reference_point(rpoint, qs_env=rtbse_env%qs_env, reference=use_mom_ref_coac, ref_point=custom_ref_point)
     370          12 :       CALL build_local_moment_matrix(rtbse_env%qs_env, moments_dbcsr_p, 1, rpoint)
     371             :       ! Copy to full matrix
     372          48 :       DO k = 1, 3
     373             :          ! Again, matrices are created from overlap template
     374          48 :          CALL copy_dbcsr_to_fm(moments_dbcsr_p(k)%matrix, rtbse_env%moments(k))
     375             :       END DO
     376             :       ! Now, repeat without reference point to get the moments for field
     377          12 :       CALL get_reference_point(rpoint, qs_env=rtbse_env%qs_env, reference=use_mom_ref_zero, ref_point=custom_ref_point)
     378          12 :       DEALLOCATE (custom_ref_point)
     379          12 :       CALL build_local_moment_matrix(rtbse_env%qs_env, moments_dbcsr_p, 1, rpoint)
     380          48 :       DO k = 1, 3
     381          48 :          CALL copy_dbcsr_to_fm(moments_dbcsr_p(k)%matrix, rtbse_env%moments_field(k))
     382             :       END DO
     383             : 
     384             :       ! Now can deallocate dbcsr matrices
     385          48 :       DO k = 1, 3
     386          36 :          CALL dbcsr_release(moments_dbcsr_p(k)%matrix)
     387          48 :          DEALLOCATE (moments_dbcsr_p(k)%matrix)
     388             :       END DO
     389          12 :       DEALLOCATE (moments_dbcsr_p)
     390             :       ! ****** END MOMENTS OPERATOR CALCULATION
     391             : 
     392             :       ! ****** START INITIAL DENSITY MATRIX CALCULATION
     393             :       ! Get the rho from fm_MOS
     394             :       ! Uses real orbitals only - no kpoints
     395          36 :       ALLOCATE (occupations(rtbse_env%n_ao))
     396             :       ! Iterate over both spins
     397          24 :       DO i = 1, rtbse_env%n_spin
     398          36 :          occupations(:) = 0.0_dp
     399          24 :          occupations(1:rtbse_env%n_occ(i)) = 1.0_dp
     400             :          ! Create real part
     401          12 :          CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(i), rtbse_env%real_workspace(1))
     402          12 :          CALL cp_fm_column_scale(rtbse_env%real_workspace(1), occupations)
     403             :          CALL parallel_gemm("N", "T", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     404             :                             1.0_dp, rtbse_env%real_workspace(1), bs_env%fm_mo_coeff_Gamma(i), &
     405          12 :                             0.0_dp, rtbse_env%real_workspace(2))
     406             :          ! Sets imaginary part to zero
     407          12 :          CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%rho(i))
     408             :          ! Save the reference value for the case of delta kick
     409          24 :          CALL cp_cfm_to_cfm(rtbse_env%rho(i), rtbse_env%rho_orig(i))
     410             :       END DO
     411          12 :       DEALLOCATE (occupations)
     412             :       ! If the restart field is provided, overwrite rho from restart
     413          12 :       IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
     414           4 :          CALL read_restart(rtbse_env)
     415             :       END IF
     416             :       ! ****** END INITIAL DENSITY MATRIX CALCULATION
     417             :       ! ****** START MOMENTS TRACE LOAD
     418             :       ! The moments are only loaded if the consistent save files exist
     419          12 :       CALL read_moments(rtbse_env)
     420             :       ! ****** END MOMENTS TRACE LOAD
     421             : 
     422             :       ! ****** START FIELD TRACE LOAD
     423             :       ! The moments are only loaded if the consistent save files exist
     424          12 :       CALL read_field(rtbse_env)
     425             :       ! ****** END FIELD TRACE LOAD
     426             : 
     427             :       ! ****** START OVERLAP + INVERSE OVERLAP CALCULATION
     428          12 :       CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, rtbse_env%S_fm)
     429          12 :       CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%S_cfm)
     430          12 :       CALL cp_fm_invert(rtbse_env%S_fm, rtbse_env%S_inv_fm)
     431             :       ! ****** END OVERLAP + INVERSE OVERLAP CALCULATION
     432             : 
     433             :       ! ****** START SINGLE PARTICLE HAMILTONIAN CALCULATION
     434          24 :       DO i = 1, rtbse_env%n_spin
     435          24 :          IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN
     436             :             ! G0W0 Hamiltonian
     437          12 :             CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(i), rtbse_env%real_workspace(1))
     438             :             ! NOTE : Gamma point is not always the zero k-point
     439             :             ! C * Lambda
     440          12 :             CALL cp_fm_column_scale(rtbse_env%real_workspace(1), bs_env%eigenval_G0W0(:, 1, i))
     441             :             ! C * Lambda * C^T
     442             :             CALL parallel_gemm("N", "T", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     443             :                                1.0_dp, rtbse_env%real_workspace(1), bs_env%fm_mo_coeff_Gamma(i), &
     444          12 :                                0.0_dp, rtbse_env%real_workspace(2))
     445             :             ! S * C * Lambda * C^T
     446             :             CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     447             :                                1.0_dp, rtbse_env%S_fm, rtbse_env%real_workspace(2), &
     448          12 :                                0.0_dp, rtbse_env%real_workspace(1))
     449             :             ! S * C * Lambda * C^T * S = H
     450             :             CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     451             :                                1.0_dp, rtbse_env%real_workspace(1), rtbse_env%S_fm, &
     452          12 :                                0.0_dp, rtbse_env%real_workspace(2))
     453          12 :             CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%ham_reference(i))
     454             :          ELSE
     455             :             ! KS Hamiltonian
     456           0 :             CALL cp_fm_to_cfm(msourcer=bs_env%fm_ks_Gamma(i), mtarget=rtbse_env%ham_reference(i))
     457             :          END IF
     458             :       END DO
     459             :       ! ****** END SINGLE PARTICLE HAMILTONIAN CALCULATION
     460             : 
     461             :       ! ****** START HARTREE POTENTIAL REFERENCE CALCULATION
     462             :       ! Calculate Coulomb RI elements, necessary for Hartree calculation
     463          12 :       CALL init_hartree(rtbse_env, rtbse_env%v_dbcsr)
     464          12 :       IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN
     465             :          ! In a non-HF calculation, copy the actual correlation part of the interaction
     466          12 :          CALL copy_fm_to_dbcsr(bs_env%fm_W_MIC_freq_zero, rtbse_env%w_dbcsr)
     467             :       ELSE
     468             :          ! In HF, correlation is set to zero
     469           0 :          CALL dbcsr_set(rtbse_env%w_dbcsr, 0.0_dp)
     470             :       END IF
     471             :       ! Add the Hartree to the screened_dbt tensor - now W = V + W^c
     472          12 :       CALL dbcsr_add(rtbse_env%w_dbcsr, rtbse_env%v_dbcsr, 1.0_dp, 1.0_dp)
     473          12 :       CALL dbt_copy_matrix_to_tensor(rtbse_env%w_dbcsr, rtbse_env%screened_dbt)
     474             :       ! Calculate the original Hartree potential
     475             :       ! Uses rho_orig - same as rho for initial run but different for continued run
     476          24 :       DO i = 1, rtbse_env%n_spin
     477          12 :          CALL get_hartree_local(rtbse_env, rtbse_env%rho_orig(i), rtbse_env%hartree_curr(i))
     478             :          ! Scaling by spin degeneracy
     479          12 :          CALL cp_fm_scale(rtbse_env%spin_degeneracy, rtbse_env%hartree_curr(i))
     480             :          ! Subtract the reference from the reference Hamiltonian
     481          12 :          CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr(i), mtarget=rtbse_env%ham_workspace(1))
     482             :          CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
     483          24 :                                    CMPLX(-1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1))
     484             :       END DO
     485             :       ! ****** END HARTREE POTENTIAL REFERENCE CALCULATION
     486             : 
     487             :       ! ****** START COHSEX REFERENCE CALCULATION
     488             :       ! Calculate the COHSEX starting energies
     489          12 :       IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN
     490             :          ! Subtract the v_xc from COH part of the self-energy, as V_xc is also not updated during the timestepping
     491          24 :          DO i = 1, rtbse_env%n_spin
     492             :             ! TODO : Allow no COH calculation for static screening
     493          12 :             CALL get_sigma(rtbse_env, rtbse_env%sigma_COH(i), -0.5_dp, rtbse_env%S_inv_fm)
     494             :             ! Copy and subtract from the complex reference hamiltonian
     495          12 :             CALL cp_fm_to_cfm(msourcer=rtbse_env%sigma_COH(i), mtarget=rtbse_env%ham_workspace(1))
     496             :             CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
     497          12 :                                       CMPLX(-1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1))
     498             :             ! Calculate exchange part - TODO : should this be applied for different spins? - TEST with O2 HF propagation?
     499             :             ! So far only closed shell tested
     500             :             ! Uses rho_orig - same as rho for initial run but different for continued run
     501          12 :             CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX(i), -1.0_dp, rtbse_env%rho_orig(i))
     502             :             ! Subtract from the complex reference Hamiltonian
     503             :             CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
     504          24 :                                       CMPLX(-1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(i))
     505             :          END DO
     506             :       ELSE
     507             :          ! KS Hamiltonian - use time-dependent Fock exchange
     508           0 :          DO i = 1, rtbse_env%n_spin
     509           0 :             CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX(i), -1.0_dp, rtbse_env%rho_orig(i))
     510             :             CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
     511           0 :                                       CMPLX(-1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(i))
     512             :          END DO
     513             :       END IF
     514             :       ! ****** END COHSEX REFERENCE CALCULATION
     515             : 
     516          12 :       CALL timestop(handle)
     517          12 :    END SUBROUTINE initialize_rtbse_env
     518             : 
     519             : ! **************************************************************************************************
     520             : !> \brief Custom reimplementation of the delta pulse routines
     521             : !> \param rtbse_env RT-BSE environment
     522             : !> \author Stepan Marek (09.24)
     523             : ! **************************************************************************************************
     524           4 :    SUBROUTINE apply_delta_pulse(rtbse_env)
     525             :       TYPE(rtbse_env_type), POINTER                :: rtbse_env
     526             :       CHARACTER(len=*), PARAMETER                  :: routineN = "apply_delta_pulse"
     527             :       REAL(kind=dp)                                :: intensity, metric
     528             :       REAL(kind=dp), DIMENSION(3)                  :: kvec
     529             :       INTEGER                                      :: i, k, handle
     530             : 
     531           4 :       CALL timeset(routineN, handle)
     532             : 
     533             :       ! Report application
     534           4 :       IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A27)') ' RTBSE| Applying delta pulse'
     535             :       ! Extra minus for the propagation of density
     536           4 :       intensity = -rtbse_env%dft_control%rtp_control%delta_pulse_scale
     537           4 :       metric = 0.0_dp
     538          16 :       kvec(:) = rtbse_env%dft_control%rtp_control%delta_pulse_direction(:)
     539           4 :       IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A38,E14.4E3,E14.4E3,E14.4E3)') &
     540           8 :          " RTBSE| Delta pulse elements (a.u.) : ", intensity*kvec(:)
     541             :       ! So far no spin dependence, but can be added by different structure of delta pulse
     542           4 :       CALL cp_fm_set_all(rtbse_env%real_workspace(1), 0.0_dp)
     543          16 :       DO k = 1, 3
     544             :          CALL cp_fm_scale_and_add(1.0_dp, rtbse_env%real_workspace(1), &
     545          16 :                                   kvec(k), rtbse_env%moments_field(k))
     546             :       END DO
     547             :       ! enforce hermiticity of the effective Hamiltonian
     548           4 :       CALL cp_fm_transpose(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
     549             :       CALL cp_fm_scale_and_add(0.5_dp, rtbse_env%real_workspace(1), &
     550           4 :                                0.5_dp, rtbse_env%real_workspace(2))
     551             :       ! Prepare the exponential/exponent for propagation
     552           4 :       IF (rtbse_env%mat_exp_method == do_bch) THEN
     553             :          ! Multiply by the S_inv matrix - in the classic ordering
     554             :          CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     555             :                             intensity, rtbse_env%S_inv_fm, rtbse_env%real_workspace(1), &
     556           4 :                             0.0_dp, rtbse_env%real_workspace(2))
     557           8 :          DO i = 1, rtbse_env%n_spin
     558             :             ! Sets real part to zero
     559           8 :             CALL cp_fm_to_cfm(msourcei=rtbse_env%real_workspace(2), mtarget=rtbse_env%ham_effective(i))
     560             :          END DO
     561           0 :       ELSE IF (rtbse_env%mat_exp_method == do_exact) THEN
     562           0 :          DO i = 1, rtbse_env%n_spin
     563           0 :             CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(1), mtarget=rtbse_env%ham_workspace(i))
     564             :             CALL cp_cfm_gexp(rtbse_env%ham_workspace(i), rtbse_env%S_cfm, rtbse_env%ham_effective(i), &
     565           0 :                              CMPLX(0.0, intensity, kind=dp), rtbse_env%rho_workspace)
     566             :          END DO
     567             :       END IF
     568             :       ! Propagate the density by the effect of the delta pulse
     569           4 :       CALL propagate_density(rtbse_env, rtbse_env%ham_effective, rtbse_env%rho, rtbse_env%rho_new)
     570           4 :       metric = rho_metric(rtbse_env%rho_new, rtbse_env%rho, rtbse_env%n_spin)
     571           4 :       IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, ('(A42,E38.8E3)')) " RTBSE| Metric difference after delta kick", metric
     572             :       ! Copy the new density to the old density
     573           8 :       DO i = 1, rtbse_env%n_spin
     574           8 :          CALL cp_cfm_to_cfm(rtbse_env%rho_new(i), rtbse_env%rho(i))
     575             :       END DO
     576             : 
     577           4 :       CALL timestop(handle)
     578           4 :    END SUBROUTINE apply_delta_pulse
     579             : ! **************************************************************************************************
     580             : !> \brief Determines the metric for the density matrix, used for convergence criterion
     581             : !> \param rho_new Array of new density matrices (one for each spin index)
     582             : !> \param rho_old Array of old density matrices (one for each spin index)
     583             : !> \param nspin Number of spin indices
     584             : !> \param workspace_opt Optionally provide external workspace to save some allocation time
     585             : ! **************************************************************************************************
     586       16938 :    FUNCTION rho_metric(rho_new, rho_old, nspin, workspace_opt) RESULT(metric)
     587             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER, INTENT(IN):: rho_new, &
     588             :                                                              rho_old
     589             :       INTEGER, INTENT(IN)                                :: nspin
     590             :       TYPE(cp_cfm_type), POINTER, OPTIONAL               :: workspace_opt
     591             :       TYPE(cp_cfm_type)                                  :: workspace
     592             :       REAL(kind=dp)                                      :: metric
     593       16938 :       REAL(kind=dp), DIMENSION(:), ALLOCATABLE           :: partial_metric
     594             :       INTEGER                                            :: j
     595             :       COMPLEX(kind=dp)                                   :: scale_factor
     596             : 
     597       50814 :       ALLOCATE (partial_metric(nspin))
     598             : 
     599             :       ! Only allocate/deallocate storage if required
     600       16938 :       IF (PRESENT(workspace_opt)) THEN
     601           0 :          workspace = workspace_opt
     602             :       ELSE
     603       16938 :          CALL cp_cfm_create(workspace, rho_new(1)%matrix_struct)
     604             :       END IF
     605       16938 :       scale_factor = 1.0
     606       33876 :       DO j = 1, nspin
     607       16938 :          CALL cp_cfm_to_cfm(rho_new(j), workspace)
     608             :          ! Get the difference in the resulting matrix
     609       16938 :          CALL cp_cfm_scale_and_add(scale_factor, workspace, -scale_factor, rho_old(j))
     610             :          ! Now, get the relevant number
     611       33876 :          partial_metric(j) = cp_cfm_norm(workspace, 'M')
     612             :       END DO
     613             :       metric = 0.0_dp
     614             :       ! For more than one spin, do Cartesian sum of the different spin norms
     615       33876 :       DO j = 1, nspin
     616       33876 :          metric = metric + partial_metric(j)*partial_metric(j)
     617             :       END DO
     618       16938 :       metric = SQRT(metric)
     619             :       ! Deallocate workspace
     620       16938 :       IF (.NOT. PRESENT(workspace_opt)) CALL cp_cfm_release(workspace)
     621       16938 :       DEALLOCATE (partial_metric)
     622       16938 :    END FUNCTION
     623             : 
     624             : ! **************************************************************************************************
     625             : !> \brief Determines the metric of the antihermitian part of the matrix
     626             : !> \param real_fm Real part of the full matrix
     627             : !> \param imag_fm Imaginary part of the full matrix
     628             : ! **************************************************************************************************
     629           0 :    SUBROUTINE antiherm_metric(real_fm, imag_fm, workspace, metric)
     630             :       TYPE(cp_fm_type), INTENT(IN)                      :: real_fm
     631             :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL            :: imag_fm
     632             :       REAL(kind=dp), INTENT(OUT)                        :: metric
     633             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER          :: workspace
     634             :       COMPLEX(kind=dp)                                  :: complex_one
     635             : 
     636             :       ! Get the complex and complex conjugate matrix
     637           0 :       IF (PRESENT(imag_fm)) THEN
     638           0 :          CALL cp_fm_to_cfm(real_fm, imag_fm, workspace(1))
     639             :       ELSE
     640           0 :          CALL cp_fm_to_cfm(msourcer=real_fm, mtarget=workspace(1))
     641             :       END IF
     642           0 :       CALL cp_cfm_transpose(workspace(1), "C", workspace(2))
     643             :       ! Subtract these, and get the metric
     644           0 :       complex_one = CMPLX(1.0, 0.0, kind=dp)
     645           0 :       CALL cp_cfm_scale_and_add(complex_one, workspace(1), -complex_one, workspace(2))
     646           0 :       metric = cp_cfm_norm(workspace(1), "M")
     647           0 :    END SUBROUTINE
     648             : 
     649             : ! **************************************************************************************************
     650             : !> \brief For Taylor and Exact exp_method, calculates the matrix exponential of the
     651             : !>        effective Hamiltonian. For BCH, calculates just the effective Hamiltonian. For other methods,
     652             : !>        aborts the execution, as they are not implemented yet.
     653             : !> \param rtbse_env Entry point of the calculation. Uses rho_workspace for Taylor and BCH. For exact,
     654             : !>                  uses complex_workspace, complex_ham, complex_s, real_eigvals and exp_eigvals.
     655             : !>                  Results are stored in ham_workspace.
     656             : ! **************************************************************************************************
     657        2432 :    SUBROUTINE ham_to_exp(rtbse_env)
     658             :       TYPE(rtbse_env_type)                               :: rtbse_env
     659             :       CHARACTER(len=*), PARAMETER                        :: routineN = "ham_to_exp"
     660             :       INTEGER                                            :: j, handle
     661        2432 :       CALL timeset(routineN, handle)
     662        4864 :       DO j = 1, rtbse_env%n_spin
     663        4864 :          IF (rtbse_env%mat_exp_method == do_bch) THEN
     664             :             ! In Taylor and BCH, we first evaluate the entire exponent and then evaluate exponential in series
     665             :             ! In order to produce correct result, need to remultiply by inverse overlap matrix
     666             :             CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     667             :                                  1.0_dp, rtbse_env%S_inv_fm, rtbse_env%ham_effective(j), &
     668        2432 :                                  0.0_dp, rtbse_env%rho_workspace(1))
     669             : 
     670             :             ! The evolution of density matrix is derived from the right multiplying term
     671             :             ! Imaginary part of the exponent = -real part of the matrix
     672        2432 :             CALL cp_cfm_scale(CMPLX(0.0, -rtbse_env%sim_dt/2, kind=dp), rtbse_env%rho_workspace(1))
     673             :             ! In BCH, exponential is not calculated explicitly, but the propagation is solved in series
     674        2432 :             CALL cp_cfm_to_cfm(rtbse_env%rho_workspace(1), rtbse_env%ham_workspace(j))
     675           0 :          ELSE IF (rtbse_env%mat_exp_method == do_exact) THEN
     676             :             CALL cp_cfm_gexp(rtbse_env%ham_effective(j), rtbse_env%S_cfm, rtbse_env%ham_workspace(j), &
     677           0 :                              CMPLX(0.0, -rtbse_env%sim_dt/2, kind=dp), rtbse_env%rho_workspace)
     678             :          ELSE
     679           0 :             CPABORT("Only BCH and Taylor matrix exponentiation implemented")
     680             :          END IF
     681             :       END DO
     682             : 
     683        2432 :       CALL timestop(handle)
     684        2432 :    END SUBROUTINE
     685             : ! **************************************************************************************************
     686             : !> \brief Updates the effective Hamiltonian, given a density matrix rho
     687             : !> \param rtbse_env Entry point of the calculation - contains current state of variables
     688             : !> \param qs_env QS env
     689             : !> \param rho Real and imaginary parts ( + spin) of the density at current time
     690             : ! **************************************************************************************************
     691        2432 :    SUBROUTINE update_effective_ham(rtbse_env, rho)
     692             :       TYPE(rtbse_env_type)                               :: rtbse_env
     693             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER           :: rho
     694             :       CHARACTER(len=*), PARAMETER                        :: routineN = "update_effective_ham"
     695             :       INTEGER                                            :: k, j, nspin, handle
     696             : 
     697        2432 :       CALL timeset(routineN, handle)
     698             :       ! Shorthand
     699        2432 :       nspin = rtbse_env%n_spin
     700             :       ! Reset the effective Hamiltonian to KS Hamiltonian + G0W0 - reference COHSEX - reference Hartree
     701        4864 :       DO j = 1, nspin
     702             :          ! Sets the imaginary part to zero
     703        4864 :          CALL cp_cfm_to_cfm(rtbse_env%ham_reference(j), rtbse_env%ham_effective(j))
     704             :       END DO
     705             :       ! Determine the field at current time
     706        2432 :       IF (rtbse_env%dft_control%apply_efield_field) THEN
     707         832 :          CALL make_field(rtbse_env%dft_control, rtbse_env%field, rtbse_env%sim_step, rtbse_env%sim_time)
     708             :       ELSE
     709             :          ! No field
     710        6400 :          rtbse_env%field(:) = 0.0_dp
     711             :       END IF
     712        4864 :       DO j = 1, nspin
     713        9728 :          DO k = 1, 3
     714             :             ! Minus sign due to charge of electrons
     715        7296 :             CALL cp_fm_to_cfm(msourcer=rtbse_env%moments_field(k), mtarget=rtbse_env%ham_workspace(1))
     716             :             CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), &
     717        9728 :                                       CMPLX(rtbse_env%field(k), 0.0, kind=dp), rtbse_env%ham_workspace(1))
     718             :          END DO
     719        2432 :          IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN
     720             :             ! Add the COH part - so far static but can be dynamic in principle throught the W updates
     721        2432 :             CALL get_sigma(rtbse_env, rtbse_env%sigma_COH(j), -0.5_dp, rtbse_env%S_inv_fm)
     722        2432 :             CALL cp_fm_to_cfm(msourcer=rtbse_env%sigma_COH(j), mtarget=rtbse_env%ham_workspace(1))
     723             :             CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), &
     724        2432 :                                       CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1))
     725             :          END IF
     726             :          ! Calculate the (S)EX part - based on provided rho
     727             :          ! iGW = - rho W
     728        2432 :          CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX(j), -1.0_dp, rho(j))
     729             :          CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), &
     730        2432 :                                    CMPLX(1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(j))
     731             :          ! Calculate Hartree potential
     732             :          ! Hartree potential is scaled by number of electrons in each MO - spin degeneracy
     733             :          CALL get_hartree_local(rtbse_env, rho(j), &
     734        2432 :                                 rtbse_env%hartree_curr(j))
     735        2432 :          CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr(j), mtarget=rtbse_env%ham_workspace(1))
     736             :          CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), &
     737        2432 :                                    CMPLX(rtbse_env%spin_degeneracy, 0.0, kind=dp), rtbse_env%ham_workspace(1))
     738             :          ! Enforce hermiticity of the effective Hamiltonian
     739             :          ! Important components without forced Hermiticity - moments matrix, sigma matrices, Hartree matrix
     740             :          ! single particle Ham
     741        2432 :          CALL cp_cfm_transpose(rtbse_env%ham_effective(j), 'C', rtbse_env%ham_workspace(1))
     742             :          CALL cp_cfm_scale_and_add(CMPLX(0.5, 0.0, kind=dp), rtbse_env%ham_effective(j), &
     743        4864 :                                    CMPLX(0.5, 0.0, kind=dp), rtbse_env%ham_workspace(1))
     744             :       END DO
     745        2432 :       CALL timestop(handle)
     746        2432 :    END SUBROUTINE update_effective_ham
     747             : 
     748             : ! **************************************************************************************************
     749             : !> \brief Does the BCH iterative determination of the exponential
     750             : !> \param propagator_matrix Matrix X which is to be exponentiated
     751             : !> \param target_matrix Matrix Y which the exponential acts upon
     752             : !> \param result_matrix Propagated matrix
     753             : !> \param workspace Matrices dedicated for work, 4 fm matrices with dimensions of X required
     754             : !> \param threshold_opt Optionally, a threshold under which the iteration is considered converged (default 1e-10)
     755             : !> \param max_iter_opt Optionally, maximum number of BCH iterations (default 20)
     756             : ! **************************************************************************************************
     757        4872 :    SUBROUTINE bch_propagate(propagator_matrix, target_matrix, result_matrix, workspace, threshold_opt, max_iter_opt)
     758             :       ! Array of complex propagator matrix X, such that
     759             :       ! the propagated matrix will follow Y' = e^X Y e^(-X), for each spin
     760             :       ! effect of e^(-X) is calculated - provide the X on the left hand side
     761             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER          :: propagator_matrix
     762             :       ! Matrix Y to be propagated into matrix Y'
     763             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER          :: target_matrix
     764             :       ! Matrix Y' is stored here on exit
     765             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER          :: result_matrix, workspace
     766             :       ! Threshold for the metric which decides when to truncate the BCH expansion
     767             :       REAL(kind=dp), OPTIONAL                           :: threshold_opt
     768             :       INTEGER, OPTIONAL                                 :: max_iter_opt
     769             :       CHARACTER(len=*), PARAMETER                       :: routineN = "bch_propagate"
     770             :       REAL(kind=dp)                                     :: threshold, prefactor, metric
     771             :       INTEGER                                           :: max_iter, i, n_spin, n_ao, k, &
     772             :                                                            w_stride, handle
     773             :       LOGICAL                                           :: converged
     774             :       CHARACTER(len=77)                                 :: error
     775             : 
     776        2436 :       CALL timeset(routineN, handle)
     777             : 
     778        2436 :       converged = .FALSE.
     779             : 
     780        2436 :       IF (PRESENT(threshold_opt)) THEN
     781        2436 :          threshold = threshold_opt
     782             :       ELSE
     783           0 :          threshold = 1.0e-10
     784             :       END IF
     785             : 
     786        2436 :       IF (PRESENT(max_iter_opt)) THEN
     787        2436 :          max_iter = max_iter_opt
     788             :       ELSE
     789             :          max_iter = 20
     790             :       END IF
     791             : 
     792        2436 :       n_spin = SIZE(target_matrix)
     793             :       n_ao = 0
     794        2436 :       CALL cp_cfm_get_info(target_matrix(1), nrow_global=n_ao)
     795        2436 :       w_stride = n_spin
     796             : 
     797             :       ! Initiate
     798        4872 :       DO i = 1, n_spin
     799        2436 :          CALL cp_cfm_to_cfm(target_matrix(i), result_matrix(i))
     800        4872 :          CALL cp_cfm_to_cfm(target_matrix(i), workspace(i))
     801             :       END DO
     802             : 
     803             :       ! Start the BCH iterations
     804             :       ! So far, no spin mixing terms
     805       15110 :       DO k = 1, max_iter
     806       15110 :          prefactor = 1.0_dp/REAL(k, kind=dp)
     807       30220 :          DO i = 1, n_spin
     808             :             CALL parallel_gemm("N", "N", n_ao, n_ao, n_ao, &
     809             :                                CMPLX(prefactor, 0.0, kind=dp), propagator_matrix(i), workspace(i), &
     810       15110 :                                CMPLX(0.0, 0.0, kind=dp), workspace(i + w_stride))
     811             :             CALL parallel_gemm("N", "C", n_ao, n_ao, n_ao, &
     812             :                                CMPLX(prefactor, 0.0, kind=dp), workspace(i), propagator_matrix(i), &
     813       15110 :                                CMPLX(1.0, 0.0, kind=dp), workspace(i + w_stride))
     814             :             ! Add to the result
     815             :             CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), result_matrix(i), &
     816       30220 :                                       CMPLX(1.0, 0.0, kind=dp), workspace(i + w_stride))
     817             :          END DO
     818       15110 :          metric = rho_metric(workspace(w_stride + 1:), workspace(1:w_stride), n_spin)
     819       15110 :          IF (metric <= threshold) THEN
     820             :             converged = .TRUE.
     821             :             EXIT
     822             :          ELSE
     823       25348 :             DO i = 1, n_spin
     824       25348 :                CALL cp_cfm_to_cfm(workspace(i + w_stride), workspace(i))
     825             :             END DO
     826             :          END IF
     827             :       END DO
     828        2436 :       IF (.NOT. converged) THEN
     829           0 :          WRITE (error, '(A35,E13.4E3,A16,E13.4E3)') "BCH did not converge, BCH Metric : ", &
     830           0 :             metric, "BCH Threshold : ", threshold
     831           0 :          CPABORT(error)
     832             :       END IF
     833             : 
     834        2436 :       CALL timestop(handle)
     835        2436 :    END SUBROUTINE bch_propagate
     836             : 
     837             : ! **************************************************************************************************
     838             : !> \brief Updates the density in rtbse_env, using the provided exponential
     839             : !>        The new density is saved to a different matrix, which enables for comparison of matrices
     840             : !> \param rtbse_env Entry point of the calculation - contains current state of variables
     841             : !> \param exponential Real and imaginary parts ( + spin) of the exponential propagator
     842             : ! **************************************************************************************************
     843        2436 :    SUBROUTINE propagate_density(rtbse_env, exponential, rho_old, rho_new)
     844             :       TYPE(rtbse_env_type)                               :: rtbse_env
     845             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER           :: exponential, &
     846             :                                                             rho_old, &
     847             :                                                             rho_new
     848             :       CHARACTER(len=*), PARAMETER                        :: routineN = "propagate_density"
     849             :       INTEGER                                            :: j, handle
     850             : 
     851        2436 :       CALL timeset(routineN, handle)
     852        2436 :       IF (rtbse_env%mat_exp_method == do_exact) THEN
     853             :          ! For these methods, exponential is explicitly constructed
     854           0 :          DO j = 1, rtbse_env%n_spin
     855             :             ! rho * (exp^dagger)
     856             :             CALL parallel_gemm("N", "C", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     857             :                                CMPLX(1.0, 0.0, kind=dp), rho_old(j), exponential(j), &
     858           0 :                                CMPLX(0.0, 0.0, kind=dp), rtbse_env%rho_workspace(1))
     859             :             ! exp * rho * (exp^dagger)
     860             :             CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     861             :                                CMPLX(1.0, 0.0, kind=dp), exponential(j), rtbse_env%rho_workspace(1), &
     862           0 :                                CMPLX(0.0, 0.0, kind=dp), rho_new(j))
     863             :          END DO
     864        2436 :       ELSE IF (rtbse_env%mat_exp_method == do_bch) THEN
     865             :          ! Same number of iterations as ETRS
     866             :          CALL bch_propagate(exponential, rho_old, rho_new, rtbse_env%rho_workspace, threshold_opt=rtbse_env%exp_accuracy, &
     867        2436 :                             max_iter_opt=rtbse_env%etrs_max_iter)
     868             :       ELSE
     869           0 :          CPABORT("Only BCH and exact matrix exponentiation implemented.")
     870             :       END IF
     871             : 
     872        2436 :       CALL timestop(handle)
     873        2436 :    END SUBROUTINE
     874             : 
     875             : ! **************************************************************************************************
     876             : !> \brief Outputs the number of electrons in the system from the density matrix
     877             : !> \note  Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3)
     878             : !> \param rtbse_env Entry point - rtbse environment
     879             : !> \param rho Density matrix in AO basis
     880             : !> \param electron_n_re Real number of electrons
     881             : !> \param electron_n_im Imaginary number of electrons, which can arise from numerical non-hermiticity
     882             : ! **************************************************************************************************
     883         608 :    SUBROUTINE get_electron_number(rtbse_env, rho, electron_n_re, electron_n_im)
     884             :       TYPE(rtbse_env_type)                               :: rtbse_env
     885             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER           :: rho
     886             :       REAL(kind=dp), INTENT(OUT)                         :: electron_n_re, electron_n_im
     887             :       COMPLEX(kind=dp)                                   :: electron_n_buffer
     888             :       INTEGER                                            :: j
     889             : 
     890         608 :       electron_n_re = 0.0_dp
     891         608 :       electron_n_im = 0.0_dp
     892         608 :       CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%rho_workspace(1))
     893        1216 :       DO j = 1, rtbse_env%n_spin
     894         608 :          CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rho(j), electron_n_buffer)
     895         608 :          electron_n_re = electron_n_re + REAL(electron_n_buffer, kind=dp)
     896        1216 :          electron_n_im = electron_n_im + REAL(AIMAG(electron_n_buffer), kind=dp)
     897             :       END DO
     898             :       ! Scale by spin degeneracy
     899         608 :       electron_n_re = electron_n_re*rtbse_env%spin_degeneracy
     900         608 :       electron_n_im = electron_n_im*rtbse_env%spin_degeneracy
     901         608 :    END SUBROUTINE get_electron_number
     902             : ! **************************************************************************************************
     903             : !> \brief Outputs the deviation from idempotence of density matrix
     904             : !> \note  Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3)
     905             : !> \param rtbse_env Entry point - rtbse environment
     906             : !> \param rho Density matrix in AO basis
     907             : !> \param electron_n_re Real number of electrons
     908             : !> \param electron_n_im Imaginary number of electrons, which can arise from numerical non-hermiticity
     909             : ! **************************************************************************************************
     910           0 :    SUBROUTINE get_idempotence_deviation(rtbse_env, rho, deviation_metric)
     911             :       TYPE(rtbse_env_type)                              :: rtbse_env
     912             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER          :: rho
     913             :       REAL(kind=dp), INTENT(OUT)                        :: deviation_metric
     914             :       COMPLEX(kind=dp)                                  :: buffer_1, buffer_2
     915             :       REAL(kind=dp)                                     :: buffer_dev
     916             :       INTEGER                                           :: j
     917             : 
     918           0 :       deviation_metric = 0.0_dp
     919           0 :       buffer_dev = 0.0_dp
     920             :       ! First, determine Tr(S * rho_re) + i Tr (S * rho_im)
     921           0 :       CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%rho_workspace(1))
     922           0 :       DO j = 1, rtbse_env%n_spin
     923           0 :          CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rho(j), buffer_1)
     924           0 :          buffer_dev = buffer_dev + REAL(ABS(buffer_1)*ABS(buffer_1), kind=dp)
     925             :       END DO
     926             :       ! Now, determine Tr(S * rho_re * S * rho_re) - Tr(S * rho_im * S * rho_im) + 2i Tr(S * rho_re * S * rho_im)
     927           0 :       DO j = 1, rtbse_env%n_spin
     928             :          ! S * rho
     929             :          CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     930             :                               1.0_dp, rtbse_env%S_fm, rho(j), &
     931           0 :                               0.0_dp, rtbse_env%rho_workspace(2))
     932             :          ! rho * S * rho
     933             :          CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
     934             :                             CMPLX(1.0, 0.0, kind=dp), rho(j), rtbse_env%rho_workspace(2), &
     935           0 :                             CMPLX(0.0, 0.0, kind=dp), rtbse_env%rho_workspace(3))
     936             :          ! Tr (S * rho * S * rho)
     937           0 :          CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rtbse_env%rho_workspace(3), buffer_2)
     938           0 :          deviation_metric = deviation_metric + REAL(ABS(buffer_2)*ABS(buffer_2), kind=dp)
     939             :       END DO
     940           0 :       deviation_metric = SQRT(deviation_metric) - SQRT(buffer_dev)
     941           0 :    END SUBROUTINE get_idempotence_deviation
     942             : 
     943             : ! **************************************************************************************************
     944             : !> \brief Calculates the self-energy by contraction of screened potential, for complex density
     945             : !> \note Can be used for both the Coulomb hole part and screened exchange part
     946             : !> \param rtbse_env Quickstep environment data, entry point of the calculation
     947             : !> \param sigma_cfm Pointer to the self-energy full matrix, which is overwritten by this routine
     948             : !> \param greens_cfm Pointer to the Green's function matrix, which is used as input data
     949             : !> \author Stepan Marek
     950             : !> \date 09.2024
     951             : ! **************************************************************************************************
     952        2444 :    SUBROUTINE get_sigma_complex(rtbse_env, sigma_cfm, prefactor_opt, greens_cfm)
     953             :       TYPE(rtbse_env_type)                               :: rtbse_env
     954             :       TYPE(cp_cfm_type)                                  :: sigma_cfm ! resulting self energy
     955             :       REAL(kind=dp), INTENT(IN), OPTIONAL                :: prefactor_opt
     956             :       TYPE(cp_cfm_type), INTENT(IN)                      :: greens_cfm ! matrix to contract with RI_W
     957             :       REAL(kind=dp)                                      :: prefactor
     958             : 
     959        2444 :       prefactor = 1.0_dp
     960        2444 :       IF (PRESENT(prefactor_opt)) prefactor = prefactor_opt
     961             : 
     962             :       ! Carry out the sigma part twice
     963             :       ! Real part
     964        2444 :       CALL cp_cfm_to_fm(msource=greens_cfm, mtargetr=rtbse_env%real_workspace(1))
     965        2444 :       CALL get_sigma(rtbse_env, rtbse_env%real_workspace(2), prefactor, rtbse_env%real_workspace(1))
     966        2444 :       CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%ham_workspace(1))
     967             :       ! Imaginary part
     968        2444 :       CALL cp_cfm_to_fm(msource=greens_cfm, mtargeti=rtbse_env%real_workspace(1))
     969        2444 :       CALL get_sigma(rtbse_env, rtbse_env%real_workspace(2), prefactor, rtbse_env%real_workspace(1))
     970        2444 :       CALL cp_fm_to_cfm(msourcei=rtbse_env%real_workspace(2), mtarget=sigma_cfm)
     971             :       ! Add the real part
     972             :       CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), sigma_cfm, &
     973        2444 :                                 CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1))
     974             : 
     975        2444 :    END SUBROUTINE get_sigma_complex
     976             : ! **************************************************************************************************
     977             : !> \brief Calculates the self-energy by contraction of screened potential, for complex density
     978             : !> \note Can be used for both the Coulomb hole part and screened exchange part
     979             : !> \param rtbse_env Quickstep environment data, entry point of the calculation
     980             : !> \param sigma_fm Pointer to the self-energy full matrix, which is overwritten by this routine
     981             : !> \param greens_fm Pointer to the Green's function matrix, which is used as input data
     982             : !> \author Stepan Marek
     983             : !> \date 09.2024
     984             : ! **************************************************************************************************
     985        7332 :    SUBROUTINE get_sigma_real(rtbse_env, sigma_fm, prefactor_opt, greens_fm)
     986             :       TYPE(rtbse_env_type)                               :: rtbse_env
     987             :       TYPE(cp_fm_type)                                   :: sigma_fm ! resulting self energy
     988             :       REAL(kind=dp), INTENT(IN), OPTIONAL                :: prefactor_opt
     989             :       TYPE(cp_fm_type), INTENT(IN)                       :: greens_fm ! matrix to contract with RI_W
     990             :       REAL(kind=dp)                                      :: prefactor
     991             : 
     992        7332 :       prefactor = 1.0_dp
     993        7332 :       IF (PRESENT(prefactor_opt)) prefactor = prefactor_opt
     994             : 
     995             :       ! Carry out the sigma part twice
     996             :       ! Convert to dbcsr
     997        7332 :       CALL copy_fm_to_dbcsr(greens_fm, rtbse_env%rho_dbcsr)
     998        7332 :       CALL get_sigma(rtbse_env, sigma_fm, prefactor, rtbse_env%rho_dbcsr)
     999             : 
    1000        7332 :    END SUBROUTINE get_sigma_real
    1001             : ! **************************************************************************************************
    1002             : !> \brief Calculates the self-energy by contraction of screened potential
    1003             : !> \note Can be used for both the Coulomb hole part and screened exchange part
    1004             : !> \param qs_env Quickstep environment data, entry point of the calculation
    1005             : !> \param greens_fm Pointer to the Green's function matrix, which is used as input data
    1006             : !> \param sigma_fm Pointer to the self-energy full matrix, which is overwritten by this routine
    1007             : !> \author Stepan Marek
    1008             : !> \date 01.2024
    1009             : ! **************************************************************************************************
    1010        7332 :    SUBROUTINE get_sigma_dbcsr(rtbse_env, sigma_fm, prefactor_opt, greens_dbcsr)
    1011             :       TYPE(rtbse_env_type)                               :: rtbse_env
    1012             :       TYPE(cp_fm_type)                                   :: sigma_fm ! resulting self energy
    1013             :       REAL(kind=dp), INTENT(IN), OPTIONAL                :: prefactor_opt
    1014             :       TYPE(dbcsr_type)                                   :: greens_dbcsr
    1015             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'get_sigma'
    1016             :       REAL(kind=dp)                                      :: prefactor
    1017             :       TYPE(dbcsr_type)                                   :: sigma_dbcsr
    1018             :       INTEGER                                            :: handle
    1019             : 
    1020        7332 :       CALL timeset(routineN, handle)
    1021             : 
    1022        7332 :       IF (PRESENT(prefactor_opt)) THEN
    1023        7332 :          prefactor = prefactor_opt
    1024             :       ELSE
    1025           0 :          prefactor = 1.0_dp
    1026             :       END IF
    1027             : 
    1028             :       ! Three-centre integrals are obtained from build_3c_integrals, from qs_tensors
    1029             :       ! These should use sparcity, while W and Sigma can be full matrices
    1030             :       ! The summation is carried out by dbt library - dbt_contract in dbt_api
    1031             :       ! The building of the tensors might be a bit hard, because it requires a lot of parallel information
    1032             :       ! Probably just use the tensors already present in bs_env? They seem to be mostly work tensors
    1033             :       ! Create by template
    1034             :       CALL dbt_contract(alpha=1.0_dp, &
    1035             :                         tensor_1=rtbse_env%screened_dbt, &
    1036             :                         tensor_2=rtbse_env%t_3c_w, &
    1037             :                         beta=0.0_dp, &
    1038             :                         tensor_3=rtbse_env%t_3c_work_RI_AO__AO, &
    1039             :                         contract_1=[2], notcontract_1=[1], map_1=[1], &
    1040        7332 :                         contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3])!,&
    1041             :       !filter_eps=bs_env%eps_filter)
    1042             :       ! t_work1 now contains B^P_(nu beta) = sum _ Q W _ (PQ) (iomega = 0) (Q| nu beta)
    1043             :       ! Next step is to convert the greens full matrix to dbcsr matrix
    1044        7332 :       CALL dbt_copy_matrix_to_tensor(greens_dbcsr, rtbse_env%greens_dbt)
    1045             :       ! Then contract it
    1046             :       ! no scaling applied - this has to be applied externally
    1047             :       CALL dbt_contract(alpha=1.0_dp, &
    1048             :                         tensor_1=rtbse_env%t_3c_work_RI_AO__AO, &
    1049             :                         tensor_2=rtbse_env%greens_dbt, &
    1050             :                         beta=0.0_dp, &
    1051             :                         tensor_3=rtbse_env%t_3c_work2_RI_AO__AO, &
    1052             :                         contract_1=[2], notcontract_1=[1, 3], map_1=[1, 3], &
    1053        7332 :                         contract_2=[2], notcontract_2=[1], map_2=[2])
    1054             :       ! workspace 2 now contains C ^ P _ (mu beta) sum _ nu B ^ P _ (nu beta) g _ (mu nu)
    1055             :       CALL dbt_contract(alpha=prefactor, &
    1056             :                         tensor_1=rtbse_env%t_3c_w, &
    1057             :                         tensor_2=rtbse_env%t_3c_work2_RI_AO__AO, &
    1058             :                         beta=0.0_dp, &
    1059             :                         tensor_3=rtbse_env%sigma_dbt, &
    1060             :                         contract_1=[1, 3], notcontract_1=[2], map_1=[1], &
    1061        7332 :                         contract_2=[1, 2], notcontract_2=[3], map_2=[2])!,&
    1062             :       !filter_eps=bs_env%eps_filter)
    1063             :       ! Finally, convert the COH tensor to matrix and then to fm matrix
    1064        7332 :       CALL dbcsr_create(sigma_dbcsr, name="sigma", template=rtbse_env%rho_dbcsr)
    1065        7332 :       CALL dbt_copy_tensor_to_matrix(rtbse_env%sigma_dbt, sigma_dbcsr)
    1066        7332 :       CALL copy_dbcsr_to_fm(sigma_dbcsr, sigma_fm)
    1067        7332 :       CALL dbcsr_release(sigma_dbcsr)
    1068             :       ! Clear workspaces - saves memory?
    1069        7332 :       CALL dbt_clear(rtbse_env%t_3c_work_RI_AO__AO)
    1070        7332 :       CALL dbt_clear(rtbse_env%t_3c_work2_RI_AO__AO)
    1071        7332 :       CALL dbt_clear(rtbse_env%sigma_dbt)
    1072        7332 :       CALL dbt_clear(rtbse_env%greens_dbt)
    1073        7332 :       CALL timestop(handle)
    1074             : 
    1075        7332 :    END SUBROUTINE get_sigma_dbcsr
    1076             : ! **************************************************************************************************
    1077             : !> \brief Creates the RI matrix and populates it with correct values
    1078             : !> \note Tensor contains Hartree elements in the auxiliary basis
    1079             : !> \param qs_env Quickstep environment - entry point of calculation
    1080             : !> \author Stepan Marek
    1081             : !> \date 01.2024
    1082             : ! **************************************************************************************************
    1083          12 :    SUBROUTINE init_hartree(rtbse_env, v_dbcsr)
    1084             :       TYPE(rtbse_env_type), POINTER, INTENT(IN)          :: rtbse_env
    1085             :       TYPE(dbcsr_type)                                   :: v_dbcsr
    1086             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1087             :       TYPE(libint_potential_type)                        :: coulomb_op
    1088             :       TYPE(cp_fm_type)                                   :: V_fm
    1089             :       TYPE(cp_fm_type)                                   :: metric_fm
    1090             :       TYPE(cp_fm_type)                                   :: metric_inv_fm, &
    1091             :                                                             work_fm
    1092             :       TYPE(dbcsr_type), DIMENSION(:), ALLOCATABLE        :: V_dbcsr_a, &
    1093          12 :                                                             metric_dbcsr
    1094             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1095          12 :          POINTER                                         :: nl_2c
    1096             : 
    1097          12 :       bs_env => rtbse_env%bs_env
    1098             : 
    1099             :       ! Allocate for bare Hartree term
    1100          24 :       ALLOCATE (V_dbcsr_a(1))
    1101          24 :       ALLOCATE (metric_dbcsr(1))
    1102          12 :       CALL dbcsr_create(V_dbcsr_a(1), name="Hartree_dbcsr", template=bs_env%mat_RI_RI%matrix)
    1103          12 :       CALL dbcsr_create(metric_dbcsr(1), name="RI_metric_dbcsr", template=bs_env%mat_RI_RI%matrix)
    1104             : 
    1105             :       ! Calculate full coulomb RI basis elements - V _ (PQ) matrix
    1106          12 :       NULLIFY (nl_2c)
    1107             :       CALL build_2c_neighbor_lists(nl_2c, bs_env%basis_set_RI, bs_env%basis_set_RI, &
    1108             :                                    coulomb_op, "Coulomb_neighbor_2c_list", rtbse_env%qs_env, &
    1109          12 :                                    sym_ij=.FALSE., molecular=.TRUE.)
    1110             :       CALL build_2c_integrals(V_dbcsr_a, bs_env%eps_filter, rtbse_env%qs_env, nl_2c, &
    1111             :                               bs_env%basis_set_RI, bs_env%basis_set_RI, coulomb_op, &
    1112          12 :                               do_kpoints=.FALSE., regularization_RI=bs_env%regularization_RI)
    1113             :       ! Calculate the RI metric elements
    1114             :       ! nl_2c is automatically rewritten (even reallocated) in this routine
    1115             :       CALL build_2c_neighbor_lists(nl_2c, bs_env%basis_set_RI, bs_env%basis_set_RI, &
    1116             :                                    bs_env%ri_metric, "Metric_neighbor_2c_list", rtbse_env%qs_env, &
    1117          12 :                                    sym_ij=.FALSE., molecular=.TRUE.)
    1118             :       CALL build_2c_integrals(metric_dbcsr, bs_env%eps_filter, rtbse_env%qs_env, nl_2c, &
    1119             :                               bs_env%basis_set_RI, bs_env%basis_set_RI, bs_env%ri_metric, &
    1120          12 :                               do_kpoints=.FALSE., regularization_RI=bs_env%regularization_RI)
    1121             :       ! nl_2c no longer needed
    1122          12 :       CALL release_neighbor_list_sets(nl_2c)
    1123          12 :       CALL cp_fm_create(metric_fm, bs_env%fm_RI_RI%matrix_struct)
    1124          12 :       CALL cp_fm_set_all(metric_fm, 0.0_dp)
    1125          12 :       CALL cp_fm_create(metric_inv_fm, bs_env%fm_RI_RI%matrix_struct)
    1126          12 :       CALL cp_fm_set_all(metric_inv_fm, 0.0_dp)
    1127          12 :       CALL cp_fm_create(work_fm, bs_env%fm_RI_RI%matrix_struct)
    1128          12 :       CALL cp_fm_set_all(work_fm, 0.0_dp)
    1129          12 :       CALL copy_dbcsr_to_fm(metric_dbcsr(1), metric_fm)
    1130          12 :       CALL cp_fm_invert(metric_fm, metric_inv_fm)
    1131          12 :       CALL cp_fm_create(V_fm, bs_env%fm_RI_RI%matrix_struct)
    1132          12 :       CALL cp_fm_set_all(V_fm, 0.0_dp)
    1133             :       ! Multiply by the inverse from each side (M^-1 is symmetric)
    1134             :       CALL cp_dbcsr_sm_fm_multiply(V_dbcsr_a(1), metric_inv_fm, &
    1135          12 :                                    work_fm, bs_env%n_RI)
    1136             :       CALL parallel_gemm("N", "N", bs_env%n_RI, bs_env%n_RI, bs_env%n_RI, &
    1137          12 :                          1.0_dp, metric_inv_fm, work_fm, 0.0_dp, V_fm)
    1138             :       ! Now, create the tensor from the matrix
    1139             :       ! First, convert full matrix to dbcsr
    1140          12 :       CALL dbcsr_clear(V_dbcsr_a(1))
    1141          12 :       CALL copy_fm_to_dbcsr(V_fm, V_dbcsr_a(1))
    1142          12 :       CALL dbcsr_create(v_dbcsr, "Hartree ri", V_dbcsr_a(1))
    1143          12 :       CALL dbcsr_copy(v_dbcsr, V_dbcsr_a(1))
    1144             :       ! Create and copy distinctly, so that unnecessary objects can be destroyed
    1145             :       ! Destroy all unnecessary matrices
    1146          12 :       CALL dbcsr_release(V_dbcsr_a(1))
    1147          12 :       CALL dbcsr_release(metric_dbcsr(1))
    1148          12 :       DEALLOCATE (V_dbcsr_a)
    1149          12 :       DEALLOCATE (metric_dbcsr)
    1150          12 :       CALL cp_fm_release(V_fm)
    1151             :       ! CALL cp_fm_release(metric_fm(1,1))
    1152          12 :       CALL cp_fm_release(metric_fm)
    1153             :       ! DEALLOCATE(metric_fm)
    1154          12 :       CALL cp_fm_release(work_fm)
    1155          12 :       CALL cp_fm_release(metric_inv_fm)
    1156          72 :    END SUBROUTINE
    1157             : ! **************************************************************************************************
    1158             : !> \brief Calculates the Hartree matrix in the atomic orbital basis, given a density matrix, in local arrays
    1159             : !>        Calculates the values for single spin species present in given rho
    1160             : !> \param qs_env Entry point
    1161             : !> \param rtbse_env Entry point of GWBSE - uses rho_dbcsr and some complex_workspace
    1162             : !> \param rho_ao Density matrix in ao basis
    1163             : !> \param v_ao Overwritten by the Hartree matrix in the atomic orbital basis
    1164             : !> \author Stepan Marek
    1165             : !> \date 02.2024
    1166             : ! **************************************************************************************************
    1167        2444 :    SUBROUTINE get_hartree_local(rtbse_env, rho_ao, v_ao)
    1168             :       TYPE(rtbse_env_type)                               :: rtbse_env
    1169             :       TYPE(cp_cfm_type), INTENT(IN)                      :: rho_ao
    1170             :       TYPE(cp_fm_type)                                   :: v_ao
    1171             :       CHARACTER(len=*), PARAMETER                        :: routineN = "get_hartree_local"
    1172             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1173             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1174             :       TYPE(dbcsr_iterator_type)                          :: iterator_matrix
    1175             :       INTEGER                                            :: i, j, k, n, nblocks, ind_1, ind_2, row_offset, col_offset, &
    1176             :                                                             row_size, col_size, j_n_AO, k_n_AO, i_n_RI, n_RI, &
    1177             :                                                             ri_offset, ind_i, handle
    1178        2444 :       REAL(kind=dp), DIMENSION(:), ALLOCATABLE           :: Pvector, Qvector
    1179        2444 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: block_matrix
    1180        2444 :       REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: int_3c
    1181             :       TYPE(dbcsr_type)                                   :: v_dbcsr_ao
    1182             : 
    1183             :       ! No memory optimisation so far - calculate all 3cs an all ranks
    1184             :       ! Importantly - dbcsr blocks are ordered by atoms - i.e. ethene with 6 atoms will have 6x6 block structure
    1185             :       ! Number of basis states on each basis set is known is post_scf_bandstructure env
    1186             : 
    1187        2444 :       CALL timeset(routineN, handle)
    1188             :       ! Get the relevant environments
    1189        2444 :       CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env, para_env=para_env)
    1190        2444 :       n_RI = rtbse_env%n_RI
    1191        2444 :       int_3c => rtbse_env%int_3c_array
    1192             : 
    1193             :       ! Allocate the Q and Pvector on each rank
    1194       61100 :       ALLOCATE (Qvector(rtbse_env%n_RI), source=0.0_dp)
    1195       58656 :       ALLOCATE (Pvector(rtbse_env%n_RI), source=0.0_dp)
    1196             : 
    1197             :       ! First step - analyze the structure of copied dbcsr matrix on all ranks
    1198        2444 :       CALL dbcsr_clear(rtbse_env%rho_dbcsr)
    1199             :       ! Only the real part of the density matrix contributes
    1200        2444 :       CALL cp_cfm_to_fm(msource=rho_ao, mtargetr=rtbse_env%real_workspace(1))
    1201        2444 :       CALL copy_fm_to_dbcsr(rtbse_env%real_workspace(1), rtbse_env%rho_dbcsr)
    1202        2444 :       nblocks = dbcsr_get_num_blocks(rtbse_env%rho_dbcsr)
    1203        2444 :       CALL dbcsr_iterator_start(iterator_matrix, rtbse_env%rho_dbcsr)
    1204        7332 :       DO n = 1, nblocks
    1205             :          CALL dbcsr_iterator_next_block(iterator_matrix, ind_1, ind_2, block_matrix, &
    1206        4888 :                                         row_offset=row_offset, col_offset=col_offset, row_size=row_size, col_size=col_size)
    1207             :          ! Now we have a block corresponding to a single atom pair
    1208        4888 :          j_n_AO = bs_env%sizes_AO(ind_1)
    1209        4888 :          k_n_AO = bs_env%sizes_AO(ind_2)
    1210             :          ! For each atom pair, we need to get contributions from all RI atoms
    1211             :          ! The allocations are as follows
    1212             :          ! Now, get the relevant 3c integrals
    1213        4888 :          ri_offset = 0
    1214       17108 :          DO ind_i = 1, bs_env%n_atom
    1215        9776 :             i_n_RI = bs_env%sizes_RI(ind_i)
    1216             :             !$OMP PARALLEL DO DEFAULT(none) PRIVATE(i,j,k) &
    1217        9776 :             !$OMP SHARED(Qvector, int_3c, block_matrix,i_n_RI,j_n_AO,k_n_AO, ri_offset, ind_1, ind_2, ind_i, bs_env)
    1218             :             DO i = 1, i_n_RI
    1219             :                DO j = 1, j_n_AO
    1220             :                   DO k = 1, k_n_AO
    1221             :                      ! Different OMP threads write to different places in memory - should be safe from data race
    1222             :                      Qvector(ri_offset + i) = Qvector(ri_offset + i) + int_3c(j + bs_env%i_ao_start_from_atom(ind_1) - 1, &
    1223             :                                                                               k + bs_env%i_ao_start_from_atom(ind_2) - 1, &
    1224             :                                                                               i + bs_env%i_RI_start_from_atom(ind_i) - 1) &
    1225             :                                               *block_matrix(j, k)
    1226             :                   END DO
    1227             :                END DO
    1228             :             END DO
    1229             :             !$OMP END PARALLEL DO
    1230       14664 :             ri_offset = ri_offset + i_n_RI
    1231             :          END DO
    1232             :       END DO
    1233        2444 :       CALL dbcsr_iterator_stop(iterator_matrix)
    1234             :       ! Now, each rank has contributions from D_jk within its scope
    1235             :       ! Need to sum over different ranks to get the total vector on all ranks
    1236        2444 :       CALL para_env%sum(Qvector)
    1237             :       ! Once this is done, Pvector is current on all ranks
    1238             :       ! Continue with V_PQ summation
    1239        2444 :       nblocks = dbcsr_get_num_blocks(rtbse_env%v_dbcsr)
    1240        2444 :       CALL dbcsr_iterator_start(iterator_matrix, rtbse_env%v_dbcsr)
    1241        7332 :       DO n = 1, nblocks
    1242             :          ! TODO : Try OMP parallelisation over different blocks - expect many more available speedup for large systems
    1243             :          CALL dbcsr_iterator_next_block(iterator_matrix, ind_1, ind_2, block_matrix, &
    1244        4888 :                                         row_offset=row_offset, col_offset=col_offset, row_size=row_size, col_size=col_size)
    1245             :          ! TODO : Better names for RI
    1246        4888 :          j_n_AO = bs_env%sizes_RI(ind_1)
    1247        4888 :          k_n_AO = bs_env%sizes_RI(ind_2)
    1248             :          ! The allocations are as follows
    1249             :          !$OMP PARALLEL DO DEFAULT(none) PRIVATE(j,k) &
    1250        7332 :          !$OMP SHARED(block_matrix, Pvector, Qvector,j_n_AO,k_n_AO,row_offset,col_offset)
    1251             :          DO j = 1, j_n_AO
    1252             :             DO k = 1, k_n_AO
    1253             :                Pvector(j + row_offset - 1) = Pvector(j + row_offset - 1) + block_matrix(j, k)*Qvector(k + col_offset - 1)
    1254             :             END DO
    1255             :          END DO
    1256             :          !$OMP END PARALLEL DO
    1257             :       END DO
    1258        2444 :       CALL dbcsr_iterator_stop(iterator_matrix)
    1259             :       ! Again, make sure that the P vector is present on all ranks
    1260        2444 :       CALL para_env%sum(Pvector)
    1261             :       ! Now, for the final trick, iterate over local blocks of v_dbcsr_ao to get the Hartree as dbcsr, then convert to fm
    1262        2444 :       CALL dbcsr_create(v_dbcsr_ao, "Hartree ao", rtbse_env%rho_dbcsr)
    1263        2444 :       CALL copy_fm_to_dbcsr(v_ao, v_dbcsr_ao)
    1264        2444 :       nblocks = dbcsr_get_num_blocks(v_dbcsr_ao)
    1265        2444 :       CALL dbcsr_iterator_start(iterator_matrix, v_dbcsr_ao)
    1266        7332 :       DO n = 1, nblocks
    1267             :          CALL dbcsr_iterator_next_block(iterator_matrix, ind_1, ind_2, block_matrix, &
    1268        4888 :                                         row_offset=row_offset, col_offset=col_offset)
    1269        4888 :          j_n_AO = bs_env%sizes_AO(ind_1)
    1270        4888 :          k_n_AO = bs_env%sizes_AO(ind_2)
    1271        4888 :          ri_offset = 0
    1272       14664 :          block_matrix(:, :) = 0.0_dp
    1273       21996 :          DO ind_i = 1, bs_env%n_atom
    1274        9776 :             i_n_RI = bs_env%sizes_RI(ind_i)
    1275             :             ! TODO : In principle, can parallelize over both directions here
    1276             :             !$OMP PARALLEL DO DEFAULT(none) PRIVATE(i,j,k) &
    1277        9776 :             !$OMP SHARED(block_matrix, Pvector, int_3c, i_n_RI, j_n_AO, k_n_AO, ri_offset, ind_1, ind_2, ind_i, bs_env)
    1278             :             DO j = 1, j_n_AO
    1279             :                DO k = 1, k_n_AO
    1280             :                   ! Add all the summed up values
    1281             :                   DO i = 1, i_n_RI
    1282             :                      block_matrix(j, k) = block_matrix(j, k) + int_3c(j + bs_env%i_ao_start_from_atom(ind_1) - 1, &
    1283             :                                                                       k + bs_env%i_ao_start_from_atom(ind_2) - 1, &
    1284             :                                                                   i + bs_env%i_RI_start_from_atom(ind_i) - 1)*Pvector(i + ri_offset)
    1285             :                   END DO
    1286             :                END DO
    1287             :             END DO
    1288             :             !$OMP END PARALLEL DO
    1289       14664 :             ri_offset = ri_offset + i_n_RI
    1290             :          END DO
    1291             :       END DO
    1292        2444 :       CALL dbcsr_iterator_stop(iterator_matrix)
    1293             :       ! Since P vector was present on all the ranks, v_dbcsr_ao has the complete Hartree result
    1294        2444 :       CALL copy_dbcsr_to_fm(v_dbcsr_ao, v_ao)
    1295        2444 :       CALL dbcsr_release(v_dbcsr_ao)
    1296        2444 :       DEALLOCATE (Qvector)
    1297        2444 :       DEALLOCATE (Pvector)
    1298             : 
    1299        2444 :       CALL timestop(handle)
    1300       12220 :    END SUBROUTINE get_hartree_local
    1301             : ! **************************************************************************************************
    1302             : !> \brief Calculates the exponential of a matrix in a generalized eigenvalue problem. Specifically,
    1303             : !>        it assumes we have a Hermitian matrix A in the eigenvalue problem AX = BXE, where B is some overlap
    1304             : !>        matrix and E is a diagonal matrix of real eigenvalues. Then, it calculates
    1305             : !>        exp(B^(-1) A) = X exp(E) X^C B
    1306             : !> \param amatrix Matrix to exponentiate
    1307             : !> \param bmatrix Overlap matrix
    1308             : !> \param exponential Exponential exp(B^(-1) A) is stored here after the routine is finished
    1309             : !> \param eig_scale_opt Optionally scale eigenvalues by a complex number before exponentiating them
    1310             : !> \param work_opt Optionally provide workspace (of size at least 4) that is used in the calculation
    1311             : !> \author Stepan Marek
    1312             : !> \date 09.2024
    1313             : ! **************************************************************************************************
    1314           0 :    SUBROUTINE cp_cfm_gexp(amatrix, bmatrix, exponential, eig_scale_opt, work_opt)
    1315             :       ! TODO : Do interface for real matrices
    1316             :       TYPE(cp_cfm_type), INTENT(IN)                      :: amatrix
    1317             :       TYPE(cp_cfm_type), INTENT(IN)                      :: bmatrix
    1318             :       TYPE(cp_cfm_type)                                  :: exponential
    1319             :       COMPLEX(kind=dp), INTENT(IN), OPTIONAL             :: eig_scale_opt
    1320             :       TYPE(cp_cfm_type), DIMENSION(:), POINTER, OPTIONAL :: work_opt
    1321             :       CHARACTER(len=*), PARAMETER                        :: routineN = "cp_cfm_gexp"
    1322             :       COMPLEX(kind=dp)                                   :: eig_scale
    1323           0 :       REAL(kind=dp), DIMENSION(:), ALLOCATABLE           :: eigenvalues
    1324           0 :       COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE        :: expvalues
    1325           0 :       TYPE(cp_cfm_type), DIMENSION(:), POINTER           :: work
    1326             :       LOGICAL                                            :: deallocate_work
    1327             :       INTEGER                                            :: nrow, i, handle
    1328             : 
    1329           0 :       CALL timeset(routineN, handle)
    1330             : 
    1331             :       ! Argument parsing and sanity checks
    1332           0 :       IF (PRESENT(eig_scale_opt)) THEN
    1333           0 :          eig_scale = eig_scale_opt
    1334             :       ELSE
    1335             :          eig_scale = CMPLX(1.0, 0.0, kind=dp)
    1336             :       END IF
    1337             : 
    1338           0 :       NULLIFY (work)
    1339           0 :       IF (PRESENT(work_opt) .AND. SIZE(work_opt) >= 4) THEN
    1340           0 :          deallocate_work = .FALSE.
    1341           0 :          work => work_opt
    1342             :       ELSE
    1343           0 :          deallocate_work = .TRUE.
    1344           0 :          ALLOCATE (work(4))
    1345             :          ! Allocate the work storage on the fly
    1346           0 :          DO i = 1, 4
    1347           0 :             CALL cp_cfm_create(work(i), amatrix%matrix_struct)
    1348             :          END DO
    1349             :       END IF
    1350             : 
    1351           0 :       nrow = amatrix%matrix_struct%nrow_global
    1352             : 
    1353           0 :       ALLOCATE (eigenvalues(nrow))
    1354           0 :       ALLOCATE (expvalues(nrow))
    1355             : 
    1356             :       ! Do not change the amatrix and bmatrix - need to copy them first
    1357           0 :       CALL cp_cfm_to_cfm(amatrix, work(1))
    1358           0 :       CALL cp_cfm_to_cfm(bmatrix, work(2))
    1359             : 
    1360             :       ! Solve the generalized eigenvalue equation
    1361           0 :       CALL cp_cfm_geeig(work(1), work(2), work(3), eigenvalues, work(4))
    1362             : 
    1363             :       ! Scale and exponentiate the eigenvalues
    1364           0 :       expvalues(:) = EXP(eigenvalues(:)*eig_scale)
    1365             : 
    1366             :       ! Copy eigenvectors to column scale them
    1367           0 :       CALL cp_cfm_to_cfm(work(3), work(1))
    1368             :       ! X * exp(E)
    1369           0 :       CALL cp_cfm_column_scale(work(1), expvalues)
    1370             : 
    1371             :       ! Carry out the remaining operations
    1372             :       ! X * exp(E) * X^C
    1373             :       CALL parallel_gemm("N", "C", nrow, nrow, nrow, &
    1374             :                          CMPLX(1.0, 0.0, kind=dp), work(1), work(3), &
    1375           0 :                          CMPLX(0.0, 0.0, kind=dp), work(2))
    1376             :       ! X * exp(E) * X^C * B
    1377             :       CALL parallel_gemm("N", "N", nrow, nrow, nrow, &
    1378             :                          CMPLX(1.0, 0.0, kind=dp), work(2), bmatrix, &
    1379           0 :                          CMPLX(0.0, 0.0, kind=dp), exponential)
    1380             : 
    1381             :       ! Deallocate work storage if necessary
    1382           0 :       IF (deallocate_work) THEN
    1383           0 :          DO i = 1, 4
    1384           0 :             CALL cp_cfm_release(work(i))
    1385             :          END DO
    1386           0 :          DEALLOCATE (work)
    1387             :       END IF
    1388             : 
    1389           0 :       DEALLOCATE (eigenvalues)
    1390           0 :       DEALLOCATE (expvalues)
    1391             : 
    1392           0 :       CALL timestop(handle)
    1393           0 :    END SUBROUTINE cp_cfm_gexp
    1394             : END MODULE rt_bse

Generated by: LCOV version 1.15