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