bse_main.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Main routines for GW + Bethe-Salpeter for computing electronic excitations
!> \par History
!>      04.2024 created [Maximilian Graml]
! **************************************************************************************************
 
MODULE bse_main
 
   USE bse_full_diag,                   ONLY: create_a,&
                                              create_b,&
                                              create_hermitian_form_of_abba,&
                                              diagonalize_a,&
                                              diagonalize_c
   USE bse_iterative,                   ONLY: do_subspace_iterations,&
                                              fill_local_3c_arrays
   USE bse_print,                       ONLY: print_bse_start_flag
   USE bse_util,                        ONLY: adapt_bse_input_params,&
                                              deallocate_matrices_bse,&
                                              estimate_bse_resources,&
                                              mult_b_with_w,&
                                              truncate_bse_matrices
   USE cp_control_types,                ONLY: dft_control_type,&
                                              tddfpt2_control_type
   USE cp_fm_types,                     ONLY: cp_fm_release,&
                                              cp_fm_type
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: debug_print_level
   USE group_dist_types,                ONLY: group_dist_d1_type
   USE input_constants,                 ONLY: bse_abba,&
                                              bse_both,&
                                              bse_fulldiag,&
                                              bse_iterdiag,&
                                              bse_screening_alpha,&
                                              bse_screening_tdhf,&
                                              bse_tda
   USE kinds,                           ONLY: dp
   USE message_passing,                 ONLY: mp_para_env_type
   USE mp2_types,                       ONLY: mp2_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main'
 
   PUBLIC :: start_bse_calculation
 
CONTAINS
 
