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