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