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