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