LCOV - code coverage report
Current view: top level - src/emd - rt_propagation_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:c5411e0) Lines: 361 365 98.9 %
Date: 2024-05-16 06:48:40 Functions: 9 9 100.0 %

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

Generated by: LCOV version 1.15