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

Generated by: LCOV version 1.15