(git:f2099e5)
Loading...
Searching...
No Matches
gw_tensor_small_cell_full_kp.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief
10!> \author Jan Wilhelm
11!> \date 05.2024
12! **************************************************************************************************
13MODULE gw_tensor_small_cell_full_kp
14 USE bibliography, ONLY: pasquier2025,&
15 cite_reference
16 USE cp_cfm_basic_linalg, ONLY: cp_cfm_get_diag
17 USE cp_cfm_types, ONLY: cp_cfm_create,&
18 cp_cfm_get_info,&
19 cp_cfm_release,&
20 cp_cfm_to_cfm,&
21 cp_cfm_type
22 USE cp_fm_types, ONLY: cp_fm_set_all,&
23 cp_fm_type
24 USE dbt_api, ONLY: dbt_clear,&
25 dbt_contract,&
26 dbt_copy,&
27 dbt_create,&
28 dbt_destroy,&
29 dbt_type
30 USE gw_utils, ONLY: add_r,&
31 analyt_conti_and_print,&
32 de_init_bs_env,&
33 get_v_tr_r,&
34 is_cell_in_index_to_cell,&
35 time_to_freq
36 USE gw_utils_communication, ONLY: fm_to_local_array,&
37 fm_to_local_tensor,&
38 local_array_to_fm,&
39 local_dbt_to_global_fm
40 USE gw_utils_fm, ONLY: cfm_contract_aba,&
41 local_add_on_diag,&
42 local_complex_power
43 USE kinds, ONLY: dp,&
44 int_8
45 USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp_small_cell
46 USE kpoint_k_r_trafo_simple, ONLY: add_kp_to_all_rs,&
47 fm_add_kp_to_all_rs,&
48 fm_rs_to_kp,&
49 rs_to_kp
50 USE machine, ONLY: m_walltime
51 USE mathconstants, ONLY: z_one,&
52 z_zero
53 USE mathlib, ONLY: gemm_square
54 USE parallel_gemm_api, ONLY: parallel_gemm
55 USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
56 USE post_scf_bandstructure_utils, ONLY: get_all_vbm_cbm_bandgaps
57 USE qs_environment_types, ONLY: qs_environment_type
58#include "./base/base_uses.f90"
59
60 IMPLICIT NONE
61
62 PRIVATE
63
64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_tensor_small_cell_full_kp'
65
66 PUBLIC :: gw_calc_tensor_small_cell_full_kp
67
68CONTAINS
69
70! **************************************************************************************************
71!> \brief Perform GW band structure calculation
72!> \param qs_env ...
73!> \param bs_env Band-structure environment containing GW parameters.
74!> \par History
75!> * 05.2024 created [Jan Wilhelm]
76! **************************************************************************************************
77 SUBROUTINE gw_calc_tensor_small_cell_full_kp(qs_env, bs_env)
78 TYPE(qs_environment_type), POINTER :: qs_env
79 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
80
81 CHARACTER(LEN=*), PARAMETER :: routinen = 'gw_calc_tensor_small_cell_full_kp'
82
83 INTEGER :: handle
84
85 CALL timeset(routinen, handle)
86
87 CALL cite_reference(pasquier2025)
88
89 ! G^occ_µλ(i|τ|,k) = sum_n^occ C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k)
90 ! G^vir_µλ(i|τ|,k) = sum_n^vir C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k)
91 ! k-point k -> cell S: G^occ/vir_µλ^S(i|τ|) = sum_k w_k G^occ/vir_µλ(i|τ|,k) e^(ikS)
92 ! χ_PQ^R(iτ) = sum_λR1νR2 [ sum_µS (µR1-S νR2 | P0) G^vir_λµ^S(i|τ|) ]
93 ! [ sum_σS (σR2-S λR1 | QR) G^occ_νσ^S(i|τ|) ]
94 CALL compute_chi(bs_env)
95
96 ! χ_PQ^R(iτ) -> χ_PQ(iω,k) -> ε_PQ(iω,k) -> W_PQ(iω,k) -> Ŵ(iω,k) = M^-1(k)*W(iω,k)*M^-1(k)
97 ! -> Ŵ_PQ^R(iτ)
98 CALL compute_w_real_space(bs_env, qs_env)
99
100 ! D_µν(k) = sum_n^occ C^*_µn(k) C_νn(k), V^tr_PQ^R = <phi_P,0|V^tr|phi_Q,R>
101 ! V^tr(k) = sum_R e^ikR V^tr^R, M(k) = sum_R e^ikR M^R, M(k) -> M^-1(k)
102 ! -> Ṽ^tr(k) = M^-1(k) * V^tr(k) * M^-1(k) -> Ṽ^tr_PQ^R = sum_k w_k e^-ikR Ṽ^tr_PQ(k)
103 ! Σ^x_λσ^R = sum_PR1νS1 [ sum_µS2 (λ0 µS1-S2 | PR1 ) D_µν^S2 ]
104 ! [ sum_QR2 (σR νS1 | QR1-R2) Ṽ^tr_PQ^R2 ]
105 CALL compute_sigma_x(bs_env, qs_env)
106
107 ! Σ^c_λσ^R(iτ) = sum_PR1νS1 [ sum_µS2 (λ0 µS1-S2 | PR1 ) G^occ/vir_µν^S2(i|τ|) ]
108 ! [ sum_QR2 (σR νS1 | QR1-R2) Ŵ_PQ^R2(iτ) ]
109 CALL compute_sigma_c(bs_env)
110
111 ! Σ^c_λσ^R(iτ,k=0) -> Σ^c_nn(ϵ,k); ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ,k) + Σ^x_nn(k) - v^xc_nn(k)
112 CALL compute_qp_energies(bs_env)
113
114 CALL de_init_bs_env(qs_env, bs_env)
115
116 CALL timestop(handle)
117
118 END SUBROUTINE gw_calc_tensor_small_cell_full_kp
119
120! **************************************************************************************************
121!> \brief ...
122!> \param bs_env ...
123! **************************************************************************************************
124 SUBROUTINE compute_chi(bs_env)
125 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
126
127 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_chi'
128
129 INTEGER :: cell_dr(3), cell_r1(3), cell_r2(3), &
130 handle, i_cell_delta_r, i_cell_r1, &
131 i_cell_r2, i_t, i_task_delta_r_local, &
132 ispin
133 LOGICAL :: cell_found
134 REAL(kind=dp) :: t1, tau
135 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: gocc_s, gvir_s, t_chi_r
136 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_gocc, t_gvir
137
138 CALL timeset(routinen, handle)
139
140 DO i_t = 1, bs_env%num_time_freq_points
141
142 CALL dbt_create_2c_r(gocc_s, bs_env%t_G, bs_env%nimages_scf_desymm)
143 CALL dbt_create_2c_r(gvir_s, bs_env%t_G, bs_env%nimages_scf_desymm)
144 CALL dbt_create_2c_r(t_chi_r, bs_env%t_chi, bs_env%nimages_scf_desymm)
145 CALL dbt_create_3c_r1_r2(t_gocc, bs_env%t_RI_AO__AO, bs_env%nimages_3c, bs_env%nimages_3c)
146 CALL dbt_create_3c_r1_r2(t_gvir, bs_env%t_RI_AO__AO, bs_env%nimages_3c, bs_env%nimages_3c)
147
148 t1 = m_walltime()
149 tau = bs_env%time_frequency_grid%imaginary_time(i_t)
150
151 DO ispin = 1, bs_env%n_spin
152
153 ! 1. compute G^occ,S(iτ) and G^vir^S(iτ) in imaginary time for cell S
154 ! Background: G^σ,S(iτ) = G^occ,S,σ(iτ) * Θ(-τ) + G^vir,S,σ(iτ) * Θ(τ), σ ∈ {↑,↓}
155 ! G^occ_µλ(i|τ|,k) = sum_n^occ C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k)
156 ! G^vir_µλ(i|τ|,k) = sum_n^vir C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k)
157 ! k-point k -> cell S: G^occ/vir_µλ^S(i|τ|) = sum_k w_k G^occ/vir_µλ(i|τ|,k) e^(ikS)
158 CALL g_occ_vir(bs_env, tau, gocc_s, ispin, occ=.true., vir=.false.)
159 CALL g_occ_vir(bs_env, tau, gvir_s, ispin, occ=.false., vir=.true.)
160
161 ! loop over ΔR = R_1 - R_2 which are local in the tensor subgroup
162 DO i_task_delta_r_local = 1, bs_env%n_tasks_Delta_R_local
163
164 IF (bs_env%skip_DR_chi(i_task_delta_r_local)) cycle
165
166 i_cell_delta_r = bs_env%task_Delta_R(i_task_delta_r_local)
167
168 DO i_cell_r2 = 1, bs_env%nimages_3c
169
170 cell_r2(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_r2)
171 cell_dr(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_delta_r)
172
173 ! R_1 = R_2 + ΔR (from ΔR = R_2 - R_1)
174 CALL add_r(cell_r2, cell_dr, bs_env%index_to_cell_3c, cell_r1, &
175 cell_found, bs_env%cell_to_index_3c, i_cell_r1)
176
177 ! 3-cells check because in M^vir_νR2,λR1,QR (step 3.): R2 is index on ν
178 IF (.NOT. cell_found) cycle
179 ! 2. M^occ/vir_λR1,νR2,P0 = sum_µS (λR1 µR2-S | P0) G^occ/vir_νµ^S(iτ)
180 CALL g_times_3c(gocc_s, t_gocc, bs_env, i_cell_r1, i_cell_r2, &
181 i_task_delta_r_local, bs_env%skip_DR_R12_S_Goccx3c_chi)
182 CALL g_times_3c(gvir_s, t_gvir, bs_env, i_cell_r2, i_cell_r1, &
183 i_task_delta_r_local, bs_env%skip_DR_R12_S_Gvirx3c_chi)
184
185 END DO ! i_cell_R2
186
187 ! 3. χ_PQ^R(iτ) = sum_λR1,νR2 M^occ_λR1,νR2,P0 M^vir_νR2,λR1,QR
188 CALL contract_m_occ_vir_to_chi(t_gocc, t_gvir, t_chi_r, bs_env, &
189 i_task_delta_r_local)
190
191 END DO ! i_cell_Delta_R_local
192
193 END DO ! ispin
194
195 CALL bs_env%para_env%sync()
196
197 CALL local_dbt_to_global_fm(t_chi_r, bs_env%fm_chi_R_t(:, i_t), bs_env%mat_RI_RI, &
198 bs_env%mat_RI_RI_tensor, bs_env)
199
200 CALL destroy_t_1d(gocc_s)
201 CALL destroy_t_1d(gvir_s)
202 CALL destroy_t_1d(t_chi_r)
203 CALL destroy_t_2d(t_gocc)
204 CALL destroy_t_2d(t_gvir)
205
206 IF (bs_env%unit_nr > 0) THEN
207 WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
208 'Computed χ^R(iτ) for time point', i_t, ' /', bs_env%num_time_freq_points, &
209 ', Execution time', m_walltime() - t1, ' s'
210 END IF
211
212 END DO ! i_t
213
214 CALL timestop(handle)
215
216 END SUBROUTINE compute_chi
217
218! *************************************************************************************************
219!> \brief ...
220!> \param R ...
221!> \param template ...
222!> \param nimages ...
223! **************************************************************************************************
224 SUBROUTINE dbt_create_2c_r(R, template, nimages)
225
226 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: r
227 TYPE(dbt_type) :: template
228 INTEGER :: nimages
229
230 CHARACTER(LEN=*), PARAMETER :: routinen = 'dbt_create_2c_R'
231
232 INTEGER :: handle, i_cell_s
233
234 CALL timeset(routinen, handle)
235
236 ALLOCATE (r(nimages))
237 DO i_cell_s = 1, nimages
238 CALL dbt_create(template, r(i_cell_s))
239 END DO
240
241 CALL timestop(handle)
242
243 END SUBROUTINE dbt_create_2c_r
244
245! **************************************************************************************************
246!> \brief ...
247!> \param t_3c_R1_R2 ...
248!> \param t_3c_template ...
249!> \param nimages_1 ...
250!> \param nimages_2 ...
251! **************************************************************************************************
252 SUBROUTINE dbt_create_3c_r1_r2(t_3c_R1_R2, t_3c_template, nimages_1, nimages_2)
253
254 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_r1_r2
255 TYPE(dbt_type) :: t_3c_template
256 INTEGER :: nimages_1, nimages_2
257
258 CHARACTER(LEN=*), PARAMETER :: routinen = 'dbt_create_3c_R1_R2'
259
260 INTEGER :: handle, i_cell, j_cell
261
262 CALL timeset(routinen, handle)
263
264 ALLOCATE (t_3c_r1_r2(nimages_1, nimages_2))
265 DO i_cell = 1, nimages_1
266 DO j_cell = 1, nimages_2
267 CALL dbt_create(t_3c_template, t_3c_r1_r2(i_cell, j_cell))
268 END DO
269 END DO
270
271 CALL timestop(handle)
272
273 END SUBROUTINE dbt_create_3c_r1_r2
274
275! **************************************************************************************************
276!> \brief ...
277!> \param t_G_S ...
278!> \param t_M ...
279!> \param bs_env ...
280!> \param i_cell_R1 ...
281!> \param i_cell_R2 ...
282!> \param i_task_Delta_R_local ...
283!> \param skip_DR_R1_S_Gx3c ...
284! **************************************************************************************************
285 SUBROUTINE g_times_3c(t_G_S, t_M, bs_env, i_cell_R1, i_cell_R2, i_task_Delta_R_local, &
286 skip_DR_R1_S_Gx3c)
287 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_g_s
288 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_m
289 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
290 INTEGER :: i_cell_r1, i_cell_r2, &
291 i_task_delta_r_local
292 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: skip_dr_r1_s_gx3c
293
294 CHARACTER(LEN=*), PARAMETER :: routinen = 'G_times_3c'
295
296 INTEGER :: handle, i_cell_r1_p_s, i_cell_s
297 INTEGER(KIND=int_8) :: flop
298 INTEGER, DIMENSION(3) :: cell_r1, cell_r1_plus_cell_s, cell_r2, &
299 cell_s
300 LOGICAL :: cell_found
301 TYPE(dbt_type) :: t_3c_int
302
303 CALL timeset(routinen, handle)
304
305 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_int)
306
307 cell_r1(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_r1)
308 cell_r2(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_r2)
309
310 DO i_cell_s = 1, bs_env%nimages_scf_desymm
311
312 IF (skip_dr_r1_s_gx3c(i_task_delta_r_local, i_cell_r1, i_cell_s)) cycle
313
314 cell_s(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_s)
315 cell_r1_plus_cell_s(1:3) = cell_r1(1:3) + cell_s(1:3)
316
317 CALL is_cell_in_index_to_cell(cell_r1_plus_cell_s, bs_env%index_to_cell_3c, cell_found)
318
319 IF (.NOT. cell_found) cycle
320
321 i_cell_r1_p_s = bs_env%cell_to_index_3c(cell_r1_plus_cell_s(1), cell_r1_plus_cell_s(2), &
322 cell_r1_plus_cell_s(3))
323
324 IF (bs_env%nblocks_3c(i_cell_r2, i_cell_r1_p_s) == 0) cycle
325
326 CALL get_t_3c_int(t_3c_int, bs_env, i_cell_r2, i_cell_r1_p_s)
327
328 CALL dbt_contract(alpha=1.0_dp, &
329 tensor_1=t_3c_int, &
330 tensor_2=t_g_s(i_cell_s), &
331 beta=1.0_dp, &
332 tensor_3=t_m(i_cell_r1, i_cell_r2), &
333 contract_1=[3], notcontract_1=[1, 2], map_1=[1, 2], &
334 contract_2=[2], notcontract_2=[1], map_2=[3], &
335 filter_eps=bs_env%eps_filter, flop=flop)
336
337 IF (flop == 0_int_8) skip_dr_r1_s_gx3c(i_task_delta_r_local, i_cell_r1, i_cell_s) = .true.
338
339 END DO
340
341 CALL dbt_destroy(t_3c_int)
342
343 CALL timestop(handle)
344
345 END SUBROUTINE g_times_3c
346
347! **************************************************************************************************
348!> \brief ...
349!> \param t_3c_int ...
350!> \param bs_env ...
351!> \param j_cell ...
352!> \param k_cell ...
353! **************************************************************************************************
354 SUBROUTINE get_t_3c_int(t_3c_int, bs_env, j_cell, k_cell)
355
356 TYPE(dbt_type) :: t_3c_int
357 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
358 INTEGER :: j_cell, k_cell
359
360 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_t_3c_int'
361
362 INTEGER :: handle
363
364 CALL timeset(routinen, handle)
365
366 CALL dbt_clear(t_3c_int)
367 IF (j_cell < k_cell) THEN
368 CALL dbt_copy(bs_env%t_3c_int(k_cell, j_cell), t_3c_int, order=[1, 3, 2])
369 ELSE
370 CALL dbt_copy(bs_env%t_3c_int(j_cell, k_cell), t_3c_int)
371 END IF
372
373 CALL timestop(handle)
374
375 END SUBROUTINE get_t_3c_int
376
377! **************************************************************************************************
378!> \brief ...
379!> \param bs_env ...
380!> \param tau ...
381!> \param G_S ...
382!> \param ispin ...
383!> \param occ ...
384!> \param vir ...
385! **************************************************************************************************
386 SUBROUTINE g_occ_vir(bs_env, tau, G_S, ispin, occ, vir)
387 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
388 REAL(kind=dp) :: tau
389 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: g_s
390 INTEGER :: ispin
391 LOGICAL :: occ, vir
392
393 CHARACTER(LEN=*), PARAMETER :: routinen = 'G_occ_vir'
394
395 INTEGER :: handle, homo, i_cell_s, ikp, j, &
396 j_col_local, n_mo, ncol_local, &
397 nimages, nkp
398 INTEGER, DIMENSION(:), POINTER :: col_indices
399 REAL(kind=dp) :: tau_e
400
401 CALL timeset(routinen, handle)
402
403 cpassert(occ .NEQV. vir)
404
405 CALL cp_cfm_get_info(matrix=bs_env%cfm_work_mo, &
406 ncol_local=ncol_local, &
407 col_indices=col_indices)
408
409 nkp = bs_env%nkp_scf_desymm
410 nimages = bs_env%nimages_scf_desymm
411 n_mo = bs_env%n_ao
412 homo = bs_env%n_occ(ispin)
413
414 DO i_cell_s = 1, bs_env%nimages_scf_desymm
415 CALL cp_fm_set_all(bs_env%fm_G_S(i_cell_s), 0.0_dp)
416 END DO
417
418 DO ikp = 1, nkp
419
420 ! get C_µn(k)
421 CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), bs_env%cfm_work_mo)
422
423 ! G^occ/vir_µλ(i|τ|,k) = sum_n^occ/vir C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k)
424 DO j_col_local = 1, ncol_local
425
426 j = col_indices(j_col_local)
427
428 ! 0.5 * |(ϵ_nk-ϵ_F)τ|
429 tau_e = abs(tau*0.5_dp*(bs_env%eigenval_scf(j, ikp, ispin) - bs_env%e_fermi(ispin)))
430
431 IF (tau_e < bs_env%stabilize_exp) THEN
432 bs_env%cfm_work_mo%local_data(:, j_col_local) = &
433 bs_env%cfm_work_mo%local_data(:, j_col_local)*exp(-tau_e)
434 ELSE
435 bs_env%cfm_work_mo%local_data(:, j_col_local) = z_zero
436 END IF
437
438 IF ((occ .AND. j > homo) .OR. (vir .AND. j <= homo)) THEN
439 bs_env%cfm_work_mo%local_data(:, j_col_local) = z_zero
440 END IF
441
442 END DO
443
444 CALL parallel_gemm(transa="N", transb="C", m=n_mo, n=n_mo, k=n_mo, alpha=z_one, &
445 matrix_a=bs_env%cfm_work_mo, matrix_b=bs_env%cfm_work_mo, &
446 beta=z_zero, matrix_c=bs_env%cfm_work_mo_2)
447
448 ! trafo k-point k -> cell S: G^occ/vir_µλ(i|τ|,k) -> G^occ/vir,S_µλ(i|τ|)
449 CALL fm_add_kp_to_all_rs(bs_env%cfm_work_mo_2, bs_env%fm_G_S, &
450 bs_env%kpoints_scf_desymm, ikp)
451
452 END DO ! ikp
453
454 ! replicate to tensor from local tensor group
455 DO i_cell_s = 1, bs_env%nimages_scf_desymm
456 CALL fm_to_local_tensor(bs_env%fm_G_S(i_cell_s), bs_env%mat_ao_ao%matrix, &
457 bs_env%mat_ao_ao_tensor%matrix, g_s(i_cell_s), bs_env)
458 END DO
459
460 CALL timestop(handle)
461
462 END SUBROUTINE g_occ_vir
463
464! **************************************************************************************************
465!> \brief ...
466!> \param t_Gocc ...
467!> \param t_Gvir ...
468!> \param t_chi_R ...
469!> \param bs_env ...
470!> \param i_task_Delta_R_local ...
471! **************************************************************************************************
472 SUBROUTINE contract_m_occ_vir_to_chi(t_Gocc, t_Gvir, t_chi_R, bs_env, i_task_Delta_R_local)
473 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_gocc, t_gvir
474 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_chi_r
475 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
476 INTEGER :: i_task_delta_r_local
477
478 CHARACTER(LEN=*), PARAMETER :: routinen = 'contract_M_occ_vir_to_chi'
479
480 INTEGER :: handle, i_cell_delta_r, i_cell_r, &
481 i_cell_r1, i_cell_r1_minus_r, &
482 i_cell_r2, i_cell_r2_minus_r
483 INTEGER(KIND=int_8) :: flop, flop_tmp
484 INTEGER, DIMENSION(3) :: cell_dr, cell_r, cell_r1, &
485 cell_r1_minus_r, cell_r2, &
486 cell_r2_minus_r
487 LOGICAL :: cell_found
488 TYPE(dbt_type) :: t_gocc_2, t_gvir_2
489
490 CALL timeset(routinen, handle)
491
492 CALL dbt_create(bs_env%t_RI__AO_AO, t_gocc_2)
493 CALL dbt_create(bs_env%t_RI__AO_AO, t_gvir_2)
494
495 flop = 0_int_8
496
497 ! χ_PQ^R(iτ) = sum_λR1,νR2 M^occ_λR1,νR2,P0 M^vir_νR2,λR1,QR
498 DO i_cell_r = 1, bs_env%nimages_scf_desymm
499
500 DO i_cell_r2 = 1, bs_env%nimages_3c
501
502 IF (bs_env%skip_DR_R_R2_MxM_chi(i_task_delta_r_local, i_cell_r2, i_cell_r)) cycle
503
504 i_cell_delta_r = bs_env%task_Delta_R(i_task_delta_r_local)
505
506 cell_r(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_r)
507 cell_r2(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_r2)
508 cell_dr(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_delta_r)
509
510 ! R_1 = R_2 + ΔR (from ΔR = R_2 - R_1)
511 CALL add_r(cell_r2, cell_dr, bs_env%index_to_cell_3c, cell_r1, &
512 cell_found, bs_env%cell_to_index_3c, i_cell_r1)
513 IF (.NOT. cell_found) cycle
514
515 ! R_1 - R
516 CALL add_r(cell_r1, -cell_r, bs_env%index_to_cell_3c, cell_r1_minus_r, &
517 cell_found, bs_env%cell_to_index_3c, i_cell_r1_minus_r)
518 IF (.NOT. cell_found) cycle
519
520 ! R_2 - R
521 CALL add_r(cell_r2, -cell_r, bs_env%index_to_cell_3c, cell_r2_minus_r, &
522 cell_found, bs_env%cell_to_index_3c, i_cell_r2_minus_r)
523 IF (.NOT. cell_found) cycle
524
525 ! reorder tensors for efficient contraction to χ_PQ^R
526 CALL dbt_copy(t_gocc(i_cell_r1, i_cell_r2), t_gocc_2, order=[1, 3, 2])
527 CALL dbt_copy(t_gvir(i_cell_r2_minus_r, i_cell_r1_minus_r), t_gvir_2)
528
529 ! χ_PQ^R(iτ) = sum_λR1,νR2 M^occ_λR1,νR2,P0 M^vir_νR2,λR1,QR
530 CALL dbt_contract(alpha=bs_env%spin_degeneracy, &
531 tensor_1=t_gocc_2, tensor_2=t_gvir_2, &
532 beta=1.0_dp, tensor_3=t_chi_r(i_cell_r), &
533 contract_1=[2, 3], notcontract_1=[1], map_1=[1], &
534 contract_2=[2, 3], notcontract_2=[1], map_2=[2], &
535 filter_eps=bs_env%eps_filter, move_data=.true., flop=flop_tmp)
536
537 IF (flop_tmp == 0_int_8) bs_env%skip_DR_R_R2_MxM_chi(i_task_delta_r_local, &
538 i_cell_r2, i_cell_r) = .true.
539
540 flop = flop + flop_tmp
541
542 END DO ! i_cell_R2
543
544 END DO ! i_cell_R
545
546 IF (flop == 0_int_8) bs_env%skip_DR_chi(i_task_delta_r_local) = .true.
547
548 ! remove all data from t_Gocc and t_Gvir to safe memory
549 DO i_cell_r1 = 1, bs_env%nimages_3c
550 DO i_cell_r2 = 1, bs_env%nimages_3c
551 CALL dbt_clear(t_gocc(i_cell_r1, i_cell_r2))
552 CALL dbt_clear(t_gvir(i_cell_r1, i_cell_r2))
553 END DO
554 END DO
555
556 CALL dbt_destroy(t_gocc_2)
557 CALL dbt_destroy(t_gvir_2)
558
559 CALL timestop(handle)
560
561 END SUBROUTINE contract_m_occ_vir_to_chi
562
563! **************************************************************************************************
564!> \brief ...
565!> \param bs_env ...
566!> \param qs_env ...
567! **************************************************************************************************
568 SUBROUTINE compute_w_real_space(bs_env, qs_env)
569 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
570 TYPE(qs_environment_type), POINTER :: qs_env
571
572 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_W_real_space'
573
574 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chi_k_w, eps_k_w, w_k_w
575 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: m_inv, m_inv_v_sqrt, v_sqrt
576 INTEGER :: handle, i_t, ikp, ikp_local, j_w, n_ri, &
577 nimages_scf_desymm
578 REAL(kind=dp) :: freq_j, t1, time_i, weight_ij
579 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: chi_r, mwm_r, w_r
580
581 CALL timeset(routinen, handle)
582
583 n_ri = bs_env%n_RI
584 nimages_scf_desymm = bs_env%nimages_scf_desymm
585
586 ALLOCATE (chi_k_w(n_ri, n_ri), eps_k_w(n_ri, n_ri), w_k_w(n_ri, n_ri))
587 ALLOCATE (chi_r(n_ri, n_ri, nimages_scf_desymm), w_r(n_ri, n_ri, nimages_scf_desymm), &
588 mwm_r(n_ri, n_ri, nimages_scf_desymm))
589
590 t1 = m_walltime()
591
592 CALL compute_minv_and_vsqrt(bs_env, qs_env, m_inv_v_sqrt, m_inv, v_sqrt)
593
594 IF (bs_env%unit_nr > 0) THEN
595 WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
596 'Computed V_PQ(k),', 'Execution time', m_walltime() - t1, ' s'
597 WRITE (bs_env%unit_nr, '(A)') ' '
598 END IF
599
600 t1 = m_walltime()
601
602 DO j_w = 1, bs_env%num_time_freq_points
603
604 ! χ_PQ^R(iτ) -> χ_PQ^R(iω_j) (which is stored in chi_R, single ω_j from j_w loop)
605 chi_r(:, :, :) = 0.0_dp
606 DO i_t = 1, bs_env%num_time_freq_points
607 freq_j = bs_env%time_frequency_grid%frequency(j_w)
608 time_i = bs_env%time_frequency_grid%imaginary_time(i_t)
609 weight_ij = bs_env%time_frequency_grid%cosine_time_to_frequency_weights(j_w, i_t)*cos(time_i*freq_j)
610
611 CALL fm_to_local_array(bs_env%fm_chi_R_t(:, i_t), chi_r, weight_ij, add=.true.)
612 END DO
613
614 ikp_local = 0
615 w_r(:, :, :) = 0.0_dp
616 DO ikp = 1, bs_env%nkp_chi_eps_W_orig_plus_extra
617
618 ! trivial parallelization over k-points
619 IF (modulo(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
620
621 ikp_local = ikp_local + 1
622
623 ! 1. χ_PQ^R(iω_j) -> χ_PQ(iω_j,k)
624 CALL rs_to_kp(chi_r, chi_k_w, bs_env%kpoints_scf_desymm%index_to_cell, &
625 bs_env%kpoints_chi_eps_W%xkp(1:3, ikp))
626
627 ! 2. remove negative eigenvalues from χ_PQ(iω,k)
628 CALL local_complex_power(chi_k_w, 1.0_dp, bs_env%eps_eigval_mat_RI)
629
630 ! 3. ε(iω_j,k_i) = I + V^0.5(k_i)*M^-1(k_i)*χ(iω_j,k_i)*M^-1(k_i)*V^0.5(k_i)
631 ! a) eps_work = V^0.5(k_i)*M^-1(k_i)*χ(iω_j,k_i)*M^-1(k_i)*V^0.5(k_i)
632 CALL gemm_square(m_inv_v_sqrt(:, :, ikp_local), 'C', chi_k_w, 'N', &
633 m_inv_v_sqrt(:, :, ikp_local), 'N', eps_k_w)
634
635 ! b) ε(iω_j,k_i) = eps_work + I
636 CALL local_add_on_diag(eps_k_w, z_one)
637
638 ! 4. W(iω_j,k_i) = V^0.5(k_i)*(ε^-1(iω_j,k_i)-I)*V^0.5(k_i)
639 ! a) Invert ε(iω_j,k_i) using its Cholesky decomposition
640 CALL local_complex_power(eps_k_w, -1.0_dp, 0.0_dp)
641
642 ! b) ε^-1(iω_j,k_i)-I
643 CALL local_add_on_diag(eps_k_w, -z_one)
644
645 ! c) W(iω_j,k_i) = V^0.5(k_i)*(ε^-1(iω_j,k_i)-I)*V^0.5(k_i)
646 CALL gemm_square(v_sqrt(:, :, ikp_local), 'N', eps_k_w, 'N', &
647 v_sqrt(:, :, ikp_local), 'C', w_k_w)
648
649 ! 5. W(iω,k_i) -> W^R(iω) = sum_k w_k e^(-ikR) W(iω,k) (k-point extrapolation here)
650 CALL add_kp_to_all_rs(w_k_w, w_r, bs_env%kpoints_chi_eps_W, ikp, &
651 index_to_cell_ext=bs_env%kpoints_scf_desymm%index_to_cell)
652
653 END DO ! ikp
654
655 CALL bs_env%para_env%sync()
656 CALL bs_env%para_env%sum(w_r)
657
658 ! 6. W^R(iω) -> W(iω,k) [k-mesh is not extrapolated for stable mult. with M^-1(k) ]
659 ! -> M^-1(k)*W(iω,k)*M^-1(k) =: Ŵ(iω,k) -> Ŵ^R(iω) (stored in MWM_R)
660 CALL mult_w_with_minv(w_r, mwm_r, bs_env, qs_env)
661
662 ! 7. Ŵ^R(iω) -> Ŵ^R(iτ) and to fully distributed fm matrix bs_env%fm_MWM_R_t
663 DO i_t = 1, bs_env%num_time_freq_points
664 freq_j = bs_env%time_frequency_grid%frequency(j_w)
665 time_i = bs_env%time_frequency_grid%imaginary_time(i_t)
666 weight_ij = bs_env%time_frequency_grid%cosine_frequency_to_time_weights(i_t, j_w)*cos(time_i*freq_j)
667 CALL local_array_to_fm(mwm_r, bs_env%fm_MWM_R_t(:, i_t), weight_ij, add=.true.)
668 END DO ! i_t
669
670 END DO ! j_w
671
672 IF (bs_env%unit_nr > 0) THEN
673 WRITE (bs_env%unit_nr, '(T2,A,T60,A,F7.1,A)') &
674 'Computed W_PQ(k,iω) for all k and τ,', 'Execution time', m_walltime() - t1, ' s'
675 WRITE (bs_env%unit_nr, '(A)') ' '
676 END IF
677
678 CALL timestop(handle)
679
680 END SUBROUTINE compute_w_real_space
681
682! **************************************************************************************************
683!> \brief ...
684!> \param bs_env ...
685!> \param qs_env ...
686!> \param M_inv_V_sqrt ...
687!> \param M_inv ...
688!> \param V_sqrt ...
689! **************************************************************************************************
690 SUBROUTINE compute_minv_and_vsqrt(bs_env, qs_env, M_inv_V_sqrt, M_inv, V_sqrt)
691 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
692 TYPE(qs_environment_type), POINTER :: qs_env
693 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: m_inv_v_sqrt, m_inv, v_sqrt
694
695 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_Minv_and_Vsqrt'
696
697 INTEGER :: handle, ikp, ikp_local, n_ri, nkp, &
698 nkp_local, nkp_orig
699 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: m_r
700
701 CALL timeset(routinen, handle)
702
703 nkp = bs_env%nkp_chi_eps_W_orig_plus_extra
704 nkp_orig = bs_env%nkp_chi_eps_W_orig
705 n_ri = bs_env%n_RI
706
707 nkp_local = 0
708 DO ikp = 1, nkp
709 ! trivial parallelization over k-points
710 IF (modulo(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
711 nkp_local = nkp_local + 1
712 END DO
713
714 ALLOCATE (m_inv_v_sqrt(n_ri, n_ri, nkp_local), m_inv(n_ri, n_ri, nkp_local), &
715 v_sqrt(n_ri, n_ri, nkp_local))
716
717 m_inv_v_sqrt(:, :, :) = z_zero
718 m_inv(:, :, :) = z_zero
719 v_sqrt(:, :, :) = z_zero
720
721 ! 1. 2c Coulomb integrals for the first "original" k-point grid
722 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
723 CALL build_2c_coulomb_matrix_kp_small_cell(v_sqrt, qs_env, bs_env%kpoints_chi_eps_W, &
724 bs_env%size_lattice_sum_V, basis_type="RI_AUX", &
725 ikp_start=1, ikp_end=nkp_orig)
726
727 ! 2. 2c Coulomb integrals for the second "extrapolation" k-point grid
728 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_extra
729 CALL build_2c_coulomb_matrix_kp_small_cell(v_sqrt, qs_env, bs_env%kpoints_chi_eps_W, &
730 bs_env%size_lattice_sum_V, basis_type="RI_AUX", &
731 ikp_start=nkp_orig + 1, ikp_end=nkp)
732
733 ! now get M^-1(k) and M^-1(k)*V^0.5(k)
734
735 ! compute M^R_PQ = <phi_P,0|V^tr(rc=3Å)|phi_Q,R> for RI metric
736 CALL get_v_tr_r(m_r, bs_env%ri_metric, bs_env%regularization_RI, bs_env, qs_env)
737
738 ikp_local = 0
739 DO ikp = 1, nkp
740
741 ! trivial parallelization
742 IF (modulo(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
743
744 ikp_local = ikp_local + 1
745
746 ! M(k) = sum_R e^ikR M^R
747 CALL rs_to_kp(m_r, m_inv(:, :, ikp_local), &
748 bs_env%kpoints_scf_desymm%index_to_cell, &
749 bs_env%kpoints_chi_eps_W%xkp(1:3, ikp))
750
751 ! invert M_PQ(k)
752 CALL local_complex_power(m_inv(:, :, ikp_local), -1.0_dp, 0.0_dp)
753
754 ! V^0.5(k)
755 CALL local_complex_power(v_sqrt(:, :, ikp_local), 0.5_dp, 0.0_dp)
756
757 ! M^-1(k)*V^0.5(k)
758 CALL gemm_square(m_inv(:, :, ikp_local), 'N', v_sqrt(:, :, ikp_local), 'C', m_inv_v_sqrt(:, :, ikp_local))
759
760 END DO ! ikp
761
762 CALL timestop(handle)
763
764 END SUBROUTINE compute_minv_and_vsqrt
765
766! **************************************************************************************************
767!> \brief ...
768!> \param W_R ...
769!> \param MWM_R ...
770!> \param bs_env ...
771!> \param qs_env ...
772! **************************************************************************************************
773 SUBROUTINE mult_w_with_minv(W_R, MWM_R, bs_env, qs_env)
774 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: w_r, mwm_r
775 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
776 TYPE(qs_environment_type), POINTER :: qs_env
777
778 CHARACTER(LEN=*), PARAMETER :: routinen = 'mult_W_with_Minv'
779
780 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: m_inv, w_k, work
781 INTEGER :: handle, ikp, n_ri
782 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: m_r
783
784 CALL timeset(routinen, handle)
785
786 ! compute M^R again
787 CALL get_v_tr_r(m_r, bs_env%ri_metric, bs_env%regularization_RI, bs_env, qs_env)
788
789 n_ri = bs_env%n_RI
790 ALLOCATE (m_inv(n_ri, n_ri), w_k(n_ri, n_ri), work(n_ri, n_ri))
791 mwm_r(:, :, :) = 0.0_dp
792
793 DO ikp = 1, bs_env%nkp_scf_desymm
794
795 ! trivial parallelization
796 IF (modulo(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
797
798 ! M(k) = sum_R e^ikR M^R
799 CALL rs_to_kp(m_r, m_inv, &
800 bs_env%kpoints_scf_desymm%index_to_cell, &
801 bs_env%kpoints_scf_desymm%xkp(1:3, ikp))
802
803 ! invert M_PQ(k)
804 CALL local_complex_power(m_inv, -1.0_dp, 0.0_dp)
805
806 ! W(k) = sum_R e^ikR W^R [only R in the supercell that is determined by the SCF k-mesh]
807 CALL rs_to_kp(w_r, w_k, &
808 bs_env%kpoints_scf_desymm%index_to_cell, &
809 bs_env%kpoints_scf_desymm%xkp(1:3, ikp))
810
811 ! Ŵ(k) = M^-1(k)*W^trunc(k)*M^-1(k)
812 CALL gemm_square(m_inv, 'N', w_k, 'N', m_inv, 'N', work)
813 w_k(:, :) = work(:, :)
814
815 ! Ŵ^R = sum_k w_k e^(-ikR) Ŵ^(k)
816 CALL add_kp_to_all_rs(w_k, mwm_r, bs_env%kpoints_scf_desymm, ikp)
817
818 END DO ! ikp
819
820 CALL bs_env%para_env%sync()
821 CALL bs_env%para_env%sum(mwm_r)
822
823 CALL timestop(handle)
824
825 END SUBROUTINE mult_w_with_minv
826
827! **************************************************************************************************
828!> \brief ...
829!> \param bs_env ...
830!> \param qs_env ...
831! **************************************************************************************************
832 SUBROUTINE compute_sigma_x(bs_env, qs_env)
833 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
834 TYPE(qs_environment_type), POINTER :: qs_env
835
836 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_Sigma_x'
837
838 INTEGER :: handle, i_task_delta_r_local, ispin
839 REAL(kind=dp) :: t1
840 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: d_s, mi_vtr_mi_r, sigma_x_r
841 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_v
842
843 CALL timeset(routinen, handle)
844
845 CALL dbt_create_2c_r(mi_vtr_mi_r, bs_env%t_W, bs_env%nimages_scf_desymm)
846 CALL dbt_create_2c_r(d_s, bs_env%t_G, bs_env%nimages_scf_desymm)
847 CALL dbt_create_2c_r(sigma_x_r, bs_env%t_G, bs_env%nimages_scf_desymm)
848 CALL dbt_create_3c_r1_r2(t_v, bs_env%t_RI_AO__AO, bs_env%nimages_3c, bs_env%nimages_3c)
849
850 t1 = m_walltime()
851
852 ! V^tr_PQ^R = <phi_P,0|V^tr|phi_Q,R>, V^tr(k) = sum_R e^ikR V^tr^R
853 ! M(k) = sum_R e^ikR M^R, M(k) -> M^-1(k) -> Ṽ^tr(k) = M^-1(k) * V^tr(k) * M^-1(k)
854 ! -> Ṽ^tr_PQ^R = sum_k w_k e^-ikR Ṽ^tr_PQ(k)
855 CALL get_minv_vtr_minv_r(mi_vtr_mi_r, bs_env, qs_env)
856
857 ! Σ^x_λσ^R = sum_PR1νS1 [ sum_µS2 (λ0 µS1-S2 | PR1 ) D_µν^S2 ]
858 ! [ sum_QR2 (σR νS1 | QR1-R2) Ṽ^tr_PQ^R2 ]
859 DO ispin = 1, bs_env%n_spin
860
861 ! compute D^S(iτ) for cell S from D_µν(k) = sum_n^occ C^*_µn(k) C_νn(k):
862 ! trafo k-point k -> cell S: D_µν^S = sum_k w_k D_µν(k) e^(ikS)
863 CALL g_occ_vir(bs_env, 0.0_dp, d_s, ispin, occ=.true., vir=.false.)
864
865 ! loop over ΔR = S_1 - R_1 which are local in the tensor subgroup
866 DO i_task_delta_r_local = 1, bs_env%n_tasks_Delta_R_local
867
868 ! M^V_σ0,νS1,PR1 = sum_QR2 ( σ0 νS1 | QR1-R2 ) Ṽ^tr_QP^R2 for i_task_local
869 CALL contract_w(t_v, mi_vtr_mi_r, bs_env, i_task_delta_r_local)
870
871 ! M^D_λ0,νS1,PR1 = sum_µS2 (λ0 µS1-S2 | PR1) D_µν^S2
872 ! Σ^x_λσ^R = sum_PR1νS1 M^D_λ0,νS1,PR1 * M^V_σR,νS1,PR1 for i_task_local, where
873 ! M^V_σR,νS1,PR1 = M^V_σ0,νS1-R,PR1-R
874 CALL contract_to_sigma(sigma_x_r, t_v, d_s, i_task_delta_r_local, bs_env, &
875 occ=.true., vir=.false., clear_t_w=.true., fill_skip=.false.)
876
877 END DO ! i_cell_Delta_R_local
878
879 CALL bs_env%para_env%sync()
880
881 CALL local_dbt_to_global_fm(sigma_x_r, bs_env%fm_Sigma_x_R, bs_env%mat_ao_ao, &
882 bs_env%mat_ao_ao_tensor, bs_env)
883
884 END DO ! ispin
885
886 IF (bs_env%unit_nr > 0) THEN
887 WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
888 'Computed Σ^x,', ' Execution time', m_walltime() - t1, ' s'
889 WRITE (bs_env%unit_nr, '(A)') ' '
890 END IF
891
892 CALL destroy_t_1d(mi_vtr_mi_r)
893 CALL destroy_t_1d(d_s)
894 CALL destroy_t_1d(sigma_x_r)
895 CALL destroy_t_2d(t_v)
896
897 CALL timestop(handle)
898
899 END SUBROUTINE compute_sigma_x
900
901! **************************************************************************************************
902!> \brief ...
903!> \param bs_env ...
904! **************************************************************************************************
905 SUBROUTINE compute_sigma_c(bs_env)
906 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
907
908 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_Sigma_c'
909
910 INTEGER :: handle, i_t, i_task_delta_r_local, ispin
911 REAL(kind=dp) :: t1, tau
912 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: gocc_s, gvir_s, sigma_c_r_neg_tau, &
913 sigma_c_r_pos_tau, w_r
914 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_w
915
916 CALL timeset(routinen, handle)
917
918 CALL dbt_create_2c_r(gocc_s, bs_env%t_G, bs_env%nimages_scf_desymm)
919 CALL dbt_create_2c_r(gvir_s, bs_env%t_G, bs_env%nimages_scf_desymm)
920 CALL dbt_create_2c_r(w_r, bs_env%t_W, bs_env%nimages_scf_desymm)
921 CALL dbt_create_3c_r1_r2(t_w, bs_env%t_RI_AO__AO, bs_env%nimages_3c, bs_env%nimages_3c)
922 CALL dbt_create_2c_r(sigma_c_r_neg_tau, bs_env%t_G, bs_env%nimages_scf_desymm)
923 CALL dbt_create_2c_r(sigma_c_r_pos_tau, bs_env%t_G, bs_env%nimages_scf_desymm)
924
925 ! Σ^c_λσ^R(iτ) = sum_PR1νS1 [ sum_µS2 (λ0 µS1-S2 | PR1 ) G^occ/vir_µν^S2(i|τ|) ]
926 ! [ sum_QR2 (σR νS1 | QR1-R2) Ŵ_PQ^R2(iτ) ]
927 DO i_t = 1, bs_env%num_time_freq_points
928
929 DO ispin = 1, bs_env%n_spin
930
931 t1 = m_walltime()
932
933 tau = bs_env%time_frequency_grid%imaginary_time(i_t)
934
935 ! G^occ_µλ(i|τ|,k) = sum_n^occ C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k), τ < 0
936 ! G^vir_µλ(i|τ|,k) = sum_n^vir C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k), τ > 0
937 ! k-point k -> cell S: G^occ/vir_µλ^S(i|τ|) = sum_k w_k G^occ/vir_µλ(i|τ|,k) e^(ikS)
938 CALL g_occ_vir(bs_env, tau, gocc_s, ispin, occ=.true., vir=.false.)
939 CALL g_occ_vir(bs_env, tau, gvir_s, ispin, occ=.false., vir=.true.)
940
941 ! write data of W^R_PQ(iτ) to W_R 2-index tensor
942 CALL fm_mwm_r_t_to_local_tensor_w_r(bs_env%fm_MWM_R_t(:, i_t), w_r, bs_env)
943
944 ! loop over ΔR = S_1 - R_1 which are local in the tensor subgroup
945 DO i_task_delta_r_local = 1, bs_env%n_tasks_Delta_R_local
946
947 IF (bs_env%skip_DR_Sigma(i_task_delta_r_local)) cycle
948
949 ! for i_task_local (i.e. fixed ΔR = S_1 - R_1) and for all τ (W(iτ) = W(-iτ)):
950 ! M^W_σ0,νS1,PR1 = sum_QR2 ( σ0 νS1 | QR1-R2 ) W(iτ)_QP^R2
951 CALL contract_w(t_w, w_r, bs_env, i_task_delta_r_local)
952
953 ! for τ < 0 and for i_task_local (i.e. fixed ΔR = S_1 - R_1):
954 ! M^G_λ0,νS1,PR1 = sum_µS2 (λ0 µS1-S2 | PR1) G^occ(i|τ|)_µν^S2
955 ! Σ^c_λσ^R(iτ) = sum_PR1νS1 M^G_λ0,νS1,PR1 * M^W_σR,νS1,PR1
956 ! where M^W_σR,νS1,PR1 = M^W_σ0,νS1-R,PR1-R
957 CALL contract_to_sigma(sigma_c_r_neg_tau, t_w, gocc_s, i_task_delta_r_local, bs_env, &
958 occ=.true., vir=.false., clear_t_w=.false., fill_skip=.false.)
959
960 ! for τ > 0: same as for τ < 0, but G^occ -> G^vir
961 CALL contract_to_sigma(sigma_c_r_pos_tau, t_w, gvir_s, i_task_delta_r_local, bs_env, &
962 occ=.false., vir=.true., clear_t_w=.true., fill_skip=.true.)
963
964 END DO ! i_cell_Delta_R_local
965
966 CALL bs_env%para_env%sync()
967
968 CALL local_dbt_to_global_fm(sigma_c_r_pos_tau, &
969 bs_env%fm_Sigma_c_R_pos_tau(:, i_t, ispin), &
970 bs_env%mat_ao_ao, bs_env%mat_ao_ao_tensor, bs_env)
971
972 CALL local_dbt_to_global_fm(sigma_c_r_neg_tau, &
973 bs_env%fm_Sigma_c_R_neg_tau(:, i_t, ispin), &
974 bs_env%mat_ao_ao, bs_env%mat_ao_ao_tensor, bs_env)
975
976 IF (bs_env%unit_nr > 0) THEN
977 WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F7.1,A)') &
978 'Computed Σ^c(iτ) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
979 ', Execution time', m_walltime() - t1, ' s'
980 END IF
981
982 END DO ! ispin
983
984 END DO ! i_t
985
986 CALL destroy_t_1d(gocc_s)
987 CALL destroy_t_1d(gvir_s)
988 CALL destroy_t_1d(w_r)
989 CALL destroy_t_1d(sigma_c_r_neg_tau)
990 CALL destroy_t_1d(sigma_c_r_pos_tau)
991 CALL destroy_t_2d(t_w)
992
993 CALL timestop(handle)
994
995 END SUBROUTINE compute_sigma_c
996
997! **************************************************************************************************
998!> \brief ...
999!> \param Mi_Vtr_Mi_R ...
1000!> \param bs_env ...
1001!> \param qs_env ...
1002! **************************************************************************************************
1003 SUBROUTINE get_minv_vtr_minv_r(Mi_Vtr_Mi_R, bs_env, qs_env)
1004 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: mi_vtr_mi_r
1005 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1006 TYPE(qs_environment_type), POINTER :: qs_env
1007
1008 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_Minv_Vtr_Minv_R'
1009
1010 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: m_kp, mi_vtr_mi_kp, v_tr_kp
1011 INTEGER :: handle, i_cell_r, ikp, n_ri, &
1012 nimages_scf, nkp_scf
1013 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: m_r, mi_vtr_mi_r_arr, v_tr_r
1014
1015 CALL timeset(routinen, handle)
1016
1017 nimages_scf = bs_env%nimages_scf_desymm
1018 nkp_scf = bs_env%kpoints_scf_desymm%nkp
1019 n_ri = bs_env%n_RI
1020
1021 CALL get_v_tr_r(v_tr_r, bs_env%trunc_coulomb, 0.0_dp, bs_env, qs_env)
1022 CALL get_v_tr_r(m_r, bs_env%ri_metric, bs_env%regularization_RI, bs_env, qs_env)
1023
1024 ALLOCATE (v_tr_kp(n_ri, n_ri), m_kp(n_ri, n_ri), &
1025 mi_vtr_mi_kp(n_ri, n_ri), mi_vtr_mi_r_arr(n_ri, n_ri, nimages_scf))
1026 mi_vtr_mi_r_arr(:, :, :) = 0.0_dp
1027
1028 DO ikp = 1, nkp_scf
1029 ! trivial parallelization
1030 IF (modulo(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) cycle
1031 ! V_tr(k) = sum_R e^ikR V_tr^R
1032 CALL rs_to_kp(v_tr_r, v_tr_kp, bs_env%kpoints_scf_desymm%index_to_cell, &
1033 bs_env%kpoints_scf_desymm%xkp(1:3, ikp))
1034 ! M(k) = sum_R e^ikR M^R
1035 CALL rs_to_kp(m_r, m_kp, bs_env%kpoints_scf_desymm%index_to_cell, &
1036 bs_env%kpoints_scf_desymm%xkp(1:3, ikp))
1037 ! M(k) -> M^-1(k)
1038 CALL local_complex_power(m_kp, -1.0_dp, 0.0_dp)
1039 ! Ṽ(k) = M^-1(k) * V_tr(k) * M^-1(k)
1040 CALL gemm_square(m_kp, 'N', v_tr_kp, 'N', m_kp, 'N', mi_vtr_mi_kp)
1041 ! Ṽ^R = sum_k w_k e^-ikR Ṽ(k)
1042 CALL add_kp_to_all_rs(mi_vtr_mi_kp, mi_vtr_mi_r_arr, bs_env%kpoints_scf_desymm, ikp)
1043 END DO
1044 CALL bs_env%para_env%sync()
1045 CALL bs_env%para_env%sum(mi_vtr_mi_r_arr)
1046
1047 ! use bs_env%fm_chi_R_t for temporary storage
1048 CALL local_array_to_fm(mi_vtr_mi_r_arr, bs_env%fm_chi_R_t(:, 1))
1049
1050 ! communicate Mi_Vtr_Mi_R to tensor format; full replication in tensor group
1051 DO i_cell_r = 1, nimages_scf
1052 CALL fm_to_local_tensor(bs_env%fm_chi_R_t(i_cell_r, 1), bs_env%mat_RI_RI%matrix, &
1053 bs_env%mat_RI_RI_tensor%matrix, mi_vtr_mi_r(i_cell_r), bs_env)
1054 END DO
1055
1056 CALL timestop(handle)
1057
1058 END SUBROUTINE get_minv_vtr_minv_r
1059
1060! **************************************************************************************************
1061!> \brief ...
1062!> \param t_1d ...
1063! **************************************************************************************************
1064 SUBROUTINE destroy_t_1d(t_1d)
1065 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_1d
1066
1067 CHARACTER(LEN=*), PARAMETER :: routinen = 'destroy_t_1d'
1068
1069 INTEGER :: handle, i
1070
1071 CALL timeset(routinen, handle)
1072
1073 DO i = 1, SIZE(t_1d)
1074 CALL dbt_destroy(t_1d(i))
1075 END DO
1076 DEALLOCATE (t_1d)
1077
1078 CALL timestop(handle)
1079
1080 END SUBROUTINE destroy_t_1d
1081
1082! **************************************************************************************************
1083!> \brief ...
1084!> \param t_2d ...
1085! **************************************************************************************************
1086 SUBROUTINE destroy_t_2d(t_2d)
1087 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_2d
1088
1089 CHARACTER(LEN=*), PARAMETER :: routinen = 'destroy_t_2d'
1090
1091 INTEGER :: handle, i, j
1092
1093 CALL timeset(routinen, handle)
1094
1095 DO i = 1, SIZE(t_2d, 1)
1096 DO j = 1, SIZE(t_2d, 2)
1097 CALL dbt_destroy(t_2d(i, j))
1098 END DO
1099 END DO
1100 DEALLOCATE (t_2d)
1101
1102 CALL timestop(handle)
1103
1104 END SUBROUTINE destroy_t_2d
1105
1106! **************************************************************************************************
1107!> \brief ...
1108!> \param t_W ...
1109!> \param W_R ...
1110!> \param bs_env ...
1111!> \param i_task_Delta_R_local ...
1112! **************************************************************************************************
1113 SUBROUTINE contract_w(t_W, W_R, bs_env, i_task_Delta_R_local)
1114 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_w
1115 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: w_r
1116 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1117 INTEGER :: i_task_delta_r_local
1118
1119 CHARACTER(LEN=*), PARAMETER :: routinen = 'contract_W'
1120
1121 INTEGER :: handle, i_cell_delta_r, i_cell_r1, &
1122 i_cell_r2, i_cell_r2_m_r1, i_cell_s1, &
1123 i_cell_s1_m_r1_p_r2
1124 INTEGER, DIMENSION(3) :: cell_dr, cell_r1, cell_r2, cell_r2_m_r1, &
1125 cell_s1, cell_s1_m_r2_p_r1
1126 LOGICAL :: cell_found
1127 TYPE(dbt_type) :: t_3c_int, t_w_tmp
1128
1129 CALL timeset(routinen, handle)
1130
1131 CALL dbt_create(bs_env%t_RI__AO_AO, t_w_tmp)
1132 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_int)
1133
1134 i_cell_delta_r = bs_env%task_Delta_R(i_task_delta_r_local)
1135
1136 DO i_cell_r1 = 1, bs_env%nimages_3c
1137
1138 cell_r1(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_r1)
1139 cell_dr(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_delta_r)
1140
1141 ! S_1 = R_1 + ΔR (from ΔR = S_1 - R_1)
1142 CALL add_r(cell_r1, cell_dr, bs_env%index_to_cell_3c, cell_s1, &
1143 cell_found, bs_env%cell_to_index_3c, i_cell_s1)
1144 IF (.NOT. cell_found) cycle
1145
1146 DO i_cell_r2 = 1, bs_env%nimages_scf_desymm
1147
1148 cell_r2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_r2)
1149
1150 ! R_2 - R_1
1151 CALL add_r(cell_r2, -cell_r1, bs_env%index_to_cell_3c, cell_r2_m_r1, &
1152 cell_found, bs_env%cell_to_index_3c, i_cell_r2_m_r1)
1153 IF (.NOT. cell_found) cycle
1154
1155 ! S_1 - R_1 + R_2
1156 CALL add_r(cell_s1, cell_r2_m_r1, bs_env%index_to_cell_3c, cell_s1_m_r2_p_r1, &
1157 cell_found, bs_env%cell_to_index_3c, i_cell_s1_m_r1_p_r2)
1158 IF (.NOT. cell_found) cycle
1159
1160 CALL get_t_3c_int(t_3c_int, bs_env, i_cell_s1_m_r1_p_r2, i_cell_r2_m_r1)
1161
1162 ! M^W_σ0,νS1,PR1 = sum_QR2 ( σ0 νS1 | QR1-R2 ) W_QP^R2
1163 ! = sum_QR2 ( σR2-R1 νS1-R1+R2 | Q0 ) W_QP^R2
1164 ! for ΔR = S_1 - R_1
1165 CALL dbt_contract(alpha=1.0_dp, &
1166 tensor_1=w_r(i_cell_r2), &
1167 tensor_2=t_3c_int, &
1168 beta=0.0_dp, &
1169 tensor_3=t_w_tmp, &
1170 contract_1=[1], notcontract_1=[2], map_1=[1], &
1171 contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3], &
1172 filter_eps=bs_env%eps_filter)
1173
1174 ! reorder tensor
1175 CALL dbt_copy(t_w_tmp, t_w(i_cell_s1, i_cell_r1), order=[1, 2, 3], &
1176 move_data=.true., summation=.true.)
1177
1178 END DO ! i_cell_R2
1179
1180 END DO ! i_cell_R1
1181
1182 CALL dbt_destroy(t_w_tmp)
1183 CALL dbt_destroy(t_3c_int)
1184
1185 CALL timestop(handle)
1186
1187 END SUBROUTINE contract_w
1188
1189! **************************************************************************************************
1190!> \brief ...
1191!> \param Sigma_R ...
1192!> \param t_W ...
1193!> \param G_S ...
1194!> \param i_task_Delta_R_local ...
1195!> \param bs_env ...
1196!> \param occ ...
1197!> \param vir ...
1198!> \param clear_t_W ...
1199!> \param fill_skip ...
1200! **************************************************************************************************
1201 SUBROUTINE contract_to_sigma(Sigma_R, t_W, G_S, i_task_Delta_R_local, bs_env, occ, vir, &
1202 clear_t_W, fill_skip)
1203 TYPE(dbt_type), DIMENSION(:) :: sigma_r
1204 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_w
1205 TYPE(dbt_type), DIMENSION(:) :: g_s
1206 INTEGER :: i_task_delta_r_local
1207 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1208 LOGICAL :: occ, vir, clear_t_w, fill_skip
1209
1210 CHARACTER(LEN=*), PARAMETER :: routinen = 'contract_to_Sigma'
1211
1212 INTEGER :: handle, handle2, i_cell_delta_r, i_cell_m_r1, i_cell_r, i_cell_r1, &
1213 i_cell_r1_minus_r, i_cell_s1, i_cell_s1_minus_r, i_cell_s1_p_s2_m_r1, i_cell_s2
1214 INTEGER(KIND=int_8) :: flop, flop_tmp
1215 INTEGER, DIMENSION(3) :: cell_dr, cell_m_r1, cell_r, cell_r1, &
1216 cell_r1_minus_r, cell_s1, &
1217 cell_s1_minus_r, cell_s1_p_s2_m_r1, &
1218 cell_s2
1219 LOGICAL :: cell_found
1220 REAL(kind=dp) :: sign_sigma
1221 TYPE(dbt_type) :: t_3c_int, t_g, t_g_2
1222
1223 CALL timeset(routinen, handle)
1224
1225 cpassert(occ .EQV. (.NOT. vir))
1226 IF (occ) sign_sigma = -1.0_dp
1227 IF (vir) sign_sigma = 1.0_dp
1228
1229 CALL dbt_create(bs_env%t_RI_AO__AO, t_g)
1230 CALL dbt_create(bs_env%t_RI_AO__AO, t_g_2)
1231 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_int)
1232
1233 i_cell_delta_r = bs_env%task_Delta_R(i_task_delta_r_local)
1234
1235 flop = 0_int_8
1236
1237 DO i_cell_r1 = 1, bs_env%nimages_3c
1238
1239 cell_r1(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_r1)
1240 cell_dr(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_delta_r)
1241
1242 ! S_1 = R_1 + ΔR (from ΔR = S_1 - R_1)
1243 CALL add_r(cell_r1, cell_dr, bs_env%index_to_cell_3c, cell_s1, cell_found, &
1244 bs_env%cell_to_index_3c, i_cell_s1)
1245 IF (.NOT. cell_found) cycle
1246
1247 DO i_cell_s2 = 1, bs_env%nimages_scf_desymm
1248
1249 IF (bs_env%skip_DR_R1_S2_Gx3c_Sigma(i_task_delta_r_local, i_cell_r1, i_cell_s2)) cycle
1250
1251 cell_s2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_s2)
1252 cell_m_r1(1:3) = -cell_r1(1:3)
1253 cell_s1_p_s2_m_r1(1:3) = cell_s1(1:3) + cell_s2(1:3) - cell_r1(1:3)
1254
1255 CALL is_cell_in_index_to_cell(cell_m_r1, bs_env%index_to_cell_3c, cell_found)
1256 IF (.NOT. cell_found) cycle
1257
1258 CALL is_cell_in_index_to_cell(cell_s1_p_s2_m_r1, bs_env%index_to_cell_3c, cell_found)
1259 IF (.NOT. cell_found) cycle
1260
1261 i_cell_m_r1 = bs_env%cell_to_index_3c(cell_m_r1(1), cell_m_r1(2), cell_m_r1(3))
1262 i_cell_s1_p_s2_m_r1 = bs_env%cell_to_index_3c(cell_s1_p_s2_m_r1(1), &
1263 cell_s1_p_s2_m_r1(2), &
1264 cell_s1_p_s2_m_r1(3))
1265
1266 CALL timeset(routinen//"_3c_x_G", handle2)
1267
1268 CALL get_t_3c_int(t_3c_int, bs_env, i_cell_m_r1, i_cell_s1_p_s2_m_r1)
1269
1270 ! M_λ0,νS1,PR1 = sum_µS2 ( λ0 µS1-S2 | PR1 ) G^occ/vir_µν^S2(i|τ|)
1271 ! = sum_µS2 ( λ-R1 µS1-S2-R1 | P0 ) G^occ/vir_µν^S2(i|τ|)
1272 ! for ΔR = S_1 - R_1
1273 CALL dbt_contract(alpha=1.0_dp, &
1274 tensor_1=g_s(i_cell_s2), &
1275 tensor_2=t_3c_int, &
1276 beta=1.0_dp, &
1277 tensor_3=t_g, &
1278 contract_1=[2], notcontract_1=[1], map_1=[3], &
1279 contract_2=[3], notcontract_2=[1, 2], map_2=[1, 2], &
1280 filter_eps=bs_env%eps_filter, flop=flop_tmp)
1281
1282 IF (flop_tmp == 0_int_8 .AND. fill_skip) THEN
1283 bs_env%skip_DR_R1_S2_Gx3c_Sigma(i_task_delta_r_local, i_cell_r1, i_cell_s2) = .true.
1284 END IF
1285
1286 CALL timestop(handle2)
1287
1288 END DO ! i_cell_S2
1289
1290 CALL dbt_copy(t_g, t_g_2, order=[1, 3, 2], move_data=.true.)
1291
1292 CALL timeset(routinen//"_contract", handle2)
1293
1294 DO i_cell_r = 1, bs_env%nimages_scf_desymm
1295
1296 IF (bs_env%skip_DR_R1_R_MxM_Sigma(i_task_delta_r_local, i_cell_r1, i_cell_r)) cycle
1297
1298 cell_r = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_r)
1299
1300 ! R_1 - R
1301 CALL add_r(cell_r1, -cell_r, bs_env%index_to_cell_3c, cell_r1_minus_r, &
1302 cell_found, bs_env%cell_to_index_3c, i_cell_r1_minus_r)
1303 IF (.NOT. cell_found) cycle
1304
1305 ! S_1 - R
1306 CALL add_r(cell_s1, -cell_r, bs_env%index_to_cell_3c, cell_s1_minus_r, &
1307 cell_found, bs_env%cell_to_index_3c, i_cell_s1_minus_r)
1308 IF (.NOT. cell_found) cycle
1309
1310 ! Σ_λσ^R = sum_PR1νS1 M^G_λ0,νS1,PR1 M^W_σR,νS1,PR1, where
1311 ! M^G_λ0,νS1,PR1 = sum_µS2 (λ0 µS1-S2 | PR1) G_µν^S2
1312 ! M^W_σR,νS1,PR1 = sum_QR2 (σR νS1 | QR1-R2) W_PQ^R2 = M^W_σ0,νS1-R,PR1-R
1313 CALL dbt_contract(alpha=sign_sigma, &
1314 tensor_1=t_g_2, &
1315 tensor_2=t_w(i_cell_s1_minus_r, i_cell_r1_minus_r), &
1316 beta=1.0_dp, &
1317 tensor_3=sigma_r(i_cell_r), &
1318 contract_1=[1, 2], notcontract_1=[3], map_1=[1], &
1319 contract_2=[1, 2], notcontract_2=[3], map_2=[2], &
1320 filter_eps=bs_env%eps_filter, flop=flop_tmp)
1321
1322 flop = flop + flop_tmp
1323
1324 IF (flop_tmp == 0_int_8 .AND. fill_skip) THEN
1325 bs_env%skip_DR_R1_R_MxM_Sigma(i_task_delta_r_local, i_cell_r1, i_cell_r) = .true.
1326 END IF
1327
1328 END DO ! i_cell_R
1329
1330 CALL dbt_clear(t_g_2)
1331
1332 CALL timestop(handle2)
1333
1334 END DO ! i_cell_R1
1335
1336 IF (vir .AND. flop == 0_int_8) bs_env%skip_DR_Sigma(i_task_delta_r_local) = .true.
1337
1338 ! release memory
1339 IF (clear_t_w) THEN
1340 DO i_cell_s1 = 1, bs_env%nimages_3c
1341 DO i_cell_r1 = 1, bs_env%nimages_3c
1342 CALL dbt_clear(t_w(i_cell_s1, i_cell_r1))
1343 END DO
1344 END DO
1345 END IF
1346
1347 CALL dbt_destroy(t_g)
1348 CALL dbt_destroy(t_g_2)
1349 CALL dbt_destroy(t_3c_int)
1350
1351 CALL timestop(handle)
1352
1353 END SUBROUTINE contract_to_sigma
1354
1355! **************************************************************************************************
1356!> \brief ...
1357!> \param fm_W_R ...
1358!> \param W_R ...
1359!> \param bs_env ...
1360! **************************************************************************************************
1361 SUBROUTINE fm_mwm_r_t_to_local_tensor_w_r(fm_W_R, W_R, bs_env)
1362 TYPE(cp_fm_type), DIMENSION(:) :: fm_w_r
1363 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: w_r
1364 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1365
1366 CHARACTER(LEN=*), PARAMETER :: routinen = 'fm_MWM_R_t_to_local_tensor_W_R'
1367
1368 INTEGER :: handle, i_cell_r
1369
1370 CALL timeset(routinen, handle)
1371
1372 ! communicate fm_W_R to tensor W_R; full replication in tensor group
1373 DO i_cell_r = 1, bs_env%nimages_scf_desymm
1374 CALL fm_to_local_tensor(fm_w_r(i_cell_r), bs_env%mat_RI_RI%matrix, &
1375 bs_env%mat_RI_RI_tensor%matrix, w_r(i_cell_r), bs_env)
1376 END DO
1377
1378 CALL timestop(handle)
1379
1380 END SUBROUTINE fm_mwm_r_t_to_local_tensor_w_r
1381
1382! **************************************************************************************************
1383!> \brief ...
1384!> \param bs_env ...
1385! **************************************************************************************************
1386 SUBROUTINE compute_qp_energies(bs_env)
1387 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1388
1389 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_QP_energies'
1390
1391 INTEGER :: handle, ikp, ispin, j_t
1392 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sigma_x_ikp_n
1393 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sigma_c_ikp_n_freq, sigma_c_ikp_n_time
1394 TYPE(cp_cfm_type) :: cfm_mo_coeff
1395
1396 CALL timeset(routinen, handle)
1397
1398 CALL cp_cfm_create(cfm_mo_coeff, bs_env%fm_s_Gamma%matrix_struct)
1399 ALLOCATE (sigma_x_ikp_n(bs_env%n_ao))
1400 ALLOCATE (sigma_c_ikp_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
1401 ALLOCATE (sigma_c_ikp_n_freq(bs_env%n_ao, bs_env%num_time_freq_points, 2))
1402
1403 DO ispin = 1, bs_env%n_spin
1404
1405 DO ikp = 1, bs_env%nkp_bs_and_DOS
1406
1407 ! 1. get C_µn(k)
1408 CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mo_coeff)
1409
1410 ! 2. Σ^x_µν(k) = sum_R Σ^x_µν^R e^ikR
1411 ! Σ^x_nn(k) = sum_µν C^*_µn(k) Σ^x_µν(k) C_νn(k)
1412 CALL trafo_to_k_and_nn(bs_env%fm_Sigma_x_R, sigma_x_ikp_n, cfm_mo_coeff, bs_env, ikp)
1413
1414 ! 3. Σ^c_µν(k,+/-i|τ_j|) = sum_R Σ^c_µν^R(+/-i|τ_j|) e^ikR
1415 ! Σ^c_nn(k,+/-i|τ_j|) = sum_µν C^*_µn(k) Σ^c_µν(k,+/-i|τ_j|) C_νn(k)
1416 DO j_t = 1, bs_env%num_time_freq_points
1417 CALL trafo_to_k_and_nn(bs_env%fm_Sigma_c_R_pos_tau(:, j_t, ispin), &
1418 sigma_c_ikp_n_time(:, j_t, 1), cfm_mo_coeff, bs_env, ikp)
1419 CALL trafo_to_k_and_nn(bs_env%fm_Sigma_c_R_neg_tau(:, j_t, ispin), &
1420 sigma_c_ikp_n_time(:, j_t, 2), cfm_mo_coeff, bs_env, ikp)
1421 END DO
1422
1423 ! 4. Σ^c_nn(k_i,iω) = ∫ from -∞ to ∞ dτ e^-iωτ Σ^c_nn(k_i,iτ)
1424 CALL time_to_freq(bs_env, sigma_c_ikp_n_time, sigma_c_ikp_n_freq, ispin)
1425
1426 ! 5. Analytic continuation Σ^c_nn(k_i,iω) -> Σ^c_nn(k_i,ϵ) and
1427 ! ϵ_nk_i^GW = ϵ_nk_i^DFT + Σ^c_nn(k_i,ϵ) + Σ^x_nn(k_i) - v^xc_nn(k_i)
1428 CALL analyt_conti_and_print(bs_env, sigma_c_ikp_n_freq, sigma_x_ikp_n, &
1429 bs_env%v_xc_n(:, ikp, ispin), &
1430 bs_env%eigenval_scf(:, ikp, ispin), ikp, ispin)
1431
1432 END DO ! ikp
1433
1434 END DO ! ispin
1435
1436 CALL get_all_vbm_cbm_bandgaps(bs_env)
1437
1438 CALL cp_cfm_release(cfm_mo_coeff)
1439
1440 CALL timestop(handle)
1441
1442 END SUBROUTINE compute_qp_energies
1443
1444! **************************************************************************************************
1445!> \brief ...
1446!> \param fm_rs ...
1447!> \param array_ikp_n ...
1448!> \param cfm_mo_coeff ...
1449!> \param bs_env ...
1450!> \param ikp ...
1451! **************************************************************************************************
1452 SUBROUTINE trafo_to_k_and_nn(fm_rs, array_ikp_n, cfm_mo_coeff, bs_env, ikp)
1453 TYPE(cp_fm_type), DIMENSION(:) :: fm_rs
1454 REAL(kind=dp), DIMENSION(:) :: array_ikp_n
1455 TYPE(cp_cfm_type) :: cfm_mo_coeff
1456 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1457 INTEGER :: ikp
1458
1459 CHARACTER(LEN=*), PARAMETER :: routinen = 'trafo_to_k_and_nn'
1460
1461 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cdiag
1462 INTEGER :: handle
1463 TYPE(cp_cfm_type) :: cfm_ikp
1464
1465 CALL timeset(routinen, handle)
1466
1467 CALL cp_cfm_create(cfm_ikp, cfm_mo_coeff%matrix_struct)
1468
1469 ! Σ_µν(k_i) = sum_R e^ik_iR Σ_µν^R
1470 CALL fm_rs_to_kp(cfm_ikp, fm_rs, bs_env%kpoints_DOS, ikp)
1471
1472 ! Σ_nm(k_i) = sum_µν C^*_µn(k_i) Σ_µν(k_i) C_νn(k_i)
1473 CALL cfm_contract_aba(cfm_mo_coeff, cfm_ikp)
1474
1475 ! get Σ_nn(k_i) which is a real quantity as Σ^x and Σ^c(iτ) is Hermitian
1476 ALLOCATE (cdiag(SIZE(array_ikp_n)))
1477 CALL cp_cfm_get_diag(cfm_ikp, cdiag)
1478 array_ikp_n = real(cdiag, kind=dp)
1479
1480 DEALLOCATE (cdiag)
1481 CALL cp_cfm_release(cfm_ikp)
1482
1483 CALL timestop(handle)
1484
1485 END SUBROUTINE trafo_to_k_and_nn
1486
1487END MODULE gw_tensor_small_cell_full_kp
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
mathlib::gemm_square
Definition mathlib.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::pasquier2025
integer, save, public pasquier2025
Definition bibliography.F:36
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_get_diag
subroutine, public cp_cfm_get_diag(matrix, diag)
returns the diagonal of a complex full matrix: diag(i)= A_{ii}. Each diagonal entry is owned by one p...
Definition cp_cfm_basic_linalg.F:1479
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
gw_tensor_small_cell_full_kp
Definition gw_tensor_small_cell_full_kp.F:13
gw_tensor_small_cell_full_kp::gw_calc_tensor_small_cell_full_kp
subroutine, public gw_calc_tensor_small_cell_full_kp(qs_env, bs_env)
Perform GW band structure calculation.
Definition gw_tensor_small_cell_full_kp.F:78
gw_utils_communication
Definition gw_utils_communication.F:13
gw_utils_communication::fm_to_local_array
subroutine, public fm_to_local_array(fm_s, array_s, weight, add)
...
Definition gw_utils_communication.F:764
gw_utils_communication::fm_to_local_tensor
subroutine, public fm_to_local_tensor(fm_global, mat_global, mat_local, tensor, bs_env, atom_ranges)
...
Definition gw_utils_communication.F:67
gw_utils_communication::local_dbt_to_global_fm
subroutine, public local_dbt_to_global_fm(t_r, fm_r, mat_global, mat_local, bs_env)
...
Definition gw_utils_communication.F:916
gw_utils_communication::local_array_to_fm
subroutine, public local_array_to_fm(array_s, fm_s, weight, add)
...
Definition gw_utils_communication.F:841
gw_utils_fm
Full-matrix operations not provided by the CP2K FM packages.
Definition gw_utils_fm.F:13
gw_utils_fm::local_add_on_diag
subroutine, public local_add_on_diag(matrix, alpha)
Add a scalar to the diagonal of a local complex matrix.
Definition gw_utils_fm.F:90
gw_utils_fm::cfm_contract_aba
subroutine, public cfm_contract_aba(matrix_a, matrix_b, matrix_c)
Computes A^H B A for complex full matrices.
Definition gw_utils_fm.F:56
gw_utils_fm::local_complex_power
subroutine, public local_complex_power(matrix, exponent, eps, cond_nr, min_ev, max_ev)
Compute a spectral power of a local complex Hermitian matrix. Eigenvalues not larger than eps are dis...
Definition gw_utils_fm.F:115
gw_utils
Definition gw_utils.F:15
gw_utils::get_v_tr_r
subroutine, public get_v_tr_r(v_tr_r, pot_type, regularization_ri, bs_env, qs_env)
...
Definition gw_utils.F:3228
gw_utils::time_to_freq
subroutine, public time_to_freq(bs_env, sigma_c_n_time, sigma_c_n_freq, ispin)
...
Definition gw_utils.F:3499
gw_utils::analyt_conti_and_print
subroutine, public analyt_conti_and_print(bs_env, sigma_c_ikp_n_freq, sigma_x_ikp_n, v_xc_ikp_n, eigenval_scf, ikp, ispin)
...
Definition gw_utils.F:3563
gw_utils::add_r
subroutine, public add_r(cell_1, cell_2, index_to_cell, cell_1_plus_2, cell_found, cell_to_index, i_cell_1_plus_2)
...
Definition gw_utils.F:2863
gw_utils::de_init_bs_env
subroutine, public de_init_bs_env(qs_env, bs_env)
Releases the memory-heavy GW intermediates that cannot be freed in bs_env_release,...
Definition gw_utils.F:3411
gw_utils::is_cell_in_index_to_cell
subroutine, public is_cell_in_index_to_cell(cell, index_to_cell, cell_found)
...
Definition gw_utils.F:2902
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
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_small_cell
subroutine, public build_2c_coulomb_matrix_kp_small_cell(v_k, qs_env, kpoints, size_lattice_sum, basis_type, ikp_start, ikp_end)
...
Definition kpoint_coulomb_2c.F:698
kpoint_k_r_trafo_simple
Implements transformations from k-space to R-space for Fortran array matrices.
Definition kpoint_k_r_trafo_simple.F:16
kpoint_k_r_trafo_simple::rs_to_kp
subroutine, public rs_to_kp(rs_real, ks_complex, index_to_cell, xkp, deriv_direction, hmat)
Integrate RS matrices (stored as Fortran array) into a kpoint matrix at given kp.
Definition kpoint_k_r_trafo_simple.F:230
kpoint_k_r_trafo_simple::fm_add_kp_to_all_rs
subroutine, public fm_add_kp_to_all_rs(cfm_kp, fm_rs, kpoints, ikp)
Adds given kpoint matrix to a single rs matrix.
Definition kpoint_k_r_trafo_simple.F:398
kpoint_k_r_trafo_simple::add_kp_to_all_rs
subroutine, public add_kp_to_all_rs(array_kp, array_rs, kpoints, ikp, index_to_cell_ext)
Adds given kpoint matrix to all rs matrices.
Definition kpoint_k_r_trafo_simple.F:359
kpoint_k_r_trafo_simple::fm_rs_to_kp
subroutine, public fm_rs_to_kp(cfm_kp, fm_rs, kpoints, ikp)
Transforms array of fm RS matrices into cfm k-space matrix, at given kpoint index.
Definition kpoint_k_r_trafo_simple.F:269
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
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::get_all_vbm_cbm_bandgaps
subroutine, public get_all_vbm_cbm_bandgaps(bs_env)
...
Definition post_scf_bandstructure_utils.F:3055
qs_environment_types
Definition qs_environment_types.F:14
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
post_scf_bandstructure_types::post_scf_bandstructure_type
Definition post_scf_bandstructure_types.F:156
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224