LCOV - code coverage report
Current view: top level - src/emd - rt_propagation_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:e4f5572) Lines: 358 361 99.2 %
Date: 2023-09-29 06:49:12 Functions: 9 9 100.0 %

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

Generated by: LCOV version 1.15