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