(git:d3d49ac)
Loading...
Searching...
No Matches
gw_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief
10!> \par History
11!> 01.2026 Maximilian Graml: add more bounds to exploit sparsity in 3c integrals, fixes
12!> \author Jan Wilhelm
13!> \date 07.2023
14! **************************************************************************************************
15MODULE gw_utils
16 USE atomic_kind_types, ONLY: atomic_kind_type,&
17 get_atomic_kind_set
18 USE basis_set_types, ONLY: get_gto_basis_set,&
19 gto_basis_set_type
20 USE bibliography, ONLY: graml2024,&
21 cite_reference
22 USE cell_types, ONLY: cell_type,&
23 pbc,&
24 scaled_to_real
25 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
26 cp_blacs_env_release,&
27 cp_blacs_env_type
28 USE cp_cfm_types, ONLY: cp_cfm_create,&
29 cp_cfm_release,&
30 cp_cfm_to_cfm,&
31 cp_cfm_to_fm,&
32 cp_cfm_type
33 USE cp_control_types, ONLY: dft_control_type
34 USE cp_dbcsr_api, ONLY: &
35 dbcsr_create, dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_p_type, &
36 dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
37 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
38 copy_fm_to_dbcsr,&
39 cp_dbcsr_dist2d_to_dist,&
40 dbcsr_allocate_matrix_set,&
41 dbcsr_deallocate_matrix_set
42 USE cp_files, ONLY: close_file,&
43 open_file
44 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
45 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
46 cp_fm_struct_release,&
47 cp_fm_struct_type
48 USE cp_fm_types, ONLY: cp_fm_create,&
49 cp_fm_get_diag,&
50 cp_fm_release,&
51 cp_fm_set_all,&
52 cp_fm_type
53 USE cp_log_handling, ONLY: cp_get_default_logger,&
54 cp_logger_type
55 USE cp_output_handling, ONLY: cp_print_key_generate_filename
56 USE dbt_api, ONLY: &
57 dbt_clear, dbt_create, dbt_destroy, dbt_filter, dbt_iterator_blocks_left, &
58 dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, &
59 dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
60 USE distribution_2d_types, ONLY: distribution_2d_type
61 USE gw_communication, ONLY: fm_to_local_array
62 USE gw_integrals, ONLY: build_3c_integral_block
63 USE input_constants, ONLY: &
64 do_potential_truncated, large_cell_gamma, large_cell_gamma_ri_rs, non_periodic_ri_rs, &
65 ri_rpa_g0w0_crossing_newton, rtp_bse_kernel_ri_ao, rtp_bse_kernel_ri_rs, rtp_method_bse, &
66 rtp_method_bse_linearized, small_cell_full_kp, 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_path_length,&
73 dp,&
74 int_8
75 USE kpoint_k_r_trafo_simple, ONLY: rs_to_kp
76 USE kpoint_types, ONLY: get_kpoint_info,&
77 kpoint_create,&
78 kpoint_type
79 USE libint_2c_3c, ONLY: libint_potential_type
80 USE libint_wrapper, ONLY: cp_libint_static_cleanup,&
81 cp_libint_static_init
82 USE machine, ONLY: m_memory,&
83 m_walltime
84 USE mathconstants, ONLY: gaussi,&
85 z_one,&
86 z_zero
87 USE mathlib, ONLY: diag_complex,&
88 gcd
89 USE message_passing, ONLY: mp_cart_type,&
90 mp_para_env_type
91 USE minimax_exp, ONLY: get_exp_minimax_coeff
92 USE minimax_exp_gw, ONLY: get_exp_minimax_coeff_gw
93 USE minimax_rpa, ONLY: get_rpa_minimax_coeff,&
94 get_rpa_minimax_coeff_larger_grid
95 USE mp2_gpw, ONLY: create_mat_munu
96 USE mp2_grids, ONLY: get_l_sq_wghts_cos_tf_t_to_w,&
97 get_l_sq_wghts_cos_tf_w_to_t,&
98 get_l_sq_wghts_sin_tf_t_to_w
99 USE mp2_ri_2c, ONLY: trunc_coulomb_for_exchange
100 USE parallel_gemm_api, ONLY: parallel_gemm
101 USE particle_methods, ONLY: get_particle_set
102 USE particle_types, ONLY: particle_type
103 USE physcon, ONLY: angstrom,&
104 evolt
105 USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
106 USE post_scf_bandstructure_utils, ONLY: rsmat_to_kp
107 USE qs_energy_types, ONLY: qs_energy_type
108 USE qs_environment_types, ONLY: get_qs_env,&
109 qs_env_part_release,&
110 qs_environment_type
111 USE qs_integral_utils, ONLY: basis_set_list_setup
112 USE qs_interactions, ONLY: init_interaction_radii_orb_basis
113 USE qs_kind_types, ONLY: get_qs_kind,&
114 qs_kind_type
115 USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
116 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
117 release_neighbor_list_sets
118 USE qs_tensors, ONLY: build_2c_integrals,&
119 build_2c_neighbor_lists,&
120 build_3c_integrals,&
121 build_3c_neighbor_lists,&
122 get_tensor_occupancy,&
123 neighbor_list_3c_destroy
124 USE qs_tensors_types, ONLY: create_2c_tensor,&
125 create_3c_tensor,&
126 distribution_3d_create,&
127 distribution_3d_type,&
128 neighbor_list_3c_type
129 USE rpa_gw, ONLY: continuation_pade
130#include "base/base_uses.f90"
131
132 IMPLICIT NONE
133
134 PRIVATE
135
136 PUBLIC :: create_and_init_bs_env_for_gw, de_init_bs_env, get_i_j_atoms, &
137 compute_xkp, time_to_freq, analyt_conti_and_print, &
138 add_r, is_cell_in_index_to_cell, get_v_tr_r, power, &
139 rtbse_resolve_rirs_flag
140
141 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils'
142
143CONTAINS
144
145! **************************************************************************************************
146!> \brief ...
147!> \param qs_env ...
148!> \param bs_env ...
149!> \param bs_sec ...
150! **************************************************************************************************
151 SUBROUTINE create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec)
152 TYPE(qs_environment_type), POINTER :: qs_env
153 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
154 TYPE(section_vals_type), POINTER :: bs_sec
155
156 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_and_init_bs_env_for_gw'
157
158 INTEGER :: handle
159
160 CALL timeset(routinen, handle)
161
162 CALL cite_reference(graml2024)
163
164 CALL read_gw_input_parameters(bs_env, bs_sec)
165
166 CALL print_header_and_input_parameters(bs_env)
167
168 CALL setup_ao_and_ri_basis_set(qs_env, bs_env)
169
170 CALL get_ri_basis_and_basis_function_indices(qs_env, bs_env)
171
172 CALL set_heuristic_parameters(bs_env, qs_env)
173
174 CALL cp_libint_static_init()
175
176 CALL setup_kpoints_chi_eps_w(bs_env, bs_env%kpoints_chi_eps_W)
177
178 IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
179 CALL setup_cells_3c(qs_env, bs_env)
180 END IF
181
182 CALL set_parallelization_parameters(qs_env, bs_env)
183
184 CALL allocate_matrices(qs_env, bs_env)
185
186 CALL compute_v_xc(qs_env, bs_env)
187
188 CALL create_tensors(qs_env, bs_env)
189
190 SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
191 CASE (large_cell_gamma, large_cell_gamma_ri_rs, non_periodic_ri_rs)
192
193 CALL allocate_gw_eigenvalues(bs_env)
194
195 CALL check_sparsity_3c(qs_env, bs_env)
196
197 CALL set_sparsity_parallelization_parameters(bs_env)
198
199 CALL check_for_restart_files(qs_env, bs_env)
200
201 CASE (small_cell_full_kp)
202
203 CALL compute_3c_integrals(qs_env, bs_env)
204
205 CALL setup_cells_delta_r(bs_env)
206
207 CALL setup_parallelization_delta_r(bs_env)
208
209 CALL allocate_matrices_small_cell_full_kp(qs_env, bs_env)
210
211 CALL trafo_v_xc_r_to_kp(qs_env, bs_env)
212
213 CALL heuristic_ri_regularization(qs_env, bs_env)
214
215 END SELECT
216
217 CALL setup_time_and_frequency_minimax_grid(bs_env)
218
219 ! free memory in qs_env; only if one is not calculating the LDOS because
220 ! we need real-space grid operations in pw_env, task_list for the LDOS
221 ! Recommendation in case of memory issues: first perform GW calculation without calculating
222 ! LDOS (to safe memor). Then, use GW restart files
223 ! in a subsequent calculation to calculate the LDOS
224 ! Marek : TODO - boolean that does not interfere with RTP init but sets this to correct value
225 IF (.NOT. bs_env%do_ldos .AND. .false.) THEN
226 CALL qs_env_part_release(qs_env)
227 END IF
228
229 CALL timestop(handle)
230
231 END SUBROUTINE create_and_init_bs_env_for_gw
232
233! **************************************************************************************************
234!> \brief Releases the memory-heavy GW intermediates that cannot be freed in bs_env_release,
235!> retaining the 3c neighbor list only when an AO-RI RT-BSE self-energy still needs it
236!> \param qs_env Quickstep environment - used to resolve the RT-BSE kernel-RI setting
237!> \param bs_env Bandstructure environment whose intermediates are released
238! **************************************************************************************************
239 SUBROUTINE de_init_bs_env(qs_env, bs_env)
240 TYPE(qs_environment_type), POINTER :: qs_env
241 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
242
243 CHARACTER(LEN=*), PARAMETER :: routinen = 'de_init_bs_env'
244
245 INTEGER :: handle
246 LOGICAL :: retain_nl_3c, rirs_kernel
247
248 CALL timeset(routinen, handle)
249 ! deallocate quantities here which:
250 ! 1. cannot be deallocated in bs_env_release due to circular dependencies
251 ! 2. consume a lot of memory and should not be kept until the quantity is
252 ! deallocated in bs_env_release
253
254 ! nl_3c feeds only the AO-RI SEX self-energy (compute_3c_integrals); the RI-RS SEX
255 ! path never reads it, and AO-RI Hartree builds its own blocks. Retain iff AO-RI SEX.
256 retain_nl_3c = .false.
257 IF (ASSOCIATED(bs_env%nl_3c%ij_list) .AND. (bs_env%rtp_method == rtp_method_bse)) THEN
258 CALL rtbse_resolve_rirs_flag(qs_env, bs_env, rirs_kernel=rirs_kernel)
259 retain_nl_3c = .NOT. rirs_kernel
260 END IF
261
262 IF (retain_nl_3c) THEN
263 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, *) "Retaining nl_3c for AO-RI RT-BSE self-energy"
264 ELSE
265 CALL neighbor_list_3c_destroy(bs_env%nl_3c)
266 END IF
267
268 CALL cp_libint_static_cleanup()
269
270 CALL timestop(handle)
271
272 END SUBROUTINE de_init_bs_env
273
274! **************************************************************************************************
275!> \brief Resolve the linRTBSE RI-RS kernel switch from the KERNEL_RI input and the GW default.
276!> \param qs_env Quickstep environment - source of the input section and the RTP method
277!> \param bs_env Bandstructure environment - provides the do_gw_ri_rs default
278!> \param rirs_kernel (optional) .TRUE. if the Hartree + SEX kernels use the RI-RS grid backend
279!> \author Maximilian Graml
280!> \note Single source of truth shared by create_rtbse_env (sets the flag) and de_init_bs_env
281!> (decides whether to retain nl_3c). KERNEL_RI=DEFAULT follows bs_env%do_gw_ri_rs;
282!> RS/AO force; forced .FALSE. for non-linearized (full) RT-BSE (warn on explicit RS).
283! **************************************************************************************************
284 SUBROUTINE rtbse_resolve_rirs_flag(qs_env, bs_env, rirs_kernel)
285 TYPE(qs_environment_type), POINTER :: qs_env
286 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
287 LOGICAL, INTENT(OUT), OPTIONAL :: rirs_kernel
288
289 INTEGER :: kernel_ri
290 LOGICAL :: my_rirs_kernel
291 TYPE(dft_control_type), POINTER :: dft_control
292 TYPE(section_vals_type), POINTER :: input
293
294 NULLIFY (dft_control, input)
295 CALL get_qs_env(qs_env, dft_control=dft_control, input=input)
296
297 CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%KERNEL_RI", &
298 i_val=kernel_ri)
299 SELECT CASE (kernel_ri)
300 CASE (rtp_bse_kernel_ri_rs)
301 my_rirs_kernel = .true.
302 CASE (rtp_bse_kernel_ri_ao)
303 my_rirs_kernel = .false.
304 CASE DEFAULT ! rtp_bse_kernel_ri_default
305 my_rirs_kernel = bs_env%do_gw_ri_rs
306 END SELECT
307
308 ! RI-RS kernels are implemented for linearized RT-BSE only; full RT-BSE always uses AO-RI.
309 IF (dft_control%rtp_control%rtp_method /= rtp_method_bse_linearized) THEN
310 IF (kernel_ri == rtp_bse_kernel_ri_rs) THEN
311 cpwarn("RI-RS kernels are implemented for linearized RT-BSE only; forcing AO")
312 END IF
313 my_rirs_kernel = .false.
314 END IF
315
316 IF (PRESENT(rirs_kernel)) rirs_kernel = my_rirs_kernel
317
318 END SUBROUTINE rtbse_resolve_rirs_flag
319
320! **************************************************************************************************
321!> \brief ...
322!> \param bs_env ...
323!> \param bs_sec ...
324! **************************************************************************************************
325 SUBROUTINE read_gw_input_parameters(bs_env, bs_sec)
326 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
327 TYPE(section_vals_type), POINTER :: bs_sec
328
329 CHARACTER(LEN=*), PARAMETER :: routinen = 'read_gw_input_parameters'
330
331 INTEGER :: handle
332 TYPE(section_vals_type), POINTER :: gw_sec
333
334 CALL timeset(routinen, handle)
335
336 NULLIFY (gw_sec)
337 gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
338
339 CALL section_vals_val_get(gw_sec, "NUM_TIME_FREQ_POINTS", i_val=bs_env%num_time_freq_points)
340 CALL section_vals_val_get(gw_sec, "EPS_FILTER", r_val=bs_env%eps_filter)
341 CALL section_vals_val_get(gw_sec, "REGULARIZATION_RI", r_val=bs_env%input_regularization_RI)
342 CALL section_vals_val_get(gw_sec, "REGULARIZATION_MINIMAX", r_val=bs_env%input_regularization_minimax)
343 CALL section_vals_val_get(gw_sec, "CUTOFF_RADIUS_RI", r_val=bs_env%ri_metric%cutoff_radius)
344 CALL section_vals_val_get(gw_sec, "MEMORY_PER_PROC", r_val=bs_env%input_memory_per_proc_GB)
345 CALL section_vals_val_get(gw_sec, "APPROX_KP_EXTRAPOL", l_val=bs_env%approx_kp_extrapol)
346 CALL section_vals_val_get(gw_sec, "SIZE_LATTICE_SUM", i_val=bs_env%size_lattice_sum_V)
347 CALL section_vals_val_get(gw_sec, "KPOINTS_W", i_vals=bs_env%nkp_grid_chi_eps_W_input)
348 CALL section_vals_val_get(gw_sec, "HEDIN_SHIFT", l_val=bs_env%do_hedin_shift)
349 CALL section_vals_val_get(gw_sec, "FREQ_MAX_FIT", r_val=bs_env%freq_max_fit)
350 CALL section_vals_val_get(gw_sec, "PRINT%PRINT_DBT_CONTRACT", l_val=bs_env%print_contract)
351 CALL section_vals_val_get(gw_sec, "PRINT%PRINT_DBT_CONTRACT_VERBOSE", l_val=bs_env%print_contract_verbose)
352 CALL section_vals_val_get(gw_sec, "TIKHONOV", r_val=bs_env%ri_rs%tikhonov)
353 CALL section_vals_val_get(gw_sec, "GRID_SELECT", i_val=bs_env%ri_rs%grid_select)
354 CALL section_vals_val_get(gw_sec, "GRID_FILE_SUFFIX", c_val=bs_env%ri_rs%grid_file_suffix)
355 CALL section_vals_val_get(gw_sec, "CUTOFF_RADIUS_RL_RI", r_val=bs_env%ri_rs%cutoff_radius_ri_rs)
356 CALL section_vals_val_get(gw_sec, "CUTOFF_RADIUS_RL_AO", r_val=bs_env%ri_rs%cutoff_radius_ri_ao)
357 CALL section_vals_val_get(gw_sec, "N_PROCS_PER_ATOM_Z_LP", i_val=bs_env%ri_rs%n_procs_per_atom_z_lp)
358 CALL section_vals_val_get(gw_sec, "N_PANELS", i_val=bs_env%ri_rs%n_panels)
359 CALL section_vals_val_get(gw_sec, "KEEP_SPARSITY_RL", l_val=bs_env%ri_rs%keep_sparsity_rirs)
360 CALL section_vals_val_get(gw_sec, "CUTOFF_RADIUS_RL_W", r_val=bs_env%ri_rs%cutoff_radius_v_w)
361 CALL section_vals_val_get(gw_sec, "CUTOFF_RADIUS_G_W", r_val=bs_env%ri_rs%cutoff_radius_g_w)
362
363 IF (bs_env%print_contract) THEN
364 bs_env%unit_nr_contract = bs_env%unit_nr
365 ELSE
366 bs_env%unit_nr_contract = 0
367 END IF
368 CALL timestop(handle)
369
370 END SUBROUTINE read_gw_input_parameters
371
372! **************************************************************************************************
373!> \brief ...
374!> \param qs_env ...
375!> \param bs_env ...
376! **************************************************************************************************
377 SUBROUTINE setup_ao_and_ri_basis_set(qs_env, bs_env)
378 TYPE(qs_environment_type), POINTER :: qs_env
379 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
380
381 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_AO_and_RI_basis_set'
382
383 INTEGER :: handle, natom, nkind
384 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
385 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
386
387 CALL timeset(routinen, handle)
388
389 CALL get_qs_env(qs_env, &
390 qs_kind_set=qs_kind_set, &
391 particle_set=particle_set, &
392 natom=natom, nkind=nkind)
393
394 ! set up basis
395 ALLOCATE (bs_env%sizes_RI(natom), bs_env%sizes_AO(natom))
396 ALLOCATE (bs_env%basis_set_RI(nkind), bs_env%basis_set_AO(nkind))
397
398 CALL basis_set_list_setup(bs_env%basis_set_RI, "RI_AUX", qs_kind_set)
399 CALL basis_set_list_setup(bs_env%basis_set_AO, "ORB", qs_kind_set)
400
401 CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_RI, &
402 basis=bs_env%basis_set_RI)
403 CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_AO, &
404 basis=bs_env%basis_set_AO)
405
406 CALL timestop(handle)
407
408 END SUBROUTINE setup_ao_and_ri_basis_set
409
410! **************************************************************************************************
411!> \brief ...
412!> \param qs_env ...
413!> \param bs_env ...
414! **************************************************************************************************
415 SUBROUTINE get_ri_basis_and_basis_function_indices(qs_env, bs_env)
416 TYPE(qs_environment_type), POINTER :: qs_env
417 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
418
419 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_RI_basis_and_basis_function_indices'
420
421 INTEGER :: handle, i_ri, iatom, ikind, iset, &
422 max_ao_bf_per_atom, n_ao_test, n_atom, &
423 n_kind, n_ri, nset, nsgf, u
424 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
425 INTEGER, DIMENSION(:), POINTER :: l_max, l_min, nsgf_set
426 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
427 TYPE(gto_basis_set_type), POINTER :: basis_set_a
428 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
429
430 CALL timeset(routinen, handle)
431
432 ! determine RI basis set size
433 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
434
435 n_kind = SIZE(qs_kind_set)
436 n_atom = bs_env%n_atom
437
438 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
439
440 DO ikind = 1, n_kind
441 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
442 basis_type="RI_AUX")
443 IF (.NOT. ASSOCIATED(basis_set_a)) THEN
444 CALL cp_abort(__location__, &
445 "At least one RI_AUX basis set was not explicitly invoked in &KIND-section.")
446 END IF
447 END DO
448
449 ALLOCATE (bs_env%i_RI_start_from_atom(n_atom))
450 ALLOCATE (bs_env%i_RI_end_from_atom(n_atom))
451 ALLOCATE (bs_env%i_ao_start_from_atom(n_atom))
452 ALLOCATE (bs_env%i_ao_end_from_atom(n_atom))
453
454 n_ri = 0
455 DO iatom = 1, n_atom
456 bs_env%i_RI_start_from_atom(iatom) = n_ri + 1
457 ikind = kind_of(iatom)
458 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
459 n_ri = n_ri + nsgf
460 bs_env%i_RI_end_from_atom(iatom) = n_ri
461 END DO
462 bs_env%n_RI = n_ri
463
464 max_ao_bf_per_atom = 0
465 n_ao_test = 0
466 DO iatom = 1, n_atom
467 bs_env%i_ao_start_from_atom(iatom) = n_ao_test + 1
468 ikind = kind_of(iatom)
469 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="ORB")
470 n_ao_test = n_ao_test + nsgf
471 bs_env%i_ao_end_from_atom(iatom) = n_ao_test
472 max_ao_bf_per_atom = max(max_ao_bf_per_atom, nsgf)
473 END DO
474 cpassert(n_ao_test == bs_env%n_ao)
475 bs_env%max_AO_bf_per_atom = max_ao_bf_per_atom
476
477 ALLOCATE (bs_env%l_RI(n_ri))
478 i_ri = 0
479 DO iatom = 1, n_atom
480 ikind = kind_of(iatom)
481
482 nset = bs_env%basis_set_RI(ikind)%gto_basis_set%nset
483 l_max => bs_env%basis_set_RI(ikind)%gto_basis_set%lmax
484 l_min => bs_env%basis_set_RI(ikind)%gto_basis_set%lmin
485 nsgf_set => bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf_set
486
487 DO iset = 1, nset
488 cpassert(l_max(iset) == l_min(iset))
489 bs_env%l_RI(i_ri + 1:i_ri + nsgf_set(iset)) = l_max(iset)
490 i_ri = i_ri + nsgf_set(iset)
491 END DO
492
493 END DO
494 cpassert(i_ri == n_ri)
495
496 u = bs_env%unit_nr
497
498 IF (u > 0) THEN
499 WRITE (u, fmt="(T2,A)") " "
500 WRITE (u, fmt="(T2,2A,T75,I8)") "Number of auxiliary Gaussian basis functions ", &
501 χε"for , , W", n_ri
502 END IF
503
504 CALL timestop(handle)
505
506 END SUBROUTINE get_ri_basis_and_basis_function_indices
507
508! **************************************************************************************************
509!> \brief ...
510!> \param bs_env ...
511!> \param kpoints ...
512! **************************************************************************************************
513 SUBROUTINE setup_kpoints_chi_eps_w(bs_env, kpoints)
514
515 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
516 TYPE(kpoint_type), POINTER :: kpoints
517
518 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_kpoints_chi_eps_W'
519
520 INTEGER :: handle, i_dim, n_dim, nkp, nkp_extra, &
521 nkp_orig, u
522 INTEGER, DIMENSION(3) :: nkp_grid, nkp_grid_extra, periodic
523 REAL(kind=dp) :: exp_s_p, n_dim_inv
524
525 CALL timeset(routinen, handle)
526
527 ! routine adapted from mp2_integrals.F
528 NULLIFY (kpoints)
529 CALL kpoint_create(kpoints)
530
531 kpoints%kp_scheme = "GENERAL"
532
533 periodic(1:3) = bs_env%periodic(1:3)
534
535 cpassert(SIZE(bs_env%nkp_grid_chi_eps_W_input) == 3)
536
537 IF (bs_env%nkp_grid_chi_eps_W_input(1) > 0 .AND. &
538 bs_env%nkp_grid_chi_eps_W_input(2) > 0 .AND. &
539 bs_env%nkp_grid_chi_eps_W_input(3) > 0) THEN
540 ! 1. k-point mesh for χ, ε, W from input
541 DO i_dim = 1, 3
542 SELECT CASE (periodic(i_dim))
543 CASE (0)
544 nkp_grid(i_dim) = 1
545 nkp_grid_extra(i_dim) = 1
546 CASE (1)
547 nkp_grid(i_dim) = bs_env%nkp_grid_chi_eps_W_input(i_dim)
548 nkp_grid_extra(i_dim) = nkp_grid(i_dim)*2
549 CASE DEFAULT
550 cpabort("Error in periodicity.")
551 END SELECT
552 END DO
553
554 ELSE IF (bs_env%nkp_grid_chi_eps_W_input(1) == -1 .AND. &
555 bs_env%nkp_grid_chi_eps_W_input(2) == -1 .AND. &
556 bs_env%nkp_grid_chi_eps_W_input(3) == -1) THEN
557 ! 2. automatic k-point mesh for χ, ε, W
558
559 DO i_dim = 1, 3
560
561 cpassert(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
562
563 SELECT CASE (periodic(i_dim))
564 CASE (0)
565 nkp_grid(i_dim) = 1
566 nkp_grid_extra(i_dim) = 1
567 CASE (1)
568 SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
569 CASE (large_cell_gamma, large_cell_gamma_ri_rs, non_periodic_ri_rs)
570 nkp_grid(i_dim) = 4
571 nkp_grid_extra(i_dim) = 6
572 CASE (small_cell_full_kp)
573 nkp_grid(i_dim) = bs_env%kpoints_scf_desymm%nkp_grid(i_dim)*4
574 nkp_grid_extra(i_dim) = bs_env%kpoints_scf_desymm%nkp_grid(i_dim)*8
575 END SELECT
576 CASE DEFAULT
577 cpabort("Error in periodicity.")
578 END SELECT
579
580 END DO
581
582 ELSE
583
584 cpabort("An error occured when setting up the k-mesh for W.")
585
586 END IF
587
588 nkp_orig = max(nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2, 1)
589
590 nkp_extra = nkp_grid_extra(1)*nkp_grid_extra(2)*nkp_grid_extra(3)/2
591
592 nkp = nkp_orig + nkp_extra
593
594 kpoints%nkp_grid(1:3) = nkp_grid(1:3)
595 kpoints%nkp = nkp
596
597 bs_env%nkp_grid_chi_eps_W_orig(1:3) = nkp_grid(1:3)
598 bs_env%nkp_grid_chi_eps_W_extra(1:3) = nkp_grid_extra(1:3)
599 bs_env%nkp_chi_eps_W_orig = nkp_orig
600 bs_env%nkp_chi_eps_W_extra = nkp_extra
601 bs_env%nkp_chi_eps_W_orig_plus_extra = nkp
602
603 ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
604 ALLOCATE (bs_env%wkp_no_extra(nkp), bs_env%wkp_s_p(nkp))
605
606 CALL compute_xkp(kpoints%xkp, 1, nkp_orig, nkp_grid)
607 CALL compute_xkp(kpoints%xkp, nkp_orig + 1, nkp, nkp_grid_extra)
608
609 n_dim = sum(periodic)
610 IF (n_dim == 0) THEN
611 ! molecules
612 kpoints%wkp(1) = 1.0_dp
613 bs_env%wkp_s_p(1) = 1.0_dp
614 bs_env%wkp_no_extra(1) = 1.0_dp
615 ELSE
616
617 n_dim_inv = 1.0_dp/real(n_dim, kind=dp)
618
619 ! k-point weights are chosen to automatically extrapolate the k-point mesh
620 CALL compute_wkp(kpoints%wkp(1:nkp_orig), nkp_orig, nkp_extra, n_dim_inv)
621 CALL compute_wkp(kpoints%wkp(nkp_orig + 1:nkp), nkp_extra, nkp_orig, n_dim_inv)
622
623 bs_env%wkp_no_extra(1:nkp_orig) = 0.0_dp
624 bs_env%wkp_no_extra(nkp_orig + 1:nkp) = 1.0_dp/real(nkp_extra, kind=dp)
625
626 IF (n_dim == 3) THEN
627 ! W_PQ(k) for an s-function P and a p-function Q diverges as 1/k at k=0
628 ! (instead of 1/k^2 for P and Q both being s-functions).
629 exp_s_p = 2.0_dp*n_dim_inv
630 CALL compute_wkp(bs_env%wkp_s_p(1:nkp_orig), nkp_orig, nkp_extra, exp_s_p)
631 CALL compute_wkp(bs_env%wkp_s_p(nkp_orig + 1:nkp), nkp_extra, nkp_orig, exp_s_p)
632 ELSE
633 bs_env%wkp_s_p(1:nkp) = bs_env%wkp_no_extra(1:nkp)
634 END IF
635
636 END IF
637
638 IF (bs_env%approx_kp_extrapol) THEN
639 bs_env%wkp_orig = 1.0_dp/real(nkp_orig, kind=dp)
640 END IF
641
642 ! heuristic parameter: how many k-points for χ, ε, and W are used simultaneously
643 ! (less simultaneous k-points: less memory, but more computational effort because of
644 ! recomputation of V(k))
645 bs_env%nkp_chi_eps_W_batch = 4
646
647 bs_env%num_chi_eps_W_batches = (bs_env%nkp_chi_eps_W_orig_plus_extra - 1)/ &
648 bs_env%nkp_chi_eps_W_batch + 1
649
650 u = bs_env%unit_nr
651
652 IF (u > 0) THEN
653 WRITE (u, fmt="(T2,A)") " "
654 WRITE (u, fmt="(T2,1A,T71,3I4)") χε"K-point mesh 1 for , , W", nkp_grid(1:3)
655 WRITE (u, fmt="(T2,2A,T71,3I4)") χε"K-point mesh 2 for , , W ", &
656 "(for k-point extrapolation of W)", nkp_grid_extra(1:3)
657 WRITE (u, fmt="(T2,A,T80,L)") "Approximate the k-point extrapolation", &
658 bs_env%approx_kp_extrapol
659 END IF
660
661 CALL timestop(handle)
662
663 END SUBROUTINE setup_kpoints_chi_eps_w
664
665! **************************************************************************************************
666!> \brief ...
667!> \param xkp ...
668!> \param ikp_start ...
669!> \param ikp_end ...
670!> \param grid ...
671! **************************************************************************************************
672 SUBROUTINE compute_xkp(xkp, ikp_start, ikp_end, grid)
673
674 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
675 INTEGER :: ikp_start, ikp_end
676 INTEGER, DIMENSION(3) :: grid
677
678 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_xkp'
679
680 INTEGER :: handle, i, ix, iy, iz
681
682 CALL timeset(routinen, handle)
683
684 i = ikp_start
685 DO ix = 1, grid(1)
686 DO iy = 1, grid(2)
687 DO iz = 1, grid(3)
688
689 IF (i > ikp_end) cycle
690
691 xkp(1, i) = real(2*ix - grid(1) - 1, kind=dp)/(2._dp*real(grid(1), kind=dp))
692 xkp(2, i) = real(2*iy - grid(2) - 1, kind=dp)/(2._dp*real(grid(2), kind=dp))
693 xkp(3, i) = real(2*iz - grid(3) - 1, kind=dp)/(2._dp*real(grid(3), kind=dp))
694 i = i + 1
695
696 END DO
697 END DO
698 END DO
699
700 CALL timestop(handle)
701
702 END SUBROUTINE compute_xkp
703
704! **************************************************************************************************
705!> \brief ...
706!> \param wkp ...
707!> \param nkp_1 ...
708!> \param nkp_2 ...
709!> \param exponent ...
710! **************************************************************************************************
711 SUBROUTINE compute_wkp(wkp, nkp_1, nkp_2, exponent)
712 REAL(kind=dp), DIMENSION(:) :: wkp
713 INTEGER :: nkp_1, nkp_2
714 REAL(kind=dp) :: exponent
715
716 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_wkp'
717
718 INTEGER :: handle
719 REAL(kind=dp) :: nkp_ratio
720
721 CALL timeset(routinen, handle)
722
723 nkp_ratio = real(nkp_2, kind=dp)/real(nkp_1, kind=dp)
724
725 wkp(:) = 1.0_dp/real(nkp_1, kind=dp)/(1.0_dp - nkp_ratio**exponent)
726
727 CALL timestop(handle)
728
729 END SUBROUTINE compute_wkp
730
731! **************************************************************************************************
732!> \brief ...
733!> \param qs_env ...
734!> \param bs_env ...
735! **************************************************************************************************
736 SUBROUTINE allocate_matrices(qs_env, bs_env)
737 TYPE(qs_environment_type), POINTER :: qs_env
738 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
739
740 CHARACTER(LEN=*), PARAMETER :: routinen = 'allocate_matrices'
741
742 INTEGER :: handle, i_t
743 TYPE(cp_blacs_env_type), POINTER :: blacs_env, blacs_env_tensor
744 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_ri_global
745 TYPE(mp_para_env_type), POINTER :: para_env
746
747 CALL timeset(routinen, handle)
748
749 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
750
751 fm_struct => bs_env%fm_ks_Gamma(1)%matrix_struct
752
753 CALL cp_fm_create(bs_env%fm_Gocc, fm_struct)
754 CALL cp_fm_create(bs_env%fm_Gvir, fm_struct)
755
756 NULLIFY (fm_struct_ri_global)
757 CALL cp_fm_struct_create(fm_struct_ri_global, context=blacs_env, nrow_global=bs_env%n_RI, &
758 ncol_global=bs_env%n_RI, para_env=para_env)
759 CALL cp_fm_create(bs_env%fm_RI_RI, fm_struct_ri_global)
760 CALL cp_fm_create(bs_env%fm_chi_Gamma_freq, fm_struct_ri_global)
761 CALL cp_fm_create(bs_env%fm_W_MIC_freq, fm_struct_ri_global)
762 IF (bs_env%approx_kp_extrapol) THEN
763 CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_extra, fm_struct_ri_global)
764 CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_no_extra, fm_struct_ri_global)
765 CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_extra, 0.0_dp)
766 CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_no_extra, 0.0_dp)
767 END IF
768 CALL cp_fm_struct_release(fm_struct_ri_global)
769
770 ! create blacs_env for subgroups of tensor operations
771 NULLIFY (blacs_env_tensor)
772 CALL cp_blacs_env_create(blacs_env=blacs_env_tensor, para_env=bs_env%para_env_tensor)
773
774 ! allocate dbcsr matrices in the tensor subgroup; actually, one only needs a small
775 ! subset of blocks in the tensor subgroup, however, all atomic blocks are allocated.
776 ! One might think of creating a dbcsr matrix with only the blocks that are needed
777 ! in the tensor subgroup
778 CALL create_mat_munu(bs_env%mat_ao_ao_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
779 blacs_env_tensor, do_ri_aux_basis=.false.)
780
781 CALL create_mat_munu(bs_env%mat_RI_RI_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
782 blacs_env_tensor, do_ri_aux_basis=.true.)
783
784 CALL create_mat_munu(bs_env%mat_RI_RI, qs_env, bs_env%eps_atom_grid_2d_mat, &
785 blacs_env, do_ri_aux_basis=.true.)
786
787 CALL cp_blacs_env_release(blacs_env_tensor)
788
789 NULLIFY (bs_env%mat_chi_Gamma_tau)
790 CALL dbcsr_allocate_matrix_set(bs_env%mat_chi_Gamma_tau, bs_env%num_time_freq_points)
791
792 DO i_t = 1, bs_env%num_time_freq_points
793 ALLOCATE (bs_env%mat_chi_Gamma_tau(i_t)%matrix)
794 CALL dbcsr_create(bs_env%mat_chi_Gamma_tau(i_t)%matrix, template=bs_env%mat_RI_RI%matrix)
795 END DO
796
797 CALL timestop(handle)
798
799 END SUBROUTINE allocate_matrices
800
801! **************************************************************************************************
802!> \brief ...
803!> \param bs_env ...
804! **************************************************************************************************
805 SUBROUTINE allocate_gw_eigenvalues(bs_env)
806 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
807
808 CHARACTER(LEN=*), PARAMETER :: routinen = 'allocate_GW_eigenvalues'
809
810 INTEGER :: handle
811
812 CALL timeset(routinen, handle)
813
814 ALLOCATE (bs_env%eigenval_G0W0(bs_env%n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
815 ALLOCATE (bs_env%eigenval_HF(bs_env%n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
816
817 CALL timestop(handle)
818
819 END SUBROUTINE allocate_gw_eigenvalues
820
821! **************************************************************************************************
822!> \brief ...
823!> \param qs_env ...
824!> \param bs_env ...
825! **************************************************************************************************
826 SUBROUTINE create_tensors(qs_env, bs_env)
827 TYPE(qs_environment_type), POINTER :: qs_env
828 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
829
830 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_tensors'
831
832 INTEGER :: handle
833
834 CALL timeset(routinen, handle)
835
836 CALL init_interaction_radii(bs_env)
837
838 ! split blocks does not improve load balancing/efficienfy for tensor contraction, so we go
839 ! with the standard atomic blocks
840 CALL create_3c_t(bs_env%t_RI_AO__AO, bs_env%para_env_tensor, "(RI AO | AO)", [1, 2], [3], &
841 bs_env%sizes_RI, bs_env%sizes_AO, &
842 create_nl_3c=.true., nl_3c=bs_env%nl_3c, qs_env=qs_env)
843 CALL create_3c_t(bs_env%t_RI__AO_AO, bs_env%para_env_tensor, "(RI | AO AO)", [1], [2, 3], &
844 bs_env%sizes_RI, bs_env%sizes_AO)
845
846 CALL create_2c_t(bs_env, bs_env%sizes_RI, bs_env%sizes_AO)
847
848 CALL timestop(handle)
849
850 END SUBROUTINE create_tensors
851
852! **************************************************************************************************
853!> \brief ...
854!> \param qs_env ...
855!> \param bs_env ...
856! **************************************************************************************************
857 SUBROUTINE check_sparsity_3c(qs_env, bs_env)
858 TYPE(qs_environment_type), POINTER :: qs_env
859 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
860
861 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_sparsity_3c'
862
863 INTEGER :: handle, n_atom_step, ri_atom
864 INTEGER(int_8) :: non_zero_elements_sum, nze
865 REAL(dp) :: max_dist_ao_atoms, occ, occupation_sum
866 REAL(kind=dp) :: t1, t2
867 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_global_array
868
869!TYPE(dbt_type) :: t_3c_global
870
871 !TYPE(neighbor_list_3c_type) :: nl_3c_global
872
873 CALL timeset(routinen, handle)
874
875 ! check the sparsity of 3c integral tensor (µν|P); calculate maximum distance between
876 ! AO atoms µ, ν where at least a single integral (µν|P) is larger than the filter threshold
877
878 ALLOCATE (t_3c_global_array(1, 1))
879 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_global_array(1, 1))
880
881 ! Allocate arrays to store min/max indices for overlap with other AO/RI functions on each atom
882 ! (Filled during loop via get_i_j_atom_ranges)
883 ALLOCATE (bs_env%min_RI_idx_from_AO_AO_atom(bs_env%n_atom, bs_env%n_atom))
884 ALLOCATE (bs_env%max_RI_idx_from_AO_AO_atom(bs_env%n_atom, bs_env%n_atom))
885 ALLOCATE (bs_env%min_AO_idx_from_RI_AO_atom(bs_env%n_atom, bs_env%n_atom))
886 ALLOCATE (bs_env%max_AO_idx_from_RI_AO_atom(bs_env%n_atom, bs_env%n_atom))
887 bs_env%min_RI_idx_from_AO_AO_atom(:, :) = bs_env%n_RI
888 bs_env%max_RI_idx_from_AO_AO_atom(:, :) = 1
889 bs_env%min_AO_idx_from_RI_AO_atom(:, :) = bs_env%n_AO
890 bs_env%max_AO_idx_from_RI_AO_atom(:, :) = 1
891
892 CALL bs_env%para_env%sync()
893 t1 = m_walltime()
894
895 occupation_sum = 0.0_dp
896 non_zero_elements_sum = 0
897 max_dist_ao_atoms = 0.0_dp
898 n_atom_step = int(sqrt(real(bs_env%n_atom, kind=dp)))
899 ! do not compute full 3c integrals at once because it may cause out of memory
900 DO ri_atom = 1, bs_env%n_atom, n_atom_step
901
902 CALL build_3c_integrals(t_3c_global_array, &
903 bs_env%eps_filter, &
904 qs_env, &
905 bs_env%nl_3c, &
906 int_eps=bs_env%eps_filter, &
907 basis_i=bs_env%basis_set_RI, &
908 basis_j=bs_env%basis_set_AO, &
909 basis_k=bs_env%basis_set_AO, &
910 bounds_i=[ri_atom, min(ri_atom + n_atom_step - 1, bs_env%n_atom)], &
911 potential_parameter=bs_env%ri_metric, &
912 desymmetrize=.false.)
913
914 CALL dbt_filter(t_3c_global_array(1, 1), bs_env%eps_filter)
915
916 CALL bs_env%para_env%sync()
917
918 CALL get_tensor_occupancy(t_3c_global_array(1, 1), nze, occ)
919 non_zero_elements_sum = non_zero_elements_sum + nze
920 occupation_sum = occupation_sum + occ
921
922 CALL get_max_dist_ao_atoms(t_3c_global_array(1, 1), max_dist_ao_atoms, qs_env)
923
924 ! Extract indices per block
925 CALL get_i_j_atom_ranges(t_3c_global_array(1, 1), bs_env)
926
927 CALL dbt_clear(t_3c_global_array(1, 1))
928
929 END DO
930
931 t2 = m_walltime()
932
933 ! Sync/max for max_dist_AO_atoms is done inside each get_max_dist_AO_atoms
934 bs_env%max_dist_AO_atoms = max_dist_ao_atoms
935 ! occupation_sum is a global quantity, also needs no sync here
936 bs_env%occupation_3c_int = occupation_sum
937
938 CALL bs_env%para_env%min(bs_env%min_RI_idx_from_AO_AO_atom)
939 CALL bs_env%para_env%max(bs_env%max_RI_idx_from_AO_AO_atom)
940 CALL bs_env%para_env%min(bs_env%min_AO_idx_from_RI_AO_atom)
941 CALL bs_env%para_env%max(bs_env%max_AO_idx_from_RI_AO_atom)
942
943 CALL dbt_destroy(t_3c_global_array(1, 1))
944 DEALLOCATE (t_3c_global_array)
945
946 IF (bs_env%unit_nr > 0) THEN
947 WRITE (bs_env%unit_nr, '(T2,A)') ''
948 WRITE (bs_env%unit_nr, '(T2,A,F27.1,A)') &
949 µν'Computed 3-center integrals (|P), execution time', t2 - t1, ' s'
950 WRITE (bs_env%unit_nr, '(T2,A,F48.3,A)') µν'Percentage of non-zero (|P)', &
951 bs_env%occupation_3c_int*100, ' %'
952 WRITE (bs_env%unit_nr, '(T2,A,F33.1,A)') µνµν'Max. distance between , in non-zero (|P)', &
953 bs_env%max_dist_AO_atoms*angstrom, ' A'
954 WRITE (bs_env%unit_nr, '(T2,2A,I20,A)') 'Required memory if storing all 3-center ', &
955 µν'integrals (|P)', int(real(non_zero_elements_sum, kind=dp)*8.0e-9_dp), ' GB'
956 END IF
957
958 CALL timestop(handle)
959
960 END SUBROUTINE check_sparsity_3c
961
962! **************************************************************************************************
963!> \brief ...
964!> \param t_3c ...
965!> \param bs_env ...
966! **************************************************************************************************
967 SUBROUTINE get_i_j_atom_ranges(t_3c, bs_env)
968 TYPE(dbt_type) :: t_3c
969 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
970
971 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_i_j_atom_ranges'
972
973 INTEGER :: handle, idx_ao_end, idx_ao_start, &
974 idx_ri_end, idx_ri_start
975 INTEGER, DIMENSION(3) :: atom_ind
976 TYPE(dbt_iterator_type) :: iter
977
978 CALL timeset(routinen, handle)
979
980 ! Loop over blocks in 3c, for given min_atom: RI_min/max index from min_atom
981!$OMP PARALLEL DEFAULT(NONE) &
982!$OMP SHARED(t_3c, bs_env) &
983!$OMP PRIVATE(iter, atom_ind, &
984!$OMP idx_RI_start, idx_RI_end, idx_AO_start, idx_AO_end)
985
986 CALL dbt_iterator_start(iter, t_3c)
987 DO WHILE (dbt_iterator_blocks_left(iter))
988 CALL dbt_iterator_next_block(iter, atom_ind)
989
990 ! Pre-fetch indices to avoid referencing 'bs_env' twice inside the ATOMIC blocks
991 idx_ri_start = bs_env%i_RI_start_from_atom(atom_ind(1))
992 idx_ri_end = bs_env%i_RI_end_from_atom(atom_ind(1))
993
994 idx_ao_start = bs_env%i_ao_start_from_atom(atom_ind(2))
995 idx_ao_end = bs_env%i_ao_end_from_atom(atom_ind(2))
996
997 ! Update values safely inside ATOMIC blocks, otherwise race conditions occur
998!$OMP ATOMIC UPDATE
999 bs_env%min_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)) = &
1000 min(bs_env%min_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)), idx_ri_start)
1001!$OMP ATOMIC UPDATE
1002 bs_env%max_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)) = &
1003 max(bs_env%max_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)), idx_ri_end)
1004
1005!$OMP ATOMIC UPDATE
1006 bs_env%min_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)) = &
1007 min(bs_env%min_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)), idx_ao_start)
1008!$OMP ATOMIC UPDATE
1009 bs_env%max_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)) = &
1010 max(bs_env%max_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)), idx_ao_end)
1011
1012 END DO
1013 CALL dbt_iterator_stop(iter)
1014!$OMP END PARALLEL
1015
1016 CALL timestop(handle)
1017
1018 END SUBROUTINE get_i_j_atom_ranges
1019
1020! **************************************************************************************************
1021!> \brief ...
1022!> \param bs_env ...
1023!> \param sizes_RI ...
1024!> \param sizes_AO ...
1025! **************************************************************************************************
1026 SUBROUTINE create_2c_t(bs_env, sizes_RI, sizes_AO)
1027 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1028 INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_ri, sizes_ao
1029
1030 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_2c_t'
1031
1032 INTEGER :: handle
1033 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2
1034 INTEGER, DIMENSION(2) :: pdims_2d
1035 TYPE(dbt_pgrid_type) :: pgrid_2d
1036
1037 CALL timeset(routinen, handle)
1038
1039 ! inspired from rpa_im_time.F / hfx_types.F
1040
1041 pdims_2d = 0
1042 CALL dbt_pgrid_create(bs_env%para_env_tensor, pdims_2d, pgrid_2d)
1043
1044 CALL create_2c_tensor(bs_env%t_G, dist_1, dist_2, pgrid_2d, sizes_ao, sizes_ao, &
1045 name="(AO | AO)")
1046 DEALLOCATE (dist_1, dist_2)
1047 CALL create_2c_tensor(bs_env%t_chi, dist_1, dist_2, pgrid_2d, sizes_ri, sizes_ri, &
1048 name="(RI | RI)")
1049 DEALLOCATE (dist_1, dist_2)
1050 CALL create_2c_tensor(bs_env%t_W, dist_1, dist_2, pgrid_2d, sizes_ri, sizes_ri, &
1051 name="(RI | RI)")
1052 DEALLOCATE (dist_1, dist_2)
1053 CALL dbt_pgrid_destroy(pgrid_2d)
1054
1055 CALL timestop(handle)
1056
1057 END SUBROUTINE create_2c_t
1058
1059! **************************************************************************************************
1060!> \brief ...
1061!> \param tensor ...
1062!> \param para_env ...
1063!> \param tensor_name ...
1064!> \param map1 ...
1065!> \param map2 ...
1066!> \param sizes_RI ...
1067!> \param sizes_AO ...
1068!> \param create_nl_3c ...
1069!> \param nl_3c ...
1070!> \param qs_env ...
1071! **************************************************************************************************
1072 SUBROUTINE create_3c_t(tensor, para_env, tensor_name, map1, map2, sizes_RI, sizes_AO, &
1073 create_nl_3c, nl_3c, qs_env)
1074 TYPE(dbt_type) :: tensor
1075 TYPE(mp_para_env_type), POINTER :: para_env
1076 CHARACTER(LEN=12) :: tensor_name
1077 INTEGER, DIMENSION(:) :: map1, map2
1078 INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_ri, sizes_ao
1079 LOGICAL, OPTIONAL :: create_nl_3c
1080 TYPE(neighbor_list_3c_type), OPTIONAL :: nl_3c
1081 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
1082
1083 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_3c_t'
1084
1085 INTEGER :: handle, nkind
1086 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_ao_1, dist_ao_2, dist_ri
1087 INTEGER, DIMENSION(3) :: pcoord, pdims, pdims_3d
1088 LOGICAL :: my_create_nl_3c
1089 TYPE(dbt_pgrid_type) :: pgrid_3d
1090 TYPE(distribution_3d_type) :: dist_3d
1091 TYPE(mp_cart_type) :: mp_comm_t3c_2
1092 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1093
1094 CALL timeset(routinen, handle)
1095
1096 pdims_3d = 0
1097 CALL dbt_pgrid_create(para_env, pdims_3d, pgrid_3d)
1098 CALL create_3c_tensor(tensor, dist_ri, dist_ao_1, dist_ao_2, &
1099 pgrid_3d, sizes_ri, sizes_ao, sizes_ao, &
1100 map1=map1, map2=map2, name=tensor_name)
1101
1102 IF (PRESENT(create_nl_3c)) THEN
1103 my_create_nl_3c = create_nl_3c
1104 ELSE
1105 my_create_nl_3c = .false.
1106 END IF
1107
1108 IF (my_create_nl_3c) THEN
1109 CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set)
1110 CALL dbt_mp_environ_pgrid(pgrid_3d, pdims, pcoord)
1111 CALL mp_comm_t3c_2%create(pgrid_3d%mp_comm_2d, 3, pdims)
1112 CALL distribution_3d_create(dist_3d, dist_ri, dist_ao_1, dist_ao_2, &
1113 nkind, particle_set, mp_comm_t3c_2, own_comm=.true.)
1114
1115 CALL build_3c_neighbor_lists(nl_3c, &
1116 qs_env%bs_env%basis_set_RI, &
1117 qs_env%bs_env%basis_set_AO, &
1118 qs_env%bs_env%basis_set_AO, &
1119 dist_3d, qs_env%bs_env%ri_metric, &
1120 "GW_3c_nl", qs_env, own_dist=.true.)
1121 END IF
1122
1123 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
1124 CALL dbt_pgrid_destroy(pgrid_3d)
1125
1126 CALL timestop(handle)
1127
1128 END SUBROUTINE create_3c_t
1129
1130! **************************************************************************************************
1131!> \brief ...
1132!> \param bs_env ...
1133! **************************************************************************************************
1134 SUBROUTINE init_interaction_radii(bs_env)
1135 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1136
1137 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_interaction_radii'
1138
1139 INTEGER :: handle, ibasis
1140 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
1141
1142 CALL timeset(routinen, handle)
1143
1144 DO ibasis = 1, SIZE(bs_env%basis_set_AO)
1145
1146 orb_basis => bs_env%basis_set_AO(ibasis)%gto_basis_set
1147 CALL init_interaction_radii_orb_basis(orb_basis, bs_env%eps_filter)
1148
1149 ri_basis => bs_env%basis_set_RI(ibasis)%gto_basis_set
1150 CALL init_interaction_radii_orb_basis(ri_basis, bs_env%eps_filter)
1151
1152 END DO
1153
1154 CALL timestop(handle)
1155
1156 END SUBROUTINE init_interaction_radii
1157
1158! **************************************************************************************************
1159!> \brief ...
1160!> \param t_3c_int ...
1161!> \param max_dist_AO_atoms ...
1162!> \param qs_env ...
1163! **************************************************************************************************
1164 SUBROUTINE get_max_dist_ao_atoms(t_3c_int, max_dist_AO_atoms, qs_env)
1165 TYPE(dbt_type) :: t_3c_int
1166 REAL(kind=dp), INTENT(INOUT) :: max_dist_ao_atoms
1167 TYPE(qs_environment_type), POINTER :: qs_env
1168
1169 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_max_dist_AO_atoms'
1170
1171 INTEGER :: atom_1, atom_2, handle, num_cells
1172 INTEGER, DIMENSION(3) :: atom_ind
1173 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1174 REAL(kind=dp) :: abs_rab
1175 REAL(kind=dp), DIMENSION(3) :: rab
1176 TYPE(cell_type), POINTER :: cell
1177 TYPE(dbt_iterator_type) :: iter
1178 TYPE(mp_para_env_type), POINTER :: para_env
1179 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1180
1181 CALL timeset(routinen, handle)
1182
1183 NULLIFY (cell, particle_set, para_env)
1184 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, para_env=para_env)
1185
1186 ! max_dist_AO_atoms is compared to earlier steps in the loop with step n_atom_step
1187 ! do not initialize/overwrite here
1188
1189! IMPORTANT: Use thread-local copy for max_dist_AO_atoms via REDUCTION to avoid race conditions
1190!$OMP PARALLEL DEFAULT(NONE) &
1191!$OMP SHARED(t_3c_int, num_cells, index_to_cell, particle_set, cell) &
1192!$OMP PRIVATE(iter, atom_ind, rab, abs_rab, atom_1, atom_2) &
1193!$OMP REDUCTION(MAX:max_dist_AO_atoms)
1194
1195 CALL dbt_iterator_start(iter, t_3c_int)
1196 DO WHILE (dbt_iterator_blocks_left(iter))
1197 CALL dbt_iterator_next_block(iter, atom_ind)
1198
1199 atom_1 = atom_ind(2)
1200 atom_2 = atom_ind(3)
1201 rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
1202 abs_rab = sqrt(rab(1)**2 + rab(2)**2 + rab(3)**2)
1203
1204 ! Reduction takes care of using a thread-local copy
1205 max_dist_ao_atoms = max(max_dist_ao_atoms, abs_rab)
1206
1207 END DO
1208 CALL dbt_iterator_stop(iter)
1209!$OMP END PARALLEL
1210
1211 CALL para_env%max(max_dist_ao_atoms)
1212
1213 CALL timestop(handle)
1214
1215 END SUBROUTINE get_max_dist_ao_atoms
1216
1217! **************************************************************************************************
1218!> \brief ...
1219!> \param bs_env ...
1220! **************************************************************************************************
1221 SUBROUTINE set_sparsity_parallelization_parameters(bs_env)
1222 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1223
1224 CHARACTER(LEN=*), PARAMETER :: routinen = 'set_sparsity_parallelization_parameters'
1225
1226 INTEGER :: handle, i_ivl, il_ivl, j_ivl, n_atom_per_il_ivl, n_atom_per_ivl, n_intervals_i, &
1227 n_intervals_inner_loop_atoms, n_intervals_j, u
1228 INTEGER(KIND=int_8) :: input_memory_per_proc
1229
1230 CALL timeset(routinen, handle)
1231
1232 ! heuristic parameter to prevent out of memory
1233 bs_env%safety_factor_memory = 0.10_dp
1234
1235 input_memory_per_proc = int(bs_env%input_memory_per_proc_GB*1.0e9_dp, kind=int_8)
1236
1237 ! choose atomic range for λ ("i_atom"), ν ("j_atom") in
1238 ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
1239 ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
1240 ! such that M and N fit into the memory
1241 n_atom_per_ivl = int(sqrt(bs_env%safety_factor_memory*input_memory_per_proc &
1242 *bs_env%group_size_tensor/24/bs_env%n_RI &
1243 /sqrt(bs_env%occupation_3c_int)))/bs_env%max_AO_bf_per_atom
1244
1245 n_intervals_i = (bs_env%n_atom_i - 1)/n_atom_per_ivl + 1
1246 n_intervals_j = (bs_env%n_atom_j - 1)/n_atom_per_ivl + 1
1247
1248 bs_env%n_atom_per_interval_ij = n_atom_per_ivl
1249 bs_env%n_intervals_i = n_intervals_i
1250 bs_env%n_intervals_j = n_intervals_j
1251
1252 ALLOCATE (bs_env%i_atom_intervals(2, n_intervals_i))
1253 ALLOCATE (bs_env%j_atom_intervals(2, n_intervals_j))
1254
1255 DO i_ivl = 1, n_intervals_i
1256 bs_env%i_atom_intervals(1, i_ivl) = (i_ivl - 1)*n_atom_per_ivl + bs_env%atoms_i(1)
1257 bs_env%i_atom_intervals(2, i_ivl) = min(i_ivl*n_atom_per_ivl + bs_env%atoms_i(1) - 1, &
1258 bs_env%atoms_i(2))
1259 END DO
1260
1261 DO j_ivl = 1, n_intervals_j
1262 bs_env%j_atom_intervals(1, j_ivl) = (j_ivl - 1)*n_atom_per_ivl + bs_env%atoms_j(1)
1263 bs_env%j_atom_intervals(2, j_ivl) = min(j_ivl*n_atom_per_ivl + bs_env%atoms_j(1) - 1, &
1264 bs_env%atoms_j(2))
1265 END DO
1266
1267 ALLOCATE (bs_env%skip_Sigma_occ(n_intervals_i, n_intervals_j))
1268 ALLOCATE (bs_env%skip_Sigma_vir(n_intervals_i, n_intervals_j))
1269 bs_env%skip_Sigma_occ(:, :) = .false.
1270 bs_env%skip_Sigma_vir(:, :) = .false.
1271 bs_env%n_skip_chi = 0
1272
1273 ALLOCATE (bs_env%skip_chi(n_intervals_i, n_intervals_j))
1274 bs_env%skip_chi(:, :) = .false.
1275 bs_env%n_skip_sigma = 0
1276
1277 ! choose atomic range for µ and σ ("inner loop (IL) atom") in
1278 ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
1279 ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
1280 n_atom_per_il_ivl = min(int(bs_env%safety_factor_memory*input_memory_per_proc &
1281 *bs_env%group_size_tensor/n_atom_per_ivl &
1282 /bs_env%max_AO_bf_per_atom &
1283 /bs_env%n_RI/8/sqrt(bs_env%occupation_3c_int) &
1284 /bs_env%max_AO_bf_per_atom), bs_env%n_atom)
1285
1286 n_intervals_inner_loop_atoms = (bs_env%n_atom - 1)/n_atom_per_il_ivl + 1
1287
1288 bs_env%n_atom_per_IL_interval = n_atom_per_il_ivl
1289 bs_env%n_intervals_inner_loop_atoms = n_intervals_inner_loop_atoms
1290
1291 ALLOCATE (bs_env%inner_loop_atom_intervals(2, n_intervals_inner_loop_atoms))
1292 DO il_ivl = 1, n_intervals_inner_loop_atoms
1293 bs_env%inner_loop_atom_intervals(1, il_ivl) = (il_ivl - 1)*n_atom_per_il_ivl + 1
1294 bs_env%inner_loop_atom_intervals(2, il_ivl) = min(il_ivl*n_atom_per_il_ivl, bs_env%n_atom)
1295 END DO
1296
1297 u = bs_env%unit_nr
1298 IF (u > 0) THEN
1299 WRITE (u, '(T2,A)') ''
1300 WRITE (u, '(T2,A,I33)') λντνλτ'Number of i and j atoms in M_P(), N_Q():', n_atom_per_ivl
1301 WRITE (u, '(T2,A,I18)') µλνµµνµλ'Number of inner loop atoms for in M_P = sum_ (|P) G_', &
1302 n_atom_per_il_ivl
1303 END IF
1304
1305 CALL timestop(handle)
1306
1307 END SUBROUTINE set_sparsity_parallelization_parameters
1308
1309! **************************************************************************************************
1310!> \brief ...
1311!> \param qs_env ...
1312!> \param bs_env ...
1313! **************************************************************************************************
1314 SUBROUTINE check_for_restart_files(qs_env, bs_env)
1315 TYPE(qs_environment_type), POINTER :: qs_env
1316 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1317
1318 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_for_restart_files'
1319
1320 CHARACTER(LEN=9) :: frmt
1321 CHARACTER(len=default_path_length) :: f_chi, f_s_n, f_s_p, f_s_x, f_w_t, &
1322 prefix, project_name, z_lp_name
1323 INTEGER :: handle, i_spin, i_t_or_w, ind, n_spin, &
1324 num_time_freq_points
1325 LOGICAL :: chi_exists, sigma_neg_time_exists, &
1326 sigma_pos_time_exists, &
1327 sigma_x_spin_exists, w_time_exists, &
1328 z_lp_exists
1329 TYPE(cp_logger_type), POINTER :: logger
1330 TYPE(section_vals_type), POINTER :: input, print_key
1331
1332 CALL timeset(routinen, handle)
1333
1334 num_time_freq_points = bs_env%num_time_freq_points
1335 n_spin = bs_env%n_spin
1336
1337 ALLOCATE (bs_env%read_chi(num_time_freq_points))
1338 ALLOCATE (bs_env%calc_chi(num_time_freq_points))
1339 ALLOCATE (bs_env%Sigma_c_exists(num_time_freq_points, n_spin))
1340
1341 CALL get_qs_env(qs_env, input=input)
1342
1343 logger => cp_get_default_logger()
1344 print_key => section_vals_get_subs_vals(input, 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART')
1345 project_name = cp_print_key_generate_filename(logger, print_key, extension="", &
1346 my_local=.false.)
1347 WRITE (prefix, '(2A)') trim(project_name), "-RESTART_"
1348 bs_env%prefix = prefix
1349
1350 bs_env%all_W_exist = .true.
1351
1352 DO i_t_or_w = 1, num_time_freq_points
1353
1354 IF (i_t_or_w < 10) THEN
1355 WRITE (frmt, '(A)') '(3A,I1,A)'
1356 WRITE (f_chi, frmt) trim(prefix), bs_env%chi_name, "_0", i_t_or_w, ".matrix"
1357 WRITE (f_w_t, frmt) trim(prefix), bs_env%W_time_name, "_0", i_t_or_w, ".matrix"
1358 ELSE IF (i_t_or_w < 100) THEN
1359 WRITE (frmt, '(A)') '(3A,I2,A)'
1360 WRITE (f_chi, frmt) trim(prefix), bs_env%chi_name, "_", i_t_or_w, ".matrix"
1361 WRITE (f_w_t, frmt) trim(prefix), bs_env%W_time_name, "_", i_t_or_w, ".matrix"
1362 ELSE
1363 cpabort('Please implement more than 99 time/frequency points.')
1364 END IF
1365
1366 INQUIRE (file=trim(f_chi), exist=chi_exists)
1367 INQUIRE (file=trim(f_w_t), exist=w_time_exists)
1368
1369 bs_env%read_chi(i_t_or_w) = chi_exists
1370 bs_env%calc_chi(i_t_or_w) = .NOT. chi_exists
1371
1372 bs_env%all_W_exist = bs_env%all_W_exist .AND. w_time_exists
1373
1374 ! the self-energy is spin-dependent
1375 DO i_spin = 1, n_spin
1376
1377 ind = i_t_or_w + (i_spin - 1)*num_time_freq_points
1378
1379 IF (ind < 10) THEN
1380 WRITE (frmt, '(A)') '(3A,I1,A)'
1381 WRITE (f_s_p, frmt) trim(prefix), bs_env%Sigma_p_name, "_0", ind, ".matrix"
1382 WRITE (f_s_n, frmt) trim(prefix), bs_env%Sigma_n_name, "_0", ind, ".matrix"
1383 ELSE IF (ind < 100) THEN
1384 WRITE (frmt, '(A)') '(3A,I2,A)'
1385 WRITE (f_s_p, frmt) trim(prefix), bs_env%Sigma_p_name, "_", ind, ".matrix"
1386 WRITE (f_s_n, frmt) trim(prefix), bs_env%Sigma_n_name, "_", ind, ".matrix"
1387 ELSE
1388 cpabort('Please implement more than 99 combined spin+freq indices.')
1389 END IF
1390
1391 INQUIRE (file=trim(f_s_p), exist=sigma_pos_time_exists)
1392 INQUIRE (file=trim(f_s_n), exist=sigma_neg_time_exists)
1393
1394 bs_env%Sigma_c_exists(i_t_or_w, i_spin) = sigma_pos_time_exists .AND. &
1395 sigma_neg_time_exists
1396
1397 END DO
1398
1399 END DO
1400
1401 ! Marek : In the RTBSE run, check also for zero frequency W
1402 IF (bs_env%rtp_method == rtp_method_bse .OR. &
1403 bs_env%rtp_method == rtp_method_bse_linearized) THEN
1404 WRITE (f_w_t, '(3A,I1,A)') trim(prefix), "W_freq_rtp", "_0", 0, ".matrix"
1405 INQUIRE (file=trim(f_w_t), exist=w_time_exists)
1406 bs_env%all_W_exist = bs_env%all_W_exist .AND. w_time_exists
1407 END IF
1408
1409 ! Check for Restart Z_lP file
1410 IF (bs_env%do_gw_ri_rs) THEN
1411 WRITE (z_lp_name, '(3A)') trim(prefix), "Z_lP", ".matrix"
1412 INQUIRE (file=trim(z_lp_name), exist=z_lp_exists)
1413 bs_env%ri_rs%Z_lP_exists = z_lp_exists
1414 END IF
1415
1416 IF (bs_env%all_W_exist) THEN
1417 bs_env%read_chi(:) = .false.
1418 bs_env%calc_chi(:) = .false.
1419 END IF
1420
1421 bs_env%Sigma_x_exists = .true.
1422 DO i_spin = 1, n_spin
1423 WRITE (f_s_x, '(3A,I1,A)') trim(prefix), bs_env%Sigma_x_name, "_0", i_spin, ".matrix"
1424 INQUIRE (file=trim(f_s_x), exist=sigma_x_spin_exists)
1425 bs_env%Sigma_x_exists = bs_env%Sigma_x_exists .AND. sigma_x_spin_exists
1426 END DO
1427
1428 ! If any restart files are read, check if the SCF converged in 1 step.
1429 ! This is important because a re-iterated SCF can lead to spurious GW results
1430 IF (any(bs_env%read_chi(:)) &
1431 .OR. any(bs_env%Sigma_c_exists) &
1432 .OR. bs_env%all_W_exist &
1433 .OR. bs_env%Sigma_x_exists &
1434 ) THEN
1435
1436 IF (qs_env%scf_env%iter_count /= 1) THEN
1437 CALL cp_warn(__location__, "SCF needed more than 1 step, "// &
1438 "which might lead to spurious GW results when using GW restart files. ")
1439 END IF
1440 END IF
1441
1442 CALL timestop(handle)
1443
1444 END SUBROUTINE check_for_restart_files
1445
1446! **************************************************************************************************
1447!> \brief ...
1448!> \param qs_env ...
1449!> \param bs_env ...
1450! **************************************************************************************************
1451 SUBROUTINE set_parallelization_parameters(qs_env, bs_env)
1452 TYPE(qs_environment_type), POINTER :: qs_env
1453 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1454
1455 CHARACTER(LEN=*), PARAMETER :: routinen = 'set_parallelization_parameters'
1456
1457 INTEGER :: color_sub, dummy_1, dummy_2, handle, &
1458 num_pe, num_t_groups, u
1459 TYPE(mp_para_env_type), POINTER :: para_env
1460
1461 CALL timeset(routinen, handle)
1462
1463 CALL get_qs_env(qs_env, para_env=para_env)
1464
1465 num_pe = para_env%num_pe
1466 ! if not already set, use all processors for the group (for large-cell GW, performance
1467 ! seems to be best for a single group with all MPI processes per group)
1468 IF (bs_env%group_size_tensor < 0 .OR. bs_env%group_size_tensor > num_pe) THEN
1469 bs_env%group_size_tensor = num_pe
1470 END IF
1471
1472 ! group_size_tensor must divide num_pe without rest; otherwise everything will be complicated
1473 IF (modulo(num_pe, bs_env%group_size_tensor) /= 0) THEN
1474 CALL find_good_group_size(num_pe, bs_env%group_size_tensor)
1475 END IF
1476
1477 ! para_env_tensor for tensor subgroups
1478 color_sub = para_env%mepos/bs_env%group_size_tensor
1479 bs_env%tensor_group_color = color_sub
1480
1481 ALLOCATE (bs_env%para_env_tensor)
1482 CALL bs_env%para_env_tensor%from_split(para_env, color_sub)
1483
1484 num_t_groups = para_env%num_pe/bs_env%group_size_tensor
1485 bs_env%num_tensor_groups = num_t_groups
1486
1487 CALL get_i_j_atoms(bs_env%atoms_i, bs_env%atoms_j, bs_env%n_atom_i, bs_env%n_atom_j, &
1488 color_sub, bs_env)
1489
1490 ALLOCATE (bs_env%atoms_i_t_group(2, num_t_groups))
1491 ALLOCATE (bs_env%atoms_j_t_group(2, num_t_groups))
1492 DO color_sub = 0, num_t_groups - 1
1493 CALL get_i_j_atoms(bs_env%atoms_i_t_group(1:2, color_sub + 1), &
1494 bs_env%atoms_j_t_group(1:2, color_sub + 1), &
1495 dummy_1, dummy_2, color_sub, bs_env)
1496 END DO
1497
1498 u = bs_env%unit_nr
1499 IF (u > 0) THEN
1500 WRITE (u, '(T2,A,I47)') 'Group size for tensor operations', bs_env%group_size_tensor
1501 IF (bs_env%group_size_tensor > 1 .AND. bs_env%n_atom < 5) THEN
1502 WRITE (u, '(T2,A)') 'The requested group size is > 1 which can lead to bad performance.'
1503 WRITE (u, '(T2,A)') 'Using more memory per MPI process might improve performance.'
1504 WRITE (u, '(T2,A)') '(Also increase MEMORY_PER_PROC when using more memory per process.)'
1505 END IF
1506 END IF
1507
1508 CALL timestop(handle)
1509
1510 END SUBROUTINE set_parallelization_parameters
1511
1512! **************************************************************************************************
1513!> \brief ...
1514!> \param num_pe ...
1515!> \param group_size ...
1516! **************************************************************************************************
1517 SUBROUTINE find_good_group_size(num_pe, group_size)
1518
1519 INTEGER :: num_pe, group_size
1520
1521 CHARACTER(LEN=*), PARAMETER :: routinen = 'find_good_group_size'
1522
1523 INTEGER :: group_size_minus, group_size_orig, &
1524 group_size_plus, handle, i_diff
1525
1526 CALL timeset(routinen, handle)
1527
1528 group_size_orig = group_size
1529
1530 DO i_diff = 1, num_pe
1531
1532 group_size_minus = group_size - i_diff
1533
1534 IF (modulo(num_pe, group_size_minus) == 0 .AND. group_size_minus > 0) THEN
1535 group_size = group_size_minus
1536 EXIT
1537 END IF
1538
1539 group_size_plus = group_size + i_diff
1540
1541 IF (modulo(num_pe, group_size_plus) == 0 .AND. group_size_plus <= num_pe) THEN
1542 group_size = group_size_plus
1543 EXIT
1544 END IF
1545
1546 END DO
1547
1548 IF (group_size_orig == group_size) cpabort("Group size error")
1549
1550 CALL timestop(handle)
1551
1552 END SUBROUTINE find_good_group_size
1553
1554! **************************************************************************************************
1555!> \brief ...
1556!> \param atoms_i ...
1557!> \param atoms_j ...
1558!> \param n_atom_i ...
1559!> \param n_atom_j ...
1560!> \param color_sub ...
1561!> \param bs_env ...
1562! **************************************************************************************************
1563 SUBROUTINE get_i_j_atoms(atoms_i, atoms_j, n_atom_i, n_atom_j, color_sub, bs_env)
1564
1565 INTEGER, DIMENSION(2) :: atoms_i, atoms_j
1566 INTEGER :: n_atom_i, n_atom_j, color_sub
1567 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1568
1569 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_i_j_atoms'
1570
1571 INTEGER :: handle, i_atoms_per_group, i_group, &
1572 ipcol, ipcol_loop, iprow, iprow_loop, &
1573 j_atoms_per_group, npcol, nprow
1574
1575 CALL timeset(routinen, handle)
1576
1577 ! create a square mesh of tensor groups for iatom and jatom; code from blacs_env_create
1578 CALL square_mesh(nprow, npcol, bs_env%num_tensor_groups)
1579
1580 i_group = 0
1581 DO ipcol_loop = 0, npcol - 1
1582 DO iprow_loop = 0, nprow - 1
1583 IF (i_group == color_sub) THEN
1584 iprow = iprow_loop
1585 ipcol = ipcol_loop
1586 END IF
1587 i_group = i_group + 1
1588 END DO
1589 END DO
1590
1591 IF (modulo(bs_env%n_atom, nprow) == 0) THEN
1592 i_atoms_per_group = bs_env%n_atom/nprow
1593 ELSE
1594 i_atoms_per_group = bs_env%n_atom/nprow + 1
1595 END IF
1596
1597 IF (modulo(bs_env%n_atom, npcol) == 0) THEN
1598 j_atoms_per_group = bs_env%n_atom/npcol
1599 ELSE
1600 j_atoms_per_group = bs_env%n_atom/npcol + 1
1601 END IF
1602
1603 atoms_i(1) = iprow*i_atoms_per_group + 1
1604 atoms_i(2) = min((iprow + 1)*i_atoms_per_group, bs_env%n_atom)
1605 n_atom_i = atoms_i(2) - atoms_i(1) + 1
1606
1607 atoms_j(1) = ipcol*j_atoms_per_group + 1
1608 atoms_j(2) = min((ipcol + 1)*j_atoms_per_group, bs_env%n_atom)
1609 n_atom_j = atoms_j(2) - atoms_j(1) + 1
1610
1611 CALL timestop(handle)
1612
1613 END SUBROUTINE get_i_j_atoms
1614
1615! **************************************************************************************************
1616!> \brief ...
1617!> \param nprow ...
1618!> \param npcol ...
1619!> \param nproc ...
1620! **************************************************************************************************
1621 SUBROUTINE square_mesh(nprow, npcol, nproc)
1622 INTEGER :: nprow, npcol, nproc
1623
1624 CHARACTER(LEN=*), PARAMETER :: routinen = 'square_mesh'
1625
1626 INTEGER :: gcd_max, handle, ipe, jpe
1627
1628 CALL timeset(routinen, handle)
1629
1630 gcd_max = -1
1631 DO ipe = 1, ceiling(sqrt(real(nproc, dp)))
1632 jpe = nproc/ipe
1633 IF (ipe*jpe /= nproc) cycle
1634 IF (gcd(ipe, jpe) >= gcd_max) THEN
1635 nprow = ipe
1636 npcol = jpe
1637 gcd_max = gcd(ipe, jpe)
1638 END IF
1639 END DO
1640
1641 CALL timestop(handle)
1642
1643 END SUBROUTINE square_mesh
1644
1645! **************************************************************************************************
1646!> \brief ...
1647!> \param bs_env ...
1648!> \param qs_env ...
1649! **************************************************************************************************
1650 SUBROUTINE set_heuristic_parameters(bs_env, qs_env)
1651 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1652 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
1653
1654 CHARACTER(LEN=*), PARAMETER :: routinen = 'set_heuristic_parameters'
1655
1656 INTEGER :: handle, u
1657 LOGICAL :: do_bvk_cell
1658
1659 CALL timeset(routinen, handle)
1660
1661 ! for generating numerically stable minimax Fourier integration weights
1662 bs_env%num_points_per_magnitude = 200
1663
1664 IF (bs_env%input_regularization_minimax > -1.0e-12_dp) THEN
1665 bs_env%regularization_minimax = bs_env%input_regularization_minimax
1666 ELSE
1667 ! for periodic systems and for 20 minimax points, we use a regularized minimax mesh
1668 ! (from experience: regularized minimax meshes converges faster for periodic systems
1669 ! and for 20 pts)
1670 IF (sum(bs_env%periodic) /= 0 .OR. bs_env%num_time_freq_points >= 20) THEN
1671 bs_env%regularization_minimax = 1.0e-6_dp
1672 ELSE
1673 bs_env%regularization_minimax = 0.0_dp
1674 END IF
1675 END IF
1676
1677 bs_env%stabilize_exp = 70.0_dp
1678 bs_env%eps_atom_grid_2d_mat = 1.0e-50_dp
1679
1680 ! use a 16-parameter Padé fit
1681 bs_env%nparam_pade = 16
1682
1683 ! resolution of the identity with the truncated Coulomb metric, cutoff radius 3 Angström
1684 bs_env%ri_metric%potential_type = do_potential_truncated
1685 bs_env%ri_metric%omega = 0.0_dp
1686 ! cutoff radius is specified in the input
1687 bs_env%ri_metric%filename = "t_c_g.dat"
1688
1689 bs_env%eps_eigval_mat_RI = 0.0_dp
1690
1691 IF (bs_env%input_regularization_RI > -1.0e-12_dp) THEN
1692 bs_env%regularization_RI = bs_env%input_regularization_RI
1693 ELSE
1694 ! default case:
1695
1696 ! 1. for periodic systems, we use the regularized resolution of the identity per default
1697 bs_env%regularization_RI = 1.0e-2_dp
1698
1699 ! 2. for molecules, no regularization is necessary
1700 IF (sum(bs_env%periodic) == 0) bs_env%regularization_RI = 0.0_dp
1701
1702 END IF
1703
1704 ! truncated Coulomb operator for exchange self-energy
1705 ! (see details in Guidon, VandeVondele, Hutter, JCTC 5, 3010 (2009) and references therein)
1706 do_bvk_cell = bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp
1707 CALL trunc_coulomb_for_exchange(qs_env, bs_env%trunc_coulomb, &
1708 rel_cutoff_trunc_coulomb_ri_x=0.5_dp, &
1709 cell_grid=bs_env%cell_grid_scf_desymm, &
1710 do_bvk_cell=do_bvk_cell)
1711
1712 ! for small-cell GW, we need more cells than normally used by the filter bs_env%eps_filter
1713 ! (in particular for computing the self-energy because of higher number of cells needed)
1714 bs_env%heuristic_filter_factor = 1.0e-4
1715
1716 u = bs_env%unit_nr
1717 IF (u > 0) THEN
1718 WRITE (u, fmt="(T2,2A,F21.1,A)") "Cutoff radius for the truncated Coulomb ", &
1719 Σ"operator in ^x:", bs_env%trunc_coulomb%cutoff_radius*angstrom, Å" "
1720 WRITE (u, fmt="(T2,2A,F15.1,A)") "Cutoff radius for the truncated Coulomb ", &
1721 "operator in RI metric:", bs_env%ri_metric%cutoff_radius*angstrom, Å" "
1722 WRITE (u, fmt="(T2,A,ES48.1)") "Regularization parameter of RI ", bs_env%regularization_RI
1723 WRITE (u, fmt="(T2,A,ES38.1)") "Regularization parameter of minimax grids", &
1724 bs_env%regularization_minimax
1725 WRITE (u, fmt="(T2,A,I53)") "Lattice sum size for V(k):", bs_env%size_lattice_sum_V
1726 END IF
1727
1728 CALL timestop(handle)
1729
1730 END SUBROUTINE set_heuristic_parameters
1731
1732! **************************************************************************************************
1733!> \brief ...
1734!> \param bs_env ...
1735! **************************************************************************************************
1736 SUBROUTINE print_header_and_input_parameters(bs_env)
1737
1738 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1739
1740 CHARACTER(LEN=*), PARAMETER :: routinen = 'print_header_and_input_parameters'
1741
1742 INTEGER :: handle, u
1743
1744 CALL timeset(routinen, handle)
1745
1746 u = bs_env%unit_nr
1747
1748 IF (u > 0) THEN
1749 WRITE (u, '(T2,A)') ' '
1750 WRITE (u, '(T2,A)') repeat('-', 79)
1751 WRITE (u, '(T2,A,A78)') '-', '-'
1752 WRITE (u, '(T2,A,A46,A32)') '-', 'GW CALCULATION', '-'
1753 WRITE (u, '(T2,A,A78)') '-', '-'
1754 WRITE (u, '(T2,A)') repeat('-', 79)
1755 WRITE (u, '(T2,A)') ' '
1756 WRITE (u, '(T2,A,I45)') 'Input: Number of time/freq. points', bs_env%num_time_freq_points
1757 WRITE (u, "(T2,A,F44.1,A)") ωΣω'Input: _max for fitting (i) (eV)', bs_env%freq_max_fit*evolt
1758 WRITE (u, '(T2,A,ES27.1)') 'Input: Filter threshold for sparse tensor operations', &
1759 bs_env%eps_filter
1760 WRITE (u, "(T2,A,L55)") 'Input: Apply Hedin shift', bs_env%do_hedin_shift
1761 WRITE (u, '(T2,A,F37.1,A)') 'Input: Available memory per MPI process', &
1762 bs_env%input_memory_per_proc_GB, ' GB'
1763 IF (bs_env%do_gw_ri_rs) THEN
1764 WRITE (u, '(A)') ' '
1765 WRITE (u, '(T2,A,ES43.2)') 'Input: RI-RS Tikhonov regularization', &
1766 bs_env%ri_rs%tikhonov
1767 IF (bs_env%ri_rs%cutoff_radius_ri_rs > 0.0_dp) THEN
1768 WRITE (u, '(T2,A,F39.2,A)') 'Input: RI-RS integration sphere cutoff', &
1769 bs_env%ri_rs%cutoff_radius_ri_rs*angstrom, Å' '
1770 END IF
1771 IF (bs_env%ri_rs%cutoff_radius_ri_ao > 0.0_dp) THEN
1772 WRITE (u, '(T2,A,F44.2,A)') 'Input: AO grid hard cutoff radius', &
1773 bs_env%ri_rs%cutoff_radius_ri_ao*angstrom, Å' '
1774 END IF
1775 WRITE (u, '(T2,A,I40)') 'Input: MPI ranks per atom in Z_lP solve', &
1776 bs_env%ri_rs%n_procs_per_atom_z_lp
1777 WRITE (u, '(T2,A,L43)') χ'Input: Keep sparsity in /G/W panels', &
1778 bs_env%ri_rs%keep_sparsity_rirs
1779 IF (bs_env%ri_rs%cutoff_radius_v_w > 0.0_dp) THEN
1780 WRITE (u, '(T2,A,F43.2,A)') 'Input: G/W panel truncation radius', &
1781 bs_env%ri_rs%cutoff_radius_v_w*angstrom, Å' '
1782 END IF
1783 IF (bs_env%ri_rs%cutoff_radius_g_w > 0.0_dp) THEN
1784 WRITE (u, '(T2,A,F40.2,A)') 'Input: G/W operator truncation radius', &
1785 bs_env%ri_rs%cutoff_radius_g_w*angstrom, Å' '
1786 END IF
1787 WRITE (u, '(A)') ' '
1788 END IF
1789 END IF
1790
1791 CALL timestop(handle)
1792
1793 END SUBROUTINE print_header_and_input_parameters
1794
1795! **************************************************************************************************
1796!> \brief ...
1797!> \param qs_env ...
1798!> \param bs_env ...
1799! **************************************************************************************************
1800 SUBROUTINE compute_v_xc(qs_env, bs_env)
1801 TYPE(qs_environment_type), POINTER :: qs_env
1802 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1803
1804 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_V_xc'
1805
1806 INTEGER :: handle, img, ispin, myfun, nimages
1807 LOGICAL :: hf_present
1808 REAL(kind=dp) :: energy_ex, energy_exc, energy_total, &
1809 myfraction
1810 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_ks_without_v_xc
1811 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
1812 TYPE(dft_control_type), POINTER :: dft_control
1813 TYPE(qs_energy_type), POINTER :: energy
1814 TYPE(section_vals_type), POINTER :: hf_section, input, xc_section
1815
1816 CALL timeset(routinen, handle)
1817
1818 CALL get_qs_env(qs_env, input=input, energy=energy, dft_control=dft_control)
1819
1820 ! previously, dft_control%nimages set to # neighbor cells, revert for Γ-only KS matrix
1821 nimages = dft_control%nimages
1822 dft_control%nimages = bs_env%nimages_scf
1823
1824 ! we need to reset XC functional, therefore, get XC input
1825 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1826 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
1827 CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=xc_none)
1828 ! IF (ASSOCIATED(section_vals_get_subs_vals(xc_section, "HF", can_return_null=.TRUE.))) THEN
1829 hf_section => section_vals_get_subs_vals(input, "DFT%XC%HF", can_return_null=.true.)
1830 hf_present = .false.
1831 IF (ASSOCIATED(hf_section)) THEN
1832 CALL section_vals_get(hf_section, explicit=hf_present)
1833 END IF
1834 IF (hf_present) THEN
1835 ! Special case for handling hfx
1836 CALL section_vals_val_get(xc_section, "HF%FRACTION", r_val=myfraction)
1837 CALL section_vals_val_set(xc_section, "HF%FRACTION", r_val=0.0_dp)
1838 END IF
1839
1840 ! save the energy before the energy gets updated
1841 energy_total = energy%total
1842 energy_exc = energy%exc
1843 energy_ex = energy%ex
1844
1845 SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
1846 CASE (large_cell_gamma, large_cell_gamma_ri_rs, non_periodic_ri_rs)
1847
1848 NULLIFY (mat_ks_without_v_xc)
1849 CALL dbcsr_allocate_matrix_set(mat_ks_without_v_xc, bs_env%n_spin)
1850
1851 DO ispin = 1, bs_env%n_spin
1852 ALLOCATE (mat_ks_without_v_xc(ispin)%matrix)
1853 IF (hf_present) THEN
1854 CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix, &
1855 matrix_type=dbcsr_type_symmetric)
1856 ELSE
1857 CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1858 END IF
1859 END DO
1860
1861 ! calculate KS-matrix without XC
1862 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., just_energy=.false., &
1863 ext_ks_matrix=mat_ks_without_v_xc)
1864
1865 DO ispin = 1, bs_env%n_spin
1866 ! transfer dbcsr matrix to fm
1867 CALL cp_fm_create(bs_env%fm_V_xc_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1868 CALL copy_dbcsr_to_fm(mat_ks_without_v_xc(ispin)%matrix, bs_env%fm_V_xc_Gamma(ispin))
1869
1870 ! v_xc = h_ks - h_ks(v_xc = 0)
1871 CALL cp_fm_scale_and_add(alpha=-1.0_dp, matrix_a=bs_env%fm_V_xc_Gamma(ispin), &
1872 beta=1.0_dp, matrix_b=bs_env%fm_ks_Gamma(ispin))
1873 END DO
1874
1875 CALL dbcsr_deallocate_matrix_set(mat_ks_without_v_xc)
1876
1877 CASE (small_cell_full_kp)
1878
1879 ! calculate KS-matrix without XC
1880 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., just_energy=.false.)
1881 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_kp)
1882
1883 ALLOCATE (bs_env%fm_V_xc_R(dft_control%nimages, bs_env%n_spin))
1884 DO ispin = 1, bs_env%n_spin
1885 DO img = 1, dft_control%nimages
1886 ! safe fm_V_xc_R in fm_matrix because saving in dbcsr matrix caused trouble...
1887 CALL copy_dbcsr_to_fm(matrix_ks_kp(ispin, img)%matrix, bs_env%fm_work_mo(1))
1888 CALL cp_fm_create(bs_env%fm_V_xc_R(img, ispin), bs_env%fm_work_mo(1)%matrix_struct, &
1889 set_zero=.true.)
1890 ! store h_ks(v_xc = 0) in fm_V_xc_R
1891 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=bs_env%fm_V_xc_R(img, ispin), &
1892 beta=1.0_dp, matrix_b=bs_env%fm_work_mo(1))
1893 END DO
1894 END DO
1895
1896 END SELECT
1897
1898 ! set back the energy
1899 energy%total = energy_total
1900 energy%exc = energy_exc
1901 energy%ex = energy_ex
1902
1903 ! set back nimages
1904 dft_control%nimages = nimages
1905
1906 ! set the DFT functional and HF fraction back
1907 CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
1908 i_val=myfun)
1909 IF (hf_present) THEN
1910 CALL section_vals_val_set(xc_section, "HF%FRACTION", &
1911 r_val=myfraction)
1912 END IF
1913
1914 IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
1915 ! calculate KS-matrix again with XC
1916 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., just_energy=.false.)
1917 DO ispin = 1, bs_env%n_spin
1918 DO img = 1, dft_control%nimages
1919 ! store h_ks in fm_work_mo
1920 CALL copy_dbcsr_to_fm(matrix_ks_kp(ispin, img)%matrix, bs_env%fm_work_mo(1))
1921 ! v_xc = h_ks - h_ks(v_xc = 0)
1922 CALL cp_fm_scale_and_add(alpha=-1.0_dp, matrix_a=bs_env%fm_V_xc_R(img, ispin), &
1923 beta=1.0_dp, matrix_b=bs_env%fm_work_mo(1))
1924 END DO
1925 END DO
1926 END IF
1927
1928 CALL timestop(handle)
1929
1930 END SUBROUTINE compute_v_xc
1931
1932! **************************************************************************************************
1933!> \brief ...
1934!> \param bs_env ...
1935! **************************************************************************************************
1936 SUBROUTINE setup_time_and_frequency_minimax_grid(bs_env)
1937 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1938
1939 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_time_and_frequency_minimax_grid'
1940
1941 INTEGER :: handle, homo, i_w, ierr, ispin, j_w, &
1942 n_mo, num_time_freq_points, u
1943 REAL(kind=dp) :: e_max, e_max_ispin, e_min, e_min_ispin, &
1944 e_range, max_error_min
1945 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: points_and_weights
1946
1947 CALL timeset(routinen, handle)
1948
1949 n_mo = bs_env%n_ao
1950 num_time_freq_points = bs_env%num_time_freq_points
1951
1952 ALLOCATE (bs_env%imag_freq_points(num_time_freq_points))
1953 ALLOCATE (bs_env%imag_time_points(num_time_freq_points))
1954 ALLOCATE (bs_env%imag_time_weights_freq_zero(num_time_freq_points))
1955 ALLOCATE (bs_env%weights_cos_t_to_w(num_time_freq_points, num_time_freq_points))
1956 ALLOCATE (bs_env%weights_cos_w_to_t(num_time_freq_points, num_time_freq_points))
1957 ALLOCATE (bs_env%weights_sin_t_to_w(num_time_freq_points, num_time_freq_points))
1958
1959 ! minimum and maximum difference between eigenvalues of unoccupied and an occupied MOs
1960 e_min = 1000.0_dp
1961 e_max = -1000.0_dp
1962 DO ispin = 1, bs_env%n_spin
1963 homo = bs_env%n_occ(ispin)
1964 SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
1965 CASE (large_cell_gamma, large_cell_gamma_ri_rs, non_periodic_ri_rs)
1966 e_min_ispin = bs_env%eigenval_scf_Gamma(homo + 1, ispin) - &
1967 bs_env%eigenval_scf_Gamma(homo, ispin)
1968 e_max_ispin = bs_env%eigenval_scf_Gamma(n_mo, ispin) - &
1969 bs_env%eigenval_scf_Gamma(1, ispin)
1970 CASE (small_cell_full_kp)
1971 e_min_ispin = minval(bs_env%eigenval_scf(homo + 1, :, ispin)) - &
1972 maxval(bs_env%eigenval_scf(homo, :, ispin))
1973 e_max_ispin = maxval(bs_env%eigenval_scf(n_mo, :, ispin)) - &
1974 minval(bs_env%eigenval_scf(1, :, ispin))
1975 END SELECT
1976 e_min = min(e_min, e_min_ispin)
1977 e_max = max(e_max, e_max_ispin)
1978 END DO
1979
1980 ! Open-shell uses ONE minimax grid for the combined [min gap, max span] over both spins (the
1981 ! superset covers each channel, so it is accurate; per-spin grids would only be more efficient).
1982 IF (bs_env%n_spin > 1) THEN
1983 CALL cp_hint(__location__, &
1984 "Open-shell GW uses one minimax grid spanning [min gap, max span] across both "// &
1985 "spin channels; raise NUM_TIME_FREQ_POINTS if QP convergence is marginal for "// &
1986 "strongly spin-asymmetric systems.")
1987 END IF
1988
1989 e_range = e_max/e_min
1990
1991 ALLOCATE (points_and_weights(2*num_time_freq_points))
1992
1993 ! frequency points
1994 IF (num_time_freq_points <= 20) THEN
1995 CALL get_rpa_minimax_coeff(num_time_freq_points, e_range, points_and_weights, ierr, .false.)
1996 ELSE
1997 CALL get_rpa_minimax_coeff_larger_grid(num_time_freq_points, e_range, points_and_weights)
1998 END IF
1999
2000 ! one needs to scale the minimax grids, see Azizi, Wilhelm, Golze, Panades-Barrueta,
2001 ! Giantomassi, Rinke, Draxl, Gonze et al., 2 publications
2002 bs_env%imag_freq_points(:) = points_and_weights(1:num_time_freq_points)*e_min
2003
2004 ! determine number of fit points in the interval [0,ω_max] for virt, or [-ω_max,0] for occ
2005 bs_env%num_freq_points_fit = 0
2006 DO i_w = 1, num_time_freq_points
2007 IF (bs_env%imag_freq_points(i_w) < bs_env%freq_max_fit) THEN
2008 bs_env%num_freq_points_fit = bs_env%num_freq_points_fit + 1
2009 END IF
2010 END DO
2011
2012 ! iω values for the analytic continuation Σ^c_n(iω,k) -> Σ^c_n(ϵ,k)
2013 ALLOCATE (bs_env%imag_freq_points_fit(bs_env%num_freq_points_fit))
2014 j_w = 0
2015 DO i_w = 1, num_time_freq_points
2016 IF (bs_env%imag_freq_points(i_w) < bs_env%freq_max_fit) THEN
2017 j_w = j_w + 1
2018 bs_env%imag_freq_points_fit(j_w) = bs_env%imag_freq_points(i_w)
2019 END IF
2020 END DO
2021
2022 ! reset the number of Padé parameters if smaller than the number of
2023 ! imaginary-frequency points for the fit
2024 IF (bs_env%num_freq_points_fit < bs_env%nparam_pade) THEN
2025 bs_env%nparam_pade = bs_env%num_freq_points_fit
2026 END IF
2027
2028 ! time points
2029 IF (num_time_freq_points <= 20) THEN
2030 CALL get_exp_minimax_coeff(num_time_freq_points, e_range, points_and_weights)
2031 ELSE
2032 CALL get_exp_minimax_coeff_gw(num_time_freq_points, e_range, points_and_weights)
2033 END IF
2034
2035 bs_env%imag_time_points(:) = points_and_weights(1:num_time_freq_points)/(2.0_dp*e_min)
2036 bs_env%imag_time_weights_freq_zero(:) = points_and_weights(num_time_freq_points + 1:)/(e_min)
2037
2038 DEALLOCATE (points_and_weights)
2039
2040 u = bs_env%unit_nr
2041 IF (u > 0) THEN
2042 WRITE (u, '(T2,A)') ''
2043 WRITE (u, '(T2,A,F55.2)') 'SCF direct band gap (eV)', e_min*evolt
2044 WRITE (u, '(T2,A,F53.2)') 'Max. SCF eigval diff. (eV)', e_max*evolt
2045 WRITE (u, '(T2,A,F55.2)') 'E-Range for minimax grid', e_range
2046 WRITE (u, '(T2,A,I27)') é'Number of Pad parameters for analytic continuation:', &
2047 bs_env%nparam_pade
2048 WRITE (u, '(T2,A)') ''
2049 END IF
2050
2051 ! in minimax grids, Fourier transforms t -> w and w -> t are split using
2052 ! e^(iwt) = cos(wt) + i sin(wt); we thus calculate weights for trafos with a cos and
2053 ! sine prefactor; details in Azizi, Wilhelm, Golze, Giantomassi, Panades-Barrueta,
2054 ! Rinke, Draxl, Gonze et al., 2 publications
2055
2056 ! cosine transform weights imaginary time to imaginary frequency
2057 CALL get_l_sq_wghts_cos_tf_t_to_w(num_time_freq_points, &
2058 bs_env%imag_time_points, &
2059 bs_env%weights_cos_t_to_w, &
2060 bs_env%imag_freq_points, &
2061 e_min, e_max, max_error_min, &
2062 bs_env%num_points_per_magnitude, &
2063 bs_env%regularization_minimax)
2064
2065 ! cosine transform weights imaginary frequency to imaginary time
2066 CALL get_l_sq_wghts_cos_tf_w_to_t(num_time_freq_points, &
2067 bs_env%imag_time_points, &
2068 bs_env%weights_cos_w_to_t, &
2069 bs_env%imag_freq_points, &
2070 e_min, e_max, max_error_min, &
2071 bs_env%num_points_per_magnitude, &
2072 bs_env%regularization_minimax)
2073
2074 ! sine transform weights imaginary time to imaginary frequency
2075 CALL get_l_sq_wghts_sin_tf_t_to_w(num_time_freq_points, &
2076 bs_env%imag_time_points, &
2077 bs_env%weights_sin_t_to_w, &
2078 bs_env%imag_freq_points, &
2079 e_min, e_max, max_error_min, &
2080 bs_env%num_points_per_magnitude, &
2081 bs_env%regularization_minimax)
2082
2083 CALL timestop(handle)
2084
2085 END SUBROUTINE setup_time_and_frequency_minimax_grid
2086
2087! **************************************************************************************************
2088!> \brief ...
2089!> \param qs_env ...
2090!> \param bs_env ...
2091! **************************************************************************************************
2092 SUBROUTINE setup_cells_3c(qs_env, bs_env)
2093
2094 TYPE(qs_environment_type), POINTER :: qs_env
2095 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2096
2097 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_cells_3c'
2098
2099 INTEGER :: atom_i, atom_j, atom_k, block_count, handle, i, i_cell_x, i_cell_x_max, &
2100 i_cell_x_min, i_size, ikind, img, j, j_cell, j_cell_max, j_cell_y, j_cell_y_max, &
2101 j_cell_y_min, j_size, k_cell, k_cell_max, k_cell_z, k_cell_z_max, k_cell_z_min, k_size, &
2102 nimage_pairs_3c, nimages_3c, nimages_3c_max, nkind, u
2103 INTEGER(KIND=int_8) :: mem_occ_per_proc
2104 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, n_other_3c_images_max
2105 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c_max, nblocks_3c_max
2106 INTEGER, DIMENSION(3) :: cell_index, n_max
2107 REAL(kind=dp) :: avail_mem_per_proc_gb, cell_dist, cell_radius_3c, dij, dik, djk, eps, &
2108 exp_min_ao, exp_min_ri, frobenius_norm, mem_3c_gb, mem_occ_per_proc_gb, radius_ao, &
2109 radius_ao_product, radius_ri
2110 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exp_ao_kind, exp_ri_kind, &
2111 radius_ao_kind, &
2112 radius_ao_product_kind, radius_ri_kind
2113 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c
2114 REAL(kind=dp), DIMENSION(3) :: rij, rik, rjk, vec_cell_j, vec_cell_k
2115 REAL(kind=dp), DIMENSION(:, :), POINTER :: exp_ao, exp_ri
2116 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2117 TYPE(cell_type), POINTER :: cell
2118 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2119
2120 CALL timeset(routinen, handle)
2121
2122 CALL get_qs_env(qs_env, nkind=nkind, atomic_kind_set=atomic_kind_set, particle_set=particle_set, cell=cell)
2123
2124 ALLOCATE (exp_ao_kind(nkind), exp_ri_kind(nkind), radius_ao_kind(nkind), &
2125 radius_ao_product_kind(nkind), radius_ri_kind(nkind))
2126
2127 exp_min_ri = 10.0_dp
2128 exp_min_ao = 10.0_dp
2129 exp_ri_kind = 10.0_dp
2130 exp_ao_kind = 10.0_dp
2131
2132 eps = bs_env%eps_filter*bs_env%heuristic_filter_factor
2133
2134 DO ikind = 1, nkind
2135
2136 CALL get_gto_basis_set(bs_env%basis_set_RI(ikind)%gto_basis_set, zet=exp_ri)
2137 CALL get_gto_basis_set(bs_env%basis_set_ao(ikind)%gto_basis_set, zet=exp_ao)
2138
2139 ! we need to remove all exponents lower than a lower bound, e.g. 1E-3, because
2140 ! for contracted basis sets, there might be exponents = 0 in zet
2141 DO i = 1, SIZE(exp_ri, 1)
2142 DO j = 1, SIZE(exp_ri, 2)
2143 IF (exp_ri(i, j) < exp_min_ri .AND. exp_ri(i, j) > 1e-3_dp) exp_min_ri = exp_ri(i, j)
2144 IF (exp_ri(i, j) < exp_ri_kind(ikind) .AND. exp_ri(i, j) > 1e-3_dp) THEN
2145 exp_ri_kind(ikind) = exp_ri(i, j)
2146 END IF
2147 END DO
2148 END DO
2149 DO i = 1, SIZE(exp_ao, 1)
2150 DO j = 1, SIZE(exp_ao, 2)
2151 IF (exp_ao(i, j) < exp_min_ao .AND. exp_ao(i, j) > 1e-3_dp) exp_min_ao = exp_ao(i, j)
2152 IF (exp_ao(i, j) < exp_ao_kind(ikind) .AND. exp_ao(i, j) > 1e-3_dp) THEN
2153 exp_ao_kind(ikind) = exp_ao(i, j)
2154 END IF
2155 END DO
2156 END DO
2157 radius_ao_kind(ikind) = sqrt(-log(eps)/exp_ao_kind(ikind))
2158 radius_ao_product_kind(ikind) = sqrt(-log(eps)/(2.0_dp*exp_ao_kind(ikind)))
2159 radius_ri_kind(ikind) = sqrt(-log(eps)/exp_ri_kind(ikind))
2160 END DO
2161
2162 radius_ao = sqrt(-log(eps)/exp_min_ao)
2163 radius_ao_product = sqrt(-log(eps)/(2.0_dp*exp_min_ao))
2164 radius_ri = sqrt(-log(eps)/exp_min_ri)
2165
2166 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
2167
2168 ! For a 3c integral (μR υS | P0) we have that cell R and cell S need to be within radius_3c
2169 cell_radius_3c = radius_ao_product + radius_ri + bs_env%ri_metric%cutoff_radius
2170
2171 n_max(1:3) = bs_env%periodic(1:3)*30
2172
2173 nimages_3c_max = 0
2174
2175 i_cell_x_min = 0
2176 i_cell_x_max = 0
2177 j_cell_y_min = 0
2178 j_cell_y_max = 0
2179 k_cell_z_min = 0
2180 k_cell_z_max = 0
2181
2182 DO i_cell_x = -n_max(1), n_max(1)
2183 DO j_cell_y = -n_max(2), n_max(2)
2184 DO k_cell_z = -n_max(3), n_max(3)
2185
2186 cell_index(1:3) = [i_cell_x, j_cell_y, k_cell_z]
2187
2188 CALL get_cell_dist(cell_index, bs_env%hmat, cell_dist)
2189
2190 IF (cell_dist < cell_radius_3c) THEN
2191 nimages_3c_max = nimages_3c_max + 1
2192 i_cell_x_min = min(i_cell_x_min, i_cell_x)
2193 i_cell_x_max = max(i_cell_x_max, i_cell_x)
2194 j_cell_y_min = min(j_cell_y_min, j_cell_y)
2195 j_cell_y_max = max(j_cell_y_max, j_cell_y)
2196 k_cell_z_min = min(k_cell_z_min, k_cell_z)
2197 k_cell_z_max = max(k_cell_z_max, k_cell_z)
2198 END IF
2199
2200 END DO
2201 END DO
2202 END DO
2203
2204 ! get index_to_cell_3c_max for the maximum possible cell range;
2205 ! compute 3c integrals later in this routine and check really which cell is needed
2206 ALLOCATE (index_to_cell_3c_max(3, nimages_3c_max))
2207
2208 img = 0
2209 DO i_cell_x = -n_max(1), n_max(1)
2210 DO j_cell_y = -n_max(2), n_max(2)
2211 DO k_cell_z = -n_max(3), n_max(3)
2212
2213 cell_index(1:3) = [i_cell_x, j_cell_y, k_cell_z]
2214
2215 CALL get_cell_dist(cell_index, bs_env%hmat, cell_dist)
2216
2217 IF (cell_dist < cell_radius_3c) THEN
2218 img = img + 1
2219 index_to_cell_3c_max(1:3, img) = cell_index(1:3)
2220 END IF
2221
2222 END DO
2223 END DO
2224 END DO
2225
2226 ! get pairs of R and S which have non-zero 3c integral (μR υS | P0)
2227 ALLOCATE (nblocks_3c_max(nimages_3c_max, nimages_3c_max))
2228 nblocks_3c_max(:, :) = 0
2229
2230 block_count = 0
2231 DO j_cell = 1, nimages_3c_max
2232 DO k_cell = 1, nimages_3c_max
2233
2234 DO atom_j = 1, bs_env%n_atom
2235 DO atom_k = 1, bs_env%n_atom
2236 DO atom_i = 1, bs_env%n_atom
2237
2238 block_count = block_count + 1
2239 IF (modulo(block_count, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
2240
2241 CALL scaled_to_real(vec_cell_j, real(index_to_cell_3c_max(1:3, j_cell), kind=dp), cell)
2242 CALL scaled_to_real(vec_cell_k, real(index_to_cell_3c_max(1:3, k_cell), kind=dp), cell)
2243
2244 rij = pbc(particle_set(atom_j)%r(:), cell) - pbc(particle_set(atom_i)%r(:), cell) + vec_cell_j(:)
2245 rjk = pbc(particle_set(atom_k)%r(:), cell) - pbc(particle_set(atom_j)%r(:), cell) &
2246 + vec_cell_k(:) - vec_cell_j(:)
2247 rik(:) = rij(:) + rjk(:)
2248 dij = norm2(rij)
2249 dik = norm2(rik)
2250 djk = norm2(rjk)
2251 IF (djk > radius_ao_kind(kind_of(atom_j)) + radius_ao_kind(kind_of(atom_k))) cycle
2252 IF (dij > radius_ao_kind(kind_of(atom_j)) + radius_ri_kind(kind_of(atom_i)) &
2253 + bs_env%ri_metric%cutoff_radius) cycle
2254 IF (dik > radius_ri_kind(kind_of(atom_i)) + radius_ao_kind(kind_of(atom_k)) &
2255 + bs_env%ri_metric%cutoff_radius) cycle
2256
2257 j_size = bs_env%i_ao_end_from_atom(atom_j) - bs_env%i_ao_start_from_atom(atom_j) + 1
2258 k_size = bs_env%i_ao_end_from_atom(atom_k) - bs_env%i_ao_start_from_atom(atom_k) + 1
2259 i_size = bs_env%i_RI_end_from_atom(atom_i) - bs_env%i_RI_start_from_atom(atom_i) + 1
2260
2261 ALLOCATE (int_3c(j_size, k_size, i_size))
2262
2263 ! compute 3-c int. ( μ(atom j) R , ν (atom k) S | P (atom i) 0 )
2264 ! ("|": truncated Coulomb operator), inside build_3c_integrals: (j k | i)
2265 CALL build_3c_integral_block(int_3c, qs_env, bs_env%ri_metric, &
2266 basis_j=bs_env%basis_set_AO, &
2267 basis_k=bs_env%basis_set_AO, &
2268 basis_i=bs_env%basis_set_RI, &
2269 cell_j=index_to_cell_3c_max(1:3, j_cell), &
2270 cell_k=index_to_cell_3c_max(1:3, k_cell), &
2271 atom_k=atom_k, atom_j=atom_j, atom_i=atom_i)
2272
2273 frobenius_norm = sqrt(sum(int_3c(:, :, :)**2))
2274
2275 DEALLOCATE (int_3c)
2276
2277 ! we use a higher threshold here to safe memory when storing the 3c integrals
2278 ! in every tensor group
2279 IF (frobenius_norm > eps) THEN
2280 nblocks_3c_max(j_cell, k_cell) = nblocks_3c_max(j_cell, k_cell) + 1
2281 END IF
2282
2283 END DO
2284 END DO
2285 END DO
2286
2287 END DO
2288 END DO
2289
2290 CALL bs_env%para_env%sum(nblocks_3c_max)
2291
2292 ALLOCATE (n_other_3c_images_max(nimages_3c_max))
2293 n_other_3c_images_max(:) = 0
2294
2295 nimages_3c = 0
2296 nimage_pairs_3c = 0
2297
2298 DO j_cell = 1, nimages_3c_max
2299 DO k_cell = 1, nimages_3c_max
2300 IF (nblocks_3c_max(j_cell, k_cell) > 0) THEN
2301 n_other_3c_images_max(j_cell) = n_other_3c_images_max(j_cell) + 1
2302 nimage_pairs_3c = nimage_pairs_3c + 1
2303 END IF
2304 END DO
2305
2306 IF (n_other_3c_images_max(j_cell) > 0) nimages_3c = nimages_3c + 1
2307
2308 END DO
2309
2310 bs_env%nimages_3c = nimages_3c
2311 ALLOCATE (bs_env%index_to_cell_3c(3, nimages_3c))
2312 ALLOCATE (bs_env%cell_to_index_3c(i_cell_x_min:i_cell_x_max, &
2313 j_cell_y_min:j_cell_y_max, &
2314 k_cell_z_min:k_cell_z_max))
2315 bs_env%cell_to_index_3c(:, :, :) = -1
2316
2317 ALLOCATE (bs_env%nblocks_3c(nimages_3c, nimages_3c))
2318 bs_env%nblocks_3c(nimages_3c, nimages_3c) = 0
2319
2320 j_cell = 0
2321 DO j_cell_max = 1, nimages_3c_max
2322 IF (n_other_3c_images_max(j_cell_max) == 0) cycle
2323 j_cell = j_cell + 1
2324 cell_index(1:3) = index_to_cell_3c_max(1:3, j_cell_max)
2325 bs_env%index_to_cell_3c(1:3, j_cell) = cell_index(1:3)
2326 bs_env%cell_to_index_3c(cell_index(1), cell_index(2), cell_index(3)) = j_cell
2327
2328 k_cell = 0
2329 DO k_cell_max = 1, nimages_3c_max
2330 IF (n_other_3c_images_max(k_cell_max) == 0) cycle
2331 k_cell = k_cell + 1
2332
2333 bs_env%nblocks_3c(j_cell, k_cell) = nblocks_3c_max(j_cell_max, k_cell_max)
2334 END DO
2335
2336 END DO
2337
2338 ! we use: 8*10^-9 GB / double precision number
2339 mem_3c_gb = real(bs_env%n_RI, kind=dp)*real(bs_env%n_ao, kind=dp)**2 &
2340 *real(nimage_pairs_3c, kind=dp)*8e-9_dp
2341
2342 CALL m_memory(mem_occ_per_proc)
2343 CALL bs_env%para_env%max(mem_occ_per_proc)
2344
2345 mem_occ_per_proc_gb = real(mem_occ_per_proc, kind=dp)/1.0e9_dp
2346
2347 ! number of processors per group that entirely stores the 3c integrals and does tensor ops
2348 avail_mem_per_proc_gb = bs_env%input_memory_per_proc_GB - mem_occ_per_proc_gb
2349
2350 ! careful: downconvering real to integer, 1.9 -> 1; thus add 1.0 for upconversion, 1.9 -> 2
2351 bs_env%group_size_tensor = max(int(mem_3c_gb/avail_mem_per_proc_gb + 1.0_dp), 1)
2352
2353 u = bs_env%unit_nr
2354
2355 IF (u > 0) THEN
2356 WRITE (u, fmt="(T2,A,F52.1,A)") "Radius of atomic orbitals", radius_ao*angstrom, Å" "
2357 WRITE (u, fmt="(T2,A,F55.1,A)") "Radius of RI functions", radius_ri*angstrom, Å" "
2358 WRITE (u, fmt="(T2,A,I47)") "Number of cells for 3c integrals", nimages_3c
2359 WRITE (u, fmt="(T2,A,I42)") "Number of cell pairs for 3c integrals", nimage_pairs_3c
2360 WRITE (u, '(T2,A)') ''
2361 WRITE (u, '(T2,A,F37.1,A)') 'Input: Available memory per MPI process', &
2362 bs_env%input_memory_per_proc_GB, ' GB'
2363 WRITE (u, '(T2,A,F35.1,A)') 'Used memory per MPI process before GW run', &
2364 mem_occ_per_proc_gb, ' GB'
2365 WRITE (u, '(T2,A,F44.1,A)') 'Memory of three-center integrals', mem_3c_gb, ' GB'
2366 END IF
2367
2368 CALL timestop(handle)
2369
2370 END SUBROUTINE setup_cells_3c
2371
2372! **************************************************************************************************
2373!> \brief ...
2374!> \param index_to_cell_1 ...
2375!> \param index_to_cell_2 ...
2376!> \param nimages_1 ...
2377!> \param nimages_2 ...
2378!> \param index_to_cell ...
2379!> \param cell_to_index ...
2380!> \param nimages ...
2381! **************************************************************************************************
2382 SUBROUTINE sum_two_r_grids(index_to_cell_1, index_to_cell_2, nimages_1, nimages_2, &
2383 index_to_cell, cell_to_index, nimages)
2384
2385 INTEGER, DIMENSION(:, :) :: index_to_cell_1, index_to_cell_2
2386 INTEGER :: nimages_1, nimages_2
2387 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
2388 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2389 INTEGER :: nimages
2390
2391 CHARACTER(LEN=*), PARAMETER :: routinen = 'sum_two_R_grids'
2392
2393 INTEGER :: handle, i_dim, img_1, img_2, nimages_max
2394 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_tmp
2395 INTEGER, DIMENSION(3) :: cell_1, cell_2, r, r_max, r_min
2396
2397 CALL timeset(routinen, handle)
2398
2399 DO i_dim = 1, 3
2400 r_min(i_dim) = minval(index_to_cell_1(i_dim, :)) + minval(index_to_cell_2(i_dim, :))
2401 r_max(i_dim) = maxval(index_to_cell_1(i_dim, :)) + maxval(index_to_cell_2(i_dim, :))
2402 END DO
2403
2404 nimages_max = (r_max(1) - r_min(1) + 1)*(r_max(2) - r_min(2) + 1)*(r_max(3) - r_min(3) + 1)
2405
2406 ALLOCATE (index_to_cell_tmp(3, nimages_max))
2407 index_to_cell_tmp(:, :) = -1
2408
2409 ALLOCATE (cell_to_index(r_min(1):r_max(1), r_min(2):r_max(2), r_min(3):r_max(3)))
2410 cell_to_index(:, :, :) = -1
2411
2412 nimages = 0
2413
2414 DO img_1 = 1, nimages_1
2415
2416 DO img_2 = 1, nimages_2
2417
2418 cell_1(1:3) = index_to_cell_1(1:3, img_1)
2419 cell_2(1:3) = index_to_cell_2(1:3, img_2)
2420
2421 r(1:3) = cell_1(1:3) + cell_2(1:3)
2422
2423 ! check whether we have found a new cell
2424 IF (cell_to_index(r(1), r(2), r(3)) == -1) THEN
2425
2426 nimages = nimages + 1
2427 cell_to_index(r(1), r(2), r(3)) = nimages
2428 index_to_cell_tmp(1:3, nimages) = r(1:3)
2429
2430 END IF
2431
2432 END DO
2433
2434 END DO
2435
2436 ALLOCATE (index_to_cell(3, nimages))
2437 index_to_cell(:, :) = index_to_cell_tmp(1:3, 1:nimages)
2438
2439 CALL timestop(handle)
2440
2441 END SUBROUTINE sum_two_r_grids
2442
2443! **************************************************************************************************
2444!> \brief ...
2445!> \param qs_env ...
2446!> \param bs_env ...
2447! **************************************************************************************************
2448 SUBROUTINE compute_3c_integrals(qs_env, bs_env)
2449
2450 TYPE(qs_environment_type), POINTER :: qs_env
2451 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2452
2453 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_3c_integrals'
2454
2455 INTEGER :: handle, j_cell, k_cell, nimages_3c
2456
2457 CALL timeset(routinen, handle)
2458
2459 nimages_3c = bs_env%nimages_3c
2460 ALLOCATE (bs_env%t_3c_int(nimages_3c, nimages_3c))
2461 DO j_cell = 1, nimages_3c
2462 DO k_cell = 1, nimages_3c
2463 CALL dbt_create(bs_env%t_RI_AO__AO, bs_env%t_3c_int(j_cell, k_cell))
2464 END DO
2465 END DO
2466
2467 CALL build_3c_integrals(bs_env%t_3c_int, &
2468 bs_env%eps_filter, &
2469 qs_env, &
2470 bs_env%nl_3c, &
2471 int_eps=bs_env%eps_filter*0.05_dp, &
2472 basis_i=bs_env%basis_set_RI, &
2473 basis_j=bs_env%basis_set_AO, &
2474 basis_k=bs_env%basis_set_AO, &
2475 potential_parameter=bs_env%ri_metric, &
2476 desymmetrize=.false., do_kpoints=.true., cell_sym=.true., &
2477 cell_to_index_ext=bs_env%cell_to_index_3c)
2478
2479 CALL bs_env%para_env%sync()
2480
2481 CALL timestop(handle)
2482
2483 END SUBROUTINE compute_3c_integrals
2484
2485! **************************************************************************************************
2486!> \brief ...
2487!> \param cell_index ...
2488!> \param hmat ...
2489!> \param cell_dist ...
2490! **************************************************************************************************
2491 SUBROUTINE get_cell_dist(cell_index, hmat, cell_dist)
2492
2493 INTEGER, DIMENSION(3) :: cell_index
2494 REAL(kind=dp) :: hmat(3, 3), cell_dist
2495
2496 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_cell_dist'
2497
2498 INTEGER :: handle, i_dim
2499 INTEGER, DIMENSION(3) :: cell_index_adj
2500 REAL(kind=dp) :: cell_dist_3(3)
2501
2502 CALL timeset(routinen, handle)
2503
2504 ! the distance of cells needs to be taken to adjacent neighbors, not
2505 ! between the center of the cells. We thus need to rescale the cell index
2506 DO i_dim = 1, 3
2507 IF (cell_index(i_dim) > 0) cell_index_adj(i_dim) = cell_index(i_dim) - 1
2508 IF (cell_index(i_dim) < 0) cell_index_adj(i_dim) = cell_index(i_dim) + 1
2509 IF (cell_index(i_dim) == 0) cell_index_adj(i_dim) = cell_index(i_dim)
2510 END DO
2511
2512 cell_dist_3(1:3) = matmul(hmat, real(cell_index_adj, kind=dp))
2513
2514 cell_dist = sqrt(abs(sum(cell_dist_3(1:3)**2)))
2515
2516 CALL timestop(handle)
2517
2518 END SUBROUTINE get_cell_dist
2519
2520! **************************************************************************************************
2521!> \brief ...
2522!> \param qs_env ...
2523!> \param bs_env ...
2524!> \param kpoints ...
2525!> \param do_print ...
2526! **************************************************************************************************
2527 SUBROUTINE setup_kpoints_scf_desymm(qs_env, bs_env, kpoints, do_print)
2528 TYPE(qs_environment_type), POINTER :: qs_env
2529 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2530 TYPE(kpoint_type), POINTER :: kpoints
2531
2532 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_kpoints_scf_desymm'
2533
2534 INTEGER :: handle, i_cell_x, i_dim, img, j_cell_y, &
2535 k_cell_z, nimages, nkp, u
2536 INTEGER, DIMENSION(3) :: cell_grid, cixd, nkp_grid
2537 TYPE(kpoint_type), POINTER :: kpoints_scf
2538
2539 LOGICAL:: do_print
2540
2541 CALL timeset(routinen, handle)
2542
2543 NULLIFY (kpoints)
2544 CALL kpoint_create(kpoints)
2545
2546 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints_scf)
2547
2548 nkp_grid(1:3) = kpoints_scf%nkp_grid(1:3)
2549 nkp = nkp_grid(1)*nkp_grid(2)*nkp_grid(3)
2550
2551 ! we need in periodic directions at least 2 k-points in the SCF
2552 DO i_dim = 1, 3
2553 IF (bs_env%periodic(i_dim) == 1) THEN
2554 cpassert(nkp_grid(i_dim) > 1)
2555 END IF
2556 END DO
2557
2558 kpoints%kp_scheme = "GENERAL"
2559 kpoints%nkp_grid(1:3) = nkp_grid(1:3)
2560 kpoints%nkp = nkp
2561 bs_env%nkp_scf_desymm = nkp
2562
2563 ALLOCATE (kpoints%xkp(1:3, nkp))
2564 CALL compute_xkp(kpoints%xkp, 1, nkp, nkp_grid)
2565
2566 ALLOCATE (kpoints%wkp(nkp))
2567 kpoints%wkp(:) = 1.0_dp/real(nkp, kind=dp)
2568
2569 ! for example 4x3x6 kpoint grid -> 3x3x5 cell grid because we need the same number of
2570 ! neighbor cells on both sides of the unit cell
2571 cell_grid(1:3) = nkp_grid(1:3) - modulo(nkp_grid(1:3) + 1, 2)
2572 ! cell index: for example for x: from -n_x/2 to +n_x/2, n_x: number of cells in x direction
2573 cixd(1:3) = cell_grid(1:3)/2
2574
2575 nimages = cell_grid(1)*cell_grid(2)*cell_grid(3)
2576
2577 bs_env%nimages_scf_desymm = nimages
2578
2579 ALLOCATE (kpoints%cell_to_index(-cixd(1):cixd(1), -cixd(2):cixd(2), -cixd(3):cixd(3)))
2580 ALLOCATE (kpoints%index_to_cell(3, nimages))
2581
2582 img = 0
2583 DO i_cell_x = -cixd(1), cixd(1)
2584 DO j_cell_y = -cixd(2), cixd(2)
2585 DO k_cell_z = -cixd(3), cixd(3)
2586 img = img + 1
2587 kpoints%cell_to_index(i_cell_x, j_cell_y, k_cell_z) = img
2588 kpoints%index_to_cell(1:3, img) = [i_cell_x, j_cell_y, k_cell_z]
2589 END DO
2590 END DO
2591 END DO
2592
2593 u = bs_env%unit_nr
2594 IF (u > 0 .AND. do_print) THEN
2595 WRITE (u, fmt="(T2,A,I49)") χΣ"Number of cells for G, , W, ", nimages
2596 END IF
2597
2598 CALL timestop(handle)
2599
2600 END SUBROUTINE setup_kpoints_scf_desymm
2601
2602! **************************************************************************************************
2603!> \brief ...
2604!> \param bs_env ...
2605! **************************************************************************************************
2606 SUBROUTINE setup_cells_delta_r(bs_env)
2607
2608 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2609
2610 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_cells_Delta_R'
2611
2612 INTEGER :: handle
2613
2614 CALL timeset(routinen, handle)
2615
2616 ! cell sums batch wise for fixed ΔR = S_1 - R_1; for example:
2617 ! Σ_λσ^R = sum_PR1νS1 M^G_λ0,νS1,PR1 M^W_σR,νS1,PR1
2618
2619 CALL sum_two_r_grids(bs_env%index_to_cell_3c, &
2620 bs_env%index_to_cell_3c, &
2621 bs_env%nimages_3c, bs_env%nimages_3c, &
2622 bs_env%index_to_cell_Delta_R, &
2623 bs_env%cell_to_index_Delta_R, &
2624 bs_env%nimages_Delta_R)
2625
2626 IF (bs_env%unit_nr > 0) THEN
2627 WRITE (bs_env%unit_nr, fmt="(T2,A,I61)") Δ"Number of cells R", bs_env%nimages_Delta_R
2628 END IF
2629
2630 CALL timestop(handle)
2631
2632 END SUBROUTINE setup_cells_delta_r
2633
2634! **************************************************************************************************
2635!> \brief ...
2636!> \param bs_env ...
2637! **************************************************************************************************
2638 SUBROUTINE setup_parallelization_delta_r(bs_env)
2639
2640 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2641
2642 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_parallelization_Delta_R'
2643
2644 INTEGER :: handle, i_cell_delta_r, i_task_local, &
2645 n_tasks_local
2646 INTEGER, ALLOCATABLE, DIMENSION(:) :: i_cell_delta_r_group, &
2647 n_tensor_ops_delta_r
2648
2649 CALL timeset(routinen, handle)
2650
2651 CALL compute_n_tensor_ops_delta_r(bs_env, n_tensor_ops_delta_r)
2652
2653 CALL compute_delta_r_dist(bs_env, n_tensor_ops_delta_r, i_cell_delta_r_group, n_tasks_local)
2654
2655 bs_env%n_tasks_Delta_R_local = n_tasks_local
2656
2657 ALLOCATE (bs_env%task_Delta_R(n_tasks_local))
2658
2659 i_task_local = 0
2660 DO i_cell_delta_r = 1, bs_env%nimages_Delta_R
2661
2662 IF (i_cell_delta_r_group(i_cell_delta_r) /= bs_env%tensor_group_color) cycle
2663
2664 i_task_local = i_task_local + 1
2665
2666 bs_env%task_Delta_R(i_task_local) = i_cell_delta_r
2667
2668 END DO
2669
2670 ALLOCATE (bs_env%skip_DR_chi(n_tasks_local))
2671 bs_env%skip_DR_chi(:) = .false.
2672 ALLOCATE (bs_env%skip_DR_Sigma(n_tasks_local))
2673 bs_env%skip_DR_Sigma(:) = .false.
2674
2675 CALL allocate_skip_3xr(bs_env%skip_DR_R12_S_Goccx3c_chi, bs_env)
2676 CALL allocate_skip_3xr(bs_env%skip_DR_R12_S_Gvirx3c_chi, bs_env)
2677 CALL allocate_skip_3xr(bs_env%skip_DR_R_R2_MxM_chi, bs_env)
2678
2679 CALL allocate_skip_3xr(bs_env%skip_DR_R1_S2_Gx3c_Sigma, bs_env)
2680 CALL allocate_skip_3xr(bs_env%skip_DR_R1_R_MxM_Sigma, bs_env)
2681
2682 CALL timestop(handle)
2683
2684 END SUBROUTINE setup_parallelization_delta_r
2685
2686! **************************************************************************************************
2687!> \brief ...
2688!> \param skip ...
2689!> \param bs_env ...
2690! **************************************************************************************************
2691 SUBROUTINE allocate_skip_3xr(skip, bs_env)
2692 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: skip
2693 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2694
2695 CHARACTER(LEN=*), PARAMETER :: routinen = 'allocate_skip_3xR'
2696
2697 INTEGER :: handle
2698
2699 CALL timeset(routinen, handle)
2700
2701 ALLOCATE (skip(bs_env%n_tasks_Delta_R_local, bs_env%nimages_3c, bs_env%nimages_scf_desymm))
2702 skip(:, :, :) = .false.
2703
2704 CALL timestop(handle)
2705
2706 END SUBROUTINE allocate_skip_3xr
2707
2708! **************************************************************************************************
2709!> \brief ...
2710!> \param bs_env ...
2711!> \param n_tensor_ops_Delta_R ...
2712!> \param i_cell_Delta_R_group ...
2713!> \param n_tasks_local ...
2714! **************************************************************************************************
2715 SUBROUTINE compute_delta_r_dist(bs_env, n_tensor_ops_Delta_R, i_cell_Delta_R_group, n_tasks_local)
2716 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2717 INTEGER, ALLOCATABLE, DIMENSION(:) :: n_tensor_ops_delta_r, &
2718 i_cell_delta_r_group
2719 INTEGER :: n_tasks_local
2720
2721 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_Delta_R_dist'
2722
2723 INTEGER :: handle, i_delta_r_max_op, i_group_min, &
2724 nimages_delta_r, u
2725 INTEGER, ALLOCATABLE, DIMENSION(:) :: n_tensor_ops_delta_r_in_group
2726
2727 CALL timeset(routinen, handle)
2728
2729 nimages_delta_r = bs_env%nimages_Delta_R
2730
2731 u = bs_env%unit_nr
2732
2733 IF (u > 0 .AND. nimages_delta_r < bs_env%num_tensor_groups) THEN
2734 WRITE (u, fmt="(T2,A,I5,A,I5,A)") "There are only ", nimages_delta_r, &
2735 " tasks to work on but there are ", bs_env%num_tensor_groups, " groups."
2736 WRITE (u, fmt="(T2,A)") "Please reduce the number of MPI processes."
2737 WRITE (u, '(T2,A)') ''
2738 END IF
2739
2740 ALLOCATE (n_tensor_ops_delta_r_in_group(bs_env%num_tensor_groups))
2741 n_tensor_ops_delta_r_in_group(:) = 0
2742 ALLOCATE (i_cell_delta_r_group(nimages_delta_r))
2743 i_cell_delta_r_group(:) = -1
2744
2745 n_tasks_local = 0
2746
2747 DO WHILE (any(n_tensor_ops_delta_r(:) /= 0))
2748
2749 ! get largest element of n_tensor_ops_Delta_R
2750 i_delta_r_max_op = maxloc(n_tensor_ops_delta_r, 1)
2751
2752 ! distribute i_Delta_R_max_op to tensor group which has currently the smallest load
2753 i_group_min = minloc(n_tensor_ops_delta_r_in_group, 1)
2754
2755 ! the tensor groups are 0-index based; but i_group_min is 1-index based
2756 i_cell_delta_r_group(i_delta_r_max_op) = i_group_min - 1
2757 n_tensor_ops_delta_r_in_group(i_group_min) = n_tensor_ops_delta_r_in_group(i_group_min) + &
2758 n_tensor_ops_delta_r(i_delta_r_max_op)
2759
2760 ! remove i_Delta_R_max_op from n_tensor_ops_Delta_R
2761 n_tensor_ops_delta_r(i_delta_r_max_op) = 0
2762
2763 IF (bs_env%tensor_group_color == i_group_min - 1) n_tasks_local = n_tasks_local + 1
2764
2765 END DO
2766
2767 CALL timestop(handle)
2768
2769 END SUBROUTINE compute_delta_r_dist
2770
2771! **************************************************************************************************
2772!> \brief ...
2773!> \param bs_env ...
2774!> \param n_tensor_ops_Delta_R ...
2775! **************************************************************************************************
2776 SUBROUTINE compute_n_tensor_ops_delta_r(bs_env, n_tensor_ops_Delta_R)
2777 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2778 INTEGER, ALLOCATABLE, DIMENSION(:) :: n_tensor_ops_delta_r
2779
2780 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_n_tensor_ops_Delta_R'
2781
2782 INTEGER :: handle, i_cell_delta_r, i_cell_r, i_cell_r1, i_cell_r1_minus_r, i_cell_r2, &
2783 i_cell_r2_m_r1, i_cell_s1, i_cell_s1_m_r1_p_r2, i_cell_s1_minus_r, i_cell_s2, &
2784 nimages_delta_r
2785 INTEGER, DIMENSION(3) :: cell_dr, cell_m_r1, cell_r, cell_r1, cell_r1_minus_r, cell_r2, &
2786 cell_r2_m_r1, cell_s1, cell_s1_m_r2_p_r1, cell_s1_minus_r, cell_s1_p_s2_m_r1, cell_s2
2787 LOGICAL :: cell_found
2788
2789 CALL timeset(routinen, handle)
2790
2791 nimages_delta_r = bs_env%nimages_Delta_R
2792
2793 ALLOCATE (n_tensor_ops_delta_r(nimages_delta_r))
2794 n_tensor_ops_delta_r(:) = 0
2795
2796 ! compute number of tensor operations for specific Delta_R
2797 DO i_cell_delta_r = 1, nimages_delta_r
2798
2799 IF (modulo(i_cell_delta_r, bs_env%num_tensor_groups) /= bs_env%tensor_group_color) cycle
2800
2801 DO i_cell_r1 = 1, bs_env%nimages_3c
2802
2803 cell_r1(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_r1)
2804 cell_dr(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_delta_r)
2805
2806 ! S_1 = R_1 + ΔR (from ΔR = S_1 - R_1)
2807 CALL add_r(cell_r1, cell_dr, bs_env%index_to_cell_3c, cell_s1, &
2808 cell_found, bs_env%cell_to_index_3c, i_cell_s1)
2809 IF (.NOT. cell_found) cycle
2810
2811 DO i_cell_r2 = 1, bs_env%nimages_scf_desymm
2812
2813 cell_r2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_r2)
2814
2815 ! R_2 - R_1
2816 CALL add_r(cell_r2, -cell_r1, bs_env%index_to_cell_3c, cell_r2_m_r1, &
2817 cell_found, bs_env%cell_to_index_3c, i_cell_r2_m_r1)
2818 IF (.NOT. cell_found) cycle
2819
2820 ! S_1 - R_1 + R_2
2821 CALL add_r(cell_s1, cell_r2_m_r1, bs_env%index_to_cell_3c, cell_s1_m_r2_p_r1, &
2822 cell_found, bs_env%cell_to_index_3c, i_cell_s1_m_r1_p_r2)
2823 IF (.NOT. cell_found) cycle
2824
2825 n_tensor_ops_delta_r(i_cell_delta_r) = n_tensor_ops_delta_r(i_cell_delta_r) + 1
2826
2827 END DO ! i_cell_R2
2828
2829 DO i_cell_s2 = 1, bs_env%nimages_scf_desymm
2830
2831 cell_s2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_s2)
2832 cell_m_r1(1:3) = -cell_r1(1:3)
2833 cell_s1_p_s2_m_r1(1:3) = cell_s1(1:3) + cell_s2(1:3) - cell_r1(1:3)
2834
2835 CALL is_cell_in_index_to_cell(cell_m_r1, bs_env%index_to_cell_3c, cell_found)
2836 IF (.NOT. cell_found) cycle
2837
2838 CALL is_cell_in_index_to_cell(cell_s1_p_s2_m_r1, bs_env%index_to_cell_3c, cell_found)
2839 IF (.NOT. cell_found) cycle
2840
2841 END DO ! i_cell_S2
2842
2843 DO i_cell_r = 1, bs_env%nimages_scf_desymm
2844
2845 cell_r = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_r)
2846
2847 ! R_1 - R
2848 CALL add_r(cell_r1, -cell_r, bs_env%index_to_cell_3c, cell_r1_minus_r, &
2849 cell_found, bs_env%cell_to_index_3c, i_cell_r1_minus_r)
2850 IF (.NOT. cell_found) cycle
2851
2852 ! S_1 - R
2853 CALL add_r(cell_s1, -cell_r, bs_env%index_to_cell_3c, cell_s1_minus_r, &
2854 cell_found, bs_env%cell_to_index_3c, i_cell_s1_minus_r)
2855 IF (.NOT. cell_found) cycle
2856
2857 END DO ! i_cell_R
2858
2859 END DO ! i_cell_R1
2860
2861 END DO ! i_cell_Delta_R
2862
2863 CALL bs_env%para_env%sum(n_tensor_ops_delta_r)
2864
2865 CALL timestop(handle)
2866
2867 END SUBROUTINE compute_n_tensor_ops_delta_r
2868
2869! **************************************************************************************************
2870!> \brief ...
2871!> \param cell_1 ...
2872!> \param cell_2 ...
2873!> \param index_to_cell ...
2874!> \param cell_1_plus_2 ...
2875!> \param cell_found ...
2876!> \param cell_to_index ...
2877!> \param i_cell_1_plus_2 ...
2878! **************************************************************************************************
2879 SUBROUTINE add_r(cell_1, cell_2, index_to_cell, cell_1_plus_2, cell_found, &
2880 cell_to_index, i_cell_1_plus_2)
2881
2882 INTEGER, DIMENSION(3) :: cell_1, cell_2
2883 INTEGER, DIMENSION(:, :) :: index_to_cell
2884 INTEGER, DIMENSION(3) :: cell_1_plus_2
2885 LOGICAL :: cell_found
2886 INTEGER, DIMENSION(:, :, :), INTENT(IN), &
2887 OPTIONAL, POINTER :: cell_to_index
2888 INTEGER, INTENT(OUT), OPTIONAL :: i_cell_1_plus_2
2889
2890 CHARACTER(LEN=*), PARAMETER :: routinen = 'add_R'
2891
2892 INTEGER :: handle
2893
2894 CALL timeset(routinen, handle)
2895
2896 cell_1_plus_2(1:3) = cell_1(1:3) + cell_2(1:3)
2897
2898 CALL is_cell_in_index_to_cell(cell_1_plus_2, index_to_cell, cell_found)
2899
2900 IF (PRESENT(i_cell_1_plus_2)) THEN
2901 IF (cell_found) THEN
2902 cpassert(PRESENT(cell_to_index))
2903 i_cell_1_plus_2 = cell_to_index(cell_1_plus_2(1), cell_1_plus_2(2), cell_1_plus_2(3))
2904 ELSE
2905 i_cell_1_plus_2 = -1000
2906 END IF
2907 END IF
2908
2909 CALL timestop(handle)
2910
2911 END SUBROUTINE add_r
2912
2913! **************************************************************************************************
2914!> \brief ...
2915!> \param cell ...
2916!> \param index_to_cell ...
2917!> \param cell_found ...
2918! **************************************************************************************************
2919 SUBROUTINE is_cell_in_index_to_cell(cell, index_to_cell, cell_found)
2920 INTEGER, DIMENSION(3) :: cell
2921 INTEGER, DIMENSION(:, :) :: index_to_cell
2922 LOGICAL :: cell_found
2923
2924 CHARACTER(LEN=*), PARAMETER :: routinen = 'is_cell_in_index_to_cell'
2925
2926 INTEGER :: handle, i_cell, nimg
2927 INTEGER, DIMENSION(3) :: cell_i
2928
2929 CALL timeset(routinen, handle)
2930
2931 nimg = SIZE(index_to_cell, 2)
2932
2933 cell_found = .false.
2934
2935 DO i_cell = 1, nimg
2936
2937 cell_i(1:3) = index_to_cell(1:3, i_cell)
2938
2939 IF (cell_i(1) == cell(1) .AND. cell_i(2) == cell(2) .AND. cell_i(3) == cell(3)) THEN
2940 cell_found = .true.
2941 END IF
2942
2943 END DO
2944
2945 CALL timestop(handle)
2946
2947 END SUBROUTINE is_cell_in_index_to_cell
2948
2949! **************************************************************************************************
2950!> \brief ...
2951!> \param qs_env ...
2952!> \param bs_env ...
2953! **************************************************************************************************
2954 SUBROUTINE allocate_matrices_small_cell_full_kp(qs_env, bs_env)
2955 TYPE(qs_environment_type), POINTER :: qs_env
2956 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2957
2958 CHARACTER(LEN=*), PARAMETER :: routinen = 'allocate_matrices_small_cell_full_kp'
2959
2960 INTEGER :: handle, i_spin, i_t, img, n_spin, &
2961 nimages_scf, num_time_freq_points
2962 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2963 TYPE(mp_para_env_type), POINTER :: para_env
2964
2965 CALL timeset(routinen, handle)
2966
2967 nimages_scf = bs_env%nimages_scf_desymm
2968 num_time_freq_points = bs_env%num_time_freq_points
2969 n_spin = bs_env%n_spin
2970
2971 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
2972
2973 ALLOCATE (bs_env%fm_G_S(nimages_scf))
2974 ALLOCATE (bs_env%fm_Sigma_x_R(nimages_scf))
2975 ALLOCATE (bs_env%fm_chi_R_t(nimages_scf, num_time_freq_points))
2976 ALLOCATE (bs_env%fm_MWM_R_t(nimages_scf, num_time_freq_points))
2977 ALLOCATE (bs_env%fm_Sigma_c_R_neg_tau(nimages_scf, num_time_freq_points, n_spin))
2978 ALLOCATE (bs_env%fm_Sigma_c_R_pos_tau(nimages_scf, num_time_freq_points, n_spin))
2979 DO img = 1, nimages_scf
2980 CALL cp_fm_create(bs_env%fm_G_S(img), bs_env%fm_work_mo(1)%matrix_struct)
2981 CALL cp_fm_create(bs_env%fm_Sigma_x_R(img), bs_env%fm_work_mo(1)%matrix_struct)
2982 DO i_t = 1, num_time_freq_points
2983 CALL cp_fm_create(bs_env%fm_chi_R_t(img, i_t), bs_env%fm_RI_RI%matrix_struct)
2984 CALL cp_fm_create(bs_env%fm_MWM_R_t(img, i_t), bs_env%fm_RI_RI%matrix_struct)
2985 CALL cp_fm_set_all(bs_env%fm_MWM_R_t(img, i_t), 0.0_dp)
2986 DO i_spin = 1, n_spin
2987 CALL cp_fm_create(bs_env%fm_Sigma_c_R_neg_tau(img, i_t, i_spin), &
2988 bs_env%fm_work_mo(1)%matrix_struct)
2989 CALL cp_fm_create(bs_env%fm_Sigma_c_R_pos_tau(img, i_t, i_spin), &
2990 bs_env%fm_work_mo(1)%matrix_struct)
2991 CALL cp_fm_set_all(bs_env%fm_Sigma_c_R_neg_tau(img, i_t, i_spin), 0.0_dp)
2992 CALL cp_fm_set_all(bs_env%fm_Sigma_c_R_pos_tau(img, i_t, i_spin), 0.0_dp)
2993 END DO
2994 END DO
2995 END DO
2996
2997 CALL timestop(handle)
2998
2999 END SUBROUTINE allocate_matrices_small_cell_full_kp
3000
3001! **************************************************************************************************
3002!> \brief ...
3003!> \param qs_env ...
3004!> \param bs_env ...
3005! **************************************************************************************************
3006 SUBROUTINE trafo_v_xc_r_to_kp(qs_env, bs_env)
3007 TYPE(qs_environment_type), POINTER :: qs_env
3008 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3009
3010 CHARACTER(LEN=*), PARAMETER :: routinen = 'trafo_V_xc_R_to_kp'
3011
3012 INTEGER :: handle, ikp, img, ispin, n_ao
3013 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf
3014 TYPE(cp_cfm_type) :: cfm_mo_coeff, cfm_tmp, cfm_v_xc
3015 TYPE(cp_fm_type) :: fm_v_xc_re
3016 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
3017 TYPE(kpoint_type), POINTER :: kpoints_scf
3018 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3019 POINTER :: sab_nl
3020
3021 CALL timeset(routinen, handle)
3022
3023 n_ao = bs_env%n_ao
3024
3025 CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, kpoints=kpoints_scf)
3026
3027 NULLIFY (sab_nl)
3028 CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
3029
3030 CALL cp_cfm_create(cfm_v_xc, bs_env%cfm_work_mo%matrix_struct)
3031 CALL cp_cfm_create(cfm_mo_coeff, bs_env%cfm_work_mo%matrix_struct)
3032 CALL cp_cfm_create(cfm_tmp, bs_env%cfm_work_mo%matrix_struct)
3033 CALL cp_fm_create(fm_v_xc_re, bs_env%cfm_work_mo%matrix_struct)
3034
3035 DO img = 1, bs_env%nimages_scf
3036 DO ispin = 1, bs_env%n_spin
3037 ! JW kind of hack because the format of matrix_ks remains dubious...
3038 CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
3039 CALL copy_fm_to_dbcsr(bs_env%fm_V_xc_R(img, ispin), matrix_ks(ispin, img)%matrix)
3040 END DO
3041 END DO
3042
3043 ALLOCATE (bs_env%v_xc_n(n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
3044
3045 DO ispin = 1, bs_env%n_spin
3046 DO ikp = 1, bs_env%nkp_bs_and_DOS
3047
3048 ! v^xc^R -> v^xc(k) (matrix_ks stores v^xc^R, see SUBROUTINE compute_V_xc)
3049 CALL rsmat_to_kp(matrix_ks, ispin, bs_env%kpoints_DOS%xkp(1:3, ikp), &
3050 cell_to_index_scf, sab_nl, bs_env, cfm_v_xc)
3051
3052 ! get C_µn(k)
3053 CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mo_coeff)
3054
3055 ! v^xc_nm(k_i) = sum_µν C^*_µn(k_i) v^xc_µν(k_i) C_νn(k_i)
3056 CALL parallel_gemm('N', 'N', n_ao, n_ao, n_ao, z_one, cfm_v_xc, cfm_mo_coeff, &
3057 z_zero, cfm_tmp)
3058 CALL parallel_gemm('C', 'N', n_ao, n_ao, n_ao, z_one, cfm_mo_coeff, cfm_tmp, &
3059 z_zero, cfm_v_xc)
3060
3061 ! get v^xc_nn(k_i) which is a real quantity as v^xc is Hermitian
3062 CALL cp_cfm_to_fm(cfm_v_xc, fm_v_xc_re)
3063 CALL cp_fm_get_diag(fm_v_xc_re, bs_env%v_xc_n(:, ikp, ispin))
3064
3065 END DO
3066
3067 END DO
3068
3069 ! just rebuild the overwritten KS matrix again
3070 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., just_energy=.false.)
3071
3072 CALL cp_cfm_release(cfm_v_xc)
3073 CALL cp_cfm_release(cfm_mo_coeff)
3074 CALL cp_cfm_release(cfm_tmp)
3075 CALL cp_fm_release(fm_v_xc_re)
3076
3077 CALL timestop(handle)
3078
3079 END SUBROUTINE trafo_v_xc_r_to_kp
3080
3081! **************************************************************************************************
3082!> \brief ...
3083!> \param qs_env ...
3084!> \param bs_env ...
3085! **************************************************************************************************
3086 SUBROUTINE heuristic_ri_regularization(qs_env, bs_env)
3087 TYPE(qs_environment_type), POINTER :: qs_env
3088 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3089
3090 CHARACTER(LEN=*), PARAMETER :: routinen = 'heuristic_RI_regularization'
3091
3092 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: m
3093 INTEGER :: handle, ikp, ikp_local, n_ri, nkp, &
3094 nkp_local, u
3095 REAL(kind=dp) :: cond_nr, cond_nr_max, max_ev, &
3096 max_ev_ikp, min_ev, min_ev_ikp
3097 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: m_r
3098
3099 CALL timeset(routinen, handle)
3100
3101 ! compute M^R_PQ = <phi_P,0|V^tr(rc)|phi_Q,R> for RI metric
3102 CALL get_v_tr_r(m_r, bs_env%ri_metric, 0.0_dp, bs_env, qs_env)
3103
3104 nkp = bs_env%nkp_chi_eps_W_orig_plus_extra
3105 n_ri = bs_env%n_RI
3106
3107 nkp_local = 0
3108 DO ikp = 1, nkp
3109 ! trivial parallelization over k-points
3110 IF (modulo(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
3111 nkp_local = nkp_local + 1
3112 END DO
3113
3114 ALLOCATE (m(n_ri, n_ri, nkp_local))
3115
3116 ikp_local = 0
3117 cond_nr_max = 0.0_dp
3118 min_ev = 1000.0_dp
3119 max_ev = -1000.0_dp
3120
3121 DO ikp = 1, nkp
3122
3123 ! trivial parallelization
3124 IF (modulo(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
3125
3126 ikp_local = ikp_local + 1
3127
3128 ! M(k) = sum_R e^ikR M^R
3129 CALL rs_to_kp(m_r, m(:, :, ikp_local), &
3130 bs_env%kpoints_scf_desymm%index_to_cell, &
3131 bs_env%kpoints_chi_eps_W%xkp(1:3, ikp))
3132
3133 ! compute condition number of M_PQ(k)
3134 CALL power(m(:, :, ikp_local), 1.0_dp, 0.0_dp, cond_nr, min_ev_ikp, max_ev_ikp)
3135
3136 IF (cond_nr > cond_nr_max) cond_nr_max = cond_nr
3137 IF (max_ev_ikp > max_ev) max_ev = max_ev_ikp
3138 IF (min_ev_ikp < min_ev) min_ev = min_ev_ikp
3139
3140 END DO ! ikp
3141
3142 CALL bs_env%para_env%max(cond_nr_max)
3143 CALL bs_env%para_env%min(min_ev)
3144 CALL bs_env%para_env%max(max_ev)
3145
3146 u = bs_env%unit_nr
3147 IF (u > 0) THEN
3148 WRITE (u, fmt="(T2,A,ES34.1)") "Min. abs. eigenvalue of RI metric matrix M(k)", min_ev
3149 WRITE (u, fmt="(T2,A,ES34.1)") "Max. abs. eigenvalue of RI metric matrix M(k)", max_ev
3150 WRITE (u, fmt="(T2,A,ES50.1)") "Max. condition number of M(k)", cond_nr_max
3151 END IF
3152
3153 CALL timestop(handle)
3154
3155 END SUBROUTINE heuristic_ri_regularization
3156
3157! **************************************************************************************************
3158!> \brief ...
3159!> \param V_tr_R ...
3160!> \param pot_type ...
3161!> \param regularization_RI ...
3162!> \param bs_env ...
3163!> \param qs_env ...
3164! **************************************************************************************************
3165 SUBROUTINE get_v_tr_r(V_tr_R, pot_type, regularization_RI, bs_env, qs_env)
3166 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: v_tr_r
3167 TYPE(libint_potential_type) :: pot_type
3168 REAL(kind=dp) :: regularization_ri
3169 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3170 TYPE(qs_environment_type), POINTER :: qs_env
3171
3172 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_V_tr_R'
3173
3174 INTEGER :: handle, img, nimages_scf_desymm
3175 INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_ri
3176 INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
3177 TYPE(cp_blacs_env_type), POINTER :: blacs_env
3178 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_v_tr_r
3179 TYPE(dbcsr_distribution_type) :: dbcsr_dist
3180 TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: mat_v_tr_r
3181 TYPE(distribution_2d_type), POINTER :: dist_2d
3182 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3183 POINTER :: sab_ri
3184 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3185 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3186
3187 CALL timeset(routinen, handle)
3188
3189 NULLIFY (sab_ri, dist_2d)
3190
3191 CALL get_qs_env(qs_env=qs_env, &
3192 blacs_env=blacs_env, &
3193 distribution_2d=dist_2d, &
3194 qs_kind_set=qs_kind_set, &
3195 particle_set=particle_set)
3196
3197 ALLOCATE (sizes_ri(bs_env%n_atom))
3198 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=bs_env%basis_set_RI)
3199 CALL build_2c_neighbor_lists(sab_ri, bs_env%basis_set_RI, bs_env%basis_set_RI, &
3200 pot_type, "2c_nl_RI", qs_env, sym_ij=.false., &
3201 dist_2d=dist_2d)
3202 CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
3203 ALLOCATE (row_bsize(SIZE(sizes_ri)))
3204 ALLOCATE (col_bsize(SIZE(sizes_ri)))
3205 row_bsize(:) = sizes_ri
3206 col_bsize(:) = sizes_ri
3207
3208 nimages_scf_desymm = bs_env%nimages_scf_desymm
3209 ALLOCATE (mat_v_tr_r(nimages_scf_desymm))
3210 CALL dbcsr_create(mat_v_tr_r(1), "(RI|RI)", dbcsr_dist, dbcsr_type_no_symmetry, &
3211 row_bsize, col_bsize)
3212 DEALLOCATE (row_bsize, col_bsize)
3213
3214 DO img = 2, nimages_scf_desymm
3215 CALL dbcsr_create(mat_v_tr_r(img), template=mat_v_tr_r(1))
3216 END DO
3217
3218 CALL build_2c_integrals(mat_v_tr_r, 0.0_dp, qs_env, sab_ri, bs_env%basis_set_RI, &
3219 bs_env%basis_set_RI, pot_type, do_kpoints=.true., &
3220 ext_kpoints=bs_env%kpoints_scf_desymm, &
3221 regularization_ri=regularization_ri)
3222
3223 ALLOCATE (fm_v_tr_r(nimages_scf_desymm))
3224 DO img = 1, nimages_scf_desymm
3225 CALL cp_fm_create(fm_v_tr_r(img), bs_env%fm_RI_RI%matrix_struct)
3226 CALL copy_dbcsr_to_fm(mat_v_tr_r(img), fm_v_tr_r(img))
3227 CALL dbcsr_release(mat_v_tr_r(img))
3228 END DO
3229
3230 IF (.NOT. ALLOCATED(v_tr_r)) THEN
3231 ALLOCATE (v_tr_r(bs_env%n_RI, bs_env%n_RI, nimages_scf_desymm))
3232 END IF
3233
3234 CALL fm_to_local_array(fm_v_tr_r, v_tr_r)
3235
3236 CALL cp_fm_release(fm_v_tr_r)
3237 CALL dbcsr_distribution_release(dbcsr_dist)
3238 CALL release_neighbor_list_sets(sab_ri)
3239
3240 CALL timestop(handle)
3241
3242 END SUBROUTINE get_v_tr_r
3243
3244! **************************************************************************************************
3245!> \brief ...
3246!> \param matrix ...
3247!> \param exponent ...
3248!> \param eps ...
3249!> \param cond_nr ...
3250!> \param min_ev ...
3251!> \param max_ev ...
3252! **************************************************************************************************
3253 SUBROUTINE power(matrix, exponent, eps, cond_nr, min_ev, max_ev)
3254 COMPLEX(KIND=dp), DIMENSION(:, :) :: matrix
3255 REAL(kind=dp) :: exponent, eps
3256 REAL(kind=dp), OPTIONAL :: cond_nr, min_ev, max_ev
3257
3258 CHARACTER(len=*), PARAMETER :: routinen = 'power'
3259
3260 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors
3261 INTEGER :: handle, i, n
3262 REAL(kind=dp) :: pos_eval
3263 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
3264
3265 CALL timeset(routinen, handle)
3266
3267 ! make matrix perfectly Hermitian
3268 matrix(:, :) = 0.5_dp*(matrix(:, :) + conjg(transpose(matrix(:, :))))
3269
3270 n = SIZE(matrix, 1)
3271 ALLOCATE (eigenvalues(n), eigenvectors(n, n))
3272 CALL diag_complex(matrix, eigenvectors, eigenvalues)
3273
3274 IF (PRESENT(cond_nr)) cond_nr = maxval(abs(eigenvalues))/minval(abs(eigenvalues))
3275 IF (PRESENT(min_ev)) min_ev = minval(abs(eigenvalues))
3276 IF (PRESENT(max_ev)) max_ev = maxval(abs(eigenvalues))
3277
3278 DO i = 1, n
3279 IF (eps < eigenvalues(i)) THEN
3280 pos_eval = (eigenvalues(i))**(0.5_dp*exponent)
3281 ELSE
3282 pos_eval = 0.0_dp
3283 END IF
3284 eigenvectors(:, i) = eigenvectors(:, i)*pos_eval
3285 END DO
3286
3287 CALL zgemm("N", "C", n, n, n, z_one, eigenvectors, n, eigenvectors, n, z_zero, matrix, n)
3288
3289 DEALLOCATE (eigenvalues, eigenvectors)
3290
3291 CALL timestop(handle)
3292
3293 END SUBROUTINE power
3294
3295! **************************************************************************************************
3296!> \brief ...
3297!> \param bs_env ...
3298!> \param Sigma_c_n_time ...
3299!> \param Sigma_c_n_freq ...
3300!> \param ispin ...
3301! **************************************************************************************************
3302 SUBROUTINE time_to_freq(bs_env, Sigma_c_n_time, Sigma_c_n_freq, ispin)
3303 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3304 REAL(kind=dp), DIMENSION(:, :, :) :: sigma_c_n_time, sigma_c_n_freq
3305 INTEGER :: ispin
3306
3307 CHARACTER(LEN=*), PARAMETER :: routinen = 'time_to_freq'
3308
3309 INTEGER :: handle, i_t, j_w, n_occ
3310 REAL(kind=dp) :: freq_j, time_i, w_cos_ij, w_sin_ij
3311 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sigma_c_n_cos_time, sigma_c_n_sin_time
3312
3313 CALL timeset(routinen, handle)
3314
3315 ALLOCATE (sigma_c_n_cos_time(bs_env%n_ao, bs_env%num_time_freq_points))
3316 ALLOCATE (sigma_c_n_sin_time(bs_env%n_ao, bs_env%num_time_freq_points))
3317
3318 sigma_c_n_cos_time(:, :) = 0.5_dp*(sigma_c_n_time(:, :, 1) + sigma_c_n_time(:, :, 2))
3319 sigma_c_n_sin_time(:, :) = 0.5_dp*(sigma_c_n_time(:, :, 1) - sigma_c_n_time(:, :, 2))
3320
3321 sigma_c_n_freq(:, :, :) = 0.0_dp
3322
3323 DO i_t = 1, bs_env%num_time_freq_points
3324
3325 DO j_w = 1, bs_env%num_time_freq_points
3326
3327 freq_j = bs_env%imag_freq_points(j_w)
3328 time_i = bs_env%imag_time_points(i_t)
3329 ! integration weights for cosine and sine transform
3330 w_cos_ij = bs_env%weights_cos_t_to_w(j_w, i_t)*cos(freq_j*time_i)
3331 w_sin_ij = bs_env%weights_sin_t_to_w(j_w, i_t)*sin(freq_j*time_i)
3332
3333 ! 1. Re(Σ^c_nn(k_i,iω)) from cosine transform
3334 sigma_c_n_freq(:, j_w, 1) = sigma_c_n_freq(:, j_w, 1) + &
3335 w_cos_ij*sigma_c_n_cos_time(:, i_t)
3336
3337 ! 2. Im(Σ^c_nn(k_i,iω)) from sine transform
3338 sigma_c_n_freq(:, j_w, 2) = sigma_c_n_freq(:, j_w, 2) + &
3339 w_sin_ij*sigma_c_n_sin_time(:, i_t)
3340
3341 END DO
3342
3343 END DO
3344
3345 ! for occupied levels, we need the correlation self-energy for negative omega.
3346 ! Therefore, weight_sin should be computed with -omega, which results in an
3347 ! additional minus for the imaginary part:
3348 n_occ = bs_env%n_occ(ispin)
3349 sigma_c_n_freq(1:n_occ, :, 2) = -sigma_c_n_freq(1:n_occ, :, 2)
3350
3351 CALL timestop(handle)
3352
3353 END SUBROUTINE time_to_freq
3354
3355! **************************************************************************************************
3356!> \brief ...
3357!> \param bs_env ...
3358!> \param Sigma_c_ikp_n_freq ...
3359!> \param Sigma_x_ikp_n ...
3360!> \param V_xc_ikp_n ...
3361!> \param eigenval_scf ...
3362!> \param ikp ...
3363!> \param ispin ...
3364! **************************************************************************************************
3365 SUBROUTINE analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_xc_ikp_n, &
3366 eigenval_scf, ikp, ispin)
3367
3368 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3369 REAL(kind=dp), DIMENSION(:, :, :) :: sigma_c_ikp_n_freq
3370 REAL(kind=dp), DIMENSION(:) :: sigma_x_ikp_n, v_xc_ikp_n, eigenval_scf
3371 INTEGER :: ikp, ispin
3372
3373 CHARACTER(LEN=*), PARAMETER :: routinen = 'analyt_conti_and_print'
3374
3375 CHARACTER(len=3) :: occ_vir
3376 CHARACTER(len=default_path_length) :: fname
3377 INTEGER :: handle, i_mo, ikp_for_print, iunit, &
3378 n_mo, nkp
3379 LOGICAL :: is_bandstruc_kpoint, print_dos_kpoints, &
3380 print_ikp
3381 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dummy, sigma_c_ikp_n_qp
3382
3383 CALL timeset(routinen, handle)
3384
3385 n_mo = bs_env%n_ao
3386 ALLOCATE (dummy(n_mo), sigma_c_ikp_n_qp(n_mo))
3387 sigma_c_ikp_n_qp(:) = 0.0_dp
3388
3389 DO i_mo = 1, n_mo
3390
3391 ! parallelization
3392 IF (modulo(i_mo, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
3393
3394 CALL continuation_pade(sigma_c_ikp_n_qp, &
3395 bs_env%imag_freq_points_fit, dummy, dummy, &
3396 sigma_c_ikp_n_freq(:, 1:bs_env%num_freq_points_fit, 1)*z_one + &
3397 sigma_c_ikp_n_freq(:, 1:bs_env%num_freq_points_fit, 2)*gaussi, &
3398 sigma_x_ikp_n(:) - v_xc_ikp_n(:), &
3399 eigenval_scf(:), eigenval_scf(:), &
3400 bs_env%do_hedin_shift, &
3401 i_mo, bs_env%n_occ(ispin), bs_env%n_vir(ispin), &
3402 bs_env%nparam_pade, bs_env%num_freq_points_fit, &
3403 ri_rpa_g0w0_crossing_newton, bs_env%n_occ(ispin), &
3404 0.0_dp, .true., .false., 1, e_fermi_ext=bs_env%e_fermi(ispin))
3405 END DO
3406
3407 CALL bs_env%para_env%sum(sigma_c_ikp_n_qp)
3408
3409 CALL correct_obvious_fitting_fails(sigma_c_ikp_n_qp, ispin, bs_env)
3410
3411 bs_env%eigenval_G0W0(:, ikp, ispin) = eigenval_scf(:) + &
3412 sigma_c_ikp_n_qp(:) + &
3413 sigma_x_ikp_n(:) - &
3414 v_xc_ikp_n(:)
3415
3416 bs_env%eigenval_HF(:, ikp, ispin) = eigenval_scf(:) + sigma_x_ikp_n(:) - v_xc_ikp_n(:)
3417
3418 ! only print eigenvalues of DOS k-points in case no bandstructure path has been given
3419 print_dos_kpoints = (bs_env%nkp_only_bs <= 0)
3420 ! in kpoints_DOS, the last nkp_only_bs are bandstructure k-points
3421 is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
3422 print_ikp = print_dos_kpoints .OR. is_bandstruc_kpoint
3423
3424 IF (bs_env%para_env%is_source() .AND. print_ikp) THEN
3425
3426 IF (print_dos_kpoints) THEN
3427 nkp = bs_env%nkp_only_DOS
3428 ikp_for_print = ikp
3429 ELSE
3430 nkp = bs_env%nkp_only_bs
3431 ikp_for_print = ikp - bs_env%nkp_only_DOS
3432 END IF
3433
3434 fname = "bandstructure_SCF_and_G0W0"
3435
3436 IF (ikp_for_print == 1 .AND. ispin == 1) THEN
3437 CALL open_file(trim(fname), unit_number=iunit, file_status="REPLACE", &
3438 file_action="WRITE")
3439 ELSE
3440 CALL open_file(trim(fname), unit_number=iunit, file_status="OLD", &
3441 file_action="WRITE", file_position="APPEND")
3442 END IF
3443
3444 WRITE (iunit, "(A)") " "
3445 WRITE (iunit, "(A10,I7,A25,3F10.4,T90,A7,I2)") "kpoint: ", ikp_for_print, "coordinate: ", &
3446 bs_env%kpoints_DOS%xkp(:, ikp), "spin: ", ispin
3447 WRITE (iunit, "(A)") " "
3448 WRITE (iunit, "(A5,A12,3A17,A16,A18)") "n", "k", ϵ"_nk^DFT (eV)", Σ"^c_nk (eV)", &
3449 Σ"^x_nk (eV)", "v_nk^xc (eV)", ϵ"_nk^G0W0 (eV)"
3450 WRITE (iunit, "(A)") " "
3451
3452 DO i_mo = 1, n_mo
3453 IF (i_mo <= bs_env%n_occ(ispin)) occ_vir = 'occ'
3454 IF (i_mo > bs_env%n_occ(ispin)) occ_vir = 'vir'
3455 WRITE (iunit, "(I5,3A,I5,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', ikp_for_print, &
3456 eigenval_scf(i_mo)*evolt, &
3457 sigma_c_ikp_n_qp(i_mo)*evolt, &
3458 sigma_x_ikp_n(i_mo)*evolt, &
3459 v_xc_ikp_n(i_mo)*evolt, &
3460 bs_env%eigenval_G0W0(i_mo, ikp, ispin)*evolt
3461 END DO
3462
3463 WRITE (iunit, "(A)") " "
3464
3465 CALL close_file(iunit)
3466
3467 END IF
3468
3469 CALL timestop(handle)
3470
3471 END SUBROUTINE analyt_conti_and_print
3472
3473! **************************************************************************************************
3474!> \brief ...
3475!> \param Sigma_c_ikp_n_qp ...
3476!> \param ispin ...
3477!> \param bs_env ...
3478! **************************************************************************************************
3479 SUBROUTINE correct_obvious_fitting_fails(Sigma_c_ikp_n_qp, ispin, bs_env)
3480 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sigma_c_ikp_n_qp
3481 INTEGER :: ispin
3482 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3483
3484 CHARACTER(LEN=*), PARAMETER :: routinen = 'correct_obvious_fitting_fails'
3485
3486 INTEGER :: handle, homo, i_mo, j_mo, &
3487 n_levels_scissor, n_mo
3488 LOGICAL :: is_occ, is_vir
3489 REAL(kind=dp) :: sum_sigma_c
3490
3491 CALL timeset(routinen, handle)
3492
3493 n_mo = bs_env%n_ao
3494 homo = bs_env%n_occ(ispin)
3495
3496 DO i_mo = 1, n_mo
3497
3498 ! if |𝚺^c| > 13 eV, we use a scissors shift
3499 IF (abs(sigma_c_ikp_n_qp(i_mo)) > 13.0_dp/evolt) THEN
3500
3501 is_occ = (i_mo <= homo)
3502 is_vir = (i_mo > homo)
3503
3504 n_levels_scissor = 0
3505 sum_sigma_c = 0.0_dp
3506
3507 ! compute scissor
3508 DO j_mo = 1, n_mo
3509
3510 ! only compute scissor from other GW levels close in energy
3511 IF (is_occ .AND. j_mo > homo) cycle
3512 IF (is_vir .AND. j_mo <= homo) cycle
3513 IF (abs(i_mo - j_mo) > 10) cycle
3514 IF (i_mo == j_mo) cycle
3515
3516 n_levels_scissor = n_levels_scissor + 1
3517 sum_sigma_c = sum_sigma_c + sigma_c_ikp_n_qp(j_mo)
3518
3519 END DO
3520
3521 ! overwrite the self-energy with scissor shift
3522 sigma_c_ikp_n_qp(i_mo) = sum_sigma_c/real(n_levels_scissor, kind=dp)
3523
3524 END IF
3525
3526 END DO ! i_mo
3527
3528 CALL timestop(handle)
3529
3530 END SUBROUTINE correct_obvious_fitting_fails
3531
3532END MODULE gw_utils
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
tensor
struct tensor_ tensor
cell_types::pbc
Definition cell_types.F:104
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:56
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::graml2024
integer, save, public graml2024
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::scaled_to_real
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:625
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:792
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:643
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_dist2d_to_dist
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
Definition cp_dbcsr_operations.F:415
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:564
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition cp_output_handling.F:754
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition distribution_2d_types.F:15
gw_communication
Definition gw_communication.F:13
gw_communication::fm_to_local_array
subroutine, public fm_to_local_array(fm_s, array_s, weight, add)
...
Definition gw_communication.F:764
gw_integrals
Utility method to build 3-center integrals for small cell GW.
Definition gw_integrals.F:11
gw_integrals::build_3c_integral_block
subroutine, public build_3c_integral_block(int_3c, qs_env, potential_parameter, basis_j, basis_k, basis_i, cell_j, cell_k, cell_i, atom_j, atom_k, atom_i, j_bf_start_from_atom, k_bf_start_from_atom, i_bf_start_from_atom)
...
Definition gw_integrals.F:461
gw_utils
Definition gw_utils.F:15
gw_utils::get_v_tr_r
subroutine, public get_v_tr_r(v_tr_r, pot_type, regularization_ri, bs_env, qs_env)
...
Definition gw_utils.F:3166
gw_utils::time_to_freq
subroutine, public time_to_freq(bs_env, sigma_c_n_time, sigma_c_n_freq, ispin)
...
Definition gw_utils.F:3303
gw_utils::rtbse_resolve_rirs_flag
subroutine, public rtbse_resolve_rirs_flag(qs_env, bs_env, rirs_kernel)
Resolve the linRTBSE RI-RS kernel switch from the KERNEL_RI input and the GW default.
Definition gw_utils.F:285
gw_utils::compute_xkp
subroutine, public compute_xkp(xkp, ikp_start, ikp_end, grid)
...
Definition gw_utils.F:673
gw_utils::analyt_conti_and_print
subroutine, public analyt_conti_and_print(bs_env, sigma_c_ikp_n_freq, sigma_x_ikp_n, v_xc_ikp_n, eigenval_scf, ikp, ispin)
...
Definition gw_utils.F:3367
gw_utils::create_and_init_bs_env_for_gw
subroutine, public create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec)
...
Definition gw_utils.F:152
gw_utils::add_r
subroutine, public add_r(cell_1, cell_2, index_to_cell, cell_1_plus_2, cell_found, cell_to_index, i_cell_1_plus_2)
...
Definition gw_utils.F:2881
gw_utils::de_init_bs_env
subroutine, public de_init_bs_env(qs_env, bs_env)
Releases the memory-heavy GW intermediates that cannot be freed in bs_env_release,...
Definition gw_utils.F:240
gw_utils::get_i_j_atoms
subroutine, public get_i_j_atoms(atoms_i, atoms_j, n_atom_i, n_atom_j, color_sub, bs_env)
...
Definition gw_utils.F:1564
gw_utils::power
subroutine, public power(matrix, exponent, eps, cond_nr, min_ev, max_ev)
...
Definition gw_utils.F:3254
gw_utils::is_cell_in_index_to_cell
subroutine, public is_cell_in_index_to_cell(cell, index_to_cell, cell_found)
...
Definition gw_utils.F:2920
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::rtp_method_bse_linearized
integer, parameter, public rtp_method_bse_linearized
Definition input_constants.F:986
input_constants::non_periodic_ri_rs
integer, parameter, public non_periodic_ri_rs
Definition input_constants.F:1124
input_constants::rtp_bse_kernel_ri_rs
integer, parameter, public rtp_bse_kernel_ri_rs
Definition input_constants.F:993
input_constants::rtp_bse_kernel_ri_ao
integer, parameter, public rtp_bse_kernel_ri_ao
Definition input_constants.F:993
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:837
input_constants::rtp_method_bse
integer, parameter, public rtp_method_bse
Definition input_constants.F:986
input_constants::small_cell_full_kp
integer, parameter, public small_cell_full_kp
Definition input_constants.F:1124
input_constants::large_cell_gamma_ri_rs
integer, parameter, public large_cell_gamma_ri_rs
Definition input_constants.F:1124
input_constants::large_cell_gamma
integer, parameter, public large_cell_gamma
Definition input_constants.F:1124
input_constants::xc_none
integer, parameter, public xc_none
Definition input_constants.F:575
input_constants::ri_rpa_g0w0_crossing_newton
integer, parameter, public ri_rpa_g0w0_crossing_newton
Definition input_constants.F:1124
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition input_section_types.F:1233
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kpoint_k_r_trafo_simple
Implements transformations from k-space to R-space for Fortran array matrices.
Definition kpoint_k_r_trafo_simple.F:16
kpoint_k_r_trafo_simple::rs_to_kp
subroutine, public rs_to_kp(rs_real, ks_complex, index_to_cell, xkp, deriv_direction, hmat)
Integrate RS matrices (stored as Fortran array) into a kpoint matrix at given kp.
Definition kpoint_k_r_trafo_simple.F:230
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:529
kpoint_types::kpoint_create
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
Definition kpoint_types.F:237
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
libint_wrapper::cp_libint_static_cleanup
subroutine, public cp_libint_static_cleanup()
Definition libint_wrapper.F:1351
libint_wrapper::cp_libint_static_init
subroutine, public cp_libint_static_init()
Definition libint_wrapper.F:1343
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:440
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diag_complex
subroutine, public diag_complex(matrix, eigenvectors, eigenvalues)
Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd.
Definition mathlib.F:1878
mathlib::gcd
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1289
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
minimax_exp_gw
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition minimax_exp_gw.F:17
minimax_exp_gw::get_exp_minimax_coeff_gw
subroutine, public get_exp_minimax_coeff_gw(k, e_range, aw)
...
Definition minimax_exp_gw.F:38
minimax_exp
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition minimax_exp.F:29
minimax_exp::get_exp_minimax_coeff
subroutine, public get_exp_minimax_coeff(k, rc, aw, mm_error, which_coeffs)
Get best minimax approximation for given input parameters. Automatically chooses the most exact set o...
Definition minimax_exp.F:127
minimax_rpa
Routines to calculate the minimax coefficients for approximating 1/x as 1/x ~ 1/pi SUM_{i}^{K} w_i x^...
Definition minimax_rpa.F:14
minimax_rpa::get_rpa_minimax_coeff_larger_grid
subroutine, public get_rpa_minimax_coeff_larger_grid(k, e_range, aw)
...
Definition minimax_rpa.F:10119
minimax_rpa::get_rpa_minimax_coeff
subroutine, public get_rpa_minimax_coeff(k, e_range, aw, ierr, print_warning)
The a_i and w_i coefficient are stored in aw such that the first 1:K elements correspond to a_i and t...
Definition minimax_rpa.F:41
mp2_gpw
Calls routines to get RI integrals and calculate total energies.
Definition mp2_gpw.F:14
mp2_gpw::create_mat_munu
subroutine, public create_mat_munu(mat_munu, qs_env, eps_grid, blacs_env_sub, do_ri_aux_basis, do_mixed_basis, group_size_prim, do_alloc_blocks_from_nbl, do_kpoints, sab_orb_sub, dbcsr_sym_type)
Encapsulate the building of dbcsr_matrix mat_munu.
Definition mp2_gpw.F:989
mp2_grids
Routines to calculate frequency and time grids (integration points and weights) for correlation metho...
Definition mp2_grids.F:14
mp2_grids::get_l_sq_wghts_cos_tf_w_to_t
subroutine, public get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, omega_tj, e_min, e_max, max_error, num_points_per_magnitude, regularization)
...
Definition mp2_grids.F:1235
mp2_grids::get_l_sq_wghts_cos_tf_t_to_w
subroutine, public get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, omega_tj, e_min, e_max, max_error, num_points_per_magnitude, regularization)
Calculate integration weights for the tau grid (in dependency of the omega node)
Definition mp2_grids.F:736
mp2_grids::get_l_sq_wghts_sin_tf_t_to_w
subroutine, public get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, omega_tj, e_min, e_max, max_error, num_points_per_magnitude, regularization)
Calculate integration weights for the tau grid (in dependency of the omega node)
Definition mp2_grids.F:872
mp2_ri_2c
Framework for 2c-integrals for RI.
Definition mp2_ri_2c.F:14
mp2_ri_2c::trunc_coulomb_for_exchange
subroutine, public trunc_coulomb_for_exchange(qs_env, trunc_coulomb, rel_cutoff_trunc_coulomb_ri_x, cell_grid, do_bvk_cell)
...
Definition mp2_ri_2c.F:1601
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:121
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
post_scf_bandstructure_types
Definition post_scf_bandstructure_types.F:13
post_scf_bandstructure_utils
Definition post_scf_bandstructure_utils.F:13
post_scf_bandstructure_utils::rsmat_to_kp
subroutine, public rsmat_to_kp(mat_rs, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_kp, imag_rs_mat)
...
Definition post_scf_bandstructure_utils.F:717
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_environment_types::qs_env_part_release
subroutine, public qs_env_part_release(qs_env)
releases part of the given qs_env in order to save memory
Definition qs_environment_types.F:1713
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_interactions
Calculate the interaction radii for the operator matrix calculation.
Definition qs_interactions.F:17
qs_interactions::init_interaction_radii_orb_basis
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Definition qs_interactions.F:439
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_methods.F:22
qs_ks_methods::qs_ks_build_kohn_sham_matrix
subroutine, public qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, print_active, ext_ks_matrix, ext_xc_section)
routine where the real calculations are made: the KS matrix is calculated
Definition qs_ks_methods.F:201
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition qs_neighbor_list_types.F:1096
qs_tensors_types
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors_types.F:12
qs_tensors_types::distribution_3d_create
subroutine, public distribution_3d_create(dist_3d, dist1, dist2, dist3, nkind, particle_set, mp_comm_3d, own_comm)
Create a 3d distribution.
Definition qs_tensors_types.F:91
qs_tensors_types::create_2c_tensor
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
Definition qs_tensors_types.F:382
qs_tensors_types::create_3c_tensor
subroutine, public create_3c_tensor(t3c, dist_1, dist_2, dist_3, pgrid, sizes_1, sizes_2, sizes_3, map1, map2, name)
...
Definition qs_tensors_types.F:335
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors.F:11
qs_tensors::build_2c_integrals
subroutine, public build_2c_integrals(t2c, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints, do_hfx_kpoints, ext_kpoints, regularization_ri)
...
Definition qs_tensors.F:3156
qs_tensors::build_2c_neighbor_lists
subroutine, public build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, sym_ij, molecular, dist_2d, pot_to_rad)
Build 2-center neighborlists adapted to different operators This mainly wraps build_neighbor_lists fo...
Definition qs_tensors.F:144
qs_tensors::build_3c_integrals
subroutine, public build_3c_integrals(t3c, filter_eps, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, int_eps, op_pos, do_kpoints, do_hfx_kpoints, desymmetrize, cell_sym, bounds_i, bounds_j, bounds_k, ri_range, img_to_ri_cell, cell_to_index_ext)
Build 3-center integral tensor.
Definition qs_tensors.F:2122
qs_tensors::neighbor_list_3c_destroy
subroutine, public neighbor_list_3c_destroy(ijk_list)
Destroy 3c neighborlist.
Definition qs_tensors.F:381
qs_tensors::get_tensor_occupancy
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Definition qs_tensors.F:3624
qs_tensors::build_3c_neighbor_lists
subroutine, public build_3c_neighbor_lists(ijk_list, basis_i, basis_j, basis_k, dist_3d, potential_parameter, name, qs_env, sym_ij, sym_jk, sym_ik, molecular, op_pos, own_dist)
Build a 3-center neighbor list.
Definition qs_tensors.F:280
rpa_gw
Routines for GW, continuous development [Jan Wilhelm].
Definition rpa_gw.F:14
rpa_gw::continuation_pade
subroutine, public continuation_pade(vec_gw_energ, vec_omega_fit_gw, z_value, m_value, vec_sigma_c_gw, vec_sigma_x_minus_vxc_gw, eigenval, eigenval_scf, do_hedin_shift, n_level_gw, gw_corr_lev_occ, gw_corr_lev_vir, nparam_pade, num_fit_points, crossing_search, homo, fermi_level_offset, do_gw_im_time, print_self_energy, count_ev_sc_gw, vec_gw_dos, dos_lower_bound, dos_precision, ndos, min_level_self_energy, max_level_self_energy, dos_eta, dos_min, dos_max, e_fermi_ext)
perform analytic continuation with pade approximation
Definition rpa_gw.F:4314
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_types.F:62
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:175
libint_2c_3c::libint_potential_type
Definition libint_2c_3c.F:67
message_passing::mp_cart_type
Definition message_passing.F:727
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
particle_types::particle_type
Definition particle_types.F:35
post_scf_bandstructure_types::post_scf_bandstructure_type
Definition post_scf_bandstructure_types.F:125
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_tensors_types::distribution_3d_type
Definition qs_tensors_types.F:57
qs_tensors_types::neighbor_list_3c_type
Definition qs_tensors_types.F:63