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