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