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_basic_linalg, ONLY: cp_cfm_get_diag
17 : USE cp_cfm_types, ONLY: cp_cfm_create,&
18 : cp_cfm_get_info,&
19 : cp_cfm_release,&
20 : cp_cfm_to_cfm,&
21 : cp_cfm_type
22 : USE cp_fm_redistribution, ONLY: cp_fm_to_replicated_array
23 : USE cp_fm_types, ONLY: cp_fm_set_all,&
24 : cp_fm_set_submatrix,&
25 : cp_fm_type
26 : USE dbt_api, ONLY: dbt_clear,&
27 : dbt_contract,&
28 : dbt_copy,&
29 : dbt_create,&
30 : dbt_destroy,&
31 : dbt_type
32 : USE dbt_matrix_redistribution, ONLY: cp_dbt_to_fm,&
33 : cp_fm_to_dbt
34 : USE gw_utils, ONLY: add_R,&
35 : analyt_conti_and_print,&
36 : de_init_bs_env,&
37 : get_V_tr_R,&
38 : is_cell_in_index_to_cell
39 : USE gw_utils_fm, ONLY: cfm_contract_ABA,&
40 : local_add_on_diag,&
41 : local_complex_power
42 : USE kinds, ONLY: dp,&
43 : int_8
44 : USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp_small_cell
45 : USE kpoint_k_r_trafo_simple, ONLY: add_kp_to_all_rs,&
46 : fm_add_kp_to_all_rs,&
47 : fm_rs_to_kp,&
48 : rs_to_kp
49 : USE machine, ONLY: m_walltime
50 : USE mathconstants, ONLY: z_one,&
51 : z_zero
52 : USE mathlib, ONLY: gemm_square
53 : USE parallel_gemm_api, ONLY: parallel_gemm
54 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
55 : USE post_scf_bandstructure_utils, ONLY: get_all_VBM_CBM_bandgaps
56 : USE qs_environment_types, ONLY: qs_environment_type
57 : USE time_frequency_grids, ONLY: negative_time,&
58 : positive_time,&
59 : transform_time_to_frequency
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 cp_dbt_to_fm(t_chi_R, bs_env%fm_chi_R_t(:, i_t), bs_env%mat_RI_RI, bs_env%mat_RI_RI_tensor, &
200 104 : bs_env%para_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 cp_fm_to_dbt(bs_env%fm_G_S(i_cell_S), bs_env%mat_ao_ao%matrix, bs_env%mat_ao_ao_tensor%matrix, &
459 : G_S(i_cell_S), bs_env%para_env, bs_env%para_env_tensor%num_pe, &
460 4320 : bs_env%eps_filter)
461 : END DO
462 :
463 432 : CALL timestop(handle)
464 :
465 432 : END SUBROUTINE G_occ_vir
466 :
467 : ! **************************************************************************************************
468 : !> \brief ...
469 : !> \param t_Gocc ...
470 : !> \param t_Gvir ...
471 : !> \param t_chi_R ...
472 : !> \param bs_env ...
473 : !> \param i_task_Delta_R_local ...
474 : ! **************************************************************************************************
475 560 : SUBROUTINE contract_M_occ_vir_to_chi(t_Gocc, t_Gvir, t_chi_R, bs_env, i_task_Delta_R_local)
476 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_Gocc, t_Gvir
477 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_chi_R
478 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
479 : INTEGER :: i_task_Delta_R_local
480 :
481 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_M_occ_vir_to_chi'
482 :
483 : INTEGER :: handle, i_cell_Delta_R, i_cell_R, &
484 : i_cell_R1, i_cell_R1_minus_R, &
485 : i_cell_R2, i_cell_R2_minus_R
486 : INTEGER(KIND=int_8) :: flop, flop_tmp
487 : INTEGER, DIMENSION(3) :: cell_DR, cell_R, cell_R1, &
488 : cell_R1_minus_R, cell_R2, &
489 : cell_R2_minus_R
490 : LOGICAL :: cell_found
491 9520 : TYPE(dbt_type) :: t_Gocc_2, t_Gvir_2
492 :
493 560 : CALL timeset(routineN, handle)
494 :
495 560 : CALL dbt_create(bs_env%t_RI__AO_AO, t_Gocc_2)
496 560 : CALL dbt_create(bs_env%t_RI__AO_AO, t_Gvir_2)
497 :
498 560 : flop = 0_int_8
499 :
500 : ! χ_PQ^R(iτ) = sum_λR1,νR2 M^occ_λR1,νR2,P0 M^vir_νR2,λR1,QR
501 5600 : DO i_cell_R = 1, bs_env%nimages_scf_desymm
502 :
503 63875 : DO i_cell_R2 = 1, bs_env%nimages_3c
504 :
505 58275 : IF (bs_env%skip_DR_R_R2_MxM_chi(i_task_Delta_R_local, i_cell_R2, i_cell_R)) CYCLE
506 :
507 49778 : i_cell_Delta_R = bs_env%task_Delta_R(i_task_Delta_R_local)
508 :
509 199112 : cell_R(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_R)
510 199112 : cell_R2(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_R2)
511 199112 : cell_DR(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_Delta_R)
512 :
513 : ! R_1 = R_2 + ΔR (from ΔR = R_2 - R_1)
514 : CALL add_R(cell_R2, cell_DR, bs_env%index_to_cell_3c, cell_R1, &
515 49778 : cell_found, bs_env%cell_to_index_3c, i_cell_R1)
516 49778 : IF (.NOT. cell_found) CYCLE
517 :
518 : ! R_1 - R
519 : CALL add_R(cell_R1, -cell_R, bs_env%index_to_cell_3c, cell_R1_minus_R, &
520 91688 : cell_found, bs_env%cell_to_index_3c, i_cell_R1_minus_R)
521 22922 : IF (.NOT. cell_found) CYCLE
522 :
523 : ! R_2 - R
524 : CALL add_R(cell_R2, -cell_R, bs_env%index_to_cell_3c, cell_R2_minus_R, &
525 51880 : cell_found, bs_env%cell_to_index_3c, i_cell_R2_minus_R)
526 12970 : IF (.NOT. cell_found) CYCLE
527 :
528 : ! reorder tensors for efficient contraction to χ_PQ^R
529 8328 : CALL dbt_copy(t_Gocc(i_cell_R1, i_cell_R2), t_Gocc_2, order=[1, 3, 2])
530 8328 : CALL dbt_copy(t_Gvir(i_cell_R2_minus_R, i_cell_R1_minus_R), t_Gvir_2)
531 :
532 : ! χ_PQ^R(iτ) = sum_λR1,νR2 M^occ_λR1,νR2,P0 M^vir_νR2,λR1,QR
533 : CALL dbt_contract(alpha=bs_env%spin_degeneracy, &
534 : tensor_1=t_Gocc_2, tensor_2=t_Gvir_2, &
535 : beta=1.0_dp, tensor_3=t_chi_R(i_cell_R), &
536 : contract_1=[2, 3], notcontract_1=[1], map_1=[1], &
537 : contract_2=[2, 3], notcontract_2=[1], map_2=[2], &
538 8328 : filter_eps=bs_env%eps_filter, move_data=.TRUE., flop=flop_tmp)
539 :
540 8328 : IF (flop_tmp == 0_int_8) bs_env%skip_DR_R_R2_MxM_chi(i_task_Delta_R_local, &
541 3485 : i_cell_R2, i_cell_R) = .TRUE.
542 :
543 71643 : flop = flop + flop_tmp
544 :
545 : END DO ! i_cell_R2
546 :
547 : END DO ! i_cell_R
548 :
549 560 : IF (flop == 0_int_8) bs_env%skip_DR_chi(i_task_Delta_R_local) = .TRUE.
550 :
551 : ! remove all data from t_Gocc and t_Gvir to safe memory
552 7035 : DO i_cell_R1 = 1, bs_env%nimages_3c
553 83510 : DO i_cell_R2 = 1, bs_env%nimages_3c
554 76475 : CALL dbt_clear(t_Gocc(i_cell_R1, i_cell_R2))
555 82950 : CALL dbt_clear(t_Gvir(i_cell_R1, i_cell_R2))
556 : END DO
557 : END DO
558 :
559 560 : CALL dbt_destroy(t_Gocc_2)
560 560 : CALL dbt_destroy(t_Gvir_2)
561 :
562 560 : CALL timestop(handle)
563 :
564 560 : END SUBROUTINE contract_M_occ_vir_to_chi
565 :
566 : ! **************************************************************************************************
567 : !> \brief ...
568 : !> \param bs_env ...
569 : !> \param qs_env ...
570 : ! **************************************************************************************************
571 16 : SUBROUTINE compute_W_real_space(bs_env, qs_env)
572 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
573 : TYPE(qs_environment_type), POINTER :: qs_env
574 :
575 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_W_real_space'
576 :
577 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chi_k_w, eps_k_w, W_k_w
578 16 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: M_inv, M_inv_V_sqrt, V_sqrt
579 : INTEGER :: handle, i_cell_R, i_t, ikp, ikp_local, &
580 : j_w, n_RI, nimages_scf_desymm
581 : REAL(KIND=dp) :: t1, weight_ij
582 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: chi_R, MWM_R, W_R
583 :
584 16 : CALL timeset(routineN, handle)
585 :
586 16 : n_RI = bs_env%n_RI
587 16 : nimages_scf_desymm = bs_env%nimages_scf_desymm
588 :
589 128 : ALLOCATE (chi_k_w(n_RI, n_RI), eps_k_w(n_RI, n_RI), W_k_w(n_RI, n_RI))
590 : ALLOCATE (chi_R(n_RI, n_RI, nimages_scf_desymm), W_R(n_RI, n_RI, nimages_scf_desymm), &
591 176 : MWM_R(n_RI, n_RI, nimages_scf_desymm))
592 :
593 16 : t1 = m_walltime()
594 :
595 16 : CALL compute_Minv_and_Vsqrt(bs_env, qs_env, M_inv_V_sqrt, M_inv, V_sqrt)
596 :
597 16 : IF (bs_env%unit_nr > 0) THEN
598 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
599 8 : 'Computed V_PQ(k),', 'Execution time', m_walltime() - t1, ' s'
600 8 : WRITE (bs_env%unit_nr, '(A)') ' '
601 : END IF
602 :
603 16 : t1 = m_walltime()
604 :
605 120 : DO j_w = 1, bs_env%num_time_freq_points
606 :
607 : ! χ_PQ^R(iτ) -> χ_PQ^R(iω_j) (which is stored in chi_R, single ω_j from j_w loop)
608 104 : chi_R(:, :, :) = 0.0_dp
609 808 : DO i_t = 1, bs_env%num_time_freq_points
610 704 : weight_ij = bs_env%time_frequency_grid%cosine_time_to_frequency(j_w, i_t)
611 :
612 808 : CALL cp_fm_to_replicated_array(bs_env%fm_chi_R_t(:, i_t), chi_R, weight_ij, add=.TRUE.)
613 : END DO
614 :
615 104 : ikp_local = 0
616 104 : W_R(:, :, :) = 0.0_dp
617 66664 : DO ikp = 1, bs_env%nkp_chi_eps_W_orig_plus_extra
618 :
619 : ! trivial parallelization over k-points
620 66560 : IF (MODULO(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
621 :
622 33280 : ikp_local = ikp_local + 1
623 :
624 : ! 1. χ_PQ^R(iω_j) -> χ_PQ(iω_j,k)
625 : CALL rs_to_kp(chi_R, chi_k_w, bs_env%kpoints_scf_desymm%index_to_cell, &
626 33280 : bs_env%kpoints_chi_eps_W%xkp(1:3, ikp))
627 :
628 : ! 2. remove negative eigenvalues from χ_PQ(iω,k)
629 33280 : CALL local_complex_power(chi_k_w, 1.0_dp, bs_env%eps_eigval_mat_RI)
630 :
631 : ! 3. ε(iω_j,k_i) = I + V^0.5(k_i)*M^-1(k_i)*χ(iω_j,k_i)*M^-1(k_i)*V^0.5(k_i)
632 : ! 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)
633 : CALL gemm_square(M_inv_V_sqrt(:, :, ikp_local), 'C', chi_k_w, 'N', &
634 33280 : M_inv_V_sqrt(:, :, ikp_local), 'N', eps_k_w)
635 :
636 : ! b) ε(iω_j,k_i) = eps_work + I
637 33280 : CALL local_add_on_diag(eps_k_w, z_one)
638 :
639 : ! 4. W(iω_j,k_i) = V^0.5(k_i)*(ε^-1(iω_j,k_i)-I)*V^0.5(k_i)
640 : ! a) Invert ε(iω_j,k_i) using its Cholesky decomposition
641 33280 : CALL local_complex_power(eps_k_w, -1.0_dp, 0.0_dp)
642 :
643 : ! b) ε^-1(iω_j,k_i)-I
644 33280 : CALL local_add_on_diag(eps_k_w, -z_one)
645 :
646 : ! c) W(iω_j,k_i) = V^0.5(k_i)*(ε^-1(iω_j,k_i)-I)*V^0.5(k_i)
647 : CALL gemm_square(V_sqrt(:, :, ikp_local), 'N', eps_k_w, 'N', &
648 33280 : V_sqrt(:, :, ikp_local), 'C', W_k_w)
649 :
650 : ! 5. W(iω,k_i) -> W^R(iω) = sum_k w_k e^(-ikR) W(iω,k) (k-point extrapolation here)
651 : CALL add_kp_to_all_rs(W_k_w, W_R, bs_env%kpoints_chi_eps_W, ikp, &
652 66664 : index_to_cell_ext=bs_env%kpoints_scf_desymm%index_to_cell)
653 :
654 : END DO ! ikp
655 :
656 104 : CALL bs_env%para_env%sync()
657 104 : CALL bs_env%para_env%sum(W_R)
658 :
659 : ! 6. W^R(iω) -> W(iω,k) [k-mesh is not extrapolated for stable mult. with M^-1(k) ]
660 : ! -> M^-1(k)*W(iω,k)*M^-1(k) =: Ŵ(iω,k) -> Ŵ^R(iω) (stored in MWM_R)
661 104 : CALL mult_W_with_Minv(W_R, MWM_R, bs_env, qs_env)
662 :
663 : ! 7. Ŵ^R(iω) -> Ŵ^R(iτ) and to fully distributed fm matrix bs_env%fm_MWM_R_t
664 824 : DO i_t = 1, bs_env%num_time_freq_points
665 704 : weight_ij = bs_env%time_frequency_grid%cosine_frequency_to_time(i_t, j_w)
666 7144 : DO i_cell_R = 1, nimages_scf_desymm
667 : CALL cp_fm_set_submatrix(bs_env%fm_MWM_R_t(i_cell_R, i_t), MWM_R(:, :, i_cell_R), &
668 7040 : alpha=weight_ij, beta=1.0_dp)
669 : END DO
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 local_complex_power(M_inv(:, :, ikp_local), -1.0_dp, 0.0_dp)
755 :
756 : ! V^0.5(k)
757 5120 : CALL local_complex_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 W_R ...
771 : !> \param MWM_R ...
772 : !> \param bs_env ...
773 : !> \param qs_env ...
774 : ! **************************************************************************************************
775 104 : SUBROUTINE mult_W_with_Minv(W_R, MWM_R, bs_env, qs_env)
776 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: W_R, MWM_R
777 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
778 : TYPE(qs_environment_type), POINTER :: qs_env
779 :
780 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mult_W_with_Minv'
781 :
782 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: M_inv, W_k, work
783 : INTEGER :: handle, ikp, n_RI
784 104 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: M_R
785 :
786 104 : CALL timeset(routineN, handle)
787 :
788 : ! compute M^R again
789 104 : CALL get_V_tr_R(M_R, bs_env%ri_metric, bs_env%regularization_RI, bs_env, qs_env)
790 :
791 104 : n_RI = bs_env%n_RI
792 832 : ALLOCATE (M_inv(n_RI, n_RI), W_k(n_RI, n_RI), work(n_RI, n_RI))
793 104 : MWM_R(:, :, :) = 0.0_dp
794 :
795 1768 : DO ikp = 1, bs_env%nkp_scf_desymm
796 :
797 : ! trivial parallelization
798 1664 : IF (MODULO(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
799 :
800 : ! M(k) = sum_R e^ikR M^R
801 : CALL rs_to_kp(M_R, M_inv, &
802 : bs_env%kpoints_scf_desymm%index_to_cell, &
803 832 : bs_env%kpoints_scf_desymm%xkp(1:3, ikp))
804 :
805 : ! invert M_PQ(k)
806 832 : CALL local_complex_power(M_inv, -1.0_dp, 0.0_dp)
807 :
808 : ! W(k) = sum_R e^ikR W^R [only R in the supercell that is determined by the SCF k-mesh]
809 : CALL rs_to_kp(W_R, W_k, &
810 : bs_env%kpoints_scf_desymm%index_to_cell, &
811 832 : bs_env%kpoints_scf_desymm%xkp(1:3, ikp))
812 :
813 : ! Ŵ(k) = M^-1(k)*W^trunc(k)*M^-1(k)
814 832 : CALL gemm_square(M_inv, 'N', W_k, 'N', M_inv, 'N', work)
815 33856 : W_k(:, :) = work(:, :)
816 :
817 : ! Ŵ^R = sum_k w_k e^(-ikR) Ŵ^(k)
818 1768 : CALL add_kp_to_all_rs(W_k, MWM_R, bs_env%kpoints_scf_desymm, ikp)
819 :
820 : END DO ! ikp
821 :
822 104 : CALL bs_env%para_env%sync()
823 104 : CALL bs_env%para_env%sum(MWM_R)
824 :
825 104 : CALL timestop(handle)
826 :
827 208 : END SUBROUTINE mult_W_with_Minv
828 :
829 : ! **************************************************************************************************
830 : !> \brief ...
831 : !> \param bs_env ...
832 : !> \param qs_env ...
833 : ! **************************************************************************************************
834 16 : SUBROUTINE compute_Sigma_x(bs_env, qs_env)
835 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
836 : TYPE(qs_environment_type), POINTER :: qs_env
837 :
838 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_x'
839 :
840 : INTEGER :: handle, i_task_Delta_R_local, ispin
841 : REAL(KIND=dp) :: t1
842 16 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: D_S, Mi_Vtr_Mi_R, Sigma_x_R
843 16 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_V
844 :
845 16 : CALL timeset(routineN, handle)
846 :
847 16 : CALL dbt_create_2c_R(Mi_Vtr_Mi_R, bs_env%t_W, bs_env%nimages_scf_desymm)
848 16 : CALL dbt_create_2c_R(D_S, bs_env%t_G, bs_env%nimages_scf_desymm)
849 16 : CALL dbt_create_2c_R(Sigma_x_R, bs_env%t_G, bs_env%nimages_scf_desymm)
850 16 : CALL dbt_create_3c_R1_R2(t_V, bs_env%t_RI_AO__AO, bs_env%nimages_3c, bs_env%nimages_3c)
851 :
852 16 : t1 = m_walltime()
853 :
854 : ! V^tr_PQ^R = <phi_P,0|V^tr|phi_Q,R>, V^tr(k) = sum_R e^ikR V^tr^R
855 : ! M(k) = sum_R e^ikR M^R, M(k) -> M^-1(k) -> Ṽ^tr(k) = M^-1(k) * V^tr(k) * M^-1(k)
856 : ! -> Ṽ^tr_PQ^R = sum_k w_k e^-ikR Ṽ^tr_PQ(k)
857 16 : CALL get_Minv_Vtr_Minv_R(Mi_Vtr_Mi_R, bs_env, qs_env)
858 :
859 : ! Σ^x_λσ^R = sum_PR1νS1 [ sum_µS2 (λ0 µS1-S2 | PR1 ) D_µν^S2 ]
860 : ! [ sum_QR2 (σR νS1 | QR1-R2) Ṽ^tr_PQ^R2 ]
861 32 : DO ispin = 1, bs_env%n_spin
862 :
863 : ! compute D^S(iτ) for cell S from D_µν(k) = sum_n^occ C^*_µn(k) C_νn(k):
864 : ! trafo k-point k -> cell S: D_µν^S = sum_k w_k D_µν(k) e^(ikS)
865 16 : CALL G_occ_vir(bs_env, 0.0_dp, D_S, ispin, occ=.TRUE., vir=.FALSE.)
866 :
867 : ! loop over ΔR = S_1 - R_1 which are local in the tensor subgroup
868 239 : DO i_task_Delta_R_local = 1, bs_env%n_tasks_Delta_R_local
869 :
870 : ! M^V_σ0,νS1,PR1 = sum_QR2 ( σ0 νS1 | QR1-R2 ) Ṽ^tr_QP^R2 for i_task_local
871 223 : CALL contract_W(t_V, Mi_Vtr_Mi_R, bs_env, i_task_Delta_R_local)
872 :
873 : ! M^D_λ0,νS1,PR1 = sum_µS2 (λ0 µS1-S2 | PR1) D_µν^S2
874 : ! Σ^x_λσ^R = sum_PR1νS1 M^D_λ0,νS1,PR1 * M^V_σR,νS1,PR1 for i_task_local, where
875 : ! M^V_σR,νS1,PR1 = M^V_σ0,νS1-R,PR1-R
876 : CALL contract_to_Sigma(Sigma_x_R, t_V, D_S, i_task_Delta_R_local, bs_env, &
877 239 : occ=.TRUE., vir=.FALSE., clear_t_W=.TRUE., fill_skip=.FALSE.)
878 :
879 : END DO ! i_cell_Delta_R_local
880 :
881 16 : CALL bs_env%para_env%sync()
882 :
883 : CALL cp_dbt_to_fm(Sigma_x_R, bs_env%fm_Sigma_x_R, bs_env%mat_ao_ao, bs_env%mat_ao_ao_tensor, &
884 32 : bs_env%para_env)
885 :
886 : END DO ! ispin
887 :
888 16 : IF (bs_env%unit_nr > 0) THEN
889 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
890 8 : 'Computed Σ^x,', ' Execution time', m_walltime() - t1, ' s'
891 8 : WRITE (bs_env%unit_nr, '(A)') ' '
892 : END IF
893 :
894 16 : CALL destroy_t_1d(Mi_Vtr_Mi_R)
895 16 : CALL destroy_t_1d(D_S)
896 16 : CALL destroy_t_1d(Sigma_x_R)
897 16 : CALL destroy_t_2d(t_V)
898 :
899 16 : CALL timestop(handle)
900 :
901 16 : END SUBROUTINE compute_Sigma_x
902 :
903 : ! **************************************************************************************************
904 : !> \brief ...
905 : !> \param bs_env ...
906 : ! **************************************************************************************************
907 16 : SUBROUTINE compute_Sigma_c(bs_env)
908 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
909 :
910 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c'
911 :
912 : INTEGER :: handle, i_t, i_task_Delta_R_local, ispin
913 : REAL(KIND=dp) :: t1, tau
914 16 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: Gocc_S, Gvir_S, Sigma_c_R_neg_tau, &
915 16 : Sigma_c_R_pos_tau, W_R
916 16 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_W
917 :
918 16 : CALL timeset(routineN, handle)
919 :
920 16 : CALL dbt_create_2c_R(Gocc_S, bs_env%t_G, bs_env%nimages_scf_desymm)
921 16 : CALL dbt_create_2c_R(Gvir_S, bs_env%t_G, bs_env%nimages_scf_desymm)
922 16 : CALL dbt_create_2c_R(W_R, bs_env%t_W, bs_env%nimages_scf_desymm)
923 16 : CALL dbt_create_3c_R1_R2(t_W, bs_env%t_RI_AO__AO, bs_env%nimages_3c, bs_env%nimages_3c)
924 16 : CALL dbt_create_2c_R(Sigma_c_R_neg_tau, bs_env%t_G, bs_env%nimages_scf_desymm)
925 16 : CALL dbt_create_2c_R(Sigma_c_R_pos_tau, bs_env%t_G, bs_env%nimages_scf_desymm)
926 :
927 : ! Σ^c_λσ^R(iτ) = sum_PR1νS1 [ sum_µS2 (λ0 µS1-S2 | PR1 ) G^occ/vir_µν^S2(i|τ|) ]
928 : ! [ sum_QR2 (σR νS1 | QR1-R2) Ŵ_PQ^R2(iτ) ]
929 120 : DO i_t = 1, bs_env%num_time_freq_points
930 :
931 224 : DO ispin = 1, bs_env%n_spin
932 :
933 104 : t1 = m_walltime()
934 :
935 104 : tau = bs_env%time_frequency_grid%imaginary_time(i_t)
936 :
937 : ! G^occ_µλ(i|τ|,k) = sum_n^occ C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k), τ < 0
938 : ! G^vir_µλ(i|τ|,k) = sum_n^vir C_µn(k)^* e^(-|(ϵ_nk-ϵ_F)τ|) C_λn(k), τ > 0
939 : ! k-point k -> cell S: G^occ/vir_µλ^S(i|τ|) = sum_k w_k G^occ/vir_µλ(i|τ|,k) e^(ikS)
940 104 : CALL G_occ_vir(bs_env, tau, Gocc_S, ispin, occ=.TRUE., vir=.FALSE.)
941 104 : CALL G_occ_vir(bs_env, tau, Gvir_S, ispin, occ=.FALSE., vir=.TRUE.)
942 :
943 : ! write data of W^R_PQ(iτ) to W_R 2-index tensor
944 104 : CALL fm_MWM_R_t_to_local_tensor_W_R(bs_env%fm_MWM_R_t(:, i_t), W_R, bs_env)
945 :
946 : ! loop over ΔR = S_1 - R_1 which are local in the tensor subgroup
947 1494 : DO i_task_Delta_R_local = 1, bs_env%n_tasks_Delta_R_local
948 :
949 1390 : IF (bs_env%skip_DR_Sigma(i_task_Delta_R_local)) CYCLE
950 :
951 : ! for i_task_local (i.e. fixed ΔR = S_1 - R_1) and for all τ (W(iτ) = W(-iτ)):
952 : ! M^W_σ0,νS1,PR1 = sum_QR2 ( σ0 νS1 | QR1-R2 ) W(iτ)_QP^R2
953 510 : CALL contract_W(t_W, W_R, bs_env, i_task_Delta_R_local)
954 :
955 : ! for τ < 0 and for i_task_local (i.e. fixed ΔR = S_1 - R_1):
956 : ! M^G_λ0,νS1,PR1 = sum_µS2 (λ0 µS1-S2 | PR1) G^occ(i|τ|)_µν^S2
957 : ! Σ^c_λσ^R(iτ) = sum_PR1νS1 M^G_λ0,νS1,PR1 * M^W_σR,νS1,PR1
958 : ! where M^W_σR,νS1,PR1 = M^W_σ0,νS1-R,PR1-R
959 : CALL contract_to_Sigma(Sigma_c_R_neg_tau, t_W, Gocc_S, i_task_Delta_R_local, bs_env, &
960 510 : occ=.TRUE., vir=.FALSE., clear_t_W=.FALSE., fill_skip=.FALSE.)
961 :
962 : ! for τ > 0: same as for τ < 0, but G^occ -> G^vir
963 : CALL contract_to_Sigma(Sigma_c_R_pos_tau, t_W, Gvir_S, i_task_Delta_R_local, bs_env, &
964 1494 : occ=.FALSE., vir=.TRUE., clear_t_W=.TRUE., fill_skip=.TRUE.)
965 :
966 : END DO ! i_cell_Delta_R_local
967 :
968 104 : CALL bs_env%para_env%sync()
969 :
970 : CALL cp_dbt_to_fm(Sigma_c_R_pos_tau, bs_env%fm_Sigma_c_R_pos_tau(:, i_t, ispin), bs_env%mat_ao_ao, &
971 104 : bs_env%mat_ao_ao_tensor, bs_env%para_env)
972 :
973 : CALL cp_dbt_to_fm(Sigma_c_R_neg_tau, bs_env%fm_Sigma_c_R_neg_tau(:, i_t, ispin), bs_env%mat_ao_ao, &
974 104 : bs_env%mat_ao_ao_tensor, bs_env%para_env)
975 :
976 208 : IF (bs_env%unit_nr > 0) THEN
977 : WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F7.1,A)') &
978 52 : 'Computed Σ^c(iτ) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
979 104 : ', Execution time', m_walltime() - t1, ' s'
980 : END IF
981 :
982 : END DO ! ispin
983 :
984 : END DO ! i_t
985 :
986 16 : CALL destroy_t_1d(Gocc_S)
987 16 : CALL destroy_t_1d(Gvir_S)
988 16 : CALL destroy_t_1d(W_R)
989 16 : CALL destroy_t_1d(Sigma_c_R_neg_tau)
990 16 : CALL destroy_t_1d(Sigma_c_R_pos_tau)
991 16 : CALL destroy_t_2d(t_W)
992 :
993 16 : CALL timestop(handle)
994 :
995 16 : END SUBROUTINE compute_Sigma_c
996 :
997 : ! **************************************************************************************************
998 : !> \brief ...
999 : !> \param Mi_Vtr_Mi_R ...
1000 : !> \param bs_env ...
1001 : !> \param qs_env ...
1002 : ! **************************************************************************************************
1003 16 : SUBROUTINE get_Minv_Vtr_Minv_R(Mi_Vtr_Mi_R, bs_env, qs_env)
1004 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: Mi_Vtr_Mi_R
1005 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1006 : TYPE(qs_environment_type), POINTER :: qs_env
1007 :
1008 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_Minv_Vtr_Minv_R'
1009 :
1010 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: M_kp, Mi_Vtr_Mi_kp, V_tr_kp
1011 : INTEGER :: handle, i_cell_R, ikp, n_RI, &
1012 : nimages_scf, nkp_scf
1013 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: M_R, Mi_Vtr_Mi_R_arr, V_tr_R
1014 :
1015 16 : CALL timeset(routineN, handle)
1016 :
1017 16 : nimages_scf = bs_env%nimages_scf_desymm
1018 16 : nkp_scf = bs_env%kpoints_scf_desymm%nkp
1019 16 : n_RI = bs_env%n_RI
1020 :
1021 16 : CALL get_V_tr_R(V_tr_R, bs_env%trunc_coulomb, 0.0_dp, bs_env, qs_env)
1022 16 : CALL get_V_tr_R(M_R, bs_env%ri_metric, bs_env%regularization_RI, bs_env, qs_env)
1023 :
1024 : ALLOCATE (V_tr_kp(n_RI, n_RI), M_kp(n_RI, n_RI), &
1025 192 : Mi_Vtr_Mi_kp(n_RI, n_RI), Mi_Vtr_Mi_R_arr(n_RI, n_RI, nimages_scf))
1026 16 : Mi_Vtr_Mi_R_arr(:, :, :) = 0.0_dp
1027 :
1028 272 : DO ikp = 1, nkp_scf
1029 : ! trivial parallelization
1030 256 : IF (MODULO(ikp, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
1031 : ! V_tr(k) = sum_R e^ikR V_tr^R
1032 : CALL rs_to_kp(V_tr_R, V_tr_kp, bs_env%kpoints_scf_desymm%index_to_cell, &
1033 128 : bs_env%kpoints_scf_desymm%xkp(1:3, ikp))
1034 : ! M(k) = sum_R e^ikR M^R
1035 : CALL rs_to_kp(M_R, M_kp, bs_env%kpoints_scf_desymm%index_to_cell, &
1036 128 : bs_env%kpoints_scf_desymm%xkp(1:3, ikp))
1037 : ! M(k) -> M^-1(k)
1038 128 : CALL local_complex_power(M_kp, -1.0_dp, 0.0_dp)
1039 : ! Ṽ(k) = M^-1(k) * V_tr(k) * M^-1(k)
1040 128 : CALL gemm_square(M_kp, 'N', V_tr_kp, 'N', M_kp, 'N', Mi_Vtr_Mi_kp)
1041 : ! Ṽ^R = sum_k w_k e^-ikR Ṽ(k)
1042 272 : CALL add_kp_to_all_rs(Mi_Vtr_Mi_kp, Mi_Vtr_Mi_R_arr, bs_env%kpoints_scf_desymm, ikp)
1043 : END DO
1044 16 : CALL bs_env%para_env%sync()
1045 16 : CALL bs_env%para_env%sum(Mi_Vtr_Mi_R_arr)
1046 :
1047 : ! use bs_env%fm_chi_R_t for temporary storage
1048 160 : DO i_cell_R = 1, nimages_scf
1049 160 : CALL cp_fm_set_submatrix(bs_env%fm_chi_R_t(i_cell_R, 1), Mi_Vtr_Mi_R_arr(:, :, i_cell_R))
1050 : END DO
1051 :
1052 : ! communicate Mi_Vtr_Mi_R to tensor format; full replication in tensor group
1053 160 : DO i_cell_R = 1, nimages_scf
1054 : CALL cp_fm_to_dbt(bs_env%fm_chi_R_t(i_cell_R, 1), bs_env%mat_RI_RI%matrix, &
1055 : bs_env%mat_RI_RI_tensor%matrix, Mi_Vtr_Mi_R(i_cell_R), bs_env%para_env, &
1056 160 : bs_env%para_env_tensor%num_pe, bs_env%eps_filter)
1057 : END DO
1058 :
1059 16 : CALL timestop(handle)
1060 :
1061 32 : END SUBROUTINE get_Minv_Vtr_Minv_R
1062 :
1063 : ! **************************************************************************************************
1064 : !> \brief ...
1065 : !> \param t_1d ...
1066 : ! **************************************************************************************************
1067 440 : SUBROUTINE destroy_t_1d(t_1d)
1068 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_1d
1069 :
1070 : CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_t_1d'
1071 :
1072 : INTEGER :: handle, i
1073 :
1074 440 : CALL timeset(routineN, handle)
1075 :
1076 4400 : DO i = 1, SIZE(t_1d)
1077 4400 : CALL dbt_destroy(t_1d(i))
1078 : END DO
1079 4400 : DEALLOCATE (t_1d)
1080 :
1081 440 : CALL timestop(handle)
1082 :
1083 440 : END SUBROUTINE destroy_t_1d
1084 :
1085 : ! **************************************************************************************************
1086 : !> \brief ...
1087 : !> \param t_2d ...
1088 : ! **************************************************************************************************
1089 240 : SUBROUTINE destroy_t_2d(t_2d)
1090 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_2d
1091 :
1092 : CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_t_2d'
1093 :
1094 : INTEGER :: handle, i, j
1095 :
1096 240 : CALL timeset(routineN, handle)
1097 :
1098 2812 : DO i = 1, SIZE(t_2d, 1)
1099 32136 : DO j = 1, SIZE(t_2d, 2)
1100 31896 : CALL dbt_destroy(t_2d(i, j))
1101 : END DO
1102 : END DO
1103 29564 : DEALLOCATE (t_2d)
1104 :
1105 240 : CALL timestop(handle)
1106 :
1107 240 : END SUBROUTINE destroy_t_2d
1108 :
1109 : ! **************************************************************************************************
1110 : !> \brief ...
1111 : !> \param t_W ...
1112 : !> \param W_R ...
1113 : !> \param bs_env ...
1114 : !> \param i_task_Delta_R_local ...
1115 : ! **************************************************************************************************
1116 733 : SUBROUTINE contract_W(t_W, W_R, bs_env, i_task_Delta_R_local)
1117 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_W
1118 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: W_R
1119 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1120 : INTEGER :: i_task_Delta_R_local
1121 :
1122 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_W'
1123 :
1124 : INTEGER :: handle, i_cell_Delta_R, i_cell_R1, &
1125 : i_cell_R2, i_cell_R2_m_R1, i_cell_S1, &
1126 : i_cell_S1_m_R1_p_R2
1127 : INTEGER, DIMENSION(3) :: cell_DR, cell_R1, cell_R2, cell_R2_m_R1, &
1128 : cell_S1, cell_S1_m_R2_p_R1
1129 : LOGICAL :: cell_found
1130 12461 : TYPE(dbt_type) :: t_3c_int, t_W_tmp
1131 :
1132 733 : CALL timeset(routineN, handle)
1133 :
1134 733 : CALL dbt_create(bs_env%t_RI__AO_AO, t_W_tmp)
1135 733 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_int)
1136 :
1137 733 : i_cell_Delta_R = bs_env%task_Delta_R(i_task_Delta_R_local)
1138 :
1139 9200 : DO i_cell_R1 = 1, bs_env%nimages_3c
1140 :
1141 33868 : cell_R1(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_R1)
1142 33868 : cell_DR(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_Delta_R)
1143 :
1144 : ! S_1 = R_1 + ΔR (from ΔR = S_1 - R_1)
1145 : CALL add_R(cell_R1, cell_DR, bs_env%index_to_cell_3c, cell_S1, &
1146 8467 : cell_found, bs_env%cell_to_index_3c, i_cell_S1)
1147 8467 : IF (.NOT. cell_found) CYCLE
1148 :
1149 50570 : DO i_cell_R2 = 1, bs_env%nimages_scf_desymm
1150 :
1151 148932 : cell_R2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_R2)
1152 :
1153 : ! R_2 - R_1
1154 : CALL add_R(cell_R2, -cell_R1, bs_env%index_to_cell_3c, cell_R2_m_R1, &
1155 148932 : cell_found, bs_env%cell_to_index_3c, i_cell_R2_m_R1)
1156 37233 : IF (.NOT. cell_found) CYCLE
1157 :
1158 : ! S_1 - R_1 + R_2
1159 : CALL add_R(cell_S1, cell_R2_m_R1, bs_env%index_to_cell_3c, cell_S1_m_R2_p_R1, &
1160 22627 : cell_found, bs_env%cell_to_index_3c, i_cell_S1_m_R1_p_R2)
1161 22627 : IF (.NOT. cell_found) CYCLE
1162 :
1163 17285 : CALL get_t_3c_int(t_3c_int, bs_env, i_cell_S1_m_R1_p_R2, i_cell_R2_m_R1)
1164 :
1165 : ! M^W_σ0,νS1,PR1 = sum_QR2 ( σ0 νS1 | QR1-R2 ) W_QP^R2
1166 : ! = sum_QR2 ( σR2-R1 νS1-R1+R2 | Q0 ) W_QP^R2
1167 : ! for ΔR = S_1 - R_1
1168 : CALL dbt_contract(alpha=1.0_dp, &
1169 : tensor_1=W_R(i_cell_R2), &
1170 : tensor_2=t_3c_int, &
1171 : beta=0.0_dp, &
1172 : tensor_3=t_W_tmp, &
1173 : contract_1=[1], notcontract_1=[2], map_1=[1], &
1174 : contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3], &
1175 17285 : filter_eps=bs_env%eps_filter)
1176 :
1177 : ! reorder tensor
1178 : CALL dbt_copy(t_W_tmp, t_W(i_cell_S1, i_cell_R1), order=[1, 2, 3], &
1179 62985 : move_data=.TRUE., summation=.TRUE.)
1180 :
1181 : END DO ! i_cell_R2
1182 :
1183 : END DO ! i_cell_R1
1184 :
1185 733 : CALL dbt_destroy(t_W_tmp)
1186 733 : CALL dbt_destroy(t_3c_int)
1187 :
1188 733 : CALL timestop(handle)
1189 :
1190 733 : END SUBROUTINE contract_W
1191 :
1192 : ! **************************************************************************************************
1193 : !> \brief ...
1194 : !> \param Sigma_R ...
1195 : !> \param t_W ...
1196 : !> \param G_S ...
1197 : !> \param i_task_Delta_R_local ...
1198 : !> \param bs_env ...
1199 : !> \param occ ...
1200 : !> \param vir ...
1201 : !> \param clear_t_W ...
1202 : !> \param fill_skip ...
1203 : ! **************************************************************************************************
1204 1243 : SUBROUTINE contract_to_Sigma(Sigma_R, t_W, G_S, i_task_Delta_R_local, bs_env, occ, vir, &
1205 : clear_t_W, fill_skip)
1206 : TYPE(dbt_type), DIMENSION(:) :: Sigma_R
1207 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_W
1208 : TYPE(dbt_type), DIMENSION(:) :: G_S
1209 : INTEGER :: i_task_Delta_R_local
1210 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1211 : LOGICAL :: occ, vir, clear_t_W, fill_skip
1212 :
1213 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_to_Sigma'
1214 :
1215 : INTEGER :: handle, handle2, i_cell_Delta_R, i_cell_m_R1, i_cell_R, i_cell_R1, &
1216 : i_cell_R1_minus_R, i_cell_S1, i_cell_S1_minus_R, i_cell_S1_p_S2_m_R1, i_cell_S2
1217 : INTEGER(KIND=int_8) :: flop, flop_tmp
1218 : INTEGER, DIMENSION(3) :: cell_DR, cell_m_R1, cell_R, cell_R1, &
1219 : cell_R1_minus_R, cell_S1, &
1220 : cell_S1_minus_R, cell_S1_p_S2_m_R1, &
1221 : cell_S2
1222 : LOGICAL :: cell_found
1223 : REAL(KIND=dp) :: sign_Sigma
1224 31075 : TYPE(dbt_type) :: t_3c_int, t_G, t_G_2
1225 :
1226 1243 : CALL timeset(routineN, handle)
1227 :
1228 1243 : CPASSERT(occ .EQV. (.NOT. vir))
1229 1243 : IF (occ) sign_Sigma = -1.0_dp
1230 1243 : IF (vir) sign_Sigma = 1.0_dp
1231 :
1232 1243 : CALL dbt_create(bs_env%t_RI_AO__AO, t_G)
1233 1243 : CALL dbt_create(bs_env%t_RI_AO__AO, t_G_2)
1234 1243 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_int)
1235 :
1236 1243 : i_cell_Delta_R = bs_env%task_Delta_R(i_task_Delta_R_local)
1237 :
1238 1243 : flop = 0_int_8
1239 :
1240 15592 : DO i_cell_R1 = 1, bs_env%nimages_3c
1241 :
1242 57396 : cell_R1(1:3) = bs_env%index_to_cell_3c(1:3, i_cell_R1)
1243 57396 : cell_DR(1:3) = bs_env%index_to_cell_Delta_R(1:3, i_cell_Delta_R)
1244 :
1245 : ! S_1 = R_1 + ΔR (from ΔR = S_1 - R_1)
1246 : CALL add_R(cell_R1, cell_DR, bs_env%index_to_cell_3c, cell_S1, cell_found, &
1247 14349 : bs_env%cell_to_index_3c, i_cell_S1)
1248 14349 : IF (.NOT. cell_found) CYCLE
1249 :
1250 73040 : DO i_cell_S2 = 1, bs_env%nimages_scf_desymm
1251 :
1252 65736 : IF (bs_env%skip_DR_R1_S2_Gx3c_Sigma(i_task_Delta_R_local, i_cell_R1, i_cell_S2)) CYCLE
1253 :
1254 207144 : cell_S2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_S2)
1255 207144 : cell_m_R1(1:3) = -cell_R1(1:3)
1256 207144 : cell_S1_p_S2_m_R1(1:3) = cell_S1(1:3) + cell_S2(1:3) - cell_R1(1:3)
1257 :
1258 51786 : CALL is_cell_in_index_to_cell(cell_m_R1, bs_env%index_to_cell_3c, cell_found)
1259 51786 : IF (.NOT. cell_found) CYCLE
1260 :
1261 38997 : CALL is_cell_in_index_to_cell(cell_S1_p_S2_m_R1, bs_env%index_to_cell_3c, cell_found)
1262 38997 : IF (.NOT. cell_found) CYCLE
1263 :
1264 21513 : i_cell_m_R1 = bs_env%cell_to_index_3c(cell_m_R1(1), cell_m_R1(2), cell_m_R1(3))
1265 : i_cell_S1_p_S2_m_R1 = bs_env%cell_to_index_3c(cell_S1_p_S2_m_R1(1), &
1266 : cell_S1_p_S2_m_R1(2), &
1267 21513 : cell_S1_p_S2_m_R1(3))
1268 :
1269 21513 : CALL timeset(routineN//"_3c_x_G", handle2)
1270 :
1271 21513 : CALL get_t_3c_int(t_3c_int, bs_env, i_cell_m_R1, i_cell_S1_p_S2_m_R1)
1272 :
1273 : ! M_λ0,νS1,PR1 = sum_µS2 ( λ0 µS1-S2 | PR1 ) G^occ/vir_µν^S2(i|τ|)
1274 : ! = sum_µS2 ( λ-R1 µS1-S2-R1 | P0 ) G^occ/vir_µν^S2(i|τ|)
1275 : ! for ΔR = S_1 - R_1
1276 : CALL dbt_contract(alpha=1.0_dp, &
1277 : tensor_1=G_S(i_cell_S2), &
1278 : tensor_2=t_3c_int, &
1279 : beta=1.0_dp, &
1280 : tensor_3=t_G, &
1281 : contract_1=[2], notcontract_1=[1], map_1=[3], &
1282 : contract_2=[3], notcontract_2=[1, 2], map_2=[1, 2], &
1283 21513 : filter_eps=bs_env%eps_filter, flop=flop_tmp)
1284 :
1285 21513 : IF (flop_tmp == 0_int_8 .AND. fill_skip) THEN
1286 2561 : bs_env%skip_DR_R1_S2_Gx3c_Sigma(i_task_Delta_R_local, i_cell_R1, i_cell_S2) = .TRUE.
1287 : END IF
1288 :
1289 94553 : CALL timestop(handle2)
1290 :
1291 : END DO ! i_cell_S2
1292 :
1293 7304 : CALL dbt_copy(t_G, t_G_2, order=[1, 3, 2], move_data=.TRUE.)
1294 :
1295 7304 : CALL timeset(routineN//"_contract", handle2)
1296 :
1297 73040 : DO i_cell_R = 1, bs_env%nimages_scf_desymm
1298 :
1299 65736 : IF (bs_env%skip_DR_R1_R_MxM_Sigma(i_task_Delta_R_local, i_cell_R1, i_cell_R)) CYCLE
1300 :
1301 189192 : cell_R = bs_env%kpoints_scf_desymm%index_to_cell(1:3, i_cell_R)
1302 :
1303 : ! R_1 - R
1304 : CALL add_R(cell_R1, -cell_R, bs_env%index_to_cell_3c, cell_R1_minus_R, &
1305 189192 : cell_found, bs_env%cell_to_index_3c, i_cell_R1_minus_R)
1306 47298 : IF (.NOT. cell_found) CYCLE
1307 :
1308 : ! S_1 - R
1309 : CALL add_R(cell_S1, -cell_R, bs_env%index_to_cell_3c, cell_S1_minus_R, &
1310 103036 : cell_found, bs_env%cell_to_index_3c, i_cell_S1_minus_R)
1311 25759 : IF (.NOT. cell_found) CYCLE
1312 :
1313 : ! Σ_λσ^R = sum_PR1νS1 M^G_λ0,νS1,PR1 M^W_σR,νS1,PR1, where
1314 : ! M^G_λ0,νS1,PR1 = sum_µS2 (λ0 µS1-S2 | PR1) G_µν^S2
1315 : ! M^W_σR,νS1,PR1 = sum_QR2 (σR νS1 | QR1-R2) W_PQ^R2 = M^W_σ0,νS1-R,PR1-R
1316 : CALL dbt_contract(alpha=sign_Sigma, &
1317 : tensor_1=t_G_2, &
1318 : tensor_2=t_W(i_cell_S1_minus_R, i_cell_R1_minus_R), &
1319 : beta=1.0_dp, &
1320 : tensor_3=Sigma_R(i_cell_R), &
1321 : contract_1=[1, 2], notcontract_1=[3], map_1=[1], &
1322 : contract_2=[1, 2], notcontract_2=[3], map_2=[2], &
1323 16108 : filter_eps=bs_env%eps_filter, flop=flop_tmp)
1324 :
1325 16108 : flop = flop + flop_tmp
1326 :
1327 23412 : IF (flop_tmp == 0_int_8 .AND. fill_skip) THEN
1328 3900 : bs_env%skip_DR_R1_R_MxM_Sigma(i_task_Delta_R_local, i_cell_R1, i_cell_R) = .TRUE.
1329 : END IF
1330 :
1331 : END DO ! i_cell_R
1332 :
1333 7304 : CALL dbt_clear(t_G_2)
1334 :
1335 30200 : CALL timestop(handle2)
1336 :
1337 : END DO ! i_cell_R1
1338 :
1339 1243 : IF (vir .AND. flop == 0_int_8) bs_env%skip_DR_Sigma(i_task_Delta_R_local) = .TRUE.
1340 :
1341 : ! release memory
1342 1243 : IF (clear_t_W) THEN
1343 9200 : DO i_cell_S1 = 1, bs_env%nimages_3c
1344 109159 : DO i_cell_R1 = 1, bs_env%nimages_3c
1345 108426 : CALL dbt_clear(t_W(i_cell_S1, i_cell_R1))
1346 : END DO
1347 : END DO
1348 : END IF
1349 :
1350 1243 : CALL dbt_destroy(t_G)
1351 1243 : CALL dbt_destroy(t_G_2)
1352 1243 : CALL dbt_destroy(t_3c_int)
1353 :
1354 1243 : CALL timestop(handle)
1355 :
1356 1243 : END SUBROUTINE contract_to_Sigma
1357 :
1358 : ! **************************************************************************************************
1359 : !> \brief ...
1360 : !> \param fm_W_R ...
1361 : !> \param W_R ...
1362 : !> \param bs_env ...
1363 : ! **************************************************************************************************
1364 104 : SUBROUTINE fm_MWM_R_t_to_local_tensor_W_R(fm_W_R, W_R, bs_env)
1365 : TYPE(cp_fm_type), DIMENSION(:) :: fm_W_R
1366 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: W_R
1367 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1368 :
1369 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_MWM_R_t_to_local_tensor_W_R'
1370 :
1371 : INTEGER :: handle, i_cell_R
1372 :
1373 104 : CALL timeset(routineN, handle)
1374 :
1375 : ! communicate fm_W_R to tensor W_R; full replication in tensor group
1376 1040 : DO i_cell_R = 1, bs_env%nimages_scf_desymm
1377 : CALL cp_fm_to_dbt(fm_W_R(i_cell_R), bs_env%mat_RI_RI%matrix, bs_env%mat_RI_RI_tensor%matrix, &
1378 : W_R(i_cell_R), bs_env%para_env, bs_env%para_env_tensor%num_pe, &
1379 1040 : bs_env%eps_filter)
1380 : END DO
1381 :
1382 104 : CALL timestop(handle)
1383 :
1384 104 : END SUBROUTINE fm_MWM_R_t_to_local_tensor_W_R
1385 :
1386 : ! **************************************************************************************************
1387 : !> \brief ...
1388 : !> \param bs_env ...
1389 : ! **************************************************************************************************
1390 16 : SUBROUTINE compute_QP_energies(bs_env)
1391 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1392 :
1393 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies'
1394 :
1395 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Sigma_c_ikp_n_freq
1396 : INTEGER :: handle, ikp, ispin, j_t
1397 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Sigma_x_ikp_n
1398 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Sigma_c_ikp_n_time
1399 : TYPE(cp_cfm_type) :: cfm_mo_coeff
1400 :
1401 16 : CALL timeset(routineN, handle)
1402 :
1403 16 : CALL cp_cfm_create(cfm_mo_coeff, bs_env%fm_s_Gamma%matrix_struct)
1404 48 : ALLOCATE (Sigma_x_ikp_n(bs_env%n_ao))
1405 80 : ALLOCATE (Sigma_c_ikp_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
1406 64 : ALLOCATE (Sigma_c_ikp_n_freq(bs_env%n_ao, bs_env%num_freq_points_fit))
1407 :
1408 32 : DO ispin = 1, bs_env%n_spin
1409 :
1410 350 : DO ikp = 1, bs_env%nkp_bs_and_DOS
1411 :
1412 : ! 1. get C_µn(k)
1413 318 : CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mo_coeff)
1414 :
1415 : ! 2. Σ^x_µν(k) = sum_R Σ^x_µν^R e^ikR
1416 : ! Σ^x_nn(k) = sum_µν C^*_µn(k) Σ^x_µν(k) C_νn(k)
1417 318 : CALL trafo_to_k_and_nn(bs_env%fm_Sigma_x_R, Sigma_x_ikp_n, cfm_mo_coeff, bs_env, ikp)
1418 :
1419 : ! 3. Σ^c_µν(k,+/-i|τ_j|) = sum_R Σ^c_µν^R(+/-i|τ_j|) e^ikR
1420 : ! Σ^c_nn(k,+/-i|τ_j|) = sum_µν C^*_µn(k) Σ^c_µν(k,+/-i|τ_j|) C_νn(k)
1421 2602 : DO j_t = 1, bs_env%num_time_freq_points
1422 : CALL trafo_to_k_and_nn(bs_env%fm_Sigma_c_R_pos_tau(:, j_t, ispin), &
1423 2284 : Sigma_c_ikp_n_time(:, j_t, positive_time), cfm_mo_coeff, bs_env, ikp)
1424 : CALL trafo_to_k_and_nn(bs_env%fm_Sigma_c_R_neg_tau(:, j_t, ispin), &
1425 2602 : Sigma_c_ikp_n_time(:, j_t, negative_time), cfm_mo_coeff, bs_env, ikp)
1426 : END DO
1427 :
1428 : ! 4. Σ^c_nn(k_i,iω) = ∫ from -∞ to ∞ dτ e^-iωτ Σ^c_nn(k_i,iτ)
1429 318 : CALL transform_time_to_frequency(bs_env%time_frequency_grid, Sigma_c_ikp_n_time, Sigma_c_ikp_n_freq)
1430 6358 : Sigma_c_ikp_n_freq(1:bs_env%n_occ(ispin), :) = CONJG(Sigma_c_ikp_n_freq(1:bs_env%n_occ(ispin), :))
1431 :
1432 : ! 5. Analytic continuation Σ^c_nn(k_i,iω) -> Σ^c_nn(k_i,ϵ) and
1433 : ! ϵ_nk_i^GW = ϵ_nk_i^DFT + Σ^c_nn(k_i,ϵ) + Σ^x_nn(k_i) - v^xc_nn(k_i)
1434 : CALL analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, &
1435 : bs_env%v_xc_n(:, ikp, ispin), &
1436 334 : bs_env%eigenval_scf(:, ikp, ispin), ikp, ispin)
1437 :
1438 : END DO ! ikp
1439 :
1440 : END DO ! ispin
1441 :
1442 16 : CALL get_all_VBM_CBM_bandgaps(bs_env)
1443 :
1444 16 : CALL cp_cfm_release(cfm_mo_coeff)
1445 :
1446 16 : CALL timestop(handle)
1447 :
1448 32 : END SUBROUTINE compute_QP_energies
1449 :
1450 : ! **************************************************************************************************
1451 : !> \brief ...
1452 : !> \param fm_rs ...
1453 : !> \param array_ikp_n ...
1454 : !> \param cfm_mo_coeff ...
1455 : !> \param bs_env ...
1456 : !> \param ikp ...
1457 : ! **************************************************************************************************
1458 4886 : SUBROUTINE trafo_to_k_and_nn(fm_rs, array_ikp_n, cfm_mo_coeff, bs_env, ikp)
1459 : TYPE(cp_fm_type), DIMENSION(:) :: fm_rs
1460 : REAL(KIND=dp), DIMENSION(:) :: array_ikp_n
1461 : TYPE(cp_cfm_type) :: cfm_mo_coeff
1462 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1463 : INTEGER :: ikp
1464 :
1465 : CHARACTER(LEN=*), PARAMETER :: routineN = 'trafo_to_k_and_nn'
1466 :
1467 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cdiag
1468 : INTEGER :: handle
1469 : TYPE(cp_cfm_type) :: cfm_ikp
1470 :
1471 4886 : CALL timeset(routineN, handle)
1472 :
1473 4886 : CALL cp_cfm_create(cfm_ikp, cfm_mo_coeff%matrix_struct)
1474 :
1475 : ! Σ_µν(k_i) = sum_R e^ik_iR Σ_µν^R
1476 4886 : CALL fm_rs_to_kp(cfm_ikp, fm_rs, bs_env%kpoints_DOS, ikp)
1477 :
1478 : ! Σ_nm(k_i) = sum_µν C^*_µn(k_i) Σ_µν(k_i) C_νn(k_i)
1479 4886 : CALL cfm_contract_ABA(cfm_mo_coeff, cfm_ikp)
1480 :
1481 : ! get Σ_nn(k_i) which is a real quantity as Σ^x and Σ^c(iτ) is Hermitian
1482 14658 : ALLOCATE (cdiag(SIZE(array_ikp_n)))
1483 4886 : CALL cp_cfm_get_diag(cfm_ikp, cdiag)
1484 63420 : array_ikp_n = REAL(cdiag, KIND=dp)
1485 :
1486 4886 : DEALLOCATE (cdiag)
1487 4886 : CALL cp_cfm_release(cfm_ikp)
1488 :
1489 4886 : CALL timestop(handle)
1490 :
1491 4886 : END SUBROUTINE trafo_to_k_and_nn
1492 :
1493 : END MODULE gw_tensor_small_cell_full_kp
|
