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