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