! **************************************************************************************************
!> \brief Main subroutine managing BSE calculations
!> \param fm_mat_S_ia_bse ...
!> \param fm_mat_S_ij_bse ...
!> \param fm_mat_S_ab_bse ...
!> \param fm_mat_Q_static_bse_gemm ...
!> \param Eigenval ...
!> \param Eigenval_scf ...
!> \param homo ...
!> \param virtual ...
!> \param dimen_RI ...
!> \param dimen_RI_red ...
!> \param bse_lev_virt ...
!> \param gd_array ...
!> \param color_sub ...
!> \param mp2_env ...
!> \param qs_env ...
!> \param mo_coeff ...
!> \param unit_nr ...
! **************************************************************************************************
   SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
                                    fm_mat_Q_static_bse_gemm, &
                                    Eigenval, Eigenval_scf, &
                                    homo, virtual, dimen_RI, dimen_RI_red, bse_lev_virt, &
                                    gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
 
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_s_ia_bse, fm_mat_s_ij_bse, &
                                                            fm_mat_s_ab_bse
      TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_mat_q_static_bse_gemm
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
         INTENT(IN)                                      :: eigenval, eigenval_scf
      INTEGER, DIMENSION(:), INTENT(IN)                  :: homo, virtual
      INTEGER, INTENT(IN)                                :: dimen_ri, dimen_ri_red, bse_lev_virt
      TYPE(group_dist_d1_type), INTENT(IN)               :: gd_array
      INTEGER, INTENT(IN)                                :: color_sub
      TYPE(mp2_type)                                     :: mp2_env
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff
      INTEGER, INTENT(IN)                                :: unit_nr
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'start_bse_calculation'
 
      INTEGER                                            :: handle, homo_red, virtual_red
      LOGICAL                                            :: my_do_abba, my_do_fulldiag, &
                                                            my_do_iterat_diag, my_do_tda
      REAL(kind=dp)                                      :: diag_runtime_est
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: eigenval_reduced
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: b_abq_bse_local, b_bar_iaq_bse_local, &
                                                            b_bar_ijq_bse_local, b_iaq_bse_local
      TYPE(cp_fm_type) :: fm_a_bse, fm_b_bse, fm_c_bse, fm_inv_sqrt_a_minus_b, fm_mat_s_ab_trunc, &
         fm_mat_s_bar_ia_bse, fm_mat_s_bar_ij_bse, fm_mat_s_ia_trunc, fm_mat_s_ij_trunc, &
         fm_sqrt_a_minus_b
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
 
      CALL timeset(routinen, handle)
 
      para_env => fm_mat_s_ia_bse%matrix_struct%para_env
 
      my_do_fulldiag = .false.
      my_do_iterat_diag = .false.
      my_do_tda = .false.
      my_do_abba = .false.
      !Method: Iterative or full diagonalization
      SELECT CASE (mp2_env%bse%bse_diag_method)
      CASE (bse_iterdiag)
         my_do_iterat_diag = .true.
         !MG: Basics of the Davidson solver are implemented, but not rigorously checked.
         cpabort("Iterative BSE not yet implemented")
      CASE (bse_fulldiag)
         my_do_fulldiag = .true.
      END SELECT
      !Approximation: TDA and/or full ABBA matrix
      SELECT CASE (mp2_env%bse%flag_tda)
      CASE (bse_tda)
         my_do_tda = .true.
      CASE (bse_abba)
         my_do_abba = .true.
      CASE (bse_both)
         my_do_tda = .true.
         my_do_abba = .true.
      END SELECT
 
      CALL print_bse_start_flag(my_do_tda, my_do_abba, unit_nr)
 
      ! Link BSE debug flag against debug print level
      logger => cp_get_default_logger()
      IF (logger%iter_info%print_level == debug_print_level) THEN
         mp2_env%bse%bse_debug_print = .true.
      END IF
 
      CALL fm_mat_s_ia_bse%matrix_struct%para_env%sync()
      ! We apply the BSE cutoffs using the DFT Eigenenergies
      ! Reduce matrices in case of energy cutoff for occupied and unoccupied in A/B-BSE-matrices
      CALL truncate_bse_matrices(fm_mat_s_ia_bse, fm_mat_s_ij_bse, fm_mat_s_ab_bse, &
                                 fm_mat_s_ia_trunc, fm_mat_s_ij_trunc, fm_mat_s_ab_trunc, &
                                 eigenval_scf(:, 1, 1), eigenval(:, 1, 1), eigenval_reduced, &
                                 homo(1), virtual(1), dimen_ri, unit_nr, &
                                 bse_lev_virt, &
                                 homo_red, virtual_red, &
                                 mp2_env)
      ! \bar{B}^P_rs = \sum_R W_PR B^R_rs where B^R_rs = \sum_T [1/sqrt(v)]_RT (T|rs)
      ! r,s: MO-index, P,R,T: RI-index
      ! B: fm_mat_S_..., W: fm_mat_Q_...
      CALL mult_b_with_w(fm_mat_s_ij_trunc, fm_mat_s_ia_trunc, fm_mat_s_bar_ia_bse, &
                         fm_mat_s_bar_ij_bse, fm_mat_q_static_bse_gemm, &
                         dimen_ri_red, homo_red, virtual_red)
 
      IF (my_do_iterat_diag) THEN
         CALL fill_local_3c_arrays(fm_mat_s_ab_trunc, fm_mat_s_ia_trunc, &
                                   fm_mat_s_bar_ia_bse, fm_mat_s_bar_ij_bse, &
                                   b_bar_ijq_bse_local, b_abq_bse_local, b_bar_iaq_bse_local, &
                                   b_iaq_bse_local, dimen_ri_red, homo_red, virtual_red, &
                                   gd_array, color_sub, para_env)
      END IF
 
      CALL adapt_bse_input_params(homo_red, virtual_red, unit_nr, mp2_env, qs_env)
 
      IF (my_do_fulldiag) THEN
         ! Quick estimate of memory consumption and runtime of diagonalizations
         CALL estimate_bse_resources(homo_red, virtual_red, unit_nr, my_do_abba, &
                                     para_env, diag_runtime_est)
         ! Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
         ! A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
         ! ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
         ! α is a spin-dependent factor
         ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
         ! W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
 
         ! For unscreened W matrix, we need fm_mat_S_ij_trunc
         IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
             mp2_env%bse%screening_method == bse_screening_alpha) THEN
            CALL create_a(fm_mat_s_ia_trunc, fm_mat_s_ij_trunc, fm_mat_s_ab_trunc, &
                          fm_a_bse, eigenval_reduced, unit_nr, &
                          homo_red, virtual_red, dimen_ri, mp2_env, &
                          para_env, qs_env)
         ELSE
            CALL create_a(fm_mat_s_ia_trunc, fm_mat_s_bar_ij_bse, fm_mat_s_ab_trunc, &
                          fm_a_bse, eigenval_reduced, unit_nr, &
                          homo_red, virtual_red, dimen_ri, mp2_env, &
                          para_env, qs_env)
         END IF
         IF (my_do_abba) THEN
            ! Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
            ! B_ia,jb = α * v_ia,jb - W_ib,aj
            ! α is a spin-dependent factor
            ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
            ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
 
            ! For unscreened W matrix, we need fm_mat_S_ia_trunc
            IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
                mp2_env%bse%screening_method == bse_screening_alpha) THEN
               CALL create_b(fm_mat_s_ia_trunc, fm_mat_s_ia_trunc, fm_b_bse, &
                             homo_red, virtual_red, dimen_ri, unit_nr, mp2_env)
            ELSE
               CALL create_b(fm_mat_s_ia_trunc, fm_mat_s_bar_ia_bse, fm_b_bse, &
                             homo_red, virtual_red, dimen_ri, unit_nr, mp2_env)
            END IF
            ! Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
            ! (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
            ! We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
            ! of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
            CALL create_hermitian_form_of_abba(fm_a_bse, fm_b_bse, fm_c_bse, &
                                               fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, &
                                               homo_red, virtual_red, unit_nr, mp2_env, diag_runtime_est)
         END IF
         CALL cp_fm_release(fm_b_bse)
 
         NULLIFY (dft_control, tddfpt_control)
         CALL get_qs_env(qs_env, dft_control=dft_control)
         tddfpt_control => dft_control%tddfpt2_control
         IF ((my_do_tda) .AND. (.NOT. tddfpt_control%do_bse)) THEN
            ! Solving the hermitian eigenvalue equation A X^n = Ω^n X^n
            CALL diagonalize_a(fm_a_bse, homo_red, virtual_red, homo(1), &
                               unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
         END IF
         ! Release to avoid faulty use of changed A matrix
         CALL cp_fm_release(fm_a_bse)
         IF (my_do_abba) THEN
            ! Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
            ! Here, the eigenvectors Z^n relate to X^n via
            ! Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
            CALL diagonalize_c(fm_c_bse, homo_red, virtual_red, homo(1), &
                               fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, &
                               unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
         END IF
         ! Release to avoid faulty use of changed C matrix
         CALL cp_fm_release(fm_c_bse)
      END IF
 
      CALL deallocate_matrices_bse(fm_mat_s_bar_ia_bse, fm_mat_s_bar_ij_bse, &
                                   fm_mat_s_ia_trunc, fm_mat_s_ij_trunc, fm_mat_s_ab_trunc, &
                                   fm_mat_q_static_bse_gemm, mp2_env)
      DEALLOCATE (eigenval_reduced)
      IF (my_do_iterat_diag) THEN
         ! Contains untested Block-Davidson algorithm
         CALL do_subspace_iterations(b_bar_ijq_bse_local, b_abq_bse_local, b_bar_iaq_bse_local, &
                                     b_iaq_bse_local, homo(1), virtual(1), mp2_env%bse%bse_spin_config, unit_nr, &
                                     eigenval(:, 1, 1), para_env, mp2_env)
         ! Deallocate local 3c-B-matrices
         DEALLOCATE (b_bar_ijq_bse_local, b_abq_bse_local, b_bar_iaq_bse_local, b_iaq_bse_local)
      END IF
 
      IF (unit_nr > 0) THEN
         WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!'
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE start_bse_calculation
 
END MODULE bse_main