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
10 : !> \par History
11 : !> 01.2026 Maximilian Graml: add more bounds to exploit sparsity in 3c integrals, fixes
12 : !> \author Jan Wilhelm
13 : !> \date 07.2023
14 : ! **************************************************************************************************
15 : MODULE gw_utils
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind_set
18 : USE basis_set_types, ONLY: get_gto_basis_set,&
19 : gto_basis_set_type
20 : USE bibliography, ONLY: Graml2024,&
21 : cite_reference
22 : USE cell_types, ONLY: cell_type,&
23 : pbc,&
24 : scaled_to_real
25 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
26 : cp_blacs_env_release,&
27 : cp_blacs_env_type
28 : USE cp_cfm_types, ONLY: cp_cfm_create,&
29 : cp_cfm_release,&
30 : cp_cfm_to_cfm,&
31 : cp_cfm_to_fm,&
32 : cp_cfm_type
33 : USE cp_control_types, ONLY: dft_control_type
34 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
35 : dbcsr_p_type,&
36 : dbcsr_set,&
37 : dbcsr_type_symmetric
38 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
39 : copy_fm_to_dbcsr,&
40 : dbcsr_allocate_matrix_set,&
41 : dbcsr_deallocate_matrix_set
42 : USE cp_files, ONLY: close_file,&
43 : open_file
44 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
45 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
46 : cp_fm_struct_release,&
47 : cp_fm_struct_type
48 : USE cp_fm_types, ONLY: cp_fm_create,&
49 : cp_fm_get_diag,&
50 : cp_fm_release,&
51 : cp_fm_set_all,&
52 : cp_fm_to_fm,&
53 : cp_fm_type
54 : USE cp_log_handling, ONLY: cp_get_default_logger,&
55 : cp_logger_type
56 : USE cp_output_handling, ONLY: cp_print_key_generate_filename
57 : USE dbt_api, ONLY: &
58 : dbt_clear, dbt_create, dbt_destroy, dbt_filter, dbt_iterator_blocks_left, &
59 : dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, &
60 : dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
61 : USE gw_auto_ri_generate_RI_basis, ONLY: generate_auto_ri_basis
62 : USE gw_auto_ri_types, ONLY: auto_ri_type
63 : USE gw_compute_V_and_W, ONLY: get_V_tr_R
64 : USE gw_utils_compute_integrals, ONLY: build_3c_integral_block
65 : USE gw_utils_fm, ONLY: cfm_contract_ABA,&
66 : fm_invert,&
67 : local_complex_power
68 : USE input_constants, ONLY: &
69 : G0W0, cohsex_qp_update, do_potential_coulomb, do_potential_truncated, evGW0, &
70 : ri_rpa_g0w0_crossing_newton, ri_rs_large_cell_Gamma, ri_rs_non_periodic, &
71 : rtp_bse_kernel_ri_ao, rtp_bse_kernel_ri_rs, rtp_method_bse, rtp_method_bse_linearized, &
72 : tensor_large_cell_Gamma, tensor_small_cell_full_kp, tensor_static_cohsex, xc_none
73 : USE input_section_types, ONLY: section_vals_get,&
74 : section_vals_get_subs_vals,&
75 : section_vals_type,&
76 : section_vals_val_get,&
77 : section_vals_val_set
78 : USE kinds, ONLY: default_path_length,&
79 : default_string_length,&
80 : dp,&
81 : int_8
82 : USE kpoint_k_r_trafo_simple, ONLY: rs_to_kp
83 : USE kpoint_types, ONLY: get_kpoint_info,&
84 : kpoint_create,&
85 : kpoint_type
86 : USE libint_wrapper, ONLY: cp_libint_static_cleanup,&
87 : cp_libint_static_init
88 : USE machine, ONLY: m_walltime
89 : USE mathlib, ONLY: gcd
90 : USE message_passing, ONLY: mp_cart_type,&
91 : mp_mem_avail_per_rank_GB,&
92 : mp_mem_used_per_rank_GB,&
93 : mp_para_env_type
94 : USE mp2_gpw, ONLY: create_mat_munu
95 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat,&
96 : trunc_coulomb_for_exchange
97 : USE particle_methods, ONLY: get_particle_set
98 : USE particle_types, ONLY: particle_type
99 : USE physcon, ONLY: angstrom,&
100 : evolt
101 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
102 : USE post_scf_bandstructure_utils, ONLY: allocate_GW_eigenvalues,&
103 : gw_flavour_label,&
104 : rsmat_to_kp
105 : USE qs_energy_types, ONLY: qs_energy_type
106 : USE qs_environment_types, ONLY: get_qs_env,&
107 : qs_env_part_release,&
108 : qs_environment_type
109 : USE qs_integral_utils, ONLY: basis_set_list_setup
110 : USE qs_interactions, ONLY: init_interaction_radii_orb_basis
111 : USE qs_kind_types, ONLY: get_qs_kind,&
112 : qs_kind_type
113 : USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
114 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
115 : USE qs_tensors, ONLY: build_3c_integrals,&
116 : build_3c_neighbor_lists,&
117 : get_tensor_occupancy,&
118 : neighbor_list_3c_destroy
119 : USE qs_tensors_types, ONLY: create_2c_tensor,&
120 : create_3c_tensor,&
121 : distribution_3d_create,&
122 : distribution_3d_type,&
123 : neighbor_list_3c_type
124 : USE ri_rs_types, ONLY: ri_rs_env,&
125 : ri_rs_grid_opt_type
126 : USE rpa_gw, ONLY: continuation_pade
127 : USE rpa_gw_im_time_util, ONLY: build_ao_propagator
128 : USE time_frequency_grids, ONLY: build_minimax_time_frequency_grid
129 : #include "base/base_uses.f90"
130 :
131 : IMPLICIT NONE
132 :
133 : PRIVATE
134 :
135 : PUBLIC :: create_and_init_bs_env_for_gw, &
136 : de_init_bs_env, &
137 : compute_xkp, &
138 : analyt_conti_and_print, &
139 : add_R, &
140 : is_cell_in_index_to_cell, &
141 : rtbse_resolve_rirs_flag, &
142 : G_occ_vir
143 :
144 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils'
145 :
146 : CONTAINS
147 :
148 : ! **************************************************************************************************
149 : !> \brief Initializes the GW environment from the input and electronic-structure data.
150 : !> \param qs_env ...
151 : !> \param bs_env Band-structure environment containing GW parameters.
152 : !> \param bs_sec BAND_STRUCTURE input section
153 : ! **************************************************************************************************
154 138 : SUBROUTINE create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec)
155 : TYPE(qs_environment_type), POINTER :: qs_env
156 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
157 : TYPE(section_vals_type), POINTER :: bs_sec
158 :
159 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_and_init_bs_env_for_gw'
160 :
161 : INTEGER :: handle
162 :
163 138 : CALL timeset(routineN, handle)
164 :
165 138 : CALL cite_reference(Graml2024)
166 :
167 138 : CALL get_parameters_from_qs_env(qs_env, bs_env)
168 :
169 138 : CALL read_gw_input_parameters(bs_env, bs_sec)
170 :
171 138 : CALL print_header_and_input_parameters(bs_env)
172 :
173 138 : CALL setup_AO_and_RI_basis_set(qs_env, bs_env)
174 :
175 138 : CALL set_heuristic_parameters(bs_env, qs_env)
176 :
177 138 : CALL cp_libint_static_init()
178 138 : IF (bs_env%auto_ri%enabled) CALL generate_auto_ri_basis(qs_env, bs_env)
179 :
180 138 : CALL get_RI_basis(qs_env, bs_env)
181 :
182 138 : CALL compute_V_xc(qs_env, bs_env)
183 :
184 138 : CALL init_interaction_radii(bs_env)
185 :
186 : ! The RI-RS drivers build their (µν|P) blocks on the fly through gw_3c_ctx
187 : ! RT-BSE with the AO-RI kernel need nl_3c it after GW, so
188 : ! keep everything when an RT-BSE run follows.
189 : IF (.NOT. bs_env%do_gw_ri_rs .OR. &
190 138 : bs_env%rtp_method == rtp_method_bse .OR. &
191 : bs_env%rtp_method == rtp_method_bse_linearized) THEN
192 112 : CALL create_tensors(qs_env, bs_env)
193 : END IF
194 :
195 138 : CALL allocate_GW_eigenvalues(bs_env)
196 :
197 260 : SELECT CASE (bs_env%gw_implementation)
198 : CASE (tensor_large_cell_Gamma, ri_rs_large_cell_Gamma, ri_rs_non_periodic, tensor_static_cohsex)
199 :
200 122 : IF (.NOT. bs_env%do_gw_ri_rs) THEN
201 82 : CALL check_sparsity_3c(qs_env, bs_env)
202 :
203 82 : CALL set_sparsity_parallelization_parameters(bs_env)
204 : END IF
205 :
206 122 : CALL check_for_restart_files(qs_env, bs_env)
207 :
208 : CASE (tensor_small_cell_full_kp)
209 :
210 16 : CALL compute_3c_integrals(qs_env, bs_env)
211 :
212 16 : CALL setup_cells_Delta_R(bs_env)
213 :
214 16 : CALL setup_parallelization_Delta_R(bs_env)
215 :
216 16 : CALL allocate_matrices_small_cell_full_kp_tensor(qs_env, bs_env)
217 :
218 16 : CALL trafo_V_xc_R_to_kp(qs_env, bs_env)
219 :
220 154 : CALL heuristic_RI_regularization(qs_env, bs_env)
221 :
222 : END SELECT
223 :
224 138 : CALL setup_time_and_frequency_minimax_grid(bs_env)
225 :
226 : ! Free memory in qs_env. SCF real-space grids, task lists, neighbor
227 : ! lists and work arrays on large systems hold GBs per rank for the whole run.
228 : ! Not required by GW. LDOS and RT-BSE still needs these information.
229 : !
230 : ! Recommendation in case of memory issues: first perform GW calculation without calculating
231 : ! LDOS (to save memory). Then, use GW restart files
232 : ! in a subsequent calculation to calculate the LDOS
233 : IF (.NOT. bs_env%do_ldos .AND. &
234 138 : bs_env%rtp_method /= rtp_method_bse .AND. &
235 : bs_env%rtp_method /= rtp_method_bse_linearized) THEN
236 64 : CALL qs_env_part_release(qs_env)
237 : END IF
238 :
239 138 : CALL timestop(handle)
240 :
241 138 : END SUBROUTINE create_and_init_bs_env_for_gw
242 :
243 : ! **************************************************************************************************
244 : !> \brief Store the QS environment parameters needed by GW.
245 : !> \param qs_env ...
246 : !> \param bs_env ...
247 : ! **************************************************************************************************
248 138 : SUBROUTINE get_parameters_from_qs_env(qs_env, bs_env)
249 : TYPE(qs_environment_type), POINTER :: qs_env
250 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
251 :
252 : TYPE(section_vals_type), POINTER :: input, rtbse_sec
253 :
254 138 : NULLIFY (input, rtbse_sec)
255 : CALL get_qs_env(qs_env, atomic_kind_set=bs_env%ri_rs%atomic_kind_set, &
256 138 : cell=bs_env%ri_rs%cell, particle_set=bs_env%ri_rs%particle_set, input=input)
257 138 : rtbse_sec => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%RTBSE")
258 138 : CALL section_vals_val_get(rtbse_sec, "CUTOFF_RADIUS_W0", r_val=bs_env%ri_rs%cutoff_radius_w0)
259 138 : END SUBROUTINE get_parameters_from_qs_env
260 :
261 : ! **************************************************************************************************
262 : !> \brief Reads the GW input parameters.
263 : !> \param bs_env ...
264 : !> \param bs_sec BAND_STRUCTURE input section
265 : ! **************************************************************************************************
266 276 : SUBROUTINE read_gw_input_parameters(bs_env, bs_sec)
267 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
268 : TYPE(section_vals_type), POINTER :: bs_sec
269 :
270 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_gw_input_parameters'
271 :
272 : INTEGER :: handle
273 : LOGICAL :: do_evgw0
274 : TYPE(section_vals_type), POINTER :: auto_ri_sec, cohsex_sec, evgw0_sec, &
275 : grid_opt_sec, gw_sec, ri_rs_sec
276 :
277 138 : CALL timeset(routineN, handle)
278 :
279 138 : NULLIFY (auto_ri_sec, cohsex_sec, evgw0_sec, grid_opt_sec, gw_sec, ri_rs_sec)
280 138 : gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
281 138 : auto_ri_sec => section_vals_get_subs_vals(gw_sec, "AUTO_RI")
282 138 : cohsex_sec => section_vals_get_subs_vals(gw_sec, "COHSEX")
283 138 : evgw0_sec => section_vals_get_subs_vals(gw_sec, "EVGW0")
284 138 : ri_rs_sec => section_vals_get_subs_vals(gw_sec, "RI_RS")
285 138 : grid_opt_sec => section_vals_get_subs_vals(ri_rs_sec, "GRID_OPTIMIZATION")
286 :
287 138 : CALL section_vals_val_get(gw_sec, "NUM_TIME_FREQ_POINTS", i_val=bs_env%num_time_freq_points)
288 138 : CALL section_vals_val_get(gw_sec, "EPS_FILTER", r_val=bs_env%eps_filter)
289 138 : CALL section_vals_val_get(gw_sec, "REGULARIZATION_RI", r_val=bs_env%input_regularization_RI)
290 138 : CALL section_vals_val_get(gw_sec, "REGULARIZATION_MINIMAX", r_val=bs_env%input_regularization_minimax)
291 138 : CALL section_vals_val_get(gw_sec, "CUTOFF_RADIUS_RI", r_val=bs_env%ri_metric%cutoff_radius)
292 138 : CALL section_vals_val_get(gw_sec, "MEMORY_PER_PROC", r_val=bs_env%input_memory_per_proc_GB)
293 138 : CALL section_vals_val_get(gw_sec, "APPROX_KP_EXTRAPOL", l_val=bs_env%approx_kp_extrapol)
294 138 : CALL section_vals_val_get(gw_sec, "SIZE_LATTICE_SUM", i_val=bs_env%size_lattice_sum_V)
295 138 : CALL section_vals_val_get(gw_sec, "KPOINTS_W", i_vals=bs_env%nkp_grid_chi_eps_W_input)
296 138 : CALL section_vals_val_get(gw_sec, "HEDIN_SHIFT", l_val=bs_env%do_hedin_shift)
297 138 : CALL section_vals_val_get(gw_sec, "FREQ_MAX_FIT", r_val=bs_env%freq_max_fit)
298 138 : CALL section_vals_val_get(cohsex_sec, "SCALE_OCC", r_val=bs_env%cohsex_scale_occ)
299 138 : CALL section_vals_val_get(cohsex_sec, "SCALE_VIRT", r_val=bs_env%cohsex_scale_virt)
300 138 : CALL read_rirs_input(bs_env%ri_rs, ri_rs_sec)
301 138 : CALL read_rirs_grid_opt(bs_env%ri_rs, grid_opt_sec)
302 138 : CALL read_evgw0_input(bs_env, evgw0_sec, do_evgw0)
303 138 : CALL read_auto_ri_input(bs_env%auto_ri, auto_ri_sec)
304 138 : CALL check_gw_input(bs_env, do_evgw0)
305 :
306 138 : CALL timestop(handle)
307 :
308 138 : END SUBROUTINE read_gw_input_parameters
309 :
310 : ! **************************************************************************************************
311 : !> \brief Reads the RI-RS input.
312 : !> \param ri_rs RI-RS calculation parameters
313 : !> \param ri_rs_sec RI_RS input section
314 : ! **************************************************************************************************
315 138 : SUBROUTINE read_rirs_input(ri_rs, ri_rs_sec)
316 : TYPE(ri_rs_env), INTENT(INOUT) :: ri_rs
317 : TYPE(section_vals_type), POINTER :: ri_rs_sec
318 :
319 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_rirs_input'
320 :
321 : INTEGER :: handle
322 :
323 138 : CALL timeset(routineN, handle)
324 :
325 138 : CALL section_vals_val_get(ri_rs_sec, "TIKHONOV", r_val=ri_rs%tikhonov)
326 138 : CALL section_vals_val_get(ri_rs_sec, "GRID_SELECT", i_val=ri_rs%grid_select)
327 138 : CALL section_vals_val_get(ri_rs_sec, "CUTOFF_RADIUS_RL_AO", r_val=ri_rs%cutoff_radius_ri_ao)
328 138 : CALL section_vals_val_get(ri_rs_sec, "CUTOFF_RADIUS_RL_W", r_val=ri_rs%cutoff_radius_v_w)
329 :
330 138 : CALL timestop(handle)
331 :
332 138 : END SUBROUTINE read_rirs_input
333 :
334 : ! **************************************************************************************************
335 : !> \brief Read the RI-RS grid optimization input.
336 : !> \param ri_rs RI-RS parameters receiving the optimized-grid settings
337 : !> \param grid_opt_sec GRID_OPTIMIZATION input section
338 : ! **************************************************************************************************
339 138 : SUBROUTINE read_rirs_grid_opt(ri_rs, grid_opt_sec)
340 : TYPE(ri_rs_env), INTENT(INOUT), TARGET :: ri_rs
341 : TYPE(section_vals_type), POINTER :: grid_opt_sec
342 :
343 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_rirs_grid_opt'
344 :
345 : INTEGER :: handle
346 : TYPE(ri_rs_grid_opt_type), POINTER :: grid_opt
347 :
348 138 : CALL timeset(routineN, handle)
349 :
350 138 : CALL section_vals_get(grid_opt_sec, explicit=ri_rs%grid_opt%enabled)
351 138 : IF (.NOT. ri_rs%grid_opt%enabled) THEN
352 128 : CALL timestop(handle)
353 128 : RETURN
354 : END IF
355 :
356 : NULLIFY (grid_opt)
357 10 : grid_opt => ri_rs%grid_opt
358 :
359 10 : CALL section_vals_val_get(grid_opt_sec, "CUTOFF_ATOMIC_CLUSTER", r_val=grid_opt%cutoff_atomic_cluster)
360 10 : CALL section_vals_val_get(grid_opt_sec, "MAX_ITER", i_val=grid_opt%max_iter)
361 10 : CALL section_vals_val_get(grid_opt_sec, "RS_AO_RATIO", r_val=grid_opt%rs_ao_ratio)
362 10 : IF (grid_opt%rs_ao_ratio <= 0.0_dp) CPABORT("RS_AO_RATIO must be positive")
363 10 : CALL timestop(handle)
364 :
365 : END SUBROUTINE read_rirs_grid_opt
366 :
367 : ! **************************************************************************************************
368 : !> \brief Reads the eigenvalue-self-consistent GW input.
369 : !> \param bs_env GW environment receiving the evGW0 convergence settings
370 : !> \param evgw0_sec EVGW0 input section
371 : !> \param do_evgw0 whether the EVGW0 section is activated
372 : ! **************************************************************************************************
373 276 : SUBROUTINE read_evgw0_input(bs_env, evgw0_sec, do_evgw0)
374 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
375 : TYPE(section_vals_type), POINTER :: evgw0_sec
376 : LOGICAL, INTENT(OUT) :: do_evgw0
377 :
378 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_evgw0_input'
379 :
380 : INTEGER :: handle
381 :
382 138 : CALL timeset(routineN, handle)
383 :
384 138 : CALL section_vals_val_get(evgw0_sec, "_SECTION_PARAMETERS_", l_val=do_evgw0)
385 138 : CALL section_vals_val_get(evgw0_sec, "MAX_ITER", i_val=bs_env%evgw0_iter)
386 138 : CALL section_vals_val_get(evgw0_sec, "EPS_ITER", r_val=bs_env%evgw0_eps_iter)
387 :
388 138 : CALL timestop(handle)
389 :
390 138 : END SUBROUTINE read_evgw0_input
391 :
392 : ! **************************************************************************************************
393 : !> \brief Reads and validates the automatic RI input.
394 : !> \param auto_ri automatic RI basis optimization parameters
395 : !> \param auto_ri_sec AUTO_RI input section
396 : ! **************************************************************************************************
397 138 : SUBROUTINE read_auto_ri_input(auto_ri, auto_ri_sec)
398 : TYPE(auto_ri_type), INTENT(INOUT) :: auto_ri
399 : TYPE(section_vals_type), POINTER :: auto_ri_sec
400 :
401 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_auto_ri_input'
402 :
403 : INTEGER :: handle
404 :
405 138 : CALL timeset(routineN, handle)
406 :
407 138 : CALL section_vals_get(auto_ri_sec, explicit=auto_ri%enabled)
408 138 : IF (.NOT. auto_ri%enabled) THEN
409 134 : CALL timestop(handle)
410 134 : RETURN
411 : END IF
412 :
413 4 : CALL section_vals_val_get(auto_ri_sec, "RI_AO_RATIO", r_val=auto_ri%ri_ao_ratio)
414 4 : CALL section_vals_val_get(auto_ri_sec, "OCC_EMPTY_FRONTIER_ORBITAL_WINDOW", r_val=auto_ri%occ_energy_window)
415 4 : CALL section_vals_val_get(auto_ri_sec, "NEIGHBOR_RADIUS", r_val=auto_ri%neighbor_radius)
416 4 : IF (.NOT. (auto_ri%ri_ao_ratio > 0.0_dp)) CPABORT("AUTO_RI%RI_AO_RATIO must be positive")
417 4 : IF (.NOT. (auto_ri%occ_energy_window >= 0.0_dp)) THEN
418 0 : CPABORT("AUTO_RI%OCC_EMPTY_FRONTIER_ORBITAL_WINDOW must be nonnegative")
419 : END IF
420 4 : IF (.NOT. (auto_ri%neighbor_radius > 0.0_dp)) CPABORT("AUTO_RI%NEIGHBOR_RADIUS must be positive")
421 :
422 4 : CALL timestop(handle)
423 :
424 : END SUBROUTINE read_auto_ri_input
425 :
426 : ! **************************************************************************************************
427 : !> \brief Checks compatibility among the GW input sections and assigns the GW flavour.
428 : !> \param bs_env ...
429 : !> \param do_evgw0 whether eigenvalue self-consistency is requested
430 : ! **************************************************************************************************
431 138 : SUBROUTINE check_gw_input(bs_env, do_evgw0)
432 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
433 : LOGICAL, INTENT(IN) :: do_evgw0
434 :
435 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_gw_input'
436 :
437 : INTEGER :: handle
438 :
439 138 : CALL timeset(routineN, handle)
440 :
441 : ! GRID_OPTIMIZATION is available only with RI-RS.
442 138 : CPASSERT(.NOT. bs_env%ri_rs%grid_opt%enabled .OR. bs_env%do_gw_ri_rs)
443 138 : IF (bs_env%ri_rs%grid_opt%enabled) THEN
444 10 : IF (bs_env%ri_rs%grid_opt%max_iter < 1) THEN
445 0 : CPABORT("GRID_OPTIMIZATION%MAX_ITER must be positive")
446 : END IF
447 10 : IF (bs_env%ri_rs%grid_opt%cutoff_atomic_cluster <= 0.0_dp) THEN
448 0 : CPABORT("GRID_OPTIMIZATION%CUTOFF_ATOMIC_CLUSTER must be positive")
449 : END IF
450 10 : IF (bs_env%ri_rs%tikhonov < 0.0_dp) THEN
451 0 : CPABORT("RI_RS%TIKHONOV must not be negative")
452 : END IF
453 10 : IF (bs_env%do_periodic) THEN
454 0 : CPABORT("GRID_OPTIMIZATION currently supports nonperiodic local environments only")
455 : END IF
456 : END IF
457 :
458 : ! evGW0 implemented so far only with RI-RS.
459 138 : CPASSERT(.NOT. do_evgw0 .OR. bs_env%do_gw_ri_rs)
460 138 : IF (bs_env%auto_ri%enabled) THEN
461 : ! AUTO_RI is available only with RI-RS.
462 4 : CPASSERT(bs_env%do_gw_ri_rs)
463 : ! AUTO_RI is implemented only for nonperiodic systems.
464 4 : CPASSERT(.NOT. bs_env%do_periodic)
465 : END IF
466 :
467 138 : bs_env%gw_flavour = G0W0
468 138 : IF (do_evgw0) THEN
469 4 : bs_env%gw_flavour = evGW0
470 : CALL cp_warn(__LOCATION__, &
471 : "evGW0 in the RI-RS GW code is experimental. The quasiparticle energies "// &
472 : "of the previous cycle are fed back into the Green's function, so the "// &
473 : "error of the analytic continuation propagates and accumulates over the "// &
474 : "cycles instead of affecting one state only, as it does in G0W0. In "// &
475 : "tests on small molecules the resulting deviation stayed below 100 meV, "// &
476 : "but this has not been established in general. Use at your own risk and "// &
477 4 : "check the convergence of the reported states.")
478 : END IF
479 :
480 138 : CALL resolve_memory_per_proc(bs_env)
481 :
482 138 : CALL timestop(handle)
483 :
484 138 : END SUBROUTINE check_gw_input
485 :
486 : ! **************************************************************************************************
487 : !> \brief Determine memory per processor available for low-sclaing tensor GW routines.
488 : !>
489 : !> MEMORY_PER_PROC = -1 (the default) : Automatically detect available memory from the system.
490 : !> MEMORY_PER_PROC > 0.0 : Override with user-defined value (in GB).
491 : !>
492 : !> \param bs_env ...
493 : ! **************************************************************************************************
494 138 : SUBROUTINE resolve_memory_per_proc(bs_env)
495 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
496 :
497 : CHARACTER(LEN=*), PARAMETER :: routineN = 'resolve_memory_per_proc'
498 : REAL(KIND=dp), PARAMETER :: detected_memory_fraction = 0.8_dp, &
499 : fallback_memory_per_proc_GB = 2.0_dp
500 :
501 : INTEGER :: handle
502 : REAL(KIND=dp) :: mem_free_per_proc_GB, mem_occ_per_proc_GB
503 :
504 138 : CALL timeset(routineN, handle)
505 :
506 138 : bs_env%auto_memory_per_proc = bs_env%input_memory_per_proc_GB <= 0.0_dp
507 :
508 138 : IF (bs_env%auto_memory_per_proc) THEN
509 :
510 132 : CALL mp_mem_avail_per_rank_GB(bs_env%para_env, mem_free_per_proc_GB)
511 132 : CALL mp_mem_used_per_rank_GB(bs_env%para_env, mem_occ_per_proc_GB)
512 :
513 : ! detected_memory_fraction is a safety net on the detected free memory
514 : bs_env%input_memory_per_proc_GB = detected_memory_fraction*mem_free_per_proc_GB &
515 132 : + mem_occ_per_proc_GB
516 :
517 : ! mp_mem_avail_per_rank_GB returns zero if the free memory cannot be detected
518 132 : IF (mem_free_per_proc_GB <= 0.0_dp) THEN
519 0 : bs_env%auto_memory_per_proc = .FALSE.
520 0 : bs_env%input_memory_per_proc_GB = fallback_memory_per_proc_GB
521 : CALL cp_warn(__LOCATION__, &
522 : "Could not detect the available memory per MPI process. Falling back "// &
523 : "to MEMORY_PER_PROC = 2 GB. Set MEMORY_PER_PROC in the GW section "// &
524 0 : "explicitly for good performance.")
525 : END IF
526 :
527 6 : ELSE IF (bs_env%do_gw_ri_rs) THEN
528 :
529 : CALL cp_warn(__LOCATION__, &
530 : "MEMORY_PER_PROC is set, but it is not used by the RI-RS GW code, which "// &
531 : "detects the available memory per MPI process automatically. The keyword "// &
532 0 : "is ignored.")
533 :
534 : END IF
535 :
536 138 : CALL timestop(handle)
537 :
538 138 : END SUBROUTINE resolve_memory_per_proc
539 :
540 : ! **************************************************************************************************
541 : !> \brief Prints the common GW input parameters and the RI-RS block when active.
542 : !> \param bs_env ...
543 : ! **************************************************************************************************
544 138 : SUBROUTINE print_header_and_input_parameters(bs_env)
545 :
546 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
547 :
548 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_header_and_input_parameters'
549 :
550 : INTEGER :: handle, u
551 :
552 138 : CALL timeset(routineN, handle)
553 :
554 138 : u = bs_env%unit_nr
555 :
556 138 : IF (u > 0) THEN
557 69 : WRITE (u, '(T2,A)') ' '
558 69 : WRITE (u, '(T2,A)') REPEAT('-', 79)
559 69 : WRITE (u, '(T2,A,A78)') '-', '-'
560 69 : WRITE (u, '(T2,A,A46,A32)') '-', 'GW CALCULATION', '-'
561 69 : WRITE (u, '(T2,A,A78)') '-', '-'
562 69 : WRITE (u, '(T2,A)') REPEAT('-', 79)
563 69 : WRITE (u, '(T2,A)') ' '
564 69 : WRITE (u, '(T2,A,I45)') 'Input: Number of time/freq. points', bs_env%num_time_freq_points
565 69 : WRITE (u, "(T2,A,F44.1,A)") 'Input: ω_max for fitting Σ(iω) (eV)', bs_env%freq_max_fit*evolt
566 69 : WRITE (u, '(T2,A,ES27.1)') 'Input: Filter threshold for sparse tensor operations', &
567 138 : bs_env%eps_filter
568 69 : WRITE (u, "(T2,A,L55)") 'Input: Apply Hedin shift', bs_env%do_hedin_shift
569 69 : IF (bs_env%auto_memory_per_proc) THEN
570 66 : WRITE (u, '(T2,A,F34.1,A)') 'Detected: Available memory per MPI process', &
571 132 : bs_env%input_memory_per_proc_GB, ' GB'
572 : ELSE
573 3 : WRITE (u, '(T2,A,F37.1,A)') 'Input: Available memory per MPI process', &
574 6 : bs_env%input_memory_per_proc_GB, ' GB'
575 : END IF
576 69 : IF (bs_env%do_gw_ri_rs) CALL print_gw_rirs_input(bs_env)
577 69 : IF (bs_env%do_static_cohsex) THEN
578 4 : IF (bs_env%cohsex_flavor == cohsex_qp_update) THEN
579 3 : WRITE (u, FMT="(T2,A,T71,A10)") "static COHSEX approximation (COHSEX)", "QP_UPDATE"
580 : ELSE
581 1 : WRITE (u, FMT="(T2,A,T71,A10)") "static COHSEX approximation (COHSEX)", "DIAGONAL"
582 : END IF
583 4 : WRITE (u, FMT="(T2,A,T74,F8.6)") "static COHSEX scaling of Σ^c, occ-occ (COHSEX%SCALE_OCC):", &
584 8 : bs_env%cohsex_scale_occ
585 4 : WRITE (u, FMT="(T2,A,T74,F8.6)") "static COHSEX scaling of Σ^c, virt-virt (COHSEX%SCALE_VIRT):", &
586 8 : bs_env%cohsex_scale_virt
587 : END IF
588 69 : WRITE (u, '(A)') ' '
589 : END IF
590 :
591 138 : CALL timestop(handle)
592 :
593 138 : END SUBROUTINE print_header_and_input_parameters
594 :
595 : ! **************************************************************************************************
596 : !> \brief Prints the input parameters specific to an RI-RS GW calculation.
597 : !> \param bs_env ...
598 : ! **************************************************************************************************
599 20 : SUBROUTINE print_gw_rirs_input(bs_env)
600 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
601 :
602 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_gw_rirs_input'
603 :
604 : INTEGER :: handle, unit_nr
605 :
606 20 : CALL timeset(routineN, handle)
607 :
608 20 : unit_nr = bs_env%unit_nr
609 :
610 20 : WRITE (unit_nr, '(A)') ' '
611 20 : WRITE (unit_nr, '(T2,A,ES43.2)') 'Input: RI-RS Tikhonov regularization', &
612 40 : bs_env%ri_rs%tikhonov
613 : CALL print_cutoff_radius_input(unit_nr, 'AO grid hard cutoff radius', &
614 20 : bs_env%ri_rs%cutoff_radius_ri_ao)
615 : CALL print_cutoff_radius_input(unit_nr, 'G/W panel truncation radius', &
616 20 : bs_env%ri_rs%cutoff_radius_v_w)
617 20 : WRITE (unit_nr, '(T2,A,A62)') 'Input: GW flavour', TRIM(gw_flavour_label(bs_env))
618 :
619 20 : IF (bs_env%ri_rs%grid_opt%enabled) THEN
620 5 : WRITE (unit_nr, '(A)') ' '
621 : WRITE (unit_nr, '(T2,A,T80,L1)') &
622 5 : 'Input: RI-RS grid optimization activated', .TRUE.
623 : WRITE (unit_nr, '(T2,A,T69,F12.1)') &
624 5 : 'Input: RS points per AO function', bs_env%ri_rs%grid_opt%rs_ao_ratio
625 : WRITE (unit_nr, '(T2,A,T73,F6.1,A)') &
626 5 : 'Input: Cutoff radius for atomic clusters', &
627 10 : bs_env%ri_rs%grid_opt%cutoff_atomic_cluster*angstrom, ' Å'
628 : WRITE (unit_nr, '(T2,A,T75,I6)') &
629 5 : 'Input: Maximum grid optimization iterations', &
630 10 : bs_env%ri_rs%grid_opt%max_iter
631 : END IF
632 :
633 20 : IF (bs_env%auto_ri%enabled) THEN
634 2 : WRITE (unit_nr, '(A)') ' '
635 : WRITE (unit_nr, '(T2,A,T80,L1)') &
636 2 : 'Input: AUTO_RI basis optimization activated:', .TRUE.
637 : WRITE (unit_nr, '(T2,A,T69,F12.1)') &
638 2 : 'Input: AUTO_RI number of RI functions per AO function:', bs_env%auto_ri%ri_ao_ratio
639 : WRITE (unit_nr, '(T2,A,T74,F4.1,A)') &
640 2 : 'Input: AUTO_RI frontier orbital window:', bs_env%auto_ri%occ_energy_window*evolt, ' eV'
641 : WRITE (unit_nr, '(T2,A,T75,F4.1,A)') &
642 2 : 'Input: AUTO_RI neighbor radius:', bs_env%auto_ri%neighbor_radius*angstrom, ' Å'
643 : END IF
644 :
645 20 : IF (bs_env%gw_flavour == evGW0) THEN
646 2 : WRITE (unit_nr, '(T2,A,I42)') 'Input: Maximum number of evGW0 cycles', &
647 4 : bs_env%evgw0_iter
648 2 : WRITE (unit_nr, '(T2,A,F42.5,A)') 'Input: evGW0 convergence threshold', &
649 4 : bs_env%evgw0_eps_iter*evolt, ' eV'
650 : END IF
651 20 : WRITE (unit_nr, '(A)') ' '
652 :
653 20 : CALL timestop(handle)
654 :
655 20 : END SUBROUTINE print_gw_rirs_input
656 :
657 : ! **************************************************************************************************
658 : !> \brief Prints a positive RI-RS cutoff radius.
659 : !> \param unit_nr output unit
660 : !> \param label description of the cutoff radius
661 : !> \param cutoff_radius cutoff radius in atomic units
662 : ! **************************************************************************************************
663 40 : SUBROUTINE print_cutoff_radius_input(unit_nr, label, cutoff_radius)
664 : INTEGER, INTENT(IN) :: unit_nr
665 : CHARACTER(LEN=*), INTENT(IN) :: label
666 : REAL(KIND=dp), INTENT(IN) :: cutoff_radius
667 :
668 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_cutoff_radius_input'
669 :
670 : INTEGER :: handle
671 :
672 40 : CALL timeset(routineN, handle)
673 :
674 40 : IF (cutoff_radius > 0.0_dp) THEN
675 : WRITE (unit_nr, '(T2,2A,T73,F6.2,A)') &
676 2 : 'Input: ', TRIM(label), cutoff_radius*angstrom, ' Å'
677 : END IF
678 :
679 40 : CALL timestop(handle)
680 :
681 40 : END SUBROUTINE print_cutoff_radius_input
682 :
683 : ! **************************************************************************************************
684 : !> \brief Creates the AO and reference RI basis sets for every atomic kind.
685 : !> \param qs_env ...
686 : !> \param bs_env ...
687 : ! **************************************************************************************************
688 138 : SUBROUTINE setup_AO_and_RI_basis_set(qs_env, bs_env)
689 : TYPE(qs_environment_type), POINTER :: qs_env
690 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
691 :
692 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_AO_and_RI_basis_set'
693 :
694 : INTEGER :: handle, natom, nkind
695 138 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
696 138 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
697 :
698 138 : CALL timeset(routineN, handle)
699 :
700 : CALL get_qs_env(qs_env, &
701 : qs_kind_set=qs_kind_set, &
702 : particle_set=particle_set, &
703 138 : natom=natom, nkind=nkind)
704 :
705 : ! set up basis
706 552 : ALLOCATE (bs_env%sizes_RI(natom), bs_env%sizes_AO(natom))
707 952 : ALLOCATE (bs_env%basis_set_RI(nkind), bs_env%basis_set_AO(nkind))
708 :
709 138 : CALL basis_set_list_setup(bs_env%basis_set_RI, "RI_AUX", qs_kind_set)
710 138 : CALL basis_set_list_setup(bs_env%basis_set_AO, "ORB", qs_kind_set)
711 :
712 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_RI, &
713 138 : basis=bs_env%basis_set_RI)
714 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_AO, &
715 138 : basis=bs_env%basis_set_AO)
716 :
717 138 : CALL timestop(handle)
718 :
719 138 : END SUBROUTINE setup_AO_and_RI_basis_set
720 :
721 : ! **************************************************************************************************
722 : !> \brief Sets internal GW heuristics that are not exposed as user input.
723 : !> \param bs_env ...
724 : !> \param qs_env ...
725 : ! **************************************************************************************************
726 138 : SUBROUTINE set_heuristic_parameters(bs_env, qs_env)
727 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
728 : TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
729 :
730 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_heuristic_parameters'
731 :
732 : INTEGER :: handle, u
733 : LOGICAL :: do_BvK_cell
734 :
735 138 : CALL timeset(routineN, handle)
736 :
737 : ! for generating numerically stable minimax Fourier integration weights
738 138 : bs_env%num_points_per_magnitude = 200
739 :
740 138 : IF (bs_env%input_regularization_minimax > -1.0E-12_dp) THEN
741 0 : bs_env%regularization_minimax = bs_env%input_regularization_minimax
742 : ELSE
743 : ! for periodic systems and for 20 minimax points, we use a regularized minimax mesh
744 : ! (from experience: regularized minimax meshes converge faster for periodic systems
745 : ! and for 20 pts)
746 138 : IF (bs_env%do_periodic .OR. bs_env%num_time_freq_points >= 20) THEN
747 110 : bs_env%regularization_minimax = 1.0E-6_dp
748 : ELSE
749 28 : bs_env%regularization_minimax = 0.0_dp
750 : END IF
751 : END IF
752 :
753 138 : bs_env%stabilize_exp = 70.0_dp
754 138 : bs_env%eps_atom_grid_2d_mat = 1.0E-50_dp
755 :
756 : ! use a 16-parameter Padé fit
757 138 : bs_env%nparam_pade = 16
758 :
759 : ! resolution of the identity with the truncated Coulomb metric
760 138 : bs_env%ri_metric%potential_type = do_potential_truncated
761 138 : bs_env%ri_metric%omega = 0.0_dp
762 : ! cutoff radius is specified in the input
763 138 : bs_env%ri_metric%filename = "t_c_g.dat"
764 :
765 138 : bs_env%eps_eigval_mat_RI = 0.0_dp
766 :
767 138 : IF (bs_env%input_regularization_RI > -1.0E-12_dp) THEN
768 0 : bs_env%regularization_RI = bs_env%input_regularization_RI
769 : ELSE
770 : ! default case:
771 :
772 : ! 1. for periodic systems, we use the regularized resolution of the identity per default
773 138 : bs_env%regularization_RI = 1.0E-2_dp
774 :
775 : ! 2. for molecules, no regularization is necessary
776 138 : IF (.NOT. bs_env%do_periodic) bs_env%regularization_RI = 0.0_dp
777 :
778 : END IF
779 :
780 : ! Coulomb operator for the exchange self-energy. Periodic systems use the truncated
781 : ! operator described by Guidon, VandeVondele, and Hutter, JCTC 5, 3010 (2009).
782 138 : IF (bs_env%do_periodic) THEN
783 32 : do_BvK_cell = bs_env%gw_implementation == tensor_small_cell_full_kp
784 : CALL trunc_coulomb_for_exchange(qs_env, bs_env%trunc_coulomb, &
785 : rel_cutoff_trunc_coulomb_ri_x=0.5_dp, &
786 : cell_grid=bs_env%cell_grid_scf_desymm, &
787 32 : do_BvK_cell=do_BvK_cell)
788 : ELSE
789 106 : bs_env%trunc_coulomb%potential_type = do_potential_coulomb
790 106 : bs_env%trunc_coulomb%cutoff_radius = -1.0_dp
791 106 : bs_env%trunc_coulomb%omega = 0.0_dp
792 106 : bs_env%trunc_coulomb%filename = ""
793 : END IF
794 :
795 : ! for small-cell GW, we need more cells than normally used by the filter bs_env%eps_filter
796 : ! (in particular for computing the self-energy because of higher number of cells needed)
797 138 : bs_env%heuristic_filter_factor = 1.0E-4_dp
798 :
799 138 : u = bs_env%unit_nr
800 138 : IF (u > 0) THEN
801 69 : IF (bs_env%do_periodic) THEN
802 16 : WRITE (u, FMT="(T2,2A,F21.1,A)") "Cutoff radius for the truncated Coulomb ", &
803 32 : "operator in Σ^x:", bs_env%trunc_coulomb%cutoff_radius*angstrom, " Å"
804 : END IF
805 69 : WRITE (u, FMT="(T2,2A,F15.1,A)") "Cutoff radius for the truncated Coulomb ", &
806 138 : "operator in RI metric:", bs_env%ri_metric%cutoff_radius*angstrom, " Å"
807 69 : WRITE (u, FMT="(T2,A,ES48.1)") "Regularization parameter of RI ", bs_env%regularization_RI
808 69 : WRITE (u, FMT="(T2,A,ES38.1)") "Regularization parameter of minimax grids", &
809 138 : bs_env%regularization_minimax
810 69 : IF (bs_env%do_periodic) THEN
811 16 : WRITE (u, FMT="(T2,A,I53)") "Lattice sum size for V(k):", bs_env%size_lattice_sum_V
812 : END IF
813 : END IF
814 :
815 138 : CALL timestop(handle)
816 :
817 138 : END SUBROUTINE set_heuristic_parameters
818 :
819 : ! **************************************************************************************************
820 : !> \brief Sets up the RI basis, its matrix distributions, and M^-1(k=0).
821 : !> \param qs_env ...
822 : !> \param bs_env ...
823 : ! **************************************************************************************************
824 138 : SUBROUTINE get_RI_basis(qs_env, bs_env)
825 : TYPE(qs_environment_type), POINTER :: qs_env
826 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
827 :
828 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_RI_basis'
829 :
830 : INTEGER :: handle
831 :
832 138 : CALL timeset(routineN, handle)
833 :
834 138 : CALL set_AO_RI_basis_function_indices(qs_env, bs_env)
835 138 : CALL setup_kpoints_chi_eps_W(bs_env, bs_env%kpoints_chi_eps_W)
836 138 : IF (bs_env%gw_implementation == tensor_small_cell_full_kp) THEN
837 16 : CALL setup_cells_3c(qs_env, bs_env)
838 : END IF
839 138 : CALL set_parallelization_parameters(qs_env, bs_env)
840 138 : CALL allocate_matrices(qs_env, bs_env)
841 138 : IF (bs_env%gw_implementation /= tensor_small_cell_full_kp) THEN
842 122 : CALL compute_Minv_Gamma(qs_env, bs_env)
843 : END IF
844 :
845 138 : CALL timestop(handle)
846 :
847 138 : END SUBROUTINE get_RI_basis
848 :
849 : ! **************************************************************************************************
850 : !> \brief Determines the AO and RI basis-function indices for every atom.
851 : !> \param qs_env ...
852 : !> \param bs_env ...
853 : ! **************************************************************************************************
854 138 : SUBROUTINE set_AO_RI_basis_function_indices(qs_env, bs_env)
855 : TYPE(qs_environment_type), POINTER :: qs_env
856 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
857 :
858 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_AO_RI_basis_function_indices'
859 :
860 : INTEGER :: handle, i_RI, iatom, ikind, iset, &
861 : max_AO_bf_per_atom, n_ao_test, n_atom, &
862 : n_kind, n_RI, nset, nsgf, u
863 138 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
864 138 : INTEGER, DIMENSION(:), POINTER :: l_max, l_min, nsgf_set
865 138 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
866 : TYPE(gto_basis_set_type), POINTER :: basis
867 138 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
868 :
869 138 : CALL timeset(routineN, handle)
870 :
871 : ! determine RI basis set size
872 138 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
873 :
874 138 : n_kind = SIZE(qs_kind_set)
875 138 : n_atom = bs_env%n_atom
876 :
877 138 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
878 :
879 338 : DO ikind = 1, n_kind
880 200 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type="RI_AUX")
881 338 : CPASSERT(ASSOCIATED(basis))
882 : END DO
883 :
884 414 : ALLOCATE (bs_env%i_RI_start_from_atom(n_atom))
885 276 : ALLOCATE (bs_env%i_RI_end_from_atom(n_atom))
886 276 : ALLOCATE (bs_env%i_ao_start_from_atom(n_atom))
887 276 : ALLOCATE (bs_env%i_ao_end_from_atom(n_atom))
888 :
889 138 : n_RI = 0
890 466 : DO iatom = 1, n_atom
891 328 : bs_env%i_RI_start_from_atom(iatom) = n_RI + 1
892 328 : n_RI = n_RI + bs_env%sizes_RI(iatom)
893 466 : bs_env%i_RI_end_from_atom(iatom) = n_RI
894 : END DO
895 138 : bs_env%n_RI = n_RI
896 :
897 138 : max_AO_bf_per_atom = 0
898 138 : n_ao_test = 0
899 466 : DO iatom = 1, n_atom
900 328 : bs_env%i_ao_start_from_atom(iatom) = n_ao_test + 1
901 328 : ikind = kind_of(iatom)
902 328 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="ORB")
903 328 : n_ao_test = n_ao_test + nsgf
904 328 : bs_env%i_ao_end_from_atom(iatom) = n_ao_test
905 466 : max_AO_bf_per_atom = MAX(max_AO_bf_per_atom, nsgf)
906 : END DO
907 138 : CPASSERT(n_ao_test == bs_env%n_ao)
908 138 : bs_env%max_AO_bf_per_atom = max_AO_bf_per_atom
909 :
910 138 : u = bs_env%unit_nr
911 4658 : ALLOCATE (bs_env%l_RI(n_RI), source=-1)
912 138 : IF (bs_env%auto_ri%enabled) THEN
913 4 : CPASSERT(bs_env%auto_ri%ready)
914 : ELSE
915 : i_RI = 0
916 450 : DO iatom = 1, n_atom
917 316 : ikind = kind_of(iatom)
918 316 : nset = bs_env%basis_set_RI(ikind)%gto_basis_set%nset
919 316 : l_max => bs_env%basis_set_RI(ikind)%gto_basis_set%lmax
920 316 : l_min => bs_env%basis_set_RI(ikind)%gto_basis_set%lmin
921 316 : nsgf_set => bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf_set
922 2188 : DO iset = 1, nset
923 1738 : CPASSERT(l_max(iset) == l_min(iset))
924 5924 : bs_env%l_RI(i_RI + 1:i_RI + nsgf_set(iset)) = l_max(iset)
925 2054 : i_RI = i_RI + nsgf_set(iset)
926 : END DO
927 : END DO
928 134 : CPASSERT(i_RI == n_RI)
929 134 : IF (u > 0) THEN
930 67 : WRITE (u, FMT="(T2,A)") " "
931 67 : WRITE (u, FMT="(T2,2A,T75,I8)") "Number of auxiliary Gaussian basis functions ", &
932 134 : "for χ, ε, W", n_RI
933 : END IF
934 : END IF
935 :
936 138 : CALL timestop(handle)
937 :
938 276 : END SUBROUTINE set_AO_RI_basis_function_indices
939 :
940 : ! **************************************************************************************************
941 : !> \brief Sets the reciprocal mesh used for χ, ε, and W.
942 : !> \param bs_env ...
943 : !> \param kpoints_chi_eps_W reciprocal mesh used for χ, ε, and W
944 : ! **************************************************************************************************
945 138 : SUBROUTINE setup_kpoints_chi_eps_W(bs_env, kpoints_chi_eps_W)
946 :
947 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
948 : TYPE(kpoint_type), POINTER :: kpoints_chi_eps_W
949 :
950 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_kpoints_chi_eps_W'
951 :
952 : INTEGER :: handle, i_dim, n_dim, nkp, nkp_extra, &
953 : nkp_orig, u
954 : INTEGER, DIMENSION(3) :: nkp_grid, nkp_grid_extra, periodic
955 : REAL(KIND=dp) :: exp_s_p, n_dim_inv
956 :
957 138 : CALL timeset(routineN, handle)
958 :
959 : ! routine adapted from mp2_integrals.F
960 138 : NULLIFY (kpoints_chi_eps_W)
961 138 : CALL kpoint_create(kpoints_chi_eps_W)
962 :
963 138 : kpoints_chi_eps_W%kp_scheme = "GENERAL"
964 :
965 552 : periodic(1:3) = bs_env%periodic(1:3)
966 :
967 138 : CPASSERT(SIZE(bs_env%nkp_grid_chi_eps_W_input) == 3)
968 :
969 : IF (bs_env%nkp_grid_chi_eps_W_input(1) > 0 .AND. &
970 138 : bs_env%nkp_grid_chi_eps_W_input(2) > 0 .AND. &
971 : bs_env%nkp_grid_chi_eps_W_input(3) > 0) THEN
972 :
973 : ! 1. k-point mesh for χ, ε, W from input
974 :
975 16 : DO i_dim = 1, 3
976 4 : SELECT CASE (periodic(i_dim))
977 : CASE (0)
978 2 : nkp_grid(i_dim) = 1
979 2 : nkp_grid_extra(i_dim) = 1
980 : CASE (1)
981 10 : nkp_grid(i_dim) = bs_env%nkp_grid_chi_eps_W_input(i_dim)
982 10 : nkp_grid_extra(i_dim) = nkp_grid(i_dim)*2
983 : CASE DEFAULT
984 12 : CPABORT("Error in periodicity.")
985 : END SELECT
986 : END DO
987 :
988 : ELSE IF (bs_env%nkp_grid_chi_eps_W_input(1) == -1 .AND. &
989 134 : bs_env%nkp_grid_chi_eps_W_input(2) == -1 .AND. &
990 : bs_env%nkp_grid_chi_eps_W_input(3) == -1) THEN
991 :
992 : ! 2. automatic k-point mesh for χ, ε, W
993 :
994 536 : DO i_dim = 1, 3
995 :
996 402 : CPASSERT(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
997 :
998 134 : SELECT CASE (periodic(i_dim))
999 : CASE (0)
1000 346 : nkp_grid(i_dim) = 1
1001 346 : nkp_grid_extra(i_dim) = 1
1002 : CASE (1)
1003 80 : SELECT CASE (bs_env%gw_implementation)
1004 : CASE (tensor_large_cell_Gamma, ri_rs_large_cell_Gamma, ri_rs_non_periodic, tensor_static_cohsex)
1005 24 : nkp_grid(i_dim) = 4
1006 24 : nkp_grid_extra(i_dim) = 6
1007 : CASE (tensor_small_cell_full_kp)
1008 32 : nkp_grid(i_dim) = bs_env%kpoints_scf_desymm%nkp_grid(i_dim)*4
1009 56 : nkp_grid_extra(i_dim) = bs_env%kpoints_scf_desymm%nkp_grid(i_dim)*8
1010 : END SELECT
1011 : CASE DEFAULT
1012 402 : CPABORT("Error in periodicity.")
1013 : END SELECT
1014 :
1015 : END DO
1016 :
1017 : ELSE
1018 :
1019 0 : CPABORT("An error occured when setting up the k-mesh for W.")
1020 :
1021 : END IF
1022 :
1023 138 : nkp_orig = MAX(nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2, 1)
1024 :
1025 138 : nkp_extra = nkp_grid_extra(1)*nkp_grid_extra(2)*nkp_grid_extra(3)/2
1026 :
1027 138 : nkp = nkp_orig + nkp_extra
1028 :
1029 552 : kpoints_chi_eps_W%nkp_grid(1:3) = nkp_grid(1:3)
1030 138 : kpoints_chi_eps_W%nkp = nkp
1031 :
1032 552 : bs_env%nkp_grid_chi_eps_W_orig(1:3) = nkp_grid(1:3)
1033 552 : bs_env%nkp_grid_chi_eps_W_extra(1:3) = nkp_grid_extra(1:3)
1034 138 : bs_env%nkp_chi_eps_W_orig = nkp_orig
1035 138 : bs_env%nkp_chi_eps_W_orig_plus_extra = nkp
1036 :
1037 690 : ALLOCATE (kpoints_chi_eps_W%xkp(3, nkp), kpoints_chi_eps_W%wkp(nkp))
1038 414 : ALLOCATE (bs_env%wkp_no_extra(nkp), bs_env%wkp_s_p(nkp))
1039 :
1040 138 : CALL compute_xkp(kpoints_chi_eps_W%xkp, 1, nkp_orig, nkp_grid)
1041 138 : CALL compute_xkp(kpoints_chi_eps_W%xkp, nkp_orig + 1, nkp, nkp_grid_extra)
1042 :
1043 552 : n_dim = SUM(periodic)
1044 138 : IF (n_dim == 0) THEN
1045 : ! molecules
1046 106 : kpoints_chi_eps_W%wkp(1) = 1.0_dp
1047 106 : bs_env%wkp_s_p(1) = 1.0_dp
1048 106 : bs_env%wkp_no_extra(1) = 1.0_dp
1049 : ELSE
1050 :
1051 32 : n_dim_inv = 1.0_dp/REAL(n_dim, KIND=dp)
1052 :
1053 : ! k-point weights are chosen to automatically extrapolate the k-point mesh
1054 32 : CALL compute_wkp(kpoints_chi_eps_W%wkp(1:nkp_orig), nkp_orig, nkp_extra, n_dim_inv)
1055 : CALL compute_wkp(kpoints_chi_eps_W%wkp(nkp_orig + 1:nkp), &
1056 32 : nkp_extra, nkp_orig, n_dim_inv)
1057 :
1058 2200 : bs_env%wkp_no_extra(1:nkp_orig) = 0.0_dp
1059 8568 : bs_env%wkp_no_extra(nkp_orig + 1:nkp) = 1.0_dp/REAL(nkp_extra, KIND=dp)
1060 :
1061 32 : IF (n_dim == 3) THEN
1062 : ! W_PQ(k) for an s-function P and a p-function Q diverges as 1/k at k=0
1063 : ! (instead of 1/k^2 for P and Q both being s-functions).
1064 2 : exp_s_p = 2.0_dp*n_dim_inv
1065 2 : CALL compute_wkp(bs_env%wkp_s_p(1:nkp_orig), nkp_orig, nkp_extra, exp_s_p)
1066 2 : CALL compute_wkp(bs_env%wkp_s_p(nkp_orig + 1:nkp), nkp_extra, nkp_orig, exp_s_p)
1067 : ELSE
1068 10662 : bs_env%wkp_s_p(1:nkp) = bs_env%wkp_no_extra(1:nkp)
1069 : END IF
1070 :
1071 : END IF
1072 :
1073 138 : IF (bs_env%approx_kp_extrapol) THEN
1074 2 : bs_env%wkp_orig = 1.0_dp/REAL(nkp_orig, KIND=dp)
1075 : END IF
1076 :
1077 : ! heuristic parameter: how many k-points for χ, ε, and W are used simultaneously
1078 : ! (less simultaneous k-points: less memory, but more computational effort because of
1079 : ! recomputation of V(k))
1080 138 : bs_env%nkp_chi_eps_W_batch = 4
1081 :
1082 : bs_env%num_chi_eps_W_batches = (bs_env%nkp_chi_eps_W_orig_plus_extra - 1)/ &
1083 138 : bs_env%nkp_chi_eps_W_batch + 1
1084 :
1085 138 : u = bs_env%unit_nr
1086 :
1087 138 : IF (u > 0 .AND. bs_env%do_periodic) THEN
1088 16 : WRITE (u, FMT="(T2,A)") " "
1089 16 : WRITE (u, FMT="(T2,1A,T71,3I4)") "K-point mesh 1 for χ, ε, W", nkp_grid(1:3)
1090 16 : WRITE (u, FMT="(T2,2A,T71,3I4)") "K-point mesh 2 for χ, ε, W ", &
1091 32 : "(for k-point extrapolation of W)", nkp_grid_extra(1:3)
1092 16 : WRITE (u, FMT="(T2,A,T80,L)") "Approximate the k-point extrapolation", &
1093 32 : bs_env%approx_kp_extrapol
1094 : END IF
1095 :
1096 138 : CALL timestop(handle)
1097 :
1098 138 : END SUBROUTINE setup_kpoints_chi_eps_W
1099 :
1100 : ! **************************************************************************************************
1101 : !> \brief ...
1102 : !> \param xkp ...
1103 : !> \param ikp_start ...
1104 : !> \param ikp_end ...
1105 : !> \param grid ...
1106 : ! **************************************************************************************************
1107 276 : SUBROUTINE compute_xkp(xkp, ikp_start, ikp_end, grid)
1108 :
1109 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1110 : INTEGER :: ikp_start, ikp_end
1111 : INTEGER, DIMENSION(3) :: grid
1112 :
1113 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_xkp'
1114 :
1115 : INTEGER :: handle, i, ix, iy, iz
1116 :
1117 276 : CALL timeset(routineN, handle)
1118 :
1119 276 : i = ikp_start
1120 704 : DO ix = 1, grid(1)
1121 4636 : DO iy = 1, grid(2)
1122 25980 : DO iz = 1, grid(3)
1123 :
1124 21620 : IF (i > ikp_end) CYCLE
1125 :
1126 10810 : xkp(1, i) = REAL(2*ix - grid(1) - 1, KIND=dp)/(2._dp*REAL(grid(1), KIND=dp))
1127 10810 : xkp(2, i) = REAL(2*iy - grid(2) - 1, KIND=dp)/(2._dp*REAL(grid(2), KIND=dp))
1128 10810 : xkp(3, i) = REAL(2*iz - grid(3) - 1, KIND=dp)/(2._dp*REAL(grid(3), KIND=dp))
1129 25552 : i = i + 1
1130 :
1131 : END DO
1132 : END DO
1133 : END DO
1134 :
1135 276 : CALL timestop(handle)
1136 :
1137 276 : END SUBROUTINE compute_xkp
1138 :
1139 : ! **************************************************************************************************
1140 : !> \brief ...
1141 : !> \param wkp ...
1142 : !> \param nkp_1 ...
1143 : !> \param nkp_2 ...
1144 : !> \param exponent ...
1145 : ! **************************************************************************************************
1146 68 : SUBROUTINE compute_wkp(wkp, nkp_1, nkp_2, exponent)
1147 : REAL(KIND=dp), DIMENSION(:) :: wkp
1148 : INTEGER :: nkp_1, nkp_2
1149 : REAL(KIND=dp) :: exponent
1150 :
1151 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_wkp'
1152 :
1153 : INTEGER :: handle
1154 : REAL(KIND=dp) :: nkp_ratio
1155 :
1156 68 : CALL timeset(routineN, handle)
1157 :
1158 68 : nkp_ratio = REAL(nkp_2, KIND=dp)/REAL(nkp_1, KIND=dp)
1159 :
1160 10844 : wkp(:) = 1.0_dp/REAL(nkp_1, KIND=dp)/(1.0_dp - nkp_ratio**exponent)
1161 :
1162 68 : CALL timestop(handle)
1163 :
1164 68 : END SUBROUTINE compute_wkp
1165 :
1166 : ! **************************************************************************************************
1167 : !> \brief ...
1168 : !> \param qs_env ...
1169 : !> \param bs_env ...
1170 : ! **************************************************************************************************
1171 16 : SUBROUTINE setup_cells_3c(qs_env, bs_env)
1172 :
1173 : TYPE(qs_environment_type), POINTER :: qs_env
1174 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1175 :
1176 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_cells_3c'
1177 :
1178 : INTEGER :: atom_i, atom_j, atom_k, block_count, handle, i, i_cell_x, i_cell_x_max, &
1179 : i_cell_x_min, i_size, ikind, img, j, j_cell, j_cell_max, j_cell_y, j_cell_y_max, &
1180 : j_cell_y_min, j_size, k_cell, k_cell_max, k_cell_z, k_cell_z_max, k_cell_z_min, k_size, &
1181 : nimage_pairs_3c, nimages_3c, nimages_3c_max, nkind, u
1182 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, n_other_3c_images_max
1183 16 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c_max, nblocks_3c_max
1184 : INTEGER, DIMENSION(3) :: cell_index, n_max
1185 : REAL(KIND=dp) :: avail_mem_per_proc_GB, cell_dist, cell_radius_3c, dij, dik, djk, eps, &
1186 : exp_min_ao, exp_min_RI, frobenius_norm, mem_3c_GB, mem_occ_per_proc_GB, radius_ao, &
1187 : radius_ao_product, radius_RI
1188 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exp_ao_kind, exp_RI_kind, &
1189 16 : radius_ao_kind, &
1190 16 : radius_ao_product_kind, radius_RI_kind
1191 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c
1192 : REAL(KIND=dp), DIMENSION(3) :: rij, rik, rjk, vec_cell_j, vec_cell_k
1193 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: exp_ao, exp_RI
1194 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1195 : TYPE(cell_type), POINTER :: cell
1196 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1197 :
1198 16 : CALL timeset(routineN, handle)
1199 :
1200 16 : CALL get_qs_env(qs_env, nkind=nkind, atomic_kind_set=atomic_kind_set, particle_set=particle_set, cell=cell)
1201 :
1202 : ALLOCATE (exp_ao_kind(nkind), exp_RI_kind(nkind), radius_ao_kind(nkind), &
1203 112 : radius_ao_product_kind(nkind), radius_RI_kind(nkind))
1204 :
1205 48 : exp_min_RI = 10.0_dp
1206 48 : exp_min_ao = 10.0_dp
1207 48 : exp_RI_kind = 10.0_dp
1208 48 : exp_AO_kind = 10.0_dp
1209 :
1210 16 : eps = bs_env%eps_filter*bs_env%heuristic_filter_factor
1211 :
1212 48 : DO ikind = 1, nkind
1213 :
1214 32 : CALL get_gto_basis_set(bs_env%basis_set_RI(ikind)%gto_basis_set, zet=exp_RI)
1215 32 : CALL get_gto_basis_set(bs_env%basis_set_ao(ikind)%gto_basis_set, zet=exp_ao)
1216 :
1217 : ! we need to remove all exponents lower than a lower bound, e.g. 1E-3, because
1218 : ! for contracted basis sets, there might be exponents = 0 in zet
1219 64 : DO i = 1, SIZE(exp_RI, 1)
1220 112 : DO j = 1, SIZE(exp_RI, 2)
1221 48 : IF (exp_RI(i, j) < exp_min_RI .AND. exp_RI(i, j) > 1E-3_dp) exp_min_RI = exp_RI(i, j)
1222 80 : IF (exp_RI(i, j) < exp_RI_kind(ikind) .AND. exp_RI(i, j) > 1E-3_dp) THEN
1223 32 : exp_RI_kind(ikind) = exp_RI(i, j)
1224 : END IF
1225 : END DO
1226 : END DO
1227 160 : DO i = 1, SIZE(exp_ao, 1)
1228 384 : DO j = 1, SIZE(exp_ao, 2)
1229 224 : IF (exp_ao(i, j) < exp_min_ao .AND. exp_ao(i, j) > 1E-3_dp) exp_min_ao = exp_ao(i, j)
1230 352 : IF (exp_ao(i, j) < exp_ao_kind(ikind) .AND. exp_ao(i, j) > 1E-3_dp) THEN
1231 96 : exp_ao_kind(ikind) = exp_ao(i, j)
1232 : END IF
1233 : END DO
1234 : END DO
1235 32 : radius_ao_kind(ikind) = SQRT(-LOG(eps)/exp_ao_kind(ikind))
1236 32 : radius_ao_product_kind(ikind) = SQRT(-LOG(eps)/(2.0_dp*exp_ao_kind(ikind)))
1237 48 : radius_RI_kind(ikind) = SQRT(-LOG(eps)/exp_RI_kind(ikind))
1238 : END DO
1239 :
1240 16 : radius_ao = SQRT(-LOG(eps)/exp_min_ao)
1241 16 : radius_ao_product = SQRT(-LOG(eps)/(2.0_dp*exp_min_ao))
1242 16 : radius_RI = SQRT(-LOG(eps)/exp_min_RI)
1243 :
1244 16 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
1245 :
1246 : ! For a 3c integral (μR υS | P0) we have that cell R and cell S need to be within radius_3c
1247 16 : cell_radius_3c = radius_ao_product + radius_RI + bs_env%ri_metric%cutoff_radius
1248 :
1249 64 : n_max(1:3) = bs_env%periodic(1:3)*30
1250 :
1251 16 : nimages_3c_max = 0
1252 :
1253 16 : i_cell_x_min = 0
1254 16 : i_cell_x_max = 0
1255 16 : j_cell_y_min = 0
1256 16 : j_cell_y_max = 0
1257 16 : k_cell_z_min = 0
1258 16 : k_cell_z_max = 0
1259 :
1260 152 : DO i_cell_x = -n_max(1), n_max(1)
1261 8448 : DO j_cell_y = -n_max(2), n_max(2)
1262 67968 : DO k_cell_z = -n_max(3), n_max(3)
1263 :
1264 238144 : cell_index(1:3) = [i_cell_x, j_cell_y, k_cell_z]
1265 :
1266 59536 : CALL get_cell_dist(cell_index, bs_env%hmat, cell_dist)
1267 :
1268 67832 : IF (cell_dist < cell_radius_3c) THEN
1269 392 : nimages_3c_max = nimages_3c_max + 1
1270 392 : i_cell_x_min = MIN(i_cell_x_min, i_cell_x)
1271 392 : i_cell_x_max = MAX(i_cell_x_max, i_cell_x)
1272 392 : j_cell_y_min = MIN(j_cell_y_min, j_cell_y)
1273 392 : j_cell_y_max = MAX(j_cell_y_max, j_cell_y)
1274 392 : k_cell_z_min = MIN(k_cell_z_min, k_cell_z)
1275 392 : k_cell_z_max = MAX(k_cell_z_max, k_cell_z)
1276 : END IF
1277 :
1278 : END DO
1279 : END DO
1280 : END DO
1281 :
1282 : ! get index_to_cell_3c_max for the maximum possible cell range;
1283 : ! compute 3c integrals later in this routine and check really which cell is needed
1284 48 : ALLOCATE (index_to_cell_3c_max(3, nimages_3c_max))
1285 :
1286 16 : img = 0
1287 152 : DO i_cell_x = -n_max(1), n_max(1)
1288 8448 : DO j_cell_y = -n_max(2), n_max(2)
1289 67968 : DO k_cell_z = -n_max(3), n_max(3)
1290 :
1291 238144 : cell_index(1:3) = [i_cell_x, j_cell_y, k_cell_z]
1292 :
1293 59536 : CALL get_cell_dist(cell_index, bs_env%hmat, cell_dist)
1294 :
1295 67832 : IF (cell_dist < cell_radius_3c) THEN
1296 392 : img = img + 1
1297 1568 : index_to_cell_3c_max(1:3, img) = cell_index(1:3)
1298 : END IF
1299 :
1300 : END DO
1301 : END DO
1302 : END DO
1303 :
1304 : ! get pairs of R and S which have non-zero 3c integral (μR υS | P0)
1305 64 : ALLOCATE (nblocks_3c_max(nimages_3c_max, nimages_3c_max))
1306 16 : nblocks_3c_max(:, :) = 0
1307 :
1308 16 : block_count = 0
1309 408 : DO j_cell = 1, nimages_3c_max
1310 10040 : DO k_cell = 1, nimages_3c_max
1311 :
1312 38038 : DO atom_j = 1, bs_env%n_atom
1313 119924 : DO atom_k = 1, bs_env%n_atom
1314 353598 : DO atom_i = 1, bs_env%n_atom
1315 :
1316 243306 : block_count = block_count + 1
1317 243306 : IF (MODULO(block_count, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
1318 :
1319 486612 : CALL scaled_to_real(vec_cell_j, REAL(index_to_cell_3c_max(1:3, j_cell), kind=dp), cell)
1320 486612 : CALL scaled_to_real(vec_cell_k, REAL(index_to_cell_3c_max(1:3, k_cell), kind=dp), cell)
1321 :
1322 486612 : rij = pbc(particle_set(atom_j)%r(:), cell) - pbc(particle_set(atom_i)%r(:), cell) + vec_cell_j(:)
1323 : rjk = pbc(particle_set(atom_k)%r(:), cell) - pbc(particle_set(atom_j)%r(:), cell) &
1324 486612 : + vec_cell_k(:) - vec_cell_j(:)
1325 486612 : rik(:) = rij(:) + rjk(:)
1326 486612 : dij = NORM2(rij)
1327 486612 : dik = NORM2(rik)
1328 486612 : djk = NORM2(rjk)
1329 121653 : IF (djk > radius_ao_kind(kind_of(atom_j)) + radius_ao_kind(kind_of(atom_k))) CYCLE
1330 38127 : IF (dij > radius_ao_kind(kind_of(atom_j)) + radius_RI_kind(kind_of(atom_i)) &
1331 : + bs_env%ri_metric%cutoff_radius) CYCLE
1332 19497 : IF (dik > radius_RI_kind(kind_of(atom_i)) + radius_ao_kind(kind_of(atom_k)) &
1333 : + bs_env%ri_metric%cutoff_radius) CYCLE
1334 :
1335 12822 : j_size = bs_env%i_ao_end_from_atom(atom_j) - bs_env%i_ao_start_from_atom(atom_j) + 1
1336 12822 : k_size = bs_env%i_ao_end_from_atom(atom_k) - bs_env%i_ao_start_from_atom(atom_k) + 1
1337 12822 : i_size = bs_env%i_RI_end_from_atom(atom_i) - bs_env%i_RI_start_from_atom(atom_i) + 1
1338 :
1339 64110 : ALLOCATE (int_3c(j_size, k_size, i_size))
1340 :
1341 : ! compute 3-c int. ( μ(atom j) R , ν (atom k) S | P (atom i) 0 )
1342 : ! ("|": truncated Coulomb operator), inside build_3c_integrals: (j k | i)
1343 : CALL build_3c_integral_block(int_3c, qs_env, bs_env%ri_metric, &
1344 : basis_j=bs_env%basis_set_AO, &
1345 : basis_k=bs_env%basis_set_AO, &
1346 : basis_i=bs_env%basis_set_RI, &
1347 : cell_j=index_to_cell_3c_max(1:3, j_cell), &
1348 : cell_k=index_to_cell_3c_max(1:3, k_cell), &
1349 12822 : atom_k=atom_k, atom_j=atom_j, atom_i=atom_i)
1350 :
1351 678851 : frobenius_norm = NORM2(int_3c)
1352 :
1353 12822 : DEALLOCATE (int_3c)
1354 :
1355 : ! we use a higher threshold here to safe memory when storing the 3c integrals
1356 : ! in every tensor group
1357 95100 : IF (frobenius_norm > eps) THEN
1358 2720 : nblocks_3c_max(j_cell, k_cell) = nblocks_3c_max(j_cell, k_cell) + 1
1359 : END IF
1360 :
1361 : END DO
1362 : END DO
1363 : END DO
1364 :
1365 : END DO
1366 : END DO
1367 :
1368 16 : CALL bs_env%para_env%sum(nblocks_3c_max)
1369 :
1370 48 : ALLOCATE (n_other_3c_images_max(nimages_3c_max))
1371 16 : n_other_3c_images_max(:) = 0
1372 :
1373 16 : nimages_3c = 0
1374 16 : nimage_pairs_3c = 0
1375 :
1376 408 : DO j_cell = 1, nimages_3c_max
1377 10024 : DO k_cell = 1, nimages_3c_max
1378 10024 : IF (nblocks_3c_max(j_cell, k_cell) > 0) THEN
1379 960 : n_other_3c_images_max(j_cell) = n_other_3c_images_max(j_cell) + 1
1380 960 : nimage_pairs_3c = nimage_pairs_3c + 1
1381 : END IF
1382 : END DO
1383 :
1384 408 : IF (n_other_3c_images_max(j_cell) > 0) nimages_3c = nimages_3c + 1
1385 :
1386 : END DO
1387 :
1388 16 : bs_env%nimages_3c = nimages_3c
1389 48 : ALLOCATE (bs_env%index_to_cell_3c(3, nimages_3c))
1390 : ALLOCATE (bs_env%cell_to_index_3c(i_cell_x_min:i_cell_x_max, &
1391 : j_cell_y_min:j_cell_y_max, &
1392 80 : k_cell_z_min:k_cell_z_max))
1393 848 : bs_env%cell_to_index_3c(:, :, :) = -1
1394 :
1395 64 : ALLOCATE (bs_env%nblocks_3c(nimages_3c, nimages_3c))
1396 16 : bs_env%nblocks_3c(nimages_3c, nimages_3c) = 0
1397 :
1398 16 : j_cell = 0
1399 408 : DO j_cell_max = 1, nimages_3c_max
1400 392 : IF (n_other_3c_images_max(j_cell_max) == 0) CYCLE
1401 178 : j_cell = j_cell + 1
1402 712 : cell_index(1:3) = index_to_cell_3c_max(1:3, j_cell_max)
1403 712 : bs_env%index_to_cell_3c(1:3, j_cell) = cell_index(1:3)
1404 178 : bs_env%cell_to_index_3c(cell_index(1), cell_index(2), cell_index(3)) = j_cell
1405 :
1406 178 : k_cell = 0
1407 4604 : DO k_cell_max = 1, nimages_3c_max
1408 4410 : IF (n_other_3c_images_max(k_cell_max) == 0) CYCLE
1409 2066 : k_cell = k_cell + 1
1410 :
1411 4802 : bs_env%nblocks_3c(j_cell, k_cell) = nblocks_3c_max(j_cell_max, k_cell_max)
1412 : END DO
1413 :
1414 : END DO
1415 :
1416 : ! we use: 8*10^-9 GB / double precision number
1417 : mem_3c_GB = REAL(bs_env%n_RI, KIND=dp)*REAL(bs_env%n_ao, KIND=dp)**2 &
1418 16 : *REAL(nimage_pairs_3c, KIND=dp)*8E-9_dp
1419 :
1420 16 : CALL mp_mem_used_per_rank_GB(bs_env%para_env, mem_occ_per_proc_GB)
1421 :
1422 : ! number of processors per group that entirely stores the 3c integrals and does tensor ops
1423 16 : avail_mem_per_proc_GB = bs_env%input_memory_per_proc_GB - mem_occ_per_proc_GB
1424 :
1425 : ! careful: downconvering real to integer, 1.9 -> 1; thus add 1.0 for upconversion, 1.9 -> 2
1426 16 : bs_env%group_size_tensor = MAX(INT(mem_3c_GB/avail_mem_per_proc_GB + 1.0_dp), 1)
1427 :
1428 16 : u = bs_env%unit_nr
1429 :
1430 16 : IF (u > 0) THEN
1431 8 : WRITE (u, FMT="(T2,A,F52.1,A)") "Radius of atomic orbitals", radius_ao*angstrom, " Å"
1432 8 : WRITE (u, FMT="(T2,A,F55.1,A)") "Radius of RI functions", radius_RI*angstrom, " Å"
1433 8 : WRITE (u, FMT="(T2,A,I47)") "Number of cells for 3c integrals", nimages_3c
1434 8 : WRITE (u, FMT="(T2,A,I42)") "Number of cell pairs for 3c integrals", nimage_pairs_3c
1435 8 : WRITE (u, '(T2,A)') ''
1436 8 : IF (bs_env%auto_memory_per_proc) THEN
1437 8 : WRITE (u, '(T2,A,F34.1,A)') 'Detected: Available memory per MPI process', &
1438 16 : bs_env%input_memory_per_proc_GB, ' GB'
1439 : ELSE
1440 0 : WRITE (u, '(T2,A,F37.1,A)') 'Input: Available memory per MPI process', &
1441 0 : bs_env%input_memory_per_proc_GB, ' GB'
1442 : END IF
1443 8 : WRITE (u, '(T2,A,F35.1,A)') 'Used memory per MPI process before GW run', &
1444 16 : mem_occ_per_proc_GB, ' GB'
1445 8 : WRITE (u, '(T2,A,F44.1,A)') 'Memory of three-center integrals', mem_3c_GB, ' GB'
1446 : END IF
1447 :
1448 16 : CALL timestop(handle)
1449 :
1450 32 : END SUBROUTINE setup_cells_3c
1451 :
1452 : ! **************************************************************************************************
1453 : !> \brief ...
1454 : !> \param cell_index ...
1455 : !> \param hmat ...
1456 : !> \param cell_dist ...
1457 : ! **************************************************************************************************
1458 119072 : SUBROUTINE get_cell_dist(cell_index, hmat, cell_dist)
1459 :
1460 : INTEGER, DIMENSION(3) :: cell_index
1461 : REAL(KIND=dp) :: hmat(3, 3), cell_dist
1462 :
1463 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_cell_dist'
1464 :
1465 : INTEGER :: handle, i_dim
1466 : INTEGER, DIMENSION(3) :: cell_index_adj
1467 : REAL(KIND=dp) :: cell_dist_3(3)
1468 :
1469 119072 : CALL timeset(routineN, handle)
1470 :
1471 : ! the distance of cells needs to be taken to adjacent neighbors, not
1472 : ! between the center of the cells. We thus need to rescale the cell index
1473 476288 : DO i_dim = 1, 3
1474 357216 : IF (cell_index(i_dim) > 0) cell_index_adj(i_dim) = cell_index(i_dim) - 1
1475 357216 : IF (cell_index(i_dim) < 0) cell_index_adj(i_dim) = cell_index(i_dim) + 1
1476 476288 : IF (cell_index(i_dim) == 0) cell_index_adj(i_dim) = cell_index(i_dim)
1477 : END DO
1478 :
1479 1905152 : cell_dist_3(1:3) = MATMUL(hmat, REAL(cell_index_adj, KIND=dp))
1480 :
1481 476288 : cell_dist = NORM2(cell_dist_3)
1482 :
1483 119072 : CALL timestop(handle)
1484 :
1485 119072 : END SUBROUTINE get_cell_dist
1486 :
1487 : ! **************************************************************************************************
1488 : !> \brief ...
1489 : !> \param qs_env ...
1490 : !> \param bs_env ...
1491 : ! **************************************************************************************************
1492 138 : SUBROUTINE set_parallelization_parameters(qs_env, bs_env)
1493 : TYPE(qs_environment_type), POINTER :: qs_env
1494 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1495 :
1496 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_parallelization_parameters'
1497 :
1498 : INTEGER :: color_sub, dummy_1, dummy_2, handle, &
1499 : num_pe, num_t_groups, u
1500 : TYPE(mp_para_env_type), POINTER :: para_env
1501 :
1502 138 : CALL timeset(routineN, handle)
1503 :
1504 138 : CALL get_qs_env(qs_env, para_env=para_env)
1505 :
1506 138 : num_pe = para_env%num_pe
1507 : ! if not already set, use all processors for the group (for large-cell GW, performance
1508 : ! seems to be best for a single group with all MPI processes per group)
1509 138 : IF (bs_env%group_size_tensor < 0 .OR. bs_env%group_size_tensor > num_pe) THEN
1510 122 : bs_env%group_size_tensor = num_pe
1511 : END IF
1512 :
1513 : ! group_size_tensor must divide num_pe without rest; otherwise everything will be complicated
1514 138 : IF (MODULO(num_pe, bs_env%group_size_tensor) /= 0) THEN
1515 0 : CALL find_good_group_size(num_pe, bs_env%group_size_tensor)
1516 : END IF
1517 :
1518 : ! para_env_tensor for tensor subgroups
1519 138 : color_sub = para_env%mepos/bs_env%group_size_tensor
1520 138 : bs_env%tensor_group_color = color_sub
1521 :
1522 138 : ALLOCATE (bs_env%para_env_tensor)
1523 138 : CALL bs_env%para_env_tensor%from_split(para_env, color_sub)
1524 :
1525 138 : num_t_groups = para_env%num_pe/bs_env%group_size_tensor
1526 138 : bs_env%num_tensor_groups = num_t_groups
1527 :
1528 : CALL get_i_j_atoms(bs_env%atoms_i, bs_env%atoms_j, bs_env%n_atom_i, bs_env%n_atom_j, &
1529 138 : color_sub, bs_env)
1530 :
1531 414 : ALLOCATE (bs_env%atoms_i_t_group(2, num_t_groups))
1532 276 : ALLOCATE (bs_env%atoms_j_t_group(2, num_t_groups))
1533 292 : DO color_sub = 0, num_t_groups - 1
1534 : CALL get_i_j_atoms(bs_env%atoms_i_t_group(1:2, color_sub + 1), &
1535 : bs_env%atoms_j_t_group(1:2, color_sub + 1), &
1536 292 : dummy_1, dummy_2, color_sub, bs_env)
1537 : END DO
1538 :
1539 138 : u = bs_env%unit_nr
1540 138 : IF (u > 0 .AND. .NOT. bs_env%do_gw_ri_rs) THEN
1541 49 : WRITE (u, '(T2,A,I47)') 'Group size for tensor operations', bs_env%group_size_tensor
1542 49 : IF (bs_env%group_size_tensor > 1 .AND. bs_env%n_atom < 5) THEN
1543 41 : WRITE (u, '(T2,A)') 'The requested group size is > 1 which can lead to bad performance.'
1544 41 : WRITE (u, '(T2,A)') 'Using more memory per MPI process might improve performance.'
1545 41 : WRITE (u, '(T2,A)') '(Also increase MEMORY_PER_PROC when using more memory per process.)'
1546 : END IF
1547 : END IF
1548 :
1549 138 : CALL timestop(handle)
1550 :
1551 138 : END SUBROUTINE set_parallelization_parameters
1552 :
1553 : ! **************************************************************************************************
1554 : !> \brief ...
1555 : !> \param num_pe ...
1556 : !> \param group_size ...
1557 : ! **************************************************************************************************
1558 0 : SUBROUTINE find_good_group_size(num_pe, group_size)
1559 :
1560 : INTEGER :: num_pe, group_size
1561 :
1562 : CHARACTER(LEN=*), PARAMETER :: routineN = 'find_good_group_size'
1563 :
1564 : INTEGER :: group_size_minus, group_size_orig, &
1565 : group_size_plus, handle, i_diff
1566 :
1567 0 : CALL timeset(routineN, handle)
1568 :
1569 0 : group_size_orig = group_size
1570 :
1571 0 : DO i_diff = 1, num_pe
1572 :
1573 0 : group_size_minus = group_size - i_diff
1574 :
1575 0 : IF (MODULO(num_pe, group_size_minus) == 0 .AND. group_size_minus > 0) THEN
1576 0 : group_size = group_size_minus
1577 0 : EXIT
1578 : END IF
1579 :
1580 0 : group_size_plus = group_size + i_diff
1581 :
1582 0 : IF (MODULO(num_pe, group_size_plus) == 0 .AND. group_size_plus <= num_pe) THEN
1583 0 : group_size = group_size_plus
1584 0 : EXIT
1585 : END IF
1586 :
1587 : END DO
1588 :
1589 0 : IF (group_size_orig == group_size) CPABORT("Group size error")
1590 :
1591 0 : CALL timestop(handle)
1592 :
1593 0 : END SUBROUTINE find_good_group_size
1594 :
1595 : ! **************************************************************************************************
1596 : !> \brief ...
1597 : !> \param atoms_i ...
1598 : !> \param atoms_j ...
1599 : !> \param n_atom_i ...
1600 : !> \param n_atom_j ...
1601 : !> \param color_sub ...
1602 : !> \param bs_env ...
1603 : ! **************************************************************************************************
1604 292 : SUBROUTINE get_i_j_atoms(atoms_i, atoms_j, n_atom_i, n_atom_j, color_sub, bs_env)
1605 :
1606 : INTEGER, DIMENSION(2) :: atoms_i, atoms_j
1607 : INTEGER :: n_atom_i, n_atom_j, color_sub
1608 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1609 :
1610 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_i_j_atoms'
1611 :
1612 : INTEGER :: handle, i_atoms_per_group, i_group, &
1613 : ipcol, ipcol_loop, iprow, iprow_loop, &
1614 : j_atoms_per_group, npcol, nprow
1615 :
1616 292 : CALL timeset(routineN, handle)
1617 :
1618 : ! create a square mesh of tensor groups for iatom and jatom; code from blacs_env_create
1619 292 : CALL square_mesh(nprow, npcol, bs_env%num_tensor_groups)
1620 :
1621 292 : i_group = 0
1622 584 : DO ipcol_loop = 0, npcol - 1
1623 924 : DO iprow_loop = 0, nprow - 1
1624 340 : IF (i_group == color_sub) THEN
1625 292 : iprow = iprow_loop
1626 292 : ipcol = ipcol_loop
1627 : END IF
1628 632 : i_group = i_group + 1
1629 : END DO
1630 : END DO
1631 :
1632 292 : IF (MODULO(bs_env%n_atom, nprow) == 0) THEN
1633 250 : i_atoms_per_group = bs_env%n_atom/nprow
1634 : ELSE
1635 42 : i_atoms_per_group = bs_env%n_atom/nprow + 1
1636 : END IF
1637 :
1638 292 : IF (MODULO(bs_env%n_atom, npcol) == 0) THEN
1639 292 : j_atoms_per_group = bs_env%n_atom/npcol
1640 : ELSE
1641 0 : j_atoms_per_group = bs_env%n_atom/npcol + 1
1642 : END IF
1643 :
1644 292 : atoms_i(1) = iprow*i_atoms_per_group + 1
1645 292 : atoms_i(2) = MIN((iprow + 1)*i_atoms_per_group, bs_env%n_atom)
1646 292 : n_atom_i = atoms_i(2) - atoms_i(1) + 1
1647 :
1648 292 : atoms_j(1) = ipcol*j_atoms_per_group + 1
1649 292 : atoms_j(2) = MIN((ipcol + 1)*j_atoms_per_group, bs_env%n_atom)
1650 292 : n_atom_j = atoms_j(2) - atoms_j(1) + 1
1651 :
1652 292 : CALL timestop(handle)
1653 :
1654 292 : END SUBROUTINE get_i_j_atoms
1655 :
1656 : ! **************************************************************************************************
1657 : !> \brief ...
1658 : !> \param nprow ...
1659 : !> \param npcol ...
1660 : !> \param nproc ...
1661 : ! **************************************************************************************************
1662 292 : SUBROUTINE square_mesh(nprow, npcol, nproc)
1663 : INTEGER :: nprow, npcol, nproc
1664 :
1665 : CHARACTER(LEN=*), PARAMETER :: routineN = 'square_mesh'
1666 :
1667 : INTEGER :: gcd_max, handle, ipe, jpe
1668 :
1669 292 : CALL timeset(routineN, handle)
1670 :
1671 292 : gcd_max = -1
1672 632 : DO ipe = 1, CEILING(SQRT(REAL(nproc, dp)))
1673 340 : jpe = nproc/ipe
1674 340 : IF (ipe*jpe /= nproc) CYCLE
1675 632 : IF (gcd(ipe, jpe) >= gcd_max) THEN
1676 340 : nprow = ipe
1677 340 : npcol = jpe
1678 340 : gcd_max = gcd(ipe, jpe)
1679 : END IF
1680 : END DO
1681 :
1682 292 : CALL timestop(handle)
1683 :
1684 292 : END SUBROUTINE square_mesh
1685 :
1686 : ! **************************************************************************************************
1687 : !> \brief ...
1688 : !> \param qs_env ...
1689 : !> \param bs_env ...
1690 : ! **************************************************************************************************
1691 138 : SUBROUTINE allocate_matrices(qs_env, bs_env)
1692 : TYPE(qs_environment_type), POINTER :: qs_env
1693 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1694 :
1695 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_matrices'
1696 :
1697 : INTEGER :: handle, i_t
1698 : TYPE(cp_blacs_env_type), POINTER :: blacs_env, blacs_env_tensor
1699 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_RI_global
1700 : TYPE(mp_para_env_type), POINTER :: para_env
1701 :
1702 138 : CALL timeset(routineN, handle)
1703 :
1704 138 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1705 :
1706 138 : fm_struct => bs_env%fm_ks_Gamma(1)%matrix_struct
1707 :
1708 138 : CALL cp_fm_create(bs_env%fm_Gocc, fm_struct)
1709 138 : CALL cp_fm_create(bs_env%fm_Gvir, fm_struct)
1710 :
1711 138 : NULLIFY (fm_struct_RI_global)
1712 : CALL cp_fm_struct_create(fm_struct_RI_global, context=blacs_env, nrow_global=bs_env%n_RI, &
1713 138 : ncol_global=bs_env%n_RI, para_env=para_env)
1714 138 : CALL cp_fm_create(bs_env%fm_RI_RI, fm_struct_RI_global)
1715 138 : CALL cp_fm_create(bs_env%fm_chi_Gamma_freq, fm_struct_RI_global)
1716 138 : CALL cp_fm_create(bs_env%fm_W_MIC_freq, fm_struct_RI_global)
1717 138 : IF (bs_env%approx_kp_extrapol) THEN
1718 2 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_extra, fm_struct_RI_global)
1719 2 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_no_extra, fm_struct_RI_global)
1720 2 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_extra, 0.0_dp)
1721 2 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_no_extra, 0.0_dp)
1722 : END IF
1723 138 : CALL cp_fm_struct_release(fm_struct_RI_global)
1724 :
1725 138 : IF (.NOT. bs_env%do_gw_ri_rs) THEN
1726 : ! create blacs_env for subgroups of tensor operations
1727 98 : NULLIFY (blacs_env_tensor)
1728 98 : CALL cp_blacs_env_create(blacs_env=blacs_env_tensor, para_env=bs_env%para_env_tensor)
1729 :
1730 : ! allocate dbcsr matrices in the tensor subgroup; actually, one only needs a small
1731 : ! subset of blocks in the tensor subgroup, however, all atomic blocks are allocated.
1732 : ! One might think of creating a dbcsr matrix with only the blocks that are needed
1733 : ! in the tensor subgroup
1734 : CALL create_mat_munu(bs_env%mat_ao_ao_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
1735 98 : blacs_env_tensor, do_ri_aux_basis=.FALSE.)
1736 :
1737 : CALL create_mat_munu(bs_env%mat_RI_RI_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
1738 98 : blacs_env_tensor, do_ri_aux_basis=.TRUE.)
1739 :
1740 98 : CALL cp_blacs_env_release(blacs_env_tensor)
1741 : END IF
1742 :
1743 : CALL create_mat_munu(bs_env%mat_RI_RI, qs_env, bs_env%eps_atom_grid_2d_mat, &
1744 : blacs_env, do_ri_aux_basis=.TRUE., &
1745 138 : custom_row_blk_sizes=bs_env%sizes_RI)
1746 :
1747 138 : NULLIFY (bs_env%mat_chi_Gamma_tau)
1748 138 : CALL dbcsr_allocate_matrix_set(bs_env%mat_chi_Gamma_tau, bs_env%num_time_freq_points)
1749 :
1750 2208 : DO i_t = 1, bs_env%num_time_freq_points
1751 2070 : ALLOCATE (bs_env%mat_chi_Gamma_tau(i_t)%matrix)
1752 2208 : CALL dbcsr_create(bs_env%mat_chi_Gamma_tau(i_t)%matrix, template=bs_env%mat_RI_RI%matrix)
1753 : END DO
1754 :
1755 138 : CALL timestop(handle)
1756 :
1757 138 : END SUBROUTINE allocate_matrices
1758 :
1759 : ! **************************************************************************************************
1760 : !> \brief Computes and stores the inverse RI metric at the Γ point.
1761 : !> \param qs_env ...
1762 : !> \param bs_env ...
1763 : ! **************************************************************************************************
1764 122 : SUBROUTINE compute_Minv_Gamma(qs_env, bs_env)
1765 : TYPE(qs_environment_type), POINTER :: qs_env
1766 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1767 :
1768 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Minv_Gamma'
1769 :
1770 : INTEGER :: handle
1771 : REAL(KIND=dp) :: eigenvalue_threshold, &
1772 : metric_regularization
1773 122 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_M
1774 :
1775 122 : CALL timeset(routineN, handle)
1776 :
1777 122 : CALL cp_fm_create(bs_env%fm_Minv_Gamma, bs_env%fm_RI_RI%matrix_struct)
1778 122 : IF (bs_env%auto_ri%enabled) THEN
1779 4 : CPASSERT(bs_env%auto_ri%ready)
1780 4 : CALL cp_fm_to_fm(bs_env%auto_ri%M_pq_inv, bs_env%fm_Minv_Gamma)
1781 : ELSE
1782 118 : metric_regularization = 0.0_dp
1783 118 : eigenvalue_threshold = 0.0_dp
1784 118 : IF (bs_env%do_gw_ri_rs) THEN
1785 36 : metric_regularization = bs_env%regularization_RI
1786 36 : eigenvalue_threshold = bs_env%eps_eigval_mat_RI
1787 : END IF
1788 : CALL RI_2c_integral_mat(qs_env, fm_M, bs_env%fm_RI_RI%matrix_struct, bs_env%n_RI, &
1789 : bs_env%ri_metric, &
1790 118 : regularization_RI=metric_regularization)
1791 118 : CALL fm_invert(fm_M(1, 1), eigenvalue_threshold, bs_env%unit_nr)
1792 118 : CALL cp_fm_to_fm(fm_M(1, 1), bs_env%fm_Minv_Gamma)
1793 118 : CALL cp_fm_release(fm_M)
1794 : END IF
1795 :
1796 122 : CALL timestop(handle)
1797 :
1798 122 : END SUBROUTINE compute_Minv_Gamma
1799 :
1800 : ! **************************************************************************************************
1801 : !> \brief ...
1802 : !> \param qs_env ...
1803 : !> \param bs_env ...
1804 : ! **************************************************************************************************
1805 276 : SUBROUTINE compute_V_xc(qs_env, bs_env)
1806 : TYPE(qs_environment_type), POINTER :: qs_env
1807 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1808 :
1809 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_xc'
1810 :
1811 : INTEGER :: handle, img, ispin, myfun, nimages
1812 : LOGICAL :: hf_present
1813 : REAL(KIND=dp) :: energy_ex, energy_exc, energy_total, &
1814 : myfraction
1815 138 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_ks_without_v_xc
1816 138 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
1817 : TYPE(dft_control_type), POINTER :: dft_control
1818 : TYPE(qs_energy_type), POINTER :: energy
1819 : TYPE(section_vals_type), POINTER :: hf_section, input, xc_section
1820 :
1821 138 : CALL timeset(routineN, handle)
1822 :
1823 138 : CALL get_qs_env(qs_env, input=input, energy=energy, dft_control=dft_control)
1824 :
1825 : ! previously, dft_control%nimages set to # neighbor cells, revert for Γ-only KS matrix
1826 138 : nimages = dft_control%nimages
1827 138 : dft_control%nimages = bs_env%nimages_scf
1828 :
1829 : ! we need to reset XC functional, therefore, get XC input
1830 138 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1831 138 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
1832 138 : CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=xc_none)
1833 138 : hf_section => section_vals_get_subs_vals(input, "DFT%XC%HF", can_return_null=.TRUE.)
1834 138 : hf_present = .FALSE.
1835 138 : IF (ASSOCIATED(hf_section)) THEN
1836 138 : CALL section_vals_get(hf_section, explicit=hf_present)
1837 : END IF
1838 138 : IF (hf_present) THEN
1839 : ! Special case for handling hfx
1840 64 : CALL section_vals_val_get(xc_section, "HF%FRACTION", r_val=myfraction)
1841 64 : CALL section_vals_val_set(xc_section, "HF%FRACTION", r_val=0.0_dp)
1842 : END IF
1843 :
1844 : ! save the energy before the energy gets updated
1845 138 : energy_total = energy%total
1846 138 : energy_exc = energy%exc
1847 138 : energy_ex = energy%ex
1848 :
1849 260 : SELECT CASE (bs_env%gw_implementation)
1850 : CASE (tensor_large_cell_Gamma, ri_rs_large_cell_Gamma, ri_rs_non_periodic, tensor_static_cohsex)
1851 :
1852 122 : NULLIFY (mat_ks_without_v_xc)
1853 122 : CALL dbcsr_allocate_matrix_set(mat_ks_without_v_xc, bs_env%n_spin)
1854 :
1855 262 : DO ispin = 1, bs_env%n_spin
1856 140 : ALLOCATE (mat_ks_without_v_xc(ispin)%matrix)
1857 262 : IF (hf_present) THEN
1858 : CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix, &
1859 76 : matrix_type=dbcsr_type_symmetric)
1860 : ELSE
1861 64 : CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1862 : END IF
1863 : END DO
1864 :
1865 : ! calculate KS-matrix without XC
1866 : CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
1867 122 : ext_ks_matrix=mat_ks_without_v_xc)
1868 :
1869 262 : DO ispin = 1, bs_env%n_spin
1870 : ! transfer dbcsr matrix to fm
1871 140 : CALL cp_fm_create(bs_env%fm_V_xc_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1872 140 : CALL copy_dbcsr_to_fm(mat_ks_without_v_xc(ispin)%matrix, bs_env%fm_V_xc_Gamma(ispin))
1873 :
1874 : ! v_xc = h_ks - h_ks(v_xc = 0)
1875 : CALL cp_fm_scale_and_add(alpha=-1.0_dp, matrix_a=bs_env%fm_V_xc_Gamma(ispin), &
1876 262 : beta=1.0_dp, matrix_b=bs_env%fm_ks_Gamma(ispin))
1877 : END DO
1878 :
1879 122 : CALL dbcsr_deallocate_matrix_set(mat_ks_without_v_xc)
1880 :
1881 : CASE (tensor_small_cell_full_kp)
1882 :
1883 : ! calculate KS-matrix without XC
1884 16 : CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
1885 16 : CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_kp)
1886 :
1887 608 : ALLOCATE (bs_env%fm_V_xc_R(dft_control%nimages, bs_env%n_spin))
1888 170 : DO ispin = 1, bs_env%n_spin
1889 560 : DO img = 1, dft_control%nimages
1890 : ! safe fm_V_xc_R in fm_matrix because saving in dbcsr matrix caused trouble...
1891 528 : CALL copy_dbcsr_to_fm(matrix_ks_kp(ispin, img)%matrix, bs_env%fm_work_mo(1))
1892 : CALL cp_fm_create(bs_env%fm_V_xc_R(img, ispin), bs_env%fm_work_mo(1)%matrix_struct, &
1893 528 : set_zero=.TRUE.)
1894 : ! store h_ks(v_xc = 0) in fm_V_xc_R
1895 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=bs_env%fm_V_xc_R(img, ispin), &
1896 544 : beta=1.0_dp, matrix_b=bs_env%fm_work_mo(1))
1897 : END DO
1898 : END DO
1899 :
1900 : END SELECT
1901 :
1902 : ! set back the energy
1903 138 : energy%total = energy_total
1904 138 : energy%exc = energy_exc
1905 138 : energy%ex = energy_ex
1906 :
1907 : ! set back nimages
1908 138 : dft_control%nimages = nimages
1909 :
1910 : ! set the DFT functional and HF fraction back
1911 : CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
1912 138 : i_val=myfun)
1913 138 : IF (hf_present) THEN
1914 : CALL section_vals_val_set(xc_section, "HF%FRACTION", &
1915 64 : r_val=myfraction)
1916 : END IF
1917 :
1918 138 : IF (bs_env%gw_implementation == tensor_small_cell_full_kp) THEN
1919 : ! calculate KS-matrix again with XC
1920 16 : CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
1921 32 : DO ispin = 1, bs_env%n_spin
1922 560 : DO img = 1, dft_control%nimages
1923 : ! store h_ks in fm_work_mo
1924 528 : CALL copy_dbcsr_to_fm(matrix_ks_kp(ispin, img)%matrix, bs_env%fm_work_mo(1))
1925 : ! v_xc = h_ks - h_ks(v_xc = 0)
1926 : CALL cp_fm_scale_and_add(alpha=-1.0_dp, matrix_a=bs_env%fm_V_xc_R(img, ispin), &
1927 544 : beta=1.0_dp, matrix_b=bs_env%fm_work_mo(1))
1928 : END DO
1929 : END DO
1930 : END IF
1931 :
1932 138 : CALL timestop(handle)
1933 :
1934 138 : END SUBROUTINE compute_V_xc
1935 :
1936 : ! **************************************************************************************************
1937 : !> \brief Sets the interaction radii used by the GW calculation.
1938 : !> \param bs_env ...
1939 : ! **************************************************************************************************
1940 138 : SUBROUTINE init_interaction_radii(bs_env)
1941 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1942 :
1943 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_interaction_radii'
1944 :
1945 : INTEGER :: handle, ibasis
1946 : TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
1947 :
1948 138 : CALL timeset(routineN, handle)
1949 :
1950 338 : DO ibasis = 1, SIZE(bs_env%basis_set_AO)
1951 :
1952 200 : orb_basis => bs_env%basis_set_AO(ibasis)%gto_basis_set
1953 200 : CALL init_interaction_radii_orb_basis(orb_basis, bs_env%eps_filter)
1954 :
1955 200 : ri_basis => bs_env%basis_set_RI(ibasis)%gto_basis_set
1956 338 : CALL init_interaction_radii_orb_basis(ri_basis, bs_env%eps_filter)
1957 :
1958 : END DO
1959 :
1960 138 : CALL timestop(handle)
1961 :
1962 138 : END SUBROUTINE init_interaction_radii
1963 :
1964 : ! **************************************************************************************************
1965 : !> \brief ...
1966 : !> \param qs_env ...
1967 : !> \param bs_env ...
1968 : ! **************************************************************************************************
1969 112 : SUBROUTINE create_tensors(qs_env, bs_env)
1970 : TYPE(qs_environment_type), POINTER :: qs_env
1971 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1972 :
1973 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_tensors'
1974 :
1975 : INTEGER :: handle
1976 :
1977 112 : CALL timeset(routineN, handle)
1978 :
1979 : ! split blocks does not improve load balancing/efficienfy for tensor contraction, so we go
1980 : ! with the standard atomic blocks
1981 : CALL create_3c_t(bs_env%t_RI_AO__AO, bs_env%para_env_tensor, "(RI AO | AO)", [1, 2], [3], &
1982 : bs_env%sizes_RI, bs_env%sizes_AO, &
1983 112 : create_nl_3c=.TRUE., nl_3c=bs_env%nl_3c, qs_env=qs_env)
1984 : CALL create_3c_t(bs_env%t_RI__AO_AO, bs_env%para_env_tensor, "(RI | AO AO)", [1], [2, 3], &
1985 112 : bs_env%sizes_RI, bs_env%sizes_AO)
1986 :
1987 112 : CALL create_2c_t(bs_env)
1988 :
1989 112 : CALL timestop(handle)
1990 :
1991 112 : END SUBROUTINE create_tensors
1992 :
1993 : ! **************************************************************************************************
1994 : !> \brief ...
1995 : !> \param tensor ...
1996 : !> \param para_env ...
1997 : !> \param tensor_name ...
1998 : !> \param map1 ...
1999 : !> \param map2 ...
2000 : !> \param sizes_RI ...
2001 : !> \param sizes_AO ...
2002 : !> \param create_nl_3c ...
2003 : !> \param nl_3c ...
2004 : !> \param qs_env ...
2005 : ! **************************************************************************************************
2006 224 : SUBROUTINE create_3c_t(tensor, para_env, tensor_name, map1, map2, sizes_RI, sizes_AO, &
2007 : create_nl_3c, nl_3c, qs_env)
2008 : TYPE(dbt_type) :: tensor
2009 : TYPE(mp_para_env_type), POINTER :: para_env
2010 : CHARACTER(LEN=12) :: tensor_name
2011 : INTEGER, DIMENSION(:) :: map1, map2
2012 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_RI, sizes_AO
2013 : LOGICAL, OPTIONAL :: create_nl_3c
2014 : TYPE(neighbor_list_3c_type), OPTIONAL :: nl_3c
2015 : TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
2016 :
2017 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_3c_t'
2018 :
2019 : INTEGER :: handle, nkind
2020 224 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_AO_1, dist_AO_2, dist_RI
2021 : INTEGER, DIMENSION(3) :: pcoord, pdims, pdims_3d
2022 : LOGICAL :: my_create_nl_3c
2023 672 : TYPE(dbt_pgrid_type) :: pgrid_3d
2024 : TYPE(distribution_3d_type) :: dist_3d
2025 224 : TYPE(mp_cart_type) :: mp_comm_t3c_2
2026 224 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2027 :
2028 224 : CALL timeset(routineN, handle)
2029 :
2030 224 : pdims_3d = 0
2031 224 : CALL dbt_pgrid_create(para_env, pdims_3d, pgrid_3d)
2032 : CALL create_3c_tensor(tensor, dist_RI, dist_AO_1, dist_AO_2, &
2033 : pgrid_3d, sizes_RI, sizes_AO, sizes_AO, &
2034 224 : map1=map1, map2=map2, name=tensor_name)
2035 :
2036 224 : IF (PRESENT(create_nl_3c)) THEN
2037 112 : my_create_nl_3c = create_nl_3c
2038 : ELSE
2039 : my_create_nl_3c = .FALSE.
2040 : END IF
2041 :
2042 112 : IF (my_create_nl_3c) THEN
2043 112 : CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set)
2044 112 : CALL dbt_mp_environ_pgrid(pgrid_3d, pdims, pcoord)
2045 112 : CALL mp_comm_t3c_2%create(pgrid_3d%mp_comm_2d, 3, pdims)
2046 : CALL distribution_3d_create(dist_3d, dist_RI, dist_AO_1, dist_AO_2, &
2047 112 : nkind, particle_set, mp_comm_t3c_2, own_comm=.TRUE.)
2048 :
2049 : CALL build_3c_neighbor_lists(nl_3c, &
2050 : qs_env%bs_env%basis_set_RI, &
2051 : qs_env%bs_env%basis_set_AO, &
2052 : qs_env%bs_env%basis_set_AO, &
2053 : dist_3d, qs_env%bs_env%ri_metric, &
2054 112 : "GW_3c_nl", qs_env, own_dist=.TRUE.)
2055 : END IF
2056 :
2057 224 : DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
2058 224 : CALL dbt_pgrid_destroy(pgrid_3d)
2059 :
2060 224 : CALL timestop(handle)
2061 :
2062 448 : END SUBROUTINE create_3c_t
2063 :
2064 : ! **************************************************************************************************
2065 : !> \brief Creates the AO-AO and RI-RI two-center tensors.
2066 : !> \param bs_env ...
2067 : ! **************************************************************************************************
2068 112 : SUBROUTINE create_2c_t(bs_env)
2069 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2070 :
2071 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_2c_t'
2072 :
2073 : INTEGER :: handle
2074 112 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2
2075 : INTEGER, DIMENSION(2) :: pdims_2d
2076 336 : TYPE(dbt_pgrid_type) :: pgrid_2d
2077 :
2078 112 : CALL timeset(routineN, handle)
2079 :
2080 : ! inspired from rpa_im_time.F / hfx_types.F
2081 :
2082 112 : pdims_2d = 0
2083 112 : CALL dbt_pgrid_create(bs_env%para_env_tensor, pdims_2d, pgrid_2d)
2084 :
2085 : CALL create_2c_tensor(bs_env%t_G, dist_1, dist_2, pgrid_2d, &
2086 : bs_env%sizes_AO, bs_env%sizes_AO, &
2087 112 : name="(AO | AO)")
2088 112 : DEALLOCATE (dist_1, dist_2)
2089 : CALL create_2c_tensor(bs_env%t_chi, dist_1, dist_2, pgrid_2d, &
2090 : bs_env%sizes_RI, bs_env%sizes_RI, &
2091 112 : name="(RI | RI)")
2092 112 : DEALLOCATE (dist_1, dist_2)
2093 : CALL create_2c_tensor(bs_env%t_W, dist_1, dist_2, pgrid_2d, &
2094 : bs_env%sizes_RI, bs_env%sizes_RI, &
2095 112 : name="(RI | RI)")
2096 112 : DEALLOCATE (dist_1, dist_2)
2097 112 : CALL dbt_pgrid_destroy(pgrid_2d)
2098 :
2099 112 : CALL timestop(handle)
2100 :
2101 112 : END SUBROUTINE create_2c_t
2102 :
2103 : ! **************************************************************************************************
2104 : !> \brief ...
2105 : !> \param qs_env ...
2106 : !> \param bs_env ...
2107 : ! **************************************************************************************************
2108 82 : SUBROUTINE check_sparsity_3c(qs_env, bs_env)
2109 : TYPE(qs_environment_type), POINTER :: qs_env
2110 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2111 :
2112 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_sparsity_3c'
2113 :
2114 : INTEGER :: handle, n_atom_step, RI_atom
2115 : INTEGER(int_8) :: non_zero_elements_sum, nze
2116 : REAL(dp) :: max_dist_AO_atoms, occ, occupation_sum
2117 : REAL(KIND=dp) :: t1, t2
2118 82 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_global_array
2119 :
2120 82 : CALL timeset(routineN, handle)
2121 :
2122 : ! check the sparsity of 3c integral tensor (µν|P); calculate maximum distance between
2123 : ! AO atoms µ, ν where at least a single integral (µν|P) is larger than the filter threshold
2124 :
2125 738 : ALLOCATE (t_3c_global_array(1, 1))
2126 82 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_global_array(1, 1))
2127 :
2128 : ! Allocate arrays to store min/max indices for overlap with other AO/RI functions on each atom
2129 : ! (Filled during loop via get_i_j_atom_ranges)
2130 328 : ALLOCATE (bs_env%min_RI_idx_from_AO_AO_atom(bs_env%n_atom, bs_env%n_atom))
2131 328 : ALLOCATE (bs_env%max_RI_idx_from_AO_AO_atom(bs_env%n_atom, bs_env%n_atom))
2132 328 : ALLOCATE (bs_env%min_AO_idx_from_RI_AO_atom(bs_env%n_atom, bs_env%n_atom))
2133 328 : ALLOCATE (bs_env%max_AO_idx_from_RI_AO_atom(bs_env%n_atom, bs_env%n_atom))
2134 670 : bs_env%min_RI_idx_from_AO_AO_atom(:, :) = bs_env%n_RI
2135 670 : bs_env%max_RI_idx_from_AO_AO_atom(:, :) = 1
2136 670 : bs_env%min_AO_idx_from_RI_AO_atom(:, :) = bs_env%n_AO
2137 670 : bs_env%max_AO_idx_from_RI_AO_atom(:, :) = 1
2138 :
2139 82 : CALL bs_env%para_env%sync()
2140 82 : t1 = m_walltime()
2141 :
2142 82 : occupation_sum = 0.0_dp
2143 82 : non_zero_elements_sum = 0
2144 82 : max_dist_AO_atoms = 0.0_dp
2145 82 : n_atom_step = INT(SQRT(REAL(bs_env%n_atom, KIND=dp)))
2146 : ! do not compute full 3c integrals at once because it may cause out of memory
2147 262 : DO RI_atom = 1, bs_env%n_atom, n_atom_step
2148 :
2149 : CALL build_3c_integrals(t_3c_global_array, &
2150 : bs_env%eps_filter, &
2151 : qs_env, &
2152 : bs_env%nl_3c, &
2153 : int_eps=bs_env%eps_filter, &
2154 : basis_i=bs_env%basis_set_RI, &
2155 : basis_j=bs_env%basis_set_AO, &
2156 : basis_k=bs_env%basis_set_AO, &
2157 : bounds_i=[RI_atom, MIN(RI_atom + n_atom_step - 1, bs_env%n_atom)], &
2158 : potential_parameter=bs_env%ri_metric, &
2159 540 : desymmetrize=.FALSE.)
2160 :
2161 180 : CALL dbt_filter(t_3c_global_array(1, 1), bs_env%eps_filter)
2162 :
2163 180 : CALL bs_env%para_env%sync()
2164 :
2165 180 : CALL get_tensor_occupancy(t_3c_global_array(1, 1), nze, occ)
2166 180 : non_zero_elements_sum = non_zero_elements_sum + nze
2167 180 : occupation_sum = occupation_sum + occ
2168 :
2169 180 : CALL get_max_dist_AO_atoms(t_3c_global_array(1, 1), max_dist_AO_atoms, qs_env)
2170 :
2171 : ! Extract indices per block
2172 180 : CALL get_i_j_atom_ranges(t_3c_global_array(1, 1), bs_env)
2173 :
2174 442 : CALL dbt_clear(t_3c_global_array(1, 1))
2175 :
2176 : END DO
2177 :
2178 82 : t2 = m_walltime()
2179 :
2180 : ! Sync/max for max_dist_AO_atoms is done inside each get_max_dist_AO_atoms
2181 82 : bs_env%max_dist_AO_atoms = max_dist_AO_atoms
2182 : ! occupation_sum is a global quantity, also needs no sync here
2183 82 : bs_env%occupation_3c_int = occupation_sum
2184 :
2185 82 : CALL bs_env%para_env%min(bs_env%min_RI_idx_from_AO_AO_atom)
2186 82 : CALL bs_env%para_env%max(bs_env%max_RI_idx_from_AO_AO_atom)
2187 82 : CALL bs_env%para_env%min(bs_env%min_AO_idx_from_RI_AO_atom)
2188 82 : CALL bs_env%para_env%max(bs_env%max_AO_idx_from_RI_AO_atom)
2189 :
2190 82 : CALL dbt_destroy(t_3c_global_array(1, 1))
2191 164 : DEALLOCATE (t_3c_global_array)
2192 :
2193 82 : IF (bs_env%unit_nr > 0) THEN
2194 41 : WRITE (bs_env%unit_nr, '(T2,A)') ''
2195 : WRITE (bs_env%unit_nr, '(T2,A,F27.1,A)') &
2196 41 : 'Computed 3-center integrals (µν|P), execution time', t2 - t1, ' s'
2197 41 : WRITE (bs_env%unit_nr, '(T2,A,F48.3,A)') 'Percentage of non-zero (µν|P)', &
2198 82 : bs_env%occupation_3c_int*100, ' %'
2199 41 : WRITE (bs_env%unit_nr, '(T2,A,F33.1,A)') 'Max. distance between µ,ν in non-zero (µν|P)', &
2200 82 : bs_env%max_dist_AO_atoms*angstrom, ' A'
2201 41 : WRITE (bs_env%unit_nr, '(T2,2A,I20,A)') 'Required memory if storing all 3-center ', &
2202 82 : 'integrals (µν|P)', INT(REAL(non_zero_elements_sum, KIND=dp)*8.0E-9_dp), ' GB'
2203 : END IF
2204 :
2205 82 : CALL timestop(handle)
2206 :
2207 164 : END SUBROUTINE check_sparsity_3c
2208 :
2209 : ! **************************************************************************************************
2210 : !> \brief ...
2211 : !> \param t_3c_int ...
2212 : !> \param max_dist_AO_atoms ...
2213 : !> \param qs_env ...
2214 : ! **************************************************************************************************
2215 180 : SUBROUTINE get_max_dist_AO_atoms(t_3c_int, max_dist_AO_atoms, qs_env)
2216 : TYPE(dbt_type) :: t_3c_int
2217 : REAL(KIND=dp), INTENT(INOUT) :: max_dist_AO_atoms
2218 : TYPE(qs_environment_type), POINTER :: qs_env
2219 :
2220 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_max_dist_AO_atoms'
2221 :
2222 : INTEGER :: atom_1, atom_2, handle, num_cells
2223 : INTEGER, DIMENSION(3) :: atom_ind
2224 180 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
2225 : REAL(KIND=dp) :: abs_rab
2226 : REAL(KIND=dp), DIMENSION(3) :: rab
2227 : TYPE(cell_type), POINTER :: cell
2228 : TYPE(dbt_iterator_type) :: iter
2229 : TYPE(mp_para_env_type), POINTER :: para_env
2230 180 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2231 :
2232 180 : CALL timeset(routineN, handle)
2233 :
2234 180 : NULLIFY (cell, particle_set, para_env)
2235 180 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, para_env=para_env)
2236 :
2237 : ! max_dist_AO_atoms is compared to earlier steps in the loop with step n_atom_step
2238 : ! do not initialize/overwrite here
2239 :
2240 : ! IMPORTANT: Use thread-local copy for max_dist_AO_atoms via REDUCTION to avoid race conditions
2241 : !$OMP PARALLEL DEFAULT(NONE) &
2242 : !$OMP SHARED(t_3c_int, num_cells, index_to_cell, particle_set, cell) &
2243 : !$OMP PRIVATE(iter, atom_ind, rab, abs_rab, atom_1, atom_2) &
2244 180 : !$OMP REDUCTION(MAX:max_dist_AO_atoms)
2245 :
2246 : CALL dbt_iterator_start(iter, t_3c_int)
2247 : DO WHILE (dbt_iterator_blocks_left(iter))
2248 : CALL dbt_iterator_next_block(iter, atom_ind)
2249 :
2250 : atom_1 = atom_ind(2)
2251 : atom_2 = atom_ind(3)
2252 : rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
2253 : abs_rab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
2254 :
2255 : ! Reduction takes care of using a thread-local copy
2256 : max_dist_AO_atoms = MAX(max_dist_AO_atoms, abs_rab)
2257 :
2258 : END DO
2259 : CALL dbt_iterator_stop(iter)
2260 : !$OMP END PARALLEL
2261 :
2262 180 : CALL para_env%max(max_dist_AO_atoms)
2263 :
2264 180 : CALL timestop(handle)
2265 :
2266 180 : END SUBROUTINE get_max_dist_AO_atoms
2267 :
2268 : ! **************************************************************************************************
2269 : !> \brief ...
2270 : !> \param t_3c ...
2271 : !> \param bs_env ...
2272 : ! **************************************************************************************************
2273 180 : SUBROUTINE get_i_j_atom_ranges(t_3c, bs_env)
2274 : TYPE(dbt_type) :: t_3c
2275 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2276 :
2277 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_i_j_atom_ranges'
2278 :
2279 : INTEGER :: handle, idx_AO_end, idx_AO_start, &
2280 : idx_RI_end, idx_RI_start
2281 : INTEGER, DIMENSION(3) :: atom_ind
2282 : TYPE(dbt_iterator_type) :: iter
2283 :
2284 180 : CALL timeset(routineN, handle)
2285 :
2286 : ! Loop over blocks in 3c, for given min_atom: RI_min/max index from min_atom
2287 : !$OMP PARALLEL DEFAULT(NONE) &
2288 : !$OMP SHARED(t_3c, bs_env) &
2289 : !$OMP PRIVATE(iter, atom_ind, &
2290 180 : !$OMP idx_RI_start, idx_RI_end, idx_AO_start, idx_AO_end)
2291 :
2292 : CALL dbt_iterator_start(iter, t_3c)
2293 : DO WHILE (dbt_iterator_blocks_left(iter))
2294 : CALL dbt_iterator_next_block(iter, atom_ind)
2295 :
2296 : ! Pre-fetch indices to avoid referencing 'bs_env' twice inside the ATOMIC blocks
2297 : idx_RI_start = bs_env%i_RI_start_from_atom(atom_ind(1))
2298 : idx_RI_end = bs_env%i_RI_end_from_atom(atom_ind(1))
2299 :
2300 : idx_AO_start = bs_env%i_ao_start_from_atom(atom_ind(2))
2301 : idx_AO_end = bs_env%i_ao_end_from_atom(atom_ind(2))
2302 :
2303 : ! Update values safely inside ATOMIC blocks, otherwise race conditions occur
2304 : !$OMP ATOMIC UPDATE
2305 : bs_env%min_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)) = &
2306 : MIN(bs_env%min_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)), idx_RI_start)
2307 : !$OMP ATOMIC UPDATE
2308 : bs_env%max_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)) = &
2309 : MAX(bs_env%max_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)), idx_RI_end)
2310 :
2311 : !$OMP ATOMIC UPDATE
2312 : bs_env%min_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)) = &
2313 : MIN(bs_env%min_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)), idx_AO_start)
2314 : !$OMP ATOMIC UPDATE
2315 : bs_env%max_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)) = &
2316 : MAX(bs_env%max_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)), idx_AO_end)
2317 :
2318 : END DO
2319 : CALL dbt_iterator_stop(iter)
2320 : !$OMP END PARALLEL
2321 :
2322 180 : CALL timestop(handle)
2323 :
2324 180 : END SUBROUTINE get_i_j_atom_ranges
2325 :
2326 : ! **************************************************************************************************
2327 : !> \brief ...
2328 : !> \param bs_env ...
2329 : ! **************************************************************************************************
2330 82 : SUBROUTINE set_sparsity_parallelization_parameters(bs_env)
2331 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2332 :
2333 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_sparsity_parallelization_parameters'
2334 :
2335 : INTEGER :: handle, i_ivl, IL_ivl, j_ivl, n_atom_per_IL_ivl, n_atom_per_ivl, n_intervals_i, &
2336 : n_intervals_inner_loop_atoms, n_intervals_j, u
2337 : INTEGER(KIND=int_8) :: input_memory_per_proc
2338 :
2339 82 : CALL timeset(routineN, handle)
2340 :
2341 : ! heuristic parameter to prevent out of memory
2342 82 : bs_env%safety_factor_memory = 0.10_dp
2343 :
2344 82 : input_memory_per_proc = INT(bs_env%input_memory_per_proc_GB*1.0E9_dp, KIND=int_8)
2345 :
2346 : ! choose atomic range for λ ("i_atom"), ν ("j_atom") in
2347 : ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
2348 : ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
2349 : ! such that M and N fit into the memory
2350 : n_atom_per_ivl = INT(SQRT(bs_env%safety_factor_memory*input_memory_per_proc &
2351 : *bs_env%group_size_tensor/24/bs_env%n_RI &
2352 82 : /SQRT(bs_env%occupation_3c_int)))/bs_env%max_AO_bf_per_atom
2353 :
2354 82 : n_intervals_i = (bs_env%n_atom_i - 1)/n_atom_per_ivl + 1
2355 82 : n_intervals_j = (bs_env%n_atom_j - 1)/n_atom_per_ivl + 1
2356 :
2357 82 : bs_env%n_intervals_i = n_intervals_i
2358 82 : bs_env%n_intervals_j = n_intervals_j
2359 :
2360 246 : ALLOCATE (bs_env%i_atom_intervals(2, n_intervals_i))
2361 246 : ALLOCATE (bs_env%j_atom_intervals(2, n_intervals_j))
2362 :
2363 164 : DO i_ivl = 1, n_intervals_i
2364 82 : bs_env%i_atom_intervals(1, i_ivl) = (i_ivl - 1)*n_atom_per_ivl + bs_env%atoms_i(1)
2365 : bs_env%i_atom_intervals(2, i_ivl) = MIN(i_ivl*n_atom_per_ivl + bs_env%atoms_i(1) - 1, &
2366 164 : bs_env%atoms_i(2))
2367 : END DO
2368 :
2369 164 : DO j_ivl = 1, n_intervals_j
2370 82 : bs_env%j_atom_intervals(1, j_ivl) = (j_ivl - 1)*n_atom_per_ivl + bs_env%atoms_j(1)
2371 : bs_env%j_atom_intervals(2, j_ivl) = MIN(j_ivl*n_atom_per_ivl + bs_env%atoms_j(1) - 1, &
2372 164 : bs_env%atoms_j(2))
2373 : END DO
2374 :
2375 328 : ALLOCATE (bs_env%skip_Sigma_occ(n_intervals_i, n_intervals_j))
2376 246 : ALLOCATE (bs_env%skip_Sigma_vir(n_intervals_i, n_intervals_j))
2377 246 : bs_env%skip_Sigma_occ(:, :) = .FALSE.
2378 246 : bs_env%skip_Sigma_vir(:, :) = .FALSE.
2379 82 : bs_env%n_skip_chi = 0
2380 :
2381 246 : ALLOCATE (bs_env%skip_chi(n_intervals_i, n_intervals_j))
2382 246 : bs_env%skip_chi(:, :) = .FALSE.
2383 82 : bs_env%n_skip_sigma = 0
2384 :
2385 : ! choose atomic range for µ and σ ("inner loop (IL) atom") in
2386 : ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
2387 : ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
2388 : n_atom_per_IL_ivl = MIN(INT(bs_env%safety_factor_memory*input_memory_per_proc &
2389 : *bs_env%group_size_tensor/n_atom_per_ivl &
2390 : /bs_env%max_AO_bf_per_atom &
2391 : /bs_env%n_RI/8/SQRT(bs_env%occupation_3c_int) &
2392 82 : /bs_env%max_AO_bf_per_atom), bs_env%n_atom)
2393 :
2394 82 : n_intervals_inner_loop_atoms = (bs_env%n_atom - 1)/n_atom_per_IL_ivl + 1
2395 :
2396 82 : bs_env%n_intervals_inner_loop_atoms = n_intervals_inner_loop_atoms
2397 :
2398 246 : ALLOCATE (bs_env%inner_loop_atom_intervals(2, n_intervals_inner_loop_atoms))
2399 164 : DO IL_ivl = 1, n_intervals_inner_loop_atoms
2400 82 : bs_env%inner_loop_atom_intervals(1, IL_ivl) = (IL_ivl - 1)*n_atom_per_IL_ivl + 1
2401 164 : bs_env%inner_loop_atom_intervals(2, IL_ivl) = MIN(IL_ivl*n_atom_per_IL_ivl, bs_env%n_atom)
2402 : END DO
2403 :
2404 82 : u = bs_env%unit_nr
2405 82 : IF (u > 0) THEN
2406 41 : WRITE (u, '(T2,A)') ''
2407 41 : WRITE (u, '(T2,A,I33)') 'Number of i and j atoms in M_λνP(τ), N_νλQ(τ):', n_atom_per_ivl
2408 41 : WRITE (u, '(T2,A,I18)') 'Number of inner loop atoms for µ in M_λνP = sum_µ (µν|P) G_µλ', &
2409 82 : n_atom_per_IL_ivl
2410 : END IF
2411 :
2412 82 : CALL timestop(handle)
2413 :
2414 82 : END SUBROUTINE set_sparsity_parallelization_parameters
2415 :
2416 : ! **************************************************************************************************
2417 : !> \brief ...
2418 : !> \param qs_env ...
2419 : !> \param bs_env ...
2420 : ! **************************************************************************************************
2421 122 : SUBROUTINE check_for_restart_files(qs_env, bs_env)
2422 : TYPE(qs_environment_type), POINTER :: qs_env
2423 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2424 :
2425 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_for_restart_files'
2426 :
2427 : CHARACTER(LEN=9) :: frmt
2428 : CHARACTER(len=default_path_length) :: f_chi, f_s_n, f_s_p, f_s_x, f_w_t, &
2429 : prefix, project_name, Z_lP_name
2430 : INTEGER :: handle, i_spin, i_t_or_w, ind, n_spin, &
2431 : num_time_freq_points
2432 : LOGICAL :: chi_exists, Sigma_neg_time_exists, &
2433 : Sigma_pos_time_exists, &
2434 : Sigma_x_spin_exists, W_time_exists, &
2435 : Z_lP_exists
2436 : TYPE(cp_logger_type), POINTER :: logger
2437 : TYPE(section_vals_type), POINTER :: input, print_key
2438 :
2439 122 : CALL timeset(routineN, handle)
2440 :
2441 122 : num_time_freq_points = bs_env%num_time_freq_points
2442 122 : n_spin = bs_env%n_spin
2443 :
2444 366 : ALLOCATE (bs_env%read_chi(num_time_freq_points))
2445 244 : ALLOCATE (bs_env%calc_chi(num_time_freq_points))
2446 488 : ALLOCATE (bs_env%Sigma_c_exists(num_time_freq_points, n_spin))
2447 :
2448 122 : CALL get_qs_env(qs_env, input=input)
2449 :
2450 122 : logger => cp_get_default_logger()
2451 122 : print_key => section_vals_get_subs_vals(input, 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART')
2452 : project_name = cp_print_key_generate_filename(logger, print_key, extension="", &
2453 122 : my_local=.FALSE.)
2454 122 : WRITE (prefix, '(2A)') TRIM(project_name), "-RESTART_"
2455 122 : bs_env%prefix = prefix
2456 :
2457 122 : bs_env%all_W_exist = .TRUE.
2458 :
2459 2088 : DO i_t_or_w = 1, num_time_freq_points
2460 :
2461 1966 : IF (i_t_or_w < 10) THEN
2462 1064 : WRITE (frmt, '(A)') '(3A,I1,A)'
2463 1064 : WRITE (f_chi, frmt) TRIM(prefix), bs_env%chi_name, "_0", i_t_or_w, ".matrix"
2464 1064 : WRITE (f_W_t, frmt) TRIM(prefix), bs_env%W_time_name, "_0", i_t_or_w, ".matrix"
2465 902 : ELSE IF (i_t_or_w < 100) THEN
2466 902 : WRITE (frmt, '(A)') '(3A,I2,A)'
2467 902 : WRITE (f_chi, frmt) TRIM(prefix), bs_env%chi_name, "_", i_t_or_w, ".matrix"
2468 902 : WRITE (f_W_t, frmt) TRIM(prefix), bs_env%W_time_name, "_", i_t_or_w, ".matrix"
2469 : ELSE
2470 0 : CPABORT('Please implement more than 99 time/frequency points.')
2471 : END IF
2472 :
2473 1966 : INQUIRE (file=TRIM(f_chi), exist=chi_exists)
2474 1966 : INQUIRE (file=TRIM(f_W_t), exist=W_time_exists)
2475 :
2476 1966 : bs_env%read_chi(i_t_or_w) = chi_exists
2477 1966 : bs_env%calc_chi(i_t_or_w) = .NOT. chi_exists
2478 :
2479 1966 : bs_env%all_W_exist = bs_env%all_W_exist .AND. W_time_exists
2480 :
2481 : ! the self-energy is spin-dependent
2482 4334 : DO i_spin = 1, n_spin
2483 :
2484 2246 : ind = i_t_or_w + (i_spin - 1)*num_time_freq_points
2485 :
2486 2246 : IF (ind < 10) THEN
2487 1064 : WRITE (frmt, '(A)') '(3A,I1,A)'
2488 1064 : WRITE (f_S_p, frmt) TRIM(prefix), bs_env%Sigma_p_name, "_0", ind, ".matrix"
2489 1064 : WRITE (f_S_n, frmt) TRIM(prefix), bs_env%Sigma_n_name, "_0", ind, ".matrix"
2490 1182 : ELSE IF (ind < 100) THEN
2491 1182 : WRITE (frmt, '(A)') '(3A,I2,A)'
2492 1182 : WRITE (f_S_p, frmt) TRIM(prefix), bs_env%Sigma_p_name, "_", ind, ".matrix"
2493 1182 : WRITE (f_S_n, frmt) TRIM(prefix), bs_env%Sigma_n_name, "_", ind, ".matrix"
2494 : ELSE
2495 0 : CPABORT('Please implement more than 99 combined spin+freq indices.')
2496 : END IF
2497 :
2498 2246 : INQUIRE (file=TRIM(f_S_p), exist=Sigma_pos_time_exists)
2499 2246 : INQUIRE (file=TRIM(f_S_n), exist=Sigma_neg_time_exists)
2500 :
2501 : bs_env%Sigma_c_exists(i_t_or_w, i_spin) = Sigma_pos_time_exists .AND. &
2502 6338 : Sigma_neg_time_exists
2503 :
2504 : END DO
2505 :
2506 : END DO
2507 :
2508 : ! Marek : In the RTBSE run, check also for zero frequency W
2509 122 : IF (bs_env%rtp_method == rtp_method_bse .OR. &
2510 : bs_env%rtp_method == rtp_method_bse_linearized) THEN
2511 72 : WRITE (f_W_t, '(3A,I1,A)') TRIM(prefix), "W_freq_rtp", "_0", 0, ".matrix"
2512 72 : INQUIRE (file=TRIM(f_W_t), exist=W_time_exists)
2513 144 : bs_env%all_W_exist = bs_env%all_W_exist .AND. W_time_exists
2514 : END IF
2515 :
2516 : ! Check for Restart Z_lP file
2517 122 : IF (bs_env%do_gw_ri_rs) THEN
2518 40 : WRITE (Z_lP_name, '(3A)') TRIM(prefix), "Z_lP", ".matrix"
2519 40 : INQUIRE (file=TRIM(Z_lP_name), exist=Z_lP_exists)
2520 40 : bs_env%ri_rs%Z_lP_exists = Z_lP_exists
2521 : END IF
2522 :
2523 122 : IF (bs_env%all_W_exist) THEN
2524 0 : bs_env%read_chi(:) = .FALSE.
2525 0 : bs_env%calc_chi(:) = .FALSE.
2526 : END IF
2527 :
2528 122 : bs_env%Sigma_x_exists = .TRUE.
2529 262 : DO i_spin = 1, n_spin
2530 140 : WRITE (f_S_x, '(3A,I1,A)') TRIM(prefix), bs_env%Sigma_x_name, "_0", i_spin, ".matrix"
2531 140 : INQUIRE (file=TRIM(f_S_x), exist=Sigma_x_spin_exists)
2532 394 : bs_env%Sigma_x_exists = bs_env%Sigma_x_exists .AND. Sigma_x_spin_exists
2533 : END DO
2534 :
2535 : ! static COHSEX: Σ^c(iω=0,k=0) restart, one file per spin; if present, χ and W are not needed
2536 122 : bs_env%Sigma_c_freq_zero_exists = bs_env%do_static_cohsex
2537 262 : DO i_spin = 1, n_spin
2538 140 : WRITE (f_S_x, '(3A,I1,A)') TRIM(prefix), bs_env%Sigma_c_freq_zero_name, "_0", i_spin, ".matrix"
2539 140 : INQUIRE (file=TRIM(f_S_x), exist=Sigma_x_spin_exists)
2540 402 : bs_env%Sigma_c_freq_zero_exists = bs_env%Sigma_c_freq_zero_exists .AND. Sigma_x_spin_exists
2541 : END DO
2542 :
2543 : ! If any restart files are read, check if the SCF converged in 1 step.
2544 : ! This is important because a re-iterated SCF can lead to spurious GW results
2545 : IF (ANY(bs_env%read_chi(:)) &
2546 : .OR. ANY(bs_env%Sigma_c_exists) &
2547 : .OR. bs_env%all_W_exist &
2548 : .OR. bs_env%Sigma_x_exists &
2549 4346 : .OR. bs_env%Sigma_c_freq_zero_exists &
2550 : ) THEN
2551 :
2552 6 : IF (qs_env%scf_env%iter_count /= 1) THEN
2553 : CALL cp_warn(__LOCATION__, "SCF needed more than 1 step, "// &
2554 6 : "which might lead to spurious GW results when using GW restart files. ")
2555 : END IF
2556 : END IF
2557 :
2558 122 : CALL timestop(handle)
2559 :
2560 122 : END SUBROUTINE check_for_restart_files
2561 :
2562 : ! **************************************************************************************************
2563 : !> \brief ...
2564 : !> \param qs_env ...
2565 : !> \param bs_env ...
2566 : ! **************************************************************************************************
2567 16 : SUBROUTINE compute_3c_integrals(qs_env, bs_env)
2568 :
2569 : TYPE(qs_environment_type), POINTER :: qs_env
2570 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2571 :
2572 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_integrals'
2573 :
2574 : INTEGER :: handle, j_cell, k_cell, nimages_3c
2575 :
2576 16 : CALL timeset(routineN, handle)
2577 :
2578 16 : nimages_3c = bs_env%nimages_3c
2579 2436 : ALLOCATE (bs_env%t_3c_int(nimages_3c, nimages_3c))
2580 194 : DO j_cell = 1, nimages_3c
2581 2260 : DO k_cell = 1, nimages_3c
2582 2244 : CALL dbt_create(bs_env%t_RI_AO__AO, bs_env%t_3c_int(j_cell, k_cell))
2583 : END DO
2584 : END DO
2585 :
2586 : CALL build_3c_integrals(bs_env%t_3c_int, &
2587 : bs_env%eps_filter, &
2588 : qs_env, &
2589 : bs_env%nl_3c, &
2590 : int_eps=bs_env%eps_filter*0.05_dp, &
2591 : basis_i=bs_env%basis_set_RI, &
2592 : basis_j=bs_env%basis_set_AO, &
2593 : basis_k=bs_env%basis_set_AO, &
2594 : potential_parameter=bs_env%ri_metric, &
2595 : desymmetrize=.FALSE., do_kpoints=.TRUE., cell_sym=.TRUE., &
2596 16 : cell_to_index_ext=bs_env%cell_to_index_3c)
2597 :
2598 16 : CALL bs_env%para_env%sync()
2599 :
2600 16 : CALL timestop(handle)
2601 :
2602 16 : END SUBROUTINE compute_3c_integrals
2603 :
2604 : ! **************************************************************************************************
2605 : !> \brief ...
2606 : !> \param bs_env ...
2607 : ! **************************************************************************************************
2608 16 : SUBROUTINE setup_cells_Delta_R(bs_env)
2609 :
2610 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2611 :
2612 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_cells_Delta_R'
2613 :
2614 : INTEGER :: handle
2615 :
2616 16 : CALL timeset(routineN, handle)
2617 :
2618 : ! cell sums batch wise for fixed ΔR = S_1 - R_1; for example:
2619 : ! Σ_λσ^R = sum_PR1νS1 M^G_λ0,νS1,PR1 M^W_σR,νS1,PR1
2620 :
2621 : CALL sum_two_R_grids(bs_env%index_to_cell_3c, &
2622 : bs_env%index_to_cell_3c, &
2623 : bs_env%nimages_3c, bs_env%nimages_3c, &
2624 : bs_env%index_to_cell_Delta_R, &
2625 : bs_env%cell_to_index_Delta_R, &
2626 16 : bs_env%nimages_Delta_R)
2627 :
2628 16 : IF (bs_env%unit_nr > 0) THEN
2629 8 : WRITE (bs_env%unit_nr, FMT="(T2,A,I61)") "Number of cells ΔR", bs_env%nimages_Delta_R
2630 : END IF
2631 :
2632 16 : CALL timestop(handle)
2633 :
2634 16 : END SUBROUTINE setup_cells_Delta_R
2635 :
2636 : ! **************************************************************************************************
2637 : !> \brief ...
2638 : !> \param index_to_cell_1 ...
2639 : !> \param index_to_cell_2 ...
2640 : !> \param nimages_1 ...
2641 : !> \param nimages_2 ...
2642 : !> \param index_to_cell ...
2643 : !> \param cell_to_index ...
2644 : !> \param nimages ...
2645 : ! **************************************************************************************************
2646 16 : SUBROUTINE sum_two_R_grids(index_to_cell_1, index_to_cell_2, nimages_1, nimages_2, &
2647 : index_to_cell, cell_to_index, nimages)
2648 :
2649 : INTEGER, DIMENSION(:, :) :: index_to_cell_1, index_to_cell_2
2650 : INTEGER :: nimages_1, nimages_2
2651 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
2652 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2653 : INTEGER :: nimages
2654 :
2655 : CHARACTER(LEN=*), PARAMETER :: routineN = 'sum_two_R_grids'
2656 :
2657 : INTEGER :: handle, i_dim, img_1, img_2, nimages_max
2658 16 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_tmp
2659 : INTEGER, DIMENSION(3) :: cell_1, cell_2, R, R_max, R_min
2660 :
2661 16 : CALL timeset(routineN, handle)
2662 :
2663 64 : DO i_dim = 1, 3
2664 1116 : R_min(i_dim) = MINVAL(index_to_cell_1(i_dim, :)) + MINVAL(index_to_cell_2(i_dim, :))
2665 1180 : R_max(i_dim) = MAXVAL(index_to_cell_1(i_dim, :)) + MAXVAL(index_to_cell_2(i_dim, :))
2666 : END DO
2667 :
2668 16 : nimages_max = (R_max(1) - R_min(1) + 1)*(R_max(2) - R_min(2) + 1)*(R_max(3) - R_min(3) + 1)
2669 :
2670 48 : ALLOCATE (index_to_cell_tmp(3, nimages_max))
2671 2176 : index_to_cell_tmp(:, :) = -1
2672 :
2673 80 : ALLOCATE (cell_to_index(R_min(1):R_max(1), R_min(2):R_max(2), R_min(3):R_max(3)))
2674 1156 : cell_to_index(:, :, :) = -1
2675 :
2676 16 : nimages = 0
2677 :
2678 194 : DO img_1 = 1, nimages_1
2679 :
2680 2260 : DO img_2 = 1, nimages_2
2681 :
2682 8264 : cell_1(1:3) = index_to_cell_1(1:3, img_1)
2683 8264 : cell_2(1:3) = index_to_cell_2(1:3, img_2)
2684 :
2685 8264 : R(1:3) = cell_1(1:3) + cell_2(1:3)
2686 :
2687 : ! check whether we have found a new cell
2688 2244 : IF (cell_to_index(R(1), R(2), R(3)) == -1) THEN
2689 :
2690 516 : nimages = nimages + 1
2691 516 : cell_to_index(R(1), R(2), R(3)) = nimages
2692 2064 : index_to_cell_tmp(1:3, nimages) = R(1:3)
2693 :
2694 : END IF
2695 :
2696 : END DO
2697 :
2698 : END DO
2699 :
2700 48 : ALLOCATE (index_to_cell(3, nimages))
2701 2080 : index_to_cell(:, :) = index_to_cell_tmp(1:3, 1:nimages)
2702 :
2703 16 : CALL timestop(handle)
2704 :
2705 32 : END SUBROUTINE sum_two_R_grids
2706 :
2707 : ! **************************************************************************************************
2708 : !> \brief ...
2709 : !> \param bs_env ...
2710 : ! **************************************************************************************************
2711 16 : SUBROUTINE setup_parallelization_Delta_R(bs_env)
2712 :
2713 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2714 :
2715 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_parallelization_Delta_R'
2716 :
2717 : INTEGER :: handle, i_cell_Delta_R, i_task_local, &
2718 : n_tasks_local
2719 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: i_cell_Delta_R_group, &
2720 16 : n_tensor_ops_Delta_R
2721 :
2722 16 : CALL timeset(routineN, handle)
2723 :
2724 16 : CALL compute_n_tensor_ops_Delta_R(bs_env, n_tensor_ops_Delta_R)
2725 :
2726 16 : CALL compute_Delta_R_dist(bs_env, n_tensor_ops_Delta_R, i_cell_Delta_R_group, n_tasks_local)
2727 :
2728 16 : bs_env%n_tasks_Delta_R_local = n_tasks_local
2729 :
2730 48 : ALLOCATE (bs_env%task_Delta_R(n_tasks_local))
2731 :
2732 16 : i_task_local = 0
2733 532 : DO i_cell_Delta_R = 1, bs_env%nimages_Delta_R
2734 :
2735 516 : IF (i_cell_Delta_R_group(i_cell_Delta_R) /= bs_env%tensor_group_color) CYCLE
2736 :
2737 223 : i_task_local = i_task_local + 1
2738 :
2739 532 : bs_env%task_Delta_R(i_task_local) = i_cell_Delta_R
2740 :
2741 : END DO
2742 :
2743 32 : ALLOCATE (bs_env%skip_DR_chi(n_tasks_local))
2744 239 : bs_env%skip_DR_chi(:) = .FALSE.
2745 32 : ALLOCATE (bs_env%skip_DR_Sigma(n_tasks_local))
2746 239 : bs_env%skip_DR_Sigma(:) = .FALSE.
2747 :
2748 16 : CALL allocate_skip_3xR(bs_env%skip_DR_R12_S_Goccx3c_chi, bs_env)
2749 16 : CALL allocate_skip_3xR(bs_env%skip_DR_R12_S_Gvirx3c_chi, bs_env)
2750 16 : CALL allocate_skip_3xR(bs_env%skip_DR_R_R2_MxM_chi, bs_env)
2751 :
2752 16 : CALL allocate_skip_3xR(bs_env%skip_DR_R1_S2_Gx3c_Sigma, bs_env)
2753 16 : CALL allocate_skip_3xR(bs_env%skip_DR_R1_R_MxM_Sigma, bs_env)
2754 :
2755 16 : CALL timestop(handle)
2756 :
2757 32 : END SUBROUTINE setup_parallelization_Delta_R
2758 :
2759 : ! **************************************************************************************************
2760 : !> \brief ...
2761 : !> \param bs_env ...
2762 : !> \param n_tensor_ops_Delta_R ...
2763 : ! **************************************************************************************************
2764 16 : SUBROUTINE compute_n_tensor_ops_Delta_R(bs_env, n_tensor_ops_Delta_R)
2765 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2766 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_tensor_ops_Delta_R
2767 :
2768 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_n_tensor_ops_Delta_R'
2769 :
2770 : INTEGER :: handle, i_cell_Delta_R, i_cell_R, i_cell_R1, i_cell_R1_minus_R, i_cell_R2, &
2771 : i_cell_R2_m_R1, i_cell_S1, i_cell_S1_m_R1_p_R2, i_cell_S1_minus_R, i_cell_S2, &
2772 : nimages_Delta_R
2773 : INTEGER, DIMENSION(3) :: cell_DR, cell_m_R1, cell_R, cell_R1, cell_R1_minus_R, cell_R2, &
2774 : cell_R2_m_R1, cell_S1, cell_S1_m_R2_p_R1, cell_S1_minus_R, cell_S1_p_S2_m_R1, cell_S2
2775 : LOGICAL :: cell_found
2776 :
2777 16 : CALL timeset(routineN, handle)
2778 :
2779 16 : nimages_Delta_R = bs_env%nimages_Delta_R
2780 :
2781 48 : ALLOCATE (n_tensor_ops_Delta_R(nimages_Delta_R))
2782 16 : n_tensor_ops_Delta_R(:) = 0
2783 :
2784 : ! compute number of tensor operations for specific Delta_R
2785 532 : DO i_cell_Delta_R = 1, nimages_Delta_R
2786 :
2787 516 : IF (MODULO(i_cell_Delta_R, bs_env%num_tensor_groups) /= bs_env%tensor_group_color) CYCLE
2788 :
2789 3279 : DO i_cell_R1 = 1, bs_env%nimages_3c
2790 :
2791 12020 : cell_R1(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_R1)
2792 12020 : cell_DR(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_Delta_R)
2793 :
2794 : ! S_1 = R_1 + ΔR (from ΔR = S_1 - R_1)
2795 : CALL add_R(cell_R1, cell_DR, bs_env%index_to_cell_3c, cell_S1, &
2796 3005 : cell_found, bs_env%cell_to_index_3c, i_cell_S1)
2797 3005 : IF (.NOT. cell_found) CYCLE
2798 :
2799 9700 : DO i_cell_R2 = 1, bs_env%nimages_scf_desymm
2800 :
2801 34920 : cell_R2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_R2)
2802 :
2803 : ! R_2 - R_1
2804 : CALL add_R(cell_R2, -cell_R1, bs_env%index_to_cell_3c, cell_R2_m_R1, &
2805 34920 : cell_found, bs_env%cell_to_index_3c, i_cell_R2_m_R1)
2806 8730 : IF (.NOT. cell_found) CYCLE
2807 :
2808 : ! S_1 - R_1 + R_2
2809 : CALL add_R(cell_S1, cell_R2_m_R1, bs_env%index_to_cell_3c, cell_S1_m_R2_p_R1, &
2810 5250 : cell_found, bs_env%cell_to_index_3c, i_cell_S1_m_R1_p_R2)
2811 5250 : IF (.NOT. cell_found) CYCLE
2812 :
2813 13208 : n_tensor_ops_Delta_R(i_cell_Delta_R) = n_tensor_ops_Delta_R(i_cell_Delta_R) + 1
2814 :
2815 : END DO ! i_cell_R2
2816 :
2817 9700 : DO i_cell_S2 = 1, bs_env%nimages_scf_desymm
2818 :
2819 34920 : cell_S2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_S2)
2820 34920 : cell_m_R1(1:3) = -cell_R1(1:3)
2821 34920 : cell_S1_p_S2_m_R1(1:3) = cell_S1(1:3) + cell_S2(1:3) - cell_R1(1:3)
2822 :
2823 8730 : CALL is_cell_in_index_to_cell(cell_m_R1, bs_env%index_to_cell_3c, cell_found)
2824 8730 : IF (.NOT. cell_found) CYCLE
2825 :
2826 7029 : CALL is_cell_in_index_to_cell(cell_S1_p_S2_m_R1, bs_env%index_to_cell_3c, cell_found)
2827 970 : IF (.NOT. cell_found) CYCLE
2828 :
2829 : END DO ! i_cell_S2
2830 :
2831 13221 : DO i_cell_R = 1, bs_env%nimages_scf_desymm
2832 :
2833 34920 : cell_R = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_R)
2834 :
2835 : ! R_1 - R
2836 : CALL add_R(cell_R1, -cell_R, bs_env%index_to_cell_3c, cell_R1_minus_R, &
2837 34920 : cell_found, bs_env%cell_to_index_3c, i_cell_R1_minus_R)
2838 8730 : IF (.NOT. cell_found) CYCLE
2839 :
2840 : ! S_1 - R
2841 : CALL add_R(cell_S1, -cell_R, bs_env%index_to_cell_3c, cell_S1_minus_R, &
2842 22764 : cell_found, bs_env%cell_to_index_3c, i_cell_S1_minus_R)
2843 3005 : IF (.NOT. cell_found) CYCLE
2844 :
2845 : END DO ! i_cell_R
2846 :
2847 : END DO ! i_cell_R1
2848 :
2849 : END DO ! i_cell_Delta_R
2850 :
2851 16 : CALL bs_env%para_env%sum(n_tensor_ops_Delta_R)
2852 :
2853 16 : CALL timestop(handle)
2854 :
2855 16 : END SUBROUTINE compute_n_tensor_ops_Delta_R
2856 :
2857 : ! **************************************************************************************************
2858 : !> \brief ...
2859 : !> \param cell_1 ...
2860 : !> \param cell_2 ...
2861 : !> \param index_to_cell ...
2862 : !> \param cell_1_plus_2 ...
2863 : !> \param cell_found ...
2864 : !> \param cell_to_index ...
2865 : !> \param i_cell_1_plus_2 ...
2866 : ! **************************************************************************************************
2867 279284 : SUBROUTINE add_R(cell_1, cell_2, index_to_cell, cell_1_plus_2, cell_found, &
2868 : cell_to_index, i_cell_1_plus_2)
2869 :
2870 : INTEGER, DIMENSION(3) :: cell_1, cell_2
2871 : INTEGER, DIMENSION(:, :) :: index_to_cell
2872 : INTEGER, DIMENSION(3) :: cell_1_plus_2
2873 : LOGICAL :: cell_found
2874 : INTEGER, DIMENSION(:, :, :), INTENT(IN), &
2875 : OPTIONAL, POINTER :: cell_to_index
2876 : INTEGER, INTENT(OUT), OPTIONAL :: i_cell_1_plus_2
2877 :
2878 : CHARACTER(LEN=*), PARAMETER :: routineN = 'add_R'
2879 :
2880 : INTEGER :: handle
2881 :
2882 279284 : CALL timeset(routineN, handle)
2883 :
2884 1117136 : cell_1_plus_2(1:3) = cell_1(1:3) + cell_2(1:3)
2885 :
2886 279284 : CALL is_cell_in_index_to_cell(cell_1_plus_2, index_to_cell, cell_found)
2887 :
2888 279284 : IF (PRESENT(i_cell_1_plus_2)) THEN
2889 279284 : IF (cell_found) THEN
2890 160334 : CPASSERT(PRESENT(cell_to_index))
2891 160334 : i_cell_1_plus_2 = cell_to_index(cell_1_plus_2(1), cell_1_plus_2(2), cell_1_plus_2(3))
2892 : ELSE
2893 118950 : i_cell_1_plus_2 = -1000
2894 : END IF
2895 : END IF
2896 :
2897 279284 : CALL timestop(handle)
2898 :
2899 279284 : END SUBROUTINE add_R
2900 :
2901 : ! **************************************************************************************************
2902 : !> \brief ...
2903 : !> \param cell ...
2904 : !> \param index_to_cell ...
2905 : !> \param cell_found ...
2906 : ! **************************************************************************************************
2907 437679 : SUBROUTINE is_cell_in_index_to_cell(cell, index_to_cell, cell_found)
2908 : INTEGER, DIMENSION(3) :: cell
2909 : INTEGER, DIMENSION(:, :) :: index_to_cell
2910 : LOGICAL :: cell_found
2911 :
2912 : CHARACTER(LEN=*), PARAMETER :: routineN = 'is_cell_in_index_to_cell'
2913 :
2914 : INTEGER :: handle, i_cell, nimg
2915 : INTEGER, DIMENSION(3) :: cell_i
2916 :
2917 437679 : CALL timeset(routineN, handle)
2918 :
2919 437679 : nimg = SIZE(index_to_cell, 2)
2920 :
2921 437679 : cell_found = .FALSE.
2922 :
2923 5604294 : DO i_cell = 1, nimg
2924 :
2925 20666460 : cell_i(1:3) = index_to_cell(1:3, i_cell)
2926 :
2927 5604294 : IF (cell_i(1) == cell(1) .AND. cell_i(2) == cell(2) .AND. cell_i(3) == cell(3)) THEN
2928 263631 : cell_found = .TRUE.
2929 : END IF
2930 :
2931 : END DO
2932 :
2933 437679 : CALL timestop(handle)
2934 :
2935 437679 : END SUBROUTINE is_cell_in_index_to_cell
2936 :
2937 : ! **************************************************************************************************
2938 : !> \brief ...
2939 : !> \param bs_env ...
2940 : !> \param n_tensor_ops_Delta_R ...
2941 : !> \param i_cell_Delta_R_group ...
2942 : !> \param n_tasks_local ...
2943 : ! **************************************************************************************************
2944 16 : SUBROUTINE compute_Delta_R_dist(bs_env, n_tensor_ops_Delta_R, i_cell_Delta_R_group, n_tasks_local)
2945 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2946 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_tensor_ops_Delta_R, &
2947 : i_cell_Delta_R_group
2948 : INTEGER :: n_tasks_local
2949 :
2950 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Delta_R_dist'
2951 :
2952 : INTEGER :: handle, i_Delta_R_max_op, i_group_min, &
2953 : nimages_Delta_R, u
2954 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_tensor_ops_Delta_R_in_group
2955 :
2956 16 : CALL timeset(routineN, handle)
2957 :
2958 16 : nimages_Delta_R = bs_env%nimages_Delta_R
2959 :
2960 16 : u = bs_env%unit_nr
2961 :
2962 16 : IF (u > 0 .AND. nimages_Delta_R < bs_env%num_tensor_groups) THEN
2963 0 : WRITE (u, FMT="(T2,A,I5,A,I5,A)") "There are only ", nimages_Delta_R, &
2964 0 : " tasks to work on but there are ", bs_env%num_tensor_groups, " groups."
2965 0 : WRITE (u, FMT="(T2,A)") "Please reduce the number of MPI processes."
2966 0 : WRITE (u, '(T2,A)') ''
2967 : END IF
2968 :
2969 48 : ALLOCATE (n_tensor_ops_Delta_R_in_group(bs_env%num_tensor_groups))
2970 16 : n_tensor_ops_Delta_R_in_group(:) = 0
2971 48 : ALLOCATE (i_cell_Delta_R_group(nimages_Delta_R))
2972 532 : i_cell_Delta_R_group(:) = -1
2973 :
2974 16 : n_tasks_local = 0
2975 :
2976 1914 : DO WHILE (ANY(n_tensor_ops_Delta_R(:) /= 0))
2977 :
2978 : ! get largest element of n_tensor_ops_Delta_R
2979 15484 : i_Delta_R_max_op = MAXLOC(n_tensor_ops_Delta_R, 1)
2980 :
2981 : ! distribute i_Delta_R_max_op to tensor group which has currently the smallest load
2982 1784 : i_group_min = MINLOC(n_tensor_ops_Delta_R_in_group, 1)
2983 :
2984 : ! the tensor groups are 0-index based; but i_group_min is 1-index based
2985 446 : i_cell_Delta_R_group(i_Delta_R_max_op) = i_group_min - 1
2986 : n_tensor_ops_Delta_R_in_group(i_group_min) = n_tensor_ops_Delta_R_in_group(i_group_min) + &
2987 446 : n_tensor_ops_Delta_R(i_Delta_R_max_op)
2988 :
2989 : ! remove i_Delta_R_max_op from n_tensor_ops_Delta_R
2990 446 : n_tensor_ops_Delta_R(i_Delta_R_max_op) = 0
2991 :
2992 462 : IF (bs_env%tensor_group_color == i_group_min - 1) n_tasks_local = n_tasks_local + 1
2993 :
2994 : END DO
2995 :
2996 16 : CALL timestop(handle)
2997 :
2998 32 : END SUBROUTINE compute_Delta_R_dist
2999 :
3000 : ! **************************************************************************************************
3001 : !> \brief ...
3002 : !> \param skip ...
3003 : !> \param bs_env ...
3004 : ! **************************************************************************************************
3005 80 : SUBROUTINE allocate_skip_3xR(skip, bs_env)
3006 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: skip
3007 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3008 :
3009 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_skip_3xR'
3010 :
3011 : INTEGER :: handle
3012 :
3013 80 : CALL timeset(routineN, handle)
3014 :
3015 400 : ALLOCATE (skip(bs_env%n_tasks_Delta_R_local, bs_env%nimages_3c, bs_env%nimages_scf_desymm))
3016 80 : skip(:, :, :) = .FALSE.
3017 :
3018 80 : CALL timestop(handle)
3019 :
3020 80 : END SUBROUTINE allocate_skip_3xR
3021 :
3022 : ! **************************************************************************************************
3023 : !> \brief ...
3024 : !> \param qs_env ...
3025 : !> \param bs_env ...
3026 : ! **************************************************************************************************
3027 16 : SUBROUTINE allocate_matrices_small_cell_full_kp_tensor(qs_env, bs_env)
3028 : TYPE(qs_environment_type), POINTER :: qs_env
3029 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3030 :
3031 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_matrices_small_cell_full_kp_tensor'
3032 :
3033 : INTEGER :: handle, i_spin, i_t, img, n_spin, &
3034 : nimages_scf, num_time_freq_points
3035 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
3036 : TYPE(mp_para_env_type), POINTER :: para_env
3037 :
3038 16 : CALL timeset(routineN, handle)
3039 :
3040 16 : nimages_scf = bs_env%nimages_scf_desymm
3041 16 : num_time_freq_points = bs_env%num_time_freq_points
3042 16 : n_spin = bs_env%n_spin
3043 :
3044 16 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
3045 :
3046 192 : ALLOCATE (bs_env%fm_G_S(nimages_scf))
3047 176 : ALLOCATE (bs_env%fm_Sigma_x_R(nimages_scf))
3048 1104 : ALLOCATE (bs_env%fm_chi_R_t(nimages_scf, num_time_freq_points))
3049 1088 : ALLOCATE (bs_env%fm_MWM_R_t(nimages_scf, num_time_freq_points))
3050 1136 : ALLOCATE (bs_env%fm_Sigma_c_R_neg_tau(nimages_scf, num_time_freq_points, n_spin))
3051 1120 : ALLOCATE (bs_env%fm_Sigma_c_R_pos_tau(nimages_scf, num_time_freq_points, n_spin))
3052 160 : DO img = 1, nimages_scf
3053 144 : CALL cp_fm_create(bs_env%fm_G_S(img), bs_env%fm_work_mo(1)%matrix_struct)
3054 144 : CALL cp_fm_create(bs_env%fm_Sigma_x_R(img), bs_env%fm_work_mo(1)%matrix_struct)
3055 1096 : DO i_t = 1, num_time_freq_points
3056 936 : CALL cp_fm_create(bs_env%fm_chi_R_t(img, i_t), bs_env%fm_RI_RI%matrix_struct)
3057 936 : CALL cp_fm_create(bs_env%fm_MWM_R_t(img, i_t), bs_env%fm_RI_RI%matrix_struct)
3058 936 : CALL cp_fm_set_all(bs_env%fm_MWM_R_t(img, i_t), 0.0_dp)
3059 2016 : DO i_spin = 1, n_spin
3060 : CALL cp_fm_create(bs_env%fm_Sigma_c_R_neg_tau(img, i_t, i_spin), &
3061 936 : bs_env%fm_work_mo(1)%matrix_struct)
3062 : CALL cp_fm_create(bs_env%fm_Sigma_c_R_pos_tau(img, i_t, i_spin), &
3063 936 : bs_env%fm_work_mo(1)%matrix_struct)
3064 936 : CALL cp_fm_set_all(bs_env%fm_Sigma_c_R_neg_tau(img, i_t, i_spin), 0.0_dp)
3065 1872 : CALL cp_fm_set_all(bs_env%fm_Sigma_c_R_pos_tau(img, i_t, i_spin), 0.0_dp)
3066 : END DO
3067 : END DO
3068 : END DO
3069 :
3070 16 : CALL timestop(handle)
3071 :
3072 16 : END SUBROUTINE allocate_matrices_small_cell_full_kp_tensor
3073 :
3074 : ! **************************************************************************************************
3075 : !> \brief ...
3076 : !> \param qs_env ...
3077 : !> \param bs_env ...
3078 : ! **************************************************************************************************
3079 16 : SUBROUTINE trafo_V_xc_R_to_kp(qs_env, bs_env)
3080 : TYPE(qs_environment_type), POINTER :: qs_env
3081 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3082 :
3083 : CHARACTER(LEN=*), PARAMETER :: routineN = 'trafo_V_xc_R_to_kp'
3084 :
3085 : INTEGER :: handle, ikp, img, ispin, n_ao
3086 16 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf
3087 : TYPE(cp_cfm_type) :: cfm_mo_coeff, cfm_V_xc
3088 : TYPE(cp_fm_type) :: fm_V_xc_re
3089 16 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
3090 : TYPE(kpoint_type), POINTER :: kpoints_scf
3091 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3092 16 : POINTER :: sab_nl
3093 :
3094 16 : CALL timeset(routineN, handle)
3095 :
3096 16 : n_ao = bs_env%n_ao
3097 :
3098 16 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, kpoints=kpoints_scf)
3099 :
3100 16 : NULLIFY (sab_nl)
3101 16 : CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
3102 :
3103 16 : CALL cp_cfm_create(cfm_V_xc, bs_env%cfm_work_mo%matrix_struct)
3104 16 : CALL cp_cfm_create(cfm_mo_coeff, bs_env%cfm_work_mo%matrix_struct)
3105 16 : CALL cp_fm_create(fm_V_xc_re, bs_env%cfm_work_mo%matrix_struct)
3106 :
3107 544 : DO img = 1, bs_env%nimages_scf
3108 1072 : DO ispin = 1, bs_env%n_spin
3109 : ! JW kind of hack because the format of matrix_ks remains dubious...
3110 528 : CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
3111 1056 : CALL copy_fm_to_dbcsr(bs_env%fm_V_xc_R(img, ispin), matrix_ks(ispin, img)%matrix)
3112 : END DO
3113 : END DO
3114 :
3115 80 : ALLOCATE (bs_env%v_xc_n(n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
3116 :
3117 32 : DO ispin = 1, bs_env%n_spin
3118 350 : DO ikp = 1, bs_env%nkp_bs_and_DOS
3119 :
3120 : ! v^xc^R -> v^xc(k) (matrix_ks stores v^xc^R, see SUBROUTINE compute_V_xc)
3121 : CALL rsmat_to_kp(matrix_ks, ispin, bs_env%kpoints_DOS%xkp(1:3, ikp), &
3122 318 : cell_to_index_scf, sab_nl, bs_env, cfm_V_xc)
3123 :
3124 : ! get C_µn(k)
3125 318 : CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mo_coeff)
3126 :
3127 : ! v^xc_nm(k_i) = sum_µν C^*_µn(k_i) v^xc_µν(k_i) C_νn(k_i)
3128 318 : CALL cfm_contract_ABA(cfm_mo_coeff, cfm_V_xc)
3129 :
3130 : ! get v^xc_nn(k_i) which is a real quantity as v^xc is Hermitian
3131 318 : CALL cp_cfm_to_fm(cfm_V_xc, fm_V_xc_re)
3132 334 : CALL cp_fm_get_diag(fm_V_xc_re, bs_env%v_xc_n(:, ikp, ispin))
3133 :
3134 : END DO
3135 :
3136 : END DO
3137 :
3138 : ! just rebuild the overwritten KS matrix again
3139 16 : CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
3140 :
3141 16 : CALL cp_cfm_release(cfm_V_xc)
3142 16 : CALL cp_cfm_release(cfm_mo_coeff)
3143 16 : CALL cp_fm_release(fm_V_xc_re)
3144 :
3145 16 : CALL timestop(handle)
3146 :
3147 32 : END SUBROUTINE trafo_V_xc_R_to_kp
3148 :
3149 : ! **************************************************************************************************
3150 : !> \brief ...
3151 : !> \param qs_env ...
3152 : !> \param bs_env ...
3153 : ! **************************************************************************************************
3154 16 : SUBROUTINE heuristic_RI_regularization(qs_env, bs_env)
3155 : TYPE(qs_environment_type), POINTER :: qs_env
3156 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3157 :
3158 : CHARACTER(LEN=*), PARAMETER :: routineN = 'heuristic_RI_regularization'
3159 :
3160 16 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: M
3161 : INTEGER :: handle, ikp, ikp_local, n_RI, nkp, &
3162 : nkp_local, u
3163 : REAL(KIND=dp) :: cond_nr, cond_nr_max, max_ev, &
3164 : max_ev_ikp, min_ev, min_ev_ikp
3165 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: M_R
3166 :
3167 16 : CALL timeset(routineN, handle)
3168 :
3169 : ! compute M^R_PQ = <phi_P,0|V^tr(rc)|phi_Q,R> for RI metric
3170 16 : CALL get_V_tr_R(M_R, bs_env%ri_metric, 0.0_dp, bs_env, qs_env)
3171 :
3172 16 : nkp = bs_env%nkp_chi_eps_W_orig_plus_extra
3173 16 : n_RI = bs_env%n_RI
3174 :
3175 16 : nkp_local = 0
3176 10256 : DO ikp = 1, nkp
3177 : ! trivial parallelization over k-points
3178 10240 : IF (MODULO(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
3179 10256 : nkp_local = nkp_local + 1
3180 : END DO
3181 :
3182 80 : ALLOCATE (M(n_RI, n_RI, nkp_local))
3183 :
3184 16 : ikp_local = 0
3185 16 : cond_nr_max = 0.0_dp
3186 16 : min_ev = 1000.0_dp
3187 16 : max_ev = -1000.0_dp
3188 :
3189 10256 : DO ikp = 1, nkp
3190 :
3191 : ! trivial parallelization
3192 10240 : IF (MODULO(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
3193 :
3194 5120 : ikp_local = ikp_local + 1
3195 :
3196 : ! M(k) = sum_R e^ikR M^R
3197 : CALL rs_to_kp(M_R, M(:, :, ikp_local), &
3198 : bs_env%kpoints_scf_desymm%index_to_cell, &
3199 5120 : bs_env%kpoints_chi_eps_W%xkp(1:3, ikp))
3200 :
3201 : ! compute condition number of M_PQ(k)
3202 5120 : CALL local_complex_power(M(:, :, ikp_local), 1.0_dp, 0.0_dp, cond_nr, min_ev_ikp, max_ev_ikp)
3203 :
3204 5120 : IF (cond_nr > cond_nr_max) cond_nr_max = cond_nr
3205 5120 : IF (max_ev_ikp > max_ev) max_ev = max_ev_ikp
3206 5136 : IF (min_ev_ikp < min_ev) min_ev = min_ev_ikp
3207 :
3208 : END DO ! ikp
3209 :
3210 16 : CALL bs_env%para_env%max(cond_nr_max)
3211 16 : CALL bs_env%para_env%min(min_ev)
3212 16 : CALL bs_env%para_env%max(max_ev)
3213 :
3214 16 : u = bs_env%unit_nr
3215 16 : IF (u > 0) THEN
3216 8 : WRITE (u, FMT="(T2,A,ES34.1)") "Min. abs. eigenvalue of RI metric matrix M(k)", min_ev
3217 8 : WRITE (u, FMT="(T2,A,ES34.1)") "Max. abs. eigenvalue of RI metric matrix M(k)", max_ev
3218 8 : WRITE (u, FMT="(T2,A,ES50.1)") "Max. condition number of M(k)", cond_nr_max
3219 : END IF
3220 :
3221 16 : CALL timestop(handle)
3222 :
3223 32 : END SUBROUTINE heuristic_RI_regularization
3224 :
3225 : ! **************************************************************************************************
3226 : !> \brief ...
3227 : !> \param bs_env ...
3228 : ! **************************************************************************************************
3229 138 : SUBROUTINE setup_time_and_frequency_minimax_grid(bs_env)
3230 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3231 :
3232 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_time_and_frequency_minimax_grid'
3233 :
3234 : INTEGER :: handle, homo, ispin, n_mo, n_top, &
3235 : num_time_freq_points, u
3236 : REAL(KIND=dp) :: E_max, E_max_ispin, E_min, E_min_ispin, &
3237 : E_range, max_error_min
3238 :
3239 138 : CALL timeset(routineN, handle)
3240 :
3241 138 : n_mo = bs_env%n_ao
3242 138 : num_time_freq_points = bs_env%num_time_freq_points
3243 :
3244 : ! minimum and maximum difference between eigenvalues of unoccupied and an occupied MOs
3245 138 : E_min = 1000.0_dp
3246 138 : E_max = -1000.0_dp
3247 294 : DO ispin = 1, bs_env%n_spin
3248 156 : homo = bs_env%n_occ(ispin)
3249 :
3250 : ! Highest index that is a real state. The canonical orthogonalization removes the
3251 : ! linearly dependent basis modes and parks them at the top of the spectrum with a
3252 : ! placeholder eigenvalue; n_mo_retained is how many real states it kept.
3253 156 : n_top = MAX(MIN(bs_env%n_mo_retained, n_mo), homo + 1)
3254 :
3255 296 : SELECT CASE (bs_env%gw_implementation)
3256 : CASE (tensor_large_cell_Gamma, ri_rs_large_cell_Gamma, ri_rs_non_periodic, tensor_static_cohsex)
3257 : E_min_ispin = bs_env%eigenval_scf_Gamma(homo + 1, ispin) - &
3258 140 : bs_env%eigenval_scf_Gamma(homo, ispin)
3259 : E_max_ispin = bs_env%eigenval_scf_Gamma(n_top, ispin) - &
3260 140 : bs_env%eigenval_scf_Gamma(1, ispin)
3261 : CASE (tensor_small_cell_full_kp)
3262 : E_min_ispin = MINVAL(bs_env%eigenval_scf(homo + 1, :, ispin)) - &
3263 652 : MAXVAL(bs_env%eigenval_scf(homo, :, ispin))
3264 : E_max_ispin = MAXVAL(bs_env%eigenval_scf(n_top, :, ispin)) - &
3265 808 : MINVAL(bs_env%eigenval_scf(1, :, ispin))
3266 : END SELECT
3267 156 : E_min = MIN(E_min, E_min_ispin)
3268 294 : E_max = MAX(E_max, E_max_ispin)
3269 : END DO
3270 :
3271 : ! Open-shell uses ONE minimax grid for the combined [min gap, max span] over both spins (the
3272 : ! superset covers each channel, so it is accurate; per-spin grids would only be more efficient).
3273 138 : IF (bs_env%n_spin > 1) THEN
3274 : CALL cp_hint(__LOCATION__, &
3275 : "Open-shell GW uses one minimax grid spanning [min gap, max span] across both "// &
3276 : "spin channels; raise NUM_TIME_FREQ_POINTS if QP convergence is marginal for "// &
3277 18 : "strongly spin-asymmetric systems.")
3278 : END IF
3279 :
3280 138 : E_range = E_max/E_min
3281 :
3282 : CALL build_minimax_time_frequency_grid(num_time_freq_points, E_min, E_max, bs_env%regularization_minimax, &
3283 : bs_env%num_points_per_magnitude, bs_env%time_frequency_grid, &
3284 : build_frequency=.TRUE., build_time=.TRUE., build_transforms=.TRUE., &
3285 : build_sine=.TRUE., time_scaling=2.0_dp, time_weight_scaling=1.0_dp, &
3286 : max_fit_error=max_error_min, print_warning=.FALSE., unit_nr=0, &
3287 138 : prefer_external_backend=.FALSE.)
3288 :
3289 : ! determine number of fit points in the interval [0,ω_max] for virt, or [-ω_max,0] for occ
3290 2208 : bs_env%num_freq_points_fit = COUNT(bs_env%time_frequency_grid%frequency < bs_env%freq_max_fit)
3291 :
3292 : ! iω values for the analytic continuation Σ^c_n(iω,k) -> Σ^c_n(ϵ,k)
3293 414 : ALLOCATE (bs_env%imag_freq_points_fit(bs_env%num_freq_points_fit))
3294 : bs_env%imag_freq_points_fit(:) = PACK(bs_env%time_frequency_grid%frequency, &
3295 2612 : bs_env%time_frequency_grid%frequency < bs_env%freq_max_fit)
3296 :
3297 : ! reset the number of Padé parameters if smaller than the number of
3298 : ! imaginary-frequency points for the fit
3299 138 : IF (bs_env%num_freq_points_fit < bs_env%nparam_pade) THEN
3300 138 : bs_env%nparam_pade = bs_env%num_freq_points_fit
3301 : END IF
3302 :
3303 138 : u = bs_env%unit_nr
3304 138 : IF (u > 0) THEN
3305 69 : WRITE (u, '(T2,A)') ''
3306 69 : WRITE (u, '(T2,A,F55.2)') 'SCF direct band gap (eV)', E_min*evolt
3307 69 : WRITE (u, '(T2,A,F53.2)') 'Max. SCF eigval diff. (eV)', E_max*evolt
3308 69 : WRITE (u, '(T2,A,F55.2)') 'E-Range for minimax grid', E_range
3309 69 : WRITE (u, '(T2,A,I27)') 'Number of Padé parameters for analytic continuation:', &
3310 138 : bs_env%nparam_pade
3311 69 : WRITE (u, '(T2,A)') ''
3312 : END IF
3313 :
3314 : ! in minimax grids, Fourier transforms t -> w and w -> t are split using
3315 : ! e^(iwt) = cos(wt) + i sin(wt); we thus calculate weights for trafos with a cos and
3316 : ! sine prefactor; details in Azizi, Wilhelm, Golze, Giantomassi, Panades-Barrueta,
3317 : ! Rinke, Draxl, Gonze et al., 2 publications
3318 :
3319 138 : CALL timestop(handle)
3320 :
3321 138 : END SUBROUTINE setup_time_and_frequency_minimax_grid
3322 :
3323 : ! **************************************************************************************************
3324 : !> \brief Releases the memory-heavy GW intermediates that cannot be freed in bs_env_release,
3325 : !> retaining the 3c neighbor list only when an AO-RI RT-BSE self-energy still needs it
3326 : !> \param qs_env ...
3327 : !> \param bs_env ...
3328 : ! **************************************************************************************************
3329 138 : SUBROUTINE de_init_bs_env(qs_env, bs_env)
3330 : TYPE(qs_environment_type), POINTER :: qs_env
3331 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3332 :
3333 : CHARACTER(LEN=*), PARAMETER :: routineN = 'de_init_bs_env'
3334 :
3335 : INTEGER :: handle
3336 : LOGICAL :: retain_nl_3c, rirs_kernel
3337 :
3338 138 : CALL timeset(routineN, handle)
3339 : ! deallocate quantities here which:
3340 : ! 1. cannot be deallocated in bs_env_release due to circular dependencies
3341 : ! 2. consume a lot of memory and should not be kept until the quantity is
3342 : ! deallocated in bs_env_release
3343 :
3344 : ! nl_3c feeds compute_3c_integrals after the GW: the AO-RI RT-BSE SEX self-energy (the RI-RS
3345 : ! path never reads it, and AO-RI Hartree builds its own blocks) and PROPERTIES%BSE (keep_3c_nl).
3346 138 : retain_nl_3c = bs_env%keep_3c_nl
3347 138 : IF (ASSOCIATED(bs_env%nl_3c%ij_list) .AND. (bs_env%rtp_method == rtp_method_bse)) THEN
3348 72 : CALL rtbse_resolve_rirs_flag(qs_env, bs_env, rirs_kernel=rirs_kernel)
3349 72 : retain_nl_3c = retain_nl_3c .OR. .NOT. rirs_kernel
3350 : END IF
3351 :
3352 66 : IF (retain_nl_3c) THEN
3353 56 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, *) "Keep nl_3c for BSE"
3354 : ELSE
3355 82 : CALL neighbor_list_3c_destroy(bs_env%nl_3c)
3356 : END IF
3357 :
3358 138 : CALL cp_libint_static_cleanup()
3359 :
3360 138 : CALL timestop(handle)
3361 :
3362 138 : END SUBROUTINE de_init_bs_env
3363 :
3364 : ! **************************************************************************************************
3365 : !> \brief Resolve the linRTBSE RI-RS kernel switch from the KERNEL_RI input and the GW default.
3366 : !> \param qs_env ...
3367 : !> \param bs_env ...
3368 : !> \param rirs_kernel (optional) .TRUE. if the Hartree + SEX kernels use the RI-RS grid backend
3369 : !> \author Maximilian Graml
3370 : !> \note Single source of truth shared by create_rtbse_env (sets the flag) and de_init_bs_env
3371 : !> (decides whether to retain nl_3c). KERNEL_RI=DEFAULT follows bs_env%do_gw_ri_rs;
3372 : !> RS/AO force; forced .FALSE. for non-linearized (full) RT-BSE (warn on explicit RS).
3373 : ! **************************************************************************************************
3374 288 : SUBROUTINE rtbse_resolve_rirs_flag(qs_env, bs_env, rirs_kernel)
3375 : TYPE(qs_environment_type), POINTER :: qs_env
3376 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3377 : LOGICAL, INTENT(OUT), OPTIONAL :: rirs_kernel
3378 :
3379 : INTEGER :: kernel_ri
3380 : LOGICAL :: my_rirs_kernel
3381 : TYPE(dft_control_type), POINTER :: dft_control
3382 : TYPE(section_vals_type), POINTER :: input
3383 :
3384 144 : NULLIFY (dft_control, input)
3385 144 : CALL get_qs_env(qs_env, dft_control=dft_control, input=input)
3386 :
3387 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%KERNEL_RI", &
3388 144 : i_val=kernel_ri)
3389 : SELECT CASE (kernel_ri)
3390 : CASE (rtp_bse_kernel_ri_rs)
3391 0 : my_rirs_kernel = .TRUE.
3392 : CASE (rtp_bse_kernel_ri_ao)
3393 0 : my_rirs_kernel = .FALSE.
3394 : CASE DEFAULT ! rtp_bse_kernel_ri_default
3395 144 : my_rirs_kernel = bs_env%do_gw_ri_rs
3396 : END SELECT
3397 :
3398 : ! RI-RS kernels are implemented for linearized RT-BSE only; full RT-BSE always uses AO-RI.
3399 144 : IF (dft_control%rtp_control%rtp_method /= rtp_method_bse_linearized) THEN
3400 28 : IF (kernel_ri == rtp_bse_kernel_ri_rs) THEN
3401 0 : CPWARN("RI-RS kernels are implemented for linearized RT-BSE only; forcing AO")
3402 : END IF
3403 : my_rirs_kernel = .FALSE.
3404 : END IF
3405 :
3406 144 : IF (PRESENT(rirs_kernel)) rirs_kernel = my_rirs_kernel
3407 :
3408 144 : END SUBROUTINE rtbse_resolve_rirs_flag
3409 :
3410 : ! **************************************************************************************************
3411 : !> \brief ...
3412 : !> \param bs_env ...
3413 : !> \param Sigma_c_ikp_n_freq Complex diagonal self-energy at the signed fit frequencies.
3414 : !> \param Sigma_x_ikp_n ...
3415 : !> \param V_xc_ikp_n ...
3416 : !> \param eigenval_scf ...
3417 : !> \param ikp ...
3418 : !> \param ispin ...
3419 : ! **************************************************************************************************
3420 482 : SUBROUTINE analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_xc_ikp_n, &
3421 482 : eigenval_scf, ikp, ispin)
3422 :
3423 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3424 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Sigma_c_ikp_n_freq
3425 : REAL(KIND=dp), DIMENSION(:) :: Sigma_x_ikp_n, V_xc_ikp_n, eigenval_scf
3426 : INTEGER :: ikp, ispin
3427 :
3428 : CHARACTER(LEN=*), PARAMETER :: routineN = 'analyt_conti_and_print'
3429 :
3430 : CHARACTER(len=3) :: occ_vir
3431 : CHARACTER(len=default_path_length) :: fname
3432 : CHARACTER(len=default_string_length) :: gw_label
3433 : INTEGER :: handle, i_mo, ikp_for_print, iunit, &
3434 : n_mo, nkp
3435 : LOGICAL :: is_bandstruc_kpoint, print_DOS_kpoints, &
3436 : print_ikp
3437 482 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dummy, eigenval_G, Sigma_c_ikp_n_qp
3438 :
3439 482 : CALL timeset(routineN, handle)
3440 :
3441 482 : n_mo = bs_env%n_ao
3442 2410 : ALLOCATE (dummy(n_mo), Sigma_c_ikp_n_qp(n_mo), eigenval_G(n_mo))
3443 482 : Sigma_c_ikp_n_qp(:) = 0.0_dp
3444 :
3445 : ! Eigenvalues of the Green's function that produced Σ^c: the DFT ones for G0W0, the
3446 : ! previous cycle's quasiparticle energies for evGW0. They set the Newton start value,
3447 : ! the Z/m linearization point and the Hedin shift; the QP equation itself stays
3448 : ! referenced to the DFT eigenvalues through the Eigenval_scf argument below.
3449 482 : IF (bs_env%gw_flavour == evGW0 .AND. ALLOCATED(bs_env%eigenval_evGW0)) THEN
3450 256 : eigenval_G(:) = bs_env%eigenval_evGW0(:, ikp, ispin)
3451 : ELSE
3452 5968 : eigenval_G(:) = eigenval_scf(:)
3453 : END IF
3454 :
3455 6224 : DO i_mo = 1, n_mo
3456 :
3457 : ! parallelization
3458 5742 : IF (MODULO(i_mo, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
3459 :
3460 : CALL continuation_pade(Sigma_c_ikp_n_qp, &
3461 : bs_env%imag_freq_points_fit, dummy, dummy, &
3462 : Sigma_c_ikp_n_freq, &
3463 : Sigma_x_ikp_n(:) - V_xc_ikp_n(:), &
3464 : eigenval_G(:), eigenval_scf(:), &
3465 : bs_env%do_hedin_shift, &
3466 : i_mo, bs_env%n_occ(ispin), bs_env%n_vir(ispin), &
3467 : bs_env%nparam_pade, bs_env%num_freq_points_fit, &
3468 : ri_rpa_g0w0_crossing_newton, bs_env%n_occ(ispin), &
3469 43459 : 0.0_dp, .TRUE., .FALSE., 1, e_fermi_ext=bs_env%e_fermi(ispin))
3470 : END DO
3471 :
3472 482 : CALL bs_env%para_env%sum(Sigma_c_ikp_n_qp)
3473 :
3474 482 : CALL correct_obvious_fitting_fails(Sigma_c_ikp_n_qp, ispin, bs_env)
3475 :
3476 : bs_env%eigenval_GW(:, ikp, ispin) = eigenval_scf(:) + &
3477 : Sigma_c_ikp_n_qp(:) + &
3478 : Sigma_x_ikp_n(:) - &
3479 6224 : V_xc_ikp_n(:)
3480 :
3481 482 : IF (ALLOCATED(bs_env%eigenval_G0W0) .AND. bs_env%evgw0_i_iter <= 1) THEN
3482 6028 : bs_env%eigenval_G0W0(:, ikp, ispin) = bs_env%eigenval_GW(:, ikp, ispin)
3483 : END IF
3484 :
3485 6224 : bs_env%eigenval_HF(:, ikp, ispin) = eigenval_scf(:) + Sigma_x_ikp_n(:) - V_xc_ikp_n(:)
3486 :
3487 : ! only print eigenvalues of DOS k-points in case no bandstructure path has been given
3488 482 : print_DOS_kpoints = (bs_env%nkp_only_bs <= 0)
3489 : ! in kpoints_DOS, the last nkp_only_bs are bandstructure k-points
3490 482 : is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
3491 482 : print_ikp = print_DOS_kpoints .OR. is_bandstruc_kpoint
3492 :
3493 482 : IF (bs_env%para_env%is_source() .AND. print_ikp) THEN
3494 :
3495 225 : IF (print_DOS_kpoints) THEN
3496 194 : nkp = bs_env%nkp_only_DOS
3497 194 : ikp_for_print = ikp
3498 : ELSE
3499 31 : nkp = bs_env%nkp_only_bs
3500 31 : ikp_for_print = ikp - bs_env%nkp_only_DOS
3501 : END IF
3502 :
3503 225 : fname = "bandstructure_SCF_and_G0W0"
3504 :
3505 : ! in an evGW0 run the spectrum of every cycle is appended, so only the very first
3506 : ! cycle replaces the file
3507 225 : IF (ikp_for_print == 1 .AND. ispin == 1 .AND. bs_env%evgw0_i_iter <= 1) THEN
3508 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", &
3509 65 : file_action="WRITE")
3510 : ELSE
3511 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="OLD", &
3512 160 : file_action="WRITE", file_position="APPEND")
3513 : END IF
3514 :
3515 225 : IF (bs_env%gw_flavour == evGW0 .AND. ikp_for_print == 1 .AND. ispin == 1) THEN
3516 9 : WRITE (iunit, "(A)") " "
3517 9 : WRITE (iunit, "(A,I0)") "evGW0 cycle: ", bs_env%evgw0_i_iter
3518 : END IF
3519 :
3520 225 : WRITE (iunit, "(A)") " "
3521 225 : WRITE (iunit, "(A10,I7,A25,3F10.4,T90,A7,I2)") "kpoint: ", ikp_for_print, "coordinate: ", &
3522 1125 : bs_env%kpoints_DOS%xkp(:, ikp), "spin: ", ispin
3523 225 : WRITE (iunit, "(A)") " "
3524 225 : gw_label = "ϵ_nk^"//TRIM(gw_flavour_label(bs_env))//" (eV)"
3525 225 : WRITE (iunit, "(A5,A12,3A17,A16,A18)") "n", "k", "ϵ_nk^DFT (eV)", "Σ^c_nk (eV)", &
3526 450 : "Σ^x_nk (eV)", "v_nk^xc (eV)", TRIM(gw_label)
3527 225 : WRITE (iunit, "(A)") " "
3528 :
3529 2952 : DO i_mo = 1, n_mo
3530 2727 : IF (i_mo <= bs_env%n_occ(ispin)) occ_vir = 'occ'
3531 2727 : IF (i_mo > bs_env%n_occ(ispin)) occ_vir = 'vir'
3532 2727 : WRITE (iunit, "(I5,3A,I5,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', ikp_for_print, &
3533 2727 : eigenval_scf(i_mo)*evolt, &
3534 2727 : Sigma_c_ikp_n_qp(i_mo)*evolt, &
3535 2727 : Sigma_x_ikp_n(i_mo)*evolt, &
3536 2727 : V_xc_ikp_n(i_mo)*evolt, &
3537 5679 : bs_env%eigenval_GW(i_mo, ikp, ispin)*evolt
3538 : END DO
3539 :
3540 225 : WRITE (iunit, "(A)") " "
3541 :
3542 225 : CALL close_file(iunit)
3543 :
3544 : END IF
3545 :
3546 482 : CALL timestop(handle)
3547 :
3548 964 : END SUBROUTINE analyt_conti_and_print
3549 :
3550 : ! **************************************************************************************************
3551 : !> \brief ...
3552 : !> \param Sigma_c_ikp_n_qp ...
3553 : !> \param ispin ...
3554 : !> \param bs_env ...
3555 : ! **************************************************************************************************
3556 482 : SUBROUTINE correct_obvious_fitting_fails(Sigma_c_ikp_n_qp, ispin, bs_env)
3557 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Sigma_c_ikp_n_qp
3558 : INTEGER :: ispin
3559 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3560 :
3561 : CHARACTER(LEN=*), PARAMETER :: routineN = 'correct_obvious_fitting_fails'
3562 :
3563 : INTEGER :: handle, homo, i_mo, j_mo, &
3564 : n_levels_scissor, n_mo
3565 : LOGICAL :: is_occ, is_vir
3566 : REAL(KIND=dp) :: sum_Sigma_c
3567 :
3568 482 : CALL timeset(routineN, handle)
3569 :
3570 482 : n_mo = bs_env%n_ao
3571 482 : homo = bs_env%n_occ(ispin)
3572 :
3573 6224 : DO i_mo = 1, n_mo
3574 :
3575 : ! if |𝚺^c| > 13 eV, we use a scissors shift
3576 6224 : IF (ABS(Sigma_c_ikp_n_qp(i_mo)) > 13.0_dp/evolt) THEN
3577 :
3578 0 : is_occ = (i_mo <= homo)
3579 0 : is_vir = (i_mo > homo)
3580 :
3581 0 : n_levels_scissor = 0
3582 0 : sum_Sigma_c = 0.0_dp
3583 :
3584 : ! compute scissor
3585 0 : DO j_mo = 1, n_mo
3586 :
3587 : ! only compute scissor from other GW levels close in energy
3588 0 : IF (is_occ .AND. j_mo > homo) CYCLE
3589 0 : IF (is_vir .AND. j_mo <= homo) CYCLE
3590 0 : IF (ABS(i_mo - j_mo) > 10) CYCLE
3591 0 : IF (i_mo == j_mo) CYCLE
3592 :
3593 0 : n_levels_scissor = n_levels_scissor + 1
3594 0 : sum_Sigma_c = sum_Sigma_c + Sigma_c_ikp_n_qp(j_mo)
3595 :
3596 : END DO
3597 :
3598 : ! overwrite the self-energy with scissor shift
3599 0 : Sigma_c_ikp_n_qp(i_mo) = sum_Sigma_c/REAL(n_levels_scissor, KIND=dp)
3600 :
3601 : END IF
3602 :
3603 : END DO ! i_mo
3604 :
3605 482 : CALL timestop(handle)
3606 :
3607 482 : END SUBROUTINE correct_obvious_fitting_fails
3608 :
3609 : ! **************************************************************************************************
3610 : !> \brief Build the real Gamma occupied or virtual one-particle propagator.
3611 : !> \param bs_env MO coefficients, eigenvalues and reusable work matrix
3612 : !> \param tau imaginary time
3613 : !> \param fm_G_Gamma output AO propagator
3614 : !> \param ispin spin channel
3615 : !> \param occ select the occupied sector
3616 : !> \param vir select the virtual sector
3617 : ! **************************************************************************************************
3618 9188 : SUBROUTINE G_occ_vir(bs_env, tau, fm_G_Gamma, ispin, occ, vir)
3619 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3620 : REAL(KIND=dp) :: tau
3621 : TYPE(cp_fm_type) :: fm_G_Gamma
3622 : INTEGER :: ispin
3623 : LOGICAL :: occ, vir
3624 :
3625 : CHARACTER(LEN=*), PARAMETER :: routineN = 'G_occ_vir'
3626 :
3627 : INTEGER :: first_mo, handle, homo, i_mo, last_mo
3628 : REAL(KIND=dp) :: tau_E
3629 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: column_factors
3630 :
3631 9188 : CALL timeset(routineN, handle)
3632 :
3633 9188 : CPASSERT(occ .NEQV. vir)
3634 :
3635 9188 : homo = bs_env%n_occ(ispin)
3636 9188 : first_mo = 1
3637 : ! Canonical orthogonalization removes the modes above n_mo_retained.
3638 9188 : last_mo = MIN(bs_env%n_ao, bs_env%n_mo_retained)
3639 9188 : IF (occ) THEN
3640 4660 : last_mo = MIN(last_mo, homo)
3641 : ELSE
3642 4528 : first_mo = homo + 1
3643 : END IF
3644 27564 : ALLOCATE (column_factors(first_mo:last_mo))
3645 59716 : DO i_mo = first_mo, last_mo
3646 50528 : tau_E = ABS(tau*0.5_dp*(bs_env%eigenval_scf_Gamma(i_mo, ispin) - bs_env%e_fermi(ispin)))
3647 59716 : IF (tau_E < bs_env%stabilize_exp) THEN
3648 50216 : column_factors(i_mo) = EXP(-tau_E)
3649 : ELSE
3650 312 : column_factors(i_mo) = 0.0_dp
3651 : END IF
3652 : END DO
3653 :
3654 : CALL build_ao_propagator(bs_env%fm_mo_coeff_Gamma(ispin), column_factors, first_mo, &
3655 9188 : bs_env%fm_work_mo(1), fm_G_Gamma)
3656 9188 : DEALLOCATE (column_factors)
3657 :
3658 9188 : CALL timestop(handle)
3659 :
3660 9188 : END SUBROUTINE G_occ_vir
3661 :
3662 : END MODULE gw_utils
|
