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