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: a GW path
10 : !> prepares the BSE environment, bse_solve runs the solver on it
11 : !> \par History
12 : !> 04.2024 created [Maximilian Graml]
13 : !> 09.2026 split into prepare_bse_env and bse_solve [Maximilian Graml]
14 : ! **************************************************************************************************
15 :
16 : MODULE bse_main
17 :
18 : USE bse_full_diag, ONLY: create_A_and_B,&
19 : create_hermitian_form_of_ABBA,&
20 : diagonalize_A,&
21 : diagonalize_C
22 : USE bse_input, ONLY: read_bse_section
23 : USE bse_iterative, ONLY: solve_bse_iteratively
24 : USE bse_matvec, ONLY: bse_memory_left_GB,&
25 : bse_resolve_memory_per_proc
26 : USE bse_print, ONLY: print_BSE_start_flag
27 : USE bse_types, ONLY: bse_env_type
28 : USE bse_util, ONLY: &
29 : adapt_BSE_input_params, adapt_BSE_screening, bse_window_from_dft, bse_window_from_mos, &
30 : estimate_BSE_resources, fulldiag_memory_estimate_GB, get_bse_spin_block_layout, &
31 : screen_slabs, truncate_BSE_matrices
32 : USE cp_control_types, ONLY: tddfpt2_control_type
33 : USE cp_fm_basic_linalg, ONLY: cp_fm_add_on_diag
34 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
35 : cp_fm_cholesky_invert
36 : USE cp_fm_types, ONLY: cp_fm_create,&
37 : cp_fm_release,&
38 : cp_fm_to_fm,&
39 : cp_fm_type
40 : USE cp_log_handling, ONLY: cp_get_default_logger,&
41 : cp_logger_type
42 : USE cp_output_handling, ONLY: debug_print_level
43 : USE input_constants, ONLY: bse_abba,&
44 : bse_both,&
45 : bse_fulldiag,&
46 : bse_iterdiag,&
47 : bse_screening_w0,&
48 : bse_tda
49 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
50 : section_vals_type
51 : USE kinds, ONLY: dp
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE mo_window, ONLY: mo_window_type
54 : USE pw_poisson_read_input, ONLY: pw_poisson_read_parameters
55 : USE pw_poisson_types, ONLY: pw_poisson_parameter_type
56 : USE qs_environment_types, ONLY: get_qs_env,&
57 : qs_environment_type
58 : USE qs_mo_types, ONLY: get_mo_set,&
59 : mo_set_type
60 : USE qs_moments, ONLY: get_multipoles_ao
61 : #include "./base/base_uses.f90"
62 :
63 : IMPLICIT NONE
64 :
65 : PRIVATE
66 :
67 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main'
68 :
69 : PUBLIC :: prepare_bse_env_from_ri_rpa_gw, prepare_bse_env_from_bandstructure_gw, bse_solve, &
70 : bse_fulldiag_memory_check, bse_fulldiag_memory_check_from_mos
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Entry for the BSE of RI_RPA%GW, which hands over the slabs of all its GW levels and the
76 : !> static polarizability Q(0): re-reads the settings, cuts the window on the DFT axis out of
77 : !> the slabs and both energy sets, and inverts 1+Q(0) for W0 screening. Starts no solver.
78 : !> \param bse_env settings re-read from its section; state filled here
79 : !> \param fm_S_ia_full B^P_ia of the GW levels per spin, the RPA copy without the energy factor
80 : !> \param fm_S_ij_full B^P_ij per spin
81 : !> \param fm_S_ab_full B^P_ab per spin
82 : !> \param fm_Q the static polarizability Q(0), untouched; the caller releases it afterwards
83 : !> \param Eigenval quasiparticle energies (level, 1, spin)
84 : !> \param Eigenval_scf DFT energies (level, 1, spin)
85 : !> \param homo occupied MOs per spin
86 : !> \param virtual virtual MOs per spin
87 : !> \param dimen_RI RI dimension of the slabs
88 : !> \param dimen_RI_red reduced RI dimension of Q(0)
89 : !> \param gw_corr_lev_occ GW-corrected occupied levels; the window may not exceed them
90 : !> \param bse_lev_virt GW-corrected virtual levels per spin
91 : !> \param qs_env source of the pointers the solver reads and of the AO multipoles
92 : !> \param unit_nr output unit, positive on the writing rank only
93 : ! **************************************************************************************************
94 112 : SUBROUTINE prepare_bse_env_from_ri_rpa_gw(bse_env, fm_S_ia_full, fm_S_ij_full, fm_S_ab_full, fm_Q, &
95 168 : Eigenval, Eigenval_scf, homo, virtual, dimen_RI, dimen_RI_red, &
96 56 : gw_corr_lev_occ, bse_lev_virt, qs_env, unit_nr)
97 :
98 : TYPE(bse_env_type), INTENT(INOUT) :: bse_env
99 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_S_ia_full, fm_S_ij_full, fm_S_ab_full
100 : TYPE(cp_fm_type), INTENT(IN) :: fm_Q
101 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval, Eigenval_scf
102 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
103 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, gw_corr_lev_occ
104 : INTEGER, DIMENSION(:), INTENT(IN) :: bse_lev_virt
105 : TYPE(qs_environment_type), POINTER :: qs_env
106 : INTEGER, INTENT(IN) :: unit_nr
107 :
108 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_bse_env_from_ri_rpa_gw'
109 :
110 : INTEGER :: handle, ispin, n_window, nspins
111 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_mo
112 56 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_dft_spin, eigenval_qp_spin
113 : TYPE(mo_window_type) :: window
114 : TYPE(section_vals_type), POINTER :: bse_section
115 :
116 56 : CALL timeset(routineN, handle)
117 :
118 : ! Re-read the settings afresh from the input, since adapt_BSE_screening here and
119 : ! adapt_BSE_input_params in bse_solve normalise them in place, and a later energy evaluation
120 : ! on the same environment runs the BSE again; the local pointer keeps the section apart from
121 : ! bse_env, whose section pointer the reader sets
122 :
123 56 : bse_section => bse_env%bse_section
124 56 : CALL read_bse_section(bse_section, bse_env)
125 56 : CALL adapt_BSE_screening(bse_env, unit_nr)
126 56 : CALL bse_resolve_memory_per_proc(bse_env%memory_per_proc_gb, fm_S_ia_full(1)%matrix_struct%para_env, unit_nr)
127 :
128 56 : nspins = SIZE(homo)
129 : ! the GW-corrected level counts are the ceiling of the window
130 56 : CALL start_bse_env(bse_env, nspins, dimen_RI, dimen_RI_red, gw_corr_lev_occ, bse_lev_virt, unit_nr)
131 :
132 : ! one absolute window on the DFT axis, then the window-sized slabs and both energy sets per spin
133 168 : ALLOCATE (n_mo(nspins))
134 124 : n_mo(:) = homo(:) + virtual(:)
135 : CALL bse_window_from_dft(Eigenval_scf(:, 1, :), n_mo, homo, bse_env%bse_cutoff_occ, &
136 56 : bse_env%bse_cutoff_empty, window, warn=.TRUE.)
137 56 : DEALLOCATE (n_mo)
138 :
139 224 : ALLOCATE (bse_env%homo_red(nspins), bse_env%virt_red(nspins), bse_env%homo_full(nspins))
140 484 : ALLOCATE (bse_env%fm_S_ia(nspins), bse_env%fm_S_ij(nspins), bse_env%fm_S_ab(nspins))
141 124 : bse_env%homo_full(:) = homo(:)
142 124 : DO ispin = 1, nspins
143 : ! the cutoff summary is printed for a closed-shell reference only
144 : CALL truncate_BSE_matrices(fm_S_ia_full(ispin), fm_S_ij_full(ispin), fm_S_ab_full(ispin), &
145 : Eigenval_scf(:, 1, ispin), Eigenval(:, 1, ispin), &
146 : eigenval_dft_spin, eigenval_qp_spin, homo(ispin), virtual(ispin), &
147 68 : dimen_RI, unit_nr, ispin, bse_env, window, nspins == 1)
148 : ! the window is one absolute range, so every spin has the same number of levels
149 68 : IF (ispin == 1) THEN
150 56 : n_window = SIZE(eigenval_qp_spin)
151 336 : ALLOCATE (bse_env%eigenval_dft(n_window, nspins), bse_env%eigenval_qp(n_window, nspins))
152 : END IF
153 1870 : bse_env%eigenval_dft(:, ispin) = eigenval_dft_spin(:)
154 1870 : bse_env%eigenval_qp(:, ispin) = eigenval_qp_spin(:)
155 124 : DEALLOCATE (eigenval_dft_spin, eigenval_qp_spin)
156 : END DO
157 :
158 56 : IF (bse_env%screening_method == bse_screening_w0) THEN
159 50 : CALL invert_static_screening(fm_Q, dimen_RI_red, bse_env%fm_eps_inv)
160 : END IF
161 :
162 56 : CALL finish_bse_env(bse_env, qs_env, fm_S_ia_full(1)%matrix_struct%para_env)
163 :
164 56 : CALL timestop(handle)
165 :
166 112 : END SUBROUTINE prepare_bse_env_from_ri_rpa_gw
167 :
168 : ! **************************************************************************************************
169 : !> \brief Entry for the BSE of PROPERTIES%BANDSTRUCTURE%GW (PROPERTIES%BSE), which hands over
170 : !> window-sized slabs and a ready inverse dielectric matrix: nothing is cut or inverted, the
171 : !> slabs are taken over as they are. Shares start_bse_env and finish_bse_env with
172 : !> prepare_bse_env_from_ri_rpa_gw. Starts no solver.
173 : !> \param bse_env settings read and screening normalised by the caller; state filled here
174 : !> \param fm_S_ia B^P_ia of the window per spin; moved into bse_env, deallocated on return
175 : !> \param fm_S_ij B^P_ij of the window per spin; moved into bse_env, deallocated on return
176 : !> \param fm_S_ab B^P_ab of the window per spin; moved into bse_env, deallocated on return
177 : !> \param fm_eps_inv the inverse dielectric matrix ε^-1_PQ the core screens the slabs with,
178 : !> \bar{B}^P_ij = sum_Q ε^-1_PQ B^Q_ij; RI x RI, copied for W0 only; stays the caller's
179 : !> May be uncreated for other screening methods.
180 : !> \param eigenval_qp quasiparticle energies of the window (level, spin)
181 : !> \param eigenval_dft DFT energies of the window (level, spin)
182 : !> \param homo_full occupied MOs per spin of the full reference
183 : !> \param homo_red occupied levels in the window per spin
184 : !> \param virt_red virtual levels in the window per spin
185 : !> \param dimen_RI RI dimension of the slabs and of fm_eps_inv
186 : !> \param qs_env source of the pointers the solver reads and of the AO multipoles
187 : !> \param unit_nr output unit, positive on the writing rank only
188 : ! **************************************************************************************************
189 6 : SUBROUTINE prepare_bse_env_from_bandstructure_gw(bse_env, fm_S_ia, fm_S_ij, fm_S_ab, fm_eps_inv, &
190 6 : eigenval_qp, eigenval_dft, homo_full, homo_red, virt_red, &
191 : dimen_RI, qs_env, unit_nr)
192 :
193 : TYPE(bse_env_type), INTENT(INOUT) :: bse_env
194 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
195 : INTENT(INOUT) :: fm_S_ia, fm_S_ij, fm_S_ab
196 : TYPE(cp_fm_type), INTENT(IN) :: fm_eps_inv
197 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eigenval_qp, eigenval_dft
198 : INTEGER, DIMENSION(:), INTENT(IN) :: homo_full, homo_red, virt_red
199 : INTEGER, INTENT(IN) :: dimen_RI
200 : TYPE(qs_environment_type), POINTER :: qs_env
201 : INTEGER, INTENT(IN) :: unit_nr
202 :
203 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_bse_env_from_bandstructure_gw'
204 :
205 : INTEGER :: handle, ispin, n_window, nspins
206 :
207 6 : CALL timeset(routineN, handle)
208 :
209 : ! the caller (bse_from_bs_env) read the settings afresh from the input; the entry of RI_RPA%GW reads them itself
210 6 : IF (.NOT. ASSOCIATED(bse_env%bse_section)) CPABORT("BSE settings not read: read_bse_section comes first")
211 :
212 6 : nspins = SIZE(homo_full)
213 6 : IF (SIZE(fm_S_ia) /= nspins .OR. SIZE(eigenval_qp, 2) /= nspins) THEN
214 0 : CPABORT("BSE needs one set of RI integrals and one set of energies per spin")
215 : END IF
216 : ! every window level carries a GW energy here, and the inverse dielectric matrix has the full RI dimension
217 14 : CALL start_bse_env(bse_env, nspins, dimen_RI, dimen_RI, MAXVAL(homo_red), virt_red, unit_nr)
218 :
219 30 : ALLOCATE (bse_env%homo_red(nspins), bse_env%virt_red(nspins), bse_env%homo_full(nspins))
220 14 : bse_env%homo_full(:) = homo_full(:)
221 14 : bse_env%homo_red(:) = homo_red(:)
222 14 : bse_env%virt_red(:) = virt_red(:)
223 6 : n_window = SIZE(eigenval_qp, 1)
224 14 : DO ispin = 1, nspins
225 : ! one absolute window for every spin: its occupied and virtual levels fill it
226 14 : IF (homo_red(ispin) + virt_red(ispin) /= n_window) THEN
227 0 : CPABORT("Occupied and virtual levels do not fill the BSE window")
228 : END IF
229 : END DO
230 36 : ALLOCATE (bse_env%eigenval_dft(n_window, nspins), bse_env%eigenval_qp(n_window, nspins))
231 142 : bse_env%eigenval_dft(:, :) = eigenval_dft(:, :)
232 142 : bse_env%eigenval_qp(:, :) = eigenval_qp(:, :)
233 :
234 : ! the slabs are window-sized already: take them over
235 6 : CALL MOVE_ALLOC(fm_S_ia, bse_env%fm_S_ia)
236 6 : CALL MOVE_ALLOC(fm_S_ij, bse_env%fm_S_ij)
237 6 : CALL MOVE_ALLOC(fm_S_ab, bse_env%fm_S_ab)
238 6 : IF (bse_env%screening_method == bse_screening_w0) THEN
239 6 : CALL cp_fm_create(bse_env%fm_eps_inv, fm_eps_inv%matrix_struct)
240 6 : CALL cp_fm_to_fm(fm_eps_inv, bse_env%fm_eps_inv)
241 : END IF
242 :
243 6 : CALL finish_bse_env(bse_env, qs_env, bse_env%fm_S_ia(1)%matrix_struct%para_env)
244 :
245 6 : CALL timestop(handle)
246 :
247 6 : END SUBROUTINE prepare_bse_env_from_bandstructure_gw
248 :
249 : ! **************************************************************************************************
250 : !> \brief The part both entries share before their window: the sizes, the run flags from the
251 : !> settings, the start banner and the debug flag
252 : !> \param bse_env settings read from its section; sizes and flags set here
253 : !> \param nspins spin channels
254 : !> \param dimen_RI RI dimension of the slabs
255 : !> \param dimen_RI_red RI dimension of the inverse dielectric matrix
256 : !> \param gw_corr_lev_occ GW-corrected occupied levels
257 : !> \param gw_corr_lev_virt GW-corrected virtual levels per spin
258 : !> \param unit_nr output unit, positive on the writing rank only
259 : ! **************************************************************************************************
260 62 : SUBROUTINE start_bse_env(bse_env, nspins, dimen_RI, dimen_RI_red, gw_corr_lev_occ, gw_corr_lev_virt, unit_nr)
261 :
262 : TYPE(bse_env_type), INTENT(INOUT) :: bse_env
263 : INTEGER, INTENT(IN) :: nspins, dimen_RI, dimen_RI_red, &
264 : gw_corr_lev_occ
265 : INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_virt
266 : INTEGER, INTENT(IN) :: unit_nr
267 :
268 : TYPE(cp_logger_type), POINTER :: logger
269 :
270 62 : bse_env%nspins = nspins
271 62 : bse_env%dimen_RI = dimen_RI
272 62 : bse_env%dimen_RI_red = dimen_RI_red
273 62 : bse_env%gw_corr_lev_occ = gw_corr_lev_occ
274 186 : ALLOCATE (bse_env%gw_corr_lev_virt(nspins))
275 138 : bse_env%gw_corr_lev_virt(:) = gw_corr_lev_virt(:)
276 :
277 : ! run flags from the settings
278 62 : bse_env%do_iterdiag = bse_env%bse_diag_method == bse_iterdiag
279 62 : bse_env%do_fulldiag = bse_env%bse_diag_method == bse_fulldiag
280 62 : bse_env%do_tda = bse_env%flag_tda == bse_tda .OR. bse_env%flag_tda == bse_both
281 62 : bse_env%do_abba = bse_env%flag_tda == bse_abba .OR. bse_env%flag_tda == bse_both
282 :
283 62 : CALL print_BSE_start_flag(bse_env%do_tda, bse_env%do_abba, unit_nr)
284 :
285 : ! the debug flag follows the debug print level
286 62 : logger => cp_get_default_logger()
287 62 : IF (logger%iter_info%print_level == debug_print_level) THEN
288 0 : bse_env%bse_debug_print = .TRUE.
289 : END IF
290 :
291 62 : END SUBROUTINE start_bse_env
292 :
293 : ! **************************************************************************************************
294 : !> \brief The part both entries share once the window and the slabs are in place: the
295 : !> QS pointers, the open-shell checks, the spin-block layout and the AO multipoles
296 : !> \param bse_env state filled here
297 : !> \param qs_env source of the pointers the solver reads and of the AO multipoles
298 : !> \param para_env_slabs the communicator the slabs live on
299 : ! **************************************************************************************************
300 62 : SUBROUTINE finish_bse_env(bse_env, qs_env, para_env_slabs)
301 :
302 : TYPE(bse_env_type), INTENT(INOUT) :: bse_env
303 : TYPE(qs_environment_type), POINTER :: qs_env
304 : TYPE(mp_para_env_type), POINTER :: para_env_slabs
305 :
306 : CHARACTER(LEN=*), PARAMETER :: routineN = 'finish_bse_env'
307 :
308 : INTEGER :: handle, nspins
309 62 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_ov, offsets
310 : TYPE(pw_poisson_parameter_type) :: poisson_params
311 : TYPE(section_vals_type), POINTER :: input, poisson_section
312 :
313 62 : CALL timeset(routineN, handle)
314 :
315 62 : nspins = bse_env%nspins
316 :
317 : ! what the solver reads from QS and the GW (pointers, none owned)
318 62 : bse_env%para_env => para_env_slabs
319 : CALL get_qs_env(qs_env, subsys=bse_env%subsys, pw_env=bse_env%pw_env, &
320 : dft_control=bse_env%dft_control, mos=bse_env%mos, &
321 62 : exstate_env=bse_env%exstate_env)
322 : ! the Poisson periodicity from the input: PROPERTIES%BSE runs after the GW released the grids
323 62 : NULLIFY (input, poisson_section)
324 62 : CALL get_qs_env(qs_env, input=input)
325 62 : poisson_section => section_vals_get_subs_vals(input, "DFT%POISSON")
326 62 : CALL pw_poisson_read_parameters(poisson_section, poisson_params)
327 248 : bse_env%n_periodic_poisson = COUNT(poisson_params%periodic == 1)
328 62 : CALL bse_env%para_env%sync()
329 :
330 62 : IF (nspins > 1) THEN
331 14 : IF (bse_env%do_iterdiag) THEN
332 0 : CPABORT("Iterative BSE is not yet available for open-shell references")
333 : END IF
334 : CALL cp_warn(__LOCATION__, &
335 : "Open-shell (UKS/LSD) BSE is a recent addition and has not been "// &
336 : "extensively validated. Verify results carefully before using them "// &
337 14 : "for production calculations.")
338 : END IF
339 :
340 248 : ALLOCATE (n_ov(nspins), offsets(nspins))
341 62 : CALL get_bse_spin_block_layout(bse_env%homo_red, bse_env%virt_red, n_ov, offsets, bse_env%n_ov_joint)
342 62 : DEALLOCATE (n_ov, offsets)
343 :
344 : ! the AO multipoles for the transition dipoles and, when asked for, the exciton descriptors;
345 : ! the open-shell path computes no descriptors
346 62 : CALL get_multipoles_ao(qs_env, 1, bse_env%matrix_dipole_ao, bse_env%multipole_ref_point)
347 62 : IF (bse_env%num_print_exc_descr /= 0 .AND. nspins == 1) THEN
348 4 : CALL get_multipoles_ao(qs_env, 2, bse_env%matrix_quadpole_ao, bse_env%multipole_ref_point)
349 : END IF
350 :
351 62 : CALL timestop(handle)
352 :
353 248 : END SUBROUTINE finish_bse_env
354 :
355 : ! **************************************************************************************************
356 : !> \brief The inverse dielectric matrix in the RI basis, [1+Q(0)]^-1_PQ on the leading dimen_RI_red
357 : !> block, P, Q ≤ dimen_RI_red, by Cholesky decomposition and inversion of a copy of Q(0)
358 : !> \param fm_Q the static polarizability Q_PQ(0), untouched
359 : !> \param dimen_RI_red reduced RI dimension of Q(0)
360 : !> \param fm_eps_inv [1+Q(0)]^-1_PQ, created here on fm_Q's structure
361 : ! **************************************************************************************************
362 200 : SUBROUTINE invert_static_screening(fm_Q, dimen_RI_red, fm_eps_inv)
363 :
364 : TYPE(cp_fm_type), INTENT(IN) :: fm_Q
365 : INTEGER, INTENT(IN) :: dimen_RI_red
366 : TYPE(cp_fm_type), INTENT(OUT) :: fm_eps_inv
367 :
368 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invert_static_screening'
369 :
370 : INTEGER :: handle, info_chol
371 : TYPE(cp_fm_type) :: fm_work
372 :
373 50 : CALL timeset(routineN, handle)
374 :
375 : ! a copy: the caller keeps Q(0), and the open-shell path inverts once for both spins
376 50 : CALL cp_fm_create(fm_eps_inv, fm_Q%matrix_struct)
377 50 : CALL cp_fm_to_fm(fm_Q, fm_eps_inv)
378 50 : CALL cp_fm_create(fm_work, fm_Q%matrix_struct)
379 :
380 : ! δ_PQ + Q_PQ(0) on the leading dimen_RI_red block
381 50 : CALL cp_fm_add_on_diag(fm_eps_inv, 1.0_dp, n_active=dimen_RI_red)
382 :
383 : ! [1+Q(0)]^-1_PQ by Cholesky, symmetrised from the stored triangle
384 50 : CALL cp_fm_cholesky_decompose(matrix=fm_eps_inv, n=dimen_RI_red, info_out=info_chol)
385 50 : IF (info_chol /= 0) THEN
386 0 : CALL cp_abort(__LOCATION__, 'Cholesky decomposition failed for static polarization in BSE')
387 : END IF
388 50 : CALL cp_fm_cholesky_invert(fm_eps_inv, n=dimen_RI_red, work=fm_work)
389 :
390 50 : CALL cp_fm_release(fm_work)
391 :
392 50 : CALL timestop(handle)
393 :
394 50 : END SUBROUTINE invert_static_screening
395 :
396 : ! **************************************************************************************************
397 : !> \brief Solves the BSE on a prepared environment: input normalisation, the screened slabs from
398 : !> [1+Q(0)]^-1, then the FULLDIAG or ITERDIAG solver on the energies USE_KS_ENERGIES selects
399 : !> \param bse_env prepared by a GW path; the screened slabs are formed here
400 : !> \param mo_coeff MO coefficients of the reference, per spin, for the properties
401 : !> \param unit_nr output unit, positive on the writing rank only
402 : ! **************************************************************************************************
403 62 : SUBROUTINE bse_solve(bse_env, mo_coeff, unit_nr)
404 :
405 : TYPE(bse_env_type), INTENT(INOUT) :: bse_env
406 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
407 : INTEGER, INTENT(IN) :: unit_nr
408 :
409 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_solve'
410 :
411 : INTEGER :: handle, ispin, n_window, nspins
412 : REAL(KIND=dp) :: diag_runtime_est
413 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_joint
414 : TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_C_BSE, &
415 : fm_inv_sqrt_A_minus_B, &
416 : fm_sqrt_A_minus_B
417 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
418 :
419 62 : CALL timeset(routineN, handle)
420 :
421 62 : nspins = bse_env%nspins
422 :
423 : ! normalisation touches no slab; n_ov_joint x 1 is the transition count on either spin path
424 62 : CALL adapt_BSE_input_params(bse_env%n_ov_joint, 1, unit_nr, bse_env)
425 :
426 : ! the energies the solver sees, flat joint layout [spin 1 window | spin 2 window]
427 62 : n_window = SIZE(bse_env%eigenval_qp, 1)
428 186 : ALLOCATE (eigenval_joint(n_window*nspins))
429 138 : DO ispin = 1, nspins
430 138 : IF (bse_env%use_ks_energies) THEN
431 116 : eigenval_joint((ispin - 1)*n_window + 1:ispin*n_window) = bse_env%eigenval_dft(:, ispin)
432 : ELSE
433 1890 : eigenval_joint((ispin - 1)*n_window + 1:ispin*n_window) = bse_env%eigenval_qp(:, ispin)
434 : END IF
435 : END DO
436 :
437 : ! Keep default-initialised containers for the common solver interfaces. Only W0 needs
438 : ! screened matrix storage; TDA needs barB_ij, and ABBA additionally needs barB_ia.
439 400 : ALLOCATE (bse_env%fm_S_bar_ia(nspins), bse_env%fm_S_bar_ij(nspins))
440 62 : IF (bse_env%screening_method == bse_screening_w0) THEN
441 126 : DO ispin = 1, nspins
442 : CALL screen_slabs(bse_env%fm_eps_inv, bse_env%fm_S_ij(ispin), bse_env%fm_S_ia(ispin), &
443 : bse_env%dimen_RI_red, bse_env%homo_red(ispin), bse_env%virt_red(ispin), &
444 126 : bse_env%fm_S_bar_ij(ispin), bse_env%fm_S_bar_ia(ispin), bse_env%do_abba)
445 : END DO
446 : END IF
447 62 : CALL cp_fm_release(bse_env%fm_eps_inv)
448 :
449 62 : IF (bse_env%do_iterdiag) THEN
450 : CALL solve_bse_iteratively(eigenval_joint, &
451 : bse_env%homo_red, bse_env%virt_red, bse_env%homo_full, &
452 : bse_env%dimen_RI, bse_env%do_tda, bse_env%do_abba, &
453 8 : bse_env, mo_coeff, unit_nr)
454 : END IF
455 :
456 62 : IF (bse_env%do_fulldiag) THEN
457 : CALL estimate_BSE_resources(bse_env%n_ov_joint, unit_nr, bse_env%do_abba, &
458 54 : bse_env%para_env, diag_runtime_est)
459 : CALL create_A_and_B(bse_env%fm_S_ia, bse_env%fm_S_ij, bse_env%fm_S_ab, &
460 : bse_env%fm_S_bar_ia, bse_env%fm_S_bar_ij, fm_A_BSE, fm_B_BSE, &
461 : bse_env%do_abba, eigenval_joint, unit_nr, bse_env%homo_red, &
462 54 : bse_env%virt_red, bse_env%dimen_RI, bse_env)
463 54 : IF (bse_env%do_abba) THEN
464 : CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, &
465 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
466 26 : unit_nr, bse_env, diag_runtime_est)
467 : END IF
468 54 : CALL cp_fm_release(fm_B_BSE)
469 54 : tddfpt_control => bse_env%dft_control%tddfpt2_control
470 : ! homo_full: per-spin absolute MO labels for the amplitude table, N_e = SUM(homo_full) for TRK
471 54 : IF (bse_env%do_tda .AND. (.NOT. tddfpt_control%do_bse)) THEN
472 : CALL diagonalize_A(fm_A_BSE, bse_env%homo_red, bse_env%virt_red, bse_env%homo_full, &
473 26 : unit_nr, diag_runtime_est, bse_env, mo_coeff)
474 : END IF
475 54 : CALL cp_fm_release(fm_A_BSE)
476 54 : IF (bse_env%do_abba) THEN
477 : CALL diagonalize_C(fm_C_BSE, bse_env%homo_red, bse_env%virt_red, bse_env%homo_full, &
478 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
479 26 : unit_nr, diag_runtime_est, bse_env, mo_coeff)
480 : END IF
481 54 : CALL cp_fm_release(fm_C_BSE)
482 : END IF
483 :
484 62 : DEALLOCATE (eigenval_joint)
485 :
486 62 : IF (unit_nr > 0) THEN
487 31 : WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!'
488 : END IF
489 :
490 62 : CALL timestop(handle)
491 :
492 62 : END SUBROUTINE bse_solve
493 :
494 : ! **************************************************************************************************
495 : !> \brief Early memory check for the full diagonalization, before any GW work: the window on the
496 : !> DFT axis gives n_ov = sum_spin n_occ n_virt, the estimate is estimate_BSE_resources', and
497 : !> the run stops when it exceeds what MEMORY_PER_PROC leaves after the memory the rank holds,
498 : !> since a full diagonalization has nothing to shrink
499 : !> \param bse_env settings; the TDA flag and MEMORY_PER_PROC, resolved here when still -1
500 : !> \param window the active MO window on the DFT axis, one absolute range for every spin
501 : !> \param homo occupied MOs per spin
502 : !> \param para_env the ranks that will hold the matrices
503 : !> \param unit_nr output unit, positive on the writing rank only
504 : ! **************************************************************************************************
505 54 : SUBROUTINE bse_fulldiag_memory_check(bse_env, window, homo, para_env, unit_nr)
506 : TYPE(bse_env_type), INTENT(INOUT) :: bse_env
507 : TYPE(mo_window_type), INTENT(IN) :: window
508 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
509 : TYPE(mp_para_env_type), POINTER :: para_env
510 : INTEGER, INTENT(IN) :: unit_nr
511 :
512 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_fulldiag_memory_check'
513 :
514 : CHARACTER(LEN=16) :: left_str, mem_str
515 : CHARACTER(LEN=80) :: remedy
516 : INTEGER :: handle, ispin, n_ov_joint, nspins
517 : REAL(KIND=dp) :: mem_est_GB, mem_left_GB
518 :
519 54 : CALL timeset(routineN, handle)
520 :
521 54 : nspins = SIZE(homo)
522 54 : n_ov_joint = 0
523 122 : DO ispin = 1, nspins
524 122 : n_ov_joint = n_ov_joint + (homo(ispin) - window%first_mo + 1)*(window%last_mo - homo(ispin))
525 : END DO
526 : mem_est_GB = fulldiag_memory_estimate_GB(n_ov_joint, bse_env%flag_tda /= bse_tda)/ &
527 54 : REAL(para_env%num_pe, KIND=dp)
528 54 : CALL bse_resolve_memory_per_proc(bse_env%memory_per_proc_gb, para_env, unit_nr)
529 54 : mem_left_GB = bse_memory_left_GB(bse_env%memory_per_proc_gb, para_env)
530 :
531 54 : IF (unit_nr > 0) THEN
532 27 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
533 27 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Early memory check for the full diagonalization'
534 27 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Transition pairs in the active window (DFT axis)', &
535 54 : n_ov_joint
536 27 : WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', 'Peak memory estimate per MPI rank (GB)', mem_est_GB
537 27 : WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', 'Memory per MPI rank for the BSE, MEMORY_PER_PROC (GB)', &
538 54 : bse_env%memory_per_proc_gb
539 27 : WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', 'Left after the memory the rank holds (GB)', mem_left_GB
540 27 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
541 : END IF
542 :
543 54 : IF (mem_est_GB > mem_left_GB) THEN
544 0 : WRITE (mem_str, '(F16.3)') mem_est_GB
545 0 : WRITE (left_str, '(F16.3)') mem_left_GB
546 : ! ITERDIAG is closed-shell only, so an open-shell run is not pointed to it
547 0 : IF (nspins == 1) THEN
548 0 : remedy = "more MPI ranks, a tighter window, ITERDIAG, a higher MEMORY_PER_PROC"
549 : ELSE
550 0 : remedy = "more MPI ranks, a tighter window, a higher MEMORY_PER_PROC"
551 : END IF
552 : CALL cp_abort(__LOCATION__, &
553 : "BSE full diagonalization: the estimate of "//TRIM(ADJUSTL(mem_str))// &
554 : " GB per MPI rank exceeds the "//TRIM(ADJUSTL(left_str))// &
555 0 : " GB left of MEMORY_PER_PROC. Try one of: "//TRIM(remedy)//".")
556 : END IF
557 :
558 54 : CALL timestop(handle)
559 :
560 54 : END SUBROUTINE bse_fulldiag_memory_check
561 :
562 : ! **************************************************************************************************
563 : !> \brief bse_fulldiag_memory_check for a caller that holds MO sets: the window on the DFT axis
564 : !> from their eigenvalues, the occupied MOs per spin from their homo
565 : !> \param bse_env settings; the window cutoffs, the TDA flag and MEMORY_PER_PROC, resolved when still -1
566 : !> \param mos the MO sets the GW will correct, for their DFT eigenvalues and occupations
567 : !> \param para_env the ranks that will hold the matrices
568 : !> \param unit_nr output unit, positive on the writing rank only
569 : ! **************************************************************************************************
570 50 : SUBROUTINE bse_fulldiag_memory_check_from_mos(bse_env, mos, para_env, unit_nr)
571 :
572 : TYPE(bse_env_type), INTENT(INOUT) :: bse_env
573 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
574 : TYPE(mp_para_env_type), POINTER :: para_env
575 : INTEGER, INTENT(IN) :: unit_nr
576 :
577 : INTEGER :: ispin
578 : INTEGER, ALLOCATABLE, DIMENSION(:) :: homo
579 : TYPE(mo_window_type) :: window
580 :
581 : ! the preparer warns about a combined window; the estimate stays silent
582 50 : CALL bse_window_from_mos(mos, bse_env%bse_cutoff_occ, bse_env%bse_cutoff_empty, window, warn=.FALSE.)
583 150 : ALLOCATE (homo(SIZE(mos)))
584 112 : DO ispin = 1, SIZE(mos)
585 112 : CALL get_mo_set(mos(ispin), homo=homo(ispin))
586 : END DO
587 50 : CALL bse_fulldiag_memory_check(bse_env, window, homo, para_env, unit_nr)
588 50 : DEALLOCATE (homo)
589 :
590 50 : END SUBROUTINE bse_fulldiag_memory_check_from_mos
591 :
592 : END MODULE bse_main
|
