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