(git:5e7fe52)
Loading...
Searching...
No Matches
gw_large_cell_gamma.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 Routines from paper [Graml2024]
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_large_cell_gamma
16 USE atomic_kind_types, ONLY: atomic_kind_type
17 USE bibliography, ONLY: graml2024,&
18 cite_reference
19 USE cell_types, ONLY: cell_type,&
20 get_cell,&
21 pbc
22 USE constants_operator, ONLY: operator_coulomb
23 USE cp_cfm_basic_linalg, ONLY: cp_cfm_uplo_to_full
24 USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
25 cp_cfm_cholesky_invert
26 USE cp_cfm_diag, ONLY: cp_cfm_geeig
27 USE cp_cfm_types, ONLY: cp_cfm_create,&
28 cp_cfm_get_info,&
29 cp_cfm_release,&
30 cp_cfm_to_cfm,&
31 cp_cfm_to_fm,&
32 cp_cfm_type,&
33 cp_fm_to_cfm
34 USE cp_dbcsr_api, ONLY: &
35 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_block_p, &
36 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
37 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_release, dbcsr_set, &
38 dbcsr_type
39 USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
40 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
41 copy_fm_to_dbcsr,&
42 dbcsr_deallocate_matrix_set
43 USE cp_files, ONLY: close_file,&
44 open_file
45 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
46 USE cp_fm_diag, ONLY: cp_fm_power
47 USE cp_fm_types, ONLY: &
48 cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_read_unformatted, cp_fm_release, &
49 cp_fm_set_all, cp_fm_to_fm, cp_fm_type, cp_fm_write_unformatted
50 USE cp_log_handling, ONLY: cp_get_default_logger,&
51 cp_logger_type
52 USE cp_output_handling, ONLY: cp_p_file,&
53 cp_print_key_should_output,&
54 cp_print_key_unit_nr
55 USE dbt_api, ONLY: dbt_clear,&
56 dbt_contract,&
57 dbt_copy,&
58 dbt_create,&
59 dbt_destroy,&
60 dbt_filter,&
61 dbt_type
62 USE gw_communication, ONLY: fm_to_local_tensor,&
63 local_dbt_to_global_mat
64 USE gw_utils, ONLY: analyt_conti_and_print,&
65 de_init_bs_env,&
66 time_to_freq
67 USE input_constants, ONLY: g0w0,&
68 rtp_method_bse,&
69 rtp_method_bse_linearized
70 USE input_section_types, ONLY: section_vals_type
71 USE kinds, ONLY: default_path_length,&
72 dp,&
73 int_8
74 USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
75 USE kpoint_types, ONLY: kpoint_type
76 USE machine, ONLY: m_walltime
77 USE mathconstants, ONLY: twopi,&
78 z_one,&
79 z_zero
80 USE message_passing, ONLY: mp_file_delete
81 USE mp2_ri_2c, ONLY: ri_2c_integral_mat
82 USE parallel_gemm_api, ONLY: parallel_gemm
83 USE particle_types, ONLY: particle_type
84 USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
85 USE post_scf_bandstructure_utils, ONLY: mic_contribution_from_ikp,&
86 cfm_ikp_from_fm_gamma,&
87 get_all_vbm_cbm_bandgaps
88 USE qs_environment_types, ONLY: get_qs_env,&
89 qs_environment_type
90 USE qs_kind_types, ONLY: qs_kind_type
91 USE qs_tensors, ONLY: build_3c_integrals
92 USE rpa_gw_kpoints_util, ONLY: cp_cfm_power
93#include "./base/base_uses.f90"
94
95 IMPLICIT NONE
96
97 PRIVATE
98
99 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_large_cell_gamma'
100
101 PUBLIC :: gw_calc_large_cell_gamma, &
102 compute_3c_integrals, g_occ_vir, fm_read, write_matrix, &
103 fill_fm_sigma_c_gamma_time, delete_unnecessary_files, &
104 multiply_fm_w_mic_time_with_minv_gamma, get_w_mic, &
105 create_fm_w_mic_time, fourier_transform_w_to_t, &
106 compute_fm_chi_gamma_freq, compute_qp_energies, fm_write
107
108CONTAINS
109
110! **************************************************************************************************
111!> \brief Perform GW band structure calculation
112!> \param qs_env ...
113!> \param bs_env ...
114!> \par History
115!> * 07.2023 created [Jan Wilhelm]
116! **************************************************************************************************
117 SUBROUTINE gw_calc_large_cell_gamma(qs_env, bs_env)
118 TYPE(qs_environment_type), POINTER :: qs_env
119 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
120
121 CHARACTER(LEN=*), PARAMETER :: routinen = 'gw_calc_large_cell_Gamma'
122
123 INTEGER :: handle
124 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_sigma_x_gamma, fm_w_mic_time
125 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_sigma_c_gamma_time
126
127 CALL timeset(routinen, handle)
128
129 CALL cite_reference(graml2024)
130
131 ! G^occ_µλ(i|τ|,k=0) = sum_n^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
132 ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
133 ! χ_PQ(iτ,k=0) = sum_λν [sum_µ (µν|P) G^occ_µλ(i|τ|)] [sum_σ (σλ|Q) G^vir_σν(i|τ|)]
134 CALL get_mat_chi_gamma_tau(bs_env, qs_env, bs_env%mat_chi_Gamma_tau)
135
136 ! χ_PQ(iτ,k=0) -> χ_PQ(iω,k) -> ε_PQ(iω,k) -> W_PQ(iω,k) -> W^MIC_PQ(iτ) -> M^-1*W^MIC*M^-1
137 CALL get_w_mic(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_w_mic_time)
138
139 ! D_µν = sum_n^occ C_µn(k=0) C_νn(k=0), V^trunc_PQ = sum_cell_R <phi_P,0|V^trunc|phi_Q,R>
140 ! Σ^x_λσ(k=0) = sum_νQ [sum_P (νσ|P) V^trunc_PQ] [sum_µ (λµ|Q) D_µν)]
141 CALL get_sigma_x(bs_env, qs_env, fm_sigma_x_gamma)
142
143 ! Σ^c_λσ(iτ,k=0) = sum_νQ [sum_P (νσ|P) W^MIC_PQ(iτ)] [sum_µ (λµ|Q) G^occ_µν(i|τ|)], τ < 0
144 ! Σ^c_λσ(iτ,k=0) = sum_νQ [sum_P (νσ|P) W^MIC_PQ(iτ)] [sum_µ (λµ|Q) G^vir_µν(i|τ|)], τ > 0
145 CALL get_sigma_c(bs_env, qs_env, fm_w_mic_time, fm_sigma_c_gamma_time)
146
147 ! Σ^c_λσ(iτ,k=0) -> Σ^c_nn(ϵ,k); ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ,k) + Σ^x_nn(k) - v^xc_nn(k)
148 CALL compute_qp_energies(bs_env, qs_env, fm_sigma_x_gamma, fm_sigma_c_gamma_time)
149
150 CALL de_init_bs_env(qs_env, bs_env)
151
152 CALL timestop(handle)
153
154 END SUBROUTINE gw_calc_large_cell_gamma
155
156! **************************************************************************************************
157!> \brief ...
158!> \param bs_env ...
159!> \param qs_env ...
160!> \param mat_chi_Gamma_tau ...
161! **************************************************************************************************
162 SUBROUTINE get_mat_chi_gamma_tau(bs_env, qs_env, mat_chi_Gamma_tau)
163 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
164 TYPE(qs_environment_type), POINTER :: qs_env
165 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
166
167 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_mat_chi_Gamma_tau'
168
169 INTEGER :: handle, i_intval_idx, i_t, inner_loop_atoms_interval_index, ispin, j_intval_idx
170 INTEGER(KIND=int_8) :: flop
171 INTEGER, DIMENSION(2) :: bounds_p, bounds_q, i_atoms, il_atoms, &
172 j_atoms
173 INTEGER, DIMENSION(2, 2) :: bounds_comb
174 LOGICAL :: dist_too_long_i, dist_too_long_j
175 REAL(kind=dp) :: t1, tau
176 TYPE(dbt_type) :: t_2c_gocc, t_2c_gvir, t_3c_for_gocc, &
177 t_3c_for_gvir, t_3c_x_gocc, &
178 t_3c_x_gocc_2, t_3c_x_gvir, &
179 t_3c_x_gvir_2
180
181 CALL timeset(routinen, handle)
182
183 DO i_t = 1, bs_env%num_time_freq_points
184
185 t1 = m_walltime()
186
187 IF (bs_env%read_chi(i_t)) THEN
188
189 CALL fm_read(bs_env%fm_RI_RI, bs_env, bs_env%chi_name, i_t)
190
191 CALL copy_fm_to_dbcsr(bs_env%fm_RI_RI, mat_chi_gamma_tau(i_t)%matrix, &
192 keep_sparsity=.false.)
193
194 IF (bs_env%unit_nr > 0) THEN
195 WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
196 'Read χ(iτ,k=0) from file for time point ', i_t, ' /', &
197 bs_env%num_time_freq_points, &
198 ', Execution time', m_walltime() - t1, ' s'
199 END IF
200
201 cycle
202
203 END IF
204
205 IF (.NOT. bs_env%calc_chi(i_t)) cycle
206
207 CALL create_tensors_chi(t_2c_gocc, t_2c_gvir, t_3c_for_gocc, t_3c_for_gvir, &
208 t_3c_x_gocc, t_3c_x_gvir, t_3c_x_gocc_2, t_3c_x_gvir_2, bs_env)
209
210 ! 1. compute G^occ and G^vir
211 ! Background: G^σ(iτ) = G^occ,σ(iτ) * Θ(-τ) + G^vir,σ(iτ) * Θ(τ), σ ∈ {↑,↓}
212 ! G^occ,σ_µλ(i|τ|,k=0) = sum_n^occ C^σ_µn(k=0) e^(-|(ϵ^σ_nk=0-ϵ_F)τ|) C^σ_λn(k=0)
213 ! G^vir,σ_µλ(i|τ|,k=0) = sum_n^vir C^σ_µn(k=0) e^(-|(ϵ^σ_nk=0-ϵ_F)τ|) C^σ_λn(k=0)
214 tau = bs_env%imag_time_points(i_t)
215
216 DO ispin = 1, bs_env%n_spin
217 CALL g_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.true., vir=.false.)
218 CALL g_occ_vir(bs_env, tau, bs_env%fm_Gvir, ispin, occ=.false., vir=.true.)
219
220 CALL fm_to_local_tensor(bs_env%fm_Gocc, bs_env%mat_ao_ao%matrix, &
221 bs_env%mat_ao_ao_tensor%matrix, t_2c_gocc, bs_env, &
222 bs_env%atoms_j_t_group)
223 CALL fm_to_local_tensor(bs_env%fm_Gvir, bs_env%mat_ao_ao%matrix, &
224 bs_env%mat_ao_ao_tensor%matrix, t_2c_gvir, bs_env, &
225 bs_env%atoms_i_t_group)
226
227 ! every group has its own range of i_atoms and j_atoms; only deal with a
228 ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
229 DO i_intval_idx = 1, bs_env%n_intervals_i
230 DO j_intval_idx = 1, bs_env%n_intervals_j
231 i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
232 j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
233
234 IF (bs_env%skip_chi(i_intval_idx, j_intval_idx)) THEN
235 ! Do that only after first timestep to avoid skips due to vanishing G
236 ! caused by gaps
237 IF (i_t == 2) THEN
238 bs_env%n_skip_chi = bs_env%n_skip_chi + 1
239 END IF
240 cycle
241 END IF
242
243 DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
244
245 il_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
246 ! Idea: Use sparsity in 3c integrals behind χ_PQ(iτ,k=0)
247 ! -> λ bounds from j_atoms -> sparse in IL_atoms through σ in
248 ! N_Qλν(iτ) = sum_σ (Qλ|σ) G^vir_νσ(i|τ|,k=0)
249 ! -> ν bounds from i_atoms -> sparse in IL_atoms through µ in
250 ! M_Pνλ(iτ) = sum_µ (Pν|µ) G^occ_λµ(i|τ|,k=0)
251 CALL check_dist(i_atoms, il_atoms, qs_env, bs_env, dist_too_long_i)
252 CALL check_dist(j_atoms, il_atoms, qs_env, bs_env, dist_too_long_j)
253 IF (.NOT. dist_too_long_i) THEN
254 ! 2. compute 3-center integrals (Pν|µ) ("|": truncated Coulomb operator)
255 CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_gocc, &
256 atoms_ao_1=i_atoms, atoms_ao_2=il_atoms)
257 ! 3. tensor operation M_Pνλ(iτ) = sum_µ (Pν|µ) G^occ_λµ(i|τ|,k=0)
258 CALL g_times_3c(t_3c_for_gocc, t_2c_gocc, t_3c_x_gocc, bs_env, &
259 j_atoms, i_atoms, il_atoms)
260 END IF
261 IF (.NOT. dist_too_long_j) THEN
262 ! 4. compute 3-center integrals (Qλ|σ) ("|": truncated Coulomb operator)
263 CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_gvir, &
264 atoms_ao_1=j_atoms, atoms_ao_2=il_atoms)
265 ! 5. tensor operation N_Qλν(iτ) = sum_σ (Qλ|σ) G^vir_νσ(i|τ|,k=0)
266 CALL g_times_3c(t_3c_for_gvir, t_2c_gvir, t_3c_x_gvir, bs_env, &
267 i_atoms, j_atoms, il_atoms)
268 END IF
269 END DO ! IL_atoms
270
271 ! 6. reorder tensors: M_Pνλ -> M_Pλν
272 CALL dbt_copy(t_3c_x_gocc, t_3c_x_gocc_2, move_data=.true., order=[1, 3, 2])
273 CALL dbt_copy(t_3c_x_gvir, t_3c_x_gvir_2, move_data=.true.)
274
275 ! 7. tensor operation χ_PQ(iτ,k=0) = sum_λν M_Pλν(iτ) N_Qλν(iτ),
276 ! Bounds:
277 ! "comb" (combined index)
278 ! -> λ bounds from j_atoms
279 ! -> ν bounds from i_atoms
280 ! P -> sparse in ν (see 3.)
281 ! Q -> sparse in λ (see 5.)
282 bounds_comb(1:2, 1) = [bs_env%i_ao_start_from_atom(j_atoms(1)), &
283 bs_env%i_ao_end_from_atom(j_atoms(2))]
284 bounds_comb(1:2, 2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
285 bs_env%i_ao_end_from_atom(i_atoms(2))]
286
287 CALL get_bounds_from_atoms(bounds_p, i_atoms, [1, bs_env%n_atom], &
288 bs_env%min_RI_idx_from_AO_AO_atom, &
289 bs_env%max_RI_idx_from_AO_AO_atom)
290 CALL get_bounds_from_atoms(bounds_q, [1, bs_env%n_atom], j_atoms, &
291 bs_env%min_RI_idx_from_AO_AO_atom, &
292 bs_env%max_RI_idx_from_AO_AO_atom)
293
294 IF (bounds_q(1) > bounds_q(2) .OR. bounds_p(1) > bounds_p(2)) THEN
295 flop = 0_int_8
296 ELSE
297 CALL dbt_contract(alpha=bs_env%spin_degeneracy, &
298 tensor_1=t_3c_x_gocc_2, tensor_2=t_3c_x_gvir_2, &
299 beta=1.0_dp, tensor_3=bs_env%t_chi, &
300 contract_1=[2, 3], notcontract_1=[1], map_1=[1], &
301 contract_2=[2, 3], notcontract_2=[1], map_2=[2], &
302 bounds_1=bounds_comb, &
303 bounds_2=bounds_p, &
304 bounds_3=bounds_q, &
305 filter_eps=bs_env%eps_filter, move_data=.false., flop=flop)
306 END IF
307 IF (flop == 0_int_8) bs_env%skip_chi(i_intval_idx, j_intval_idx) = .true.
308
309 END DO ! j_atoms
310 END DO ! i_atoms
311 END DO ! ispin
312
313 ! 8. communicate data of χ_PQ(iτ,k=0) in tensor bs_env%t_chi (which local in the
314 ! subgroup) to the global dbcsr matrix mat_chi_Gamma_tau (which stores
315 ! χ_PQ(iτ,k=0) for all time points)
316 CALL local_dbt_to_global_mat(bs_env%t_chi, bs_env%mat_RI_RI_tensor%matrix, &
317 mat_chi_gamma_tau(i_t)%matrix, bs_env%para_env)
318
319 CALL write_matrix(mat_chi_gamma_tau(i_t)%matrix, i_t, bs_env%chi_name, &
320 bs_env%fm_RI_RI, qs_env)
321
322 CALL destroy_tensors_chi(t_2c_gocc, t_2c_gvir, t_3c_for_gocc, t_3c_for_gvir, &
323 t_3c_x_gocc, t_3c_x_gvir, t_3c_x_gocc_2, t_3c_x_gvir_2)
324
325 IF (bs_env%unit_nr > 0) THEN
326 WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F10.1,A)') &
327 'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
328 ', Execution time', m_walltime() - t1, ' s'
329 END IF
330
331 END DO ! i_t
332
333 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
334
335 CALL timestop(handle)
336
337 END SUBROUTINE get_mat_chi_gamma_tau
338
339! **************************************************************************************************
340!> \brief ...
341!> \param fm ...
342!> \param bs_env ...
343!> \param mat_name ...
344!> \param idx ...
345! **************************************************************************************************
346 SUBROUTINE fm_read(fm, bs_env, mat_name, idx)
347 TYPE(cp_fm_type) :: fm
348 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
349 CHARACTER(LEN=*) :: mat_name
350 INTEGER :: idx
351
352 CHARACTER(LEN=*), PARAMETER :: routinen = 'fm_read'
353
354 CHARACTER(LEN=default_path_length) :: f_chi
355 INTEGER :: handle, unit_nr
356
357 CALL timeset(routinen, handle)
358
359 unit_nr = -1
360 IF (bs_env%para_env%is_source()) THEN
361
362 IF (idx < 10) THEN
363 WRITE (f_chi, '(3A,I1,A)') trim(bs_env%prefix), trim(mat_name), "_0", idx, ".matrix"
364 ELSE IF (idx < 100) THEN
365 WRITE (f_chi, '(3A,I2,A)') trim(bs_env%prefix), trim(mat_name), "_", idx, ".matrix"
366 ELSE
367 cpabort('Please implement more than 99 time/frequency points.')
368 END IF
369
370 CALL open_file(file_name=trim(f_chi), file_action="READ", file_form="UNFORMATTED", &
371 file_position="REWIND", file_status="OLD", unit_number=unit_nr)
372
373 END IF
374
375 CALL cp_fm_read_unformatted(fm, unit_nr)
376
377 IF (bs_env%para_env%is_source()) CALL close_file(unit_number=unit_nr)
378
379 CALL timestop(handle)
380
381 END SUBROUTINE fm_read
382
383! **************************************************************************************************
384!> \brief ...
385!> \param t_2c_Gocc ...
386!> \param t_2c_Gvir ...
387!> \param t_3c_for_Gocc ...
388!> \param t_3c_for_Gvir ...
389!> \param t_3c_x_Gocc ...
390!> \param t_3c_x_Gvir ...
391!> \param t_3c_x_Gocc_2 ...
392!> \param t_3c_x_Gvir_2 ...
393!> \param bs_env ...
394! **************************************************************************************************
395 SUBROUTINE create_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
396 t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2, bs_env)
397
398 TYPE(dbt_type) :: t_2c_gocc, t_2c_gvir, t_3c_for_gocc, &
399 t_3c_for_gvir, t_3c_x_gocc, &
400 t_3c_x_gvir, t_3c_x_gocc_2, &
401 t_3c_x_gvir_2
402 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
403
404 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_tensors_chi'
405
406 INTEGER :: handle
407
408 CALL timeset(routinen, handle)
409
410 CALL dbt_create(bs_env%t_G, t_2c_gocc, name="Gocc 2c (AO|AO)")
411 CALL dbt_create(bs_env%t_G, t_2c_gvir, name="Gvir 2c (AO|AO)")
412 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_gocc, name="Gocc 3c (RI AO|AO)")
413 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_gvir, name="Gvir 3c (RI AO|AO)")
414 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_gocc, name="xGocc 3c (RI AO|AO)")
415 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_gvir, name="xGvir 3c (RI AO|AO)")
416 CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_gocc_2, name="x2Gocc 3c (RI AO|AO)")
417 CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_gvir_2, name="x2Gvir 3c (RI AO|AO)")
418
419 CALL timestop(handle)
420
421 END SUBROUTINE create_tensors_chi
422
423! **************************************************************************************************
424!> \brief ...
425!> \param t_2c_Gocc ...
426!> \param t_2c_Gvir ...
427!> \param t_3c_for_Gocc ...
428!> \param t_3c_for_Gvir ...
429!> \param t_3c_x_Gocc ...
430!> \param t_3c_x_Gvir ...
431!> \param t_3c_x_Gocc_2 ...
432!> \param t_3c_x_Gvir_2 ...
433! **************************************************************************************************
434 SUBROUTINE destroy_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
435 t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2)
436 TYPE(dbt_type) :: t_2c_gocc, t_2c_gvir, t_3c_for_gocc, &
437 t_3c_for_gvir, t_3c_x_gocc, &
438 t_3c_x_gvir, t_3c_x_gocc_2, &
439 t_3c_x_gvir_2
440
441 CHARACTER(LEN=*), PARAMETER :: routinen = 'destroy_tensors_chi'
442
443 INTEGER :: handle
444
445 CALL timeset(routinen, handle)
446
447 CALL dbt_destroy(t_2c_gocc)
448 CALL dbt_destroy(t_2c_gvir)
449 CALL dbt_destroy(t_3c_for_gocc)
450 CALL dbt_destroy(t_3c_for_gvir)
451 CALL dbt_destroy(t_3c_x_gocc)
452 CALL dbt_destroy(t_3c_x_gvir)
453 CALL dbt_destroy(t_3c_x_gocc_2)
454 CALL dbt_destroy(t_3c_x_gvir_2)
455
456 CALL timestop(handle)
457
458 END SUBROUTINE destroy_tensors_chi
459
460! **************************************************************************************************
461!> \brief ...
462!> \param matrix ...
463!> \param matrix_index ...
464!> \param matrix_name ...
465!> \param fm ...
466!> \param qs_env ...
467! **************************************************************************************************
468 SUBROUTINE write_matrix(matrix, matrix_index, matrix_name, fm, qs_env)
469 TYPE(dbcsr_type) :: matrix
470 INTEGER :: matrix_index
471 CHARACTER(LEN=*) :: matrix_name
472 TYPE(cp_fm_type), INTENT(IN), POINTER :: fm
473 TYPE(qs_environment_type), POINTER :: qs_env
474
475 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_matrix'
476
477 INTEGER :: handle
478
479 CALL timeset(routinen, handle)
480
481 CALL cp_fm_set_all(fm, 0.0_dp)
482
483 CALL copy_dbcsr_to_fm(matrix, fm)
484
485 CALL fm_write(fm, matrix_index, matrix_name, qs_env)
486
487 CALL timestop(handle)
488
489 END SUBROUTINE write_matrix
490
491! **************************************************************************************************
492!> \brief ...
493!> \param fm ...
494!> \param matrix_index ...
495!> \param matrix_name ...
496!> \param qs_env ...
497! **************************************************************************************************
498 SUBROUTINE fm_write(fm, matrix_index, matrix_name, qs_env)
499 TYPE(cp_fm_type) :: fm
500 INTEGER :: matrix_index
501 CHARACTER(LEN=*) :: matrix_name
502 TYPE(qs_environment_type), POINTER :: qs_env
503
504 CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
505 routinen = 'fm_write'
506
507 CHARACTER(LEN=default_path_length) :: filename
508 INTEGER :: handle, unit_nr
509 TYPE(cp_logger_type), POINTER :: logger
510 TYPE(section_vals_type), POINTER :: input
511
512 CALL timeset(routinen, handle)
513
514 CALL get_qs_env(qs_env, input=input)
515
516 logger => cp_get_default_logger()
517
518 IF (btest(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
519
520 IF (matrix_index < 10) THEN
521 WRITE (filename, '(3A,I1)') "RESTART_", matrix_name, "_0", matrix_index
522 ELSE IF (matrix_index < 100) THEN
523 WRITE (filename, '(3A,I2)') "RESTART_", matrix_name, "_", matrix_index
524 ELSE
525 cpabort('Please implement more than 99 time/frequency points.')
526 END IF
527
528 unit_nr = cp_print_key_unit_nr(logger, input, key, extension=".matrix", &
529 file_form="UNFORMATTED", middle_name=trim(filename), &
530 file_position="REWIND", file_action="WRITE")
531
532 CALL cp_fm_write_unformatted(fm, unit_nr)
533 IF (unit_nr > 0) THEN
534 CALL close_file(unit_nr)
535 END IF
536 END IF
537
538 CALL timestop(handle)
539
540 END SUBROUTINE fm_write
541
542! **************************************************************************************************
543!> \brief ...
544!> \param bs_env ...
545!> \param tau ...
546!> \param fm_G_Gamma ...
547!> \param ispin ...
548!> \param occ ...
549!> \param vir ...
550! **************************************************************************************************
551 SUBROUTINE g_occ_vir(bs_env, tau, fm_G_Gamma, ispin, occ, vir)
552 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
553 REAL(kind=dp) :: tau
554 TYPE(cp_fm_type) :: fm_g_gamma
555 INTEGER :: ispin
556 LOGICAL :: occ, vir
557
558 CHARACTER(LEN=*), PARAMETER :: routinen = 'G_occ_vir'
559
560 INTEGER :: handle, homo, i_row_local, j_col, &
561 j_col_local, n_mo, ncol_local, &
562 nrow_local
563 INTEGER, DIMENSION(:), POINTER :: col_indices
564 REAL(kind=dp) :: tau_e
565
566 CALL timeset(routinen, handle)
567
568 cpassert(occ .NEQV. vir)
569
570 CALL cp_fm_get_info(matrix=bs_env%fm_work_mo(1), &
571 nrow_local=nrow_local, &
572 ncol_local=ncol_local, &
573 col_indices=col_indices)
574
575 n_mo = bs_env%n_ao
576 homo = bs_env%n_occ(ispin)
577
578 CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(ispin), bs_env%fm_work_mo(1))
579
580 DO i_row_local = 1, nrow_local
581 DO j_col_local = 1, ncol_local
582
583 j_col = col_indices(j_col_local)
584
585 tau_e = abs(tau*0.5_dp*(bs_env%eigenval_scf_Gamma(j_col, ispin) - bs_env%e_fermi(ispin)))
586
587 IF (tau_e < bs_env%stabilize_exp) THEN
588 bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = &
589 bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local)*exp(-tau_e)
590 ELSE
591 bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = 0.0_dp
592 END IF
593
594 IF ((occ .AND. j_col > homo) .OR. (vir .AND. j_col <= homo)) THEN
595 bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = 0.0_dp
596 END IF
597
598 END DO
599 END DO
600
601 CALL parallel_gemm(transa="N", transb="T", m=n_mo, n=n_mo, k=n_mo, alpha=1.0_dp, &
602 matrix_a=bs_env%fm_work_mo(1), matrix_b=bs_env%fm_work_mo(1), &
603 beta=0.0_dp, matrix_c=fm_g_gamma)
604
605 CALL timestop(handle)
606
607 END SUBROUTINE g_occ_vir
608
609! **************************************************************************************************
610!> \brief ...
611!> \param qs_env ...
612!> \param bs_env ...
613!> \param t_3c ...
614!> \param atoms_AO_1 ...
615!> \param atoms_AO_2 ...
616!> \param atoms_RI ...
617! **************************************************************************************************
618 SUBROUTINE compute_3c_integrals(qs_env, bs_env, t_3c, atoms_AO_1, atoms_AO_2, atoms_RI)
619 TYPE(qs_environment_type), POINTER :: qs_env
620 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
621 TYPE(dbt_type) :: t_3c
622 INTEGER, DIMENSION(2), OPTIONAL :: atoms_ao_1, atoms_ao_2, atoms_ri
623
624 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_3c_integrals'
625
626 INTEGER :: handle
627 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_array
628
629 CALL timeset(routinen, handle)
630
631 ! free memory (not clear whether memory has been freed previously)
632 CALL dbt_clear(t_3c)
633
634 ALLOCATE (t_3c_array(1, 1))
635 CALL dbt_create(t_3c, t_3c_array(1, 1))
636
637 CALL build_3c_integrals(t_3c_array, &
638 bs_env%eps_filter, &
639 qs_env, &
640 bs_env%nl_3c, &
641 int_eps=bs_env%eps_filter, &
642 basis_i=bs_env%basis_set_RI, &
643 basis_j=bs_env%basis_set_AO, &
644 basis_k=bs_env%basis_set_AO, &
645 potential_parameter=bs_env%ri_metric, &
646 bounds_i=atoms_ri, &
647 bounds_j=atoms_ao_1, &
648 bounds_k=atoms_ao_2, &
649 desymmetrize=.false.)
650
651 CALL dbt_filter(t_3c_array(1, 1), bs_env%eps_filter)
652
653 CALL dbt_copy(t_3c_array(1, 1), t_3c, move_data=.true.)
654
655 CALL dbt_destroy(t_3c_array(1, 1))
656 DEALLOCATE (t_3c_array)
657
658 CALL timestop(handle)
659
660 END SUBROUTINE compute_3c_integrals
661
662! **************************************************************************************************
663!> \brief ...
664!> \param t_3c_for_G ...
665!> \param t_G ...
666!> \param t_M ...
667!> \param bs_env ...
668!> \param atoms_AO_1 ...
669!> \param atoms_AO_2 ...
670!> \param atoms_IL ...
671! **************************************************************************************************
672 SUBROUTINE g_times_3c(t_3c_for_G, t_G, t_M, bs_env, atoms_AO_1, atoms_AO_2, atoms_IL)
673 TYPE(dbt_type) :: t_3c_for_g, t_g, t_m
674 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
675 INTEGER, DIMENSION(2) :: atoms_ao_1, atoms_ao_2, atoms_il
676
677 CHARACTER(LEN=*), PARAMETER :: routinen = 'G_times_3c'
678
679 INTEGER :: handle
680 INTEGER(KIND=int_8) :: flop
681 INTEGER, DIMENSION(2) :: bounds_ao_1, bounds_il
682 INTEGER, DIMENSION(2, 2) :: bounds_comb
683
684 CALL timeset(routinen, handle)
685
686 ! Bounds reduce needed memory and therefore scaling behavior
687 ! Operations are of the form, e.g, M_Pνλ = sum_µ (Pν|µ) G_λµ
688 ! "comb" (combined index)
689 ! -> P sparse in ν and µ
690 ! -> λ bounds from j_atoms (via atoms_AO_1)
691 ! µ bounds from inner loop "IL" indices and sparse in P and ν
692 ! ν bounds from i_atoms (via atoms_AO_2) and sparse in P and µ
693
694 ! µ index
695 CALL get_bounds_from_atoms(bounds_il, [1, bs_env%n_atom], atoms_ao_2, &
696 bs_env%min_AO_idx_from_RI_AO_atom, &
697 bs_env%max_AO_idx_from_RI_AO_atom, &
698 atoms_3=atoms_il, &
699 indices_3_start=bs_env%i_ao_start_from_atom, &
700 indices_3_end=bs_env%i_ao_end_from_atom)
701
702 ! P index
703 CALL get_bounds_from_atoms(bounds_comb(:, 1), atoms_il, atoms_ao_2, &
704 bs_env%min_RI_idx_from_AO_AO_atom, &
705 bs_env%max_RI_idx_from_AO_AO_atom)
706
707 ! ν index
708 CALL get_bounds_from_atoms(bounds_comb(:, 2), [1, bs_env%n_atom], atoms_il, &
709 bs_env%min_AO_idx_from_RI_AO_atom, &
710 bs_env%max_AO_idx_from_RI_AO_atom, &
711 atoms_3=atoms_ao_2, &
712 indices_3_start=bs_env%i_ao_start_from_atom, &
713 indices_3_end=bs_env%i_ao_end_from_atom)
714
715 ! λ index
716 bounds_ao_1(1:2) = [bs_env%i_ao_start_from_atom(atoms_ao_1(1)), &
717 bs_env%i_ao_end_from_atom(atoms_ao_1(2))]
718
719 IF (bounds_il(1) > bounds_il(2) .OR. bounds_comb(1, 2) > bounds_comb(2, 2)) THEN
720 flop = 0_int_8
721 ELSE
722 CALL dbt_contract(alpha=1.0_dp, &
723 tensor_1=t_3c_for_g, &
724 tensor_2=t_g, &
725 beta=1.0_dp, &
726 tensor_3=t_m, &
727 contract_1=[3], notcontract_1=[1, 2], map_1=[1, 2], &
728 contract_2=[2], notcontract_2=[1], map_2=[3], &
729 bounds_1=bounds_il, &
730 bounds_2=bounds_comb, &
731 bounds_3=bounds_ao_1, &
732 flop=flop, &
733 filter_eps=bs_env%eps_filter)
734 END IF
735
736 CALL dbt_clear(t_3c_for_g)
737
738 CALL timestop(handle)
739
740 END SUBROUTINE g_times_3c
741
742! **************************************************************************************************
743!> \brief ...
744!> \param atoms_1 ...
745!> \param atoms_2 ...
746!> \param qs_env ...
747!> \param bs_env ...
748!> \param dist_too_long ...
749! **************************************************************************************************
750 SUBROUTINE check_dist(atoms_1, atoms_2, qs_env, bs_env, dist_too_long)
751 INTEGER, DIMENSION(2) :: atoms_1, atoms_2
752 TYPE(qs_environment_type), POINTER :: qs_env
753 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
754 LOGICAL :: dist_too_long
755
756 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_dist'
757
758 INTEGER :: atom_1, atom_2, handle
759 REAL(dp) :: abs_rab, min_dist_ao_atoms
760 REAL(kind=dp), DIMENSION(3) :: rab
761 TYPE(cell_type), POINTER :: cell
762 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
763
764 CALL timeset(routinen, handle)
765
766 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
767
768 min_dist_ao_atoms = huge(1.0_dp)
769 DO atom_1 = atoms_1(1), atoms_1(2)
770 DO atom_2 = atoms_2(1), atoms_2(2)
771 rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
772
773 abs_rab = sqrt(rab(1)**2 + rab(2)**2 + rab(3)**2)
774
775 min_dist_ao_atoms = min(min_dist_ao_atoms, abs_rab)
776 END DO
777 END DO
778
779 dist_too_long = (min_dist_ao_atoms > bs_env%max_dist_AO_atoms)
780
781 CALL timestop(handle)
782
783 END SUBROUTINE check_dist
784
785! **************************************************************************************************
786!> \brief ...
787!> \param bs_env ...
788!> \param qs_env ...
789!> \param mat_chi_Gamma_tau ...
790!> \param fm_W_MIC_time ...
791! **************************************************************************************************
792 SUBROUTINE get_w_mic(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
793 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
794 TYPE(qs_environment_type), POINTER :: qs_env
795 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
796 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
797
798 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_W_MIC'
799
800 INTEGER :: handle
801
802 CALL timeset(routinen, handle)
803
804 IF (bs_env%all_W_exist) THEN
805 CALL read_w_mic_time(bs_env, mat_chi_gamma_tau, fm_w_mic_time)
806 ELSE
807 CALL compute_w_mic(bs_env, qs_env, mat_chi_gamma_tau, fm_w_mic_time)
808 END IF
809
810 CALL timestop(handle)
811
812 END SUBROUTINE get_w_mic
813
814! **************************************************************************************************
815!> \brief ...
816!> \param bs_env ...
817!> \param qs_env ...
818!> \param fm_V_kp ...
819!> \param ikp_batch ...
820! **************************************************************************************************
821 SUBROUTINE compute_v_k_by_lattice_sum(bs_env, qs_env, fm_V_kp, ikp_batch)
822 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
823 TYPE(qs_environment_type), POINTER :: qs_env
824 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_v_kp
825 INTEGER :: ikp_batch
826
827 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_V_k_by_lattice_sum'
828
829 INTEGER :: handle, ikp, ikp_end, ikp_start, &
830 nkp_chi_eps_w_batch, re_im
831 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
832 TYPE(cell_type), POINTER :: cell
833 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_v_kp
834 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
835 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
836
837 CALL timeset(routinen, handle)
838
839 nkp_chi_eps_w_batch = bs_env%nkp_chi_eps_W_batch
840
841 ikp_start = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + 1
842 ikp_end = min(ikp_batch*bs_env%nkp_chi_eps_W_batch, bs_env%kpoints_chi_eps_W%nkp)
843
844 NULLIFY (mat_v_kp)
845 ALLOCATE (mat_v_kp(ikp_start:ikp_end, 2))
846
847 DO re_im = 1, 2
848 DO ikp = ikp_start, ikp_end
849 NULLIFY (mat_v_kp(ikp, re_im)%matrix)
850 ALLOCATE (mat_v_kp(ikp, re_im)%matrix)
851 CALL dbcsr_create(mat_v_kp(ikp, re_im)%matrix, template=bs_env%mat_RI_RI%matrix)
852 CALL dbcsr_reserve_all_blocks(mat_v_kp(ikp, re_im)%matrix)
853 CALL dbcsr_set(mat_v_kp(ikp, re_im)%matrix, 0.0_dp)
854 END DO ! ikp
855 END DO ! re_im
856
857 CALL get_qs_env(qs_env=qs_env, &
858 particle_set=particle_set, &
859 cell=cell, &
860 qs_kind_set=qs_kind_set, &
861 atomic_kind_set=atomic_kind_set)
862
863 IF (ikp_end <= bs_env%nkp_chi_eps_W_orig) THEN
864
865 ! 1. 2c Coulomb integrals for the first "original" k-point grid
866 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
867
868 ELSE IF (ikp_start > bs_env%nkp_chi_eps_W_orig .AND. &
869 ikp_end <= bs_env%nkp_chi_eps_W_orig_plus_extra) THEN
870
871 ! 2. 2c Coulomb integrals for the second "extrapolation" k-point grid
872 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_extra
873
874 ELSE
875
876 cpabort("Error with k-point parallelization.")
877
878 END IF
879
880 CALL build_2c_coulomb_matrix_kp(mat_v_kp, &
881 bs_env%kpoints_chi_eps_W, &
882 basis_type="RI_AUX", &
883 cell=cell, &
884 particle_set=particle_set, &
885 qs_kind_set=qs_kind_set, &
886 atomic_kind_set=atomic_kind_set, &
887 size_lattice_sum=bs_env%size_lattice_sum_V, &
888 operator_type=operator_coulomb, &
889 ikp_start=ikp_start, &
890 ikp_end=ikp_end)
891
892 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
893
894 ALLOCATE (fm_v_kp(ikp_start:ikp_end, 2))
895 DO re_im = 1, 2
896 DO ikp = ikp_start, ikp_end
897 CALL cp_fm_create(fm_v_kp(ikp, re_im), bs_env%fm_RI_RI%matrix_struct)
898 CALL copy_dbcsr_to_fm(mat_v_kp(ikp, re_im)%matrix, fm_v_kp(ikp, re_im))
899 CALL dbcsr_deallocate_matrix(mat_v_kp(ikp, re_im)%matrix)
900 END DO
901 END DO
902 DEALLOCATE (mat_v_kp)
903
904 CALL timestop(handle)
905
906 END SUBROUTINE compute_v_k_by_lattice_sum
907
908! **************************************************************************************************
909!> \brief ...
910!> \param bs_env ...
911!> \param qs_env ...
912!> \param fm_V_kp ...
913!> \param cfm_V_sqrt_ikp ...
914!> \param cfm_M_inv_V_sqrt_ikp ...
915!> \param ikp ...
916! **************************************************************************************************
917 SUBROUTINE compute_minvvsqrt_vsqrt(bs_env, qs_env, fm_V_kp, cfm_V_sqrt_ikp, &
918 cfm_M_inv_V_sqrt_ikp, ikp)
919 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
920 TYPE(qs_environment_type), POINTER :: qs_env
921 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_v_kp
922 TYPE(cp_cfm_type) :: cfm_v_sqrt_ikp, cfm_m_inv_v_sqrt_ikp
923 INTEGER :: ikp
924
925 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_MinvVsqrt_Vsqrt'
926
927 INTEGER :: handle, info, n_ri
928 TYPE(cp_cfm_type) :: cfm_m_inv_ikp, cfm_work
929 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_m_ikp
930
931 CALL timeset(routinen, handle)
932
933 n_ri = bs_env%n_RI
934
935 ! get here M(k) and write it to fm_M_ikp
936 CALL ri_2c_integral_mat(qs_env, fm_m_ikp, fm_v_kp(ikp, 1), &
937 n_ri, bs_env%ri_metric, do_kpoints=.true., &
938 kpoints=bs_env%kpoints_chi_eps_W, &
939 regularization_ri=bs_env%regularization_RI, ikp_ext=ikp, &
940 do_build_cell_index=(ikp == 1))
941
942 IF (ikp == 1) THEN
943 CALL cp_cfm_create(cfm_v_sqrt_ikp, fm_v_kp(ikp, 1)%matrix_struct)
944 CALL cp_cfm_create(cfm_m_inv_v_sqrt_ikp, fm_v_kp(ikp, 1)%matrix_struct)
945 END IF
946 CALL cp_cfm_create(cfm_m_inv_ikp, fm_v_kp(ikp, 1)%matrix_struct)
947
948 CALL cp_fm_to_cfm(fm_m_ikp(1, 1), fm_m_ikp(1, 2), cfm_m_inv_ikp)
949 CALL cp_fm_to_cfm(fm_v_kp(ikp, 1), fm_v_kp(ikp, 2), cfm_v_sqrt_ikp)
950
951 CALL cp_fm_release(fm_m_ikp)
952
953 CALL cp_cfm_create(cfm_work, fm_v_kp(ikp, 1)%matrix_struct)
954
955 ! M(k) -> M^-1(k)
956 CALL cp_cfm_to_cfm(cfm_m_inv_ikp, cfm_work)
957 CALL cp_cfm_cholesky_decompose(matrix=cfm_m_inv_ikp, n=n_ri, info_out=info)
958 IF (info == 0) THEN
959 ! successful Cholesky decomposition
960 CALL cp_cfm_cholesky_invert(cfm_m_inv_ikp)
961 ! symmetrize the result
962 CALL cp_cfm_uplo_to_full(cfm_m_inv_ikp)
963 ELSE
964 ! Cholesky decomposition not successful: use expensive diagonalization
965 CALL cp_cfm_power(cfm_work, threshold=bs_env%eps_eigval_mat_RI, exponent=-1.0_dp)
966 CALL cp_cfm_to_cfm(cfm_work, cfm_m_inv_ikp)
967 END IF
968
969 ! V(k) -> L(k) with L^H(k)*L(k) = V(k) [L(k) can be just considered to be V^0.5(k)]
970 CALL cp_cfm_to_cfm(cfm_v_sqrt_ikp, cfm_work)
971 CALL cp_cfm_cholesky_decompose(matrix=cfm_v_sqrt_ikp, n=n_ri, info_out=info)
972 IF (info == 0) THEN
973 ! successful Cholesky decomposition
974 CALL clean_lower_part(cfm_v_sqrt_ikp)
975 ELSE
976 ! Cholesky decomposition not successful: use expensive diagonalization
977 CALL cp_cfm_power(cfm_work, threshold=0.0_dp, exponent=0.5_dp)
978 CALL cp_cfm_to_cfm(cfm_work, cfm_v_sqrt_ikp)
979 END IF
980 CALL cp_cfm_release(cfm_work)
981
982 ! get M^-1(k)*V^0.5(k)
983 CALL parallel_gemm("N", "C", n_ri, n_ri, n_ri, z_one, cfm_m_inv_ikp, cfm_v_sqrt_ikp, &
984 z_zero, cfm_m_inv_v_sqrt_ikp)
985
986 CALL cp_cfm_release(cfm_m_inv_ikp)
987
988 CALL timestop(handle)
989
990 END SUBROUTINE compute_minvvsqrt_vsqrt
991
992! **************************************************************************************************
993!> \brief ...
994!> \param bs_env ...
995!> \param mat_chi_Gamma_tau ...
996!> \param fm_W_MIC_time ...
997! **************************************************************************************************
998 SUBROUTINE read_w_mic_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
999 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1000 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
1001 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
1002
1003 CHARACTER(LEN=*), PARAMETER :: routinen = 'read_W_MIC_time'
1004
1005 INTEGER :: handle, i_t
1006 REAL(kind=dp) :: t1
1007
1008 CALL timeset(routinen, handle)
1009
1010 CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
1011 CALL create_fm_w_mic_time(bs_env, fm_w_mic_time)
1012
1013 DO i_t = 1, bs_env%num_time_freq_points
1014
1015 t1 = m_walltime()
1016
1017 CALL fm_read(fm_w_mic_time(i_t), bs_env, bs_env%W_time_name, i_t)
1018
1019 IF (bs_env%unit_nr > 0) THEN
1020 WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
1021 'Read W^MIC(iτ) from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1022 ', Execution time', m_walltime() - t1, ' s'
1023 END IF
1024
1025 END DO
1026
1027 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1028
1029 ! Marek : Reading of the W(w=0) potential for RTP
1030 ! TODO : is the condition bs_env%all_W_exist sufficient for reading?
1031 ! This block builds
1032 ! bs_env%fm_W_MIC_freq_zero specifically for RT-BSE consumption (read by
1033 ! rt_bse_linearized.F initialize_cohsex_selfenergy and by
1034 ! rt_bse_ri_rs.F rt_bse_ri_rs_ensure_W0_grid). RT-BSE-specific compute
1035 ! embedded in GW; left here because moving it would require keeping
1036 ! fm_W_MIC_time alive past compute_W_MIC.
1037 IF (bs_env%rtp_method == rtp_method_bse .OR. &
1038 bs_env%rtp_method == rtp_method_bse_linearized) THEN
1039 CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
1040 t1 = m_walltime()
1041 CALL fm_read(bs_env%fm_W_MIC_freq_zero, bs_env, "W_freq_rtp", 0)
1042 IF (bs_env%unit_nr > 0) THEN
1043 WRITE (bs_env%unit_nr, '(T2,A,I3,A,I3,A,F10.1,A)') &
1044 'Read W^MIC(f=0) from file for freq. point ', 1, ' /', 1, &
1045 ', Execution time', m_walltime() - t1, ' s'
1046 END IF
1047 END IF
1048
1049 CALL timestop(handle)
1050
1051 END SUBROUTINE read_w_mic_time
1052
1053! **************************************************************************************************
1054!> \brief ...
1055!> \param bs_env ...
1056!> \param qs_env ...
1057!> \param mat_chi_Gamma_tau ...
1058!> \param fm_W_MIC_time ...
1059! **************************************************************************************************
1060 SUBROUTINE compute_w_mic(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
1061 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1062 TYPE(qs_environment_type), POINTER :: qs_env
1063 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
1064 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
1065
1066 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_W_MIC'
1067
1068 INTEGER :: handle, i_t, ikp, ikp_batch, &
1069 ikp_in_batch, j_w
1070 REAL(kind=dp) :: t1
1071 TYPE(cp_cfm_type) :: cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp
1072 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_v_kp
1073
1074 CALL timeset(routinen, handle)
1075
1076 CALL create_fm_w_mic_time(bs_env, fm_w_mic_time)
1077
1078 DO ikp_batch = 1, bs_env%num_chi_eps_W_batches
1079
1080 t1 = m_walltime()
1081
1082 ! Compute V_PQ(k) = sum_R e^(ikR) <phi_P, cell 0 | 1/r | phi_Q, cell R>
1083 CALL compute_v_k_by_lattice_sum(bs_env, qs_env, fm_v_kp, ikp_batch)
1084
1085 DO ikp_in_batch = 1, bs_env%nkp_chi_eps_W_batch
1086
1087 ikp = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + ikp_in_batch
1088
1089 IF (ikp > bs_env%nkp_chi_eps_W_orig_plus_extra) cycle
1090
1091 CALL compute_minvvsqrt_vsqrt(bs_env, qs_env, fm_v_kp, &
1092 cfm_v_sqrt_ikp, cfm_m_inv_v_sqrt_ikp, ikp)
1093
1094 CALL bs_env%para_env%sync()
1095 CALL cp_fm_release(fm_v_kp(ikp, 1))
1096 CALL cp_fm_release(fm_v_kp(ikp, 2))
1097
1098 DO j_w = 1, bs_env%num_time_freq_points
1099
1100 ! check if we need this (ikp, ω_j) combination for approximate k-point extrapolation
1101 IF (bs_env%approx_kp_extrapol .AND. j_w > 1 .AND. &
1102 ikp > bs_env%nkp_chi_eps_W_orig) cycle
1103
1104 CALL compute_fm_w_mic_freq_j(bs_env, qs_env, bs_env%fm_W_MIC_freq, j_w, ikp, &
1105 mat_chi_gamma_tau, cfm_m_inv_v_sqrt_ikp, &
1106 cfm_v_sqrt_ikp)
1107
1108 ! Fourier trafo from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
1109 CALL fourier_transform_w_to_t(bs_env, fm_w_mic_time, bs_env%fm_W_MIC_freq, j_w)
1110
1111 END DO ! ω_j
1112
1113 END DO ! ikp_in_batch
1114
1115 DEALLOCATE (fm_v_kp)
1116
1117 IF (bs_env%unit_nr > 0) THEN
1118 WRITE (bs_env%unit_nr, '(T2,A,I12,A,I3,A,F10.1,A)') &
1119 'Computed W(iτ,k) for k-point batch', &
1120 ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
1121 ', Execution time', m_walltime() - t1, ' s'
1122 END IF
1123
1124 END DO ! ikp_batch
1125
1126 IF (bs_env%approx_kp_extrapol) THEN
1127 CALL apply_extrapol_factor(bs_env, fm_w_mic_time)
1128 END IF
1129
1130 ! M^-1(k=0)*W^MIC(iτ)*M^-1(k=0)
1131 CALL multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_w_mic_time)
1132
1133 DO i_t = 1, bs_env%num_time_freq_points
1134 CALL fm_write(fm_w_mic_time(i_t), i_t, bs_env%W_time_name, qs_env)
1135 END DO
1136
1137 CALL cp_cfm_release(cfm_m_inv_v_sqrt_ikp)
1138 CALL cp_cfm_release(cfm_v_sqrt_ikp)
1139 CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
1140
1141 ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0
1142 ! Same RT-BSE coupling as read_W_MIC_time.
1143 IF (bs_env%rtp_method == rtp_method_bse .OR. &
1144 bs_env%rtp_method == rtp_method_bse_linearized) THEN
1145 t1 = m_walltime()
1146 CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
1147 ! Set to zero
1148 CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
1149 ! Sum over all times
1150 DO i_t = 1, bs_env%num_time_freq_points
1151 ! Add the relevant structure with correct weight
1152 CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
1153 bs_env%imag_time_weights_freq_zero(i_t), fm_w_mic_time(i_t))
1154 END DO
1155 ! Done, save to file
1156 CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
1157 ! Report calculation
1158 IF (bs_env%unit_nr > 0) THEN
1159 WRITE (bs_env%unit_nr, '(T2,A,I11,A,I3,A,F10.1,A)') &
1160 'Computed W(f=0,k) for k-point batch', &
1161 1, ' /', 1, &
1162 ', Execution time', m_walltime() - t1, ' s'
1163 END IF
1164 END IF
1165
1166 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1167
1168 CALL timestop(handle)
1169
1170 END SUBROUTINE compute_w_mic
1171
1172! **************************************************************************************************
1173!> \brief ...
1174!> \param bs_env ...
1175!> \param qs_env ...
1176!> \param fm_W_MIC_freq_j ...
1177!> \param j_w ...
1178!> \param ikp ...
1179!> \param mat_chi_Gamma_tau ...
1180!> \param cfm_M_inv_V_sqrt_ikp ...
1181!> \param cfm_V_sqrt_ikp ...
1182! **************************************************************************************************
1183 SUBROUTINE compute_fm_w_mic_freq_j(bs_env, qs_env, fm_W_MIC_freq_j, j_w, ikp, mat_chi_Gamma_tau, &
1184 cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
1185 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1186 TYPE(qs_environment_type), POINTER :: qs_env
1187 TYPE(cp_fm_type) :: fm_w_mic_freq_j
1188 INTEGER :: j_w, ikp
1189 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
1190 TYPE(cp_cfm_type) :: cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp
1191
1192 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_fm_W_MIC_freq_j'
1193
1194 INTEGER :: handle
1195 TYPE(cp_cfm_type) :: cfm_chi_ikp_freq_j, cfm_w_ikp_freq_j
1196
1197 CALL timeset(routinen, handle)
1198
1199 ! 1. Fourier transformation of χ_PQ(iτ,k=0) to χ_PQ(iω_j,k=0)
1200 CALL compute_fm_chi_gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_gamma_tau)
1201
1202 CALL cp_fm_set_all(fm_w_mic_freq_j, 0.0_dp)
1203
1204 ! 2. Get χ_PQ(iω_j,k_i) from χ_PQ(iω_j,k=0) using the minimum image convention
1205 CALL cfm_ikp_from_fm_gamma(cfm_chi_ikp_freq_j, bs_env%fm_chi_Gamma_freq, &
1206 ikp, qs_env, bs_env%kpoints_chi_eps_W, "RI_AUX")
1207
1208 ! 3. Remove all negative eigenvalues from χ_PQ(iω_j,k_i)
1209 CALL cp_cfm_power(cfm_chi_ikp_freq_j, threshold=0.0_dp, exponent=1.0_dp)
1210
1211 ! 4. ε(iω_j,k_i) = Id - V^0.5(k_i)*M^-1(k_i)*χ(iω_j,k_i)*M^-1(k_i)*V^0.5(k_i)
1212 ! W(iω_j,k_i) = V^0.5(k_i)*(ε^-1(iω_j,k_i)-Id)*V^0.5(k_i)
1213 CALL compute_cfm_w_ikp_freq_j(bs_env, cfm_chi_ikp_freq_j, cfm_v_sqrt_ikp, &
1214 cfm_m_inv_v_sqrt_ikp, cfm_w_ikp_freq_j)
1215
1216 ! 5. k-point integration W_PQ(iω_j, k_i) to W_PQ^MIC(iω_j)
1217 SELECT CASE (bs_env%approx_kp_extrapol)
1218 CASE (.false.)
1219 ! default: standard k-point extrapolation
1220 CALL mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, cfm_w_ikp_freq_j, ikp, &
1221 bs_env%kpoints_chi_eps_W, "RI_AUX")
1222 CASE (.true.)
1223 ! for approximate kpoint extrapolation: get W_PQ^MIC(iω_1) with and without k-point
1224 ! extrapolation to compute the extrapolation factor f_PQ for every PQ-matrix element,
1225 ! f_PQ = (W_PQ^MIC(iω_1) with extrapolation) / (W_PQ^MIC(iω_1) without extrapolation)
1226
1227 ! for ω_1, we compute the k-point extrapolated result using all k-points
1228 IF (j_w == 1) THEN
1229
1230 ! k-point extrapolated
1231 CALL mic_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_extra, &
1232 cfm_w_ikp_freq_j, ikp, bs_env%kpoints_chi_eps_W, &
1233 "RI_AUX")
1234 ! non-kpoint extrapolated
1235 IF (ikp <= bs_env%nkp_chi_eps_W_orig) THEN
1236 CALL mic_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_no_extra, &
1237 cfm_w_ikp_freq_j, ikp, bs_env%kpoints_chi_eps_W, &
1238 "RI_AUX", wkp_ext=bs_env%wkp_orig)
1239 END IF
1240
1241 END IF
1242
1243 ! for all ω_j, we need to compute W^MIC without k-point extrpolation
1244 IF (ikp <= bs_env%nkp_chi_eps_W_orig) THEN
1245 CALL mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, cfm_w_ikp_freq_j, &
1246 ikp, bs_env%kpoints_chi_eps_W, "RI_AUX", &
1247 wkp_ext=bs_env%wkp_orig)
1248 END IF
1249 END SELECT
1250
1251 CALL cp_cfm_release(cfm_w_ikp_freq_j)
1252
1253 CALL timestop(handle)
1254
1255 END SUBROUTINE compute_fm_w_mic_freq_j
1256
1257! **************************************************************************************************
1258!> \brief ...
1259!> \param cfm_mat ...
1260! **************************************************************************************************
1261 SUBROUTINE clean_lower_part(cfm_mat)
1262 TYPE(cp_cfm_type) :: cfm_mat
1263
1264 CHARACTER(LEN=*), PARAMETER :: routinen = 'clean_lower_part'
1265
1266 INTEGER :: handle, i_row, j_col, j_global, &
1267 ncol_local, nrow_local
1268 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1269
1270 CALL timeset(routinen, handle)
1271
1272 CALL cp_cfm_get_info(matrix=cfm_mat, &
1273 nrow_local=nrow_local, ncol_local=ncol_local, &
1274 row_indices=row_indices, col_indices=col_indices)
1275
1276 DO j_col = 1, ncol_local
1277 j_global = col_indices(j_col)
1278 DO i_row = 1, nrow_local
1279 IF (j_global < row_indices(i_row)) cfm_mat%local_data(i_row, j_col) = z_zero
1280 END DO
1281 END DO
1282
1283 CALL timestop(handle)
1284
1285 END SUBROUTINE clean_lower_part
1286
1287! **************************************************************************************************
1288!> \brief ...
1289!> \param bs_env ...
1290!> \param fm_W_MIC_time ...
1291! **************************************************************************************************
1292 SUBROUTINE apply_extrapol_factor(bs_env, fm_W_MIC_time)
1293 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1294 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
1295
1296 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_extrapol_factor'
1297
1298 INTEGER :: handle, i, i_t, j, ncol_local, nrow_local
1299 REAL(kind=dp) :: extrapol_factor, w_extra_1, w_no_extra_1
1300
1301 CALL timeset(routinen, handle)
1302
1303 CALL cp_fm_get_info(matrix=fm_w_mic_time(1), nrow_local=nrow_local, ncol_local=ncol_local)
1304
1305 DO i_t = 1, bs_env%num_time_freq_points
1306 DO j = 1, ncol_local
1307 DO i = 1, nrow_local
1308
1309 w_extra_1 = bs_env%fm_W_MIC_freq_1_extra%local_data(i, j)
1310 w_no_extra_1 = bs_env%fm_W_MIC_freq_1_no_extra%local_data(i, j)
1311
1312 IF (abs(w_no_extra_1) > 1.0e-13) THEN
1313 extrapol_factor = abs(w_extra_1/w_no_extra_1)
1314 ELSE
1315 extrapol_factor = 1.0_dp
1316 END IF
1317
1318 ! reset extrapolation factor if it is very large
1319 IF (extrapol_factor > 10.0_dp) extrapol_factor = 1.0_dp
1320
1321 fm_w_mic_time(i_t)%local_data(i, j) = fm_w_mic_time(i_t)%local_data(i, j) &
1322 *extrapol_factor
1323 END DO
1324 END DO
1325 END DO
1326
1327 CALL timestop(handle)
1328
1329 END SUBROUTINE apply_extrapol_factor
1330
1331! **************************************************************************************************
1332!> \brief ...
1333!> \param bs_env ...
1334!> \param fm_chi_Gamma_freq ...
1335!> \param j_w ...
1336!> \param mat_chi_Gamma_tau ...
1337! **************************************************************************************************
1338 SUBROUTINE compute_fm_chi_gamma_freq(bs_env, fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
1339 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1340 TYPE(cp_fm_type) :: fm_chi_gamma_freq
1341 INTEGER :: j_w
1342 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
1343
1344 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_fm_chi_Gamma_freq'
1345
1346 INTEGER :: handle, i_t
1347 REAL(kind=dp) :: freq_j, time_i, weight_ij
1348
1349 CALL timeset(routinen, handle)
1350
1351 CALL dbcsr_set(bs_env%mat_RI_RI%matrix, 0.0_dp)
1352
1353 freq_j = bs_env%imag_freq_points(j_w)
1354
1355 DO i_t = 1, bs_env%num_time_freq_points
1356
1357 time_i = bs_env%imag_time_points(i_t)
1358 weight_ij = bs_env%weights_cos_t_to_w(j_w, i_t)
1359
1360 ! actual Fourier transform
1361 CALL dbcsr_add(bs_env%mat_RI_RI%matrix, mat_chi_gamma_tau(i_t)%matrix, &
1362 1.0_dp, cos(time_i*freq_j)*weight_ij)
1363
1364 END DO
1365
1366 CALL copy_dbcsr_to_fm(bs_env%mat_RI_RI%matrix, fm_chi_gamma_freq)
1367
1368 CALL timestop(handle)
1369
1370 END SUBROUTINE compute_fm_chi_gamma_freq
1371
1372! **************************************************************************************************
1373!> \brief ...
1374!> \param mat_ikp_re ...
1375!> \param mat_ikp_im ...
1376!> \param mat_Gamma ...
1377!> \param kpoints ...
1378!> \param ikp ...
1379!> \param qs_env ...
1380! **************************************************************************************************
1381 SUBROUTINE mat_ikp_from_mat_gamma(mat_ikp_re, mat_ikp_im, mat_Gamma, kpoints, ikp, qs_env)
1382 TYPE(dbcsr_type) :: mat_ikp_re, mat_ikp_im, mat_gamma
1383 TYPE(kpoint_type), POINTER :: kpoints
1384 INTEGER :: ikp
1385 TYPE(qs_environment_type), POINTER :: qs_env
1386
1387 CHARACTER(LEN=*), PARAMETER :: routinen = 'mat_ikp_from_mat_Gamma'
1388
1389 INTEGER :: col, handle, i_cell, j_cell, num_cells, &
1390 row
1391 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1392 LOGICAL :: f, i_cell_is_the_minimum_image_cell
1393 REAL(kind=dp) :: abs_rab_cell_i, abs_rab_cell_j, arg
1394 REAL(kind=dp), DIMENSION(3) :: cell_vector, cell_vector_j, rab_cell_i, &
1395 rab_cell_j
1396 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1397 REAL(kind=dp), DIMENSION(:, :), POINTER :: block_im, block_re, data_block
1398 TYPE(cell_type), POINTER :: cell
1399 TYPE(dbcsr_iterator_type) :: iter
1400 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1401
1402 CALL timeset(routinen, handle)
1403
1404 ! get the same blocks in mat_ikp_re and mat_ikp_im as in mat_Gamma
1405 CALL dbcsr_copy(mat_ikp_re, mat_gamma)
1406 CALL dbcsr_copy(mat_ikp_im, mat_gamma)
1407 CALL dbcsr_set(mat_ikp_re, 0.0_dp)
1408 CALL dbcsr_set(mat_ikp_im, 0.0_dp)
1409
1410 NULLIFY (cell, particle_set)
1411 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1412 CALL get_cell(cell=cell, h=hmat)
1413
1414 index_to_cell => kpoints%index_to_cell
1415
1416 num_cells = SIZE(index_to_cell, 2)
1417
1418 DO i_cell = 1, num_cells
1419
1420 CALL dbcsr_iterator_start(iter, mat_gamma)
1421 DO WHILE (dbcsr_iterator_blocks_left(iter))
1422 CALL dbcsr_iterator_next_block(iter, row, col, data_block)
1423
1424 cell_vector(1:3) = matmul(hmat, real(index_to_cell(1:3, i_cell), dp))
1425
1426 rab_cell_i(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
1427 (pbc(particle_set(col)%r(1:3), cell) + cell_vector(1:3))
1428 abs_rab_cell_i = sqrt(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
1429
1430 ! minimum image convention
1431 i_cell_is_the_minimum_image_cell = .true.
1432 DO j_cell = 1, num_cells
1433 cell_vector_j(1:3) = matmul(hmat, real(index_to_cell(1:3, j_cell), dp))
1434 rab_cell_j(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
1435 (pbc(particle_set(col)%r(1:3), cell) + cell_vector_j(1:3))
1436 abs_rab_cell_j = sqrt(rab_cell_j(1)**2 + rab_cell_j(2)**2 + rab_cell_j(3)**2)
1437
1438 IF (abs_rab_cell_i > abs_rab_cell_j + 1.0e-6_dp) THEN
1439 i_cell_is_the_minimum_image_cell = .false.
1440 END IF
1441 END DO
1442
1443 IF (i_cell_is_the_minimum_image_cell) THEN
1444 NULLIFY (block_re, block_im)
1445 CALL dbcsr_get_block_p(matrix=mat_ikp_re, row=row, col=col, block=block_re, found=f)
1446 CALL dbcsr_get_block_p(matrix=mat_ikp_im, row=row, col=col, block=block_im, found=f)
1447 cpassert(all(abs(block_re) < 1.0e-10_dp))
1448 cpassert(all(abs(block_im) < 1.0e-10_dp))
1449
1450 arg = real(index_to_cell(1, i_cell), dp)*kpoints%xkp(1, ikp) + &
1451 REAL(index_to_cell(2, i_cell), dp)*kpoints%xkp(2, ikp) + &
1452 REAL(index_to_cell(3, i_cell), dp)*kpoints%xkp(3, ikp)
1453
1454 block_re(:, :) = cos(twopi*arg)*data_block(:, :)
1455 block_im(:, :) = sin(twopi*arg)*data_block(:, :)
1456 END IF
1457
1458 END DO
1459 CALL dbcsr_iterator_stop(iter)
1460
1461 END DO
1462
1463 CALL timestop(handle)
1464
1465 END SUBROUTINE mat_ikp_from_mat_gamma
1466
1467! **************************************************************************************************
1468!> \brief ...
1469!> \param bs_env ...
1470!> \param cfm_chi_ikp_freq_j ...
1471!> \param cfm_V_sqrt_ikp ...
1472!> \param cfm_M_inv_V_sqrt_ikp ...
1473!> \param cfm_W_ikp_freq_j ...
1474! **************************************************************************************************
1475 SUBROUTINE compute_cfm_w_ikp_freq_j(bs_env, cfm_chi_ikp_freq_j, cfm_V_sqrt_ikp, &
1476 cfm_M_inv_V_sqrt_ikp, cfm_W_ikp_freq_j)
1477
1478 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1479 TYPE(cp_cfm_type) :: cfm_chi_ikp_freq_j, cfm_v_sqrt_ikp, &
1480 cfm_m_inv_v_sqrt_ikp, cfm_w_ikp_freq_j
1481
1482 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_cfm_W_ikp_freq_j'
1483
1484 INTEGER :: handle, info, n_ri
1485 TYPE(cp_cfm_type) :: cfm_eps_ikp_freq_j, cfm_work
1486
1487 CALL timeset(routinen, handle)
1488
1489 CALL cp_cfm_create(cfm_work, cfm_chi_ikp_freq_j%matrix_struct)
1490 n_ri = bs_env%n_RI
1491
1492 ! 1. ε(iω_j,k) = Id - V^0.5(k)*M^-1(k)*χ(iω_j,k)*M^-1(k)*V^0.5(k)
1493
1494 ! 1. a) work = χ(iω_j,k)*M^-1(k)*V^0.5(k)
1495 CALL parallel_gemm('N', 'N', n_ri, n_ri, n_ri, z_one, &
1496 cfm_chi_ikp_freq_j, cfm_m_inv_v_sqrt_ikp, z_zero, cfm_work)
1497 CALL cp_cfm_release(cfm_chi_ikp_freq_j)
1498
1499 ! 1. b) eps_work = V^0.5(k)*M^-1(k)*work
1500 CALL cp_cfm_create(cfm_eps_ikp_freq_j, cfm_work%matrix_struct)
1501 CALL parallel_gemm('C', 'N', n_ri, n_ri, n_ri, z_one, &
1502 cfm_m_inv_v_sqrt_ikp, cfm_work, z_zero, cfm_eps_ikp_freq_j)
1503
1504 ! 1. c) ε(iω_j,k) = eps_work - Id
1505 CALL cfm_add_on_diag(cfm_eps_ikp_freq_j, z_one)
1506
1507 ! 2. W(iω_j,k) = V^0.5(k)*(ε^-1(iω_j,k)-Id)*V^0.5(k)
1508
1509 ! 2. a) Cholesky decomposition of ε(iω_j,k) as preparation for inversion
1510 CALL cp_cfm_cholesky_decompose(matrix=cfm_eps_ikp_freq_j, n=n_ri, info_out=info)
1511 cpassert(info == 0)
1512
1513 ! 2. b) Inversion of ε(iω_j,k) using its Cholesky decomposition
1514 CALL cp_cfm_cholesky_invert(cfm_eps_ikp_freq_j)
1515 CALL cp_cfm_uplo_to_full(cfm_eps_ikp_freq_j)
1516
1517 ! 2. c) ε^-1(iω_j,k)-Id
1518 CALL cfm_add_on_diag(cfm_eps_ikp_freq_j, -z_one)
1519
1520 ! 2. d) work = (ε^-1(iω_j,k)-Id)*V^0.5(k)
1521 CALL parallel_gemm('N', 'N', n_ri, n_ri, n_ri, z_one, cfm_eps_ikp_freq_j, cfm_v_sqrt_ikp, &
1522 z_zero, cfm_work)
1523
1524 ! 2. e) W(iw,k) = V^0.5(k)*work
1525 CALL cp_cfm_create(cfm_w_ikp_freq_j, cfm_work%matrix_struct)
1526 CALL parallel_gemm('C', 'N', n_ri, n_ri, n_ri, z_one, cfm_v_sqrt_ikp, cfm_work, &
1527 z_zero, cfm_w_ikp_freq_j)
1528
1529 CALL cp_cfm_release(cfm_work)
1530 CALL cp_cfm_release(cfm_eps_ikp_freq_j)
1531
1532 CALL timestop(handle)
1533
1534 END SUBROUTINE compute_cfm_w_ikp_freq_j
1535
1536! **************************************************************************************************
1537!> \brief ...
1538!> \param cfm ...
1539!> \param alpha ...
1540! **************************************************************************************************
1541 SUBROUTINE cfm_add_on_diag(cfm, alpha)
1542
1543 TYPE(cp_cfm_type) :: cfm
1544 COMPLEX(KIND=dp) :: alpha
1545
1546 CHARACTER(LEN=*), PARAMETER :: routinen = 'cfm_add_on_diag'
1547
1548 INTEGER :: handle, i_row, j_col, j_global, &
1549 ncol_local, nrow_local
1550 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1551
1552 CALL timeset(routinen, handle)
1553
1554 CALL cp_cfm_get_info(matrix=cfm, &
1555 nrow_local=nrow_local, &
1556 ncol_local=ncol_local, &
1557 row_indices=row_indices, &
1558 col_indices=col_indices)
1559
1560 ! add 1 on the diagonal
1561 DO j_col = 1, ncol_local
1562 j_global = col_indices(j_col)
1563 DO i_row = 1, nrow_local
1564 IF (j_global == row_indices(i_row)) THEN
1565 cfm%local_data(i_row, j_col) = cfm%local_data(i_row, j_col) + alpha
1566 END IF
1567 END DO
1568 END DO
1569
1570 CALL timestop(handle)
1571
1572 END SUBROUTINE cfm_add_on_diag
1573
1574! **************************************************************************************************
1575!> \brief ...
1576!> \param bs_env ...
1577!> \param fm_W_MIC_time ...
1578! **************************************************************************************************
1579 SUBROUTINE create_fm_w_mic_time(bs_env, fm_W_MIC_time)
1580 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1581 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
1582
1583 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_fm_W_MIC_time'
1584
1585 INTEGER :: handle, i_t
1586
1587 CALL timeset(routinen, handle)
1588
1589 ALLOCATE (fm_w_mic_time(bs_env%num_time_freq_points))
1590 DO i_t = 1, bs_env%num_time_freq_points
1591 CALL cp_fm_create(fm_w_mic_time(i_t), bs_env%fm_RI_RI%matrix_struct, set_zero=.true.)
1592 END DO
1593
1594 CALL timestop(handle)
1595
1596 END SUBROUTINE create_fm_w_mic_time
1597
1598! **************************************************************************************************
1599!> \brief ...
1600!> \param bs_env ...
1601!> \param fm_W_MIC_time ...
1602!> \param fm_W_MIC_freq_j ...
1603!> \param j_w ...
1604! **************************************************************************************************
1605 SUBROUTINE fourier_transform_w_to_t(bs_env, fm_W_MIC_time, fm_W_MIC_freq_j, j_w)
1606 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1607 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
1608 TYPE(cp_fm_type) :: fm_w_mic_freq_j
1609 INTEGER :: j_w
1610
1611 CHARACTER(LEN=*), PARAMETER :: routinen = 'Fourier_transform_w_to_t'
1612
1613 INTEGER :: handle, i_t
1614 REAL(kind=dp) :: freq_j, time_i, weight_ij
1615
1616 CALL timeset(routinen, handle)
1617
1618 freq_j = bs_env%imag_freq_points(j_w)
1619
1620 DO i_t = 1, bs_env%num_time_freq_points
1621
1622 time_i = bs_env%imag_time_points(i_t)
1623 weight_ij = bs_env%weights_cos_w_to_t(i_t, j_w)
1624
1625 ! actual Fourier transform
1626 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_w_mic_time(i_t), &
1627 beta=weight_ij*cos(time_i*freq_j), matrix_b=fm_w_mic_freq_j)
1628
1629 END DO
1630
1631 CALL timestop(handle)
1632
1633 END SUBROUTINE fourier_transform_w_to_t
1634
1635! **************************************************************************************************
1636!> \brief ...
1637!> \param bs_env ...
1638!> \param qs_env ...
1639!> \param fm_W_MIC_time ...
1640! **************************************************************************************************
1641 SUBROUTINE multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_W_MIC_time)
1642 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1643 TYPE(qs_environment_type), POINTER :: qs_env
1644 TYPE(cp_fm_type), DIMENSION(:) :: fm_w_mic_time
1645
1646 CHARACTER(LEN=*), PARAMETER :: routinen = 'multiply_fm_W_MIC_time_with_Minv_Gamma'
1647
1648 INTEGER :: handle, i_t, n_ri, ndep
1649 TYPE(cp_fm_type) :: fm_work
1650 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_minv_gamma
1651
1652 CALL timeset(routinen, handle)
1653
1654 n_ri = bs_env%n_RI
1655
1656 CALL cp_fm_create(fm_work, fm_w_mic_time(1)%matrix_struct)
1657
1658 ! compute Gamma-only RI-metric matrix M(k=0); no regularization
1659 CALL ri_2c_integral_mat(qs_env, fm_minv_gamma, fm_w_mic_time(1), n_ri, &
1660 bs_env%ri_metric, do_kpoints=.false.)
1661
1662 CALL cp_fm_power(fm_minv_gamma(1, 1), fm_work, -1.0_dp, 0.0_dp, ndep)
1663
1664 ! M^-1(k=0)*W^MIC(iτ)*M^-1(k=0)
1665 DO i_t = 1, SIZE(fm_w_mic_time)
1666
1667 CALL parallel_gemm('N', 'N', n_ri, n_ri, n_ri, 1.0_dp, fm_minv_gamma(1, 1), &
1668 fm_w_mic_time(i_t), 0.0_dp, fm_work)
1669
1670 CALL parallel_gemm('N', 'N', n_ri, n_ri, n_ri, 1.0_dp, fm_work, &
1671 fm_minv_gamma(1, 1), 0.0_dp, fm_w_mic_time(i_t))
1672
1673 END DO
1674
1675 CALL cp_fm_release(fm_work)
1676 CALL cp_fm_release(fm_minv_gamma)
1677
1678 CALL timestop(handle)
1679
1680 END SUBROUTINE multiply_fm_w_mic_time_with_minv_gamma
1681
1682! **************************************************************************************************
1683!> \brief ...
1684!> \param bs_env ...
1685!> \param qs_env ...
1686!> \param fm_Sigma_x_Gamma ...
1687! **************************************************************************************************
1688 SUBROUTINE get_sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
1689 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1690 TYPE(qs_environment_type), POINTER :: qs_env
1691 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_sigma_x_gamma
1692
1693 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_Sigma_x'
1694
1695 INTEGER :: handle, ispin
1696
1697 CALL timeset(routinen, handle)
1698
1699 ALLOCATE (fm_sigma_x_gamma(bs_env%n_spin))
1700 DO ispin = 1, bs_env%n_spin
1701 CALL cp_fm_create(fm_sigma_x_gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1702 END DO
1703
1704 IF (bs_env%Sigma_x_exists) THEN
1705 DO ispin = 1, bs_env%n_spin
1706 CALL fm_read(fm_sigma_x_gamma(ispin), bs_env, bs_env%Sigma_x_name, ispin)
1707 END DO
1708 ELSE
1709 CALL compute_sigma_x(bs_env, qs_env, fm_sigma_x_gamma)
1710 END IF
1711
1712 CALL timestop(handle)
1713
1714 END SUBROUTINE get_sigma_x
1715
1716! **************************************************************************************************
1717!> \brief ...
1718!> \param bs_env ...
1719!> \param qs_env ...
1720!> \param fm_Sigma_x_Gamma ...
1721! **************************************************************************************************
1722 SUBROUTINE compute_sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
1723 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1724 TYPE(qs_environment_type), POINTER :: qs_env
1725 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_sigma_x_gamma
1726
1727 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_Sigma_x'
1728
1729 INTEGER :: handle, i_intval_idx, ispin, j_intval_idx
1730 INTEGER, DIMENSION(2) :: i_atoms, j_atoms
1731 REAL(kind=dp) :: t1
1732 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_vtr_gamma
1733 TYPE(dbcsr_type) :: mat_sigma_x_gamma
1734 TYPE(dbt_type) :: t_2c_d, t_2c_sigma_x, t_2c_v, t_3c_x_v
1735
1736 CALL timeset(routinen, handle)
1737
1738 t1 = m_walltime()
1739
1740 CALL dbt_create(bs_env%t_G, t_2c_d)
1741 CALL dbt_create(bs_env%t_W, t_2c_v)
1742 CALL dbt_create(bs_env%t_G, t_2c_sigma_x)
1743 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_v)
1744 CALL dbcsr_create(mat_sigma_x_gamma, template=bs_env%mat_ao_ao%matrix)
1745
1746 ! 1. Compute truncated Coulomb operator matrix V^tr(k=0) (cutoff rad: cellsize/2)
1747 CALL ri_2c_integral_mat(qs_env, fm_vtr_gamma, bs_env%fm_RI_RI, bs_env%n_RI, &
1748 bs_env%trunc_coulomb, do_kpoints=.false.)
1749
1750 ! 2. Compute M^-1(k=0) and get M^-1(k=0)*V^tr(k=0)*M^-1(k=0)
1751 CALL multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_vtr_gamma(:, 1))
1752
1753 DO ispin = 1, bs_env%n_spin
1754
1755 ! 3. Compute density matrix D_µν
1756 CALL g_occ_vir(bs_env, 0.0_dp, bs_env%fm_work_mo(2), ispin, occ=.true., vir=.false.)
1757
1758 CALL fm_to_local_tensor(bs_env%fm_work_mo(2), bs_env%mat_ao_ao%matrix, &
1759 bs_env%mat_ao_ao_tensor%matrix, t_2c_d, bs_env, &
1760 bs_env%atoms_i_t_group)
1761
1762 CALL fm_to_local_tensor(fm_vtr_gamma(1, 1), bs_env%mat_RI_RI%matrix, &
1763 bs_env%mat_RI_RI_tensor%matrix, t_2c_v, bs_env, &
1764 bs_env%atoms_j_t_group)
1765
1766 ! every group has its own range of i_atoms and j_atoms; only deal with a
1767 ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
1768 DO i_intval_idx = 1, bs_env%n_intervals_i
1769 DO j_intval_idx = 1, bs_env%n_intervals_j
1770 i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
1771 j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
1772
1773 ! 4. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
1774 ! 5. M_Qνσ(iτ) = sum_P (νσ|P) (M^-1(k=0)*V^tr(k=0)*M^-1(k=0))_QP(iτ)
1775 CALL compute_3c_and_contract_w(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_v, t_2c_v)
1776
1777 ! 6. tensor operations with D and computation of Σ^x
1778 ! Σ^x_λσ(k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) D_νµ
1779 CALL contract_to_sigma(t_2c_d, t_3c_x_v, t_2c_sigma_x, i_atoms, j_atoms, &
1780 qs_env, bs_env, occ=.true., vir=.false.)
1781
1782 END DO ! j_atoms
1783 END DO ! i_atoms
1784
1785 CALL local_dbt_to_global_mat(t_2c_sigma_x, bs_env%mat_ao_ao_tensor%matrix, &
1786 mat_sigma_x_gamma, bs_env%para_env)
1787
1788 CALL write_matrix(mat_sigma_x_gamma, ispin, bs_env%Sigma_x_name, &
1789 bs_env%fm_work_mo(1), qs_env)
1790
1791 CALL copy_dbcsr_to_fm(mat_sigma_x_gamma, fm_sigma_x_gamma(ispin))
1792
1793 END DO ! ispin
1794
1795 IF (bs_env%unit_nr > 0) THEN
1796 WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
1797 'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
1798 WRITE (bs_env%unit_nr, '(A)') ' '
1799 END IF
1800
1801 CALL dbcsr_release(mat_sigma_x_gamma)
1802 CALL dbt_destroy(t_2c_d)
1803 CALL dbt_destroy(t_2c_v)
1804 CALL dbt_destroy(t_2c_sigma_x)
1805 CALL dbt_destroy(t_3c_x_v)
1806 CALL cp_fm_release(fm_vtr_gamma)
1807
1808 CALL timestop(handle)
1809
1810 END SUBROUTINE compute_sigma_x
1811
1812! **************************************************************************************************
1813!> \brief ...
1814!> \param bs_env ...
1815!> \param qs_env ...
1816!> \param fm_W_MIC_time ...
1817!> \param fm_Sigma_c_Gamma_time ...
1818! **************************************************************************************************
1819 SUBROUTINE get_sigma_c(bs_env, qs_env, fm_W_MIC_time, fm_Sigma_c_Gamma_time)
1820 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1821 TYPE(qs_environment_type), POINTER :: qs_env
1822 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
1823 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_sigma_c_gamma_time
1824
1825 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_Sigma_c'
1826
1827 INTEGER :: handle, i_intval_idx, i_t, ispin, &
1828 j_intval_idx, read_write_index
1829 INTEGER, DIMENSION(2) :: i_atoms, j_atoms
1830 REAL(kind=dp) :: t1, tau
1831 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_sigma_neg_tau, mat_sigma_pos_tau
1832 TYPE(dbt_type) :: t_2c_gocc, t_2c_gvir, &
1833 t_2c_sigma_neg_tau, &
1834 t_2c_sigma_pos_tau, t_2c_w, t_3c_x_w
1835
1836 CALL timeset(routinen, handle)
1837
1838 CALL create_mat_for_sigma_c(bs_env, t_2c_gocc, t_2c_gvir, t_2c_w, t_2c_sigma_neg_tau, &
1839 t_2c_sigma_pos_tau, t_3c_x_w, &
1840 mat_sigma_neg_tau, mat_sigma_pos_tau)
1841
1842 DO i_t = 1, bs_env%num_time_freq_points
1843
1844 DO ispin = 1, bs_env%n_spin
1845
1846 t1 = m_walltime()
1847
1848 read_write_index = i_t + (ispin - 1)*bs_env%num_time_freq_points
1849
1850 ! read self-energy from restart
1851 IF (bs_env%Sigma_c_exists(i_t, ispin)) THEN
1852 CALL fm_read(bs_env%fm_work_mo(1), bs_env, bs_env%Sigma_p_name, read_write_index)
1853 CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_sigma_pos_tau(i_t, ispin)%matrix, &
1854 keep_sparsity=.false.)
1855 CALL fm_read(bs_env%fm_work_mo(1), bs_env, bs_env%Sigma_n_name, read_write_index)
1856 CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_sigma_neg_tau(i_t, ispin)%matrix, &
1857 keep_sparsity=.false.)
1858 IF (bs_env%unit_nr > 0) THEN
1859 WRITE (bs_env%unit_nr, '(T2,2A,I3,A,I3,A,F10.1,A)') 'Read Σ^c(iτ,k=0) ', &
1860 'from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1861 ', Execution time', m_walltime() - t1, ' s'
1862 END IF
1863
1864 cycle
1865
1866 END IF
1867
1868 tau = bs_env%imag_time_points(i_t)
1869
1870 CALL g_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.true., vir=.false.)
1871 CALL g_occ_vir(bs_env, tau, bs_env%fm_Gvir, ispin, occ=.false., vir=.true.)
1872
1873 ! fm G^occ, G^vir and W to local tensor
1874 CALL fm_to_local_tensor(bs_env%fm_Gocc, bs_env%mat_ao_ao%matrix, &
1875 bs_env%mat_ao_ao_tensor%matrix, t_2c_gocc, bs_env, &
1876 bs_env%atoms_i_t_group)
1877 CALL fm_to_local_tensor(bs_env%fm_Gvir, bs_env%mat_ao_ao%matrix, &
1878 bs_env%mat_ao_ao_tensor%matrix, t_2c_gvir, bs_env, &
1879 bs_env%atoms_i_t_group)
1880 CALL fm_to_local_tensor(fm_w_mic_time(i_t), bs_env%mat_RI_RI%matrix, &
1881 bs_env%mat_RI_RI_tensor%matrix, t_2c_w, bs_env, &
1882 bs_env%atoms_j_t_group)
1883
1884 ! every group has its own range of i_atoms and j_atoms; only deal with a
1885 ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
1886 DO i_intval_idx = 1, bs_env%n_intervals_i
1887 DO j_intval_idx = 1, bs_env%n_intervals_j
1888 i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
1889 j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
1890
1891 IF (bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx) .AND. &
1892 bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) THEN
1893 ! Do that only after first timestep to avoid skips due to vanishing G
1894 ! caused by gaps
1895 IF (i_t == 2) THEN
1896 bs_env%n_skip_sigma = bs_env%n_skip_sigma + 1
1897 END IF
1898 cycle
1899 END IF
1900
1901 ! 1. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
1902 ! 2. tensor operation M_Qνσ(iτ) = sum_P (νσ|P) W^MIC_QP(iτ)
1903 CALL compute_3c_and_contract_w(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_w, t_2c_w)
1904
1905 ! 3. Σ_λσ(iτ,k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) G^occ_νµ(i|τ|) for τ < 0
1906 ! (recall M_Qνσ(iτ) = M_Qνσ(-iτ) because W^MIC_PQ(iτ) = W^MIC_PQ(-iτ) )
1907 CALL contract_to_sigma(t_2c_gocc, t_3c_x_w, t_2c_sigma_neg_tau, i_atoms, j_atoms, &
1908 qs_env, bs_env, occ=.true., vir=.false., &
1909 can_skip=bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx))
1910
1911 ! Σ_λσ(iτ,k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) G^vir_νµ(i|τ|) for τ > 0
1912 CALL contract_to_sigma(t_2c_gvir, t_3c_x_w, t_2c_sigma_pos_tau, i_atoms, j_atoms, &
1913 qs_env, bs_env, occ=.false., vir=.true., &
1914 can_skip=bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx))
1915
1916 END DO ! j_atoms
1917 END DO ! i_atoms
1918
1919 ! 4. communicate data tensor t_2c_Sigma (which is local in the subgroup)
1920 ! to the global dbcsr matrix mat_Sigma_pos/neg_tau (which stores Σ for all iτ)
1921 CALL local_dbt_to_global_mat(t_2c_sigma_neg_tau, bs_env%mat_ao_ao_tensor%matrix, &
1922 mat_sigma_neg_tau(i_t, ispin)%matrix, bs_env%para_env)
1923 CALL local_dbt_to_global_mat(t_2c_sigma_pos_tau, bs_env%mat_ao_ao_tensor%matrix, &
1924 mat_sigma_pos_tau(i_t, ispin)%matrix, bs_env%para_env)
1925
1926 CALL write_matrix(mat_sigma_pos_tau(i_t, ispin)%matrix, read_write_index, &
1927 bs_env%Sigma_p_name, bs_env%fm_work_mo(1), qs_env)
1928 CALL write_matrix(mat_sigma_neg_tau(i_t, ispin)%matrix, read_write_index, &
1929 bs_env%Sigma_n_name, bs_env%fm_work_mo(1), qs_env)
1930
1931 IF (bs_env%unit_nr > 0) THEN
1932 WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F10.1,A)') &
1933 'Computed Σ^c(iτ,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1934 ', Execution time', m_walltime() - t1, ' s'
1935 END IF
1936
1937 END DO ! ispin
1938
1939 END DO ! i_t
1940
1941 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1942
1943 CALL fill_fm_sigma_c_gamma_time(fm_sigma_c_gamma_time, bs_env, &
1944 mat_sigma_pos_tau, mat_sigma_neg_tau)
1945
1946 CALL print_skipping(bs_env)
1947
1948 CALL destroy_mat_sigma_c(t_2c_gocc, t_2c_gvir, t_2c_w, t_2c_sigma_neg_tau, &
1949 t_2c_sigma_pos_tau, t_3c_x_w, fm_w_mic_time, &
1950 mat_sigma_neg_tau, mat_sigma_pos_tau)
1951
1952 CALL delete_unnecessary_files(bs_env)
1953
1954 CALL timestop(handle)
1955
1956 END SUBROUTINE get_sigma_c
1957
1958! **************************************************************************************************
1959!> \brief ...
1960!> \param bs_env ...
1961!> \param t_2c_Gocc ...
1962!> \param t_2c_Gvir ...
1963!> \param t_2c_W ...
1964!> \param t_2c_Sigma_neg_tau ...
1965!> \param t_2c_Sigma_pos_tau ...
1966!> \param t_3c_x_W ...
1967!> \param mat_Sigma_neg_tau ...
1968!> \param mat_Sigma_pos_tau ...
1969! **************************************************************************************************
1970 SUBROUTINE create_mat_for_sigma_c(bs_env, t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
1971 t_2c_Sigma_pos_tau, t_3c_x_W, &
1972 mat_Sigma_neg_tau, mat_Sigma_pos_tau)
1973
1974 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1975 TYPE(dbt_type) :: t_2c_gocc, t_2c_gvir, t_2c_w, &
1976 t_2c_sigma_neg_tau, &
1977 t_2c_sigma_pos_tau, t_3c_x_w
1978 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_sigma_neg_tau, mat_sigma_pos_tau
1979
1980 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_mat_for_Sigma_c'
1981
1982 INTEGER :: handle, i_t, ispin
1983
1984 CALL timeset(routinen, handle)
1985
1986 CALL dbt_create(bs_env%t_G, t_2c_gocc)
1987 CALL dbt_create(bs_env%t_G, t_2c_gvir)
1988 CALL dbt_create(bs_env%t_W, t_2c_w)
1989 CALL dbt_create(bs_env%t_G, t_2c_sigma_neg_tau)
1990 CALL dbt_create(bs_env%t_G, t_2c_sigma_pos_tau)
1991 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_w)
1992
1993 NULLIFY (mat_sigma_neg_tau, mat_sigma_pos_tau)
1994 ALLOCATE (mat_sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1995 ALLOCATE (mat_sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1996
1997 DO ispin = 1, bs_env%n_spin
1998 DO i_t = 1, bs_env%num_time_freq_points
1999 ALLOCATE (mat_sigma_neg_tau(i_t, ispin)%matrix)
2000 ALLOCATE (mat_sigma_pos_tau(i_t, ispin)%matrix)
2001 CALL dbcsr_create(mat_sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
2002 CALL dbcsr_create(mat_sigma_pos_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
2003 END DO
2004 END DO
2005
2006 CALL timestop(handle)
2007
2008 END SUBROUTINE create_mat_for_sigma_c
2009
2010! **************************************************************************************************
2011!> \brief ...
2012!> \param qs_env ...
2013!> \param bs_env ...
2014!> \param i_atoms ...
2015!> \param j_atoms ...
2016!> \param t_3c_x_W ...
2017!> \param t_2c_W ...
2018! **************************************************************************************************
2019 SUBROUTINE compute_3c_and_contract_w(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_W, t_2c_W)
2020
2021 TYPE(qs_environment_type), POINTER :: qs_env
2022 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2023 INTEGER, DIMENSION(2) :: i_atoms, j_atoms
2024 TYPE(dbt_type) :: t_3c_x_w, t_2c_w
2025
2026 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_3c_and_contract_W'
2027
2028 INTEGER :: handle, ri_intval_idx
2029 INTEGER(KIND=int_8) :: flop
2030 INTEGER, DIMENSION(2) :: bounds_p, bounds_q, ri_atoms
2031 INTEGER, DIMENSION(2, 2) :: bounds_ao
2032 TYPE(dbt_type) :: t_3c_for_w, t_3c_x_w_tmp
2033
2034 CALL timeset(routinen, handle)
2035
2036 CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_w_tmp)
2037 CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_for_w)
2038
2039 ! final layout will be: M_Qνσ(iτ) = sum_P (P|νσ) W^MIC_QP(iτ)
2040 ! Bounds:
2041 ! "AO"
2042 ! -> ν (AO_1 in compute_3c_integrals) bounds from i_atoms and sparse in σ and P
2043 ! -> σ (AO_2 in compute_3c_integrals) sparse in ν and P
2044 ! Q bounds from j_atoms
2045 ! P bounds from inner loop indices and sparse in ν and σ
2046
2047 bounds_q(1:2) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
2048 bs_env%i_RI_end_from_atom(j_atoms(2))]
2049
2050 DO ri_intval_idx = 1, bs_env%n_intervals_inner_loop_atoms
2051 ri_atoms = bs_env%inner_loop_atom_intervals(1:2, ri_intval_idx)
2052
2053 CALL get_bounds_from_atoms(bounds_p, i_atoms, [1, bs_env%n_atom], &
2054 bs_env%min_RI_idx_from_AO_AO_atom, &
2055 bs_env%max_RI_idx_from_AO_AO_atom, &
2056 atoms_3=ri_atoms, &
2057 indices_3_start=bs_env%i_RI_start_from_atom, &
2058 indices_3_end=bs_env%i_RI_end_from_atom)
2059
2060 ! σ
2061 CALL get_bounds_from_atoms(bounds_ao(:, 2), ri_atoms, i_atoms, &
2062 bs_env%min_AO_idx_from_RI_AO_atom, &
2063 bs_env%max_AO_idx_from_RI_AO_atom)
2064 ! ν
2065 CALL get_bounds_from_atoms(bounds_ao(:, 1), ri_atoms, [1, bs_env%n_atom], &
2066 bs_env%min_AO_idx_from_RI_AO_atom, &
2067 bs_env%max_AO_idx_from_RI_AO_atom, &
2068 atoms_3=i_atoms, &
2069 indices_3_start=bs_env%i_ao_start_from_atom, &
2070 indices_3_end=bs_env%i_ao_end_from_atom)
2071
2072 IF (bounds_p(1) > bounds_p(2) .OR. bounds_ao(1, 2) > bounds_ao(2, 2)) THEN
2073 cycle
2074 END IF
2075
2076 ! 1. compute 3-center integrals (P|µν) ("|": truncated Coulomb operator)
2077 CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_w, &
2078 atoms_ao_1=i_atoms, atoms_ri=ri_atoms)
2079
2080 ! 2. tensor operation M_Qνσ(iτ) = sum_P W^MIC_QP(iτ) (P|νσ)
2081 CALL dbt_contract(alpha=1.0_dp, &
2082 tensor_1=t_2c_w, &
2083 tensor_2=t_3c_for_w, &
2084 beta=1.0_dp, &
2085 tensor_3=t_3c_x_w_tmp, &
2086 contract_1=[2], notcontract_1=[1], map_1=[1], &
2087 contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3], &
2088 bounds_1=bounds_p, &
2089 bounds_2=bounds_q, &
2090 bounds_3=bounds_ao, &
2091 flop=flop, &
2092 move_data=.false., &
2093 filter_eps=bs_env%eps_filter)
2094
2095 END DO ! RI_atoms
2096
2097 ! 3. reorder tensor
2098 CALL dbt_copy(t_3c_x_w_tmp, t_3c_x_w, order=[1, 2, 3], move_data=.true.)
2099
2100 CALL dbt_destroy(t_3c_x_w_tmp)
2101 CALL dbt_destroy(t_3c_for_w)
2102
2103 CALL timestop(handle)
2104
2105 END SUBROUTINE compute_3c_and_contract_w
2106
2107! **************************************************************************************************
2108!> \brief ...
2109!> \param t_2c_G ...
2110!> \param t_3c_x_W ...
2111!> \param t_2c_Sigma ...
2112!> \param i_atoms ...
2113!> \param j_atoms ...
2114!> \param qs_env ...
2115!> \param bs_env ...
2116!> \param occ ...
2117!> \param vir ...
2118!> \param can_skip ...
2119! **************************************************************************************************
2120 SUBROUTINE contract_to_sigma(t_2c_G, t_3c_x_W, t_2c_Sigma, i_atoms, j_atoms, qs_env, bs_env, &
2121 occ, vir, can_skip)
2122 TYPE(dbt_type) :: t_2c_g, t_3c_x_w, t_2c_sigma
2123 INTEGER, DIMENSION(2) :: i_atoms, j_atoms
2124 TYPE(qs_environment_type), POINTER :: qs_env
2125 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2126 LOGICAL :: occ, vir
2127 LOGICAL, OPTIONAL :: can_skip
2128
2129 CHARACTER(LEN=*), PARAMETER :: routinen = 'contract_to_Sigma'
2130
2131 INTEGER :: handle, inner_loop_atoms_interval_index
2132 INTEGER(KIND=int_8) :: flop
2133 INTEGER, DIMENSION(2) :: bounds_lambda, bounds_mu, bounds_nu, &
2134 bounds_sigma, il_atoms
2135 INTEGER, DIMENSION(2, 2) :: bounds_comb
2136 REAL(kind=dp) :: sign_sigma
2137 TYPE(dbt_type) :: t_3c_for_g, t_3c_x_g, t_3c_x_g_2
2138
2139 CALL timeset(routinen, handle)
2140
2141 cpassert(occ .EQV. (.NOT. vir))
2142 IF (occ) sign_sigma = -1.0_dp
2143 IF (vir) sign_sigma = 1.0_dp
2144
2145 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_g)
2146 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_g)
2147 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_g_2)
2148
2149 ! Here, in the first step e.g., is computed: N_Qλν = sum_µ (Qλ|µ) G_νµ
2150 ! Afterwards e.g., is computed: Σ_λσ = sum_νQ M_Qνσ N_Qνλ (after reordering)
2151 ! Bounds:
2152 ! "comb" (combined index)
2153 ! -> Q bounds from j_atoms and sparse in λ
2154 ! -> λ (AO_1 in compute_3c_integrals) sparse in Q and µ
2155 ! µ (AO_2 in compute_3c_integrals) bounds from inner loop "IL" indices and sparse in Q and λ
2156 ! ν bounds from i_atoms
2157 ! σ sparse in ν
2158
2159 ! ν
2160 bounds_nu(1:2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
2161 bs_env%i_ao_end_from_atom(i_atoms(2))]
2162
2163 DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
2164 il_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
2165
2166 ! µ
2167 CALL get_bounds_from_atoms(bounds_mu, j_atoms, [1, bs_env%n_atom], &
2168 bs_env%min_AO_idx_from_RI_AO_atom, &
2169 bs_env%max_AO_idx_from_RI_AO_atom, &
2170 atoms_3=il_atoms, &
2171 indices_3_start=bs_env%i_ao_start_from_atom, &
2172 indices_3_end=bs_env%i_ao_end_from_atom)
2173
2174 ! Q
2175 CALL get_bounds_from_atoms(bounds_comb(:, 1), il_atoms, [1, bs_env%n_atom], &
2176 bs_env%min_RI_idx_from_AO_AO_atom, &
2177 bs_env%max_RI_idx_from_AO_AO_atom, &
2178 atoms_3=j_atoms, &
2179 indices_3_start=bs_env%i_RI_start_from_atom, &
2180 indices_3_end=bs_env%i_RI_end_from_atom)
2181
2182 ! λ
2183 CALL get_bounds_from_atoms(bounds_comb(:, 2), j_atoms, il_atoms, &
2184 bs_env%min_AO_idx_from_RI_AO_atom, &
2185 bs_env%max_AO_idx_from_RI_AO_atom)
2186
2187 IF (bounds_mu(1) > bounds_mu(2) .OR. bounds_comb(1, 1) > bounds_comb(2, 1) .OR. &
2188 bounds_comb(1, 2) > bounds_comb(2, 2)) THEN
2189 cycle
2190 END IF
2191
2192 CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_g, &
2193 atoms_ri=j_atoms, atoms_ao_2=il_atoms)
2194
2195 CALL dbt_contract(alpha=1.0_dp, &
2196 tensor_1=t_2c_g, &
2197 tensor_2=t_3c_for_g, &
2198 beta=1.0_dp, &
2199 tensor_3=t_3c_x_g, &
2200 contract_1=[2], notcontract_1=[1], map_1=[3], &
2201 contract_2=[3], notcontract_2=[1, 2], map_2=[1, 2], &
2202 bounds_1=bounds_mu, &
2203 bounds_2=bounds_nu, &
2204 bounds_3=bounds_comb, &
2205 flop=flop, &
2206 move_data=.false., &
2207 filter_eps=bs_env%eps_filter)
2208 END DO ! IL_atoms
2209
2210 ! Reordering: N_Qλν -> N_Qνλ
2211 CALL dbt_copy(t_3c_x_g, t_3c_x_g_2, order=[1, 3, 2], move_data=.true.)
2212
2213 ! Here, the last contraction is done, e.g., Σ_λσ = sum_νQ M_Qνσ N_Qνλ
2214 ! Bounds as above, new "comb" with upper ingredients
2215 bounds_comb(1:2, 1) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
2216 bs_env%i_RI_end_from_atom(j_atoms(2))]
2217 bounds_comb(1:2, 2) = bounds_nu(1:2)
2218
2219 CALL get_bounds_from_atoms(bounds_lambda, j_atoms, [1, bs_env%n_atom], &
2220 bs_env%min_AO_idx_from_RI_AO_atom, &
2221 bs_env%max_AO_idx_from_RI_AO_atom)
2222 CALL get_bounds_from_atoms(bounds_sigma, [1, bs_env%n_atom], i_atoms, &
2223 bs_env%min_AO_idx_from_RI_AO_atom, &
2224 bs_env%max_AO_idx_from_RI_AO_atom)
2225
2226 IF (bounds_sigma(1) > bounds_sigma(2) .OR. bounds_lambda(1) > bounds_lambda(2)) THEN
2227 flop = 0_int_8
2228 ELSE
2229 CALL dbt_contract(alpha=sign_sigma, &
2230 tensor_1=t_3c_x_w, &
2231 tensor_2=t_3c_x_g_2, &
2232 beta=1.0_dp, &
2233 tensor_3=t_2c_sigma, &
2234 contract_1=[1, 2], notcontract_1=[3], map_1=[1], &
2235 contract_2=[1, 2], notcontract_2=[3], map_2=[2], &
2236 bounds_1=bounds_comb, &
2237 bounds_2=bounds_sigma, &
2238 bounds_3=bounds_lambda, &
2239 filter_eps=bs_env%eps_filter, move_data=.false., flop=flop)
2240 END IF
2241
2242 IF (PRESENT(can_skip)) THEN
2243 IF (flop == 0_int_8) can_skip = .true.
2244 END IF
2245
2246 CALL dbt_destroy(t_3c_for_g)
2247 CALL dbt_destroy(t_3c_x_g)
2248 CALL dbt_destroy(t_3c_x_g_2)
2249
2250 CALL timestop(handle)
2251
2252 END SUBROUTINE contract_to_sigma
2253
2254! **************************************************************************************************
2255!> \brief ...
2256!> \param fm_Sigma_c_Gamma_time ...
2257!> \param bs_env ...
2258!> \param mat_Sigma_pos_tau ...
2259!> \param mat_Sigma_neg_tau ...
2260! **************************************************************************************************
2261 SUBROUTINE fill_fm_sigma_c_gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
2262 mat_Sigma_pos_tau, mat_Sigma_neg_tau)
2263
2264 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_sigma_c_gamma_time
2265 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2266 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_sigma_pos_tau, mat_sigma_neg_tau
2267
2268 CHARACTER(LEN=*), PARAMETER :: routinen = 'fill_fm_Sigma_c_Gamma_time'
2269
2270 INTEGER :: handle, i_t, ispin, pos_neg
2271
2272 CALL timeset(routinen, handle)
2273
2274 ALLOCATE (fm_sigma_c_gamma_time(bs_env%num_time_freq_points, 2, bs_env%n_spin))
2275 DO ispin = 1, bs_env%n_spin
2276 DO i_t = 1, bs_env%num_time_freq_points
2277 DO pos_neg = 1, 2
2278 CALL cp_fm_create(fm_sigma_c_gamma_time(i_t, pos_neg, ispin), &
2279 bs_env%fm_s_Gamma%matrix_struct)
2280 END DO
2281 CALL copy_dbcsr_to_fm(mat_sigma_pos_tau(i_t, ispin)%matrix, &
2282 fm_sigma_c_gamma_time(i_t, 1, ispin))
2283 CALL copy_dbcsr_to_fm(mat_sigma_neg_tau(i_t, ispin)%matrix, &
2284 fm_sigma_c_gamma_time(i_t, 2, ispin))
2285 END DO
2286 END DO
2287
2288 CALL timestop(handle)
2289
2290 END SUBROUTINE fill_fm_sigma_c_gamma_time
2291
2292! **************************************************************************************************
2293!> \brief ...
2294!> \param bs_env ...
2295! **************************************************************************************************
2296 SUBROUTINE print_skipping(bs_env)
2297
2298 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2299
2300 CHARACTER(LEN=*), PARAMETER :: routinen = 'print_skipping'
2301
2302 INTEGER :: handle, n_pairs
2303
2304 CALL timeset(routinen, handle)
2305
2306 n_pairs = bs_env%n_intervals_i*bs_env%n_intervals_j*bs_env%n_spin
2307
2308 CALL bs_env%para_env_tensor%sum(bs_env%n_skip_sigma)
2309 CALL bs_env%para_env_tensor%sum(bs_env%n_skip_chi)
2310 CALL bs_env%para_env_tensor%sum(n_pairs)
2311
2312 IF (bs_env%unit_nr > 0) THEN
2313 WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
2314 'Sparsity of Σ^c(iτ,k=0): Percentage of skipped atom pairs:', &
2315 REAL(100*bs_env%n_skip_sigma, kind=dp)/real(n_pairs, kind=dp), ' %'
2316 WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
2317 'Sparsity of χ(iτ,k=0): Percentage of skipped atom pairs:', &
2318 REAL(100*bs_env%n_skip_chi, kind=dp)/real(n_pairs, kind=dp), ' %'
2319 END IF
2320
2321 CALL timestop(handle)
2322
2323 END SUBROUTINE print_skipping
2324
2325! **************************************************************************************************
2326!> \brief ...
2327!> \param t_2c_Gocc ...
2328!> \param t_2c_Gvir ...
2329!> \param t_2c_W ...
2330!> \param t_2c_Sigma_neg_tau ...
2331!> \param t_2c_Sigma_pos_tau ...
2332!> \param t_3c_x_W ...
2333!> \param fm_W_MIC_time ...
2334!> \param mat_Sigma_neg_tau ...
2335!> \param mat_Sigma_pos_tau ...
2336! **************************************************************************************************
2337 SUBROUTINE destroy_mat_sigma_c(t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
2338 t_2c_Sigma_pos_tau, t_3c_x_W, fm_W_MIC_time, &
2339 mat_Sigma_neg_tau, mat_Sigma_pos_tau)
2340
2341 TYPE(dbt_type) :: t_2c_gocc, t_2c_gvir, t_2c_w, &
2342 t_2c_sigma_neg_tau, &
2343 t_2c_sigma_pos_tau, t_3c_x_w
2344 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
2345 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_sigma_neg_tau, mat_sigma_pos_tau
2346
2347 CHARACTER(LEN=*), PARAMETER :: routinen = 'destroy_mat_Sigma_c'
2348
2349 INTEGER :: handle
2350
2351 CALL timeset(routinen, handle)
2352
2353 CALL dbt_destroy(t_2c_gocc)
2354 CALL dbt_destroy(t_2c_gvir)
2355 CALL dbt_destroy(t_2c_w)
2356 CALL dbt_destroy(t_2c_sigma_neg_tau)
2357 CALL dbt_destroy(t_2c_sigma_pos_tau)
2358 CALL dbt_destroy(t_3c_x_w)
2359 CALL cp_fm_release(fm_w_mic_time)
2360 CALL dbcsr_deallocate_matrix_set(mat_sigma_neg_tau)
2361 CALL dbcsr_deallocate_matrix_set(mat_sigma_pos_tau)
2362
2363 CALL timestop(handle)
2364
2365 END SUBROUTINE destroy_mat_sigma_c
2366
2367! **************************************************************************************************
2368!> \brief ...
2369!> \param bs_env ...
2370! **************************************************************************************************
2371 SUBROUTINE delete_unnecessary_files(bs_env)
2372 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2373
2374 CHARACTER(LEN=*), PARAMETER :: routinen = 'delete_unnecessary_files'
2375
2376 CHARACTER(LEN=default_path_length) :: f_chi, f_w_t, prefix
2377 INTEGER :: handle, i_t
2378
2379 CALL timeset(routinen, handle)
2380
2381 prefix = bs_env%prefix
2382
2383 DO i_t = 1, bs_env%num_time_freq_points
2384
2385 IF (i_t < 10) THEN
2386 WRITE (f_chi, '(3A,I1,A)') trim(prefix), bs_env%chi_name, "_00", i_t, ".matrix"
2387 WRITE (f_w_t, '(3A,I1,A)') trim(prefix), bs_env%W_time_name, "_00", i_t, ".matrix"
2388 ELSE IF (i_t < 100) THEN
2389 WRITE (f_chi, '(3A,I2,A)') trim(prefix), bs_env%chi_name, "_0", i_t, ".matrix"
2390 WRITE (f_w_t, '(3A,I2,A)') trim(prefix), bs_env%W_time_name, "_0", i_t, ".matrix"
2391 ELSE
2392 cpabort('Please implement more than 99 time/frequency points.')
2393 END IF
2394
2395 CALL safe_delete(f_chi, bs_env)
2396 CALL safe_delete(f_w_t, bs_env)
2397
2398 END DO
2399
2400 CALL timestop(handle)
2401
2402 END SUBROUTINE delete_unnecessary_files
2403
2404! **************************************************************************************************
2405!> \brief ...
2406!> \param filename ...
2407!> \param bs_env ...
2408! **************************************************************************************************
2409 SUBROUTINE safe_delete(filename, bs_env)
2410 CHARACTER(LEN=*) :: filename
2411 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2412
2413 CHARACTER(LEN=*), PARAMETER :: routinen = 'safe_delete'
2414
2415 INTEGER :: handle
2416 LOGICAL :: file_exists
2417
2418 CALL timeset(routinen, handle)
2419
2420 IF (bs_env%para_env%mepos == 0) THEN
2421
2422 INQUIRE (file=trim(filename), exist=file_exists)
2423 IF (file_exists) CALL mp_file_delete(trim(filename))
2424
2425 END IF
2426
2427 CALL timestop(handle)
2428
2429 END SUBROUTINE safe_delete
2430
2431! **************************************************************************************************
2432!> \brief ...
2433!> \param bs_env ...
2434!> \param qs_env ...
2435!> \param fm_Sigma_x_Gamma ...
2436!> \param fm_Sigma_c_Gamma_time ...
2437! **************************************************************************************************
2438 SUBROUTINE compute_qp_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
2439
2440 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2441 TYPE(qs_environment_type), POINTER :: qs_env
2442 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_sigma_x_gamma
2443 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_sigma_c_gamma_time
2444
2445 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_QP_energies'
2446
2447 INTEGER :: handle, ikp, ispin, j_t
2448 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sigma_x_ikp_n, v_xc_ikp_n
2449 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sigma_c_ikp_n_freq, sigma_c_ikp_n_time
2450 TYPE(cp_cfm_type) :: cfm_ks_ikp, cfm_mos_ikp, cfm_s_ikp, &
2451 cfm_sigma_x_ikp, cfm_work_ikp
2452
2453 CALL timeset(routinen, handle)
2454
2455 CALL cp_cfm_create(cfm_mos_ikp, bs_env%fm_s_Gamma%matrix_struct)
2456 CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_s_Gamma%matrix_struct)
2457 ! JW TODO: fully distribute these arrays at given time; also eigenvalues in bs_env
2458 ALLOCATE (v_xc_ikp_n(bs_env%n_ao), sigma_x_ikp_n(bs_env%n_ao))
2459 ALLOCATE (sigma_c_ikp_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
2460 ALLOCATE (sigma_c_ikp_n_freq(bs_env%n_ao, bs_env%num_time_freq_points, 2))
2461
2462 DO ispin = 1, bs_env%n_spin
2463
2464 DO ikp = 1, bs_env%nkp_bs_and_DOS
2465
2466 ! 1. get H^KS_µν(k_i) from H^KS_µν(k=0)
2467 CALL cfm_ikp_from_fm_gamma(cfm_ks_ikp, bs_env%fm_ks_Gamma(ispin), &
2468 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2469
2470 ! 2. get S_µν(k_i) from S_µν(k=0)
2471 CALL cfm_ikp_from_fm_gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
2472 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2473
2474 ! 3. Diagonalize (Roothaan-Hall): H_KS(k_i)*C(k_i) = S(k_i)*C(k_i)*ϵ(k_i)
2475 CALL cp_cfm_geeig(cfm_ks_ikp, cfm_s_ikp, cfm_mos_ikp, &
2476 bs_env%eigenval_scf(:, ikp, ispin), cfm_work_ikp)
2477
2478 ! 4. V^xc_µν(k=0) -> V^xc_µν(k_i) -> V^xc_nn(k_i)
2479 CALL to_ikp_and_mo(v_xc_ikp_n, bs_env%fm_V_xc_Gamma(ispin), &
2480 ikp, qs_env, bs_env, cfm_mos_ikp)
2481
2482 ! 5. Σ^x_µν(k=0) -> Σ^x_µν(k_i) -> Σ^x_nn(k_i)
2483 CALL to_ikp_and_mo(sigma_x_ikp_n, fm_sigma_x_gamma(ispin), &
2484 ikp, qs_env, bs_env, cfm_mos_ikp)
2485
2486 ! 6. Σ^c_µν(k=0,+/-i|τ_j|) -> Σ^c_µν(k_i,+/-i|τ_j|) -> Σ^c_nn(k_i,+/-i|τ_j|)
2487 DO j_t = 1, bs_env%num_time_freq_points
2488 CALL to_ikp_and_mo(sigma_c_ikp_n_time(:, j_t, 1), &
2489 fm_sigma_c_gamma_time(j_t, 1, ispin), &
2490 ikp, qs_env, bs_env, cfm_mos_ikp)
2491 CALL to_ikp_and_mo(sigma_c_ikp_n_time(:, j_t, 2), &
2492 fm_sigma_c_gamma_time(j_t, 2, ispin), &
2493 ikp, qs_env, bs_env, cfm_mos_ikp)
2494 END DO
2495
2496 ! 7. Σ^c_nn(k_i,iτ) -> Σ^c_nn(k_i,iω)
2497 CALL time_to_freq(bs_env, sigma_c_ikp_n_time, sigma_c_ikp_n_freq, ispin)
2498
2499 ! 8. Analytic continuation Σ^c_nn(k_i,iω) -> Σ^c_nn(k_i,ϵ) and
2500 ! ϵ_nk_i^GW = ϵ_nk_i^DFT + Σ^c_nn(k_i,ϵ) + Σ^x_nn(k_i) - v^xc_nn(k_i)
2501 CALL analyt_conti_and_print(bs_env, sigma_c_ikp_n_freq, sigma_x_ikp_n, v_xc_ikp_n, &
2502 bs_env%eigenval_scf(:, ikp, ispin), ikp, ispin)
2503
2504 END DO ! ikp_DOS
2505
2506 END DO ! ispin
2507
2508 CALL get_all_vbm_cbm_bandgaps(bs_env)
2509
2510 ! Σ^x is releases here in case of G0W0
2511 IF (bs_env%gw_flavour == g0w0) CALL cp_fm_release(fm_sigma_x_gamma)
2512 CALL cp_fm_release(fm_sigma_c_gamma_time)
2513 CALL cp_cfm_release(cfm_ks_ikp)
2514 CALL cp_cfm_release(cfm_s_ikp)
2515 CALL cp_cfm_release(cfm_mos_ikp)
2516 CALL cp_cfm_release(cfm_work_ikp)
2517 CALL cp_cfm_release(cfm_sigma_x_ikp)
2518
2519 CALL timestop(handle)
2520
2521 END SUBROUTINE compute_qp_energies
2522
2523! **************************************************************************************************
2524!> \brief ...
2525!> \param array_ikp_n ...
2526!> \param fm_Gamma ...
2527!> \param ikp ...
2528!> \param qs_env ...
2529!> \param bs_env ...
2530!> \param cfm_mos_ikp ...
2531! **************************************************************************************************
2532 SUBROUTINE to_ikp_and_mo(array_ikp_n, fm_Gamma, ikp, qs_env, bs_env, cfm_mos_ikp)
2533
2534 REAL(kind=dp), DIMENSION(:) :: array_ikp_n
2535 TYPE(cp_fm_type) :: fm_gamma
2536 INTEGER :: ikp
2537 TYPE(qs_environment_type), POINTER :: qs_env
2538 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2539 TYPE(cp_cfm_type) :: cfm_mos_ikp
2540
2541 CHARACTER(LEN=*), PARAMETER :: routinen = 'to_ikp_and_mo'
2542
2543 INTEGER :: handle
2544 TYPE(cp_fm_type) :: fm_ikp_mo_re
2545
2546 CALL timeset(routinen, handle)
2547
2548 CALL cp_fm_create(fm_ikp_mo_re, fm_gamma%matrix_struct)
2549
2550 CALL fm_gamma_ao_to_cfm_ikp_mo(fm_gamma, fm_ikp_mo_re, ikp, qs_env, bs_env, cfm_mos_ikp)
2551
2552 CALL cp_fm_get_diag(fm_ikp_mo_re, array_ikp_n)
2553
2554 CALL cp_fm_release(fm_ikp_mo_re)
2555
2556 CALL timestop(handle)
2557
2558 END SUBROUTINE to_ikp_and_mo
2559
2560! **************************************************************************************************
2561!> \brief ...
2562!> \param fm_Gamma ...
2563!> \param fm_ikp_mo_re ...
2564!> \param ikp ...
2565!> \param qs_env ...
2566!> \param bs_env ...
2567!> \param cfm_mos_ikp ...
2568! **************************************************************************************************
2569 SUBROUTINE fm_gamma_ao_to_cfm_ikp_mo(fm_Gamma, fm_ikp_mo_re, ikp, qs_env, bs_env, cfm_mos_ikp)
2570 TYPE(cp_fm_type) :: fm_gamma, fm_ikp_mo_re
2571 INTEGER :: ikp
2572 TYPE(qs_environment_type), POINTER :: qs_env
2573 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2574 TYPE(cp_cfm_type) :: cfm_mos_ikp
2575
2576 CHARACTER(LEN=*), PARAMETER :: routinen = 'fm_Gamma_ao_to_cfm_ikp_mo'
2577
2578 INTEGER :: handle, nmo
2579 TYPE(cp_cfm_type) :: cfm_ikp_ao, cfm_ikp_mo, cfm_tmp
2580
2581 CALL timeset(routinen, handle)
2582
2583 CALL cp_cfm_create(cfm_ikp_ao, fm_gamma%matrix_struct)
2584 CALL cp_cfm_create(cfm_ikp_mo, fm_gamma%matrix_struct)
2585 CALL cp_cfm_create(cfm_tmp, fm_gamma%matrix_struct)
2586
2587 ! get cfm_µν(k_i) from fm_µν(k=0)
2588 CALL cfm_ikp_from_fm_gamma(cfm_ikp_ao, fm_gamma, ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2589
2590 nmo = bs_env%n_ao
2591 CALL parallel_gemm('N', 'N', nmo, nmo, nmo, z_one, cfm_ikp_ao, cfm_mos_ikp, z_zero, cfm_tmp)
2592 CALL parallel_gemm('C', 'N', nmo, nmo, nmo, z_one, cfm_mos_ikp, cfm_tmp, z_zero, cfm_ikp_mo)
2593
2594 CALL cp_cfm_to_fm(cfm_ikp_mo, fm_ikp_mo_re)
2595
2596 CALL cp_cfm_release(cfm_ikp_mo)
2597 CALL cp_cfm_release(cfm_ikp_ao)
2598 CALL cp_cfm_release(cfm_tmp)
2599
2600 CALL timestop(handle)
2601
2602 END SUBROUTINE fm_gamma_ao_to_cfm_ikp_mo
2603
2604! **************************************************************************************************
2605!> \brief Computes bounds (AO or RI) for given atom intervals atoms_1 and atoms_2 from indices_min
2606!> and indices_max and returns them in bounds_out.
2607!> In case, atoms_3 and indices_3 are given, the bounds are computed as the intersection
2608!> \param bounds_out Bounds to be computed
2609!> \param atoms_1 First atom interval
2610!> \param atoms_2 Second atom interval
2611!> \param indices_min Minimum indices for each atom pair (typically from bs_env,
2612!> computed in get_i_j_atom_ranges in gw_utils.F, e.g. bs_env%min_RI_idx_from_AO_AO_atom)
2613!> \param indices_max Maximum indices for each atom pair (typically from bs_env,
2614!> computed in get_i_j_atom_ranges in gw_utils.F)
2615!> \param atoms_3 (Optional) Third atom interval for intersection
2616!> \param indices_3_start (Optional) Indices for third atom interval for intersection
2617!> \param indices_3_end (Optional) Indices for third atom interval for intersection
2618! **************************************************************************************************
2619 SUBROUTINE get_bounds_from_atoms(bounds_out, atoms_1, atoms_2, indices_min, indices_max, &
2620 atoms_3, indices_3_start, indices_3_end)
2621
2622 INTEGER, DIMENSION(2), INTENT(OUT) :: bounds_out
2623 INTEGER, DIMENSION(2), INTENT(IN) :: atoms_1, atoms_2
2624 INTEGER, DIMENSION(:, :) :: indices_min, indices_max
2625 INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: atoms_3
2626 INTEGER, DIMENSION(:), OPTIONAL :: indices_3_start, indices_3_end
2627
2628 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_bounds_from_atoms'
2629
2630 INTEGER :: handle, i_at, j_at
2631
2632 CALL timeset(routinen, handle)
2633 bounds_out(1) = huge(0)
2634 bounds_out(2) = -1
2635 !Loop over all atoms in the two intervals and find min/max indices
2636 DO i_at = atoms_1(1), atoms_1(2)
2637 DO j_at = atoms_2(1), atoms_2(2)
2638 bounds_out(1) = min(bounds_out(1), indices_min(i_at, j_at))
2639 bounds_out(2) = max(bounds_out(2), indices_max(i_at, j_at))
2640 END DO
2641 END DO
2642
2643 IF (PRESENT(atoms_3) .AND. PRESENT(indices_3_start) .AND. PRESENT(indices_3_end)) THEN
2644 bounds_out(1) = max(bounds_out(1), indices_3_start(atoms_3(1)))
2645 bounds_out(2) = min(bounds_out(2), indices_3_end(atoms_3(2)))
2646 END IF
2647
2648 CALL timestop(handle)
2649
2650 END SUBROUTINE get_bounds_from_atoms
2651
2652END MODULE gw_large_cell_gamma
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:161
cell_types::pbc
Definition cell_types.F:104
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:56
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
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
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
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::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
constants_operator
constants for the different operators of the 2c-integrals
Definition constants_operator.F:14
constants_operator::operator_coulomb
integer, parameter, public operator_coulomb
Definition constants_operator.F:22
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_uplo_to_full
subroutine, public cp_cfm_uplo_to_full(matrix, workspace, uplo)
...
Definition cp_cfm_basic_linalg.F:1350
cp_cfm_cholesky
various cholesky decomposition related routines
Definition cp_cfm_cholesky.F:13
cp_cfm_cholesky::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_cholesky.F:53
cp_cfm_cholesky::cp_cfm_cholesky_invert
subroutine, public cp_cfm_cholesky_invert(matrix, n, info_out)
Used to replace Cholesky decomposition by the inverse.
Definition cp_cfm_cholesky.F:124
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work, lowest_subset)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:329
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:844
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:631
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:792
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
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_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_reserve_all_blocks
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
Definition cp_dbcsr_contrib.F:295
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:504
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_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:1144
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_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition cp_fm_types.F:2132
cp_fm_types::cp_fm_read_unformatted
subroutine, public cp_fm_read_unformatted(fm, unit)
...
Definition cp_fm_types.F:2363
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_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
gw_communication
Definition gw_communication.F:13
gw_communication::fm_to_local_tensor
subroutine, public fm_to_local_tensor(fm_global, mat_global, mat_local, tensor, bs_env, atom_ranges)
...
Definition gw_communication.F:67
gw_communication::local_dbt_to_global_mat
subroutine, public local_dbt_to_global_mat(tensor, mat_tensor, mat_global, para_env)
...
Definition gw_communication.F:110
gw_large_cell_gamma
Routines from paper [Graml2024].
Definition gw_large_cell_gamma.F:15
gw_large_cell_gamma::compute_fm_chi_gamma_freq
subroutine, public compute_fm_chi_gamma_freq(bs_env, fm_chi_gamma_freq, j_w, mat_chi_gamma_tau)
...
Definition gw_large_cell_gamma.F:1339
gw_large_cell_gamma::delete_unnecessary_files
subroutine, public delete_unnecessary_files(bs_env)
...
Definition gw_large_cell_gamma.F:2372
gw_large_cell_gamma::fill_fm_sigma_c_gamma_time
subroutine, public fill_fm_sigma_c_gamma_time(fm_sigma_c_gamma_time, bs_env, mat_sigma_pos_tau, mat_sigma_neg_tau)
...
Definition gw_large_cell_gamma.F:2263
gw_large_cell_gamma::fm_write
subroutine, public fm_write(fm, matrix_index, matrix_name, qs_env)
...
Definition gw_large_cell_gamma.F:499
gw_large_cell_gamma::gw_calc_large_cell_gamma
subroutine, public gw_calc_large_cell_gamma(qs_env, bs_env)
Perform GW band structure calculation.
Definition gw_large_cell_gamma.F:118
gw_large_cell_gamma::compute_qp_energies
subroutine, public compute_qp_energies(bs_env, qs_env, fm_sigma_x_gamma, fm_sigma_c_gamma_time)
...
Definition gw_large_cell_gamma.F:2439
gw_large_cell_gamma::multiply_fm_w_mic_time_with_minv_gamma
subroutine, public multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_w_mic_time)
...
Definition gw_large_cell_gamma.F:1642
gw_large_cell_gamma::fm_read
subroutine, public fm_read(fm, bs_env, mat_name, idx)
...
Definition gw_large_cell_gamma.F:347
gw_large_cell_gamma::get_w_mic
subroutine, public get_w_mic(bs_env, qs_env, mat_chi_gamma_tau, fm_w_mic_time)
...
Definition gw_large_cell_gamma.F:793
gw_large_cell_gamma::write_matrix
subroutine, public write_matrix(matrix, matrix_index, matrix_name, fm, qs_env)
...
Definition gw_large_cell_gamma.F:469
gw_large_cell_gamma::g_occ_vir
subroutine, public g_occ_vir(bs_env, tau, fm_g_gamma, ispin, occ, vir)
...
Definition gw_large_cell_gamma.F:552
gw_large_cell_gamma::create_fm_w_mic_time
subroutine, public create_fm_w_mic_time(bs_env, fm_w_mic_time)
...
Definition gw_large_cell_gamma.F:1580
gw_large_cell_gamma::compute_3c_integrals
subroutine, public compute_3c_integrals(qs_env, bs_env, t_3c, atoms_ao_1, atoms_ao_2, atoms_ri)
...
Definition gw_large_cell_gamma.F:619
gw_large_cell_gamma::fourier_transform_w_to_t
subroutine, public fourier_transform_w_to_t(bs_env, fm_w_mic_time, fm_w_mic_freq_j, j_w)
...
Definition gw_large_cell_gamma.F:1606
gw_utils
Definition gw_utils.F:15
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:3439
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:3503
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:282
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:1000
input_constants::rtp_method_bse
integer, parameter, public rtp_method_bse
Definition input_constants.F:1000
input_constants::g0w0
integer, parameter, public g0w0
Definition input_constants.F:1138
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
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kpoint_coulomb_2c
Routines to compute the Coulomb integral V_(alpha beta)(k) for a k-point k using lattice summation in...
Definition kpoint_coulomb_2c.F:15
kpoint_coulomb_2c::build_2c_coulomb_matrix_kp
subroutine, public build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, atomic_kind_set, size_lattice_sum, operator_type, ikp_start, ikp_end)
...
Definition kpoint_coulomb_2c.F:76
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
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:143
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
message_passing::mp_file_delete
subroutine, public mp_file_delete(filepath, info)
Deletes a file. Auxiliary routine to emulate 'replace' action for mp_file_open. Only the master proce...
Definition message_passing.F:3953
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_kpoints, fm_matrix_l, dimen_ri, ri_metric, do_kpoints, kpoints, put_mat_ks_env, regularization_ri, ikp_ext, do_build_cell_index)
...
Definition mp2_ri_2c.F:564
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
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::cfm_ikp_from_fm_gamma
subroutine, public cfm_ikp_from_fm_gamma(cfm_ikp, fm_gamma, ikp, qs_env, kpoints, basis_type)
...
Definition post_scf_bandstructure_utils.F:2535
post_scf_bandstructure_utils::get_all_vbm_cbm_bandgaps
subroutine, public get_all_vbm_cbm_bandgaps(bs_env)
...
Definition post_scf_bandstructure_utils.F:3039
post_scf_bandstructure_utils::mic_contribution_from_ikp
subroutine, public mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, cfm_w_ikp_freq_j, ikp, kpoints, basis_type, wkp_ext)
...
Definition post_scf_bandstructure_utils.F:2664
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:539
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors.F:11
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
rpa_gw_kpoints_util
Routines treating GW and RPA calculations with kpoints.
Definition rpa_gw_kpoints_util.F:13
rpa_gw_kpoints_util::cp_cfm_power
subroutine, public cp_cfm_power(matrix, threshold, exponent, min_eigval)
...
Definition rpa_gw_kpoints_util.F:274
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
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_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
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:187
particle_types::particle_type
Definition particle_types.F:35
post_scf_bandstructure_types::post_scf_bandstructure_type
Definition post_scf_bandstructure_types.F:154
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190