Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : determine_bse_combined_window,&
27 : estimate_BSE_resources,&
28 : get_bse_spin_block_layout,&
29 : mult_B_with_W,&
30 : truncate_BSE_matrices
31 : USE cp_control_types, ONLY: dft_control_type,&
32 : tddfpt2_control_type
33 : USE cp_fm_types, ONLY: cp_fm_create,&
34 : cp_fm_release,&
35 : cp_fm_to_fm,&
36 : cp_fm_type
37 : USE cp_log_handling, ONLY: cp_get_default_logger,&
38 : cp_logger_type
39 : USE cp_output_handling, ONLY: debug_print_level
40 : USE group_dist_types, ONLY: group_dist_d1_type
41 : USE input_constants, ONLY: bse_abba,&
42 : bse_both,&
43 : bse_fulldiag,&
44 : bse_iterdiag,&
45 : bse_screening_alpha,&
46 : bse_screening_tdhf,&
47 : bse_tda
48 : USE kinds, ONLY: dp
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE mp2_types, ONLY: mp2_type
51 : USE qs_environment_types, ONLY: get_qs_env,&
52 : qs_environment_type
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main'
60 :
61 : PUBLIC :: start_bse_calculation
62 :
63 : CONTAINS
64 :
65 : ! **************************************************************************************************
66 : !> \brief Main subroutine managing BSE calculations
67 : !> \param fm_mat_S_ia_bse ...
68 : !> \param fm_mat_S_ij_bse ...
69 : !> \param fm_mat_S_ab_bse ...
70 : !> \param fm_mat_Q_static_bse_gemm ...
71 : !> \param Eigenval ...
72 : !> \param Eigenval_scf ...
73 : !> \param homo ...
74 : !> \param virtual ...
75 : !> \param dimen_RI ...
76 : !> \param dimen_RI_red ...
77 : !> \param bse_lev_virt ...
78 : !> \param gd_array ...
79 : !> \param color_sub ...
80 : !> \param mp2_env ...
81 : !> \param qs_env ...
82 : !> \param mo_coeff ...
83 : !> \param unit_nr ...
84 : ! **************************************************************************************************
85 42 : SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
86 : fm_mat_Q_static_bse_gemm, &
87 : Eigenval, Eigenval_scf, &
88 84 : homo, virtual, dimen_RI, dimen_RI_red, bse_lev_virt, &
89 42 : gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
90 :
91 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, &
92 : fm_mat_S_ab_bse
93 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_Q_static_bse_gemm
94 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
95 : INTENT(IN) :: Eigenval, Eigenval_scf
96 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
97 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
98 : INTEGER, DIMENSION(:), INTENT(IN) :: bse_lev_virt
99 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
100 : INTEGER, INTENT(IN) :: color_sub
101 : TYPE(mp2_type) :: mp2_env
102 : TYPE(qs_environment_type), POINTER :: qs_env
103 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
104 : INTEGER, INTENT(IN) :: unit_nr
105 :
106 : CHARACTER(LEN=*), PARAMETER :: routineN = 'start_bse_calculation'
107 :
108 : INTEGER :: first_active_mo, handle, ispin, &
109 : last_active_mo, n_ov_joint, nspins
110 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: homo_red_arr, n_ov_arr, offsets_arr, &
111 42 : virt_red_arr
112 : LOGICAL :: my_do_abba, my_do_fulldiag, &
113 : my_do_iterat_diag, my_do_tda
114 : REAL(KIND=dp) :: diag_runtime_est
115 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval_reduced, Eigenval_reduced_1, &
116 42 : Eigenval_reduced_2, &
117 42 : Eigenval_reduced_joint
118 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_abQ_bse_local, B_bar_iaQ_bse_local, &
119 42 : B_bar_ijQ_bse_local, B_iaQ_bse_local
120 : TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_C_BSE, &
121 : fm_inv_sqrt_A_minus_B, fm_Q_copy, &
122 : fm_sqrt_A_minus_B
123 42 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_S_ab_trunc_arr, &
124 42 : fm_mat_S_bar_ia_bse_arr, fm_mat_S_bar_ij_bse_arr, fm_mat_S_ia_trunc_arr, &
125 42 : fm_mat_S_ij_trunc_arr
126 : TYPE(cp_logger_type), POINTER :: logger
127 : TYPE(dft_control_type), POINTER :: dft_control
128 : TYPE(mp_para_env_type), POINTER :: para_env
129 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
130 :
131 42 : CALL timeset(routineN, handle)
132 :
133 42 : nspins = SIZE(homo)
134 42 : para_env => fm_mat_S_ia_bse(1)%matrix_struct%para_env
135 :
136 42 : my_do_fulldiag = .FALSE.
137 42 : my_do_iterat_diag = .FALSE.
138 42 : my_do_tda = .FALSE.
139 42 : my_do_abba = .FALSE.
140 : !Method: Iterative or full diagonalization
141 42 : SELECT CASE (mp2_env%bse%bse_diag_method)
142 : CASE (bse_iterdiag)
143 0 : my_do_iterat_diag = .TRUE.
144 : !MG: Basics of the Davidson solver are implemented, but not rigorously checked.
145 0 : CPABORT("Iterative BSE not yet implemented")
146 : CASE (bse_fulldiag)
147 42 : my_do_fulldiag = .TRUE.
148 : END SELECT
149 : !Approximation: TDA and/or full ABBA matrix
150 62 : SELECT CASE (mp2_env%bse%flag_tda)
151 : CASE (bse_tda)
152 20 : my_do_tda = .TRUE.
153 : CASE (bse_abba)
154 22 : my_do_abba = .TRUE.
155 : CASE (bse_both)
156 0 : my_do_tda = .TRUE.
157 42 : my_do_abba = .TRUE.
158 : END SELECT
159 :
160 42 : CALL print_BSE_start_flag(my_do_tda, my_do_abba, unit_nr)
161 :
162 : ! Link BSE debug flag against debug print level
163 42 : logger => cp_get_default_logger()
164 42 : IF (logger%iter_info%print_level == debug_print_level) THEN
165 0 : mp2_env%bse%bse_debug_print = .TRUE.
166 : END IF
167 :
168 42 : CALL fm_mat_S_ia_bse(1)%matrix_struct%para_env%sync()
169 :
170 168 : ALLOCATE (homo_red_arr(nspins), virt_red_arr(nspins))
171 360 : ALLOCATE (fm_mat_S_ia_trunc_arr(nspins), fm_mat_S_ij_trunc_arr(nspins), fm_mat_S_ab_trunc_arr(nspins))
172 226 : ALLOCATE (fm_mat_S_bar_ia_bse_arr(nspins), fm_mat_S_bar_ij_bse_arr(nspins))
173 :
174 42 : IF (nspins > 1) THEN
175 : CALL cp_warn(__LOCATION__, &
176 : "Open-shell (UKS/LSD) BSE is a recent addition and has not been "// &
177 : "extensively validated. Verify results carefully before using them "// &
178 8 : "for production calculations.")
179 : ! === Open-shell path: combined-window cutoff, per-spin truncation, joint A (+B for ABBA) ===
180 : ! Determine union of per-spin active MO windows
181 : CALL determine_bse_combined_window(Eigenval_scf(:, 1, :), homo, virtual, &
182 : mp2_env%bse%bse_cutoff_occ, &
183 : mp2_env%bse%bse_cutoff_empty, &
184 8 : first_active_mo, last_active_mo)
185 24 : DO ispin = 1, nspins
186 : CALL truncate_BSE_matrices(fm_mat_S_ia_bse(ispin), fm_mat_S_ij_bse(ispin), &
187 : fm_mat_S_ab_bse(ispin), &
188 : fm_mat_S_ia_trunc_arr(ispin), fm_mat_S_ij_trunc_arr(ispin), &
189 : fm_mat_S_ab_trunc_arr(ispin), &
190 : Eigenval_scf(:, 1, ispin), Eigenval(:, 1, ispin), &
191 : Eigenval_reduced, homo(ispin), virtual(ispin), dimen_RI, &
192 : unit_nr, bse_lev_virt(ispin), homo_red_arr(ispin), &
193 : virt_red_arr(ispin), mp2_env, &
194 : homo_incl_in=first_active_mo, &
195 16 : virt_incl_in=last_active_mo - homo(ispin))
196 16 : IF (ispin == 1) THEN
197 24 : ALLOCATE (Eigenval_reduced_1(SIZE(Eigenval_reduced)))
198 300 : Eigenval_reduced_1(:) = Eigenval_reduced(:)
199 : ELSE
200 24 : ALLOCATE (Eigenval_reduced_2(SIZE(Eigenval_reduced)))
201 300 : Eigenval_reduced_2(:) = Eigenval_reduced(:)
202 : END IF
203 24 : DEALLOCATE (Eigenval_reduced)
204 : END DO
205 : ! Flat eigenvalue layout: [sigma=1 levels, sigma=2 levels]
206 24 : ALLOCATE (Eigenval_reduced_joint(SIZE(Eigenval_reduced_1) + SIZE(Eigenval_reduced_2)))
207 300 : Eigenval_reduced_joint(1:SIZE(Eigenval_reduced_1)) = Eigenval_reduced_1
208 300 : Eigenval_reduced_joint(SIZE(Eigenval_reduced_1) + 1:) = Eigenval_reduced_2
209 8 : DEALLOCATE (Eigenval_reduced_1, Eigenval_reduced_2)
210 :
211 24 : ALLOCATE (n_ov_arr(nspins), offsets_arr(nspins))
212 8 : CALL get_bse_spin_block_layout(homo_red_arr, virt_red_arr, n_ov_arr, offsets_arr, n_ov_joint)
213 :
214 8 : CALL adapt_BSE_input_params(n_ov_joint, 1, unit_nr, mp2_env, qs_env)
215 :
216 : ! W: mult_B_with_W modifies Q in-place (Cholesky); copy the original Q for each spin call
217 24 : DO ispin = 1, nspins
218 16 : CALL cp_fm_create(fm_Q_copy, fm_mat_Q_static_bse_gemm%matrix_struct)
219 16 : CALL cp_fm_to_fm(fm_mat_Q_static_bse_gemm, fm_Q_copy)
220 : CALL mult_B_with_W(fm_mat_S_ij_trunc_arr(ispin), fm_mat_S_ia_trunc_arr(ispin), &
221 : fm_mat_S_bar_ia_bse_arr(ispin), fm_mat_S_bar_ij_bse_arr(ispin), &
222 16 : fm_Q_copy, dimen_RI_red, homo_red_arr(ispin), virt_red_arr(ispin))
223 40 : CALL cp_fm_release(fm_Q_copy)
224 : END DO
225 8 : CALL cp_fm_release(fm_mat_Q_static_bse_gemm)
226 :
227 8 : IF (my_do_fulldiag) THEN
228 8 : CALL estimate_BSE_resources(n_ov_joint, unit_nr, my_do_abba, para_env, diag_runtime_est)
229 8 : IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
230 : mp2_env%bse%screening_method == bse_screening_alpha) THEN
231 : CALL create_A(fm_mat_S_ia_trunc_arr, fm_mat_S_ij_trunc_arr, fm_mat_S_ab_trunc_arr, &
232 : fm_A_BSE, Eigenval_reduced_joint, unit_nr, &
233 0 : homo_red_arr, virt_red_arr, dimen_RI, mp2_env, para_env, qs_env)
234 : ELSE
235 : CALL create_A(fm_mat_S_ia_trunc_arr, fm_mat_S_bar_ij_bse_arr, fm_mat_S_ab_trunc_arr, &
236 : fm_A_BSE, Eigenval_reduced_joint, unit_nr, &
237 8 : homo_red_arr, virt_red_arr, dimen_RI, mp2_env, para_env, qs_env)
238 : END IF
239 8 : IF (my_do_abba) THEN
240 4 : IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
241 : mp2_env%bse%screening_method == bse_screening_alpha) THEN
242 : CALL create_B(fm_mat_S_ia_trunc_arr, fm_mat_S_ia_trunc_arr, fm_B_BSE, &
243 0 : homo_red_arr, virt_red_arr, dimen_RI, unit_nr, mp2_env)
244 : ELSE
245 : CALL create_B(fm_mat_S_ia_trunc_arr, fm_mat_S_bar_ia_bse_arr, fm_B_BSE, &
246 4 : homo_red_arr, virt_red_arr, dimen_RI, unit_nr, mp2_env)
247 : END IF
248 : CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, &
249 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
250 4 : unit_nr, mp2_env, diag_runtime_est)
251 4 : CALL cp_fm_release(fm_B_BSE)
252 : END IF
253 8 : NULLIFY (dft_control, tddfpt_control)
254 8 : CALL get_qs_env(qs_env, dft_control=dft_control)
255 8 : tddfpt_control => dft_control%tddfpt2_control
256 : ! 4th arg (homo_irred) = full per-spin occupied counts: per-spin absolute-MO labels for
257 : ! the amplitude table, and N_e = SUM(homo_irred) for the TRK print.
258 8 : IF (my_do_tda .AND. (.NOT. tddfpt_control%do_bse)) THEN
259 : CALL diagonalize_A(fm_A_BSE, homo_red_arr, virt_red_arr, homo, &
260 4 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
261 : END IF
262 8 : CALL cp_fm_release(fm_A_BSE)
263 8 : IF (my_do_abba) THEN
264 : CALL diagonalize_C(fm_C_BSE, homo_red_arr, virt_red_arr, homo, &
265 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
266 4 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
267 4 : CALL cp_fm_release(fm_C_BSE)
268 : END IF
269 : END IF
270 :
271 24 : DO ispin = 1, nspins
272 16 : CALL cp_fm_release(fm_mat_S_bar_ia_bse_arr(ispin))
273 16 : CALL cp_fm_release(fm_mat_S_bar_ij_bse_arr(ispin))
274 16 : CALL cp_fm_release(fm_mat_S_ia_trunc_arr(ispin))
275 16 : CALL cp_fm_release(fm_mat_S_ij_trunc_arr(ispin))
276 24 : CALL cp_fm_release(fm_mat_S_ab_trunc_arr(ispin))
277 : END DO
278 8 : IF (mp2_env%bse%do_nto_analysis) DEALLOCATE (mp2_env%bse%bse_nto_state_list_final)
279 16 : DEALLOCATE (Eigenval_reduced_joint, n_ov_arr, offsets_arr)
280 :
281 : ELSE
282 : ! === Closed-shell n_spin=1 path (bit-identical) ===
283 : CALL truncate_BSE_matrices(fm_mat_S_ia_bse(1), fm_mat_S_ij_bse(1), fm_mat_S_ab_bse(1), &
284 : fm_mat_S_ia_trunc_arr(1), fm_mat_S_ij_trunc_arr(1), &
285 : fm_mat_S_ab_trunc_arr(1), &
286 : Eigenval_scf(:, 1, 1), Eigenval(:, 1, 1), Eigenval_reduced, &
287 : homo(1), virtual(1), dimen_RI, unit_nr, &
288 34 : bse_lev_virt(1), homo_red_arr(1), virt_red_arr(1), mp2_env)
289 : CALL mult_B_with_W(fm_mat_S_ij_trunc_arr(1), fm_mat_S_ia_trunc_arr(1), &
290 : fm_mat_S_bar_ia_bse_arr(1), fm_mat_S_bar_ij_bse_arr(1), &
291 34 : fm_mat_Q_static_bse_gemm, dimen_RI_red, homo_red_arr(1), virt_red_arr(1))
292 :
293 34 : IF (my_do_iterat_diag) THEN
294 : CALL fill_local_3c_arrays(fm_mat_S_ab_trunc_arr(1), fm_mat_S_ia_trunc_arr(1), &
295 : fm_mat_S_bar_ia_bse_arr(1), fm_mat_S_bar_ij_bse_arr(1), &
296 : B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
297 : B_iaQ_bse_local, dimen_RI_red, homo_red_arr(1), &
298 0 : virt_red_arr(1), gd_array, color_sub, para_env)
299 : END IF
300 :
301 34 : CALL adapt_BSE_input_params(homo_red_arr(1), virt_red_arr(1), unit_nr, mp2_env, qs_env)
302 :
303 34 : IF (my_do_fulldiag) THEN
304 34 : n_ov_joint = homo_red_arr(1)*virt_red_arr(1)
305 34 : CALL estimate_BSE_resources(n_ov_joint, unit_nr, my_do_abba, para_env, diag_runtime_est)
306 34 : IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
307 : mp2_env%bse%screening_method == bse_screening_alpha) THEN
308 : CALL create_A(fm_mat_S_ia_trunc_arr, fm_mat_S_ij_trunc_arr, fm_mat_S_ab_trunc_arr, &
309 : fm_A_BSE, Eigenval_reduced, unit_nr, &
310 4 : homo_red_arr, virt_red_arr, dimen_RI, mp2_env, para_env, qs_env)
311 : ELSE
312 : CALL create_A(fm_mat_S_ia_trunc_arr, fm_mat_S_bar_ij_bse_arr, fm_mat_S_ab_trunc_arr, &
313 : fm_A_BSE, Eigenval_reduced, unit_nr, &
314 30 : homo_red_arr, virt_red_arr, dimen_RI, mp2_env, para_env, qs_env)
315 : END IF
316 34 : IF (my_do_abba) THEN
317 18 : IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
318 : mp2_env%bse%screening_method == bse_screening_alpha) THEN
319 : CALL create_B(fm_mat_S_ia_trunc_arr, fm_mat_S_ia_trunc_arr, fm_B_BSE, &
320 4 : homo_red_arr, virt_red_arr, dimen_RI, unit_nr, mp2_env)
321 : ELSE
322 : CALL create_B(fm_mat_S_ia_trunc_arr, fm_mat_S_bar_ia_bse_arr, fm_B_BSE, &
323 14 : homo_red_arr, virt_red_arr, dimen_RI, unit_nr, mp2_env)
324 : END IF
325 : CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, &
326 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
327 18 : unit_nr, mp2_env, diag_runtime_est)
328 : END IF
329 34 : CALL cp_fm_release(fm_B_BSE)
330 :
331 34 : NULLIFY (dft_control, tddfpt_control)
332 34 : CALL get_qs_env(qs_env, dft_control=dft_control)
333 34 : tddfpt_control => dft_control%tddfpt2_control
334 34 : IF (my_do_tda .AND. (.NOT. tddfpt_control%do_bse)) THEN
335 : CALL diagonalize_A(fm_A_BSE, homo_red_arr, virt_red_arr, homo, &
336 14 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
337 : END IF
338 34 : CALL cp_fm_release(fm_A_BSE)
339 34 : IF (my_do_abba) THEN
340 : CALL diagonalize_C(fm_C_BSE, homo_red_arr, virt_red_arr, homo, &
341 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
342 18 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
343 : END IF
344 34 : CALL cp_fm_release(fm_C_BSE)
345 : END IF
346 :
347 : CALL deallocate_matrices_bse(fm_mat_S_bar_ia_bse_arr(1), fm_mat_S_bar_ij_bse_arr(1), &
348 : fm_mat_S_ia_trunc_arr(1), fm_mat_S_ij_trunc_arr(1), &
349 34 : fm_mat_S_ab_trunc_arr(1), fm_mat_Q_static_bse_gemm, mp2_env)
350 34 : DEALLOCATE (Eigenval_reduced)
351 34 : IF (my_do_iterat_diag) THEN
352 : CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
353 : B_iaQ_bse_local, homo(1), virtual(1), &
354 : mp2_env%bse%bse_spin_config, unit_nr, &
355 0 : Eigenval(:, 1, 1), para_env, mp2_env)
356 0 : DEALLOCATE (B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, B_iaQ_bse_local)
357 : END IF
358 :
359 : END IF
360 :
361 42 : DEALLOCATE (homo_red_arr, virt_red_arr)
362 42 : DEALLOCATE (fm_mat_S_ia_trunc_arr, fm_mat_S_ij_trunc_arr, fm_mat_S_ab_trunc_arr)
363 42 : DEALLOCATE (fm_mat_S_bar_ia_bse_arr, fm_mat_S_bar_ij_bse_arr)
364 :
365 42 : IF (unit_nr > 0) THEN
366 21 : WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!'
367 : END IF
368 :
369 42 : CALL timestop(handle)
370 :
371 84 : END SUBROUTINE start_bse_calculation
372 :
373 : END MODULE bse_main
|
