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