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