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