LCOV - code coverage report
Current view: top level - src/emd - rt_propagator_init.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 194 199 97.5 %
Date: 2023-03-30 11:55:16 Functions: 5 5 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 that prepare rtp and EMD
      10             : !> \author Florian Schiffmann (02.09)
      11             : ! **************************************************************************************************
      12             : MODULE rt_propagator_init
      13             : 
      14             :    USE arnoldi_api,                     ONLY: arnoldi_extremal
      15             :    USE cp_control_types,                ONLY: dft_control_type,&
      16             :                                               rtp_control_type
      17             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      18             :                                               copy_fm_to_dbcsr,&
      19             :                                               cp_dbcsr_plus_fm_fm_t
      20             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      21             :                                               cp_fm_scale,&
      22             :                                               cp_fm_upper_to_full
      23             :    USE cp_fm_diag,                      ONLY: cp_fm_syevd
      24             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      25             :                                               cp_fm_get_info,&
      26             :                                               cp_fm_release,&
      27             :                                               cp_fm_set_all,&
      28             :                                               cp_fm_to_fm,&
      29             :                                               cp_fm_type
      30             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      31             :                                               cp_logger_get_default_unit_nr,&
      32             :                                               cp_logger_type
      33             :    USE dbcsr_api,                       ONLY: &
      34             :         dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_filter, dbcsr_multiply, &
      35             :         dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type
      36             :    USE input_constants,                 ONLY: do_arnoldi,&
      37             :                                               do_bch,&
      38             :                                               do_cn,&
      39             :                                               do_em,&
      40             :                                               do_etrs,&
      41             :                                               do_pade,&
      42             :                                               do_taylor
      43             :    USE iterate_matrix,                  ONLY: matrix_sqrt_Newton_Schulz
      44             :    USE kinds,                           ONLY: dp
      45             :    USE matrix_exp,                      ONLY: get_nsquare_norder
      46             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      47             :    USE qs_environment_types,            ONLY: get_qs_env,&
      48             :                                               qs_environment_type
      49             :    USE qs_mo_types,                     ONLY: mo_set_type
      50             :    USE rt_make_propagators,             ONLY: compute_exponential,&
      51             :                                               compute_exponential_sparse,&
      52             :                                               propagate_arnoldi
      53             :    USE rt_propagation_methods,          ONLY: calc_SinvH,&
      54             :                                               put_data_to_history,&
      55             :                                               s_matrices_create
      56             :    USE rt_propagation_types,            ONLY: get_rtp,&
      57             :                                               rt_prop_release_mos,&
      58             :                                               rt_prop_type
      59             : #include "../base/base_uses.f90"
      60             : 
      61             :    IMPLICIT NONE
      62             : 
      63             :    PRIVATE
      64             : 
      65             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagator_init'
      66             : 
      67             :    PUBLIC :: init_propagators, &
      68             :              rt_initialize_rho_from_mos
      69             : 
      70             : CONTAINS
      71             : 
      72             : ! **************************************************************************************************
      73             : !> \brief prepares the initial matrices for the propagators
      74             : !> \param qs_env ...
      75             : !> \author Florian Schiffmann (02.09)
      76             : ! **************************************************************************************************
      77             : 
      78         356 :    SUBROUTINE init_propagators(qs_env)
      79             : 
      80             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      81             : 
      82             :       INTEGER                                            :: i, imat, unit_nr
      83             :       REAL(KIND=dp)                                      :: dt, prefac
      84         178 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new, mos_next, mos_old
      85             :       TYPE(cp_logger_type), POINTER                      :: logger
      86         178 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: exp_H_new, exp_H_old, matrix_ks, &
      87         178 :                                                             matrix_ks_im, propagator_matrix, &
      88         178 :                                                             rho_old, s_mat
      89             :       TYPE(dft_control_type), POINTER                    :: dft_control
      90             :       TYPE(rt_prop_type), POINTER                        :: rtp
      91             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
      92             : 
      93             :       CALL get_qs_env(qs_env, &
      94             :                       rtp=rtp, &
      95             :                       dft_control=dft_control, &
      96             :                       matrix_s=s_mat, &
      97             :                       matrix_ks=matrix_ks, &
      98         178 :                       matrix_ks_im=matrix_ks_im)
      99             : 
     100         178 :       rtp_control => dft_control%rtp_control
     101             :       CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, exp_H_new=exp_H_new, &
     102         178 :                    propagator_matrix=propagator_matrix, dt=dt)
     103         178 :       CALL s_matrices_create(s_mat, rtp)
     104         178 :       CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
     105         654 :       DO i = 1, SIZE(exp_H_old)
     106         654 :          CALL dbcsr_copy(exp_H_old(i)%matrix, exp_H_new(i)%matrix)
     107             :       END DO
     108             :       ! use the fact that CN propagator is a first order pade approximation on the EM propagator
     109         178 :       IF (rtp_control%propagator == do_cn) THEN
     110          12 :          rtp%orders(1, :) = 0; rtp%orders(2, :) = 1; rtp_control%mat_exp = do_pade; rtp_control%propagator = do_em
     111         174 :       ELSE IF (rtp_control%mat_exp == do_pade .OR. rtp_control%mat_exp == do_taylor) THEN
     112         142 :          IF (rtp%linear_scaling) THEN
     113          76 :             CALL get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
     114             :          ELSE
     115          66 :             CALL get_maxabs_eigval(rtp, s_mat, matrix_ks, rtp_control)
     116             :          END IF
     117             :       END IF
     118         178 :       IF (rtp_control%mat_exp == do_pade .AND. rtp%linear_scaling) THEN
     119             :          ! get a useful output_unit
     120           0 :          logger => cp_get_default_logger()
     121           0 :          IF (logger%para_env%is_source()) THEN
     122           0 :             unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     123           0 :             WRITE (unit_nr, *) "linear_scaling currently does not support pade exponentials, switching to taylor"
     124             :          END IF
     125           0 :          rtp_control%mat_exp = do_taylor
     126             :       END IF
     127             :       ! We have no clue yet about next H so we use initial H for t and t+dt
     128             :       ! Due to different nature of the propagator the prefactor has to be adopted
     129         350 :       SELECT CASE (rtp_control%propagator)
     130             :       CASE (do_etrs)
     131         172 :          prefac = -0.5_dp*dt
     132             :       CASE (do_em)
     133         178 :          prefac = -1.0_dp*dt
     134             :       END SELECT
     135         654 :       DO imat = 1, SIZE(exp_H_new)
     136         476 :          CALL dbcsr_copy(propagator_matrix(imat)%matrix, exp_H_new(imat)%matrix)
     137         654 :          CALL dbcsr_scale(propagator_matrix(imat)%matrix, prefac)
     138             :       END DO
     139             : 
     140             :       ! For ETRS this bit could be avoided but it drastically simplifies the workflow afterwards.
     141             :       ! If we compute the half propagated mos/exponential already here, we ensure everything is done
     142             :       ! with the correct S matrix and all information as during RTP/EMD are computed.
     143             :       ! Therefore we might accept to compute an unnesscesary exponential but understand the code afterwards
     144         178 :       IF (rtp_control%propagator == do_etrs) THEN
     145         172 :          IF (rtp_control%mat_exp == do_arnoldi) THEN
     146          18 :             rtp%iter = 0
     147          18 :             CALL propagate_arnoldi(rtp, rtp_control)
     148          18 :             CALL get_rtp(rtp=rtp, mos_new=mos_new, mos_next=mos_next)
     149          66 :             DO imat = 1, SIZE(mos_new)
     150          66 :                CALL cp_fm_to_fm(mos_new(imat), mos_next(imat))
     151             :             END DO
     152         154 :          ELSEIF (rtp_control%mat_exp == do_bch) THEN
     153             :          ELSE
     154         140 :             IF (rtp%linear_scaling) THEN
     155          76 :                CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp_control, rtp)
     156             :             ELSE
     157          64 :                CALL compute_exponential(exp_H_new, propagator_matrix, rtp_control, rtp)
     158             :             END IF
     159         524 :             DO imat = 1, SIZE(exp_H_new)
     160         524 :                CALL dbcsr_copy(exp_H_old(imat)%matrix, exp_H_new(imat)%matrix)
     161             :             END DO
     162             :          END IF
     163             :       END IF
     164             : 
     165         178 :       IF (rtp%linear_scaling) THEN
     166          90 :          CALL get_rtp(rtp=rtp, rho_old=rho_old)
     167             :       ELSE
     168          88 :          CALL get_rtp(rtp=rtp, mos_old=mos_old)
     169             :       END IF
     170         178 :       CALL put_data_to_history(rtp, mos=mos_old, s_mat=s_mat, ihist=1, rho=rho_old)
     171             : 
     172         178 :    END SUBROUTINE init_propagators
     173             : 
     174             : ! **************************************************************************************************
     175             : !> \brief gets an estimate for the 2-norm of KS (diagnaliztion of KS) and
     176             : !>        calculates the order and number of squaring steps for Taylor or
     177             : !>        Pade matrix exponential
     178             : !> \param rtp ...
     179             : !> \param s_mat ...
     180             : !> \param matrix_ks ...
     181             : !> \param rtp_control ...
     182             : !> \author Florian Schiffmann (02.09)
     183             : ! **************************************************************************************************
     184             : 
     185          66 :    SUBROUTINE get_maxabs_eigval(rtp, s_mat, matrix_ks, rtp_control)
     186             :       TYPE(rt_prop_type), POINTER                        :: rtp
     187             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: s_mat, matrix_ks
     188             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     189             : 
     190             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'get_maxabs_eigval'
     191             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     192             : 
     193             :       INTEGER                                            :: handle, ispin, method, ndim
     194             :       LOGICAL                                            :: emd
     195             :       REAL(dp)                                           :: max_eval, min_eval, norm2, scale, t
     196             :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: eigval_H
     197             :       TYPE(cp_fm_type)                                   :: eigvec_H, H_fm, S_half, S_inv_fm, &
     198             :                                                             S_minus_half, tmp, tmp_mat_H
     199             :       TYPE(dbcsr_type), POINTER                          :: S_inv
     200             : 
     201          66 :       CALL timeset(routineN, handle)
     202             : 
     203          66 :       CALL get_rtp(rtp=rtp, S_inv=S_inv, dt=t)
     204             : 
     205             :       CALL cp_fm_create(S_inv_fm, &
     206             :                         matrix_struct=rtp%ao_ao_fmstruct, &
     207          66 :                         name="S_inv")
     208          66 :       CALL copy_dbcsr_to_fm(S_inv, S_inv_fm)
     209             : 
     210             :       CALL cp_fm_create(S_half, &
     211             :                         matrix_struct=rtp%ao_ao_fmstruct, &
     212          66 :                         name="S_half")
     213             : 
     214             :       CALL cp_fm_create(S_minus_half, &
     215             :                         matrix_struct=rtp%ao_ao_fmstruct, &
     216          66 :                         name="S_minus_half")
     217             : 
     218             :       CALL cp_fm_create(H_fm, &
     219             :                         matrix_struct=rtp%ao_ao_fmstruct, &
     220          66 :                         name="RTP_H_FM")
     221             : 
     222             :       CALL cp_fm_create(tmp_mat_H, &
     223             :                         matrix_struct=rtp%ao_ao_fmstruct, &
     224          66 :                         name="TMP_H")
     225             : 
     226          66 :       ndim = S_inv_fm%matrix_struct%nrow_global
     227          66 :       scale = 1.0_dp
     228          66 :       IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
     229          66 :       t = -t/scale
     230             : 
     231             :       ! Create the overlap matrices
     232             : 
     233             :       CALL cp_fm_create(tmp, &
     234             :                         matrix_struct=rtp%ao_ao_fmstruct, &
     235          66 :                         name="tmp_mat")
     236             : 
     237             :       CALL cp_fm_create(eigvec_H, &
     238             :                         matrix_struct=rtp%ao_ao_fmstruct, &
     239          66 :                         name="tmp_EVEC")
     240             : 
     241         198 :       ALLOCATE (eigval_H(ndim))
     242          66 :       CALL copy_dbcsr_to_fm(s_mat(1)%matrix, tmp)
     243          66 :       CALL cp_fm_upper_to_full(tmp, eigvec_H)
     244             : 
     245          66 :       CALL cp_fm_syevd(tmp, eigvec_H, eigval_H)
     246             : 
     247         962 :       eigval_H(:) = one/eigval_H(:)
     248          66 :       CALL backtransform_matrix(eigval_H, eigvec_H, S_inv_fm)
     249         962 :       eigval_H(:) = SQRT(eigval_H(:))
     250          66 :       CALL backtransform_matrix(eigval_H, eigvec_H, S_minus_half)
     251         962 :       eigval_H(:) = one/eigval_H(:)
     252          66 :       CALL backtransform_matrix(eigval_H, eigvec_H, S_half)
     253          66 :       CALL cp_fm_release(eigvec_H)
     254          66 :       CALL cp_fm_release(tmp)
     255             : 
     256          66 :       IF (rtp_control%mat_exp == do_taylor) method = 1
     257          66 :       IF (rtp_control%mat_exp == do_pade) method = 2
     258          66 :       emd = (.NOT. rtp_control%fixed_ions)
     259             : 
     260         148 :       DO ispin = 1, SIZE(matrix_ks)
     261             : 
     262          82 :          CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, H_fm)
     263          82 :          CALL cp_fm_upper_to_full(H_fm, tmp_mat_H)
     264          82 :          CALL cp_fm_scale(t, H_fm)
     265             : 
     266             :          CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, H_fm, S_minus_half, zero, &
     267          82 :                             tmp_mat_H)
     268             :          CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, S_minus_half, tmp_mat_H, zero, &
     269          82 :                             H_fm)
     270             : 
     271          82 :          CALL cp_fm_syevd(H_fm, tmp_mat_H, eigval_H)
     272        1324 :          min_eval = MINVAL(eigval_H)
     273        1324 :          max_eval = MAXVAL(eigval_H)
     274          82 :          norm2 = 2.0_dp*MAX(ABS(min_eval), ABS(max_eval))
     275             :          CALL get_nsquare_norder(norm2, rtp%orders(1, ispin), rtp%orders(2, ispin), &
     276         148 :                                  rtp_control%eps_exp, method, emd)
     277             :       END DO
     278             : 
     279          66 :       DEALLOCATE (eigval_H)
     280             : 
     281          66 :       CALL copy_fm_to_dbcsr(S_inv_fm, S_inv)
     282          66 :       CALL cp_fm_release(S_inv_fm)
     283          66 :       CALL cp_fm_release(S_half)
     284          66 :       CALL cp_fm_release(S_minus_half)
     285          66 :       CALL cp_fm_release(H_fm)
     286          66 :       CALL cp_fm_release(tmp_mat_H)
     287             : 
     288          66 :       CALL timestop(handle)
     289             : 
     290          66 :    END SUBROUTINE get_maxabs_eigval
     291             : 
     292             : ! **************************************************************************************************
     293             : !> \brief gets an estimate for the 2-norm of KS (diagnaliztion of KS) and
     294             : !>        calculates the order and number of squaring steps for Taylor or
     295             : !>        Pade matrix exponential. Based on the full matrix code.
     296             : !> \param rtp ...
     297             : !> \param s_mat ...
     298             : !> \param matrix_ks ...
     299             : !> \param rtp_control ...
     300             : !> \author Samuel Andermatt (02.14)
     301             : ! **************************************************************************************************
     302             : 
     303         152 :    SUBROUTINE get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
     304             :       TYPE(rt_prop_type), POINTER                        :: rtp
     305             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: s_mat, matrix_ks
     306             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     307             : 
     308             :       CHARACTER(len=*), PARAMETER :: routineN = 'get_maxabs_eigval_sparse'
     309             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     310             : 
     311             :       INTEGER                                            :: handle, ispin, method
     312             :       LOGICAL                                            :: converged, emd
     313             :       REAL(dp)                                           :: max_ev, min_ev, norm2, scale, t
     314             :       TYPE(dbcsr_type), POINTER                          :: s_half, s_minus_half, tmp, tmp2
     315             : 
     316          76 :       CALL timeset(routineN, handle)
     317             : 
     318          76 :       CALL get_rtp(rtp=rtp, dt=t)
     319             : 
     320             :       NULLIFY (s_half)
     321          76 :       ALLOCATE (s_half)
     322          76 :       CALL dbcsr_create(s_half, template=s_mat(1)%matrix)
     323             :       NULLIFY (s_minus_half)
     324          76 :       ALLOCATE (s_minus_half)
     325          76 :       CALL dbcsr_create(s_minus_half, template=s_mat(1)%matrix)
     326             :       NULLIFY (tmp)
     327          76 :       ALLOCATE (tmp)
     328          76 :       CALL dbcsr_create(tmp, template=s_mat(1)%matrix, matrix_type="N")
     329             :       NULLIFY (tmp2)
     330          76 :       ALLOCATE (tmp2)
     331          76 :       CALL dbcsr_create(tmp2, template=s_mat(1)%matrix, matrix_type="N")
     332          76 :       scale = 1.0_dp
     333          76 :       IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
     334          76 :       t = -t/scale
     335          76 :       emd = (.NOT. rtp_control%fixed_ions)
     336             : 
     337          76 :       IF (rtp_control%mat_exp == do_taylor) method = 1
     338          76 :       IF (rtp_control%mat_exp == do_pade) method = 2
     339             :       CALL matrix_sqrt_Newton_Schulz(s_half, s_minus_half, s_mat(1)%matrix, rtp%filter_eps, &
     340          76 :                                      rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
     341         188 :       DO ispin = 1, SIZE(matrix_ks)
     342             :          CALL dbcsr_multiply("N", "N", t, matrix_ks(ispin)%matrix, s_minus_half, zero, tmp, &
     343         112 :                              filter_eps=rtp%filter_eps)
     344             :          CALL dbcsr_multiply("N", "N", one, s_minus_half, tmp, zero, tmp2, &
     345         112 :                              filter_eps=rtp%filter_eps)
     346             :          CALL arnoldi_extremal(tmp2, max_ev, min_ev, threshold=rtp%lanzcos_threshold, &
     347         112 :                                max_iter=rtp%lanzcos_max_iter, converged=converged)
     348         112 :          norm2 = 2.0_dp*MAX(ABS(min_ev), ABS(max_ev))
     349             :          CALL get_nsquare_norder(norm2, rtp%orders(1, ispin), rtp%orders(2, ispin), &
     350         188 :                                  rtp_control%eps_exp, method, emd)
     351             :       END DO
     352             : 
     353          76 :       CALL dbcsr_deallocate_matrix(s_half)
     354          76 :       CALL dbcsr_deallocate_matrix(s_minus_half)
     355          76 :       CALL dbcsr_deallocate_matrix(tmp)
     356          76 :       CALL dbcsr_deallocate_matrix(tmp2)
     357             : 
     358          76 :       CALL timestop(handle)
     359             : 
     360          76 :    END SUBROUTINE
     361             : 
     362             : ! **************************************************************************************************
     363             : !> \brief Is still left from diagonalization, should be removed later but is
     364             : !>  still needed for the guess for the pade/Taylor method
     365             : !> \param Eval ...
     366             : !> \param eigenvec ...
     367             : !> \param matrix ...
     368             : !> \author Florian Schiffmann (02.09)
     369             : ! **************************************************************************************************
     370             : 
     371         198 :    SUBROUTINE backtransform_matrix(Eval, eigenvec, matrix)
     372             : 
     373             :       REAL(dp), DIMENSION(:), INTENT(in)                 :: Eval
     374             :       TYPE(cp_fm_type), INTENT(IN)                       :: eigenvec, matrix
     375             : 
     376             :       CHARACTER(len=*), PARAMETER :: routineN = 'backtransform_matrix'
     377             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     378             : 
     379             :       INTEGER                                            :: handle, i, j, l, ncol_local, ndim, &
     380             :                                                             nrow_local
     381         198 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     382             :       TYPE(cp_fm_type)                                   :: tmp
     383             : 
     384         198 :       CALL timeset(routineN, handle)
     385             :       CALL cp_fm_create(tmp, &
     386             :                         matrix_struct=matrix%matrix_struct, &
     387         198 :                         name="TMP_BT")
     388             :       CALL cp_fm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
     389         198 :                           row_indices=row_indices, col_indices=col_indices)
     390             : 
     391         198 :       ndim = matrix%matrix_struct%nrow_global
     392             : 
     393         198 :       CALL cp_fm_set_all(tmp, zero, zero)
     394        2886 :       DO i = 1, ncol_local
     395        2688 :          l = col_indices(i)
     396       25056 :          DO j = 1, nrow_local
     397       24858 :             tmp%local_data(j, i) = eigenvec%local_data(j, i)*Eval(l)
     398             :          END DO
     399             :       END DO
     400             :       CALL parallel_gemm("N", "T", ndim, ndim, ndim, one, tmp, eigenvec, zero, &
     401         198 :                          matrix)
     402             : 
     403         198 :       CALL cp_fm_release(tmp)
     404         198 :       CALL timestop(handle)
     405             : 
     406         198 :    END SUBROUTINE backtransform_matrix
     407             : 
     408             : ! **************************************************************************************************
     409             : !> \brief Computes the density matrix from the mos
     410             : !>        Update: Initialized the density matrix from complex mos (for
     411             : !>        instance after delta kick)
     412             : !> \param rtp ...
     413             : !> \param mos ...
     414             : !> \param mos_old ...
     415             : !> \author Samuel Andermatt (08.15)
     416             : !> \author Guillaume Le Breton (01.23)
     417             : ! **************************************************************************************************
     418             : 
     419          64 :    SUBROUTINE rt_initialize_rho_from_mos(rtp, mos, mos_old)
     420             : 
     421             :       TYPE(rt_prop_type), POINTER                        :: rtp
     422             :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     423             :       TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER  :: mos_old
     424             : 
     425             :       INTEGER                                            :: im, ispin, ncol, re
     426             :       REAL(KIND=dp)                                      :: alpha
     427             :       TYPE(cp_fm_type)                                   :: mos_occ, mos_occ_im
     428          64 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_new, rho_old
     429             : 
     430          64 :       CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
     431             : 
     432          64 :       IF (PRESENT(mos_old)) THEN
     433             :          ! Used the mos from delta kick. Initialize both real and im part
     434         106 :          DO ispin = 1, SIZE(mos_old)/2
     435          66 :             re = 2*ispin - 1; im = 2*ispin
     436          66 :             CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
     437          66 :             CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
     438             :             CALL cp_fm_create(mos_occ, &
     439             :                               matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
     440          66 :                               name="mos_occ")
     441             :             !Real part of rho
     442          66 :             CALL cp_fm_to_fm(mos_old(re), mos_occ)
     443          66 :             IF (mos(ispin)%uniform_occupation) THEN
     444          58 :                alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
     445         280 :                CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
     446             :                CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
     447             :                                           matrix_v=mos_occ, &
     448             :                                           ncol=ncol, &
     449          58 :                                           alpha=alpha, keep_sparsity=.FALSE.)
     450             :             ELSE
     451           8 :                alpha = 1.0_dp
     452         116 :                CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
     453             :                CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
     454             :                                           matrix_v=mos_old(re), &
     455             :                                           matrix_g=mos_occ, &
     456             :                                           ncol=ncol, &
     457           8 :                                           alpha=alpha, keep_sparsity=.FALSE.)
     458             :             END IF
     459             : 
     460             :             ! Add complex part of MOs, i*i=-1
     461          66 :             CALL cp_fm_to_fm(mos_old(im), mos_occ)
     462          66 :             IF (mos(ispin)%uniform_occupation) THEN
     463          58 :                alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
     464         280 :                CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
     465             :                CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
     466             :                                           matrix_v=mos_occ, &
     467             :                                           ncol=ncol, &
     468          58 :                                           alpha=alpha, keep_sparsity=.FALSE.)
     469             :             ELSE
     470           8 :                alpha = 1.0_dp
     471         116 :                CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
     472             :                CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
     473             :                                           matrix_v=mos_old(im), &
     474             :                                           matrix_g=mos_occ, &
     475             :                                           ncol=ncol, &
     476           8 :                                           alpha=alpha, keep_sparsity=.FALSE.)
     477             :             END IF
     478          66 :             CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
     479          66 :             CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
     480             : 
     481             :             ! Imaginary part of rho
     482             :             CALL cp_fm_create(mos_occ_im, &
     483             :                               matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
     484          66 :                               name="mos_occ_im")
     485          66 :             alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
     486          66 :             CALL cp_fm_to_fm(mos_old(re), mos_occ)
     487         396 :             CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
     488          66 :             CALL cp_fm_to_fm(mos_old(im), mos_occ_im)
     489         396 :             CALL cp_fm_column_scale(mos_occ_im, mos(ispin)%occupation_numbers/alpha)
     490             :             CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(im)%matrix, &
     491             :                                        matrix_v=mos_occ_im, &
     492             :                                        matrix_g=mos_occ, &
     493             :                                        ncol=ncol, &
     494             :                                        alpha=2.0_dp*alpha, &
     495          66 :                                        symmetry_mode=-1, keep_sparsity=.FALSE.)
     496             : 
     497          66 :             CALL dbcsr_filter(rho_old(im)%matrix, rtp%filter_eps)
     498          66 :             CALL dbcsr_copy(rho_new(im)%matrix, rho_old(im)%matrix)
     499          66 :             CALL cp_fm_release(mos_occ_im)
     500         172 :             CALL cp_fm_release(mos_occ)
     501             :          END DO
     502             :          ! Release the mos used to apply the delta kick, no longer required
     503          40 :          CALL rt_prop_release_mos(rtp)
     504             :       ELSE
     505          50 :          DO ispin = 1, SIZE(mos)
     506          26 :             re = 2*ispin - 1
     507          26 :             CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
     508          26 :             CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
     509             : 
     510             :             CALL cp_fm_create(mos_occ, &
     511             :                               matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
     512          26 :                               name="mos_occ")
     513          26 :             CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
     514          26 :             IF (mos(ispin)%uniform_occupation) THEN
     515          22 :                alpha = 3.0_dp - REAL(SIZE(mos), dp)
     516         110 :                CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
     517             :                CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
     518             :                                           matrix_v=mos_occ, &
     519             :                                           ncol=ncol, &
     520          22 :                                           alpha=alpha, keep_sparsity=.FALSE.)
     521             :             ELSE
     522           4 :                alpha = 1.0_dp
     523          88 :                CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
     524             :                CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
     525             :                                           matrix_v=mos(ispin)%mo_coeff, &
     526             :                                           matrix_g=mos_occ, &
     527             :                                           ncol=ncol, &
     528           4 :                                           alpha=alpha, keep_sparsity=.FALSE.)
     529             :             END IF
     530          26 :             CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
     531          26 :             CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
     532          76 :             CALL cp_fm_release(mos_occ)
     533             :          END DO
     534             :       END IF
     535             : 
     536          64 :    END SUBROUTINE rt_initialize_rho_from_mos
     537             : 
     538             : END MODULE rt_propagator_init

Generated by: LCOV version 1.15