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