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 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_print, ONLY: print_BSE_start_flag
24 : USE bse_util, ONLY: adapt_BSE_input_params,&
25 : deallocate_matrices_bse,&
26 : estimate_BSE_resources,&
27 : mult_B_with_W,&
28 : truncate_BSE_matrices
29 : USE cp_fm_types, ONLY: cp_fm_release,&
30 : cp_fm_type
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: debug_print_level
34 : USE group_dist_types, ONLY: group_dist_d1_type
35 : USE input_constants, ONLY: bse_abba,&
36 : bse_both,&
37 : bse_fulldiag,&
38 : bse_iterdiag,&
39 : bse_screening_alpha,&
40 : bse_screening_tdhf,&
41 : bse_tda
42 : USE kinds, ONLY: dp
43 : USE message_passing, ONLY: mp_para_env_type
44 : USE mp2_types, ONLY: mp2_type
45 : USE qs_environment_types, ONLY: qs_environment_type
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 :
50 : PRIVATE
51 :
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main'
53 :
54 : PUBLIC :: start_bse_calculation
55 :
56 : CONTAINS
57 :
58 : ! **************************************************************************************************
59 : !> \brief Main subroutine managing BSE calculations
60 : !> \param fm_mat_S_ia_bse ...
61 : !> \param fm_mat_S_ij_bse ...
62 : !> \param fm_mat_S_ab_bse ...
63 : !> \param fm_mat_Q_static_bse_gemm ...
64 : !> \param Eigenval ...
65 : !> \param Eigenval_scf ...
66 : !> \param homo ...
67 : !> \param virtual ...
68 : !> \param dimen_RI ...
69 : !> \param dimen_RI_red ...
70 : !> \param bse_lev_virt ...
71 : !> \param gd_array ...
72 : !> \param color_sub ...
73 : !> \param mp2_env ...
74 : !> \param qs_env ...
75 : !> \param mo_coeff ...
76 : !> \param unit_nr ...
77 : ! **************************************************************************************************
78 30 : SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
79 : fm_mat_Q_static_bse_gemm, &
80 : Eigenval, Eigenval_scf, &
81 30 : homo, virtual, dimen_RI, dimen_RI_red, bse_lev_virt, &
82 30 : gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
83 :
84 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, &
85 : fm_mat_S_ab_bse
86 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_Q_static_bse_gemm
87 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
88 : INTENT(IN) :: Eigenval, Eigenval_scf
89 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
90 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, bse_lev_virt
91 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
92 : INTEGER, INTENT(IN) :: color_sub
93 : TYPE(mp2_type) :: mp2_env
94 : TYPE(qs_environment_type), POINTER :: qs_env
95 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
96 : INTEGER, INTENT(IN) :: unit_nr
97 :
98 : CHARACTER(LEN=*), PARAMETER :: routineN = 'start_bse_calculation'
99 :
100 : INTEGER :: handle, homo_red, virtual_red
101 : LOGICAL :: my_do_abba, my_do_fulldiag, &
102 : my_do_iterat_diag, my_do_tda
103 : REAL(KIND=dp) :: diag_runtime_est
104 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval_reduced
105 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_abQ_bse_local, B_bar_iaQ_bse_local, &
106 30 : B_bar_ijQ_bse_local, B_iaQ_bse_local
107 : TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_C_BSE, fm_inv_sqrt_A_minus_B, fm_mat_S_ab_trunc, &
108 : fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, &
109 : fm_sqrt_A_minus_B
110 : TYPE(cp_logger_type), POINTER :: logger
111 : TYPE(mp_para_env_type), POINTER :: para_env
112 :
113 30 : CALL timeset(routineN, handle)
114 :
115 30 : para_env => fm_mat_S_ia_bse%matrix_struct%para_env
116 :
117 30 : my_do_fulldiag = .FALSE.
118 30 : my_do_iterat_diag = .FALSE.
119 30 : my_do_tda = .FALSE.
120 30 : my_do_abba = .FALSE.
121 : !Method: Iterative or full diagonalization
122 30 : SELECT CASE (mp2_env%bse%bse_diag_method)
123 : CASE (bse_iterdiag)
124 0 : my_do_iterat_diag = .TRUE.
125 : !MG: Basics of the Davidson solver are implemented, but not rigorously checked.
126 0 : CPABORT("Iterative BSE not yet implemented")
127 : CASE (bse_fulldiag)
128 30 : my_do_fulldiag = .TRUE.
129 : END SELECT
130 : !Approximation: TDA and/or full ABBA matrix
131 42 : SELECT CASE (mp2_env%bse%flag_tda)
132 : CASE (bse_tda)
133 12 : my_do_tda = .TRUE.
134 : CASE (bse_abba)
135 18 : my_do_abba = .TRUE.
136 : CASE (bse_both)
137 0 : my_do_tda = .TRUE.
138 30 : my_do_abba = .TRUE.
139 : END SELECT
140 :
141 30 : CALL print_BSE_start_flag(my_do_tda, my_do_abba, unit_nr)
142 :
143 : ! Link BSE debug flag against debug print level
144 30 : logger => cp_get_default_logger()
145 30 : IF (logger%iter_info%print_level == debug_print_level) THEN
146 0 : mp2_env%bse%bse_debug_print = .TRUE.
147 : END IF
148 :
149 30 : CALL fm_mat_S_ia_bse%matrix_struct%para_env%sync()
150 : ! We apply the BSE cutoffs using the DFT Eigenenergies
151 : ! Reduce matrices in case of energy cutoff for occupied and unoccupied in A/B-BSE-matrices
152 : CALL truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
153 : fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
154 : Eigenval_scf(:, 1, 1), Eigenval(:, 1, 1), Eigenval_reduced, &
155 : homo(1), virtual(1), dimen_RI, unit_nr, &
156 : bse_lev_virt, &
157 : homo_red, virtual_red, &
158 30 : mp2_env)
159 : ! \bar{B}^P_rs = \sum_R W_PR B^R_rs where B^R_rs = \sum_T [1/sqrt(v)]_RT (T|rs)
160 : ! r,s: MO-index, P,R,T: RI-index
161 : ! B: fm_mat_S_..., W: fm_mat_Q_...
162 : CALL mult_B_with_W(fm_mat_S_ij_trunc, fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, &
163 : fm_mat_S_bar_ij_bse, fm_mat_Q_static_bse_gemm, &
164 30 : dimen_RI_red, homo_red, virtual_red)
165 :
166 30 : IF (my_do_iterat_diag) THEN
167 : CALL fill_local_3c_arrays(fm_mat_S_ab_trunc, fm_mat_S_ia_trunc, &
168 : fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
169 : B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
170 : B_iaQ_bse_local, dimen_RI_red, homo_red, virtual_red, &
171 : gd_array, color_sub, para_env)
172 : END IF
173 :
174 30 : CALL adapt_BSE_input_params(homo_red, virtual_red, unit_nr, mp2_env, qs_env)
175 :
176 30 : IF (my_do_fulldiag) THEN
177 : ! Quick estimate of memory consumption and runtime of diagonalizations
178 : CALL estimate_BSE_resources(homo_red, virtual_red, unit_nr, my_do_abba, &
179 30 : para_env, diag_runtime_est)
180 : ! Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
181 : ! A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
182 : ! ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
183 : ! α is a spin-dependent factor
184 : ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
185 : ! W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
186 :
187 : ! For unscreened W matrix, we need fm_mat_S_ij_trunc
188 30 : IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
189 : mp2_env%bse%screening_method == bse_screening_alpha) THEN
190 : CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
191 : fm_A_BSE, Eigenval_reduced, unit_nr, &
192 : homo_red, virtual_red, dimen_RI, mp2_env, &
193 4 : para_env)
194 : ELSE
195 : CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_bar_ij_bse, fm_mat_S_ab_trunc, &
196 : fm_A_BSE, Eigenval_reduced, unit_nr, &
197 : homo_red, virtual_red, dimen_RI, mp2_env, &
198 26 : para_env)
199 : END IF
200 30 : IF (my_do_abba) THEN
201 : ! Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
202 : ! B_ia,jb = α * v_ia,jb - W_ib,aj
203 : ! α is a spin-dependent factor
204 : ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
205 : ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
206 :
207 : ! For unscreened W matrix, we need fm_mat_S_ia_trunc
208 18 : IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
209 : mp2_env%bse%screening_method == bse_screening_alpha) THEN
210 : CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_ia_trunc, fm_B_BSE, &
211 4 : homo_red, virtual_red, dimen_RI, unit_nr, mp2_env)
212 : ELSE
213 : CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, fm_B_BSE, &
214 14 : homo_red, virtual_red, dimen_RI, unit_nr, mp2_env)
215 : END IF
216 : ! Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
217 : ! (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
218 : ! We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
219 : ! of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
220 : CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, &
221 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
222 18 : homo_red, virtual_red, unit_nr, mp2_env, diag_runtime_est)
223 : END IF
224 30 : CALL cp_fm_release(fm_B_BSE)
225 30 : IF (my_do_tda) THEN
226 : ! Solving the hermitian eigenvalue equation A X^n = Ω^n X^n
227 : CALL diagonalize_A(fm_A_BSE, homo_red, virtual_red, homo(1), &
228 12 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
229 : END IF
230 : ! Release to avoid faulty use of changed A matrix
231 30 : CALL cp_fm_release(fm_A_BSE)
232 30 : IF (my_do_abba) THEN
233 : ! Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
234 : ! Here, the eigenvectors Z^n relate to X^n via
235 : ! Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
236 : CALL diagonalize_C(fm_C_BSE, homo_red, virtual_red, homo(1), &
237 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
238 18 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
239 : END IF
240 : ! Release to avoid faulty use of changed C matrix
241 30 : CALL cp_fm_release(fm_C_BSE)
242 : END IF
243 :
244 : CALL deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
245 : fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
246 30 : fm_mat_Q_static_bse_gemm, mp2_env)
247 30 : DEALLOCATE (Eigenval_reduced)
248 30 : IF (my_do_iterat_diag) THEN
249 : ! Contains untested Block-Davidson algorithm
250 : CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
251 : B_iaQ_bse_local, homo(1), virtual(1), mp2_env%bse%bse_spin_config, unit_nr, &
252 0 : Eigenval(:, 1, 1), para_env, mp2_env)
253 : ! Deallocate local 3c-B-matrices
254 0 : DEALLOCATE (B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, B_iaQ_bse_local)
255 : END IF
256 :
257 30 : IF (unit_nr > 0) THEN
258 15 : WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!'
259 : END IF
260 :
261 30 : CALL timestop(handle)
262 :
263 60 : END SUBROUTINE start_bse_calculation
264 :
265 : END MODULE bse_main
|
