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