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