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