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: gto_basis_set_type
17 : USE bibliography, ONLY: Graml2024,&
18 : cite_reference
19 : USE cell_types, ONLY: cell_type,&
20 : pbc
21 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
22 : cp_blacs_env_release,&
23 : cp_blacs_env_type
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26 : dbcsr_allocate_matrix_set,&
27 : dbcsr_deallocate_matrix_set
28 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
29 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
30 : cp_fm_struct_release,&
31 : cp_fm_struct_type
32 : USE cp_fm_types, ONLY: cp_fm_create,&
33 : cp_fm_set_all
34 : USE cp_log_handling, ONLY: cp_get_default_logger,&
35 : cp_logger_type
36 : USE cp_output_handling, ONLY: cp_print_key_generate_filename
37 : USE dbcsr_api, ONLY: dbcsr_create,&
38 : dbcsr_p_type
39 : USE dbt_api, ONLY: &
40 : dbt_clear, dbt_create, dbt_destroy, dbt_filter, dbt_iterator_blocks_left, &
41 : dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, &
42 : dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
43 : USE input_constants, ONLY: do_potential_truncated,&
44 : xc_none
45 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
46 : section_vals_type,&
47 : section_vals_val_get,&
48 : section_vals_val_set
49 : USE kinds, ONLY: default_string_length,&
50 : dp,&
51 : int_8
52 : USE kpoint_methods, ONLY: kpoint_init_cell_index
53 : USE kpoint_types, ONLY: kpoint_create,&
54 : kpoint_type
55 : USE libint_wrapper, ONLY: cp_libint_static_cleanup,&
56 : cp_libint_static_init
57 : USE machine, ONLY: m_memory,&
58 : m_walltime
59 : USE mathlib, ONLY: gcd
60 : USE message_passing, ONLY: mp_cart_type,&
61 : mp_para_env_type
62 : USE minimax_exp, ONLY: get_exp_minimax_coeff
63 : USE minimax_exp_gw, ONLY: get_exp_minimax_coeff_gw
64 : USE minimax_rpa, ONLY: get_rpa_minimax_coeff,&
65 : get_rpa_minimax_coeff_larger_grid
66 : USE mp2_gpw, ONLY: create_mat_munu
67 : USE mp2_grids, ONLY: get_l_sq_wghts_cos_tf_t_to_w,&
68 : get_l_sq_wghts_cos_tf_w_to_t,&
69 : get_l_sq_wghts_sin_tf_t_to_w
70 : USE mp2_ri_2c, ONLY: setup_trunc_coulomb_pot_for_exchange_self_energy
71 : USE particle_methods, ONLY: get_particle_set
72 : USE particle_types, ONLY: particle_type
73 : USE physcon, ONLY: angstrom,&
74 : evolt
75 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
76 : USE qs_energy_types, ONLY: qs_energy_type
77 : USE qs_environment_types, ONLY: get_qs_env,&
78 : qs_env_part_release,&
79 : qs_environment_type
80 : USE qs_integral_utils, ONLY: basis_set_list_setup
81 : USE qs_interactions, ONLY: init_interaction_radii_orb_basis
82 : USE qs_kind_types, ONLY: get_qs_kind,&
83 : qs_kind_type
84 : USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
85 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
86 : USE qs_tensors, ONLY: build_3c_integrals,&
87 : build_3c_neighbor_lists,&
88 : get_tensor_occupancy,&
89 : neighbor_list_3c_destroy
90 : USE qs_tensors_types, ONLY: create_2c_tensor,&
91 : create_3c_tensor,&
92 : distribution_3d_create,&
93 : distribution_3d_type,&
94 : neighbor_list_3c_type
95 : #include "base/base_uses.f90"
96 :
97 : IMPLICIT NONE
98 :
99 : PRIVATE
100 :
101 : PUBLIC :: create_and_init_bs_env_for_gw, de_init_bs_env, get_i_j_atoms, &
102 : kpoint_init_cell_index_simple, compute_xkp
103 :
104 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils'
105 :
106 : CONTAINS
107 :
108 : ! **************************************************************************************************
109 : !> \brief ...
110 : !> \param qs_env ...
111 : !> \param bs_env ...
112 : !> \param bs_sec ...
113 : ! **************************************************************************************************
114 18 : SUBROUTINE create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec)
115 : TYPE(qs_environment_type), POINTER :: qs_env
116 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
117 : TYPE(section_vals_type), POINTER :: bs_sec
118 :
119 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_and_init_bs_env_for_gw'
120 :
121 : INTEGER :: handle
122 :
123 18 : CALL timeset(routineN, handle)
124 :
125 18 : CALL cite_reference(Graml2024)
126 :
127 18 : CALL read_gw_input_parameters(bs_env, bs_sec)
128 :
129 18 : CALL print_header_and_input_parameters(bs_env)
130 :
131 18 : CALL setup_AO_and_RI_basis_set(qs_env, bs_env)
132 :
133 18 : CALL get_RI_basis_and_basis_function_indices(qs_env, bs_env)
134 :
135 18 : CALL setup_kpoints_chi_eps_W(qs_env, bs_env, bs_env%kpoints_chi_eps_W)
136 :
137 18 : CALL set_heuristic_parameters(bs_env, qs_env)
138 :
139 18 : CALL set_parallelization_parameters(qs_env, bs_env)
140 :
141 18 : CALL allocate_and_fill_matrices_and_arrays(qs_env, bs_env)
142 :
143 18 : CALL cp_libint_static_init()
144 :
145 18 : CALL create_tensors(qs_env, bs_env)
146 :
147 18 : CALL set_sparsity_parallelization_parameters(bs_env)
148 :
149 18 : CALL check_for_restart_files(qs_env, bs_env)
150 :
151 18 : CALL compute_fm_V_xc_Gamma(qs_env, bs_env)
152 :
153 18 : CALL setup_time_and_frequency_minimax_grid(bs_env)
154 :
155 : ! free memory in qs_env; only if one is not calculating the LDOS because
156 : ! we need real-space grid operations in pw_env, task_list for the LDOS
157 : ! Recommendation in case of memory issues: first perform GW calculation without calculating
158 : ! LDOS (to safe memor). Then, use GW restart files
159 : ! in a subsequent calculation to calculate the LDOS
160 18 : IF (.NOT. bs_env%do_ldos) THEN
161 14 : CALL qs_env_part_release(qs_env)
162 : END IF
163 :
164 18 : CALL timestop(handle)
165 :
166 18 : END SUBROUTINE create_and_init_bs_env_for_gw
167 :
168 : ! **************************************************************************************************
169 : !> \brief ...
170 : !> \param bs_env ...
171 : ! **************************************************************************************************
172 18 : SUBROUTINE de_init_bs_env(bs_env)
173 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
174 :
175 : CHARACTER(LEN=*), PARAMETER :: routineN = 'de_init_bs_env'
176 :
177 : INTEGER :: handle
178 :
179 18 : CALL timeset(routineN, handle)
180 : ! deallocate quantities here which:
181 : ! 1. cannot be deallocated in bs_env_release due to circular dependencies
182 : ! 2. consume a lot of memory and should not be kept until the quantity is
183 : ! deallocated in bs_env_release
184 :
185 18 : IF (ASSOCIATED(bs_env%nl_3c%ij_list)) CALL neighbor_list_3c_destroy(bs_env%nl_3c)
186 :
187 18 : CALL cp_libint_static_cleanup()
188 :
189 18 : CALL timestop(handle)
190 :
191 18 : END SUBROUTINE de_init_bs_env
192 :
193 : ! **************************************************************************************************
194 : !> \brief ...
195 : !> \param bs_env ...
196 : !> \param bs_sec ...
197 : ! **************************************************************************************************
198 18 : SUBROUTINE read_gw_input_parameters(bs_env, bs_sec)
199 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
200 : TYPE(section_vals_type), POINTER :: bs_sec
201 :
202 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_gw_input_parameters'
203 :
204 : INTEGER :: handle
205 : TYPE(section_vals_type), POINTER :: gw_sec
206 :
207 18 : CALL timeset(routineN, handle)
208 :
209 18 : NULLIFY (gw_sec)
210 18 : gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
211 :
212 18 : CALL section_vals_val_get(gw_sec, "NUM_TIME_FREQ_POINTS", i_val=bs_env%num_time_freq_points)
213 18 : CALL section_vals_val_get(gw_sec, "EPS_FILTER", r_val=bs_env%eps_filter)
214 18 : CALL section_vals_val_get(gw_sec, "MEMORY_PER_PROC", r_val=bs_env%input_memory_per_proc_GB)
215 18 : CALL section_vals_val_get(gw_sec, "APPROX_KP_EXTRAPOL", l_val=bs_env%approx_kp_extrapol)
216 :
217 18 : CALL timestop(handle)
218 :
219 18 : END SUBROUTINE read_gw_input_parameters
220 :
221 : ! **************************************************************************************************
222 : !> \brief ...
223 : !> \param qs_env ...
224 : !> \param bs_env ...
225 : ! **************************************************************************************************
226 18 : SUBROUTINE setup_AO_and_RI_basis_set(qs_env, bs_env)
227 : TYPE(qs_environment_type), POINTER :: qs_env
228 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
229 :
230 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_AO_and_RI_basis_set'
231 :
232 : INTEGER :: handle, natom, nkind
233 18 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
234 18 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
235 :
236 18 : CALL timeset(routineN, handle)
237 :
238 : CALL get_qs_env(qs_env, &
239 : qs_kind_set=qs_kind_set, &
240 : particle_set=particle_set, &
241 18 : natom=natom, nkind=nkind)
242 :
243 : ! set up basis
244 90 : ALLOCATE (bs_env%sizes_RI(natom), bs_env%sizes_AO(natom))
245 162 : ALLOCATE (bs_env%basis_set_RI(nkind), bs_env%basis_set_AO(nkind))
246 :
247 18 : CALL basis_set_list_setup(bs_env%basis_set_RI, "RI_AUX", qs_kind_set)
248 18 : CALL basis_set_list_setup(bs_env%basis_set_AO, "ORB", qs_kind_set)
249 :
250 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_RI, &
251 18 : basis=bs_env%basis_set_RI)
252 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_AO, &
253 18 : basis=bs_env%basis_set_AO)
254 :
255 18 : CALL timestop(handle)
256 :
257 18 : END SUBROUTINE setup_AO_and_RI_basis_set
258 :
259 : ! **************************************************************************************************
260 : !> \brief ...
261 : !> \param qs_env ...
262 : !> \param bs_env ...
263 : ! **************************************************************************************************
264 18 : SUBROUTINE get_RI_basis_and_basis_function_indices(qs_env, bs_env)
265 : TYPE(qs_environment_type), POINTER :: qs_env
266 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
267 :
268 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_RI_basis_and_basis_function_indices'
269 :
270 : INTEGER :: handle, i_RI, iatom, ikind, iset, &
271 : max_AO_bf_per_atom, n_ao_test, n_atom, &
272 : n_kind, n_RI, nset, nsgf, u
273 18 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
274 18 : INTEGER, DIMENSION(:), POINTER :: l_max, l_min, nsgf_set
275 18 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
276 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
277 18 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
278 :
279 18 : CALL timeset(routineN, handle)
280 :
281 : ! determine RI basis set size
282 18 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
283 :
284 18 : n_kind = SIZE(qs_kind_set)
285 18 : n_atom = bs_env%n_atom
286 :
287 18 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
288 :
289 54 : DO ikind = 1, n_kind
290 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
291 36 : basis_type="RI_AUX")
292 54 : CPASSERT(ASSOCIATED(basis_set_a))
293 : END DO
294 :
295 54 : ALLOCATE (bs_env%i_RI_start_from_atom(n_atom))
296 36 : ALLOCATE (bs_env%i_RI_end_from_atom(n_atom))
297 36 : ALLOCATE (bs_env%i_ao_start_from_atom(n_atom))
298 36 : ALLOCATE (bs_env%i_ao_end_from_atom(n_atom))
299 :
300 18 : n_RI = 0
301 92 : DO iatom = 1, n_atom
302 74 : bs_env%i_RI_start_from_atom(iatom) = n_RI + 1
303 74 : ikind = kind_of(iatom)
304 74 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
305 74 : n_RI = n_RI + nsgf
306 92 : bs_env%i_RI_end_from_atom(iatom) = n_RI
307 : END DO
308 18 : bs_env%n_RI = n_RI
309 :
310 18 : max_AO_bf_per_atom = 0
311 18 : n_ao_test = 0
312 92 : DO iatom = 1, n_atom
313 74 : bs_env%i_ao_start_from_atom(iatom) = n_ao_test + 1
314 74 : ikind = kind_of(iatom)
315 74 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="ORB")
316 74 : n_ao_test = n_ao_test + nsgf
317 74 : bs_env%i_ao_end_from_atom(iatom) = n_ao_test
318 92 : max_AO_bf_per_atom = MAX(max_AO_bf_per_atom, nsgf)
319 : END DO
320 18 : CPASSERT(n_ao_test == bs_env%n_ao)
321 18 : bs_env%max_AO_bf_per_atom = max_AO_bf_per_atom
322 :
323 54 : ALLOCATE (bs_env%l_RI(n_RI))
324 18 : i_RI = 0
325 92 : DO iatom = 1, n_atom
326 74 : ikind = kind_of(iatom)
327 :
328 74 : nset = bs_env%basis_set_RI(ikind)%gto_basis_set%nset
329 74 : l_max => bs_env%basis_set_RI(ikind)%gto_basis_set%lmax
330 74 : l_min => bs_env%basis_set_RI(ikind)%gto_basis_set%lmin
331 74 : nsgf_set => bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf_set
332 :
333 202 : DO iset = 1, nset
334 110 : CPASSERT(l_max(iset) == l_min(iset))
335 296 : bs_env%l_RI(i_RI + 1:i_RI + nsgf_set(iset)) = l_max(iset)
336 184 : i_RI = i_RI + nsgf_set(iset)
337 : END DO
338 :
339 : END DO
340 18 : CPASSERT(i_RI == n_RI)
341 :
342 18 : u = bs_env%unit_nr
343 :
344 18 : IF (u > 0) THEN
345 9 : WRITE (UNIT=u, FMT="(T2,A)") " "
346 9 : WRITE (UNIT=u, FMT="(T2,2A,T75,I8)") "Number of auxiliary Gaussian basis functions ", &
347 18 : "for χ, ε, W", n_RI
348 : END IF
349 :
350 18 : CALL timestop(handle)
351 :
352 36 : END SUBROUTINE get_RI_basis_and_basis_function_indices
353 :
354 : ! **************************************************************************************************
355 : !> \brief ...
356 : !> \param qs_env ...
357 : !> \param bs_env ...
358 : !> \param kpoints ...
359 : ! **************************************************************************************************
360 18 : SUBROUTINE setup_kpoints_chi_eps_W(qs_env, bs_env, kpoints)
361 :
362 : TYPE(qs_environment_type), POINTER :: qs_env
363 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
364 : TYPE(kpoint_type), POINTER :: kpoints
365 :
366 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_kpoints_chi_eps_W'
367 :
368 : INTEGER :: handle, i_dim, n_dim, nkp, nkp_extra, &
369 : nkp_orig, u
370 : INTEGER, DIMENSION(3) :: nkp_grid, nkp_grid_extra, periodic
371 : REAL(KIND=dp) :: exp_s_p, n_dim_inv
372 :
373 18 : CALL timeset(routineN, handle)
374 :
375 : ! routine adapted from mp2_integrals.F
376 18 : NULLIFY (kpoints)
377 18 : CALL kpoint_create(kpoints)
378 :
379 18 : kpoints%kp_scheme = "GENERAL"
380 :
381 72 : periodic(1:3) = bs_env%periodic(1:3)
382 :
383 72 : DO i_dim = 1, 3
384 :
385 54 : CPASSERT(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
386 :
387 54 : IF (periodic(i_dim) == 1) nkp_grid(i_dim) = 4
388 54 : IF (periodic(i_dim) == 0) nkp_grid(i_dim) = 1
389 :
390 54 : IF (periodic(i_dim) == 1) THEN
391 : ! only even k-meshes in periodic direction implemented
392 36 : CPASSERT(MODULO(nkp_grid(i_dim), 2) == 0)
393 : END IF
394 54 : IF (periodic(i_dim) == 0) THEN
395 : ! single k-kpoint in non-periodic direction needed
396 18 : CPASSERT(nkp_grid(i_dim) == 1)
397 : END IF
398 :
399 54 : IF (periodic(i_dim) == 1) nkp_grid_extra(i_dim) = nkp_grid(i_dim) + 2
400 72 : IF (periodic(i_dim) == 0) nkp_grid_extra(i_dim) = 1
401 :
402 : END DO
403 :
404 18 : nkp_orig = MAX(nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2, 1)
405 :
406 18 : nkp_extra = nkp_grid_extra(1)*nkp_grid_extra(2)*nkp_grid_extra(3)/2
407 :
408 18 : nkp = nkp_orig + nkp_extra
409 :
410 72 : kpoints%nkp_grid(1:3) = nkp_grid(1:3)
411 18 : kpoints%nkp = nkp
412 :
413 72 : bs_env%nkp_grid_chi_eps_W_orig(1:3) = nkp_grid(1:3)
414 72 : bs_env%nkp_grid_chi_eps_W_extra(1:3) = nkp_grid_extra(1:3)
415 18 : bs_env%nkp_chi_eps_W_orig = nkp_orig
416 18 : bs_env%nkp_chi_eps_W_extra = nkp_extra
417 18 : bs_env%nkp_chi_eps_W_orig_plus_extra = nkp
418 :
419 90 : ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
420 90 : ALLOCATE (bs_env%wkp_no_extra(nkp), bs_env%wkp_s_p(nkp))
421 :
422 18 : CALL compute_xkp(kpoints%xkp, 1, nkp_orig, nkp_grid)
423 18 : CALL compute_xkp(kpoints%xkp, nkp_orig + 1, nkp, nkp_grid_extra)
424 :
425 72 : n_dim = SUM(periodic)
426 18 : IF (n_dim == 0) THEN
427 : ! molecules
428 0 : kpoints%wkp(1) = 1.0_dp
429 0 : bs_env%wkp_s_p(1) = 1.0_dp
430 0 : bs_env%wkp_no_extra(1) = 1.0_dp
431 : ELSE
432 :
433 18 : n_dim_inv = 1.0_dp/REAL(n_dim, KIND=dp)
434 :
435 : ! k-point weights are chosen to automatically extrapolate the k-point mesh
436 18 : CALL compute_wkp(kpoints%wkp(1:nkp_orig), nkp_orig, nkp_extra, n_dim_inv)
437 18 : CALL compute_wkp(kpoints%wkp(nkp_orig + 1:nkp), nkp_extra, nkp_orig, n_dim_inv)
438 :
439 162 : bs_env%wkp_no_extra(1:nkp_orig) = 0.0_dp
440 342 : bs_env%wkp_no_extra(nkp_orig + 1:nkp) = 1.0_dp/REAL(nkp_extra, KIND=dp)
441 :
442 18 : IF (n_dim == 3) THEN
443 : ! W_PQ(k) for an s-function P and a p-function Q diverges as 1/k at k=0
444 : ! (instead of 1/k^2 for P and Q both being s-functions).
445 0 : exp_s_p = 2.0_dp*n_dim_inv
446 0 : CALL compute_wkp(bs_env%wkp_s_p(1:nkp_orig), nkp_orig, nkp_extra, exp_s_p)
447 0 : CALL compute_wkp(bs_env%wkp_s_p(nkp_orig + 1:nkp), nkp_extra, nkp_orig, exp_s_p)
448 : ELSE
449 486 : bs_env%wkp_s_p(1:nkp) = bs_env%wkp_no_extra(1:nkp)
450 : END IF
451 :
452 : END IF
453 :
454 18 : IF (bs_env%approx_kp_extrapol) THEN
455 2 : bs_env%wkp_orig = 1.0_dp/REAL(nkp_orig, KIND=dp)
456 : END IF
457 :
458 18 : CALL kpoint_init_cell_index_simple(kpoints, qs_env)
459 :
460 : ! heuristic parameter: how many k-points for χ, ε, and W are used simultaneously
461 : ! (less simultaneous k-points: less memory, but more computational effort because of
462 : ! recomputation of V(k))
463 18 : bs_env%nkp_chi_eps_W_batch = 4
464 :
465 : bs_env%num_chi_eps_W_batches = (bs_env%nkp_chi_eps_W_orig_plus_extra - 1)/ &
466 18 : bs_env%nkp_chi_eps_W_batch + 1
467 :
468 18 : u = bs_env%unit_nr
469 :
470 18 : IF (u > 0) THEN
471 9 : WRITE (UNIT=u, FMT="(T2,A)") " "
472 9 : WRITE (UNIT=u, FMT="(T2,1A,T71,3I4)") "K-point mesh 1 for χ, ε, W", nkp_grid(1:3)
473 9 : WRITE (UNIT=u, FMT="(T2,2A,T71,3I4)") "K-point mesh 2 for χ, ε, W ", &
474 18 : "(for k-point extrapolation of W)", nkp_grid_extra(1:3)
475 9 : WRITE (UNIT=u, FMT="(T2,A,T80,L)") "Approximate the k-point extrapolation", &
476 18 : bs_env%approx_kp_extrapol
477 : END IF
478 :
479 18 : CALL timestop(handle)
480 :
481 18 : END SUBROUTINE setup_kpoints_chi_eps_W
482 :
483 : ! **************************************************************************************************
484 : !> \brief ...
485 : !> \param kpoints ...
486 : !> \param qs_env ...
487 : ! **************************************************************************************************
488 36 : SUBROUTINE kpoint_init_cell_index_simple(kpoints, qs_env)
489 :
490 : TYPE(kpoint_type), POINTER :: kpoints
491 : TYPE(qs_environment_type), POINTER :: qs_env
492 :
493 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_init_cell_index_simple'
494 :
495 : INTEGER :: handle
496 : TYPE(dft_control_type), POINTER :: dft_control
497 : TYPE(mp_para_env_type), POINTER :: para_env
498 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
499 36 : POINTER :: sab_orb
500 :
501 36 : CALL timeset(routineN, handle)
502 :
503 36 : NULLIFY (dft_control, para_env, sab_orb)
504 36 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, dft_control=dft_control, sab_orb=sab_orb)
505 36 : CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, dft_control)
506 :
507 36 : CALL timestop(handle)
508 :
509 36 : END SUBROUTINE kpoint_init_cell_index_simple
510 :
511 : ! **************************************************************************************************
512 : !> \brief ...
513 : !> \param xkp ...
514 : !> \param ikp_start ...
515 : !> \param ikp_end ...
516 : !> \param grid ...
517 : ! **************************************************************************************************
518 54 : SUBROUTINE compute_xkp(xkp, ikp_start, ikp_end, grid)
519 :
520 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
521 : INTEGER :: ikp_start, ikp_end
522 : INTEGER, DIMENSION(3) :: grid
523 :
524 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_xkp'
525 :
526 : INTEGER :: handle, i, ix, iy, iz
527 :
528 54 : CALL timeset(routineN, handle)
529 :
530 54 : i = ikp_start
531 216 : DO ix = 1, grid(1)
532 956 : DO iy = 1, grid(2)
533 1902 : DO iz = 1, grid(3)
534 :
535 1000 : IF (i > ikp_end) CYCLE
536 :
537 500 : xkp(1, i) = REAL(2*ix - grid(1) - 1, KIND=dp)/(2._dp*REAL(grid(1), KIND=dp))
538 500 : xkp(2, i) = REAL(2*iy - grid(2) - 1, KIND=dp)/(2._dp*REAL(grid(2), KIND=dp))
539 500 : xkp(3, i) = REAL(2*iz - grid(3) - 1, KIND=dp)/(2._dp*REAL(grid(3), KIND=dp))
540 1740 : i = i + 1
541 :
542 : END DO
543 : END DO
544 : END DO
545 :
546 54 : CALL timestop(handle)
547 :
548 54 : END SUBROUTINE compute_xkp
549 :
550 : ! **************************************************************************************************
551 : !> \brief ...
552 : !> \param wkp ...
553 : !> \param nkp_1 ...
554 : !> \param nkp_2 ...
555 : !> \param exponent ...
556 : ! **************************************************************************************************
557 36 : SUBROUTINE compute_wkp(wkp, nkp_1, nkp_2, exponent)
558 : REAL(KIND=dp), DIMENSION(:) :: wkp
559 : INTEGER :: nkp_1, nkp_2
560 : REAL(KIND=dp) :: exponent
561 :
562 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_wkp'
563 :
564 : INTEGER :: handle
565 : REAL(KIND=dp) :: nkp_ratio
566 :
567 36 : CALL timeset(routineN, handle)
568 :
569 36 : nkp_ratio = REAL(nkp_2, KIND=dp)/REAL(nkp_1, KIND=dp)
570 :
571 504 : wkp(:) = 1.0_dp/REAL(nkp_1, KIND=dp)/(1.0_dp - nkp_ratio**exponent)
572 :
573 36 : CALL timestop(handle)
574 :
575 36 : END SUBROUTINE
576 :
577 : ! **************************************************************************************************
578 : !> \brief ...
579 : !> \param qs_env ...
580 : !> \param bs_env ...
581 : ! **************************************************************************************************
582 18 : SUBROUTINE allocate_and_fill_matrices_and_arrays(qs_env, bs_env)
583 : TYPE(qs_environment_type), POINTER :: qs_env
584 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
585 :
586 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_and_fill_matrices_and_arrays'
587 :
588 : INTEGER :: handle, i_t, num_time_freq_points
589 : TYPE(cp_blacs_env_type), POINTER :: blacs_env, blacs_env_tensor
590 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_RI_global
591 18 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
592 : TYPE(mp_para_env_type), POINTER :: para_env
593 :
594 18 : CALL timeset(routineN, handle)
595 :
596 18 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env, matrix_s=matrix_s)
597 :
598 18 : fm_struct => bs_env%fm_ks_Gamma(1)%matrix_struct
599 :
600 18 : CALL cp_fm_create(bs_env%fm_Gocc, fm_struct)
601 18 : CALL cp_fm_create(bs_env%fm_Gvir, fm_struct)
602 18 : CALL cp_fm_create(bs_env%fm_h_G0W0_Gamma, fm_struct)
603 :
604 18 : NULLIFY (fm_struct_RI_global)
605 : CALL cp_fm_struct_create(fm_struct_RI_global, context=blacs_env, nrow_global=bs_env%n_RI, &
606 18 : ncol_global=bs_env%n_RI, para_env=para_env)
607 18 : CALL cp_fm_create(bs_env%fm_RI_RI, fm_struct_RI_global)
608 18 : CALL cp_fm_create(bs_env%fm_chi_Gamma_freq, fm_struct_RI_global)
609 18 : CALL cp_fm_create(bs_env%fm_W_MIC_freq, fm_struct_RI_global)
610 18 : IF (bs_env%approx_kp_extrapol) THEN
611 2 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_extra, fm_struct_RI_global)
612 2 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_no_extra, fm_struct_RI_global)
613 2 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_extra, 0.0_dp)
614 2 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_no_extra, 0.0_dp)
615 : END IF
616 18 : CALL cp_fm_struct_release(fm_struct_RI_global)
617 :
618 90 : ALLOCATE (bs_env%eigenval_G0W0(bs_env%n_ao, bs_env%nkp_DOS, bs_env%n_spin))
619 :
620 18 : num_time_freq_points = bs_env%num_time_freq_points
621 :
622 54 : ALLOCATE (bs_env%imag_freq_points(num_time_freq_points))
623 54 : ALLOCATE (bs_env%imag_time_points(num_time_freq_points))
624 72 : ALLOCATE (bs_env%weights_cos_t_to_w(num_time_freq_points, num_time_freq_points))
625 72 : ALLOCATE (bs_env%weights_cos_w_to_t(num_time_freq_points, num_time_freq_points))
626 72 : ALLOCATE (bs_env%weights_sin_t_to_w(num_time_freq_points, num_time_freq_points))
627 :
628 : ! create blacs_env for subgroups of tensor operations
629 18 : NULLIFY (blacs_env_tensor)
630 18 : CALL cp_blacs_env_create(blacs_env=blacs_env_tensor, para_env=bs_env%para_env_tensor)
631 :
632 : ! allocate dbcsr matrices in the tensor subgroup; actually, one only needs a small
633 : ! subset of blocks in the tensor subgroup, however, all atomic blocks are allocated.
634 : ! One might think of creating a dbcsr matrix with only the blocks that are needed
635 : ! in the tensor subgroup
636 : CALL create_mat_munu(bs_env%mat_ao_ao_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
637 18 : blacs_env_tensor, do_ri_aux_basis=.FALSE.)
638 :
639 : CALL create_mat_munu(bs_env%mat_RI_RI_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
640 18 : blacs_env_tensor, do_ri_aux_basis=.TRUE.)
641 :
642 : CALL create_mat_munu(bs_env%mat_RI_RI, qs_env, bs_env%eps_atom_grid_2d_mat, &
643 18 : blacs_env, do_ri_aux_basis=.TRUE.)
644 :
645 18 : CALL cp_blacs_env_release(blacs_env_tensor)
646 :
647 18 : NULLIFY (bs_env%mat_chi_Gamma_tau)
648 18 : CALL dbcsr_allocate_matrix_set(bs_env%mat_chi_Gamma_tau, bs_env%num_time_freq_points)
649 :
650 190 : DO i_t = 1, bs_env%num_time_freq_points
651 172 : ALLOCATE (bs_env%mat_chi_Gamma_tau(i_t)%matrix)
652 190 : CALL dbcsr_create(bs_env%mat_chi_Gamma_tau(i_t)%matrix, template=bs_env%mat_RI_RI%matrix)
653 : END DO
654 :
655 18 : CALL timestop(handle)
656 :
657 18 : END SUBROUTINE allocate_and_fill_matrices_and_arrays
658 :
659 : ! **************************************************************************************************
660 : !> \brief ...
661 : !> \param qs_env ...
662 : !> \param bs_env ...
663 : ! **************************************************************************************************
664 18 : SUBROUTINE create_tensors(qs_env, bs_env)
665 : TYPE(qs_environment_type), POINTER :: qs_env
666 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
667 :
668 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_tensors'
669 :
670 : INTEGER :: handle, n_atom_step, RI_atom
671 : INTEGER(int_8) :: mem, non_zero_elements_sum, nze
672 : REAL(dp) :: max_dist_AO_atoms, occ, occupation_sum
673 126 : TYPE(dbt_type) :: t_3c_global
674 18 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_global_array
675 : TYPE(neighbor_list_3c_type) :: nl_3c_global
676 :
677 18 : CALL timeset(routineN, handle)
678 :
679 : ! be careful: routine needs bs_env%eps_3c_int which is set in set_heuristic_parameters
680 18 : CALL init_interaction_radii(bs_env)
681 :
682 : ! split blocks does not improve load balancing/efficienfy for tensor contraction, so we go
683 : ! with the standard atomic blocks
684 : CALL create_3c_t(bs_env%t_RI_AO__AO, bs_env%para_env_tensor, "(RI AO | AO)", [1, 2], [3], &
685 : bs_env%sizes_RI, bs_env%sizes_AO, &
686 18 : create_nl_3c=.TRUE., nl_3c=bs_env%nl_3c, qs_env=qs_env)
687 : CALL create_3c_t(bs_env%t_RI__AO_AO, bs_env%para_env_tensor, "(RI | AO AO)", [1], [2, 3], &
688 18 : bs_env%sizes_RI, bs_env%sizes_AO)
689 :
690 18 : CALL create_2c_t(bs_env, bs_env%sizes_RI, bs_env%sizes_AO)
691 :
692 : ! check the sparsity of 3c integral tensor (µν|P); calculate maximum distance between
693 : ! AO atoms µ, ν where at least a single integral (µν|P) is larger than the filter threshold
694 : CALL create_3c_t(t_3c_global, bs_env%para_env, "(RI AO | AO)", [1, 2], [3], &
695 : bs_env%sizes_RI, bs_env%sizes_AO, &
696 18 : create_nl_3c=.TRUE., nl_3c=nl_3c_global, qs_env=qs_env)
697 :
698 18 : CALL m_memory(mem)
699 18 : CALL bs_env%para_env%max(mem)
700 :
701 162 : ALLOCATE (t_3c_global_array(1, 1))
702 18 : CALL dbt_create(t_3c_global, t_3c_global_array(1, 1))
703 :
704 18 : CALL bs_env%para_env%sync()
705 18 : bs_env%t1 = m_walltime()
706 :
707 18 : occupation_sum = 0.0_dp
708 18 : non_zero_elements_sum = 0
709 18 : max_dist_AO_atoms = 0.0_dp
710 18 : n_atom_step = INT(SQRT(REAL(bs_env%n_atom, KIND=dp)))
711 : ! do not compute full 3c integrals at once because it may cause out of memory
712 62 : DO RI_atom = 1, bs_env%n_atom, n_atom_step
713 :
714 : CALL build_3c_integrals(t_3c_global_array, &
715 : bs_env%eps_filter, &
716 : qs_env, &
717 : nl_3c_global, &
718 : int_eps=bs_env%eps_3c_int, &
719 : basis_i=bs_env%basis_set_RI, &
720 : basis_j=bs_env%basis_set_AO, &
721 : basis_k=bs_env%basis_set_AO, &
722 : bounds_i=[RI_atom, MIN(RI_atom + n_atom_step - 1, bs_env%n_atom)], &
723 : potential_parameter=bs_env%ri_metric, &
724 132 : desymmetrize=.FALSE.)
725 :
726 44 : CALL dbt_filter(t_3c_global_array(1, 1), bs_env%eps_filter)
727 :
728 44 : CALL bs_env%para_env%sync()
729 :
730 44 : CALL get_tensor_occupancy(t_3c_global_array(1, 1), nze, occ)
731 44 : non_zero_elements_sum = non_zero_elements_sum + nze
732 44 : occupation_sum = occupation_sum + occ
733 :
734 44 : CALL get_max_dist_AO_atoms(t_3c_global_array(1, 1), max_dist_AO_atoms, qs_env)
735 :
736 106 : CALL dbt_clear(t_3c_global_array(1, 1))
737 :
738 : END DO
739 :
740 18 : bs_env%t2 = m_walltime()
741 :
742 18 : bs_env%occupation_3c_int = occupation_sum
743 18 : bs_env%max_dist_AO_atoms = max_dist_AO_atoms
744 :
745 18 : CALL dbt_destroy(t_3c_global)
746 18 : CALL dbt_destroy(t_3c_global_array(1, 1))
747 36 : DEALLOCATE (t_3c_global_array)
748 :
749 18 : CALL neighbor_list_3c_destroy(nl_3c_global)
750 :
751 18 : IF (bs_env%unit_nr > 0) THEN
752 9 : WRITE (bs_env%unit_nr, '(T2,A)') ''
753 : WRITE (bs_env%unit_nr, '(T2,A,F27.1,A)') &
754 9 : 'Computed 3-center integrals (µν|P), execution time', bs_env%t2 - bs_env%t1, ' s'
755 9 : WRITE (bs_env%unit_nr, '(T2,A,F48.3,A)') 'Percentage of non-zero (µν|P)', &
756 18 : occupation_sum*100, ' %'
757 9 : WRITE (bs_env%unit_nr, '(T2,A,F33.1,A)') 'Max. distance between µ,ν in non-zero (µν|P)', &
758 18 : max_dist_AO_atoms*angstrom, ' A'
759 9 : WRITE (bs_env%unit_nr, '(T2,2A,I20,A)') 'Required memory if storing all 3-center ', &
760 18 : 'integrals (µν|P)', INT(REAL(non_zero_elements_sum, KIND=dp)*8.0E-9_dp), ' GB'
761 : END IF
762 :
763 18 : CALL timestop(handle)
764 :
765 72 : END SUBROUTINE create_tensors
766 :
767 : ! **************************************************************************************************
768 : !> \brief ...
769 : !> \param bs_env ...
770 : !> \param sizes_RI ...
771 : !> \param sizes_AO ...
772 : ! **************************************************************************************************
773 18 : SUBROUTINE create_2c_t(bs_env, sizes_RI, sizes_AO)
774 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
775 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_RI, sizes_AO
776 :
777 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_2c_t'
778 :
779 : INTEGER :: handle
780 18 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2
781 : INTEGER, DIMENSION(2) :: pdims_2d
782 54 : TYPE(dbt_pgrid_type) :: pgrid_2d
783 :
784 18 : CALL timeset(routineN, handle)
785 :
786 : ! inspired from rpa_im_time.F / hfx_types.F
787 :
788 18 : pdims_2d = 0
789 18 : CALL dbt_pgrid_create(bs_env%para_env_tensor, pdims_2d, pgrid_2d)
790 :
791 : CALL create_2c_tensor(bs_env%t_G, dist_1, dist_2, pgrid_2d, sizes_AO, sizes_AO, &
792 18 : name="(AO | AO)")
793 18 : DEALLOCATE (dist_1, dist_2)
794 : CALL create_2c_tensor(bs_env%t_chi, dist_1, dist_2, pgrid_2d, sizes_RI, sizes_RI, &
795 18 : name="(RI | RI)")
796 18 : DEALLOCATE (dist_1, dist_2)
797 : CALL create_2c_tensor(bs_env%t_W, dist_1, dist_2, pgrid_2d, sizes_RI, sizes_RI, &
798 18 : name="(RI | RI)")
799 18 : DEALLOCATE (dist_1, dist_2)
800 18 : CALL dbt_pgrid_destroy(pgrid_2d)
801 :
802 18 : CALL timestop(handle)
803 :
804 18 : END SUBROUTINE create_2c_t
805 :
806 : ! **************************************************************************************************
807 : !> \brief ...
808 : !> \param tensor ...
809 : !> \param para_env ...
810 : !> \param tensor_name ...
811 : !> \param map1 ...
812 : !> \param map2 ...
813 : !> \param sizes_RI ...
814 : !> \param sizes_AO ...
815 : !> \param create_nl_3c ...
816 : !> \param nl_3c ...
817 : !> \param qs_env ...
818 : ! **************************************************************************************************
819 54 : SUBROUTINE create_3c_t(tensor, para_env, tensor_name, map1, map2, sizes_RI, sizes_AO, &
820 : create_nl_3c, nl_3c, qs_env)
821 : TYPE(dbt_type) :: tensor
822 : TYPE(mp_para_env_type), POINTER :: para_env
823 : CHARACTER(LEN=12) :: tensor_name
824 : INTEGER, DIMENSION(:) :: map1, map2
825 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_RI, sizes_AO
826 : LOGICAL, OPTIONAL :: create_nl_3c
827 : TYPE(neighbor_list_3c_type), OPTIONAL :: nl_3c
828 : TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
829 :
830 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_3c_t'
831 :
832 : INTEGER :: handle, nkind
833 54 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_AO_1, dist_AO_2, dist_RI
834 : INTEGER, DIMENSION(3) :: pcoord, pdims, pdims_3d
835 : LOGICAL :: my_create_nl_3c
836 162 : TYPE(dbt_pgrid_type) :: pgrid_3d
837 : TYPE(distribution_3d_type) :: dist_3d
838 54 : TYPE(mp_cart_type) :: mp_comm_t3c_2
839 54 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
840 :
841 54 : CALL timeset(routineN, handle)
842 :
843 54 : pdims_3d = 0
844 54 : CALL dbt_pgrid_create(para_env, pdims_3d, pgrid_3d)
845 : CALL create_3c_tensor(tensor, dist_RI, dist_AO_1, dist_AO_2, &
846 : pgrid_3d, sizes_RI, sizes_AO, sizes_AO, &
847 54 : map1=map1, map2=map2, name=tensor_name)
848 :
849 54 : IF (PRESENT(create_nl_3c)) THEN
850 36 : my_create_nl_3c = create_nl_3c
851 : ELSE
852 : my_create_nl_3c = .FALSE.
853 : END IF
854 :
855 36 : IF (my_create_nl_3c) THEN
856 36 : CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set)
857 36 : CALL dbt_mp_environ_pgrid(pgrid_3d, pdims, pcoord)
858 36 : CALL mp_comm_t3c_2%create(pgrid_3d%mp_comm_2d, 3, pdims)
859 : CALL distribution_3d_create(dist_3d, dist_RI, dist_AO_1, dist_AO_2, &
860 36 : nkind, particle_set, mp_comm_t3c_2, own_comm=.TRUE.)
861 :
862 : CALL build_3c_neighbor_lists(nl_3c, &
863 : qs_env%bs_env%basis_set_RI, &
864 : qs_env%bs_env%basis_set_AO, &
865 : qs_env%bs_env%basis_set_AO, &
866 : dist_3d, qs_env%bs_env%ri_metric, &
867 36 : "GW_3c_nl", qs_env, own_dist=.TRUE.)
868 : END IF
869 :
870 54 : DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
871 54 : CALL dbt_pgrid_destroy(pgrid_3d)
872 :
873 54 : CALL timestop(handle)
874 :
875 108 : END SUBROUTINE create_3c_t
876 :
877 : ! **************************************************************************************************
878 : !> \brief ...
879 : !> \param bs_env ...
880 : ! **************************************************************************************************
881 18 : SUBROUTINE init_interaction_radii(bs_env)
882 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
883 :
884 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_interaction_radii'
885 :
886 : INTEGER :: handle, ibasis
887 : TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
888 :
889 18 : CALL timeset(routineN, handle)
890 :
891 54 : DO ibasis = 1, SIZE(bs_env%basis_set_AO)
892 :
893 36 : orb_basis => bs_env%basis_set_AO(ibasis)%gto_basis_set
894 36 : CALL init_interaction_radii_orb_basis(orb_basis, bs_env%eps_3c_int)
895 :
896 36 : ri_basis => bs_env%basis_set_RI(ibasis)%gto_basis_set
897 54 : CALL init_interaction_radii_orb_basis(ri_basis, bs_env%eps_3c_int)
898 :
899 : END DO
900 :
901 18 : CALL timestop(handle)
902 :
903 18 : END SUBROUTINE init_interaction_radii
904 :
905 : ! **************************************************************************************************
906 : !> \brief ...
907 : !> \param t_3c_int ...
908 : !> \param max_dist_AO_atoms ...
909 : !> \param qs_env ...
910 : ! **************************************************************************************************
911 44 : SUBROUTINE get_max_dist_AO_atoms(t_3c_int, max_dist_AO_atoms, qs_env)
912 : TYPE(dbt_type) :: t_3c_int
913 : REAL(KIND=dp) :: max_dist_AO_atoms
914 : TYPE(qs_environment_type), POINTER :: qs_env
915 :
916 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_max_dist_AO_atoms'
917 :
918 : INTEGER :: atom_1, atom_2, handle, num_cells
919 : INTEGER, DIMENSION(3) :: atom_ind
920 44 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
921 : REAL(KIND=dp) :: abs_rab
922 : REAL(KIND=dp), DIMENSION(3) :: rab
923 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
924 : TYPE(cell_type), POINTER :: cell
925 : TYPE(dbt_iterator_type) :: iter
926 : TYPE(mp_para_env_type), POINTER :: para_env
927 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
928 :
929 44 : CALL timeset(routineN, handle)
930 :
931 44 : NULLIFY (cell, particle_set, para_env)
932 44 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, para_env=para_env)
933 :
934 : !$OMP PARALLEL DEFAULT(NONE) &
935 : !$OMP SHARED(t_3c_int, max_dist_AO_atoms, num_cells, index_to_cell, hmat, particle_set, cell) &
936 44 : !$OMP PRIVATE(iter,atom_ind,rab, abs_rab, atom_1, atom_2)
937 : CALL dbt_iterator_start(iter, t_3c_int)
938 : DO WHILE (dbt_iterator_blocks_left(iter))
939 : CALL dbt_iterator_next_block(iter, atom_ind)
940 :
941 : atom_1 = atom_ind(2)
942 : atom_2 = atom_ind(3)
943 :
944 : rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
945 :
946 : abs_rab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
947 :
948 : max_dist_AO_atoms = MAX(max_dist_AO_atoms, abs_rab)
949 :
950 : END DO
951 : CALL dbt_iterator_stop(iter)
952 : !$OMP END PARALLEL
953 :
954 44 : CALL para_env%max(max_dist_AO_atoms)
955 :
956 44 : CALL timestop(handle)
957 :
958 44 : END SUBROUTINE get_max_dist_AO_atoms
959 :
960 : ! **************************************************************************************************
961 : !> \brief ...
962 : !> \param bs_env ...
963 : ! **************************************************************************************************
964 18 : SUBROUTINE set_sparsity_parallelization_parameters(bs_env)
965 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
966 :
967 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_sparsity_parallelization_parameters'
968 :
969 : INTEGER :: handle, i_ivl, IL_ivl, j_ivl, n_atom_per_IL_ivl, n_atom_per_ivl, n_intervals_i, &
970 : n_intervals_inner_loop_atoms, n_intervals_j, u
971 :
972 18 : CALL timeset(routineN, handle)
973 :
974 : ! heuristic parameter to prevent out of memory
975 18 : bs_env%safety_factor_memory = 0.10_dp
976 :
977 : ! choose atomic range for λ ("i_atom"), ν ("j_atom") in
978 : ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
979 : ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
980 : ! such that M and N fit into the memory
981 : n_atom_per_ivl = INT(SQRT(bs_env%safety_factor_memory*bs_env%input_memory_per_proc &
982 : *bs_env%group_size_tensor/24/bs_env%n_RI &
983 18 : /SQRT(bs_env%occupation_3c_int)))/bs_env%max_AO_bf_per_atom
984 :
985 18 : n_intervals_i = (bs_env%n_atom_i - 1)/n_atom_per_ivl + 1
986 18 : n_intervals_j = (bs_env%n_atom_j - 1)/n_atom_per_ivl + 1
987 :
988 18 : bs_env%n_atom_per_interval_ij = n_atom_per_ivl
989 18 : bs_env%n_intervals_i = n_intervals_i
990 18 : bs_env%n_intervals_j = n_intervals_j
991 :
992 54 : ALLOCATE (bs_env%i_atom_intervals(2, n_intervals_i))
993 54 : ALLOCATE (bs_env%j_atom_intervals(2, n_intervals_j))
994 :
995 36 : DO i_ivl = 1, n_intervals_i
996 18 : bs_env%i_atom_intervals(1, i_ivl) = (i_ivl - 1)*n_atom_per_ivl + bs_env%atoms_i(1)
997 : bs_env%i_atom_intervals(2, i_ivl) = MIN(i_ivl*n_atom_per_ivl + bs_env%atoms_i(1) - 1, &
998 36 : bs_env%atoms_i(2))
999 : END DO
1000 :
1001 36 : DO j_ivl = 1, n_intervals_j
1002 18 : bs_env%j_atom_intervals(1, j_ivl) = (j_ivl - 1)*n_atom_per_ivl + bs_env%atoms_j(1)
1003 : bs_env%j_atom_intervals(2, j_ivl) = MIN(j_ivl*n_atom_per_ivl + bs_env%atoms_j(1) - 1, &
1004 36 : bs_env%atoms_j(2))
1005 : END DO
1006 :
1007 72 : ALLOCATE (bs_env%skip_Sigma_occ(n_intervals_i, n_intervals_j))
1008 54 : ALLOCATE (bs_env%skip_Sigma_vir(n_intervals_i, n_intervals_j))
1009 54 : bs_env%skip_Sigma_occ(:, :) = .FALSE.
1010 54 : bs_env%skip_Sigma_vir(:, :) = .FALSE.
1011 :
1012 : ! choose atomic range for µ and σ ("inner loop (IL) atom") in
1013 : ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
1014 : ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
1015 : n_atom_per_IL_ivl = MIN(INT(bs_env%safety_factor_memory*bs_env%input_memory_per_proc &
1016 : *bs_env%group_size_tensor/n_atom_per_ivl &
1017 : /bs_env%max_AO_bf_per_atom &
1018 : /bs_env%n_RI/8/SQRT(bs_env%occupation_3c_int) &
1019 18 : /bs_env%max_AO_bf_per_atom), bs_env%n_atom)
1020 :
1021 18 : n_intervals_inner_loop_atoms = (bs_env%n_atom - 1)/n_atom_per_IL_ivl + 1
1022 :
1023 18 : bs_env%n_atom_per_IL_interval = n_atom_per_IL_ivl
1024 18 : bs_env%n_intervals_inner_loop_atoms = n_intervals_inner_loop_atoms
1025 :
1026 54 : ALLOCATE (bs_env%inner_loop_atom_intervals(2, n_intervals_inner_loop_atoms))
1027 36 : DO IL_ivl = 1, n_intervals_inner_loop_atoms
1028 18 : bs_env%inner_loop_atom_intervals(1, IL_ivl) = (IL_ivl - 1)*n_atom_per_IL_ivl + 1
1029 36 : bs_env%inner_loop_atom_intervals(2, IL_ivl) = MIN(IL_ivl*n_atom_per_IL_ivl, bs_env%n_atom)
1030 : END DO
1031 :
1032 18 : u = bs_env%unit_nr
1033 18 : IF (u > 0) THEN
1034 9 : WRITE (u, '(T2,A)') ''
1035 9 : WRITE (u, '(T2,A,I33)') 'Number of i and j atoms in M_λνP(τ), N_νλQ(τ):', n_atom_per_ivl
1036 9 : WRITE (u, '(T2,A,I18)') 'Number of inner loop atoms for µ in M_λνP = sum_µ (µν|P) G_µλ', &
1037 18 : n_atom_per_IL_ivl
1038 : END IF
1039 :
1040 18 : CALL timestop(handle)
1041 :
1042 18 : END SUBROUTINE set_sparsity_parallelization_parameters
1043 :
1044 : ! **************************************************************************************************
1045 : !> \brief ...
1046 : !> \param qs_env ...
1047 : !> \param bs_env ...
1048 : ! **************************************************************************************************
1049 18 : SUBROUTINE check_for_restart_files(qs_env, bs_env)
1050 : TYPE(qs_environment_type), POINTER :: qs_env
1051 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1052 :
1053 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_for_restart_files'
1054 :
1055 : CHARACTER(LEN=9) :: frmt
1056 : CHARACTER(LEN=default_string_length) :: f_chi, f_S_n, f_S_p, f_S_x, f_W_t, &
1057 : prefix, project_name
1058 : INTEGER :: handle, i_spin, i_t_or_w, ind, n_spin, &
1059 : num_time_freq_points
1060 : LOGICAL :: chi_exists, Sigma_neg_time_exists, &
1061 : Sigma_pos_time_exists, &
1062 : Sigma_x_spin_exists, W_time_exists
1063 : TYPE(cp_logger_type), POINTER :: logger
1064 : TYPE(section_vals_type), POINTER :: input, print_key
1065 :
1066 18 : CALL timeset(routineN, handle)
1067 :
1068 18 : num_time_freq_points = bs_env%num_time_freq_points
1069 18 : n_spin = bs_env%n_spin
1070 :
1071 54 : ALLOCATE (bs_env%read_chi(num_time_freq_points))
1072 36 : ALLOCATE (bs_env%calc_chi(num_time_freq_points))
1073 72 : ALLOCATE (bs_env%Sigma_c_exists(num_time_freq_points, n_spin))
1074 :
1075 18 : CALL get_qs_env(qs_env, input=input)
1076 :
1077 18 : logger => cp_get_default_logger()
1078 18 : print_key => section_vals_get_subs_vals(input, 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART')
1079 : project_name = cp_print_key_generate_filename(logger, print_key, extension="", &
1080 18 : my_local=.FALSE.)
1081 18 : WRITE (prefix, '(2A)') TRIM(project_name), "-RESTART_"
1082 18 : bs_env%prefix = prefix
1083 :
1084 18 : bs_env%all_W_exist = .TRUE.
1085 :
1086 190 : DO i_t_or_w = 1, num_time_freq_points
1087 :
1088 172 : IF (i_t_or_w < 10) THEN
1089 156 : WRITE (frmt, '(A)') '(3A,I1,A)'
1090 156 : WRITE (f_chi, frmt) TRIM(prefix), bs_env%chi_name, "_0", i_t_or_w, ".matrix"
1091 156 : WRITE (f_W_t, frmt) TRIM(prefix), bs_env%W_time_name, "_0", i_t_or_w, ".matrix"
1092 16 : ELSE IF (i_t_or_w < 100) THEN
1093 16 : WRITE (frmt, '(A)') '(3A,I2,A)'
1094 16 : WRITE (f_chi, frmt) TRIM(prefix), bs_env%chi_name, "_", i_t_or_w, ".matrix"
1095 16 : WRITE (f_W_t, frmt) TRIM(prefix), bs_env%W_time_name, "_", i_t_or_w, ".matrix"
1096 : ELSE
1097 0 : CPABORT('Please implement more than 99 time/frequency points.')
1098 : END IF
1099 :
1100 172 : INQUIRE (file=TRIM(f_chi), exist=chi_exists)
1101 172 : INQUIRE (file=TRIM(f_W_t), exist=W_time_exists)
1102 :
1103 172 : bs_env%read_chi(i_t_or_w) = chi_exists
1104 172 : bs_env%calc_chi(i_t_or_w) = .NOT. chi_exists
1105 :
1106 172 : bs_env%all_W_exist = bs_env%all_W_exist .AND. W_time_exists
1107 :
1108 : ! the self-energy is spin-dependent
1109 422 : DO i_spin = 1, n_spin
1110 :
1111 232 : ind = i_t_or_w + (i_spin - 1)*num_time_freq_points
1112 :
1113 232 : IF (ind < 10) THEN
1114 156 : WRITE (frmt, '(A)') '(3A,I1,A)'
1115 156 : WRITE (f_S_p, frmt) TRIM(prefix), bs_env%Sigma_p_name, "_0", ind, ".matrix"
1116 156 : WRITE (f_S_n, frmt) TRIM(prefix), bs_env%Sigma_n_name, "_0", ind, ".matrix"
1117 76 : ELSE IF (i_t_or_w < 100) THEN
1118 76 : WRITE (frmt, '(A)') '(3A,I2,A)'
1119 76 : WRITE (f_S_p, frmt) TRIM(prefix), bs_env%Sigma_p_name, "_", ind, ".matrix"
1120 76 : WRITE (f_S_n, frmt) TRIM(prefix), bs_env%Sigma_n_name, "_", ind, ".matrix"
1121 : END IF
1122 :
1123 232 : INQUIRE (file=TRIM(f_S_p), exist=Sigma_pos_time_exists)
1124 232 : INQUIRE (file=TRIM(f_S_n), exist=Sigma_neg_time_exists)
1125 :
1126 : bs_env%Sigma_c_exists(i_t_or_w, i_spin) = Sigma_pos_time_exists .AND. &
1127 556 : Sigma_neg_time_exists
1128 :
1129 : END DO
1130 :
1131 : END DO
1132 :
1133 18 : IF (bs_env%all_W_exist) THEN
1134 44 : bs_env%read_chi(:) = .FALSE.
1135 44 : bs_env%calc_chi(:) = .FALSE.
1136 : END IF
1137 :
1138 18 : bs_env%Sigma_x_exists = .TRUE.
1139 42 : DO i_spin = 1, n_spin
1140 24 : WRITE (f_S_x, '(3A,I1,A)') TRIM(prefix), bs_env%Sigma_x_name, "_0", i_spin, ".matrix"
1141 24 : INQUIRE (file=TRIM(f_S_x), exist=Sigma_x_spin_exists)
1142 58 : bs_env%Sigma_x_exists = bs_env%Sigma_x_exists .AND. Sigma_x_spin_exists
1143 : END DO
1144 :
1145 18 : CALL timestop(handle)
1146 :
1147 18 : END SUBROUTINE check_for_restart_files
1148 :
1149 : ! **************************************************************************************************
1150 : !> \brief ...
1151 : !> \param qs_env ...
1152 : !> \param bs_env ...
1153 : ! **************************************************************************************************
1154 18 : SUBROUTINE set_parallelization_parameters(qs_env, bs_env)
1155 : TYPE(qs_environment_type), POINTER :: qs_env
1156 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1157 :
1158 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_parallelization_parameters'
1159 :
1160 : INTEGER :: color_sub, dummy_1, dummy_2, handle, &
1161 : num_pe, num_t_groups, u
1162 : INTEGER(KIND=int_8) :: mem
1163 : TYPE(mp_para_env_type), POINTER :: para_env
1164 :
1165 18 : CALL timeset(routineN, handle)
1166 :
1167 18 : CALL get_qs_env(qs_env, para_env=para_env)
1168 :
1169 18 : num_pe = para_env%num_pe
1170 : ! use all processors for the group (in principle, number could be changed, but performance
1171 : ! seems to be best for a single group with all MPI processes per group)
1172 18 : bs_env%group_size_tensor = num_pe
1173 :
1174 : ! group_size_tensor must divide num_pe without rest; otherwise everything will be complicated
1175 18 : IF (MODULO(num_pe, bs_env%group_size_tensor) .NE. 0) THEN
1176 0 : CALL find_good_group_size(num_pe, bs_env%group_size_tensor)
1177 : END IF
1178 :
1179 : ! para_env_tensor for tensor subgroups
1180 18 : color_sub = para_env%mepos/bs_env%group_size_tensor
1181 18 : bs_env%tensor_group_color = color_sub
1182 :
1183 18 : ALLOCATE (bs_env%para_env_tensor)
1184 18 : CALL bs_env%para_env_tensor%from_split(para_env, color_sub)
1185 :
1186 18 : num_t_groups = para_env%num_pe/bs_env%group_size_tensor
1187 18 : bs_env%num_tensor_groups = num_t_groups
1188 :
1189 : CALL get_i_j_atoms(bs_env%atoms_i, bs_env%atoms_j, bs_env%n_atom_i, bs_env%n_atom_j, &
1190 18 : color_sub, bs_env)
1191 :
1192 54 : ALLOCATE (bs_env%atoms_i_t_group(2, num_t_groups))
1193 54 : ALLOCATE (bs_env%atoms_j_t_group(2, num_t_groups))
1194 36 : DO color_sub = 0, num_t_groups - 1
1195 : CALL get_i_j_atoms(bs_env%atoms_i_t_group(1:2, color_sub + 1), &
1196 : bs_env%atoms_j_t_group(1:2, color_sub + 1), &
1197 36 : dummy_1, dummy_2, color_sub, bs_env)
1198 : END DO
1199 :
1200 18 : CALL m_memory(mem)
1201 18 : CALL bs_env%para_env%max(mem)
1202 :
1203 18 : bs_env%input_memory_per_proc = INT(bs_env%input_memory_per_proc_GB*1.0E9_dp, KIND=int_8)
1204 :
1205 18 : u = bs_env%unit_nr
1206 18 : IF (u > 0) THEN
1207 9 : WRITE (u, '(T2,A)') ''
1208 9 : WRITE (u, '(T2,A,I47)') 'Group size for tensor operations', bs_env%group_size_tensor
1209 9 : WRITE (u, '(T2,A)') ''
1210 9 : WRITE (u, '(T2,A,F37.1,A)') 'Input: Available memory per MPI process', &
1211 18 : bs_env%input_memory_per_proc_GB, ' GB'
1212 9 : WRITE (u, '(T2,A,F35.1,A)') 'Used memory per MPI process before GW run', &
1213 18 : REAL(mem, KIND=dp)/1.E9_dp, ' GB'
1214 : END IF
1215 :
1216 18 : CALL timestop(handle)
1217 :
1218 36 : END SUBROUTINE set_parallelization_parameters
1219 :
1220 : ! **************************************************************************************************
1221 : !> \brief ...
1222 : !> \param num_pe ...
1223 : !> \param group_size ...
1224 : ! **************************************************************************************************
1225 0 : SUBROUTINE find_good_group_size(num_pe, group_size)
1226 :
1227 : INTEGER :: num_pe, group_size
1228 :
1229 : CHARACTER(LEN=*), PARAMETER :: routineN = 'find_good_group_size'
1230 :
1231 : INTEGER :: group_size_minus, group_size_orig, &
1232 : group_size_plus, handle, i_diff
1233 :
1234 0 : CALL timeset(routineN, handle)
1235 :
1236 0 : group_size_orig = group_size
1237 :
1238 0 : DO i_diff = 1, num_pe
1239 :
1240 0 : group_size_minus = group_size - i_diff
1241 :
1242 0 : IF (MODULO(num_pe, group_size_minus) == 0 .AND. group_size_minus > 0) THEN
1243 0 : group_size = group_size_minus
1244 0 : EXIT
1245 : END IF
1246 :
1247 0 : group_size_plus = group_size + i_diff
1248 :
1249 0 : IF (MODULO(num_pe, group_size_plus) == 0 .AND. group_size_plus <= num_pe) THEN
1250 0 : group_size = group_size_plus
1251 0 : EXIT
1252 : END IF
1253 :
1254 : END DO
1255 :
1256 0 : IF (group_size_orig == group_size) CPABORT("Group size error")
1257 :
1258 0 : CALL timestop(handle)
1259 :
1260 0 : END SUBROUTINE find_good_group_size
1261 :
1262 : ! **************************************************************************************************
1263 : !> \brief ...
1264 : !> \param atoms_i ...
1265 : !> \param atoms_j ...
1266 : !> \param n_atom_i ...
1267 : !> \param n_atom_j ...
1268 : !> \param color_sub ...
1269 : !> \param bs_env ...
1270 : ! **************************************************************************************************
1271 36 : SUBROUTINE get_i_j_atoms(atoms_i, atoms_j, n_atom_i, n_atom_j, color_sub, bs_env)
1272 :
1273 : INTEGER, DIMENSION(2) :: atoms_i, atoms_j
1274 : INTEGER :: n_atom_i, n_atom_j, color_sub
1275 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1276 :
1277 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_i_j_atoms'
1278 :
1279 : INTEGER :: handle, i_atoms_per_group, i_group, &
1280 : ipcol, ipcol_loop, iprow, iprow_loop, &
1281 : j_atoms_per_group, npcol, nprow
1282 :
1283 36 : CALL timeset(routineN, handle)
1284 :
1285 : ! create a square mesh of tensor groups for iatom and jatom; code from blacs_env_create
1286 36 : CALL square_mesh(nprow, npcol, bs_env%num_tensor_groups)
1287 :
1288 36 : i_group = 0
1289 72 : DO ipcol_loop = 0, npcol - 1
1290 108 : DO iprow_loop = 0, nprow - 1
1291 36 : IF (i_group == color_sub) THEN
1292 36 : iprow = iprow_loop
1293 36 : ipcol = ipcol_loop
1294 : END IF
1295 72 : i_group = i_group + 1
1296 : END DO
1297 : END DO
1298 :
1299 36 : IF (MODULO(bs_env%n_atom, nprow) == 0) THEN
1300 36 : i_atoms_per_group = bs_env%n_atom/nprow
1301 : ELSE
1302 0 : i_atoms_per_group = bs_env%n_atom/nprow + 1
1303 : END IF
1304 :
1305 36 : IF (MODULO(bs_env%n_atom, npcol) == 0) THEN
1306 36 : j_atoms_per_group = bs_env%n_atom/npcol
1307 : ELSE
1308 0 : j_atoms_per_group = bs_env%n_atom/npcol + 1
1309 : END IF
1310 :
1311 36 : atoms_i(1) = iprow*i_atoms_per_group + 1
1312 36 : atoms_i(2) = MIN((iprow + 1)*i_atoms_per_group, bs_env%n_atom)
1313 36 : n_atom_i = atoms_i(2) - atoms_i(1) + 1
1314 :
1315 36 : atoms_j(1) = ipcol*j_atoms_per_group + 1
1316 36 : atoms_j(2) = MIN((ipcol + 1)*j_atoms_per_group, bs_env%n_atom)
1317 36 : n_atom_j = atoms_j(2) - atoms_j(1) + 1
1318 :
1319 36 : CALL timestop(handle)
1320 :
1321 36 : END SUBROUTINE get_i_j_atoms
1322 :
1323 : ! **************************************************************************************************
1324 : !> \brief ...
1325 : !> \param nprow ...
1326 : !> \param npcol ...
1327 : !> \param nproc ...
1328 : ! **************************************************************************************************
1329 36 : SUBROUTINE square_mesh(nprow, npcol, nproc)
1330 : INTEGER :: nprow, npcol, nproc
1331 :
1332 : CHARACTER(LEN=*), PARAMETER :: routineN = 'square_mesh'
1333 :
1334 : INTEGER :: gcd_max, handle, ipe, jpe
1335 :
1336 36 : CALL timeset(routineN, handle)
1337 :
1338 36 : gcd_max = -1
1339 72 : DO ipe = 1, CEILING(SQRT(REAL(nproc, dp)))
1340 36 : jpe = nproc/ipe
1341 36 : IF (ipe*jpe .NE. nproc) CYCLE
1342 72 : IF (gcd(ipe, jpe) >= gcd_max) THEN
1343 36 : nprow = ipe
1344 36 : npcol = jpe
1345 36 : gcd_max = gcd(ipe, jpe)
1346 : END IF
1347 : END DO
1348 :
1349 36 : CALL timestop(handle)
1350 :
1351 36 : END SUBROUTINE square_mesh
1352 :
1353 : ! **************************************************************************************************
1354 : !> \brief ...
1355 : !> \param bs_env ...
1356 : !> \param qs_env ...
1357 : ! **************************************************************************************************
1358 18 : SUBROUTINE set_heuristic_parameters(bs_env, qs_env)
1359 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1360 : TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
1361 :
1362 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_heuristic_parameters'
1363 :
1364 : INTEGER :: handle
1365 :
1366 18 : CALL timeset(routineN, handle)
1367 :
1368 : ! use the same threshold for computing 3-center integrals (µν|P) as for filtering
1369 : ! tensor operations
1370 18 : bs_env%eps_3c_int = bs_env%eps_filter
1371 :
1372 : ! Determines number of cells used for summing the cells R in the Coulomb matrix,
1373 : ! V_PQ(k) = \sum_R <P,cell=0 | 1/r | Q,cell=R>. SIZE_LATTICE_SUM_V 3 gives
1374 : ! good convergence
1375 18 : bs_env%size_lattice_sum_V = 3
1376 :
1377 : ! for generating numerically stable minimax Fourier integration weights
1378 18 : bs_env%num_points_per_magnitude = 200
1379 :
1380 : ! for periodic systems and for 20 minimax points, we use a regularized minimax mesh
1381 : ! (from experience: regularized minimax meshes converges faster for periodic systems
1382 : ! and for 20 pts)
1383 72 : IF (SUM(bs_env%periodic) .NE. 0 .OR. bs_env%num_time_freq_points == 20) THEN
1384 18 : bs_env%regularization_minimax = 1.0E-6_dp
1385 : ELSE
1386 0 : bs_env%regularization_minimax = 0.0_dp
1387 : END IF
1388 :
1389 18 : bs_env%stabilize_exp = 70.0_dp
1390 18 : bs_env%eps_atom_grid_2d_mat = 1.0E-50_dp
1391 :
1392 : ! only use interval ω in [0, 27.211 eV] (1 Hartree = 27.211 eV) for virt, and ω in
1393 : ! [-27.211 eV, 0] for occ for use in analytic continuation of
1394 : ! self-energy Σ^c_n(iω,k) -> Σ^c_n(ϵ,k)
1395 18 : bs_env%freq_max_fit = 1.0_dp
1396 :
1397 : ! use a 16-parameter Padé fit
1398 18 : bs_env%nparam_pade = 16
1399 :
1400 : ! minimum block size for tensor operations, taken from MP2/RPA input
1401 18 : bs_env%min_block_size = 5
1402 :
1403 : ! resolution of the identity with the truncated Coulomb metric, cutoff radius 3 Angström
1404 18 : bs_env%ri_metric%potential_type = do_potential_truncated
1405 18 : bs_env%ri_metric%omega = 0.0_dp
1406 : ! cutoff radius = 3 Angström
1407 18 : bs_env%ri_metric%cutoff_radius = 3.0_dp/angstrom
1408 18 : bs_env%ri_metric%filename = "t_c_g.dat"
1409 :
1410 18 : bs_env%eps_eigval_mat_RI = 0.0_dp
1411 :
1412 : ! for periodic systems, we use the regularized resolution of the identity
1413 72 : IF (SUM(bs_env%periodic) == 0) THEN
1414 0 : bs_env%regularization_RI = 0.0_dp
1415 : ELSE
1416 18 : bs_env%regularization_RI = 1.0E-2_dp
1417 : END IF
1418 :
1419 : ! truncated Coulomb operator for exchange self-energy
1420 : ! (see details in Guidon, VandeVondele, Hutter, JCTC 5, 3010 (2009) and references therein)
1421 : CALL setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env, bs_env%trunc_coulomb, &
1422 18 : rel_cutoff_trunc_coulomb_ri_x=0.5_dp)
1423 :
1424 18 : CALL timestop(handle)
1425 :
1426 18 : END SUBROUTINE set_heuristic_parameters
1427 :
1428 : ! **************************************************************************************************
1429 : !> \brief ...
1430 : !> \param bs_env ...
1431 : ! **************************************************************************************************
1432 18 : SUBROUTINE print_header_and_input_parameters(bs_env)
1433 :
1434 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1435 :
1436 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_header_and_input_parameters'
1437 :
1438 : INTEGER :: handle, u
1439 :
1440 18 : CALL timeset(routineN, handle)
1441 :
1442 18 : u = bs_env%unit_nr
1443 :
1444 18 : IF (u > 0) THEN
1445 9 : WRITE (u, *) ' '
1446 9 : WRITE (u, '(T2,2A)') '------------------------------------------------', &
1447 18 : '-------------------------------'
1448 9 : WRITE (u, '(T2,2A)') '- ', &
1449 18 : ' -'
1450 9 : WRITE (u, '(T2,2A)') '- GW CALCULATION ', &
1451 18 : ' -'
1452 9 : WRITE (u, '(T2,2A)') '- ', &
1453 18 : ' -'
1454 9 : WRITE (u, '(T2,2A)') '------------------------------------------------', &
1455 18 : '-------------------------------'
1456 9 : WRITE (u, '(T2,A)') ' '
1457 9 : WRITE (u, '(T2,A,I45)') 'Input: Number of time/freq. points', bs_env%num_time_freq_points
1458 9 : WRITE (u, '(T2,A,ES27.1)') 'Input: Filter threshold for sparse tensor operations', &
1459 18 : bs_env%eps_filter
1460 : END IF
1461 :
1462 18 : CALL timestop(handle)
1463 :
1464 18 : END SUBROUTINE print_header_and_input_parameters
1465 :
1466 : ! **************************************************************************************************
1467 : !> \brief ...
1468 : !> \param qs_env ...
1469 : !> \param bs_env ...
1470 : ! **************************************************************************************************
1471 36 : SUBROUTINE compute_fm_V_xc_Gamma(qs_env, bs_env)
1472 : TYPE(qs_environment_type), POINTER :: qs_env
1473 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1474 :
1475 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_V_xc_Gamma'
1476 :
1477 : INTEGER :: handle, ispin, myfun, nimages
1478 : REAL(KIND=dp) :: energy_ex, energy_exc, energy_total
1479 18 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_ks_without_v_xc
1480 : TYPE(dft_control_type), POINTER :: dft_control
1481 : TYPE(qs_energy_type), POINTER :: energy
1482 : TYPE(section_vals_type), POINTER :: input, xc_section
1483 :
1484 18 : CALL timeset(routineN, handle)
1485 :
1486 18 : CALL get_qs_env(qs_env, input=input, energy=energy, dft_control=dft_control)
1487 :
1488 : ! previously, dft_control%nimages set to # neighbor cells, revert for Γ-only KS matrix
1489 18 : nimages = dft_control%nimages
1490 18 : dft_control%nimages = 1
1491 :
1492 : ! we need to reset XC functional, therefore, get XC input
1493 18 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1494 18 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
1495 18 : CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=xc_none)
1496 :
1497 : ! save the energy before the energy gets updated
1498 18 : energy_total = energy%total
1499 18 : energy_exc = energy%exc
1500 18 : energy_ex = energy%ex
1501 :
1502 18 : NULLIFY (mat_ks_without_v_xc)
1503 18 : CALL dbcsr_allocate_matrix_set(mat_ks_without_v_xc, bs_env%n_spin)
1504 :
1505 42 : DO ispin = 1, bs_env%n_spin
1506 24 : ALLOCATE (mat_ks_without_v_xc(ispin)%matrix)
1507 42 : CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1508 : END DO
1509 :
1510 : ! calculate KS-matrix without XC
1511 : CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
1512 18 : ext_ks_matrix=mat_ks_without_v_xc)
1513 :
1514 42 : DO ispin = 1, bs_env%n_spin
1515 : ! transfer dbcsr matrix to fm
1516 24 : CALL cp_fm_create(bs_env%fm_V_xc_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1517 24 : CALL copy_dbcsr_to_fm(mat_ks_without_v_xc(ispin)%matrix, bs_env%fm_V_xc_Gamma(ispin))
1518 :
1519 : ! finally compute the xc potential in the AO basis
1520 : CALL cp_fm_scale_and_add(alpha=-1.0_dp, matrix_a=bs_env%fm_V_xc_Gamma(ispin), &
1521 42 : beta=1.0_dp, matrix_b=bs_env%fm_ks_Gamma(ispin))
1522 : END DO
1523 :
1524 : ! set back the energy
1525 18 : energy%total = energy_total
1526 18 : energy%exc = energy_exc
1527 18 : energy%ex = energy_ex
1528 :
1529 : ! set back nimages
1530 18 : dft_control%nimages = nimages
1531 :
1532 : ! set the DFT functional and HF fraction back
1533 : CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
1534 18 : i_val=myfun)
1535 :
1536 18 : CALL dbcsr_deallocate_matrix_set(mat_ks_without_v_xc)
1537 :
1538 18 : CALL timestop(handle)
1539 :
1540 18 : END SUBROUTINE compute_fm_V_xc_Gamma
1541 :
1542 : ! **************************************************************************************************
1543 : !> \brief ...
1544 : !> \param bs_env ...
1545 : ! **************************************************************************************************
1546 18 : SUBROUTINE setup_time_and_frequency_minimax_grid(bs_env)
1547 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1548 :
1549 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_time_and_frequency_minimax_grid'
1550 :
1551 : INTEGER :: handle, homo, i_w, ierr, j_w, n_mo, &
1552 : num_time_freq_points, u
1553 : REAL(KIND=dp) :: E_max, E_min, E_range, max_error_min
1554 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: points_and_weights
1555 :
1556 18 : CALL timeset(routineN, handle)
1557 :
1558 18 : homo = bs_env%n_occ(1)
1559 18 : n_mo = bs_env%n_ao
1560 18 : num_time_freq_points = bs_env%num_time_freq_points
1561 :
1562 : ! minimum and maximum difference between eigenvalues of unoccupied and an occupied MOs
1563 : E_min = MINVAL(bs_env%eigenval_scf_Gamma(homo + 1, :)) - &
1564 102 : MAXVAL(bs_env%eigenval_scf_Gamma(homo, :))
1565 : E_max = MAXVAL(bs_env%eigenval_scf_Gamma(n_mo, :)) - &
1566 102 : MINVAL(bs_env%eigenval_scf_Gamma(1, :))
1567 :
1568 18 : E_range = E_max/E_min
1569 :
1570 54 : ALLOCATE (points_and_weights(2*num_time_freq_points))
1571 :
1572 : ! frequency points
1573 18 : IF (num_time_freq_points .LE. 20) THEN
1574 18 : CALL get_rpa_minimax_coeff(num_time_freq_points, E_range, points_and_weights, ierr, .FALSE.)
1575 : ELSE
1576 0 : CALL get_rpa_minimax_coeff_larger_grid(num_time_freq_points, E_range, points_and_weights)
1577 : END IF
1578 :
1579 : ! one needs to scale the minimax grids, see Azizi, Wilhelm, Golze, Panades-Barrueta,
1580 : ! Giantomassi, Rinke, Draxl, Gonze et al., 2 publications
1581 190 : bs_env%imag_freq_points(:) = points_and_weights(1:num_time_freq_points)*E_min
1582 :
1583 : ! determine number of fit points in the interval [0,ω_max] for virt, or [-ω_max,0] for occ
1584 18 : bs_env%num_freq_points_fit = 0
1585 190 : DO i_w = 1, bs_env%num_time_freq_points
1586 190 : IF (bs_env%imag_freq_points(i_w) < bs_env%freq_max_fit) THEN
1587 82 : bs_env%num_freq_points_fit = bs_env%num_freq_points_fit + 1
1588 : END IF
1589 : END DO
1590 :
1591 : ! iω values for the analytic continuation Σ^c_n(iω,k) -> Σ^c_n(ϵ,k)
1592 54 : ALLOCATE (bs_env%imag_freq_points_fit(bs_env%num_freq_points_fit))
1593 18 : j_w = 0
1594 190 : DO i_w = 1, bs_env%num_time_freq_points
1595 190 : IF (bs_env%imag_freq_points(i_w) < bs_env%freq_max_fit) THEN
1596 82 : j_w = j_w + 1
1597 82 : bs_env%imag_freq_points_fit(j_w) = bs_env%imag_freq_points(i_w)
1598 : END IF
1599 : END DO
1600 :
1601 : ! reset the number of Padé parameters if smaller than the number of
1602 : ! imaginary-frequency points for the fit
1603 18 : IF (bs_env%num_freq_points_fit < bs_env%nparam_pade) THEN
1604 18 : bs_env%nparam_pade = bs_env%num_freq_points_fit
1605 : END IF
1606 :
1607 : ! time points
1608 18 : IF (num_time_freq_points .LE. 20) THEN
1609 18 : CALL get_exp_minimax_coeff(num_time_freq_points, E_range, points_and_weights)
1610 : ELSE
1611 0 : CALL get_exp_minimax_coeff_gw(num_time_freq_points, E_range, points_and_weights)
1612 : END IF
1613 :
1614 190 : bs_env%imag_time_points(:) = points_and_weights(1:num_time_freq_points)/(2.0_dp*E_min)
1615 :
1616 18 : DEALLOCATE (points_and_weights)
1617 :
1618 : ! in minimax grids, Fourier transforms t -> w and w -> t are split using
1619 : ! e^(iwt) = cos(wt) + i sin(wt); we thus calculate weights for trafos with a cos and
1620 : ! sine prefactor; details in Azizi, Wilhelm, Golze, Giantomassi, Panades-Barrueta,
1621 : ! Rinke, Draxl, Gonze et al., 2 publications
1622 :
1623 : ! cosine transform weights imaginary time to imaginary frequency
1624 : CALL get_l_sq_wghts_cos_tf_t_to_w(num_time_freq_points, &
1625 : bs_env%imag_time_points, &
1626 : bs_env%weights_cos_t_to_w, &
1627 : bs_env%imag_freq_points, &
1628 : E_min, E_max, max_error_min, &
1629 : bs_env%num_points_per_magnitude, &
1630 18 : bs_env%regularization_minimax)
1631 :
1632 : ! cosine transform weights imaginary frequency to imaginary time
1633 : CALL get_l_sq_wghts_cos_tf_w_to_t(num_time_freq_points, &
1634 : bs_env%imag_time_points, &
1635 : bs_env%weights_cos_w_to_t, &
1636 : bs_env%imag_freq_points, &
1637 : E_min, E_max, max_error_min, &
1638 : bs_env%num_points_per_magnitude, &
1639 18 : bs_env%regularization_minimax)
1640 :
1641 : ! sine transform weights imaginary time to imaginary frequency
1642 : CALL get_l_sq_wghts_sin_tf_t_to_w(num_time_freq_points, &
1643 : bs_env%imag_time_points, &
1644 : bs_env%weights_sin_t_to_w, &
1645 : bs_env%imag_freq_points, &
1646 : E_min, E_max, max_error_min, &
1647 : bs_env%num_points_per_magnitude, &
1648 18 : bs_env%regularization_minimax)
1649 :
1650 18 : u = bs_env%unit_nr
1651 18 : IF (u > 0) THEN
1652 9 : WRITE (u, '(T2,A)') ''
1653 9 : WRITE (u, '(T2,A,F44.2)') 'SCF direct band gap at Γ-point (eV)', E_min*evolt
1654 9 : WRITE (u, '(T2,A,F42.2)') 'Max. SCF eigval diff. at Γ-point (eV)', E_max*evolt
1655 9 : WRITE (u, '(T2,A,F55.2)') 'E-Range for minimax grid', E_range
1656 9 : WRITE (u, '(T2,A,I27)') 'Number of Padé parameters for analytic continuation:', &
1657 18 : bs_env%nparam_pade
1658 9 : WRITE (u, '(T2,A)') ''
1659 : END IF
1660 18 : CALL timestop(handle)
1661 :
1662 36 : END SUBROUTINE setup_time_and_frequency_minimax_grid
1663 :
1664 : END MODULE gw_utils
|
