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

Generated by: LCOV version 1.15