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