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