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