LCOV - code coverage report
Current view: top level - src - ls_matrix_exp.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:f7c1873) Lines: 204 213 95.8 %
Date: 2024-03-26 06:29:31 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 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             : 
      14             : MODULE ls_matrix_exp
      15             : 
      16             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      17             :    USE dbcsr_api,                       ONLY: &
      18             :         dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, &
      19             :         dbcsr_filter, dbcsr_frobenius_norm, dbcsr_multiply, dbcsr_p_type, dbcsr_scale, dbcsr_set, &
      20             :         dbcsr_transposed, dbcsr_type, dbcsr_type_complex_8
      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             : 
      63        3946 :    SUBROUTINE cp_complex_dbcsr_gemm_3(transa, transb, alpha, A_re, A_im, &
      64             :                                       B_re, B_im, beta, C_re, C_im, filter_eps)
      65             :       CHARACTER(LEN=1), INTENT(IN)                       :: transa, transb
      66             :       REAL(KIND=dp), INTENT(IN)                          :: alpha
      67             :       TYPE(dbcsr_type), INTENT(IN)                       :: A_re, A_im, B_re, B_im
      68             :       REAL(KIND=dp), INTENT(IN)                          :: beta
      69             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: C_re, C_im
      70             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: filter_eps
      71             : 
      72             :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_complex_dbcsr_gemm_3'
      73             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
      74             : 
      75             :       CHARACTER(LEN=1)                                   :: transa2, transb2
      76             :       INTEGER                                            :: handle
      77             :       REAL(KIND=dp)                                      :: alpha2, alpha3, alpha4
      78             :       TYPE(dbcsr_type), POINTER                          :: a_plus_b, ac, bd, c_plus_d
      79             : 
      80        3946 :       CALL timeset(routineN, handle)
      81             :       !A complex matrix matrix multiplication can be done with only three multiplications
      82             :       !(a+ib)*(c+id)=ac-bd+i((a+b)*(c+d) - ac - bd)
      83             :       !A_re=a, A_im=b, B_re=c, B_im=d
      84             : 
      85        3946 :       alpha2 = -alpha
      86        3946 :       alpha3 = alpha
      87        3946 :       alpha4 = alpha
      88             : 
      89        3946 :       IF (transa == "C") THEN
      90           0 :          alpha2 = -alpha2
      91           0 :          alpha3 = -alpha3
      92           0 :          transa2 = "T"
      93             :       ELSE
      94        3946 :          transa2 = transa
      95             :       END IF
      96        3946 :       IF (transb == "C") THEN
      97        1072 :          alpha2 = -alpha2
      98        1072 :          alpha4 = -alpha4
      99        1072 :          transb2 = "T"
     100             :       ELSE
     101        2874 :          transb2 = transb
     102             :       END IF
     103             : 
     104             :       !create the work matrices
     105             :       NULLIFY (ac)
     106        3946 :       ALLOCATE (ac)
     107        3946 :       CALL dbcsr_create(ac, template=A_re, matrix_type="N")
     108             :       NULLIFY (bd)
     109        3946 :       ALLOCATE (bd)
     110        3946 :       CALL dbcsr_create(bd, template=A_re, matrix_type="N")
     111             :       NULLIFY (a_plus_b)
     112        3946 :       ALLOCATE (a_plus_b)
     113        3946 :       CALL dbcsr_create(a_plus_b, template=A_re, matrix_type="N")
     114             :       NULLIFY (c_plus_d)
     115        3946 :       ALLOCATE (c_plus_d)
     116        3946 :       CALL dbcsr_create(c_plus_d, template=A_re, matrix_type="N")
     117             : 
     118             :       !Do the neccesarry multiplications
     119        3946 :       CALL dbcsr_multiply(transa2, transb2, alpha, A_re, B_re, zero, ac, filter_eps=filter_eps)
     120        3946 :       CALL dbcsr_multiply(transa2, transb2, alpha2, A_im, B_im, zero, bd, filter_eps=filter_eps)
     121             : 
     122        3946 :       CALL dbcsr_add(a_plus_b, A_re, zero, alpha)
     123        3946 :       CALL dbcsr_add(a_plus_b, A_im, one, alpha3)
     124        3946 :       CALL dbcsr_add(c_plus_d, B_re, zero, alpha)
     125        3946 :       CALL dbcsr_add(c_plus_d, B_im, one, alpha4)
     126             : 
     127             :       !this can already be written into C_im
     128             :       !now both matrixes have been scaled which means we currently multiplied by alpha squared
     129        3946 :       CALL dbcsr_multiply(transa2, transb2, one/alpha, a_plus_b, c_plus_d, beta, C_im, filter_eps=filter_eps)
     130             : 
     131             :       !now add up all the terms into the result
     132        3946 :       CALL dbcsr_add(C_re, ac, beta, one)
     133             :       !the minus sign was already taken care of at the definition of alpha2
     134        3946 :       CALL dbcsr_add(C_re, bd, one, one)
     135        3946 :       CALL dbcsr_filter(C_re, filter_eps)
     136             : 
     137        3946 :       CALL dbcsr_add(C_im, ac, one, -one)
     138             :       !the minus sign was already taken care of at the definition of alpha2
     139        3946 :       CALL dbcsr_add(C_im, bd, one, one)
     140        3946 :       CALL dbcsr_filter(C_im, filter_eps)
     141             : 
     142             :       !Deallocate the work matrices
     143        3946 :       CALL dbcsr_deallocate_matrix(ac)
     144        3946 :       CALL dbcsr_deallocate_matrix(bd)
     145        3946 :       CALL dbcsr_deallocate_matrix(a_plus_b)
     146        3946 :       CALL dbcsr_deallocate_matrix(c_plus_d)
     147             : 
     148        3946 :       CALL timestop(handle)
     149             : 
     150        3946 :    END SUBROUTINE
     151             : 
     152             : ! **************************************************************************************************
     153             : !> \brief specialized subroutine for purely imaginary matrix exponentials
     154             : !> \param exp_H ...
     155             : !> \param im_matrix ...
     156             : !> \param nsquare ...
     157             : !> \param ntaylor ...
     158             : !> \param filter_eps ...
     159             : !> \author Samuel Andermatt (02.2014)
     160             : ! **************************************************************************************************
     161             : 
     162         708 :    SUBROUTINE taylor_only_imaginary_dbcsr(exp_H, im_matrix, nsquare, ntaylor, filter_eps)
     163             : 
     164             :       TYPE(dbcsr_p_type), DIMENSION(2)                   :: exp_H
     165             :       TYPE(dbcsr_type), POINTER                          :: im_matrix
     166             :       INTEGER, INTENT(in)                                :: nsquare, ntaylor
     167             :       REAL(KIND=dp), INTENT(in)                          :: filter_eps
     168             : 
     169             :       CHARACTER(len=*), PARAMETER :: routineN = 'taylor_only_imaginary_dbcsr'
     170             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     171             : 
     172             :       INTEGER                                            :: handle, i, nloop
     173             :       REAL(KIND=dp)                                      :: square_fac, Tfac, tmp
     174             :       TYPE(dbcsr_type), POINTER                          :: T1, T2, Tres_im, Tres_re
     175             : 
     176         708 :       CALL timeset(routineN, handle)
     177             : 
     178             :       !The divider that comes from the scaling and squaring procedure
     179         708 :       square_fac = 1.0_dp/(2.0_dp**REAL(nsquare, dp))
     180             : 
     181             :       !Allocate work matrices
     182             :       NULLIFY (T1)
     183         708 :       ALLOCATE (T1)
     184         708 :       CALL dbcsr_create(T1, template=im_matrix, matrix_type="N")
     185             :       NULLIFY (T2)
     186         708 :       ALLOCATE (T2)
     187         708 :       CALL dbcsr_create(T2, template=im_matrix, matrix_type="N")
     188             :       NULLIFY (Tres_re)
     189         708 :       ALLOCATE (Tres_re)
     190         708 :       CALL dbcsr_create(Tres_re, template=im_matrix, matrix_type="N")
     191             :       NULLIFY (Tres_im)
     192         708 :       ALLOCATE (Tres_im)
     193         708 :       CALL dbcsr_create(Tres_im, template=im_matrix, matrix_type="N")
     194             : 
     195             :       !Create the unit matrices
     196         708 :       CALL dbcsr_set(T1, zero)
     197         708 :       CALL dbcsr_add_on_diag(T1, one)
     198         708 :       CALL dbcsr_set(Tres_re, zero)
     199         708 :       CALL dbcsr_add_on_diag(Tres_re, one)
     200         708 :       CALL dbcsr_set(Tres_im, zero)
     201             : 
     202         708 :       nloop = CEILING(REAL(ntaylor, dp)/2.0_dp)
     203             :       !the inverse of the prefactor in the taylor series
     204         708 :       tmp = 1.0_dp
     205        3516 :       DO i = 1, nloop
     206        2808 :          CALL dbcsr_scale(T1, 1.0_dp/(REAL(i, dp)*2.0_dp - 1.0_dp))
     207        2808 :          CALL dbcsr_filter(T1, filter_eps)
     208             :          CALL dbcsr_multiply("N", "N", square_fac, im_matrix, T1, zero, &
     209        2808 :                              T2, filter_eps=filter_eps)
     210        2808 :          Tfac = one
     211        2808 :          IF (MOD(i, 2) == 0) Tfac = -Tfac
     212        2808 :          CALL dbcsr_add(Tres_im, T2, one, Tfac)
     213        2808 :          CALL dbcsr_scale(T2, 1.0_dp/(REAL(i, dp)*2.0_dp))
     214        2808 :          CALL dbcsr_filter(T2, filter_eps)
     215             :          CALL dbcsr_multiply("N", "N", square_fac, im_matrix, T2, zero, &
     216        2808 :                              T1, filter_eps=filter_eps)
     217        2808 :          Tfac = one
     218        2808 :          IF (MOD(i, 2) == 1) Tfac = -Tfac
     219        3516 :          CALL dbcsr_add(Tres_re, T1, one, Tfac)
     220             :       END DO
     221             : 
     222             :       !Square the matrices, due to the scaling and squaring procedure
     223         708 :       IF (nsquare .GT. 0) THEN
     224           0 :          DO i = 1, nsquare
     225             :             CALL cp_complex_dbcsr_gemm_3("N", "N", one, Tres_re, Tres_im, &
     226             :                                          Tres_re, Tres_im, zero, exp_H(1)%matrix, exp_H(2)%matrix, &
     227           0 :                                          filter_eps=filter_eps)
     228           0 :             CALL dbcsr_copy(Tres_re, exp_H(1)%matrix)
     229           0 :             CALL dbcsr_copy(Tres_im, exp_H(2)%matrix)
     230             :          END DO
     231             :       ELSE
     232         708 :          CALL dbcsr_copy(exp_H(1)%matrix, Tres_re)
     233         708 :          CALL dbcsr_copy(exp_H(2)%matrix, Tres_im)
     234             :       END IF
     235         708 :       CALL dbcsr_deallocate_matrix(T1)
     236         708 :       CALL dbcsr_deallocate_matrix(T2)
     237         708 :       CALL dbcsr_deallocate_matrix(Tres_re)
     238         708 :       CALL dbcsr_deallocate_matrix(Tres_im)
     239             : 
     240         708 :       CALL timestop(handle)
     241             : 
     242         708 :    END SUBROUTINE taylor_only_imaginary_dbcsr
     243             : 
     244             : ! **************************************************************************************************
     245             : !> \brief subroutine for general complex matrix exponentials
     246             : !>        on input a separate dbcsr_type for real and complex part
     247             : !>        on output a size 2 dbcsr_p_type, first element is the real part of
     248             : !>        the exponential second the imaginary
     249             : !> \param exp_H ...
     250             : !> \param re_part ...
     251             : !> \param im_part ...
     252             : !> \param nsquare ...
     253             : !> \param ntaylor ...
     254             : !> \param filter_eps ...
     255             : !> \author Samuel Andermatt (02.2014)
     256             : ! **************************************************************************************************
     257         264 :    SUBROUTINE taylor_full_complex_dbcsr(exp_H, re_part, im_part, nsquare, ntaylor, filter_eps)
     258             :       TYPE(dbcsr_p_type), DIMENSION(2)                   :: exp_H
     259             :       TYPE(dbcsr_type), POINTER                          :: re_part, im_part
     260             :       INTEGER, INTENT(in)                                :: nsquare, ntaylor
     261             :       REAL(KIND=dp), INTENT(in)                          :: filter_eps
     262             : 
     263             :       CHARACTER(len=*), PARAMETER :: routineN = 'taylor_full_complex_dbcsr'
     264             :       COMPLEX(KIND=dp), PARAMETER                        :: one = (1.0_dp, 0.0_dp), &
     265             :                                                             zero = (0.0_dp, 0.0_dp)
     266             : 
     267             :       COMPLEX(KIND=dp)                                   :: square_fac
     268             :       INTEGER                                            :: handle, i
     269             :       TYPE(dbcsr_type), POINTER                          :: T1, T2, T3, Tres
     270             : 
     271         264 :       CALL timeset(routineN, handle)
     272             : 
     273             :       !The divider that comes from the scaling and squaring procedure
     274         264 :       square_fac = CMPLX(1.0_dp/(2.0_dp**REAL(nsquare, dp)), 0.0_dp, KIND=dp)
     275             : 
     276             :       !Allocate work matrices
     277             :       NULLIFY (T1)
     278         264 :       ALLOCATE (T1)
     279             :       CALL dbcsr_create(T1, template=re_part, matrix_type="N", &
     280         264 :                         data_type=dbcsr_type_complex_8)
     281             :       NULLIFY (T2)
     282         264 :       ALLOCATE (T2)
     283             :       CALL dbcsr_create(T2, template=re_part, matrix_type="N", &
     284         264 :                         data_type=dbcsr_type_complex_8)
     285             :       NULLIFY (T3)
     286         264 :       ALLOCATE (T3)
     287             :       CALL dbcsr_create(T3, template=re_part, matrix_type="N", &
     288         264 :                         data_type=dbcsr_type_complex_8)
     289             :       NULLIFY (Tres)
     290         264 :       ALLOCATE (Tres)
     291             :       CALL dbcsr_create(Tres, template=re_part, matrix_type="N", &
     292         264 :                         data_type=dbcsr_type_complex_8)
     293             : 
     294             :       !Fuse the input matrices to a single complex matrix
     295         264 :       CALL dbcsr_copy(T3, re_part)
     296         264 :       CALL dbcsr_copy(Tres, im_part) !will later on be set back to zero
     297         264 :       CALL dbcsr_scale(Tres, CMPLX(0.0_dp, 1.0_dp, KIND=dp))
     298         264 :       CALL dbcsr_add(T3, Tres, one, one)
     299             : 
     300             :       !Create the unit matrices
     301         264 :       CALL dbcsr_set(T1, zero)
     302         264 :       CALL dbcsr_add_on_diag(T1, one)
     303         264 :       CALL dbcsr_set(Tres, zero)
     304         264 :       CALL dbcsr_add_on_diag(Tres, one)
     305             : 
     306        2522 :       DO i = 1, ntaylor
     307        2258 :          CALL dbcsr_scale(T1, one/CMPLX(i*1.0_dp, 0.0_dp, KIND=dp))
     308        2258 :          CALL dbcsr_filter(T1, filter_eps)
     309             :          CALL dbcsr_multiply("N", "N", square_fac, T1, T3, &
     310        2258 :                              zero, T2, filter_eps=filter_eps)
     311        2258 :          CALL dbcsr_add(Tres, T2, one, one)
     312        2522 :          CALL dbcsr_copy(T1, T2)
     313             :       END DO
     314             : 
     315         264 :       IF (nsquare .GT. 0) THEN
     316         970 :          DO i = 1, nsquare
     317             :             CALL dbcsr_multiply("N", "N", one, Tres, Tres, zero, &
     318         758 :                                 T2, filter_eps=filter_eps)
     319         970 :             CALL dbcsr_copy(Tres, T2)
     320             :          END DO
     321             :       END IF
     322             : 
     323         264 :       CALL dbcsr_copy(exp_H(1)%matrix, Tres, keep_imaginary=.FALSE.)
     324         264 :       CALL dbcsr_scale(Tres, CMPLX(0.0_dp, -1.0_dp, KIND=dp))
     325         264 :       CALL dbcsr_copy(exp_H(2)%matrix, Tres, keep_imaginary=.FALSE.)
     326             : 
     327         264 :       CALL dbcsr_deallocate_matrix(T1)
     328         264 :       CALL dbcsr_deallocate_matrix(T2)
     329         264 :       CALL dbcsr_deallocate_matrix(T3)
     330         264 :       CALL dbcsr_deallocate_matrix(Tres)
     331             : 
     332         264 :       CALL timestop(handle)
     333             : 
     334         264 :    END SUBROUTINE taylor_full_complex_dbcsr
     335             : 
     336             : ! **************************************************************************************************
     337             : !> \brief  The Baker-Campbell-Hausdorff expansion for a purely imaginary exponent (e.g. rtp)
     338             : !>         Works for a non unitary propagator, because the density matrix is hermitian
     339             : !> \param propagator The exponent of the matrix exponential
     340             : !> \param density_re Real part of the density matrix
     341             : !> \param density_im Imaginary part of the density matrix
     342             : !> \param filter_eps The filtering threshold for all matrices
     343             : !> \param filter_eps_small ...
     344             : !> \param eps_exp The accuracy of the exponential
     345             : !> \author Samuel Andermatt (02.2014)
     346             : ! **************************************************************************************************
     347             : 
     348         112 :    SUBROUTINE bch_expansion_imaginary_propagator(propagator, density_re, density_im, filter_eps, filter_eps_small, eps_exp)
     349             :       TYPE(dbcsr_type), POINTER                          :: propagator, density_re, density_im
     350             :       REAL(KIND=dp), INTENT(in)                          :: filter_eps, filter_eps_small, eps_exp
     351             : 
     352             :       CHARACTER(len=*), PARAMETER :: routineN = 'bch_expansion_imaginary_propagator'
     353             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     354             : 
     355             :       INTEGER                                            :: handle, i, unit_nr
     356             :       LOGICAL                                            :: convergence
     357             :       REAL(KIND=dp)                                      :: alpha, max_alpha, prefac
     358             :       TYPE(dbcsr_type), POINTER                          :: comm, comm2, tmp, tmp2
     359             : 
     360         112 :       CALL timeset(routineN, handle)
     361             : 
     362         112 :       unit_nr = cp_logger_get_default_io_unit()
     363             : 
     364             :       NULLIFY (tmp)
     365         112 :       ALLOCATE (tmp)
     366         112 :       CALL dbcsr_create(tmp, template=propagator)
     367             :       NULLIFY (tmp2)
     368         112 :       ALLOCATE (tmp2)
     369         112 :       CALL dbcsr_create(tmp2, template=propagator)
     370             :       NULLIFY (comm)
     371         112 :       ALLOCATE (comm)
     372         112 :       CALL dbcsr_create(comm, template=propagator)
     373             :       NULLIFY (comm2)
     374         112 :       ALLOCATE (comm2)
     375         112 :       CALL dbcsr_create(comm2, template=propagator)
     376             : 
     377         112 :       CALL dbcsr_copy(tmp, density_re)
     378         112 :       CALL dbcsr_copy(tmp2, density_im)
     379             : 
     380         112 :       convergence = .FALSE.
     381         494 :       DO i = 1, 20
     382         494 :          prefac = one/i
     383             :          CALL dbcsr_multiply("N", "N", -prefac, propagator, tmp2, zero, comm, &
     384         494 :                              filter_eps=filter_eps_small)
     385             :          CALL dbcsr_multiply("N", "N", prefac, propagator, tmp, zero, comm2, &
     386         494 :                              filter_eps=filter_eps_small)
     387         494 :          CALL dbcsr_transposed(tmp, comm)
     388         494 :          CALL dbcsr_transposed(tmp2, comm2)
     389         494 :          CALL dbcsr_add(comm, tmp, one, one)
     390         494 :          CALL dbcsr_add(comm2, tmp2, one, -one)
     391         494 :          CALL dbcsr_add(density_re, comm, one, one)
     392         494 :          CALL dbcsr_add(density_im, comm2, one, one)
     393         494 :          CALL dbcsr_copy(tmp, comm)
     394         494 :          CALL dbcsr_copy(tmp2, comm2)
     395             :          !check for convergence
     396         494 :          max_alpha = zero
     397         494 :          alpha = dbcsr_frobenius_norm(comm)
     398         494 :          IF (alpha > max_alpha) max_alpha = alpha
     399         494 :          alpha = dbcsr_frobenius_norm(comm2)
     400         494 :          IF (alpha > max_alpha) max_alpha = alpha
     401         494 :          IF (max_alpha < eps_exp) convergence = .TRUE.
     402         382 :          IF (convergence) THEN
     403         112 :             IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="((T3,A,I2,A))") &
     404          56 :                "BCH converged after ", i, " steps"
     405             :             EXIT
     406             :          END IF
     407             :       END DO
     408             : 
     409         112 :       CALL dbcsr_filter(density_re, filter_eps)
     410         112 :       CALL dbcsr_filter(density_im, filter_eps)
     411             : 
     412         112 :       IF (.NOT. convergence) &
     413           0 :          CPWARN("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
     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             : 
     437         130 :    SUBROUTINE bch_expansion_complex_propagator(propagator_re, propagator_im, density_re, density_im, filter_eps, &
     438             :                                                filter_eps_small, eps_exp)
     439             :       TYPE(dbcsr_type), POINTER                          :: propagator_re, propagator_im, &
     440             :                                                             density_re, density_im
     441             :       REAL(KIND=dp), INTENT(in)                          :: filter_eps, filter_eps_small, eps_exp
     442             : 
     443             :       CHARACTER(len=*), PARAMETER :: routineN = 'bch_expansion_complex_propagator'
     444             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     445             : 
     446             :       INTEGER                                            :: handle, i, unit_nr
     447             :       LOGICAL                                            :: convergence
     448             :       REAL(KIND=dp)                                      :: alpha, max_alpha, prefac
     449             :       TYPE(dbcsr_type), POINTER                          :: comm, comm2, tmp, tmp2
     450             : 
     451         130 :       CALL timeset(routineN, handle)
     452             : 
     453         130 :       unit_nr = cp_logger_get_default_io_unit()
     454             : 
     455             :       NULLIFY (tmp)
     456         130 :       ALLOCATE (tmp)
     457         130 :       CALL dbcsr_create(tmp, template=propagator_re)
     458             :       NULLIFY (tmp2)
     459         130 :       ALLOCATE (tmp2)
     460         130 :       CALL dbcsr_create(tmp2, template=propagator_re)
     461             :       NULLIFY (comm)
     462         130 :       ALLOCATE (comm)
     463         130 :       CALL dbcsr_create(comm, template=propagator_re)
     464             :       NULLIFY (comm2)
     465         130 :       ALLOCATE (comm2)
     466         130 :       CALL dbcsr_create(comm2, template=propagator_re)
     467             : 
     468         130 :       CALL dbcsr_copy(tmp, density_re)
     469         130 :       CALL dbcsr_copy(tmp2, density_im)
     470             : 
     471         130 :       convergence = .FALSE.
     472             : 
     473        1254 :       DO i = 1, 20
     474        1254 :          prefac = one/i
     475             :          CALL cp_complex_dbcsr_gemm_3("N", "N", prefac, propagator_re, propagator_im, &
     476        1254 :                                       tmp, tmp2, zero, comm, comm2, filter_eps=filter_eps_small)
     477        1254 :          CALL dbcsr_transposed(tmp, comm)
     478        1254 :          CALL dbcsr_transposed(tmp2, comm2)
     479        1254 :          CALL dbcsr_add(comm, tmp, one, one)
     480        1254 :          CALL dbcsr_add(comm2, tmp2, one, -one)
     481        1254 :          CALL dbcsr_add(density_re, comm, one, one)
     482        1254 :          CALL dbcsr_add(density_im, comm2, one, one)
     483        1254 :          CALL dbcsr_copy(tmp, comm)
     484        1254 :          CALL dbcsr_copy(tmp2, comm2)
     485             :          !check for convergence
     486        1254 :          max_alpha = zero
     487        1254 :          alpha = dbcsr_frobenius_norm(comm)
     488        1254 :          IF (alpha > max_alpha) max_alpha = alpha
     489        1254 :          alpha = dbcsr_frobenius_norm(comm2)
     490        1254 :          IF (alpha > max_alpha) max_alpha = alpha
     491        1254 :          IF (max_alpha < eps_exp) convergence = .TRUE.
     492        1124 :          IF (convergence) THEN
     493         130 :             IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="((T3,A,I2,A))") &
     494          65 :                "BCH converged after ", i, " steps"
     495             :             EXIT
     496             :          END IF
     497             :       END DO
     498             : 
     499         130 :       CALL dbcsr_filter(density_re, filter_eps)
     500         130 :       CALL dbcsr_filter(density_im, filter_eps)
     501             : 
     502         130 :       IF (.NOT. convergence) &
     503           0 :          CPWARN("BCH method did not converge ")
     504             : 
     505         130 :       CALL dbcsr_deallocate_matrix(tmp)
     506         130 :       CALL dbcsr_deallocate_matrix(tmp2)
     507         130 :       CALL dbcsr_deallocate_matrix(comm)
     508         130 :       CALL dbcsr_deallocate_matrix(comm2)
     509             : 
     510         130 :       CALL timestop(handle)
     511             : 
     512         130 :    END SUBROUTINE
     513             : 
     514             : END MODULE ls_matrix_exp

Generated by: LCOV version 1.15