LCOV - code coverage report
Current view: top level - src - ls_matrix_exp.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 202 209 96.7 %
Date: 2025-02-18 08:24:35 Functions: 5 5 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines for calculating a complex matrix exponential with dbcsr matrices.
      10             : !>        Based on the code in matrix_exp.F from Florian Schiffmann
      11             : !> \author Samuel Andermatt (02.14)
      12             : ! **************************************************************************************************
      13             : MODULE ls_matrix_exp
      14             : 
      15             :    USE cp_dbcsr_api,                    ONLY: &
      16             :         dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, &
      17             :         dbcsr_filter, dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, &
      18             :         dbcsr_transposed, dbcsr_type
      19             :    USE cp_dbcsr_contrib,                ONLY: dbcsr_frobenius_norm
      20             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      21             :    USE kinds,                           ONLY: dp
      22             : #include "./base/base_uses.f90"
      23             : 
      24             :    IMPLICIT NONE
      25             : 
      26             :    PRIVATE
      27             : 
      28             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ls_matrix_exp'
      29             : 
      30             :    PUBLIC :: taylor_only_imaginary_dbcsr, &
      31             :              taylor_full_complex_dbcsr, &
      32             :              cp_complex_dbcsr_gemm_3, &
      33             :              bch_expansion_imaginary_propagator, &
      34             :              bch_expansion_complex_propagator
      35             : 
      36             : CONTAINS
      37             : 
      38             : ! **************************************************************************************************
      39             : !> \brief Convenience function. Computes the matrix multiplications needed
      40             : !>        for the multiplication of complex sparse matrices.
      41             : !>        C = beta * C + alpha * ( A  ** transa ) * ( B ** transb )
      42             : !> \param transa : 'N' -> normal   'T' -> transpose
      43             : !>      alpha,beta :: can be 0.0_dp and 1.0_dp
      44             : !> \param transb ...
      45             : !> \param alpha ...
      46             : !> \param A_re m x k matrix ( ! for transa = 'N'), real part
      47             : !> \param A_im m x k matrix ( ! for transa = 'N'), imaginary part
      48             : !> \param B_re k x n matrix ( ! for transb = 'N'), real part
      49             : !> \param B_im k x n matrix ( ! for transb = 'N'), imaginary part
      50             : !> \param beta ...
      51             : !> \param C_re m x n matrix, real part
      52             : !> \param C_im m x n matrix, imaginary part
      53             : !> \param filter_eps ...
      54             : !> \author Samuel Andermatt
      55             : !> \note
      56             : !>      C should have no overlap with A, B
      57             : !>      This subroutine uses three real matrix multiplications instead of two complex
      58             : !>      This reduces the amount of flops and memory bandwidth by 25%, but for memory bandwidth
      59             : !>      true complex algebra is still superior (one third less bandwidth needed)
      60             : !>      limited cases matrix multiplications
      61             : ! **************************************************************************************************
      62        6962 :    SUBROUTINE cp_complex_dbcsr_gemm_3(transa, transb, alpha, A_re, A_im, &
      63             :                                       B_re, B_im, beta, C_re, C_im, filter_eps)
      64             :       CHARACTER(LEN=1), INTENT(IN)                       :: transa, transb
      65             :       REAL(KIND=dp), INTENT(IN)                          :: alpha
      66             :       TYPE(dbcsr_type), INTENT(IN)                       :: A_re, A_im, B_re, B_im
      67             :       REAL(KIND=dp), INTENT(IN)                          :: beta
      68             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: C_re, C_im
      69             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: filter_eps
      70             : 
      71             :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_complex_dbcsr_gemm_3'
      72             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
      73             : 
      74             :       CHARACTER(LEN=1)                                   :: transa2, transb2
      75             :       INTEGER                                            :: handle
      76             :       REAL(KIND=dp)                                      :: alpha2, alpha3, alpha4
      77             :       TYPE(dbcsr_type), POINTER                          :: a_plus_b, ac, bd, c_plus_d
      78             : 
      79        6962 :       CALL timeset(routineN, handle)
      80             :       !A complex matrix matrix multiplication can be done with only three multiplications
      81             :       !(a+ib)*(c+id)=ac-bd+i((a+b)*(c+d) - ac - bd)
      82             :       !A_re=a, A_im=b, B_re=c, B_im=d
      83             : 
      84        6962 :       alpha2 = -alpha
      85        6962 :       alpha3 = alpha
      86        6962 :       alpha4 = alpha
      87             : 
      88        6962 :       IF (transa == "C") THEN
      89           0 :          alpha2 = -alpha2
      90           0 :          alpha3 = -alpha3
      91           0 :          transa2 = "T"
      92             :       ELSE
      93        6962 :          transa2 = transa
      94             :       END IF
      95        6962 :       IF (transb == "C") THEN
      96        1072 :          alpha2 = -alpha2
      97        1072 :          alpha4 = -alpha4
      98        1072 :          transb2 = "T"
      99             :       ELSE
     100        5890 :          transb2 = transb
     101             :       END IF
     102             : 
     103             :       !create the work matrices
     104             :       NULLIFY (ac)
     105        6962 :       ALLOCATE (ac)
     106        6962 :       CALL dbcsr_create(ac, template=A_re, matrix_type="N")
     107             :       NULLIFY (bd)
     108        6962 :       ALLOCATE (bd)
     109        6962 :       CALL dbcsr_create(bd, template=A_re, matrix_type="N")
     110             :       NULLIFY (a_plus_b)
     111        6962 :       ALLOCATE (a_plus_b)
     112        6962 :       CALL dbcsr_create(a_plus_b, template=A_re, matrix_type="N")
     113             :       NULLIFY (c_plus_d)
     114        6962 :       ALLOCATE (c_plus_d)
     115        6962 :       CALL dbcsr_create(c_plus_d, template=A_re, matrix_type="N")
     116             : 
     117             :       !Do the neccesarry multiplications
     118        6962 :       CALL dbcsr_multiply(transa2, transb2, alpha, A_re, B_re, zero, ac, filter_eps=filter_eps)
     119        6962 :       CALL dbcsr_multiply(transa2, transb2, alpha2, A_im, B_im, zero, bd, filter_eps=filter_eps)
     120             : 
     121        6962 :       CALL dbcsr_add(a_plus_b, A_re, zero, alpha)
     122        6962 :       CALL dbcsr_add(a_plus_b, A_im, one, alpha3)
     123        6962 :       CALL dbcsr_add(c_plus_d, B_re, zero, alpha)
     124        6962 :       CALL dbcsr_add(c_plus_d, B_im, one, alpha4)
     125             : 
     126             :       !this can already be written into C_im
     127             :       !now both matrixes have been scaled which means we currently multiplied by alpha squared
     128        6962 :       CALL dbcsr_multiply(transa2, transb2, one/alpha, a_plus_b, c_plus_d, beta, C_im, filter_eps=filter_eps)
     129             : 
     130             :       !now add up all the terms into the result
     131        6962 :       CALL dbcsr_add(C_re, ac, beta, one)
     132             :       !the minus sign was already taken care of at the definition of alpha2
     133        6962 :       CALL dbcsr_add(C_re, bd, one, one)
     134        6962 :       CALL dbcsr_filter(C_re, filter_eps)
     135             : 
     136        6962 :       CALL dbcsr_add(C_im, ac, one, -one)
     137             :       !the minus sign was already taken care of at the definition of alpha2
     138        6962 :       CALL dbcsr_add(C_im, bd, one, one)
     139        6962 :       CALL dbcsr_filter(C_im, filter_eps)
     140             : 
     141             :       !Deallocate the work matrices
     142        6962 :       CALL dbcsr_deallocate_matrix(ac)
     143        6962 :       CALL dbcsr_deallocate_matrix(bd)
     144        6962 :       CALL dbcsr_deallocate_matrix(a_plus_b)
     145        6962 :       CALL dbcsr_deallocate_matrix(c_plus_d)
     146             : 
     147        6962 :       CALL timestop(handle)
     148             : 
     149        6962 :    END SUBROUTINE cp_complex_dbcsr_gemm_3
     150             : 
     151             : ! **************************************************************************************************
     152             : !> \brief specialized subroutine for purely imaginary matrix exponentials
     153             : !> \param exp_H ...
     154             : !> \param im_matrix ...
     155             : !> \param nsquare ...
     156             : !> \param ntaylor ...
     157             : !> \param filter_eps ...
     158             : !> \author Samuel Andermatt (02.2014)
     159             : ! **************************************************************************************************
     160         708 :    SUBROUTINE taylor_only_imaginary_dbcsr(exp_H, im_matrix, nsquare, ntaylor, filter_eps)
     161             : 
     162             :       TYPE(dbcsr_p_type), DIMENSION(2)                   :: exp_H
     163             :       TYPE(dbcsr_type), POINTER                          :: im_matrix
     164             :       INTEGER, INTENT(in)                                :: nsquare, ntaylor
     165             :       REAL(KIND=dp), INTENT(in)                          :: filter_eps
     166             : 
     167             :       CHARACTER(len=*), PARAMETER :: routineN = 'taylor_only_imaginary_dbcsr'
     168             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     169             : 
     170             :       INTEGER                                            :: handle, i, nloop
     171             :       REAL(KIND=dp)                                      :: square_fac, Tfac, tmp
     172             :       TYPE(dbcsr_type), POINTER                          :: T1, T2, Tres_im, Tres_re
     173             : 
     174         708 :       CALL timeset(routineN, handle)
     175             : 
     176             :       !The divider that comes from the scaling and squaring procedure
     177         708 :       square_fac = 1.0_dp/(2.0_dp**REAL(nsquare, dp))
     178             : 
     179             :       !Allocate work matrices
     180             :       NULLIFY (T1)
     181         708 :       ALLOCATE (T1)
     182         708 :       CALL dbcsr_create(T1, template=im_matrix, matrix_type="N")
     183             :       NULLIFY (T2)
     184         708 :       ALLOCATE (T2)
     185         708 :       CALL dbcsr_create(T2, template=im_matrix, matrix_type="N")
     186             :       NULLIFY (Tres_re)
     187         708 :       ALLOCATE (Tres_re)
     188         708 :       CALL dbcsr_create(Tres_re, template=im_matrix, matrix_type="N")
     189             :       NULLIFY (Tres_im)
     190         708 :       ALLOCATE (Tres_im)
     191         708 :       CALL dbcsr_create(Tres_im, template=im_matrix, matrix_type="N")
     192             : 
     193             :       !Create the unit matrices
     194         708 :       CALL dbcsr_set(T1, zero)
     195         708 :       CALL dbcsr_add_on_diag(T1, one)
     196         708 :       CALL dbcsr_set(Tres_re, zero)
     197         708 :       CALL dbcsr_add_on_diag(Tres_re, one)
     198         708 :       CALL dbcsr_set(Tres_im, zero)
     199             : 
     200         708 :       nloop = CEILING(REAL(ntaylor, dp)/2.0_dp)
     201             :       !the inverse of the prefactor in the taylor series
     202         708 :       tmp = 1.0_dp
     203        3516 :       DO i = 1, nloop
     204        2808 :          CALL dbcsr_scale(T1, 1.0_dp/(REAL(i, dp)*2.0_dp - 1.0_dp))
     205        2808 :          CALL dbcsr_filter(T1, filter_eps)
     206             :          CALL dbcsr_multiply("N", "N", square_fac, im_matrix, T1, zero, &
     207        2808 :                              T2, filter_eps=filter_eps)
     208        2808 :          Tfac = one
     209        2808 :          IF (MOD(i, 2) == 0) Tfac = -Tfac
     210        2808 :          CALL dbcsr_add(Tres_im, T2, one, Tfac)
     211        2808 :          CALL dbcsr_scale(T2, 1.0_dp/(REAL(i, dp)*2.0_dp))
     212        2808 :          CALL dbcsr_filter(T2, filter_eps)
     213             :          CALL dbcsr_multiply("N", "N", square_fac, im_matrix, T2, zero, &
     214        2808 :                              T1, filter_eps=filter_eps)
     215        2808 :          Tfac = one
     216        2808 :          IF (MOD(i, 2) == 1) Tfac = -Tfac
     217        3516 :          CALL dbcsr_add(Tres_re, T1, one, Tfac)
     218             :       END DO
     219             : 
     220             :       !Square the matrices, due to the scaling and squaring procedure
     221         708 :       IF (nsquare .GT. 0) THEN
     222           0 :          DO i = 1, nsquare
     223             :             CALL cp_complex_dbcsr_gemm_3("N", "N", one, Tres_re, Tres_im, &
     224             :                                          Tres_re, Tres_im, zero, exp_H(1)%matrix, exp_H(2)%matrix, &
     225           0 :                                          filter_eps=filter_eps)
     226           0 :             CALL dbcsr_copy(Tres_re, exp_H(1)%matrix)
     227           0 :             CALL dbcsr_copy(Tres_im, exp_H(2)%matrix)
     228             :          END DO
     229             :       ELSE
     230         708 :          CALL dbcsr_copy(exp_H(1)%matrix, Tres_re)
     231         708 :          CALL dbcsr_copy(exp_H(2)%matrix, Tres_im)
     232             :       END IF
     233         708 :       CALL dbcsr_deallocate_matrix(T1)
     234         708 :       CALL dbcsr_deallocate_matrix(T2)
     235         708 :       CALL dbcsr_deallocate_matrix(Tres_re)
     236         708 :       CALL dbcsr_deallocate_matrix(Tres_im)
     237             : 
     238         708 :       CALL timestop(handle)
     239             : 
     240         708 :    END SUBROUTINE taylor_only_imaginary_dbcsr
     241             : 
     242             : ! **************************************************************************************************
     243             : !> \brief subroutine for general complex matrix exponentials
     244             : !>        on input a separate dbcsr_type for real and complex part
     245             : !>        on output a size 2 dbcsr_p_type, first element is the real part of
     246             : !>        the exponential second the imaginary
     247             : !> \param exp_H ...
     248             : !> \param re_part ...
     249             : !> \param im_part ...
     250             : !> \param nsquare ...
     251             : !> \param ntaylor ...
     252             : !> \param filter_eps ...
     253             : !> \author Samuel Andermatt (02.2014)
     254             : ! **************************************************************************************************
     255         264 :    SUBROUTINE taylor_full_complex_dbcsr(exp_H, re_part, im_part, nsquare, ntaylor, filter_eps)
     256             :       TYPE(dbcsr_p_type), DIMENSION(2)                   :: exp_H
     257             :       TYPE(dbcsr_type), POINTER                          :: re_part, im_part
     258             :       INTEGER, INTENT(in)                                :: nsquare, ntaylor
     259             :       REAL(KIND=dp), INTENT(in)                          :: filter_eps
     260             : 
     261             :       CHARACTER(len=*), PARAMETER :: routineN = 'taylor_full_complex_dbcsr'
     262             : 
     263             :       INTEGER                                            :: handle, i
     264             :       REAL(KIND=dp)                                      :: square_fac
     265             :       TYPE(dbcsr_type)                                   :: T1_im, T1_re, T2_im, T2_re
     266             :       TYPE(dbcsr_type), POINTER                          :: Tres_im, Tres_re
     267             : 
     268         264 :       CALL timeset(routineN, handle)
     269             : 
     270             :       ! convenient aliases for result matrices
     271         264 :       Tres_re => exp_H(1)%matrix
     272         264 :       Tres_im => exp_H(2)%matrix
     273             : 
     274             :       ! The divider that comes from the scaling and squaring procedure
     275         264 :       square_fac = 1.0_dp/REAL(2**nsquare, dp)
     276             : 
     277             :       ! Allocate work matrices
     278         264 :       CALL dbcsr_create(T1_re, template=re_part, matrix_type="N")
     279         264 :       CALL dbcsr_create(T1_im, template=re_part, matrix_type="N")
     280         264 :       CALL dbcsr_create(T2_re, template=re_part, matrix_type="N")
     281         264 :       CALL dbcsr_create(T2_im, template=re_part, matrix_type="N")
     282             : 
     283             :       ! T1 = identity
     284         264 :       CALL dbcsr_set(T1_re, 0.0_dp)
     285         264 :       CALL dbcsr_set(T1_im, 0.0_dp)
     286         264 :       CALL dbcsr_add_on_diag(T1_re, 1.0_dp)
     287             : 
     288             :       ! Tres = identity
     289         264 :       CALL dbcsr_set(Tres_re, 0.0_dp)
     290         264 :       CALL dbcsr_set(Tres_im, 0.0_dp)
     291         264 :       CALL dbcsr_add_on_diag(Tres_re, 1.0_dp)
     292             : 
     293        2522 :       DO i = 1, ntaylor
     294             :          ! T1 = T1 / i
     295        2258 :          CALL dbcsr_scale(T1_re, 1.0_dp/REAL(i, dp))
     296        2258 :          CALL dbcsr_scale(T1_im, 1.0_dp/REAL(i, dp))
     297        2258 :          CALL dbcsr_filter(T1_re, filter_eps)
     298        2258 :          CALL dbcsr_filter(T1_im, filter_eps)
     299             : 
     300             :          ! T2 = square_fac * T1 * input
     301             :          CALL cp_complex_dbcsr_gemm_3("N", "N", alpha=square_fac, beta=0.0_dp, &
     302             :                                       A_re=T1_re, A_im=T1_im, &
     303             :                                       B_re=re_part, B_im=im_part, &
     304        2258 :                                       C_re=T2_re, C_im=T2_im, filter_eps=filter_eps)
     305             : 
     306             :          ! Tres = Tres + T2
     307        2258 :          CALL dbcsr_add(Tres_re, T2_re, 1.0_dp, 1.0_dp)
     308        2258 :          CALL dbcsr_add(Tres_im, T2_im, 1.0_dp, 1.0_dp)
     309             : 
     310             :          ! T1 = T2
     311        2258 :          CALL dbcsr_copy(T1_re, T2_re)
     312        2522 :          CALL dbcsr_copy(T1_im, T2_im)
     313             :       END DO
     314             : 
     315         264 :       IF (nsquare .GT. 0) THEN
     316         970 :          DO i = 1, nsquare
     317             :             ! T2 = Tres * Tres
     318             :             CALL cp_complex_dbcsr_gemm_3("N", "N", alpha=1.0_dp, beta=0.0_dp, &
     319             :                                          A_re=Tres_re, A_im=Tres_im, &
     320             :                                          B_re=Tres_re, B_im=Tres_im, &
     321         758 :                                          C_re=T2_re, C_im=T2_im, filter_eps=filter_eps)
     322             : 
     323             :             ! Tres = T2
     324         758 :             CALL dbcsr_copy(Tres_re, T2_re)
     325         970 :             CALL dbcsr_copy(Tres_im, T2_im)
     326             :          END DO
     327             :       END IF
     328             : 
     329         264 :       CALL dbcsr_release(T1_re)
     330         264 :       CALL dbcsr_release(T1_im)
     331         264 :       CALL dbcsr_release(T2_re)
     332         264 :       CALL dbcsr_release(T2_im)
     333             : 
     334         264 :       CALL timestop(handle)
     335             : 
     336         264 :    END SUBROUTINE taylor_full_complex_dbcsr
     337             : 
     338             : ! **************************************************************************************************
     339             : !> \brief  The Baker-Campbell-Hausdorff expansion for a purely imaginary exponent (e.g. rtp)
     340             : !>         Works for a non unitary propagator, because the density matrix is hermitian
     341             : !> \param propagator The exponent of the matrix exponential
     342             : !> \param density_re Real part of the density matrix
     343             : !> \param density_im Imaginary part of the density matrix
     344             : !> \param filter_eps The filtering threshold for all matrices
     345             : !> \param filter_eps_small ...
     346             : !> \param eps_exp The accuracy of the exponential
     347             : !> \author Samuel Andermatt (02.2014)
     348             : ! **************************************************************************************************
     349         112 :    SUBROUTINE bch_expansion_imaginary_propagator(propagator, density_re, density_im, filter_eps, filter_eps_small, eps_exp)
     350             :       TYPE(dbcsr_type), POINTER                          :: propagator, density_re, density_im
     351             :       REAL(KIND=dp), INTENT(in)                          :: filter_eps, filter_eps_small, eps_exp
     352             : 
     353             :       CHARACTER(len=*), PARAMETER :: routineN = 'bch_expansion_imaginary_propagator'
     354             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     355             : 
     356             :       INTEGER                                            :: handle, i, unit_nr
     357             :       LOGICAL                                            :: convergence
     358             :       REAL(KIND=dp)                                      :: alpha, max_alpha, prefac
     359             :       TYPE(dbcsr_type), POINTER                          :: comm, comm2, tmp, tmp2
     360             : 
     361         112 :       CALL timeset(routineN, handle)
     362             : 
     363         112 :       unit_nr = cp_logger_get_default_io_unit()
     364             : 
     365             :       NULLIFY (tmp)
     366         112 :       ALLOCATE (tmp)
     367         112 :       CALL dbcsr_create(tmp, template=propagator)
     368             :       NULLIFY (tmp2)
     369         112 :       ALLOCATE (tmp2)
     370         112 :       CALL dbcsr_create(tmp2, template=propagator)
     371             :       NULLIFY (comm)
     372         112 :       ALLOCATE (comm)
     373         112 :       CALL dbcsr_create(comm, template=propagator)
     374             :       NULLIFY (comm2)
     375         112 :       ALLOCATE (comm2)
     376         112 :       CALL dbcsr_create(comm2, template=propagator)
     377             : 
     378         112 :       CALL dbcsr_copy(tmp, density_re)
     379         112 :       CALL dbcsr_copy(tmp2, density_im)
     380             : 
     381         112 :       convergence = .FALSE.
     382         494 :       DO i = 1, 20
     383         494 :          prefac = one/i
     384             :          CALL dbcsr_multiply("N", "N", -prefac, propagator, tmp2, zero, comm, &
     385         494 :                              filter_eps=filter_eps_small)
     386             :          CALL dbcsr_multiply("N", "N", prefac, propagator, tmp, zero, comm2, &
     387         494 :                              filter_eps=filter_eps_small)
     388         494 :          CALL dbcsr_transposed(tmp, comm)
     389         494 :          CALL dbcsr_transposed(tmp2, comm2)
     390         494 :          CALL dbcsr_add(comm, tmp, one, one)
     391         494 :          CALL dbcsr_add(comm2, tmp2, one, -one)
     392         494 :          CALL dbcsr_add(density_re, comm, one, one)
     393         494 :          CALL dbcsr_add(density_im, comm2, one, one)
     394         494 :          CALL dbcsr_copy(tmp, comm)
     395         494 :          CALL dbcsr_copy(tmp2, comm2)
     396             :          !check for convergence
     397         494 :          max_alpha = zero
     398         494 :          alpha = dbcsr_frobenius_norm(comm)
     399         494 :          IF (alpha > max_alpha) max_alpha = alpha
     400         494 :          alpha = dbcsr_frobenius_norm(comm2)
     401         494 :          IF (alpha > max_alpha) max_alpha = alpha
     402         494 :          IF (max_alpha < eps_exp) convergence = .TRUE.
     403         382 :          IF (convergence) THEN
     404         112 :             IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="((T3,A,I2,A))") &
     405          56 :                "BCH converged after ", i, " steps"
     406             :             EXIT
     407             :          END IF
     408             :       END DO
     409             : 
     410         112 :       CALL dbcsr_filter(density_re, filter_eps)
     411         112 :       CALL dbcsr_filter(density_im, filter_eps)
     412             : 
     413         112 :       CPWARN_IF(.NOT. convergence, "BCH method did not converge")
     414             : 
     415         112 :       CALL dbcsr_deallocate_matrix(tmp)
     416         112 :       CALL dbcsr_deallocate_matrix(tmp2)
     417         112 :       CALL dbcsr_deallocate_matrix(comm)
     418         112 :       CALL dbcsr_deallocate_matrix(comm2)
     419             : 
     420         112 :       CALL timestop(handle)
     421             : 
     422         112 :    END SUBROUTINE bch_expansion_imaginary_propagator
     423             : 
     424             : ! **************************************************************************************************
     425             : !> \brief  The Baker-Campbell-Hausdorff expansion for a complex exponent (e.g. rtp)
     426             : !>         Works for a non unitary propagator, because the density matrix is hermitian
     427             : !> \param propagator_re Real part of the exponent
     428             : !> \param propagator_im Imaginary part of the exponent
     429             : !> \param density_re Real part of the density matrix
     430             : !> \param density_im Imaginary part of the density matrix
     431             : !> \param filter_eps The filtering threshold for all matrices
     432             : !> \param filter_eps_small ...
     433             : !> \param eps_exp The accuracy of the exponential
     434             : !> \author Samuel Andermatt (02.2014)
     435             : ! **************************************************************************************************
     436         130 :    SUBROUTINE bch_expansion_complex_propagator(propagator_re, propagator_im, density_re, density_im, &
     437             :                                                filter_eps, filter_eps_small, eps_exp)
     438             :       TYPE(dbcsr_type), POINTER                          :: propagator_re, propagator_im, &
     439             :                                                             density_re, density_im
     440             :       REAL(KIND=dp), INTENT(in)                          :: filter_eps, filter_eps_small, eps_exp
     441             : 
     442             :       CHARACTER(len=*), PARAMETER :: routineN = 'bch_expansion_complex_propagator'
     443             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     444             : 
     445             :       INTEGER                                            :: handle, i, unit_nr
     446             :       LOGICAL                                            :: convergence
     447             :       REAL(KIND=dp)                                      :: alpha, max_alpha, prefac
     448             :       TYPE(dbcsr_type), POINTER                          :: comm, comm2, tmp, tmp2
     449             : 
     450         130 :       CALL timeset(routineN, handle)
     451             : 
     452         130 :       unit_nr = cp_logger_get_default_io_unit()
     453             : 
     454             :       NULLIFY (tmp)
     455         130 :       ALLOCATE (tmp)
     456         130 :       CALL dbcsr_create(tmp, template=propagator_re)
     457             :       NULLIFY (tmp2)
     458         130 :       ALLOCATE (tmp2)
     459         130 :       CALL dbcsr_create(tmp2, template=propagator_re)
     460             :       NULLIFY (comm)
     461         130 :       ALLOCATE (comm)
     462         130 :       CALL dbcsr_create(comm, template=propagator_re)
     463             :       NULLIFY (comm2)
     464         130 :       ALLOCATE (comm2)
     465         130 :       CALL dbcsr_create(comm2, template=propagator_re)
     466             : 
     467         130 :       CALL dbcsr_copy(tmp, density_re)
     468         130 :       CALL dbcsr_copy(tmp2, density_im)
     469             : 
     470         130 :       convergence = .FALSE.
     471             : 
     472        1254 :       DO i = 1, 20
     473        1254 :          prefac = one/i
     474             :          CALL cp_complex_dbcsr_gemm_3("N", "N", prefac, propagator_re, propagator_im, &
     475        1254 :                                       tmp, tmp2, zero, comm, comm2, filter_eps=filter_eps_small)
     476        1254 :          CALL dbcsr_transposed(tmp, comm)
     477        1254 :          CALL dbcsr_transposed(tmp2, comm2)
     478        1254 :          CALL dbcsr_add(comm, tmp, one, one)
     479        1254 :          CALL dbcsr_add(comm2, tmp2, one, -one)
     480        1254 :          CALL dbcsr_add(density_re, comm, one, one)
     481        1254 :          CALL dbcsr_add(density_im, comm2, one, one)
     482        1254 :          CALL dbcsr_copy(tmp, comm)
     483        1254 :          CALL dbcsr_copy(tmp2, comm2)
     484             :          !check for convergence
     485        1254 :          max_alpha = zero
     486        1254 :          alpha = dbcsr_frobenius_norm(comm)
     487        1254 :          IF (alpha > max_alpha) max_alpha = alpha
     488        1254 :          alpha = dbcsr_frobenius_norm(comm2)
     489        1254 :          IF (alpha > max_alpha) max_alpha = alpha
     490        1254 :          IF (max_alpha < eps_exp) convergence = .TRUE.
     491        1124 :          IF (convergence) THEN
     492         130 :             IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="((T3,A,I2,A))") &
     493          65 :                "BCH converged after ", i, " steps"
     494             :             EXIT
     495             :          END IF
     496             :       END DO
     497             : 
     498         130 :       CALL dbcsr_filter(density_re, filter_eps)
     499         130 :       CALL dbcsr_filter(density_im, filter_eps)
     500             : 
     501         130 :       CPWARN_IF(.NOT. convergence, "BCH method did not converge")
     502             : 
     503         130 :       CALL dbcsr_deallocate_matrix(tmp)
     504         130 :       CALL dbcsr_deallocate_matrix(tmp2)
     505         130 :       CALL dbcsr_deallocate_matrix(comm)
     506         130 :       CALL dbcsr_deallocate_matrix(comm2)
     507             : 
     508         130 :       CALL timestop(handle)
     509             : 
     510         130 :    END SUBROUTINE bch_expansion_complex_propagator
     511             : 
     512             : END MODULE ls_matrix_exp

Generated by: LCOV version 1.15