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