LCOV - code coverage report
Current view: top level - src - bse_main.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b1f098b) Lines: 42 47 89.4 %
Date: 2024-05-05 06:30:09 Functions: 1 1 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 Main routines for GW + Bethe-Salpeter for computing electronic excitations
      10             : !> \par History
      11             : !>      04.2024 created [Maximilian Graml]
      12             : ! **************************************************************************************************
      13             : 
      14             : MODULE bse_main
      15             : 
      16             :    USE bse_full_diag,                   ONLY: create_A,&
      17             :                                               create_B,&
      18             :                                               create_hermitian_form_of_ABBA,&
      19             :                                               diagonalize_A,&
      20             :                                               diagonalize_C
      21             :    USE bse_iterative,                   ONLY: do_subspace_iterations,&
      22             :                                               fill_local_3c_arrays
      23             :    USE bse_util,                        ONLY: deallocate_matrices_bse,&
      24             :                                               estimate_BSE_resources,&
      25             :                                               mult_B_with_W,&
      26             :                                               print_BSE_start_flag,&
      27             :                                               truncate_BSE_matrices
      28             :    USE cp_fm_types,                     ONLY: cp_fm_release,&
      29             :                                               cp_fm_type
      30             :    USE group_dist_types,                ONLY: group_dist_d1_type
      31             :    USE input_constants,                 ONLY: bse_abba,&
      32             :                                               bse_both,&
      33             :                                               bse_fulldiag,&
      34             :                                               bse_iterdiag,&
      35             :                                               bse_tda
      36             :    USE kinds,                           ONLY: dp
      37             :    USE message_passing,                 ONLY: mp_para_env_type
      38             :    USE mp2_types,                       ONLY: mp2_type
      39             : #include "./base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main'
      46             : 
      47             :    PUBLIC :: start_bse_calculation
      48             : 
      49             : CONTAINS
      50             : 
      51             : ! **************************************************************************************************
      52             : !> \brief Main subroutine managing BSE calculations
      53             : !> \param fm_mat_S_ia_bse ...
      54             : !> \param fm_mat_S_ij_bse ...
      55             : !> \param fm_mat_S_ab_bse ...
      56             : !> \param fm_mat_Q_static_bse ...
      57             : !> \param fm_mat_Q_static_bse_gemm ...
      58             : !> \param Eigenval ...
      59             : !> \param homo ...
      60             : !> \param virtual ...
      61             : !> \param dimen_RI ...
      62             : !> \param dimen_RI_red ...
      63             : !> \param gd_array ...
      64             : !> \param color_sub ...
      65             : !> \param mp2_env ...
      66             : !> \param unit_nr ...
      67             : ! **************************************************************************************************
      68           4 :    SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
      69             :                                     fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm, &
      70           4 :                                     Eigenval, homo, virtual, dimen_RI, dimen_RI_red, &
      71             :                                     gd_array, color_sub, mp2_env, unit_nr)
      72             : 
      73             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, &
      74             :                                                             fm_mat_S_ab_bse
      75             :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_mat_Q_static_bse, &
      76             :                                                             fm_mat_Q_static_bse_gemm
      77             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
      78             :          INTENT(IN)                                      :: Eigenval
      79             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo, virtual
      80             :       INTEGER, INTENT(IN)                                :: dimen_RI, dimen_RI_red
      81             :       TYPE(group_dist_d1_type), INTENT(IN)               :: gd_array
      82             :       INTEGER, INTENT(IN)                                :: color_sub
      83             :       TYPE(mp2_type)                                     :: mp2_env
      84             :       INTEGER, INTENT(IN)                                :: unit_nr
      85             : 
      86             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'start_bse_calculation'
      87             : 
      88             :       INTEGER                                            :: handle, homo_red, virtual_red
      89             :       LOGICAL                                            :: my_do_abba, my_do_fulldiag, &
      90             :                                                             my_do_iterat_diag, my_do_tda
      91             :       REAL(KIND=dp)                                      :: diag_runtime_est
      92           4 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: Eigenval_reduced
      93           4 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: B_abQ_bse_local, B_bar_iaQ_bse_local, &
      94           4 :                                                             B_bar_ijQ_bse_local, B_iaQ_bse_local
      95             :       TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_C_BSE, fm_inv_sqrt_A_minus_B, fm_mat_S_ab_trunc, &
      96             :          fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, &
      97             :          fm_sqrt_A_minus_B
      98             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      99             : 
     100           4 :       CALL timeset(routineN, handle)
     101             : 
     102           4 :       para_env => fm_mat_S_ia_bse%matrix_struct%para_env
     103             : 
     104           4 :       my_do_fulldiag = .FALSE.
     105           4 :       my_do_iterat_diag = .FALSE.
     106           4 :       my_do_tda = .FALSE.
     107           4 :       my_do_abba = .FALSE.
     108             :       !Method: Iterative or full diagonalization
     109           4 :       SELECT CASE (mp2_env%ri_g0w0%bse_diag_method)
     110             :       CASE (bse_iterdiag)
     111           0 :          my_do_iterat_diag = .TRUE.
     112             :          !MG: Basics of the Davidson solver are implemented, but not rigorously checked.
     113           0 :          CPABORT("Iterative BSE not yet implemented")
     114             :       CASE (bse_fulldiag)
     115           4 :          my_do_fulldiag = .TRUE.
     116             :       END SELECT
     117             :       !Approximation: TDA and/or full ABBA matrix
     118           6 :       SELECT CASE (mp2_env%ri_g0w0%bse_approx)
     119             :       CASE (bse_tda)
     120           2 :          my_do_tda = .TRUE.
     121             :       CASE (bse_abba)
     122           2 :          my_do_abba = .TRUE.
     123             :       CASE (bse_both)
     124           0 :          my_do_tda = .TRUE.
     125           4 :          my_do_abba = .TRUE.
     126             :       END SELECT
     127             : 
     128           4 :       CALL print_BSE_start_flag(my_do_tda, my_do_abba, unit_nr)
     129             : 
     130           4 :       CALL fm_mat_S_ia_bse%matrix_struct%para_env%sync()
     131             : 
     132             :       ! Reduce matrices in case of energy cutoff for occupied and unoccupied in A/B-BSE-matrices
     133             :       CALL truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
     134             :                                  fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
     135             :                                  Eigenval(:, 1, 1), Eigenval_reduced, &
     136             :                                  homo(1), virtual(1), dimen_RI, unit_nr, &
     137             :                                  homo_red, virtual_red, &
     138           4 :                                  mp2_env)
     139             :       ! \bar{B}^P_rs = \sum_R W_PR B^R_rs where B^R_rs = \sum_T [1/sqrt(v)]_RT (T|rs)
     140             :       ! r,s: MO-index, P,R,T: RI-index
     141             :       ! B: fm_mat_S_..., W: fm_mat_Q_...
     142             :       CALL mult_B_with_W(fm_mat_S_ij_trunc, fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, &
     143             :                          fm_mat_S_bar_ij_bse, fm_mat_Q_static_bse_gemm, &
     144           4 :                          dimen_RI_red, homo_red, virtual_red)
     145             : 
     146           4 :       IF (my_do_iterat_diag) THEN
     147             :          CALL fill_local_3c_arrays(fm_mat_S_ab_trunc, fm_mat_S_ia_trunc, &
     148             :                                    fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
     149             :                                    B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
     150             :                                    B_iaQ_bse_local, dimen_RI_red, homo_red, virtual_red, &
     151             :                                    gd_array, color_sub, para_env)
     152             :       END IF
     153             : 
     154           4 :       IF (my_do_fulldiag) THEN
     155             :          ! Quick estimate of memory consumption and runtime of diagonalizations
     156             :          CALL estimate_BSE_resources(homo_red, virtual_red, unit_nr, my_do_abba, &
     157           4 :                                      para_env, diag_runtime_est)
     158             :          ! Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
     159             :          ! A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
     160             :          ! ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
     161             :          ! α is a spin-dependent factor
     162             :          ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
     163             :          ! W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
     164             :          CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_bar_ij_bse, fm_mat_S_ab_trunc, &
     165             :                        fm_A_BSE, Eigenval_reduced, unit_nr, &
     166             :                        homo_red, virtual_red, dimen_RI, mp2_env, &
     167           4 :                        para_env)
     168           4 :          IF (my_do_abba) THEN
     169             :             ! Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
     170             :             ! B_ia,jb = α * v_ia,jb - W_ib,aj
     171             :             ! α is a spin-dependent factor
     172             :             ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
     173             :             ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
     174             :             CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, fm_B_BSE, &
     175           2 :                           homo_red, virtual_red, dimen_RI, unit_nr, mp2_env)
     176             :             ! Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
     177             :             ! (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
     178             :             ! We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
     179             :             ! of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
     180             :             CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, &
     181             :                                                fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
     182           2 :                                                homo_red, virtual_red, unit_nr, mp2_env, diag_runtime_est)
     183             :          END IF
     184           4 :          CALL cp_fm_release(fm_B_BSE)
     185           4 :          IF (my_do_tda) THEN
     186             :             ! Solving the hermitian eigenvalue equation A X^n = Ω^n X^n
     187             :             CALL diagonalize_A(fm_A_BSE, homo_red, virtual_red, homo(1), &
     188           2 :                                unit_nr, diag_runtime_est, mp2_env)
     189             :          END IF
     190             :          ! Release to avoid faulty use of changed A matrix
     191           4 :          CALL cp_fm_release(fm_A_BSE)
     192           4 :          IF (my_do_abba) THEN
     193             :             ! Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
     194             :             ! Here, the eigenvectors Z^n relate to X^n via
     195             :             ! Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
     196             :             CALL diagonalize_C(fm_C_BSE, homo_red, virtual_red, homo(1), &
     197             :                                fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
     198           2 :                                unit_nr, diag_runtime_est, mp2_env)
     199             :          END IF
     200             :          ! Release to avoid faulty use of changed C matrix
     201           4 :          CALL cp_fm_release(fm_C_BSE)
     202             :       END IF
     203             : 
     204             :       CALL deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
     205             :                                    fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
     206           4 :                                    fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm)
     207           4 :       DEALLOCATE (Eigenval_reduced)
     208           4 :       IF (my_do_iterat_diag) THEN
     209             :          ! Contains untested Block-Davidson algorithm
     210             :          CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
     211             :                                      B_iaQ_bse_local, homo(1), virtual(1), mp2_env%ri_g0w0%bse_spin_config, unit_nr, &
     212           0 :                                      Eigenval(:, 1, 1), para_env, mp2_env)
     213             :          ! Deallocate local 3c-B-matrices
     214           0 :          DEALLOCATE (B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, B_iaQ_bse_local)
     215             :       END IF
     216             : 
     217           4 :       IF (unit_nr > 0) THEN
     218           2 :          WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!'
     219             :       END IF
     220             : 
     221           4 :       CALL timestop(handle)
     222             : 
     223           8 :    END SUBROUTINE start_bse_calculation
     224             : 
     225             : END MODULE bse_main

Generated by: LCOV version 1.15