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