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