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