LCOV - code coverage report
Current view: top level - src/emd - rt_propagation_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 361 365 98.9 %
Date: 2025-02-18 08:24:35 Functions: 9 9 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines for propagating the orbitals
      10             : !> \author Florian Schiffmann (02.09)
      11             : ! **************************************************************************************************
      12             : MODULE rt_propagation_methods
      13             :    USE bibliography,                    ONLY: Kolafa2004,&
      14             :                                               Kuhne2007,&
      15             :                                               cite_reference
      16             :    USE cell_types,                      ONLY: cell_type
      17             :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_triangular_multiply
      18             :    USE cp_cfm_cholesky,                 ONLY: cp_cfm_cholesky_decompose
      19             :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      20             :                                               cp_cfm_release,&
      21             :                                               cp_cfm_type
      22             :    USE cp_control_types,                ONLY: dft_control_type,&
      23             :                                               rtp_control_type
      24             :    USE cp_dbcsr_api,                    ONLY: &
      25             :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
      26             :         dbcsr_filter, dbcsr_get_block_p, dbcsr_init_p, dbcsr_iterator_blocks_left, &
      27             :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      28             :         dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, &
      29             :         dbcsr_type, dbcsr_type_antisymmetric
      30             :    USE cp_dbcsr_cholesky,               ONLY: cp_dbcsr_cholesky_decompose,&
      31             :                                               cp_dbcsr_cholesky_invert
      32             :    USE cp_dbcsr_contrib,                ONLY: dbcsr_frobenius_norm
      33             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      34             :                                               dbcsr_allocate_matrix_set,&
      35             :                                               dbcsr_deallocate_matrix_set
      36             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add
      37             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      38             :                                               cp_fm_struct_double,&
      39             :                                               cp_fm_struct_release,&
      40             :                                               cp_fm_struct_type
      41             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      42             :                                               cp_fm_get_info,&
      43             :                                               cp_fm_release,&
      44             :                                               cp_fm_to_fm,&
      45             :                                               cp_fm_type
      46             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      47             :                                               cp_logger_get_default_unit_nr,&
      48             :                                               cp_logger_type,&
      49             :                                               cp_to_string
      50             :    USE efield_utils,                    ONLY: efield_potential_lengh_gauge
      51             :    USE input_constants,                 ONLY: do_arnoldi,&
      52             :                                               do_bch,&
      53             :                                               do_em,&
      54             :                                               do_pade,&
      55             :                                               do_taylor
      56             :    USE iterate_matrix,                  ONLY: matrix_sqrt_Newton_Schulz
      57             :    USE kinds,                           ONLY: dp
      58             :    USE ls_matrix_exp,                   ONLY: cp_complex_dbcsr_gemm_3
      59             :    USE mathlib,                         ONLY: binomial
      60             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      61             :    USE qs_energy_init,                  ONLY: qs_energies_init
      62             :    USE qs_energy_types,                 ONLY: qs_energy_type
      63             :    USE qs_environment_types,            ONLY: get_qs_env,&
      64             :                                               qs_environment_type
      65             :    USE qs_ks_methods,                   ONLY: qs_ks_update_qs_env
      66             :    USE qs_ks_types,                     ONLY: set_ks_env
      67             :    USE rt_make_propagators,             ONLY: propagate_arnoldi,&
      68             :                                               propagate_bch,&
      69             :                                               propagate_exp,&
      70             :                                               propagate_exp_density
      71             :    USE rt_propagation_output,           ONLY: report_density_occupation,&
      72             :                                               rt_convergence,&
      73             :                                               rt_convergence_density
      74             :    USE rt_propagation_types,            ONLY: get_rtp,&
      75             :                                               rt_prop_type
      76             :    USE rt_propagation_utils,            ONLY: calc_S_derivs,&
      77             :                                               calc_update_rho,&
      78             :                                               calc_update_rho_sparse
      79             :    USE rt_propagation_velocity_gauge,   ONLY: update_vector_potential,&
      80             :                                               velocity_gauge_ks_matrix
      81             : #include "../base/base_uses.f90"
      82             : 
      83             :    IMPLICIT NONE
      84             : 
      85             :    PRIVATE
      86             : 
      87             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_methods'
      88             : 
      89             :    PUBLIC :: propagation_step, &
      90             :              s_matrices_create, &
      91             :              calc_sinvH, &
      92             :              put_data_to_history
      93             : 
      94             : CONTAINS
      95             : 
      96             : ! **************************************************************************************************
      97             : !> \brief performs a single propagation step a(t+Dt)=U(t+Dt,t)*a(0)
      98             : !>        and calculates the new exponential
      99             : !> \param qs_env ...
     100             : !> \param rtp ...
     101             : !> \param rtp_control ...
     102             : !> \author Florian Schiffmann (02.09)
     103             : ! **************************************************************************************************
     104             : 
     105        2294 :    SUBROUTINE propagation_step(qs_env, rtp, rtp_control)
     106             : 
     107             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     108             :       TYPE(rt_prop_type), POINTER                        :: rtp
     109             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     110             : 
     111             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'propagation_step'
     112             : 
     113             :       INTEGER                                            :: aspc_order, handle, i, im, re, unit_nr
     114             :       TYPE(cell_type), POINTER                           :: cell
     115        2294 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: delta_mos, mos_new
     116             :       TYPE(cp_logger_type), POINTER                      :: logger
     117        2294 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: delta_P, H_last_iter, ks_mix, ks_mix_im, &
     118        2294 :                                                             matrix_ks, matrix_ks_im, matrix_s, &
     119        2294 :                                                             rho_new
     120             :       TYPE(dft_control_type), POINTER                    :: dft_control
     121             : 
     122        2294 :       CALL timeset(routineN, handle)
     123             : 
     124        2294 :       logger => cp_get_default_logger()
     125        2294 :       IF (logger%para_env%is_source()) THEN
     126        1147 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     127             :       ELSE
     128             :          unit_nr = -1
     129             :       END IF
     130             : 
     131        2294 :       NULLIFY (cell, delta_P, rho_new, delta_mos, mos_new)
     132        2294 :       NULLIFY (ks_mix, ks_mix_im)
     133             :       ! get everything needed and set some values
     134        2294 :       CALL get_qs_env(qs_env, cell=cell, matrix_s=matrix_s, dft_control=dft_control)
     135             : 
     136        2294 :       IF (rtp%iter == 1) THEN
     137         614 :          CALL qs_energies_init(qs_env, .FALSE.)
     138             :          !the above recalculates matrix_s, but matrix not changed if ions are fixed
     139         614 :          IF (rtp_control%fixed_ions) CALL set_ks_env(qs_env%ks_env, s_mstruct_changed=.FALSE.)
     140             : 
     141             :          ! add additional terms to matrix_h and matrix_h_im in the case of applied electric field,
     142             :          ! either in the lengh or velocity gauge.
     143             :          ! should be called  after qs_energies_init and before qs_ks_update_qs_env
     144         614 :          IF (dft_control%apply_efield_field) THEN
     145         216 :             IF (ANY(cell%perd(1:3) /= 0)) &
     146           0 :                CPABORT("Length gauge (efield) and periodicity are not compatible")
     147          54 :             CALL efield_potential_lengh_gauge(qs_env)
     148         560 :          ELSE IF (rtp_control%velocity_gauge) THEN
     149          32 :             IF (dft_control%apply_vector_potential) &
     150          32 :                CALL update_vector_potential(qs_env, dft_control)
     151          32 :             CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.)
     152             :          END IF
     153             : 
     154         614 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     155         614 :          IF (.NOT. rtp_control%fixed_ions) THEN
     156         270 :             CALL s_matrices_create(matrix_s, rtp)
     157             :          END IF
     158         614 :          rtp%delta_iter = 100.0_dp
     159         614 :          rtp%mixing_factor = 1.0_dp
     160         614 :          rtp%mixing = .FALSE.
     161         614 :          aspc_order = rtp_control%aspc_order
     162         614 :          CALL aspc_extrapolate(rtp, matrix_s, aspc_order)
     163         614 :          IF (rtp%linear_scaling) THEN
     164         182 :             CALL calc_update_rho_sparse(qs_env)
     165             :          ELSE
     166         432 :             CALL calc_update_rho(qs_env)
     167             :          END IF
     168         614 :          CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
     169             :       END IF
     170        2294 :       IF (.NOT. rtp_control%fixed_ions) THEN
     171        1150 :          CALL calc_S_derivs(qs_env)
     172             :       END IF
     173        2294 :       rtp%converged = .FALSE.
     174             : 
     175        2294 :       IF (rtp%linear_scaling) THEN
     176             :          ! keep temporary copy of the starting density matrix to check for convergence
     177         816 :          CALL get_rtp(rtp=rtp, rho_new=rho_new)
     178         816 :          NULLIFY (delta_P)
     179         816 :          CALL dbcsr_allocate_matrix_set(delta_P, SIZE(rho_new))
     180        2940 :          DO i = 1, SIZE(rho_new)
     181        2124 :             CALL dbcsr_init_p(delta_P(i)%matrix)
     182        2124 :             CALL dbcsr_create(delta_P(i)%matrix, template=rho_new(i)%matrix)
     183        2940 :             CALL dbcsr_copy(delta_P(i)%matrix, rho_new(i)%matrix)
     184             :          END DO
     185             :       ELSE
     186             :          ! keep temporary copy of the starting mos to check for convergence
     187        1478 :          CALL get_rtp(rtp=rtp, mos_new=mos_new)
     188        8214 :          ALLOCATE (delta_mos(SIZE(mos_new)))
     189        5258 :          DO i = 1, SIZE(mos_new)
     190             :             CALL cp_fm_create(delta_mos(i), &
     191             :                               matrix_struct=mos_new(i)%matrix_struct, &
     192        3780 :                               name="delta_mos"//TRIM(ADJUSTL(cp_to_string(i))))
     193        5258 :             CALL cp_fm_to_fm(mos_new(i), delta_mos(i))
     194             :          END DO
     195             :       END IF
     196             : 
     197             :       CALL get_qs_env(qs_env, &
     198             :                       matrix_ks=matrix_ks, &
     199        2294 :                       matrix_ks_im=matrix_ks_im)
     200             : 
     201        2294 :       CALL get_rtp(rtp=rtp, H_last_iter=H_last_iter)
     202        2294 :       IF (rtp%mixing) THEN
     203          96 :          IF (unit_nr > 0) THEN
     204          48 :             WRITE (unit_nr, '(t3,a,2f16.8)') "Mixing the Hamiltonians to improve robustness, mixing factor: ", rtp%mixing_factor
     205             :          END IF
     206          96 :          CALL dbcsr_allocate_matrix_set(ks_mix, SIZE(matrix_ks))
     207          96 :          CALL dbcsr_allocate_matrix_set(ks_mix_im, SIZE(matrix_ks))
     208         192 :          DO i = 1, SIZE(matrix_ks)
     209          96 :             CALL dbcsr_init_p(ks_mix(i)%matrix)
     210          96 :             CALL dbcsr_create(ks_mix(i)%matrix, template=matrix_ks(1)%matrix)
     211          96 :             CALL dbcsr_init_p(ks_mix_im(i)%matrix)
     212         192 :             CALL dbcsr_create(ks_mix_im(i)%matrix, template=matrix_ks(1)%matrix, matrix_type=dbcsr_type_antisymmetric)
     213             :          END DO
     214         192 :          DO i = 1, SIZE(matrix_ks)
     215          96 :             re = 2*i - 1
     216          96 :             im = 2*i
     217          96 :             CALL dbcsr_add(ks_mix(i)%matrix, matrix_ks(i)%matrix, 0.0_dp, rtp%mixing_factor)
     218          96 :             CALL dbcsr_add(ks_mix(i)%matrix, H_last_iter(re)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor)
     219         192 :             IF (rtp%propagate_complex_ks) THEN
     220           0 :                CALL dbcsr_add(ks_mix_im(i)%matrix, matrix_ks_im(i)%matrix, 0.0_dp, rtp%mixing_factor)
     221           0 :                CALL dbcsr_add(ks_mix_im(i)%matrix, H_last_iter(im)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor)
     222             :             END IF
     223             :          END DO
     224          96 :          CALL calc_SinvH(rtp, ks_mix, ks_mix_im, rtp_control)
     225         192 :          DO i = 1, SIZE(matrix_ks)
     226          96 :             re = 2*i - 1
     227          96 :             im = 2*i
     228          96 :             CALL dbcsr_copy(H_last_iter(re)%matrix, ks_mix(i)%matrix)
     229         192 :             IF (rtp%propagate_complex_ks) THEN
     230           0 :                CALL dbcsr_copy(H_last_iter(im)%matrix, ks_mix_im(i)%matrix)
     231             :             END IF
     232             :          END DO
     233          96 :          CALL dbcsr_deallocate_matrix_set(ks_mix)
     234          96 :          CALL dbcsr_deallocate_matrix_set(ks_mix_im)
     235             :       ELSE
     236        2198 :          CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
     237        5054 :          DO i = 1, SIZE(matrix_ks)
     238        2856 :             re = 2*i - 1
     239        2856 :             im = 2*i
     240        2856 :             CALL dbcsr_copy(H_last_iter(re)%matrix, matrix_ks(i)%matrix)
     241        5054 :             IF (rtp%propagate_complex_ks) THEN
     242         422 :                CALL dbcsr_copy(H_last_iter(im)%matrix, matrix_ks_im(i)%matrix)
     243             :             END IF
     244             :          END DO
     245             :       END IF
     246             : 
     247        2294 :       CALL compute_propagator_matrix(rtp, rtp_control%propagator)
     248             : 
     249        3974 :       SELECT CASE (rtp_control%mat_exp)
     250             :       CASE (do_pade, do_taylor)
     251        1680 :          IF (rtp%linear_scaling) THEN
     252         626 :             CALL propagate_exp_density(rtp, rtp_control)
     253         626 :             CALL calc_update_rho_sparse(qs_env)
     254             :          ELSE
     255        1054 :             CALL propagate_exp(rtp, rtp_control)
     256        1054 :             CALL calc_update_rho(qs_env)
     257             :          END IF
     258             :       CASE (do_arnoldi)
     259         424 :          CALL propagate_arnoldi(rtp, rtp_control)
     260         424 :          CALL calc_update_rho(qs_env)
     261             :       CASE (do_bch)
     262         190 :          CALL propagate_bch(rtp, rtp_control)
     263        2484 :          CALL calc_update_rho_sparse(qs_env)
     264             :       END SELECT
     265        2294 :       CALL step_finalize(qs_env, rtp_control, delta_mos, delta_P)
     266        2294 :       IF (rtp%linear_scaling) THEN
     267         816 :          CALL dbcsr_deallocate_matrix_set(delta_P)
     268             :       ELSE
     269        1478 :          CALL cp_fm_release(delta_mos)
     270             :       END IF
     271             : 
     272        2294 :       CALL timestop(handle)
     273             : 
     274        2294 :    END SUBROUTINE propagation_step
     275             : 
     276             : ! **************************************************************************************************
     277             : !> \brief Performs all the stuff to finish the step:
     278             : !>        convergence checks
     279             : !>        copying stuff into right place for the next step
     280             : !>        updating the history for extrapolation
     281             : !> \param qs_env ...
     282             : !> \param rtp_control ...
     283             : !> \param delta_mos ...
     284             : !> \param delta_P ...
     285             : !> \author Florian Schiffmann (02.09)
     286             : ! **************************************************************************************************
     287             : 
     288        2294 :    SUBROUTINE step_finalize(qs_env, rtp_control, delta_mos, delta_P)
     289             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     290             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     291             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: delta_mos
     292             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: delta_P
     293             : 
     294             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'step_finalize'
     295             : 
     296             :       INTEGER                                            :: handle, i, ihist
     297        2294 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new, mos_old
     298        2294 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: exp_H_new, exp_H_old, matrix_ks, &
     299        2294 :                                                             matrix_ks_im, rho_new, rho_old, s_mat
     300             :       TYPE(qs_energy_type), POINTER                      :: energy
     301             :       TYPE(rt_prop_type), POINTER                        :: rtp
     302             : 
     303        2294 :       CALL timeset(routineN, handle)
     304             : 
     305             :       CALL get_qs_env(qs_env=qs_env, rtp=rtp, matrix_s=s_mat, &
     306        2294 :                       matrix_ks=matrix_ks, matrix_ks_im=matrix_ks_im, energy=energy)
     307        2294 :       CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, exp_H_new=exp_H_new)
     308             : 
     309        2294 :       IF (rtp_control%sc_check_start .LT. rtp%iter) THEN
     310        2294 :          rtp%delta_iter_old = rtp%delta_iter
     311        2294 :          IF (rtp%linear_scaling) THEN
     312         816 :             CALL rt_convergence_density(rtp, delta_P, rtp%delta_iter)
     313             :          ELSE
     314        1478 :             CALL rt_convergence(rtp, s_mat(1)%matrix, delta_mos, rtp%delta_iter)
     315             :          END IF
     316        2294 :          rtp%converged = (rtp%delta_iter .LT. rtp_control%eps_ener)
     317             :          !Apply mixing if scf loop is not converging
     318             : 
     319             :          !It would be better to redo the current step with mixixng,
     320             :          !but currently the decision is made to use mixing from the next step on
     321        2294 :          IF (rtp_control%sc_check_start .LT. rtp%iter + 1) THEN
     322        2294 :             IF (rtp%delta_iter/rtp%delta_iter_old > 0.9) THEN
     323           6 :                rtp%mixing_factor = MAX(rtp%mixing_factor/2.0_dp, 0.125_dp)
     324           6 :                rtp%mixing = .TRUE.
     325             :             END IF
     326             :          END IF
     327             :       END IF
     328             : 
     329        2294 :       IF (rtp%converged) THEN
     330         614 :          IF (rtp%linear_scaling) THEN
     331         182 :             CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
     332             :             CALL purify_mcweeny_complex_nonorth(rho_new, s_mat, rtp%filter_eps, rtp%filter_eps_small, &
     333         182 :                                                 rtp_control%mcweeny_max_iter, rtp_control%mcweeny_eps)
     334         182 :             IF (rtp_control%mcweeny_max_iter > 0) CALL calc_update_rho_sparse(qs_env)
     335         182 :             CALL report_density_occupation(rtp%filter_eps, rho_new)
     336         686 :             DO i = 1, SIZE(rho_new)
     337         686 :                CALL dbcsr_copy(rho_old(i)%matrix, rho_new(i)%matrix)
     338             :             END DO
     339             :          ELSE
     340         432 :             CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
     341        1500 :             DO i = 1, SIZE(mos_new)
     342        1500 :                CALL cp_fm_to_fm(mos_new(i), mos_old(i))
     343             :             END DO
     344             :          END IF
     345         614 :          IF (rtp_control%propagator == do_em) CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
     346        2186 :          DO i = 1, SIZE(exp_H_new)
     347        2186 :             CALL dbcsr_copy(exp_H_old(i)%matrix, exp_H_new(i)%matrix)
     348             :          END DO
     349         614 :          ihist = MOD(rtp%istep, rtp_control%aspc_order) + 1
     350         614 :          IF (rtp_control%fixed_ions) THEN
     351         344 :             CALL put_data_to_history(rtp, rho=rho_new, mos=mos_new, ihist=ihist)
     352             :          ELSE
     353         270 :             CALL put_data_to_history(rtp, rho=rho_new, mos=mos_new, s_mat=s_mat, ihist=ihist)
     354             :          END IF
     355             :       END IF
     356             : 
     357        2294 :       rtp%energy_new = energy%total
     358             : 
     359        2294 :       CALL timestop(handle)
     360             : 
     361        2294 :    END SUBROUTINE step_finalize
     362             : 
     363             : ! **************************************************************************************************
     364             : !> \brief computes the propagator matrix for EM/ETRS, RTP/EMD
     365             : !> \param rtp ...
     366             : !> \param propagator ...
     367             : !> \author Florian Schiffmann (02.09)
     368             : ! **************************************************************************************************
     369             : 
     370        4588 :    SUBROUTINE compute_propagator_matrix(rtp, propagator)
     371             :       TYPE(rt_prop_type), POINTER                        :: rtp
     372             :       INTEGER                                            :: propagator
     373             : 
     374             :       CHARACTER(len=*), PARAMETER :: routineN = 'compute_propagator_matrix'
     375             : 
     376             :       INTEGER                                            :: handle, i
     377             :       REAL(Kind=dp)                                      :: dt, prefac
     378        2294 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: exp_H_new, exp_H_old, propagator_matrix
     379             : 
     380        2294 :       CALL timeset(routineN, handle)
     381             :       CALL get_rtp(rtp=rtp, exp_H_new=exp_H_new, exp_H_old=exp_H_old, &
     382        2294 :                    propagator_matrix=propagator_matrix, dt=dt)
     383             : 
     384        2294 :       prefac = -0.5_dp*dt
     385             : 
     386        8198 :       DO i = 1, SIZE(exp_H_new)
     387        5904 :          CALL dbcsr_add(propagator_matrix(i)%matrix, exp_H_new(i)%matrix, 0.0_dp, prefac)
     388        5904 :          IF (propagator == do_em) &
     389        2358 :             CALL dbcsr_add(propagator_matrix(i)%matrix, exp_H_old(i)%matrix, 1.0_dp, prefac)
     390             :       END DO
     391             : 
     392        2294 :       CALL timestop(handle)
     393             : 
     394        2294 :    END SUBROUTINE compute_propagator_matrix
     395             : 
     396             : ! **************************************************************************************************
     397             : !> \brief computes S_inv*H, if needed Sinv*B and S_inv*H_imag and store these quantities to the
     398             : !> \brief exp_H for the real and imag part (for RTP and EMD)
     399             : !> \param rtp ...
     400             : !> \param matrix_ks ...
     401             : !> \param matrix_ks_im ...
     402             : !> \param rtp_control ...
     403             : !> \author Florian Schiffmann (02.09)
     404             : ! **************************************************************************************************
     405             : 
     406        2504 :    SUBROUTINE calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
     407             :       TYPE(rt_prop_type), POINTER                        :: rtp
     408             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_ks_im
     409             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     410             : 
     411             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'calc_SinvH'
     412             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     413             : 
     414             :       INTEGER                                            :: handle, im, ispin, re
     415        2504 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: exp_H, SinvB, SinvH, SinvH_imag
     416             :       TYPE(dbcsr_type)                                   :: matrix_ks_nosym
     417             :       TYPE(dbcsr_type), POINTER                          :: B_mat, S_inv
     418             : 
     419        2504 :       CALL timeset(routineN, handle)
     420        2504 :       CALL get_rtp(rtp=rtp, S_inv=S_inv, exp_H_new=exp_H)
     421        5734 :       DO ispin = 1, SIZE(matrix_ks)
     422        3230 :          re = ispin*2 - 1
     423        3230 :          im = ispin*2
     424        3230 :          CALL dbcsr_set(exp_H(re)%matrix, zero)
     425        5734 :          CALL dbcsr_set(exp_H(im)%matrix, zero)
     426             :       END DO
     427        2504 :       CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks(1)%matrix, matrix_type="N")
     428             : 
     429             :       ! Real part of S_inv x H -> imag part of exp_H
     430        5734 :       DO ispin = 1, SIZE(matrix_ks)
     431        3230 :          re = ispin*2 - 1
     432        3230 :          im = ispin*2
     433        3230 :          CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, matrix_ks_nosym)
     434             :          CALL dbcsr_multiply("N", "N", one, S_inv, matrix_ks_nosym, zero, exp_H(im)%matrix, &
     435        3230 :                              filter_eps=rtp%filter_eps)
     436        5734 :          IF (.NOT. rtp_control%fixed_ions) THEN
     437        1472 :             CALL get_rtp(rtp=rtp, SinvH=SinvH)
     438        1472 :             CALL dbcsr_copy(SinvH(ispin)%matrix, exp_H(im)%matrix)
     439             :          END IF
     440             :       END DO
     441             : 
     442             :       ! Imag part of S_inv x H -> real part of exp_H
     443        2504 :       IF (rtp%propagate_complex_ks) THEN
     444         910 :          DO ispin = 1, SIZE(matrix_ks)
     445         466 :             re = ispin*2 - 1
     446         466 :             im = ispin*2
     447         466 :             CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
     448         466 :             CALL dbcsr_desymmetrize(matrix_ks_im(ispin)%matrix, matrix_ks_nosym)
     449             :             ! - SinvH_imag is added to exp_H(re)%matrix
     450             :             CALL dbcsr_multiply("N", "N", -one, S_inv, matrix_ks_nosym, zero, exp_H(re)%matrix, &
     451         466 :                                 filter_eps=rtp%filter_eps)
     452         910 :             IF (.NOT. rtp_control%fixed_ions) THEN
     453         266 :                CALL get_rtp(rtp=rtp, SinvH_imag=SinvH_imag)
     454             :                ! -SinvH_imag is saved
     455         266 :                CALL dbcsr_copy(SinvH_imag(ispin)%matrix, exp_H(re)%matrix)
     456             :             END IF
     457             :          END DO
     458             :       END IF
     459             :       ! EMD case: the real part of exp_H should be updated with B
     460        2504 :       IF (.NOT. rtp_control%fixed_ions) THEN
     461        1230 :          CALL get_rtp(rtp=rtp, B_mat=B_mat, SinvB=SinvB)
     462        1230 :          CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
     463        1230 :          CALL dbcsr_multiply("N", "N", one, S_inv, B_mat, zero, matrix_ks_nosym, filter_eps=rtp%filter_eps)
     464        2702 :          DO ispin = 1, SIZE(matrix_ks)
     465        1472 :             re = ispin*2 - 1
     466        1472 :             im = ispin*2
     467             :             ! + SinvB is added to exp_H(re)%matrix
     468        1472 :             CALL dbcsr_add(exp_H(re)%matrix, matrix_ks_nosym, 1.0_dp, 1.0_dp)
     469             :             ! + SinvB is saved
     470        2702 :             CALL dbcsr_copy(SinvB(ispin)%matrix, matrix_ks_nosym)
     471             :          END DO
     472             :       END IF
     473             :       ! Otherwise no real part for exp_H
     474             : 
     475        2504 :       CALL dbcsr_release(matrix_ks_nosym)
     476        2504 :       CALL timestop(handle)
     477             : 
     478        2504 :    END SUBROUTINE calc_SinvH
     479             : 
     480             : ! **************************************************************************************************
     481             : !> \brief calculates the needed overlap-like matrices
     482             : !>        depending on the way the exponential is calculated, only S^-1 is needed
     483             : !> \param s_mat ...
     484             : !> \param rtp ...
     485             : !> \author Florian Schiffmann (02.09)
     486             : ! **************************************************************************************************
     487             : 
     488         468 :    SUBROUTINE s_matrices_create(s_mat, rtp)
     489             : 
     490             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: s_mat
     491             :       TYPE(rt_prop_type), POINTER                        :: rtp
     492             : 
     493             :       CHARACTER(len=*), PARAMETER                        :: routineN = 's_matrices_create'
     494             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     495             : 
     496             :       INTEGER                                            :: handle
     497             :       TYPE(dbcsr_type), POINTER                          :: S_half, S_inv, S_minus_half
     498             : 
     499         468 :       CALL timeset(routineN, handle)
     500             : 
     501         468 :       CALL get_rtp(rtp=rtp, S_inv=S_inv)
     502             : 
     503         468 :       IF (rtp%linear_scaling) THEN
     504         136 :          CALL get_rtp(rtp=rtp, S_half=S_half, S_minus_half=S_minus_half)
     505             :          CALL matrix_sqrt_Newton_Schulz(S_half, S_minus_half, s_mat(1)%matrix, rtp%filter_eps, &
     506         136 :                                         rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
     507             :          CALL dbcsr_multiply("N", "N", one, S_minus_half, S_minus_half, zero, S_inv, &
     508         136 :                              filter_eps=rtp%filter_eps)
     509             :       ELSE
     510         332 :          CALL dbcsr_copy(S_inv, s_mat(1)%matrix)
     511             :          CALL cp_dbcsr_cholesky_decompose(S_inv, para_env=rtp%ao_ao_fmstruct%para_env, &
     512         332 :                                           blacs_env=rtp%ao_ao_fmstruct%context)
     513             :          CALL cp_dbcsr_cholesky_invert(S_inv, para_env=rtp%ao_ao_fmstruct%para_env, &
     514         332 :                                        blacs_env=rtp%ao_ao_fmstruct%context, uplo_to_full=.TRUE.)
     515             :       END IF
     516             : 
     517         468 :       CALL timestop(handle)
     518         468 :    END SUBROUTINE s_matrices_create
     519             : 
     520             : ! **************************************************************************************************
     521             : !> \brief Calculates the frobenius norm of a complex matrix represented by two real matrices
     522             : !> \param frob_norm ...
     523             : !> \param mat_re ...
     524             : !> \param mat_im ...
     525             : !> \author Samuel Andermatt (04.14)
     526             : ! **************************************************************************************************
     527             : 
     528         568 :    SUBROUTINE complex_frobenius_norm(frob_norm, mat_re, mat_im)
     529             : 
     530             :       REAL(KIND=dp), INTENT(out)                         :: frob_norm
     531             :       TYPE(dbcsr_type), POINTER                          :: mat_re, mat_im
     532             : 
     533             :       CHARACTER(len=*), PARAMETER :: routineN = 'complex_frobenius_norm'
     534             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     535             : 
     536             :       INTEGER                                            :: col_atom, handle, row_atom
     537             :       LOGICAL                                            :: found
     538         284 :       REAL(dp), DIMENSION(:, :), POINTER                 :: block_values, block_values2
     539             :       TYPE(dbcsr_iterator_type)                          :: iter
     540             :       TYPE(dbcsr_type), POINTER                          :: tmp
     541             : 
     542         284 :       CALL timeset(routineN, handle)
     543             : 
     544             :       NULLIFY (tmp)
     545         284 :       ALLOCATE (tmp)
     546         284 :       CALL dbcsr_create(tmp, template=mat_re)
     547             :       !make sure the tmp has the same sparsity pattern as the real and the complex part combined
     548         284 :       CALL dbcsr_add(tmp, mat_re, zero, one)
     549         284 :       CALL dbcsr_add(tmp, mat_im, zero, one)
     550         284 :       CALL dbcsr_set(tmp, zero)
     551             :       !calculate the hadamard product
     552         284 :       CALL dbcsr_iterator_start(iter, tmp)
     553        1804 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     554        1520 :          CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
     555        1520 :          CALL dbcsr_get_block_p(mat_re, row_atom, col_atom, block_values2, found=found)
     556        1520 :          IF (found) THEN
     557      264788 :             block_values = block_values2*block_values2
     558             :          END IF
     559        1520 :          CALL dbcsr_get_block_p(mat_im, row_atom, col_atom, block_values2, found=found)
     560        1520 :          IF (found) THEN
     561      250308 :             block_values = block_values + block_values2*block_values2
     562             :          END IF
     563      133438 :          block_values = SQRT(block_values)
     564             :       END DO
     565         284 :       CALL dbcsr_iterator_stop(iter)
     566         284 :       frob_norm = dbcsr_frobenius_norm(tmp)
     567             : 
     568         284 :       CALL dbcsr_deallocate_matrix(tmp)
     569             : 
     570         284 :       CALL timestop(handle)
     571             : 
     572         284 :    END SUBROUTINE complex_frobenius_norm
     573             : 
     574             : ! **************************************************************************************************
     575             : !> \brief Does McWeeny for complex matrices in the non-orthogonal basis
     576             : !> \param P ...
     577             : !> \param s_mat ...
     578             : !> \param eps ...
     579             : !> \param eps_small ...
     580             : !> \param max_iter ...
     581             : !> \param threshold ...
     582             : !> \author Samuel Andermatt (04.14)
     583             : ! **************************************************************************************************
     584             : 
     585         182 :    SUBROUTINE purify_mcweeny_complex_nonorth(P, s_mat, eps, eps_small, max_iter, threshold)
     586             : 
     587             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: P, s_mat
     588             :       REAL(KIND=dp), INTENT(in)                          :: eps, eps_small
     589             :       INTEGER, INTENT(in)                                :: max_iter
     590             :       REAL(KIND=dp), INTENT(in)                          :: threshold
     591             : 
     592             :       CHARACTER(len=*), PARAMETER :: routineN = 'purify_mcweeny_complex_nonorth'
     593             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     594             : 
     595             :       INTEGER                                            :: handle, i, im, imax, ispin, re, unit_nr
     596             :       REAL(KIND=dp)                                      :: frob_norm
     597             :       TYPE(cp_logger_type), POINTER                      :: logger
     598         182 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: PS, PSP, tmp
     599             : 
     600         182 :       CALL timeset(routineN, handle)
     601             : 
     602         182 :       logger => cp_get_default_logger()
     603         182 :       IF (logger%para_env%is_source()) THEN
     604          91 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     605             :       ELSE
     606             :          unit_nr = -1
     607             :       END IF
     608             : 
     609         182 :       NULLIFY (tmp, PS, PSP)
     610         182 :       CALL dbcsr_allocate_matrix_set(tmp, SIZE(P))
     611         182 :       CALL dbcsr_allocate_matrix_set(PSP, SIZE(P))
     612         182 :       CALL dbcsr_allocate_matrix_set(PS, SIZE(P))
     613         686 :       DO i = 1, SIZE(P)
     614         504 :          CALL dbcsr_init_p(PS(i)%matrix)
     615         504 :          CALL dbcsr_create(PS(i)%matrix, template=P(1)%matrix)
     616         504 :          CALL dbcsr_init_p(PSP(i)%matrix)
     617         504 :          CALL dbcsr_create(PSP(i)%matrix, template=P(1)%matrix)
     618         504 :          CALL dbcsr_init_p(tmp(i)%matrix)
     619         686 :          CALL dbcsr_create(tmp(i)%matrix, template=P(1)%matrix)
     620             :       END DO
     621         182 :       IF (SIZE(P) == 2) THEN
     622         112 :          CALL dbcsr_scale(P(1)%matrix, one/2)
     623         112 :          CALL dbcsr_scale(P(2)%matrix, one/2)
     624             :       END IF
     625         434 :       DO ispin = 1, SIZE(P)/2
     626         252 :          re = 2*ispin - 1
     627         252 :          im = 2*ispin
     628         252 :          imax = MAX(max_iter, 1) !if max_iter is 0 then only the deviation from idempotency needs to be calculated
     629         516 :          DO i = 1, imax
     630             :             CALL dbcsr_multiply("N", "N", one, P(re)%matrix, s_mat(1)%matrix, &
     631         284 :                                 zero, PS(re)%matrix, filter_eps=eps_small)
     632             :             CALL dbcsr_multiply("N", "N", one, P(im)%matrix, s_mat(1)%matrix, &
     633         284 :                                 zero, PS(im)%matrix, filter_eps=eps_small)
     634             :             CALL cp_complex_dbcsr_gemm_3("N", "N", one, PS(re)%matrix, PS(im)%matrix, &
     635             :                                          P(re)%matrix, P(im)%matrix, zero, PSP(re)%matrix, PSP(im)%matrix, &
     636         284 :                                          filter_eps=eps_small)
     637         284 :             CALL dbcsr_copy(tmp(re)%matrix, PSP(re)%matrix)
     638         284 :             CALL dbcsr_copy(tmp(im)%matrix, PSP(im)%matrix)
     639         284 :             CALL dbcsr_add(tmp(re)%matrix, P(re)%matrix, 1.0_dp, -1.0_dp)
     640         284 :             CALL dbcsr_add(tmp(im)%matrix, P(im)%matrix, 1.0_dp, -1.0_dp)
     641         284 :             CALL complex_frobenius_norm(frob_norm, tmp(re)%matrix, tmp(im)%matrix)
     642         284 :             IF (unit_nr .GT. 0) WRITE (unit_nr, '(t3,a,2f16.8)') "Deviation from idempotency: ", frob_norm
     643         800 :             IF (frob_norm .GT. threshold .AND. max_iter > 0) THEN
     644         264 :                CALL dbcsr_copy(P(re)%matrix, PSP(re)%matrix)
     645         264 :                CALL dbcsr_copy(P(im)%matrix, PSP(im)%matrix)
     646             :                CALL cp_complex_dbcsr_gemm_3("N", "N", -2.0_dp, PS(re)%matrix, PS(im)%matrix, &
     647             :                                             PSP(re)%matrix, PSP(im)%matrix, 3.0_dp, P(re)%matrix, P(im)%matrix, &
     648         264 :                                             filter_eps=eps_small)
     649         264 :                CALL dbcsr_filter(P(re)%matrix, eps)
     650         264 :                CALL dbcsr_filter(P(im)%matrix, eps)
     651             :                !make sure P is exactly hermitian
     652         264 :                CALL dbcsr_transposed(tmp(re)%matrix, P(re)%matrix)
     653         264 :                CALL dbcsr_add(P(re)%matrix, tmp(re)%matrix, one/2, one/2)
     654         264 :                CALL dbcsr_transposed(tmp(im)%matrix, P(im)%matrix)
     655         264 :                CALL dbcsr_add(P(im)%matrix, tmp(im)%matrix, one/2, -one/2)
     656             :             ELSE
     657             :                EXIT
     658             :             END IF
     659             :          END DO
     660             :          !make sure P is hermitian
     661         252 :          CALL dbcsr_transposed(tmp(re)%matrix, P(re)%matrix)
     662         252 :          CALL dbcsr_add(P(re)%matrix, tmp(re)%matrix, one/2, one/2)
     663         252 :          CALL dbcsr_transposed(tmp(im)%matrix, P(im)%matrix)
     664         434 :          CALL dbcsr_add(P(im)%matrix, tmp(im)%matrix, one/2, -one/2)
     665             :       END DO
     666         182 :       IF (SIZE(P) == 2) THEN
     667         112 :          CALL dbcsr_scale(P(1)%matrix, one*2)
     668         112 :          CALL dbcsr_scale(P(2)%matrix, one*2)
     669             :       END IF
     670         182 :       CALL dbcsr_deallocate_matrix_set(tmp)
     671         182 :       CALL dbcsr_deallocate_matrix_set(PS)
     672         182 :       CALL dbcsr_deallocate_matrix_set(PSP)
     673             : 
     674         182 :       CALL timestop(handle)
     675             : 
     676         182 :    END SUBROUTINE purify_mcweeny_complex_nonorth
     677             : 
     678             : ! **************************************************************************************************
     679             : !> \brief ...
     680             : !> \param rtp ...
     681             : !> \param matrix_s ...
     682             : !> \param aspc_order ...
     683             : ! **************************************************************************************************
     684         614 :    SUBROUTINE aspc_extrapolate(rtp, matrix_s, aspc_order)
     685             :       TYPE(rt_prop_type), POINTER                        :: rtp
     686             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     687             :       INTEGER, INTENT(in)                                :: aspc_order
     688             : 
     689             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'aspc_extrapolate'
     690             :       COMPLEX(KIND=dp), PARAMETER                        :: cone = (1.0_dp, 0.0_dp), &
     691             :                                                             czero = (0.0_dp, 0.0_dp)
     692             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     693             : 
     694             :       INTEGER                                            :: handle, i, iaspc, icol_local, ihist, &
     695             :                                                             imat, k, kdbl, n, naspc, ncol_local, &
     696             :                                                             nmat
     697             :       REAL(KIND=dp)                                      :: alpha
     698             :       TYPE(cp_cfm_type)                                  :: cfm_tmp, cfm_tmp1, csc
     699             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct, matrix_struct_new
     700             :       TYPE(cp_fm_type)                                   :: fm_tmp, fm_tmp1, fm_tmp2
     701         614 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     702         614 :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: mo_hist
     703         614 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_new, s_hist
     704         614 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_hist
     705             : 
     706         614 :       NULLIFY (rho_hist)
     707         614 :       CALL timeset(routineN, handle)
     708         614 :       CALL cite_reference(Kolafa2004)
     709         614 :       CALL cite_reference(Kuhne2007)
     710             : 
     711         614 :       IF (rtp%linear_scaling) THEN
     712         182 :          CALL get_rtp(rtp=rtp, rho_new=rho_new)
     713             :       ELSE
     714         432 :          CALL get_rtp(rtp=rtp, mos_new=mos_new)
     715             :       END IF
     716             : 
     717         614 :       naspc = MIN(rtp%istep, aspc_order)
     718         614 :       IF (rtp%linear_scaling) THEN
     719         182 :          nmat = SIZE(rho_new)
     720         182 :          rho_hist => rtp%history%rho_history
     721         686 :          DO imat = 1, nmat
     722        1454 :             DO iaspc = 1, naspc
     723             :                alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
     724         768 :                        binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
     725         768 :                ihist = MOD(rtp%istep - iaspc, aspc_order) + 1
     726        1272 :                IF (iaspc == 1) THEN
     727         504 :                   CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, zero, alpha)
     728             :                ELSE
     729         264 :                   CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, one, alpha)
     730             :                END IF
     731             :             END DO
     732             :          END DO
     733             :       ELSE
     734         432 :          mo_hist => rtp%history%mo_history
     735         432 :          nmat = SIZE(mos_new)
     736        1500 :          DO imat = 1, nmat
     737        3876 :             DO iaspc = 1, naspc
     738             :                alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
     739        2376 :                        binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
     740        2376 :                ihist = MOD(rtp%istep - iaspc, aspc_order) + 1
     741        3444 :                IF (iaspc == 1) THEN
     742        1068 :                   CALL cp_fm_scale_and_add(zero, mos_new(imat), alpha, mo_hist(imat, ihist))
     743             :                ELSE
     744        1308 :                   CALL cp_fm_scale_and_add(one, mos_new(imat), alpha, mo_hist(imat, ihist))
     745             :                END IF
     746             :             END DO
     747             :          END DO
     748             : 
     749         432 :          mo_hist => rtp%history%mo_history
     750         432 :          s_hist => rtp%history%s_history
     751         966 :          DO i = 1, SIZE(mos_new)/2
     752         534 :             NULLIFY (matrix_struct, matrix_struct_new)
     753             : 
     754             :             CALL cp_fm_struct_double(matrix_struct, &
     755             :                                      mos_new(2*i)%matrix_struct, &
     756             :                                      mos_new(2*i)%matrix_struct%context, &
     757         534 :                                      .TRUE., .FALSE.)
     758             : 
     759         534 :             CALL cp_fm_create(fm_tmp, matrix_struct)
     760         534 :             CALL cp_fm_create(fm_tmp1, matrix_struct)
     761         534 :             CALL cp_fm_create(fm_tmp2, mos_new(2*i)%matrix_struct)
     762         534 :             CALL cp_cfm_create(cfm_tmp, mos_new(2*i)%matrix_struct)
     763         534 :             CALL cp_cfm_create(cfm_tmp1, mos_new(2*i)%matrix_struct)
     764             : 
     765         534 :             CALL cp_fm_get_info(fm_tmp, ncol_global=kdbl)
     766             : 
     767             :             CALL cp_fm_get_info(mos_new(2*i), &
     768             :                                 nrow_global=n, &
     769             :                                 ncol_global=k, &
     770         534 :                                 ncol_local=ncol_local)
     771             : 
     772             :             CALL cp_fm_struct_create(matrix_struct_new, &
     773             :                                      template_fmstruct=mos_new(2*i)%matrix_struct, &
     774             :                                      nrow_global=k, &
     775         534 :                                      ncol_global=k)
     776         534 :             CALL cp_cfm_create(csc, matrix_struct_new)
     777             : 
     778         534 :             CALL cp_fm_struct_release(matrix_struct_new)
     779         534 :             CALL cp_fm_struct_release(matrix_struct)
     780             : 
     781             :             ! first the most recent
     782             : 
     783             : ! reorthogonalize vectors
     784        2292 :             DO icol_local = 1, ncol_local
     785       13980 :                fm_tmp%local_data(:, icol_local) = mos_new(2*i - 1)%local_data(:, icol_local)
     786       14514 :                fm_tmp%local_data(:, icol_local + ncol_local) = mos_new(2*i)%local_data(:, icol_local)
     787             :             END DO
     788             : 
     789         534 :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, fm_tmp, fm_tmp1, kdbl)
     790             : 
     791        2292 :             DO icol_local = 1, ncol_local
     792             :                cfm_tmp%local_data(:, icol_local) = CMPLX(fm_tmp1%local_data(:, icol_local), &
     793       13980 :                                                          fm_tmp1%local_data(:, icol_local + ncol_local), dp)
     794             :                cfm_tmp1%local_data(:, icol_local) = CMPLX(mos_new(2*i - 1)%local_data(:, icol_local), &
     795       14514 :                                                           mos_new(2*i)%local_data(:, icol_local), dp)
     796             :             END DO
     797         534 :             CALL parallel_gemm('C', 'N', k, k, n, cone, cfm_tmp1, cfm_tmp, czero, csc)
     798         534 :             CALL cp_cfm_cholesky_decompose(csc)
     799         534 :             CALL cp_cfm_triangular_multiply(csc, cfm_tmp1, n_cols=k, side='R', invert_tr=.TRUE.)
     800        2292 :             DO icol_local = 1, ncol_local
     801       13980 :                mos_new(2*i - 1)%local_data(:, icol_local) = REAL(cfm_tmp1%local_data(:, icol_local), dp)
     802       14514 :                mos_new(2*i)%local_data(:, icol_local) = AIMAG(cfm_tmp1%local_data(:, icol_local))
     803             :             END DO
     804             : 
     805             : ! deallocate work matrices
     806         534 :             CALL cp_cfm_release(csc)
     807         534 :             CALL cp_fm_release(fm_tmp)
     808         534 :             CALL cp_fm_release(fm_tmp1)
     809         534 :             CALL cp_fm_release(fm_tmp2)
     810         534 :             CALL cp_cfm_release(cfm_tmp)
     811        2034 :             CALL cp_cfm_release(cfm_tmp1)
     812             :          END DO
     813             : 
     814             :       END IF
     815             : 
     816         614 :       CALL timestop(handle)
     817             : 
     818         614 :    END SUBROUTINE aspc_extrapolate
     819             : 
     820             : ! **************************************************************************************************
     821             : !> \brief ...
     822             : !> \param rtp ...
     823             : !> \param mos ...
     824             : !> \param rho ...
     825             : !> \param s_mat ...
     826             : !> \param ihist ...
     827             : ! **************************************************************************************************
     828         812 :    SUBROUTINE put_data_to_history(rtp, mos, rho, s_mat, ihist)
     829             :       TYPE(rt_prop_type), POINTER                        :: rtp
     830             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos
     831             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho
     832             :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
     833             :          POINTER                                         :: s_mat
     834             :       INTEGER                                            :: ihist
     835             : 
     836             :       INTEGER                                            :: i
     837             : 
     838         812 :       IF (rtp%linear_scaling) THEN
     839        1032 :          DO i = 1, SIZE(rho)
     840        1032 :             CALL dbcsr_copy(rtp%history%rho_history(i, ihist)%matrix, rho(i)%matrix)
     841             :          END DO
     842             :       ELSE
     843        1884 :          DO i = 1, SIZE(mos)
     844        1884 :             CALL cp_fm_to_fm(mos(i), rtp%history%mo_history(i, ihist))
     845             :          END DO
     846         540 :          IF (PRESENT(s_mat)) THEN
     847         332 :             IF (ASSOCIATED(rtp%history%s_history(ihist)%matrix)) THEN ! the sparsity might be different
     848             :                ! (future struct:check)
     849         124 :                CALL dbcsr_deallocate_matrix(rtp%history%s_history(ihist)%matrix)
     850             :             END IF
     851         332 :             ALLOCATE (rtp%history%s_history(ihist)%matrix)
     852         332 :             CALL dbcsr_copy(rtp%history%s_history(ihist)%matrix, s_mat(1)%matrix)
     853             :          END IF
     854             :       END IF
     855             : 
     856         812 :    END SUBROUTINE put_data_to_history
     857             : 
     858             : END MODULE rt_propagation_methods

Generated by: LCOV version 1.15