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 GW using RI-RS Approximation for molecules
10 : !> \par History
11 : !> 04.2026 created [Ritaj Tyagi]
12 : ! **************************************************************************************************
13 :
14 : MODULE gw_ri_rs_non_periodic
15 : USE cp_dbcsr_api, ONLY: dbcsr_get_info
16 : USE cp_fm_diag, ONLY: cp_fm_geeig
17 : USE cp_fm_types, ONLY: cp_fm_create,&
18 : cp_fm_get_diag,&
19 : cp_fm_release,&
20 : cp_fm_to_fm,&
21 : cp_fm_type
22 : USE gw_compute_V_and_W, ONLY: compute_W
23 : USE gw_ri_rs_chi_sigma_kernels, ONLY: compute_Sigma_c,&
24 : compute_Sigma_x,&
25 : get_mat_chi_Gamma_tau
26 : USE gw_tensor_large_cell_Gamma, ONLY: delete_unnecessary_files
27 : USE gw_utils, ONLY: analyt_conti_and_print,&
28 : de_init_bs_env
29 : USE gw_utils_fm, ONLY: fm_contract_ABA
30 : USE input_constants, ONLY: G0W0,&
31 : evGW0
32 : USE kinds, ONLY: dp,&
33 : int_8,&
34 : max_line_length
35 : USE message_passing, ONLY: mp_mem_avail_per_rank_GB,&
36 : mp_print_mem_per_rank
37 : USE physcon, ONLY: evolt
38 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
39 : USE post_scf_bandstructure_utils, ONLY: get_all_VBM_CBM_bandgaps
40 : USE qs_environment_types, ONLY: qs_environment_type
41 : USE ri_rs_matrices, ONLY: mask_grid_blocks_near_panel,&
42 : panel_mem_estimate_GB,&
43 : panel_template_elems,&
44 : resolve_grid_panels
45 : USE ri_rs_setup, ONLY: ensure_ri_rs_representation
46 : USE time_frequency_grids, ONLY: negative_time,&
47 : positive_time,&
48 : transform_time_to_frequency
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_ri_rs_non_periodic'
56 :
57 : PUBLIC :: gw_calc_ri_rs_non_periodic
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief GW calculation using RI-RS formalism for molecules
63 : !> \param qs_env ...
64 : !> \param bs_env Band-structure environment containing GW parameters.
65 : ! **************************************************************************************************
66 40 : SUBROUTINE gw_calc_ri_rs_non_periodic(qs_env, bs_env)
67 :
68 : TYPE(qs_environment_type), POINTER :: qs_env
69 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
70 :
71 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_ri_rs_non_periodic'
72 :
73 : INTEGER :: handle
74 40 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma, fm_W_time
75 :
76 40 : CALL timeset(routineN, handle)
77 :
78 : ! Build the common grid, phi and Z representation with the system's boundary conditions.
79 40 : CALL ensure_ri_rs_representation(qs_env, bs_env)
80 40 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
81 40 : CALL print_ri_rs_memory_estimate(bs_env)
82 :
83 : CALL mp_print_mem_per_rank(bs_env%para_env, bs_env%unit_nr, &
84 40 : label='Memory per MPI process after computing Z_lP:')
85 :
86 : ! ==============================================================================
87 : ! 4. Polarizability matrix χ on the imaginary-time grid
88 : ! G^occ_µλ(i|τ|) = Σ_n^occ C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
89 : ! G^vir_µλ(i|τ|) = Σ_n^vir C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
90 : ! G^occ_ll'(i|τ|) = Σ_µν ϕ_µ(r_l) G^occ_µν ϕ_ν(r_l') (G^vir analogous)
91 : ! χ_ll'(iτ) = G^occ_ll'(i|τ|) ∘ G^vir_ll'(i|τ|) (element-wise)
92 : ! χ_PQ(iτ) = Σ_ll' Z_lP χ_ll'(iτ) Z_l'Q
93 : ! ==============================================================================
94 40 : CALL get_mat_chi_Gamma_tau(bs_env, bs_env%mat_chi_Gamma_tau, bs_env%ri_rs%cutoff_radius_v_w)
95 :
96 : ! ==============================================================================
97 : ! 5. Screened Coulomb interaction W (RI basis)
98 : ! χ_PQ(iτ) -> χ_PQ(iω) -> ε_PQ(iω) -> W_PQ(iω) -> W_PQ(iτ)
99 : ! ==============================================================================
100 40 : CALL compute_W(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_time)
101 :
102 : ! ==============================================================================
103 : ! 6. Exact-exchange self-energy Σ^x
104 : ! D_µν = Σ_n^occ C_µn C_νn (density matrix)
105 : ! D_ll' = Σ_µν ϕ_µ(r_l) D_µν ϕ_ν(r_l')
106 : ! V^tr_ll' = Σ_PQ Z_lP V^tr_PQ Z_l'Q (truncated Coulomb)
107 : ! Σ^x_ll' = D_ll' ∘ V^tr_ll'
108 : ! Σ^x_λσ(k=0) = -Σ_ll' ϕ_λ(r_l) Σ^x_ll' ϕ_σ(r_l')
109 : ! ==============================================================================
110 40 : CALL compute_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma, bs_env%ri_rs%cutoff_radius_v_w)
111 :
112 : ! ==============================================================================
113 : ! 7. Correlation self-energy Σ^c and quasiparticle energies, iterated
114 : ! until eigenvalue self-consistency if the &EVGW0 section is given and
115 : ! done in a single pass for G0W0.
116 : !
117 : ! (a) W_ll'(iτ) = Σ_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
118 : ! Σ^c_ll'(iτ) = -G^occ_ll'(i|τ|) ∘ W_ll'(iτ), τ < 0
119 : ! Σ^c_ll'(iτ) = G^vir_ll'(i|τ|) ∘ W_ll'(iτ), τ > 0
120 : ! Σ^c_λσ(iτ) = Σ_ll' ϕ_λ(r_l) Σ^c_ll'(iτ) ϕ_σ(r_l')
121 : ! (b) Σ^c_λσ(iτ) -> Σ^c_nn(ϵ)
122 : ! ϵ_n^GW,(i) = ϵ_n^DFT + Σ^c_nn[G^(i-1),W](ϵ_n^GW,(i)) + Σ^x_nn - v^xc_nn
123 : ! ==============================================================================
124 40 : CALL compute_Sigma_c_and_QP_energies(bs_env, fm_W_time, fm_Sigma_x_Gamma)
125 :
126 40 : CALL de_init_bs_env(qs_env, bs_env)
127 :
128 40 : CALL timestop(handle)
129 :
130 40 : END SUBROUTINE gw_calc_ri_rs_non_periodic
131 :
132 : ! **************************************************************************************************
133 : !> \brief Correlation self-energy and quasiparticle energies, iterated to eigenvalue
134 : !> self-consistency in G (evGW0).
135 : !>
136 : !> For G0W0 this runs once with G^(0) built from the DFT eigenvalues. For evGW0 the
137 : !> Green's function is rebuilt from the quasiparticle energies of all states until the
138 : !> quasiparticle HOMO, LUMO and HOMO-LUMO gap change by less than EPS_ITER between two
139 : !> cycles, or MAX_ITER cycles are spent. W, Σ^x and everything computed before this
140 : !> routine stay frozen: they either do not depend on the eigenvalues at all (Σ^x is G
141 : !> at τ = 0, a pure density matrix) or are held fixed by construction in GW0. The
142 : !> imaginary-time grid is fixed as well, since W(iτ) lives on it.
143 : !>
144 : !> \param bs_env ...
145 : !> \param fm_W_time ...
146 : !> \param fm_Sigma_x_Gamma ...
147 : ! **************************************************************************************************
148 40 : SUBROUTINE compute_Sigma_c_and_QP_energies(bs_env, fm_W_time, fm_Sigma_x_Gamma)
149 :
150 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
151 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time, fm_Sigma_x_Gamma
152 :
153 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c_and_QP_energies'
154 :
155 : INTEGER :: handle, i_iter, n_iter
156 : LOGICAL :: converged
157 40 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenval_scf_Gamma_dft
158 : REAL(KIND=dp), DIMENSION(2) :: e_fermi_dft
159 : REAL(KIND=dp), DIMENSION(3, 2) :: band_prev
160 40 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
161 :
162 40 : CALL timeset(routineN, handle)
163 :
164 : converged = .FALSE.
165 :
166 40 : CALL init_evGW0(bs_env, n_iter, band_prev, eigenval_scf_Gamma_dft, e_fermi_dft)
167 :
168 : ! evGW0 self-consistency loop; for G0W0, loop is terminated after one iteration
169 54 : DO i_iter = 1, n_iter
170 :
171 54 : bs_env%evgw0_i_iter = i_iter
172 :
173 : ! W_ll'(iτ) = Σ_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
174 : ! Σ^c_ll'(iτ) = -G^occ_ll'(i|τ|) ∘ W_ll'(iτ), τ < 0
175 : ! Σ^c_ll'(iτ) = G^vir_ll'(i|τ|) ∘ W_ll'(iτ), τ > 0
176 : ! Σ^c_λσ(iτ) = Σ_ll' ϕ_λ(r_l) Σ^c_ll'(iτ) ϕ_σ(r_l')
177 : CALL compute_Sigma_c(bs_env, fm_W_time, fm_Sigma_c_Gamma_time, &
178 54 : bs_env%ri_rs%cutoff_radius_v_w, keep_w_and_files=.TRUE.)
179 :
180 : ! Σ^c_λσ(iτ) -> Σ^c_nn(ϵ)
181 : ! ϵ_n^GW,(i) = ϵ_n^DFT + Σ^c_nn[G^(i-1),W](ϵ_n^GW,(i)) + Σ^x_nn - v^xc_nn
182 54 : CALL compute_QP_energies(bs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
183 :
184 54 : IF (bs_env%gw_flavour == G0W0) EXIT
185 :
186 18 : CALL print_evGW0_band_edges(bs_env, band_prev, i_iter, n_iter, converged)
187 :
188 18 : IF (i_iter == n_iter .OR. converged) EXIT
189 :
190 : ! eigenvalues ϵ_n to be updated in G:
191 : ! G^occ_µλ(i|τ|) = Σ_n^occ C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
192 : ! G^vir_µλ(i|τ|) = Σ_n^vir C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
193 72 : CALL update_eigenvalues_G(bs_env)
194 :
195 : END DO
196 :
197 40 : CALL cp_fm_release(fm_W_time)
198 :
199 40 : CALL reset_and_clean_bs_env(bs_env, eigenval_scf_Gamma_dft, e_fermi_dft, fm_Sigma_x_Gamma)
200 :
201 40 : CALL delete_unnecessary_files(bs_env)
202 :
203 40 : CALL timestop(handle)
204 :
205 40 : END SUBROUTINE compute_Sigma_c_and_QP_energies
206 :
207 : ! **************************************************************************************************
208 : !> \brief Sets up the evGW0 eigenvalue self-consistency loop: the cycle count, the eigenvalues
209 : !> the first Green's function is built from, and the DFT reference that the loop
210 : !> overwrites. A G0W0 run reduces to a single cycle and needs none of it.
211 : !> \param bs_env ...
212 : !> \param n_iter ...
213 : !> \param band_prev ...
214 : !> \param eigenval_scf_Gamma_dft ...
215 : !> \param e_fermi_dft ...
216 : ! **************************************************************************************************
217 40 : SUBROUTINE init_evGW0(bs_env, n_iter, band_prev, eigenval_scf_Gamma_dft, e_fermi_dft)
218 :
219 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
220 : INTEGER, INTENT(OUT) :: n_iter
221 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(OUT) :: band_prev
222 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
223 : INTENT(OUT) :: eigenval_scf_Gamma_dft
224 : REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: e_fermi_dft
225 :
226 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_evGW0'
227 :
228 : INTEGER :: handle
229 :
230 40 : CALL timeset(routineN, handle)
231 :
232 40 : n_iter = 1
233 40 : band_prev(:, :) = 0.0_dp
234 40 : e_fermi_dft(:) = 0.0_dp
235 :
236 40 : IF (bs_env%gw_flavour == evGW0) THEN
237 4 : n_iter = bs_env%evgw0_iter
238 :
239 : ! eigenvalues currently in G; the first cycle starts is a plain G0W0 step
240 70 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_scf(:, :, :)
241 : ! the loop overwrites these; post-GW printing expects the DFT values back
242 76 : ALLOCATE (eigenval_scf_Gamma_dft, SOURCE=bs_env%eigenval_scf_Gamma)
243 12 : e_fermi_dft(:) = bs_env%e_fermi(:)
244 : END IF
245 :
246 40 : CALL timestop(handle)
247 :
248 40 : END SUBROUTINE init_evGW0
249 :
250 : ! **************************************************************************************************
251 : !> \brief Feeds the quasiparticle energies of the current evGW0 cycle back into the Green's
252 : !> function used by the next one, and re-centres the Fermi level between the new
253 : !> quasiparticle HOMO and LUMO.
254 : !> \param bs_env ...
255 : ! **************************************************************************************************
256 14 : SUBROUTINE update_eigenvalues_G(bs_env)
257 :
258 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
259 :
260 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_eigenvalues_G'
261 :
262 : INTEGER :: handle, i_mo, ispin
263 :
264 14 : CALL timeset(routineN, handle)
265 :
266 : ! record this cycle's evGW0 result; it is also what the next Green's function is built from
267 230 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_GW(:, :, :)
268 :
269 34 : DO ispin = 1, bs_env%n_spin
270 : ! Update all physical states; exclude linear-dependency placeholders.
271 196 : DO i_mo = 1, bs_env%n_mo_retained
272 196 : bs_env%eigenval_scf_Gamma(i_mo, ispin) = bs_env%eigenval_GW(i_mo, 1, ispin)
273 : END DO
274 : bs_env%e_fermi(ispin) = &
275 : 0.5_dp*(bs_env%eigenval_GW(bs_env%n_occ(ispin), 1, ispin) + &
276 34 : bs_env%eigenval_GW(bs_env%n_occ(ispin) + 1, 1, ispin))
277 : END DO
278 :
279 14 : CALL timestop(handle)
280 :
281 14 : END SUBROUTINE update_eigenvalues_G
282 :
283 : ! **************************************************************************************************
284 : !> \brief Restores the DFT reference that the evGW0 loop overwrote and cleanup
285 : !> \param bs_env ...
286 : !> \param eigenval_scf_Gamma_dft ...
287 : !> \param e_fermi_dft ...
288 : !> \param fm_Sigma_x_Gamma ..
289 : ! **************************************************************************************************
290 40 : SUBROUTINE reset_and_clean_bs_env(bs_env, eigenval_scf_Gamma_dft, e_fermi_dft, fm_Sigma_x_Gamma)
291 :
292 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
293 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
294 : INTENT(INOUT) :: eigenval_scf_Gamma_dft
295 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: e_fermi_dft
296 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
297 :
298 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_and_clean_bs_env'
299 :
300 : INTEGER :: handle
301 :
302 40 : CALL timeset(routineN, handle)
303 :
304 40 : IF (bs_env%gw_flavour == evGW0) THEN
305 : ! Update array with final evGW0 result
306 70 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_GW(:, :, :)
307 : ! put the DFT reference back for the post-GW DOS/band-edge printing
308 64 : bs_env%eigenval_scf_Gamma(:, :) = eigenval_scf_Gamma_dft(:, :)
309 12 : bs_env%e_fermi(:) = e_fermi_dft(:)
310 4 : DEALLOCATE (eigenval_scf_Gamma_dft)
311 4 : CALL cp_fm_release(fm_Sigma_x_Gamma)
312 : END IF
313 :
314 40 : CALL timestop(handle)
315 :
316 40 : END SUBROUTINE reset_and_clean_bs_env
317 :
318 : ! **************************************************************************************************
319 : !> \brief Print the quasiparticle HOMO, LUMO and HOMO-LUMO gap of the current evGW0 cycle
320 : !> \param bs_env ...
321 : !> \param band_prev ...
322 : !> \param i_iter ...
323 : !> \param n_iter ...
324 : !> \param converged ...
325 : ! **************************************************************************************************
326 18 : SUBROUTINE print_evGW0_band_edges(bs_env, band_prev, i_iter, n_iter, converged)
327 :
328 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
329 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(INOUT) :: band_prev
330 : INTEGER, INTENT(IN) :: i_iter, n_iter
331 : LOGICAL, INTENT(OUT) :: converged
332 :
333 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_evGW0_band_edges'
334 :
335 : INTEGER :: handle, homo, ispin, u
336 : REAL(KIND=dp) :: max_delta
337 : REAL(KIND=dp), DIMENSION(3) :: band
338 :
339 18 : CALL timeset(routineN, handle)
340 :
341 18 : u = bs_env%unit_nr
342 :
343 18 : converged = (i_iter > 1)
344 18 : max_delta = 0.0_dp
345 :
346 18 : IF (u > 0) THEN
347 9 : WRITE (u, '(A)') ' '
348 9 : WRITE (u, '(T2,A)') REPEAT('-', 79)
349 9 : WRITE (u, '(T2,A,I4,A,I4)') 'evGW0 cycle', i_iter, ' /', n_iter
350 9 : WRITE (u, '(T2,A)') REPEAT('-', 79)
351 : END IF
352 :
353 44 : DO ispin = 1, bs_env%n_spin
354 :
355 26 : homo = bs_env%n_occ(ispin)
356 26 : band(1) = bs_env%eigenval_GW(homo, 1, ispin)
357 26 : band(2) = bs_env%eigenval_GW(homo + 1, 1, ispin)
358 26 : band(3) = band(2) - band(1)
359 :
360 26 : IF (i_iter > 1) THEN
361 80 : max_delta = MAX(max_delta, MAXVAL(ABS(band(:) - band_prev(:, ispin))))
362 : END IF
363 :
364 26 : IF (u > 0) THEN
365 13 : IF (bs_env%n_spin == 2) WRITE (u, '(T2,A,I0)') 'Spin ', ispin
366 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 HOMO (eV)', band(1)*evolt
367 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 LUMO (eV)', band(2)*evolt
368 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 HOMO-LUMO gap (eV)', band(3)*evolt
369 : END IF
370 :
371 122 : band_prev(:, ispin) = band(:)
372 :
373 : END DO
374 :
375 18 : IF (i_iter > 1) THEN
376 14 : converged = (max_delta < bs_env%evgw0_eps_iter)
377 21 : IF (u > 0) WRITE (u, '(T2,A,T61,F20.6)') 'Max. change to previous cycle (eV)', &
378 14 : max_delta*evolt
379 : END IF
380 :
381 11 : IF (u > 0) WRITE (u, '(T2,A)') REPEAT('-', 79)
382 :
383 18 : IF (converged) THEN
384 4 : IF (u > 0) THEN
385 2 : WRITE (u, '(A)') ' '
386 : WRITE (u, '(T2,A,I4,A)') &
387 2 : 'evGW0 eigenvalue self-consistency reached in', i_iter, ' cycles.'
388 2 : WRITE (u, '(A)') ' '
389 : END IF
390 14 : ELSE IF (i_iter == n_iter) THEN
391 : CALL cp_warn(__LOCATION__, &
392 : "The evGW0 eigenvalue self-consistency cycle did not converge "// &
393 : "within MAX_ITER cycles. The reported quasiparticle energies are "// &
394 0 : "those of the last cycle.")
395 : END IF
396 :
397 18 : CALL timestop(handle)
398 :
399 18 : END SUBROUTINE print_evGW0_band_edges
400 :
401 : ! **************************************************************************************************
402 : !> \brief Estimates and prints per-process memory requirements for the RI-RS GW calculation.
403 : !> \param bs_env ...
404 : ! **************************************************************************************************
405 40 : SUBROUTINE print_ri_rs_memory_estimate(bs_env)
406 :
407 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
408 :
409 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_ri_rs_memory_estimate'
410 :
411 : CHARACTER(LEN=max_line_length) :: msg
412 : INTEGER :: handle, ipan, n_grid_total, n_procs, &
413 : n_RI, pan_rows, pan_width
414 : INTEGER(KIND=int_8) :: nze_tmpl
415 40 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
416 40 : INTEGER, DIMENSION(:), POINTER :: r_blk_sizes
417 : LOGICAL :: use_cutoff
418 40 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: grid_used
419 : REAL(KIND=dp) :: mem_avail_GB, mem_pan_GB, mem_panels_GB, &
420 : mem_Z_lP_GB
421 :
422 40 : CALL timeset(routineN, handle)
423 :
424 40 : n_grid_total = bs_env%ri_rs%n_grid_points
425 40 : CALL dbcsr_get_info(bs_env%ri_rs%mat_phi_mu_l, row_blk_size=r_blk_sizes)
426 570 : CPASSERT(SUM(r_blk_sizes) == n_grid_total)
427 40 : n_RI = bs_env%n_RI
428 40 : n_procs = bs_env%para_env%num_pe
429 :
430 : ! Z_lP upper bound: dense n_grid × n_RI, distributed evenly across all ranks.
431 : ! The actual sparse Z_lP is smaller due to the per-atom locality cutoff.
432 : mem_Z_lP_GB = REAL(n_grid_total, dp)*REAL(n_RI, dp)*8.0_dp/ &
433 40 : REAL(n_procs, dp)*1.0E-9_dp
434 :
435 : ! Peak panel memory during Σ^c: two G panels (A_occ, A_vir) + one W panel plus the
436 : ! grid × RI / grid × AO intermediates. With the CUTOFF_RADIUS_RL_W restriction the panel
437 : ! matrices only allocate the geo-template blocks, so use the same nze-aware model as the
438 : ! panel planner (panel_mem_estimate_GB); without the cutoff, dense panel_rows × n_grid.
439 : ! Plus the n_RI × n_RI W_aux matrix. Distributed over n_procs ranks.
440 : use_cutoff = bs_env%ri_rs%cutoff_radius_v_w > 0.0_dp .AND. &
441 40 : ALLOCATED(bs_env%ri_rs%chunk_centroids)
442 : CALL resolve_grid_panels(bs_env%ri_rs, bs_env%para_env, bs_env%dbcsr_msg_elem_limit, &
443 40 : pan_first, pan_last, bs_env%ri_rs%cutoff_radius_v_w)
444 40 : mem_panels_GB = 0.0_dp
445 80 : DO ipan = 1, SIZE(pan_first)
446 570 : pan_rows = SUM(r_blk_sizes(pan_first(ipan):pan_last(ipan)))
447 40 : IF (use_cutoff) THEN
448 : CALL mask_grid_blocks_near_panel(bs_env%ri_rs%chunk_centroids, pan_first(ipan), &
449 : pan_last(ipan), bs_env%ri_rs%cutoff_radius_v_w, &
450 2 : grid_used)
451 44 : pan_width = SUM(r_blk_sizes, MASK=grid_used)
452 : CALL panel_template_elems(r_blk_sizes, bs_env%ri_rs%chunk_centroids, &
453 : grid_used, pan_first(ipan), pan_last(ipan), &
454 2 : bs_env%ri_rs%cutoff_radius_v_w, nze_tmpl)
455 : ELSE
456 38 : pan_width = n_grid_total
457 38 : nze_tmpl = INT(pan_rows, int_8)*INT(pan_width, int_8)
458 : END IF
459 : CALL panel_mem_estimate_GB(nze_tmpl, pan_rows, pan_width, n_grid_total, &
460 40 : n_RI, bs_env%n_ao, n_procs, mem_pan_GB)
461 80 : mem_panels_GB = MAX(mem_panels_GB, mem_pan_GB)
462 : END DO
463 : mem_panels_GB = mem_panels_GB + &
464 40 : REAL(n_RI, dp)*REAL(n_RI, dp)*8.0_dp/REAL(n_procs, dp)*1.0E-9_dp
465 :
466 : ! Available memory per process = node MemLikelyFree / ranks-per-node, min across ranks
467 : ! (0 on non-Linux => warnings suppressed below).
468 40 : CALL mp_mem_avail_per_rank_GB(bs_env%para_env, mem_avail_GB)
469 :
470 40 : IF (bs_env%unit_nr > 0) THEN
471 20 : WRITE (bs_env%unit_nr, '(A)') ' '
472 20 : WRITE (bs_env%unit_nr, '(T2,A)') 'RI-RS memory estimate per MPI process:'
473 : WRITE (bs_env%unit_nr, '(T4,A,F37.2,A)') &
474 20 : 'Available memory per process (system)', mem_avail_GB, ' GB'
475 : WRITE (bs_env%unit_nr, '(T4,A,F18.2,A)') &
476 20 : 'Required for Z_lP (dense upper bound; actual is sparser)', mem_Z_lP_GB, ' GB'
477 : WRITE (bs_env%unit_nr, '(T4,A,F25.2,A)') &
478 20 : 'Required for χ, W, Σ panels (peak per panel step)', mem_panels_GB, ' GB'
479 20 : IF (bs_env%ri_rs%Z_lP_exists .AND. .NOT. bs_env%auto_ri%enabled) THEN
480 1 : WRITE (bs_env%unit_nr, '(T4,A)') 'Z_lP fitting skipped (matrix read from restart)'
481 : END IF
482 :
483 : END IF
484 :
485 40 : IF (mem_avail_GB > 0.0_dp .AND. mem_Z_lP_GB > mem_avail_GB) THEN
486 : WRITE (msg, '(A,F0.2,A,F0.2,A)') &
487 0 : "The estimated memory for Z_lP, ", mem_Z_lP_GB, " GB per process, exceeds the "// &
488 0 : "available ", mem_avail_GB, " GB. Z_lP (n_grid x n_RI) is distributed across all "// &
489 0 : "MPI ranks, so add nodes or use fewer MPI ranks per node."
490 0 : CPWARN(TRIM(msg))
491 : END IF
492 :
493 40 : IF (mem_avail_GB > 0.0_dp .AND. mem_panels_GB > mem_avail_GB) THEN
494 : WRITE (msg, '(A,F0.2,A,F0.2,A)') &
495 0 : "The estimated peak memory of the chi/W/Sigma panels, ", mem_panels_GB, &
496 0 : " GB per process, exceeds the available ", mem_avail_GB, " GB. Panel memory "// &
497 : "scales roughly as 3*panel_size*n_grid/n_procs, so add nodes, use fewer MPI ranks "// &
498 0 : "per node."
499 0 : CPWARN(TRIM(msg))
500 : END IF
501 :
502 40 : CALL timestop(handle)
503 :
504 80 : END SUBROUTINE print_ri_rs_memory_estimate
505 :
506 : ! **************************************************************************************************
507 : !> \brief Σ^c_λσ(iτ) -> Σ^c_nn(ϵ) and the quasi-particle levels of the non-periodic RI-RS path,
508 : !> ϵ_n^GW = ϵ_n^DFT + Σ^c_nn(ϵ_n^GW) + Σ^x_nn - v^xc_nn.
509 : !> \param bs_env ...
510 : !> \param fm_Sigma_x_Gamma ...
511 : !> \param fm_Sigma_c_Gamma_time ...
512 : ! **************************************************************************************************
513 54 : SUBROUTINE compute_QP_energies(bs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
514 :
515 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
516 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
517 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
518 :
519 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies'
520 :
521 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Sigma_c_n_freq
522 : INTEGER :: handle, ispin, j_t
523 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Sigma_x_n, V_xc_n
524 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Sigma_c_n_time
525 : TYPE(cp_fm_type) :: fm_ks, fm_mos, fm_s, fm_work
526 :
527 54 : CALL timeset(routineN, handle)
528 :
529 54 : CALL cp_fm_create(fm_ks, bs_env%fm_s_Gamma%matrix_struct)
530 54 : CALL cp_fm_create(fm_s, bs_env%fm_s_Gamma%matrix_struct)
531 54 : CALL cp_fm_create(fm_mos, bs_env%fm_s_Gamma%matrix_struct)
532 54 : CALL cp_fm_create(fm_work, bs_env%fm_s_Gamma%matrix_struct)
533 :
534 216 : ALLOCATE (V_xc_n(bs_env%n_ao), Sigma_x_n(bs_env%n_ao))
535 270 : ALLOCATE (Sigma_c_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
536 216 : ALLOCATE (Sigma_c_n_freq(bs_env%n_ao, bs_env%num_freq_points_fit))
537 :
538 120 : DO ispin = 1, bs_env%n_spin
539 :
540 : ! 1. Roothaan-Hall H^KS_µν C_νn = S_µν C_νn ϵ_n
541 66 : CALL cp_fm_to_fm(bs_env%fm_ks_Gamma(ispin), fm_ks)
542 66 : CALL cp_fm_to_fm(bs_env%fm_s_Gamma, fm_s)
543 66 : CALL cp_fm_geeig(fm_ks, fm_s, fm_mos, bs_env%eigenval_scf(:, 1, ispin), fm_work)
544 :
545 : ! 2. v^xc_µν -> v^xc_nn and Σ^x_µν -> Σ^x_nn
546 66 : CALL to_Gamma_and_mo_real(V_xc_n, bs_env%fm_V_xc_Gamma(ispin), fm_mos)
547 66 : CALL to_Gamma_and_mo_real(Sigma_x_n, fm_Sigma_x_Gamma(ispin), fm_mos)
548 :
549 : ! 3. Σ^c_µν(+/-i|τ_j|) -> Σ^c_nn(+/-i|τ_j|)
550 896 : DO j_t = 1, bs_env%num_time_freq_points
551 : CALL to_Gamma_and_mo_real(Sigma_c_n_time(:, j_t, positive_time), &
552 830 : fm_Sigma_c_Gamma_time(j_t, positive_time, ispin), fm_mos)
553 : CALL to_Gamma_and_mo_real(Sigma_c_n_time(:, j_t, negative_time), &
554 896 : fm_Sigma_c_Gamma_time(j_t, negative_time, ispin), fm_mos)
555 : END DO
556 :
557 : ! 4. Σ^c_nn(iτ) -> Σ^c_nn(iω)
558 66 : CALL transform_time_to_frequency(bs_env%time_frequency_grid, Sigma_c_n_time, Sigma_c_n_freq)
559 818 : Sigma_c_n_freq(1:bs_env%n_occ(ispin), :) = CONJG(Sigma_c_n_freq(1:bs_env%n_occ(ispin), :))
560 :
561 : ! 5. Analytic continuation Σ^c_nn(iω) -> Σ^c_nn(ϵ) and the QP levels
562 : CALL analyt_conti_and_print(bs_env, Sigma_c_n_freq, Sigma_x_n, V_xc_n, &
563 120 : bs_env%eigenval_scf(:, 1, ispin), 1, ispin)
564 :
565 : END DO ! ispin
566 :
567 54 : CALL get_all_VBM_CBM_bandgaps(bs_env)
568 :
569 54 : IF (bs_env%gw_flavour == G0W0) CALL cp_fm_release(fm_Sigma_x_Gamma)
570 54 : CALL cp_fm_release(fm_Sigma_c_Gamma_time)
571 :
572 54 : CALL cp_fm_release(fm_ks)
573 54 : CALL cp_fm_release(fm_s)
574 54 : CALL cp_fm_release(fm_mos)
575 54 : CALL cp_fm_release(fm_work)
576 :
577 54 : CALL timestop(handle)
578 :
579 108 : END SUBROUTINE compute_QP_energies
580 :
581 : ! **************************************************************************************************
582 : !> \brief AO -> MO transform of a Γ-point matrix
583 : !> \param array_n ...
584 : !> \param fm_Gamma ...
585 : !> \param fm_mos ...
586 : ! **************************************************************************************************
587 1792 : SUBROUTINE to_Gamma_and_mo_real(array_n, fm_Gamma, fm_mos)
588 :
589 : REAL(KIND=dp), DIMENSION(:) :: array_n
590 : TYPE(cp_fm_type) :: fm_Gamma, fm_mos
591 :
592 : CHARACTER(LEN=*), PARAMETER :: routineN = 'to_Gamma_and_mo_real'
593 :
594 : INTEGER :: handle
595 : TYPE(cp_fm_type) :: fm_mo
596 :
597 1792 : CALL timeset(routineN, handle)
598 :
599 1792 : CALL cp_fm_create(fm_mo, fm_Gamma%matrix_struct)
600 :
601 : ! A_nn' = Σ_μν C_μn A_μν C_νn'
602 1792 : CALL fm_contract_ABA(fm_mos, fm_Gamma, fm_mo)
603 :
604 1792 : CALL cp_fm_get_diag(fm_mo, array_n)
605 :
606 1792 : CALL cp_fm_release(fm_mo)
607 :
608 1792 : CALL timestop(handle)
609 :
610 1792 : END SUBROUTINE to_Gamma_and_mo_real
611 :
612 : END MODULE gw_ri_rs_non_periodic
|
