Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines from paper [Graml2024]
10 : !> \par History
11 : !> 01.2026 Maximilian Graml: add more bounds to exploit sparsity in 3c integrals, fixes
12 : !> \author Jan Wilhelm
13 : !> \date 07.2023
14 : ! **************************************************************************************************
15 : MODULE gw_tensor_large_cell_Gamma
16 : USE bibliography, ONLY: Graml2024,&
17 : cite_reference
18 : USE cell_types, ONLY: cell_type,&
19 : pbc
20 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_get_diag
21 : USE cp_cfm_diag, ONLY: cp_cfm_geeig
22 : USE cp_cfm_types, ONLY: cp_cfm_create,&
23 : cp_cfm_release,&
24 : cp_cfm_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
26 : dbcsr_p_type,&
27 : dbcsr_release,&
28 : dbcsr_type
29 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
30 : copy_fm_to_dbcsr,&
31 : dbcsr_deallocate_matrix_set
32 : USE cp_fm_types, ONLY: cp_fm_create,&
33 : cp_fm_release,&
34 : cp_fm_type
35 : USE dbt_api, ONLY: dbt_clear,&
36 : dbt_contract,&
37 : dbt_copy,&
38 : dbt_create,&
39 : dbt_destroy,&
40 : dbt_filter,&
41 : dbt_type
42 : USE dbt_matrix_redistribution, ONLY: cp_dbt_to_dbcsr,&
43 : cp_fm_to_dbt
44 : USE gw_compute_V_and_W, ONLY: fm_read,&
45 : fm_write,&
46 : get_W_MIC
47 : USE gw_utils, ONLY: G_occ_vir,&
48 : analyt_conti_and_print,&
49 : de_init_bs_env
50 : USE gw_utils_fm, ONLY: cfm_contract_ABA,&
51 : fm_contract_ABA
52 : USE input_constants, ONLY: G0W0
53 : USE kinds, ONLY: default_path_length,&
54 : dp,&
55 : int_8
56 : USE machine, ONLY: m_walltime
57 : USE message_passing, ONLY: mp_file_delete
58 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat
59 : USE particle_types, ONLY: particle_type
60 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
61 : USE post_scf_bandstructure_utils, ONLY: cfm_ikp_from_fm_Gamma,&
62 : get_all_VBM_CBM_bandgaps
63 : USE qs_environment_types, ONLY: get_qs_env,&
64 : qs_environment_type
65 : USE qs_tensors, ONLY: build_3c_integrals
66 : USE time_frequency_grids, ONLY: negative_time,&
67 : positive_time,&
68 : transform_time_to_frequency
69 : #include "./base/base_uses.f90"
70 :
71 : IMPLICIT NONE
72 :
73 : PRIVATE
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_tensor_large_cell_Gamma'
76 :
77 : PUBLIC :: gw_calc_tensor_large_cell_Gamma, &
78 : compute_3c_integrals, &
79 : G_occ_vir, &
80 : delete_unnecessary_files, &
81 : compute_QP_energies, &
82 : get_mat_chi_Gamma_tau, &
83 : get_Sigma_x, &
84 : create_mat_for_Sigma_c, &
85 : destroy_mat_Sigma_c, &
86 : compute_Sigma_c_occ_vir, &
87 : print_skipping
88 :
89 : CONTAINS
90 :
91 : ! **************************************************************************************************
92 : !> \brief Perform GW band structure calculation
93 : !> \param qs_env ...
94 : !> \param bs_env Band-structure environment containing GW parameters.
95 : !> \par History
96 : !> * 07.2023 created [Jan Wilhelm]
97 : ! **************************************************************************************************
98 74 : SUBROUTINE gw_calc_tensor_large_cell_Gamma(qs_env, bs_env)
99 : TYPE(qs_environment_type), POINTER :: qs_env
100 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
101 :
102 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_tensor_large_cell_Gamma'
103 :
104 : INTEGER :: handle
105 74 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma, fm_W_MIC_time
106 74 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
107 :
108 74 : CALL timeset(routineN, handle)
109 :
110 74 : CALL cite_reference(Graml2024)
111 :
112 : ! G^occ_µλ(i|τ|,k=0) = sum_n^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
113 : ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
114 : ! χ_PQ(iτ,k=0) = sum_λν [sum_µ (µν|P) G^occ_µλ(i|τ|)] [sum_σ (σλ|Q) G^vir_σν(i|τ|)]
115 74 : CALL get_mat_chi_Gamma_tau(bs_env, qs_env, bs_env%mat_chi_Gamma_tau)
116 :
117 : ! χ_PQ(iτ,k=0) -> χ_PQ(iω,k) -> ε_PQ(iω,k) -> W_PQ(iω,k) -> W^MIC_PQ(iτ) -> M^-1*W^MIC*M^-1
118 74 : CALL get_W_MIC(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_MIC_time)
119 :
120 : ! D_µν = sum_n^occ C_µn(k=0) C_νn(k=0), V^trunc_PQ = sum_cell_R <phi_P,0|V^trunc|phi_Q,R>
121 : ! Σ^x_λσ(k=0) = sum_νQ [sum_P (νσ|P) V^trunc_PQ] [sum_µ (λµ|Q) D_µν)]
122 74 : CALL get_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
123 :
124 : ! Σ^c_λσ(iτ,k=0) = sum_νQ [sum_P (νσ|P) W^MIC_PQ(iτ)] [sum_µ (λµ|Q) G^occ_µν(i|τ|)], τ < 0
125 : ! Σ^c_λσ(iτ,k=0) = sum_νQ [sum_P (νσ|P) W^MIC_PQ(iτ)] [sum_µ (λµ|Q) G^vir_µν(i|τ|)], τ > 0
126 74 : CALL get_Sigma_c(bs_env, qs_env, fm_W_MIC_time, fm_Sigma_c_Gamma_time)
127 :
128 : ! Σ^c_λσ(iτ,k=0) -> Σ^c_nn(ϵ,k); ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ,k) + Σ^x_nn(k) - v^xc_nn(k)
129 74 : CALL compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
130 :
131 74 : CALL de_init_bs_env(qs_env, bs_env)
132 :
133 74 : CALL timestop(handle)
134 :
135 74 : END SUBROUTINE gw_calc_tensor_large_cell_Gamma
136 :
137 : ! **************************************************************************************************
138 : !> \brief ...
139 : !> \param bs_env ...
140 : !> \param qs_env ...
141 : !> \param mat_chi_Gamma_tau ...
142 : ! **************************************************************************************************
143 82 : SUBROUTINE get_mat_chi_Gamma_tau(bs_env, qs_env, mat_chi_Gamma_tau)
144 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
145 : TYPE(qs_environment_type), POINTER :: qs_env
146 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
147 :
148 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_chi_Gamma_tau'
149 :
150 : INTEGER :: handle, i_intval_idx, i_t, inner_loop_atoms_interval_index, ispin, j_intval_idx
151 : INTEGER(KIND=int_8) :: flop
152 : INTEGER, DIMENSION(2) :: bounds_P, bounds_Q, i_atoms, IL_atoms, &
153 : j_atoms
154 : INTEGER, DIMENSION(2, 2) :: bounds_comb
155 : LOGICAL :: dist_too_long_i, dist_too_long_j
156 : REAL(KIND=dp) :: t1, tau
157 2050 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
158 1394 : t_3c_for_Gvir, t_3c_x_Gocc, &
159 1394 : t_3c_x_Gocc_2, t_3c_x_Gvir, &
160 738 : t_3c_x_Gvir_2
161 :
162 82 : CALL timeset(routineN, handle)
163 :
164 1498 : DO i_t = 1, bs_env%num_time_freq_points
165 :
166 1416 : t1 = m_walltime()
167 :
168 1416 : IF (bs_env%read_chi(i_t)) THEN
169 :
170 0 : CALL fm_read(bs_env%fm_RI_RI, bs_env, bs_env%chi_name, i_t)
171 :
172 : CALL copy_fm_to_dbcsr(bs_env%fm_RI_RI, mat_chi_Gamma_tau(i_t)%matrix, &
173 0 : keep_sparsity=.FALSE.)
174 :
175 0 : IF (bs_env%unit_nr > 0) THEN
176 : WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
177 0 : 'Read χ(iτ,k=0) from file for time point ', i_t, ' /', &
178 0 : bs_env%num_time_freq_points, &
179 0 : ', Execution time', m_walltime() - t1, ' s'
180 : END IF
181 :
182 : CYCLE
183 :
184 : END IF
185 :
186 1416 : IF (.NOT. bs_env%calc_chi(i_t)) CYCLE
187 :
188 : CALL create_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
189 1416 : t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2, bs_env)
190 :
191 : ! 1. compute G^occ and G^vir
192 : ! Background: G^σ(iτ) = G^occ,σ(iτ) * Θ(-τ) + G^vir,σ(iτ) * Θ(τ), σ ∈ {↑,↓}
193 : ! G^occ,σ_µλ(i|τ|,k=0) = sum_n^occ C^σ_µn(k=0) e^(-|(ϵ^σ_nk=0-ϵ_F)τ|) C^σ_λn(k=0)
194 : ! G^vir,σ_µλ(i|τ|,k=0) = sum_n^vir C^σ_µn(k=0) e^(-|(ϵ^σ_nk=0-ϵ_F)τ|) C^σ_λn(k=0)
195 1416 : tau = bs_env%time_frequency_grid%imaginary_time(i_t)
196 :
197 3032 : DO ispin = 1, bs_env%n_spin
198 1616 : CALL G_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.TRUE., vir=.FALSE.)
199 1616 : CALL G_occ_vir(bs_env, tau, bs_env%fm_Gvir, ispin, occ=.FALSE., vir=.TRUE.)
200 :
201 : CALL cp_fm_to_dbt(bs_env%fm_Gocc, bs_env%mat_ao_ao%matrix, bs_env%mat_ao_ao_tensor%matrix, t_2c_Gocc, &
202 : bs_env%para_env, bs_env%para_env_tensor%num_pe, bs_env%eps_filter, &
203 1616 : bs_env%atoms_j_t_group)
204 : CALL cp_fm_to_dbt(bs_env%fm_Gvir, bs_env%mat_ao_ao%matrix, bs_env%mat_ao_ao_tensor%matrix, t_2c_Gvir, &
205 : bs_env%para_env, bs_env%para_env_tensor%num_pe, bs_env%eps_filter, &
206 1616 : bs_env%atoms_i_t_group)
207 :
208 : ! every group has its own range of i_atoms and j_atoms; only deal with a
209 : ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
210 4648 : DO i_intval_idx = 1, bs_env%n_intervals_i
211 4848 : DO j_intval_idx = 1, bs_env%n_intervals_j
212 4848 : i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
213 4848 : j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
214 :
215 1616 : IF (bs_env%skip_chi(i_intval_idx, j_intval_idx)) THEN
216 : ! Do that only after first timestep to avoid skips due to vanishing G
217 : ! caused by gaps
218 36 : IF (i_t == 2) THEN
219 0 : bs_env%n_skip_chi = bs_env%n_skip_chi + 1
220 : END IF
221 : CYCLE
222 : END IF
223 :
224 3160 : DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
225 :
226 4740 : IL_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
227 : ! Idea: Use sparsity in 3c integrals behind χ_PQ(iτ,k=0)
228 : ! -> λ bounds from j_atoms -> sparse in IL_atoms through σ in
229 : ! N_Qλν(iτ) = sum_σ (Qλ|σ) G^vir_νσ(i|τ|,k=0)
230 : ! -> ν bounds from i_atoms -> sparse in IL_atoms through µ in
231 : ! M_Pνλ(iτ) = sum_µ (Pν|µ) G^occ_λµ(i|τ|,k=0)
232 1580 : CALL check_dist(i_atoms, IL_atoms, qs_env, bs_env, dist_too_long_i)
233 1580 : CALL check_dist(j_atoms, IL_atoms, qs_env, bs_env, dist_too_long_j)
234 1580 : IF (.NOT. dist_too_long_i) THEN
235 : ! 2. compute 3-center integrals (Pν|µ) ("|": truncated Coulomb operator)
236 : CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gocc, &
237 1580 : atoms_AO_1=i_atoms, atoms_AO_2=IL_atoms)
238 : ! 3. tensor operation M_Pνλ(iτ) = sum_µ (Pν|µ) G^occ_λµ(i|τ|,k=0)
239 : CALL G_times_3c(t_3c_for_Gocc, t_2c_Gocc, t_3c_x_Gocc, bs_env, &
240 1580 : j_atoms, i_atoms, IL_atoms)
241 : END IF
242 3160 : IF (.NOT. dist_too_long_j) THEN
243 : ! 4. compute 3-center integrals (Qλ|σ) ("|": truncated Coulomb operator)
244 : CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gvir, &
245 1580 : atoms_AO_1=j_atoms, atoms_AO_2=IL_atoms)
246 : ! 5. tensor operation N_Qλν(iτ) = sum_σ (Qλ|σ) G^vir_νσ(i|τ|,k=0)
247 : CALL G_times_3c(t_3c_for_Gvir, t_2c_Gvir, t_3c_x_Gvir, bs_env, &
248 1580 : i_atoms, j_atoms, IL_atoms)
249 : END IF
250 : END DO ! IL_atoms
251 :
252 : ! 6. reorder tensors: M_Pνλ -> M_Pλν
253 1580 : CALL dbt_copy(t_3c_x_Gocc, t_3c_x_Gocc_2, move_data=.TRUE., order=[1, 3, 2])
254 1580 : CALL dbt_copy(t_3c_x_Gvir, t_3c_x_Gvir_2, move_data=.TRUE.)
255 :
256 : ! 7. tensor operation χ_PQ(iτ,k=0) = sum_λν M_Pλν(iτ) N_Qλν(iτ),
257 : ! Bounds:
258 : ! "comb" (combined index)
259 : ! -> λ bounds from j_atoms
260 : ! -> ν bounds from i_atoms
261 : ! P -> sparse in ν (see 3.)
262 : ! Q -> sparse in λ (see 5.)
263 : bounds_comb(1:2, 1) = [bs_env%i_ao_start_from_atom(j_atoms(1)), &
264 4740 : bs_env%i_ao_end_from_atom(j_atoms(2))]
265 : bounds_comb(1:2, 2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
266 4740 : bs_env%i_ao_end_from_atom(i_atoms(2))]
267 :
268 : CALL get_bounds_from_atoms(bounds_P, i_atoms, [1, bs_env%n_atom], &
269 : bs_env%min_RI_idx_from_AO_AO_atom, &
270 4740 : bs_env%max_RI_idx_from_AO_AO_atom)
271 : CALL get_bounds_from_atoms(bounds_Q, [1, bs_env%n_atom], j_atoms, &
272 : bs_env%min_RI_idx_from_AO_AO_atom, &
273 4740 : bs_env%max_RI_idx_from_AO_AO_atom)
274 :
275 1580 : IF (bounds_Q(1) > bounds_Q(2) .OR. bounds_P(1) > bounds_P(2)) THEN
276 0 : flop = 0_int_8
277 : ELSE
278 : CALL dbt_contract(alpha=bs_env%spin_degeneracy, &
279 : tensor_1=t_3c_x_Gocc_2, tensor_2=t_3c_x_Gvir_2, &
280 : beta=1.0_dp, tensor_3=bs_env%t_chi, &
281 : contract_1=[2, 3], notcontract_1=[1], map_1=[1], &
282 : contract_2=[2, 3], notcontract_2=[1], map_2=[2], &
283 : bounds_1=bounds_comb, &
284 : bounds_2=bounds_P, &
285 : bounds_3=bounds_Q, &
286 1580 : filter_eps=bs_env%eps_filter, move_data=.FALSE., flop=flop)
287 : END IF
288 3196 : IF (flop == 0_int_8) bs_env%skip_chi(i_intval_idx, j_intval_idx) = .TRUE.
289 :
290 : END DO ! j_atoms
291 : END DO ! i_atoms
292 : END DO ! ispin
293 :
294 : ! 8. communicate data of χ_PQ(iτ,k=0) in tensor bs_env%t_chi (which local in the
295 : ! subgroup) to the global dbcsr matrix mat_chi_Gamma_tau (which stores
296 : ! χ_PQ(iτ,k=0) for all time points)
297 : CALL cp_dbt_to_dbcsr(bs_env%t_chi, bs_env%mat_RI_RI_tensor%matrix, &
298 1416 : mat_chi_Gamma_tau(i_t)%matrix, bs_env%para_env)
299 :
300 : CALL fm_write(bs_env%fm_RI_RI, i_t, bs_env%chi_name, qs_env, &
301 1416 : matrix=mat_chi_Gamma_tau(i_t)%matrix)
302 :
303 : CALL destroy_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
304 1416 : t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2)
305 :
306 1498 : IF (bs_env%unit_nr > 0) THEN
307 : WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F10.1,A)') &
308 708 : 'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
309 1416 : ', Execution time', m_walltime() - t1, ' s'
310 : END IF
311 :
312 : END DO ! i_t
313 :
314 82 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
315 :
316 82 : CALL timestop(handle)
317 :
318 82 : END SUBROUTINE get_mat_chi_Gamma_tau
319 :
320 : ! **************************************************************************************************
321 : !> \brief ...
322 : !> \param t_2c_Gocc ...
323 : !> \param t_2c_Gvir ...
324 : !> \param t_3c_for_Gocc ...
325 : !> \param t_3c_for_Gvir ...
326 : !> \param t_3c_x_Gocc ...
327 : !> \param t_3c_x_Gvir ...
328 : !> \param t_3c_x_Gocc_2 ...
329 : !> \param t_3c_x_Gvir_2 ...
330 : !> \param bs_env ...
331 : ! **************************************************************************************************
332 1416 : SUBROUTINE create_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
333 : t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2, bs_env)
334 :
335 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
336 : t_3c_for_Gvir, t_3c_x_Gocc, &
337 : t_3c_x_Gvir, t_3c_x_Gocc_2, &
338 : t_3c_x_Gvir_2
339 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
340 :
341 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_tensors_chi'
342 :
343 : INTEGER :: handle
344 :
345 1416 : CALL timeset(routineN, handle)
346 :
347 1416 : CALL dbt_create(bs_env%t_G, t_2c_Gocc, name="Gocc 2c (AO|AO)")
348 1416 : CALL dbt_create(bs_env%t_G, t_2c_Gvir, name="Gvir 2c (AO|AO)")
349 1416 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gocc, name="Gocc 3c (RI AO|AO)")
350 1416 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gvir, name="Gvir 3c (RI AO|AO)")
351 1416 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gocc, name="xGocc 3c (RI AO|AO)")
352 1416 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gvir, name="xGvir 3c (RI AO|AO)")
353 1416 : CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gocc_2, name="x2Gocc 3c (RI AO|AO)")
354 1416 : CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gvir_2, name="x2Gvir 3c (RI AO|AO)")
355 :
356 1416 : CALL timestop(handle)
357 :
358 1416 : END SUBROUTINE create_tensors_chi
359 :
360 : ! **************************************************************************************************
361 : !> \brief ...
362 : !> \param t_2c_Gocc ...
363 : !> \param t_2c_Gvir ...
364 : !> \param t_3c_for_Gocc ...
365 : !> \param t_3c_for_Gvir ...
366 : !> \param t_3c_x_Gocc ...
367 : !> \param t_3c_x_Gvir ...
368 : !> \param t_3c_x_Gocc_2 ...
369 : !> \param t_3c_x_Gvir_2 ...
370 : ! **************************************************************************************************
371 1416 : SUBROUTINE destroy_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
372 : t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2)
373 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
374 : t_3c_for_Gvir, t_3c_x_Gocc, &
375 : t_3c_x_Gvir, t_3c_x_Gocc_2, &
376 : t_3c_x_Gvir_2
377 :
378 : CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_tensors_chi'
379 :
380 : INTEGER :: handle
381 :
382 1416 : CALL timeset(routineN, handle)
383 :
384 1416 : CALL dbt_destroy(t_2c_Gocc)
385 1416 : CALL dbt_destroy(t_2c_Gvir)
386 1416 : CALL dbt_destroy(t_3c_for_Gocc)
387 1416 : CALL dbt_destroy(t_3c_for_Gvir)
388 1416 : CALL dbt_destroy(t_3c_x_Gocc)
389 1416 : CALL dbt_destroy(t_3c_x_Gvir)
390 1416 : CALL dbt_destroy(t_3c_x_Gocc_2)
391 1416 : CALL dbt_destroy(t_3c_x_Gvir_2)
392 :
393 1416 : CALL timestop(handle)
394 :
395 1416 : END SUBROUTINE destroy_tensors_chi
396 :
397 : ! **************************************************************************************************
398 : !> \brief ...
399 : !> \param qs_env ...
400 : !> \param bs_env ...
401 : !> \param t_3c ...
402 : !> \param atoms_AO_1 ...
403 : !> \param atoms_AO_2 ...
404 : !> \param atoms_RI ...
405 : ! **************************************************************************************************
406 7692 : SUBROUTINE compute_3c_integrals(qs_env, bs_env, t_3c, atoms_AO_1, atoms_AO_2, atoms_RI)
407 : TYPE(qs_environment_type), POINTER :: qs_env
408 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
409 : TYPE(dbt_type) :: t_3c
410 : INTEGER, DIMENSION(2), OPTIONAL :: atoms_AO_1, atoms_AO_2, atoms_RI
411 :
412 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_integrals'
413 :
414 : INTEGER :: handle
415 7692 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_array
416 :
417 7692 : CALL timeset(routineN, handle)
418 :
419 : ! free memory (not clear whether memory has been freed previously)
420 7692 : CALL dbt_clear(t_3c)
421 :
422 84612 : ALLOCATE (t_3c_array(1, 1))
423 7692 : CALL dbt_create(t_3c, t_3c_array(1, 1))
424 :
425 : CALL build_3c_integrals(t_3c_array, &
426 : bs_env%eps_filter, &
427 : qs_env, &
428 : bs_env%nl_3c, &
429 : int_eps=bs_env%eps_filter, &
430 : basis_i=bs_env%basis_set_RI, &
431 : basis_j=bs_env%basis_set_AO, &
432 : basis_k=bs_env%basis_set_AO, &
433 : potential_parameter=bs_env%ri_metric, &
434 : bounds_i=atoms_RI, &
435 : bounds_j=atoms_AO_1, &
436 : bounds_k=atoms_AO_2, &
437 7692 : desymmetrize=.FALSE.)
438 :
439 7692 : CALL dbt_filter(t_3c_array(1, 1), bs_env%eps_filter)
440 :
441 7692 : CALL dbt_copy(t_3c_array(1, 1), t_3c, move_data=.TRUE.)
442 :
443 7692 : CALL dbt_destroy(t_3c_array(1, 1))
444 15384 : DEALLOCATE (t_3c_array)
445 :
446 7692 : CALL timestop(handle)
447 :
448 15384 : END SUBROUTINE compute_3c_integrals
449 :
450 : ! **************************************************************************************************
451 : !> \brief ...
452 : !> \param t_3c_for_G ...
453 : !> \param t_G ...
454 : !> \param t_M ...
455 : !> \param bs_env ...
456 : !> \param atoms_AO_1 ...
457 : !> \param atoms_AO_2 ...
458 : !> \param atoms_IL ...
459 : ! **************************************************************************************************
460 3160 : SUBROUTINE G_times_3c(t_3c_for_G, t_G, t_M, bs_env, atoms_AO_1, atoms_AO_2, atoms_IL)
461 : TYPE(dbt_type) :: t_3c_for_G, t_G, t_M
462 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
463 : INTEGER, DIMENSION(2) :: atoms_AO_1, atoms_AO_2, atoms_IL
464 :
465 : CHARACTER(LEN=*), PARAMETER :: routineN = 'G_times_3c'
466 :
467 : INTEGER :: handle
468 : INTEGER(KIND=int_8) :: flop
469 : INTEGER, DIMENSION(2) :: bounds_ao_1, bounds_IL
470 : INTEGER, DIMENSION(2, 2) :: bounds_comb
471 :
472 3160 : CALL timeset(routineN, handle)
473 :
474 : ! Bounds reduce needed memory and therefore scaling behavior
475 : ! Operations are of the form, e.g, M_Pνλ = sum_µ (Pν|µ) G_λµ
476 : ! "comb" (combined index)
477 : ! -> P sparse in ν and µ
478 : ! -> λ bounds from j_atoms (via atoms_AO_1)
479 : ! µ bounds from inner loop "IL" indices and sparse in P and ν
480 : ! ν bounds from i_atoms (via atoms_AO_2) and sparse in P and µ
481 :
482 : ! µ index
483 : CALL get_bounds_from_atoms(bounds_IL, [1, bs_env%n_atom], atoms_AO_2, &
484 : bs_env%min_AO_idx_from_RI_AO_atom, &
485 : bs_env%max_AO_idx_from_RI_AO_atom, &
486 : atoms_3=atoms_IL, &
487 : indices_3_start=bs_env%i_ao_start_from_atom, &
488 9480 : indices_3_end=bs_env%i_ao_end_from_atom)
489 :
490 : ! P index
491 : CALL get_bounds_from_atoms(bounds_comb(:, 1), atoms_IL, atoms_AO_2, &
492 : bs_env%min_RI_idx_from_AO_AO_atom, &
493 3160 : bs_env%max_RI_idx_from_AO_AO_atom)
494 :
495 : ! ν index
496 : CALL get_bounds_from_atoms(bounds_comb(:, 2), [1, bs_env%n_atom], atoms_IL, &
497 : bs_env%min_AO_idx_from_RI_AO_atom, &
498 : bs_env%max_AO_idx_from_RI_AO_atom, &
499 : atoms_3=atoms_AO_2, &
500 : indices_3_start=bs_env%i_ao_start_from_atom, &
501 9480 : indices_3_end=bs_env%i_ao_end_from_atom)
502 :
503 : ! λ index
504 : bounds_ao_1(1:2) = [bs_env%i_ao_start_from_atom(atoms_AO_1(1)), &
505 9480 : bs_env%i_ao_end_from_atom(atoms_AO_1(2))]
506 :
507 3160 : IF (bounds_IL(1) > bounds_IL(2) .OR. bounds_comb(1, 2) > bounds_comb(2, 2)) THEN
508 : flop = 0_int_8
509 : ELSE
510 : CALL dbt_contract(alpha=1.0_dp, &
511 : tensor_1=t_3c_for_G, &
512 : tensor_2=t_G, &
513 : beta=1.0_dp, &
514 : tensor_3=t_M, &
515 : contract_1=[3], notcontract_1=[1, 2], map_1=[1, 2], &
516 : contract_2=[2], notcontract_2=[1], map_2=[3], &
517 : bounds_1=bounds_IL, &
518 : bounds_2=bounds_comb, &
519 : bounds_3=bounds_ao_1, &
520 : flop=flop, &
521 3160 : filter_eps=bs_env%eps_filter)
522 : END IF
523 :
524 3160 : CALL dbt_clear(t_3c_for_G)
525 :
526 3160 : CALL timestop(handle)
527 :
528 3160 : END SUBROUTINE G_times_3c
529 :
530 : ! **************************************************************************************************
531 : !> \brief ...
532 : !> \param atoms_1 ...
533 : !> \param atoms_2 ...
534 : !> \param qs_env ...
535 : !> \param bs_env ...
536 : !> \param dist_too_long ...
537 : ! **************************************************************************************************
538 3160 : SUBROUTINE check_dist(atoms_1, atoms_2, qs_env, bs_env, dist_too_long)
539 : INTEGER, DIMENSION(2) :: atoms_1, atoms_2
540 : TYPE(qs_environment_type), POINTER :: qs_env
541 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
542 : LOGICAL :: dist_too_long
543 :
544 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_dist'
545 :
546 : INTEGER :: atom_1, atom_2, handle
547 : REAL(dp) :: abs_rab, min_dist_AO_atoms
548 : REAL(KIND=dp), DIMENSION(3) :: rab
549 : TYPE(cell_type), POINTER :: cell
550 3160 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
551 :
552 3160 : CALL timeset(routineN, handle)
553 :
554 3160 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
555 :
556 3160 : min_dist_AO_atoms = HUGE(1.0_dp)
557 9936 : DO atom_1 = atoms_1(1), atoms_1(2)
558 24856 : DO atom_2 = atoms_2(1), atoms_2(2)
559 14920 : rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
560 :
561 14920 : abs_rab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
562 :
563 21696 : min_dist_AO_atoms = MIN(min_dist_AO_atoms, abs_rab)
564 : END DO
565 : END DO
566 :
567 3160 : dist_too_long = (min_dist_AO_atoms > bs_env%max_dist_AO_atoms)
568 :
569 3160 : CALL timestop(handle)
570 :
571 3160 : END SUBROUTINE check_dist
572 :
573 : ! **************************************************************************************************
574 : !> \brief ...
575 : !> \param bs_env ...
576 : !> \param qs_env ...
577 : !> \param fm_Sigma_x_Gamma ...
578 : ! **************************************************************************************************
579 82 : SUBROUTINE get_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
580 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
581 : TYPE(qs_environment_type), POINTER :: qs_env
582 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
583 :
584 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_Sigma_x'
585 :
586 : INTEGER :: handle, ispin
587 :
588 82 : CALL timeset(routineN, handle)
589 :
590 340 : ALLOCATE (fm_Sigma_x_Gamma(bs_env%n_spin))
591 176 : DO ispin = 1, bs_env%n_spin
592 176 : CALL cp_fm_create(fm_Sigma_x_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
593 : END DO
594 :
595 82 : IF (bs_env%Sigma_x_exists) THEN
596 14 : DO ispin = 1, bs_env%n_spin
597 14 : CALL fm_read(fm_Sigma_x_Gamma(ispin), bs_env, bs_env%Sigma_x_name, ispin)
598 : END DO
599 : ELSE
600 76 : CALL compute_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
601 : END IF
602 :
603 82 : CALL timestop(handle)
604 :
605 82 : END SUBROUTINE get_Sigma_x
606 :
607 : ! **************************************************************************************************
608 : !> \brief ...
609 : !> \param bs_env ...
610 : !> \param qs_env ...
611 : !> \param fm_Sigma_x_Gamma ...
612 : ! **************************************************************************************************
613 76 : SUBROUTINE compute_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
614 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
615 : TYPE(qs_environment_type), POINTER :: qs_env
616 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
617 :
618 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_x'
619 :
620 : INTEGER :: handle, i_intval_idx, ispin, j_intval_idx
621 : INTEGER, DIMENSION(2) :: i_atoms, j_atoms
622 : REAL(KIND=dp) :: t1
623 76 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Vtr_Gamma
624 : TYPE(dbcsr_type) :: mat_Sigma_x_Gamma
625 2508 : TYPE(dbt_type) :: t_2c_D, t_2c_Sigma_x, t_2c_V, t_3c_x_V
626 :
627 76 : CALL timeset(routineN, handle)
628 :
629 76 : t1 = m_walltime()
630 :
631 76 : CALL dbt_create(bs_env%t_G, t_2c_D)
632 76 : CALL dbt_create(bs_env%t_W, t_2c_V)
633 76 : CALL dbt_create(bs_env%t_G, t_2c_Sigma_x)
634 76 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_V)
635 76 : CALL dbcsr_create(mat_Sigma_x_Gamma, template=bs_env%mat_ao_ao%matrix)
636 :
637 : ! 1. Compute truncated Coulomb operator matrix V^tr(k=0) (cutoff rad: cellsize/2)
638 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI%matrix_struct, bs_env%n_RI, &
639 76 : bs_env%trunc_coulomb)
640 :
641 : ! M^-1(k=0) V^tr M^-1(k=0) -> fm_Vtr_Gamma
642 76 : CALL fm_contract_ABA(bs_env%fm_Minv_Gamma, fm_Vtr_Gamma(:, 1))
643 :
644 162 : DO ispin = 1, bs_env%n_spin
645 :
646 : ! 3. Compute density matrix D_µν
647 86 : CALL G_occ_vir(bs_env, 0.0_dp, bs_env%fm_work_mo(2), ispin, occ=.TRUE., vir=.FALSE.)
648 :
649 : CALL cp_fm_to_dbt(bs_env%fm_work_mo(2), bs_env%mat_ao_ao%matrix, bs_env%mat_ao_ao_tensor%matrix, t_2c_D, &
650 : bs_env%para_env, bs_env%para_env_tensor%num_pe, bs_env%eps_filter, &
651 86 : bs_env%atoms_i_t_group)
652 :
653 : CALL cp_fm_to_dbt(fm_Vtr_Gamma(1, 1), bs_env%mat_RI_RI%matrix, bs_env%mat_RI_RI_tensor%matrix, t_2c_V, &
654 : bs_env%para_env, bs_env%para_env_tensor%num_pe, bs_env%eps_filter, &
655 86 : bs_env%atoms_j_t_group)
656 :
657 : ! every group has its own range of i_atoms and j_atoms; only deal with a
658 : ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
659 172 : DO i_intval_idx = 1, bs_env%n_intervals_i
660 258 : DO j_intval_idx = 1, bs_env%n_intervals_j
661 258 : i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
662 258 : j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
663 :
664 : ! 4. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
665 : ! 5. M_Qνσ(iτ) = sum_P (νσ|P) (M^-1(k=0)*V^tr(k=0)*M^-1(k=0))_QP(iτ)
666 86 : CALL compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_V, t_2c_V)
667 :
668 : ! 6. tensor operations with D and computation of Σ^x
669 : ! Σ^x_λσ(k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) D_νµ
670 : CALL contract_to_Sigma(t_2c_D, t_3c_x_V, t_2c_Sigma_x, i_atoms, j_atoms, &
671 172 : qs_env, bs_env, occ=.TRUE., vir=.FALSE.)
672 :
673 : END DO ! j_atoms
674 : END DO ! i_atoms
675 :
676 : CALL cp_dbt_to_dbcsr(t_2c_Sigma_x, bs_env%mat_ao_ao_tensor%matrix, &
677 86 : mat_Sigma_x_Gamma, bs_env%para_env)
678 :
679 : CALL fm_write(bs_env%fm_work_mo(1), ispin, bs_env%Sigma_x_name, qs_env, &
680 86 : matrix=mat_Sigma_x_Gamma)
681 :
682 162 : CALL copy_dbcsr_to_fm(mat_Sigma_x_Gamma, fm_Sigma_x_Gamma(ispin))
683 :
684 : END DO ! ispin
685 :
686 76 : IF (bs_env%unit_nr > 0) THEN
687 : WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
688 38 : 'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
689 38 : WRITE (bs_env%unit_nr, '(A)') ' '
690 : END IF
691 :
692 76 : CALL dbcsr_release(mat_Sigma_x_Gamma)
693 76 : CALL dbt_destroy(t_2c_D)
694 76 : CALL dbt_destroy(t_2c_V)
695 76 : CALL dbt_destroy(t_2c_Sigma_x)
696 76 : CALL dbt_destroy(t_3c_x_V)
697 76 : CALL cp_fm_release(fm_Vtr_Gamma)
698 :
699 76 : CALL timestop(handle)
700 :
701 152 : END SUBROUTINE compute_Sigma_x
702 :
703 : ! **************************************************************************************************
704 : !> \brief ...
705 : !> \param bs_env ...
706 : !> \param qs_env ...
707 : !> \param fm_W_MIC_time ...
708 : !> \param fm_Sigma_c_Gamma_time ...
709 : ! **************************************************************************************************
710 74 : SUBROUTINE get_Sigma_c(bs_env, qs_env, fm_W_MIC_time, fm_Sigma_c_Gamma_time)
711 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
712 : TYPE(qs_environment_type), POINTER :: qs_env
713 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
714 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
715 :
716 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_Sigma_c'
717 :
718 : INTEGER :: handle, i_t, ispin, read_write_index
719 : REAL(KIND=dp) :: t1, tau
720 74 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
721 1258 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, &
722 666 : t_2c_Sigma_neg_tau, &
723 1850 : t_2c_Sigma_pos_tau, t_2c_W, t_3c_x_W
724 :
725 74 : CALL timeset(routineN, handle)
726 :
727 : CALL create_mat_for_Sigma_c(bs_env, t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
728 : t_2c_Sigma_pos_tau, t_3c_x_W, &
729 74 : mat_Sigma_neg_tau, mat_Sigma_pos_tau)
730 :
731 1438 : DO i_t = 1, bs_env%num_time_freq_points
732 :
733 3002 : DO ispin = 1, bs_env%n_spin
734 :
735 1564 : t1 = m_walltime()
736 :
737 1564 : read_write_index = i_t + (ispin - 1)*bs_env%num_time_freq_points
738 :
739 : ! read self-energy from restart
740 1564 : IF (bs_env%Sigma_c_exists(i_t, ispin)) THEN
741 120 : CALL fm_read(bs_env%fm_work_mo(1), bs_env, bs_env%Sigma_p_name, read_write_index)
742 : CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_Sigma_pos_tau(i_t, ispin)%matrix, &
743 120 : keep_sparsity=.FALSE.)
744 120 : CALL fm_read(bs_env%fm_work_mo(1), bs_env, bs_env%Sigma_n_name, read_write_index)
745 : CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_Sigma_neg_tau(i_t, ispin)%matrix, &
746 120 : keep_sparsity=.FALSE.)
747 120 : IF (bs_env%unit_nr > 0) THEN
748 60 : WRITE (bs_env%unit_nr, '(T2,2A,I3,A,I3,A,F10.1,A)') 'Read Σ^c(iτ,k=0) ', &
749 60 : 'from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
750 120 : ', Execution time', m_walltime() - t1, ' s'
751 : END IF
752 :
753 : CYCLE
754 :
755 : END IF
756 :
757 1444 : tau = bs_env%time_frequency_grid%imaginary_time(i_t)
758 :
759 : CALL compute_Sigma_c_occ_vir(bs_env, qs_env, tau, ispin, fm_W_MIC_time(i_t), &
760 : t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
761 : t_2c_Sigma_pos_tau, t_3c_x_W, &
762 : mat_Sigma_neg_tau(i_t, ispin)%matrix, &
763 : mat_Sigma_pos_tau(i_t, ispin)%matrix, &
764 1444 : count_skipped=(i_t == 2))
765 :
766 : CALL fm_write(bs_env%fm_work_mo(1), read_write_index, bs_env%Sigma_p_name, qs_env, &
767 1444 : matrix=mat_Sigma_pos_tau(i_t, ispin)%matrix)
768 : CALL fm_write(bs_env%fm_work_mo(1), read_write_index, bs_env%Sigma_n_name, qs_env, &
769 1444 : matrix=mat_Sigma_neg_tau(i_t, ispin)%matrix)
770 :
771 2808 : IF (bs_env%unit_nr > 0) THEN
772 : WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F10.1,A)') &
773 722 : 'Computed Σ^c(iτ,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
774 1444 : ', Execution time', m_walltime() - t1, ' s'
775 : END IF
776 :
777 : END DO ! ispin
778 :
779 : END DO ! i_t
780 :
781 74 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
782 :
783 : CALL fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
784 74 : mat_Sigma_pos_tau, mat_Sigma_neg_tau)
785 :
786 74 : CALL print_skipping(bs_env)
787 :
788 : CALL destroy_mat_Sigma_c(t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
789 : t_2c_Sigma_pos_tau, t_3c_x_W, fm_W_MIC_time, &
790 74 : mat_Sigma_neg_tau, mat_Sigma_pos_tau)
791 :
792 74 : CALL delete_unnecessary_files(bs_env)
793 :
794 74 : CALL timestop(handle)
795 :
796 148 : END SUBROUTINE get_Sigma_c
797 :
798 : ! **************************************************************************************************
799 : !> \brief ...
800 : !> \brief Occupied and virtual parts of the correlation self-energy at one imaginary time τ
801 : !> Σ^neg_λσ(iτ) = -sum_νQ [sum_P (νσ|P) W_PQ] [sum_µ (λµ|Q) G^occ_µν(i|τ|)]
802 : !> Σ^pos_λσ(iτ) = sum_νQ [sum_P (νσ|P) W_PQ] [sum_µ (λµ|Q) G^vir_µν(i|τ|)]
803 : !> (G0W0: W = W^MIC(iτ); static COHSEX: W = W^MIC(iω=0))
804 : !> \param bs_env ...
805 : !> \param qs_env ...
806 : !> \param tau imaginary time τ
807 : !> \param ispin ...
808 : !> \param fm_W W^MIC_PQ
809 : !> \param t_2c_Gocc ...
810 : !> \param t_2c_Gvir ...
811 : !> \param t_2c_W ...
812 : !> \param t_2c_Sigma_neg work tensor for Σ^neg
813 : !> \param t_2c_Sigma_pos work tensor for Σ^pos
814 : !> \param t_3c_x_W ...
815 : !> \param mat_Sigma_neg Σ^neg_λσ(iτ) on exit
816 : !> \param mat_Sigma_pos Σ^pos_λσ(iτ) on exit
817 : !> \param count_skipped count the skipped atom-interval pairs (for print_skipping)
818 : ! **************************************************************************************************
819 1452 : SUBROUTINE compute_Sigma_c_occ_vir(bs_env, qs_env, tau, ispin, fm_W, t_2c_Gocc, t_2c_Gvir, &
820 : t_2c_W, t_2c_Sigma_neg, t_2c_Sigma_pos, t_3c_x_W, &
821 : mat_Sigma_neg, mat_Sigma_pos, count_skipped)
822 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
823 : TYPE(qs_environment_type), POINTER :: qs_env
824 : REAL(KIND=dp) :: tau
825 : INTEGER :: ispin
826 : TYPE(cp_fm_type) :: fm_W
827 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_2c_W, &
828 : t_2c_Sigma_neg, t_2c_Sigma_pos, &
829 : t_3c_x_W
830 : TYPE(dbcsr_type) :: mat_Sigma_neg, mat_Sigma_pos
831 : LOGICAL :: count_skipped
832 :
833 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c_occ_vir'
834 :
835 : INTEGER :: handle, i_intval_idx, j_intval_idx
836 : INTEGER, DIMENSION(2) :: i_atoms, j_atoms
837 :
838 1452 : CALL timeset(routineN, handle)
839 :
840 1452 : CALL G_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.TRUE., vir=.FALSE.)
841 1452 : CALL G_occ_vir(bs_env, tau, bs_env%fm_Gvir, ispin, occ=.FALSE., vir=.TRUE.)
842 :
843 : ! fm G^occ, G^vir and W to local tensor
844 : CALL cp_fm_to_dbt(bs_env%fm_Gocc, bs_env%mat_ao_ao%matrix, bs_env%mat_ao_ao_tensor%matrix, t_2c_Gocc, &
845 : bs_env%para_env, bs_env%para_env_tensor%num_pe, bs_env%eps_filter, &
846 1452 : bs_env%atoms_i_t_group)
847 : CALL cp_fm_to_dbt(bs_env%fm_Gvir, bs_env%mat_ao_ao%matrix, bs_env%mat_ao_ao_tensor%matrix, t_2c_Gvir, &
848 : bs_env%para_env, bs_env%para_env_tensor%num_pe, bs_env%eps_filter, &
849 1452 : bs_env%atoms_i_t_group)
850 : CALL cp_fm_to_dbt(fm_W, bs_env%mat_RI_RI%matrix, bs_env%mat_RI_RI_tensor%matrix, t_2c_W, &
851 : bs_env%para_env, bs_env%para_env_tensor%num_pe, bs_env%eps_filter, &
852 1452 : bs_env%atoms_j_t_group)
853 :
854 : ! every group has its own range of i_atoms and j_atoms; only deal with a
855 : ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
856 2904 : DO i_intval_idx = 1, bs_env%n_intervals_i
857 4356 : DO j_intval_idx = 1, bs_env%n_intervals_j
858 4356 : i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
859 4356 : j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
860 :
861 1452 : IF (bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx) .AND. &
862 : bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) THEN
863 : ! Do that only after first timestep to avoid skips due to vanishing G
864 : ! caused by gaps
865 18 : IF (count_skipped) THEN
866 0 : bs_env%n_skip_sigma = bs_env%n_skip_sigma + 1
867 : END IF
868 : CYCLE
869 : END IF
870 :
871 : ! 1. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
872 : ! 2. tensor operation M_Qνσ(iτ) = sum_P (νσ|P) W^MIC_QP(iτ)
873 1434 : CALL compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_W, t_2c_W)
874 :
875 : ! 3. Σ_λσ(iτ,k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) G^occ_νµ(i|τ|) for τ < 0
876 : ! (recall M_Qνσ(iτ) = M_Qνσ(-iτ) because W^MIC_PQ(iτ) = W^MIC_PQ(-iτ) )
877 : CALL contract_to_Sigma(t_2c_Gocc, t_3c_x_W, t_2c_Sigma_neg, i_atoms, j_atoms, &
878 : qs_env, bs_env, occ=.TRUE., vir=.FALSE., &
879 1434 : can_skip=bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx))
880 :
881 : ! Σ_λσ(iτ,k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) G^vir_νµ(i|τ|) for τ > 0
882 : CALL contract_to_Sigma(t_2c_Gvir, t_3c_x_W, t_2c_Sigma_pos, i_atoms, j_atoms, &
883 : qs_env, bs_env, occ=.FALSE., vir=.TRUE., &
884 2886 : can_skip=bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx))
885 :
886 : END DO ! j_atoms
887 : END DO ! i_atoms
888 :
889 : ! 4. communicate data tensor t_2c_Sigma (which is local in the subgroup)
890 : ! to the global dbcsr matrix mat_Sigma_pos/neg_tau (which stores Σ for all iτ)
891 : CALL cp_dbt_to_dbcsr(t_2c_Sigma_neg, bs_env%mat_ao_ao_tensor%matrix, &
892 1452 : mat_Sigma_neg, bs_env%para_env)
893 : CALL cp_dbt_to_dbcsr(t_2c_Sigma_pos, bs_env%mat_ao_ao_tensor%matrix, &
894 1452 : mat_Sigma_pos, bs_env%para_env)
895 :
896 1452 : CALL timestop(handle)
897 :
898 1452 : END SUBROUTINE compute_Sigma_c_occ_vir
899 :
900 : ! **************************************************************************************************
901 : !> \brief ...
902 : !> \param bs_env ...
903 : !> \param t_2c_Gocc ...
904 : !> \param t_2c_Gvir ...
905 : !> \param t_2c_W ...
906 : !> \param t_2c_Sigma_neg_tau ...
907 : !> \param t_2c_Sigma_pos_tau ...
908 : !> \param t_3c_x_W ...
909 : !> \param mat_Sigma_neg_tau ...
910 : !> \param mat_Sigma_pos_tau ...
911 : ! **************************************************************************************************
912 82 : SUBROUTINE create_mat_for_Sigma_c(bs_env, t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
913 : t_2c_Sigma_pos_tau, t_3c_x_W, &
914 : mat_Sigma_neg_tau, mat_Sigma_pos_tau)
915 :
916 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
917 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_2c_W, &
918 : t_2c_Sigma_neg_tau, &
919 : t_2c_Sigma_pos_tau, t_3c_x_W
920 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
921 :
922 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_mat_for_Sigma_c'
923 :
924 : INTEGER :: handle, i_t, ispin
925 :
926 82 : CALL timeset(routineN, handle)
927 :
928 82 : CALL dbt_create(bs_env%t_G, t_2c_Gocc)
929 82 : CALL dbt_create(bs_env%t_G, t_2c_Gvir)
930 82 : CALL dbt_create(bs_env%t_W, t_2c_W)
931 82 : CALL dbt_create(bs_env%t_G, t_2c_Sigma_neg_tau)
932 82 : CALL dbt_create(bs_env%t_G, t_2c_Sigma_pos_tau)
933 82 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_W)
934 :
935 82 : NULLIFY (mat_Sigma_neg_tau, mat_Sigma_pos_tau)
936 2038 : ALLOCATE (mat_Sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
937 2038 : ALLOCATE (mat_Sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
938 :
939 176 : DO ispin = 1, bs_env%n_spin
940 1792 : DO i_t = 1, bs_env%num_time_freq_points
941 1616 : ALLOCATE (mat_Sigma_neg_tau(i_t, ispin)%matrix)
942 1616 : ALLOCATE (mat_Sigma_pos_tau(i_t, ispin)%matrix)
943 1616 : CALL dbcsr_create(mat_Sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
944 1710 : CALL dbcsr_create(mat_Sigma_pos_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
945 : END DO
946 : END DO
947 :
948 82 : CALL timestop(handle)
949 :
950 82 : END SUBROUTINE create_mat_for_Sigma_c
951 :
952 : ! **************************************************************************************************
953 : !> \brief ...
954 : !> \param qs_env ...
955 : !> \param bs_env ...
956 : !> \param i_atoms ...
957 : !> \param j_atoms ...
958 : !> \param t_3c_x_W ...
959 : !> \param t_2c_W ...
960 : ! **************************************************************************************************
961 1520 : SUBROUTINE compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_W, t_2c_W)
962 :
963 : TYPE(qs_environment_type), POINTER :: qs_env
964 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
965 : INTEGER, DIMENSION(2) :: i_atoms, j_atoms
966 : TYPE(dbt_type) :: t_3c_x_W, t_2c_W
967 :
968 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_and_contract_W'
969 :
970 : INTEGER :: handle, RI_intval_idx
971 : INTEGER(KIND=int_8) :: flop
972 : INTEGER, DIMENSION(2) :: bounds_P, bounds_Q, RI_atoms
973 : INTEGER, DIMENSION(2, 2) :: bounds_ao
974 25840 : TYPE(dbt_type) :: t_3c_for_W, t_3c_x_W_tmp
975 :
976 1520 : CALL timeset(routineN, handle)
977 :
978 1520 : CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_W_tmp)
979 1520 : CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_for_W)
980 :
981 : ! final layout will be: M_Qνσ(iτ) = sum_P (P|νσ) W^MIC_QP(iτ)
982 : ! Bounds:
983 : ! "AO"
984 : ! -> ν (AO_1 in compute_3c_integrals) bounds from i_atoms and sparse in σ and P
985 : ! -> σ (AO_2 in compute_3c_integrals) sparse in ν and P
986 : ! Q bounds from j_atoms
987 : ! P bounds from inner loop indices and sparse in ν and σ
988 :
989 : bounds_Q(1:2) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
990 4560 : bs_env%i_RI_end_from_atom(j_atoms(2))]
991 :
992 3040 : DO RI_intval_idx = 1, bs_env%n_intervals_inner_loop_atoms
993 4560 : RI_atoms = bs_env%inner_loop_atom_intervals(1:2, RI_intval_idx)
994 :
995 : CALL get_bounds_from_atoms(bounds_P, i_atoms, [1, bs_env%n_atom], &
996 : bs_env%min_RI_idx_from_AO_AO_atom, &
997 : bs_env%max_RI_idx_from_AO_AO_atom, &
998 : atoms_3=RI_atoms, &
999 : indices_3_start=bs_env%i_RI_start_from_atom, &
1000 4560 : indices_3_end=bs_env%i_RI_end_from_atom)
1001 :
1002 : ! σ
1003 : CALL get_bounds_from_atoms(bounds_ao(:, 2), RI_atoms, i_atoms, &
1004 : bs_env%min_AO_idx_from_RI_AO_atom, &
1005 1520 : bs_env%max_AO_idx_from_RI_AO_atom)
1006 : ! ν
1007 : CALL get_bounds_from_atoms(bounds_ao(:, 1), RI_atoms, [1, bs_env%n_atom], &
1008 : bs_env%min_AO_idx_from_RI_AO_atom, &
1009 : bs_env%max_AO_idx_from_RI_AO_atom, &
1010 : atoms_3=i_atoms, &
1011 : indices_3_start=bs_env%i_ao_start_from_atom, &
1012 4560 : indices_3_end=bs_env%i_ao_end_from_atom)
1013 :
1014 1520 : IF (bounds_P(1) > bounds_P(2) .OR. bounds_ao(1, 2) > bounds_ao(2, 2)) THEN
1015 : CYCLE
1016 : END IF
1017 :
1018 : ! 1. compute 3-center integrals (P|µν) ("|": truncated Coulomb operator)
1019 : CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_W, &
1020 1520 : atoms_AO_1=i_atoms, atoms_RI=RI_atoms)
1021 :
1022 : ! 2. tensor operation M_Qνσ(iτ) = sum_P W^MIC_QP(iτ) (P|νσ)
1023 : CALL dbt_contract(alpha=1.0_dp, &
1024 : tensor_1=t_2c_W, &
1025 : tensor_2=t_3c_for_W, &
1026 : beta=1.0_dp, &
1027 : tensor_3=t_3c_x_W_tmp, &
1028 : contract_1=[2], notcontract_1=[1], map_1=[1], &
1029 : contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3], &
1030 : bounds_1=bounds_P, &
1031 : bounds_2=bounds_Q, &
1032 : bounds_3=bounds_ao, &
1033 : flop=flop, &
1034 : move_data=.FALSE., &
1035 3040 : filter_eps=bs_env%eps_filter)
1036 :
1037 : END DO ! RI_atoms
1038 :
1039 : ! 3. reorder tensor
1040 1520 : CALL dbt_copy(t_3c_x_W_tmp, t_3c_x_W, order=[1, 2, 3], move_data=.TRUE.)
1041 :
1042 1520 : CALL dbt_destroy(t_3c_x_W_tmp)
1043 1520 : CALL dbt_destroy(t_3c_for_W)
1044 :
1045 1520 : CALL timestop(handle)
1046 :
1047 1520 : END SUBROUTINE compute_3c_and_contract_W
1048 :
1049 : ! **************************************************************************************************
1050 : !> \brief ...
1051 : !> \param t_2c_G ...
1052 : !> \param t_3c_x_W ...
1053 : !> \param t_2c_Sigma ...
1054 : !> \param i_atoms ...
1055 : !> \param j_atoms ...
1056 : !> \param qs_env ...
1057 : !> \param bs_env ...
1058 : !> \param occ ...
1059 : !> \param vir ...
1060 : !> \param can_skip ...
1061 : ! **************************************************************************************************
1062 2954 : SUBROUTINE contract_to_Sigma(t_2c_G, t_3c_x_W, t_2c_Sigma, i_atoms, j_atoms, qs_env, bs_env, &
1063 : occ, vir, can_skip)
1064 : TYPE(dbt_type) :: t_2c_G, t_3c_x_W, t_2c_Sigma
1065 : INTEGER, DIMENSION(2) :: i_atoms, j_atoms
1066 : TYPE(qs_environment_type), POINTER :: qs_env
1067 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1068 : LOGICAL :: occ, vir
1069 : LOGICAL, OPTIONAL :: can_skip
1070 :
1071 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_to_Sigma'
1072 :
1073 : INTEGER :: handle, inner_loop_atoms_interval_index
1074 : INTEGER(KIND=int_8) :: flop
1075 : INTEGER, DIMENSION(2) :: bounds_lambda, bounds_mu, bounds_nu, &
1076 : bounds_sigma, IL_atoms
1077 : INTEGER, DIMENSION(2, 2) :: bounds_comb
1078 : REAL(KIND=dp) :: sign_Sigma
1079 73850 : TYPE(dbt_type) :: t_3c_for_G, t_3c_x_G, t_3c_x_G_2
1080 :
1081 2954 : CALL timeset(routineN, handle)
1082 :
1083 2954 : CPASSERT(occ .EQV. (.NOT. vir))
1084 2954 : IF (occ) sign_Sigma = -1.0_dp
1085 2954 : IF (vir) sign_Sigma = 1.0_dp
1086 :
1087 2954 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_G)
1088 2954 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_G)
1089 2954 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_G_2)
1090 :
1091 : ! Here, in the first step e.g., is computed: N_Qλν = sum_µ (Qλ|µ) G_νµ
1092 : ! Afterwards e.g., is computed: Σ_λσ = sum_νQ M_Qνσ N_Qνλ (after reordering)
1093 : ! Bounds:
1094 : ! "comb" (combined index)
1095 : ! -> Q bounds from j_atoms and sparse in λ
1096 : ! -> λ (AO_1 in compute_3c_integrals) sparse in Q and µ
1097 : ! µ (AO_2 in compute_3c_integrals) bounds from inner loop "IL" indices and sparse in Q and λ
1098 : ! ν bounds from i_atoms
1099 : ! σ sparse in ν
1100 :
1101 : ! ν
1102 : bounds_nu(1:2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
1103 8862 : bs_env%i_ao_end_from_atom(i_atoms(2))]
1104 :
1105 5908 : DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
1106 8862 : IL_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
1107 :
1108 : ! µ
1109 : CALL get_bounds_from_atoms(bounds_mu, j_atoms, [1, bs_env%n_atom], &
1110 : bs_env%min_AO_idx_from_RI_AO_atom, &
1111 : bs_env%max_AO_idx_from_RI_AO_atom, &
1112 : atoms_3=IL_atoms, &
1113 : indices_3_start=bs_env%i_ao_start_from_atom, &
1114 8862 : indices_3_end=bs_env%i_ao_end_from_atom)
1115 :
1116 : ! Q
1117 : CALL get_bounds_from_atoms(bounds_comb(:, 1), IL_atoms, [1, bs_env%n_atom], &
1118 : bs_env%min_RI_idx_from_AO_AO_atom, &
1119 : bs_env%max_RI_idx_from_AO_AO_atom, &
1120 : atoms_3=j_atoms, &
1121 : indices_3_start=bs_env%i_RI_start_from_atom, &
1122 8862 : indices_3_end=bs_env%i_RI_end_from_atom)
1123 :
1124 : ! λ
1125 : CALL get_bounds_from_atoms(bounds_comb(:, 2), j_atoms, IL_atoms, &
1126 : bs_env%min_AO_idx_from_RI_AO_atom, &
1127 2954 : bs_env%max_AO_idx_from_RI_AO_atom)
1128 :
1129 2954 : IF (bounds_mu(1) > bounds_mu(2) .OR. bounds_comb(1, 1) > bounds_comb(2, 1) .OR. &
1130 : bounds_comb(1, 2) > bounds_comb(2, 2)) THEN
1131 : CYCLE
1132 : END IF
1133 :
1134 : CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_G, &
1135 2954 : atoms_RI=j_atoms, atoms_AO_2=IL_atoms)
1136 :
1137 : CALL dbt_contract(alpha=1.0_dp, &
1138 : tensor_1=t_2c_G, &
1139 : tensor_2=t_3c_for_G, &
1140 : beta=1.0_dp, &
1141 : tensor_3=t_3c_x_G, &
1142 : contract_1=[2], notcontract_1=[1], map_1=[3], &
1143 : contract_2=[3], notcontract_2=[1, 2], map_2=[1, 2], &
1144 : bounds_1=bounds_mu, &
1145 : bounds_2=bounds_nu, &
1146 : bounds_3=bounds_comb, &
1147 : flop=flop, &
1148 : move_data=.FALSE., &
1149 5908 : filter_eps=bs_env%eps_filter)
1150 : END DO ! IL_atoms
1151 :
1152 : ! Reordering: N_Qλν -> N_Qνλ
1153 2954 : CALL dbt_copy(t_3c_x_G, t_3c_x_G_2, order=[1, 3, 2], move_data=.TRUE.)
1154 :
1155 : ! Here, the last contraction is done, e.g., Σ_λσ = sum_νQ M_Qνσ N_Qνλ
1156 : ! Bounds as above, new "comb" with upper ingredients
1157 : bounds_comb(1:2, 1) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
1158 8862 : bs_env%i_RI_end_from_atom(j_atoms(2))]
1159 8862 : bounds_comb(1:2, 2) = bounds_nu(1:2)
1160 :
1161 : CALL get_bounds_from_atoms(bounds_lambda, j_atoms, [1, bs_env%n_atom], &
1162 : bs_env%min_AO_idx_from_RI_AO_atom, &
1163 8862 : bs_env%max_AO_idx_from_RI_AO_atom)
1164 : CALL get_bounds_from_atoms(bounds_sigma, [1, bs_env%n_atom], i_atoms, &
1165 : bs_env%min_AO_idx_from_RI_AO_atom, &
1166 8862 : bs_env%max_AO_idx_from_RI_AO_atom)
1167 :
1168 2954 : IF (bounds_sigma(1) > bounds_sigma(2) .OR. bounds_lambda(1) > bounds_lambda(2)) THEN
1169 0 : flop = 0_int_8
1170 : ELSE
1171 : CALL dbt_contract(alpha=sign_Sigma, &
1172 : tensor_1=t_3c_x_W, &
1173 : tensor_2=t_3c_x_G_2, &
1174 : beta=1.0_dp, &
1175 : tensor_3=t_2c_Sigma, &
1176 : contract_1=[1, 2], notcontract_1=[3], map_1=[1], &
1177 : contract_2=[1, 2], notcontract_2=[3], map_2=[2], &
1178 : bounds_1=bounds_comb, &
1179 : bounds_2=bounds_sigma, &
1180 : bounds_3=bounds_lambda, &
1181 2954 : filter_eps=bs_env%eps_filter, move_data=.FALSE., flop=flop)
1182 : END IF
1183 :
1184 2954 : IF (PRESENT(can_skip)) THEN
1185 2868 : IF (flop == 0_int_8) can_skip = .TRUE.
1186 : END IF
1187 :
1188 2954 : CALL dbt_destroy(t_3c_for_G)
1189 2954 : CALL dbt_destroy(t_3c_x_G)
1190 2954 : CALL dbt_destroy(t_3c_x_G_2)
1191 :
1192 2954 : CALL timestop(handle)
1193 :
1194 2954 : END SUBROUTINE contract_to_Sigma
1195 :
1196 : ! **************************************************************************************************
1197 : !> \brief Copy the positive- and negative-time self-energy matrices to the FM time array.
1198 : !> \param fm_Sigma_c_Gamma_time allocated output indexed by time, time sector and spin
1199 : !> \param bs_env time grid, spin count and AO layout
1200 : !> \param mat_Sigma_pos_tau positive-time self-energy matrices
1201 : !> \param mat_Sigma_neg_tau negative-time self-energy matrices
1202 : ! **************************************************************************************************
1203 74 : SUBROUTINE fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
1204 : mat_Sigma_pos_tau, mat_Sigma_neg_tau)
1205 :
1206 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
1207 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1208 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_pos_tau, mat_Sigma_neg_tau
1209 :
1210 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_fm_Sigma_c_Gamma_time'
1211 :
1212 : INTEGER :: handle, i_t, ispin, pos_neg
1213 :
1214 74 : CALL timeset(routineN, handle)
1215 :
1216 3682 : ALLOCATE (fm_Sigma_c_Gamma_time(bs_env%num_time_freq_points, 2, bs_env%n_spin))
1217 160 : DO ispin = 1, bs_env%n_spin
1218 1724 : DO i_t = 1, bs_env%num_time_freq_points
1219 4692 : DO pos_neg = 1, 2
1220 : CALL cp_fm_create(fm_Sigma_c_Gamma_time(i_t, pos_neg, ispin), &
1221 4692 : bs_env%fm_s_Gamma%matrix_struct)
1222 : END DO
1223 : CALL copy_dbcsr_to_fm(mat_Sigma_pos_tau(i_t, ispin)%matrix, &
1224 1564 : fm_Sigma_c_Gamma_time(i_t, 1, ispin))
1225 : CALL copy_dbcsr_to_fm(mat_Sigma_neg_tau(i_t, ispin)%matrix, &
1226 1650 : fm_Sigma_c_Gamma_time(i_t, 2, ispin))
1227 : END DO
1228 : END DO
1229 :
1230 74 : CALL timestop(handle)
1231 :
1232 74 : END SUBROUTINE fill_fm_Sigma_c_Gamma_time
1233 :
1234 : ! **************************************************************************************************
1235 : !> \brief ...
1236 : !> \param bs_env ...
1237 : ! **************************************************************************************************
1238 82 : SUBROUTINE print_skipping(bs_env)
1239 :
1240 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1241 :
1242 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_skipping'
1243 :
1244 : INTEGER :: handle, n_pairs
1245 :
1246 82 : CALL timeset(routineN, handle)
1247 :
1248 82 : n_pairs = bs_env%n_intervals_i*bs_env%n_intervals_j*bs_env%n_spin
1249 :
1250 82 : CALL bs_env%para_env_tensor%sum(bs_env%n_skip_sigma)
1251 82 : CALL bs_env%para_env_tensor%sum(bs_env%n_skip_chi)
1252 82 : CALL bs_env%para_env_tensor%sum(n_pairs)
1253 :
1254 82 : IF (bs_env%unit_nr > 0) THEN
1255 : WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
1256 41 : 'Sparsity of Σ^c(iτ,k=0): Percentage of skipped atom pairs:', &
1257 82 : REAL(100*bs_env%n_skip_sigma, KIND=dp)/REAL(n_pairs, KIND=dp), ' %'
1258 : WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
1259 41 : 'Sparsity of χ(iτ,k=0): Percentage of skipped atom pairs:', &
1260 82 : REAL(100*bs_env%n_skip_chi, KIND=dp)/REAL(n_pairs, KIND=dp), ' %'
1261 : END IF
1262 :
1263 82 : CALL timestop(handle)
1264 :
1265 82 : END SUBROUTINE print_skipping
1266 :
1267 : ! **************************************************************************************************
1268 : !> \brief ...
1269 : !> \param t_2c_Gocc ...
1270 : !> \param t_2c_Gvir ...
1271 : !> \param t_2c_W ...
1272 : !> \param t_2c_Sigma_neg_tau ...
1273 : !> \param t_2c_Sigma_pos_tau ...
1274 : !> \param t_3c_x_W ...
1275 : !> \param fm_W_MIC_time ...
1276 : !> \param mat_Sigma_neg_tau ...
1277 : !> \param mat_Sigma_pos_tau ...
1278 : ! **************************************************************************************************
1279 82 : SUBROUTINE destroy_mat_Sigma_c(t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
1280 : t_2c_Sigma_pos_tau, t_3c_x_W, fm_W_MIC_time, &
1281 : mat_Sigma_neg_tau, mat_Sigma_pos_tau)
1282 :
1283 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_2c_W, &
1284 : t_2c_Sigma_neg_tau, &
1285 : t_2c_Sigma_pos_tau, t_3c_x_W
1286 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1287 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
1288 :
1289 : CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_mat_Sigma_c'
1290 :
1291 : INTEGER :: handle
1292 :
1293 82 : CALL timeset(routineN, handle)
1294 :
1295 82 : CALL dbt_destroy(t_2c_Gocc)
1296 82 : CALL dbt_destroy(t_2c_Gvir)
1297 82 : CALL dbt_destroy(t_2c_W)
1298 82 : CALL dbt_destroy(t_2c_Sigma_neg_tau)
1299 82 : CALL dbt_destroy(t_2c_Sigma_pos_tau)
1300 82 : CALL dbt_destroy(t_3c_x_W)
1301 82 : CALL cp_fm_release(fm_W_MIC_time)
1302 82 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_neg_tau)
1303 82 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_pos_tau)
1304 :
1305 82 : CALL timestop(handle)
1306 :
1307 82 : END SUBROUTINE destroy_mat_Sigma_c
1308 :
1309 : ! **************************************************************************************************
1310 : !> \brief ...
1311 : !> \param bs_env ...
1312 : ! **************************************************************************************************
1313 122 : SUBROUTINE delete_unnecessary_files(bs_env)
1314 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1315 :
1316 : CHARACTER(LEN=*), PARAMETER :: routineN = 'delete_unnecessary_files'
1317 :
1318 : CHARACTER(LEN=default_path_length) :: f_chi, f_W_t, prefix
1319 : INTEGER :: handle, i_t
1320 :
1321 122 : CALL timeset(routineN, handle)
1322 :
1323 122 : prefix = bs_env%prefix
1324 :
1325 2088 : DO i_t = 1, bs_env%num_time_freq_points
1326 :
1327 1966 : IF (i_t < 10) THEN
1328 1064 : WRITE (f_chi, '(3A,I1,A)') TRIM(prefix), bs_env%chi_name, "_0", i_t, ".matrix"
1329 1064 : WRITE (f_W_t, '(3A,I1,A)') TRIM(prefix), bs_env%W_time_name, "_0", i_t, ".matrix"
1330 902 : ELSE IF (i_t < 100) THEN
1331 902 : WRITE (f_chi, '(3A,I2,A)') TRIM(prefix), bs_env%chi_name, "_", i_t, ".matrix"
1332 902 : WRITE (f_W_t, '(3A,I2,A)') TRIM(prefix), bs_env%W_time_name, "_", i_t, ".matrix"
1333 : ELSE
1334 0 : CPABORT('Please implement more than 99 time/frequency points.')
1335 : END IF
1336 :
1337 1966 : CALL safe_delete(f_chi, bs_env)
1338 2088 : CALL safe_delete(f_W_t, bs_env)
1339 :
1340 : END DO
1341 :
1342 122 : CALL timestop(handle)
1343 :
1344 122 : END SUBROUTINE delete_unnecessary_files
1345 :
1346 : ! **************************************************************************************************
1347 : !> \brief ...
1348 : !> \param filename ...
1349 : !> \param bs_env ...
1350 : ! **************************************************************************************************
1351 3932 : SUBROUTINE safe_delete(filename, bs_env)
1352 : CHARACTER(LEN=*) :: filename
1353 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1354 :
1355 : CHARACTER(LEN=*), PARAMETER :: routineN = 'safe_delete'
1356 :
1357 : INTEGER :: handle
1358 : LOGICAL :: file_exists
1359 :
1360 3932 : CALL timeset(routineN, handle)
1361 :
1362 3932 : IF (bs_env%para_env%mepos == 0) THEN
1363 :
1364 1966 : INQUIRE (file=TRIM(filename), exist=file_exists)
1365 1966 : IF (file_exists) CALL mp_file_delete(TRIM(filename))
1366 :
1367 : END IF
1368 :
1369 3932 : CALL timestop(handle)
1370 :
1371 3932 : END SUBROUTINE safe_delete
1372 :
1373 : ! **************************************************************************************************
1374 : !> \brief ...
1375 : !> \param bs_env ...
1376 : !> \param qs_env ...
1377 : !> \param fm_Sigma_x_Gamma ...
1378 : !> \param fm_Sigma_c_Gamma_time ...
1379 : ! **************************************************************************************************
1380 74 : SUBROUTINE compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
1381 :
1382 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1383 : TYPE(qs_environment_type), POINTER :: qs_env
1384 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
1385 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
1386 :
1387 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies'
1388 :
1389 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Sigma_c_ikp_n_freq
1390 : INTEGER :: handle, ikp, ispin, j_t
1391 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Sigma_x_ikp_n, V_xc_ikp_n
1392 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Sigma_c_ikp_n_time
1393 : TYPE(cp_cfm_type) :: cfm_ks_ikp, cfm_mos_ikp, cfm_s_ikp, &
1394 : cfm_work_ikp
1395 :
1396 74 : CALL timeset(routineN, handle)
1397 :
1398 74 : CALL cp_cfm_create(cfm_mos_ikp, bs_env%fm_s_Gamma%matrix_struct)
1399 74 : CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_s_Gamma%matrix_struct)
1400 : ! JW TODO: fully distribute these arrays at given time; also eigenvalues in bs_env
1401 296 : ALLOCATE (V_xc_ikp_n(bs_env%n_ao), Sigma_x_ikp_n(bs_env%n_ao))
1402 370 : ALLOCATE (Sigma_c_ikp_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
1403 296 : ALLOCATE (Sigma_c_ikp_n_freq(bs_env%n_ao, bs_env%num_freq_points_fit))
1404 :
1405 160 : DO ispin = 1, bs_env%n_spin
1406 :
1407 258 : DO ikp = 1, bs_env%nkp_bs_and_DOS
1408 :
1409 : ! 1. get H^KS_µν(k_i) from H^KS_µν(k=0)
1410 : CALL cfm_ikp_from_fm_Gamma(cfm_ks_ikp, bs_env%fm_ks_Gamma(ispin), &
1411 98 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1412 :
1413 : ! 2. get S_µν(k_i) from S_µν(k=0)
1414 : CALL cfm_ikp_from_fm_Gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
1415 98 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1416 :
1417 : ! 3. Diagonalize (Roothaan-Hall): H_KS(k_i)*C(k_i) = S(k_i)*C(k_i)*ϵ(k_i)
1418 : CALL cp_cfm_geeig(cfm_ks_ikp, cfm_s_ikp, cfm_mos_ikp, &
1419 98 : bs_env%eigenval_scf(:, ikp, ispin), cfm_work_ikp)
1420 :
1421 : ! 4. V^xc_µν(k=0) -> V^xc_µν(k_i) -> V^xc_nn(k_i)
1422 : CALL to_ikp_and_mo(V_xc_ikp_n, bs_env%fm_V_xc_Gamma(ispin), &
1423 98 : ikp, qs_env, bs_env, cfm_mos_ikp, cfm_ks_ikp, cfm_work_ikp)
1424 :
1425 : ! 5. Σ^x_µν(k=0) -> Σ^x_µν(k_i) -> Σ^x_nn(k_i)
1426 : CALL to_ikp_and_mo(Sigma_x_ikp_n, fm_Sigma_x_Gamma(ispin), &
1427 98 : ikp, qs_env, bs_env, cfm_mos_ikp, cfm_ks_ikp, cfm_work_ikp)
1428 :
1429 : ! 6. Σ^c_µν(k=0,+/-i|τ_j|) -> Σ^c_µν(k_i,+/-i|τ_j|) -> Σ^c_nn(k_i,+/-i|τ_j|)
1430 1766 : DO j_t = 1, bs_env%num_time_freq_points
1431 : CALL to_ikp_and_mo(Sigma_c_ikp_n_time(:, j_t, positive_time), &
1432 : fm_Sigma_c_Gamma_time(j_t, positive_time, ispin), &
1433 1668 : ikp, qs_env, bs_env, cfm_mos_ikp, cfm_ks_ikp, cfm_work_ikp)
1434 : CALL to_ikp_and_mo(Sigma_c_ikp_n_time(:, j_t, negative_time), &
1435 : fm_Sigma_c_Gamma_time(j_t, negative_time, ispin), &
1436 1766 : ikp, qs_env, bs_env, cfm_mos_ikp, cfm_ks_ikp, cfm_work_ikp)
1437 : END DO
1438 :
1439 : ! 7. Σ^c_nn(k_i,iτ) -> Σ^c_nn(k_i,iω)
1440 98 : CALL transform_time_to_frequency(bs_env%time_frequency_grid, Sigma_c_ikp_n_time, Sigma_c_ikp_n_freq)
1441 1580 : Sigma_c_ikp_n_freq(1:bs_env%n_occ(ispin), :) = CONJG(Sigma_c_ikp_n_freq(1:bs_env%n_occ(ispin), :))
1442 :
1443 : ! 8. Analytic continuation Σ^c_nn(k_i,iω) -> Σ^c_nn(k_i,ϵ) and
1444 : ! ϵ_nk_i^GW = ϵ_nk_i^DFT + Σ^c_nn(k_i,ϵ) + Σ^x_nn(k_i) - v^xc_nn(k_i)
1445 : CALL analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_xc_ikp_n, &
1446 184 : bs_env%eigenval_scf(:, ikp, ispin), ikp, ispin)
1447 :
1448 : END DO ! ikp_DOS
1449 :
1450 : END DO ! ispin
1451 :
1452 74 : CALL get_all_VBM_CBM_bandgaps(bs_env)
1453 :
1454 : ! Σ^x is releases here in case of G0W0
1455 74 : IF (bs_env%gw_flavour == G0W0) CALL cp_fm_release(fm_Sigma_x_Gamma)
1456 74 : CALL cp_fm_release(fm_Sigma_c_Gamma_time)
1457 74 : CALL cp_cfm_release(cfm_ks_ikp)
1458 74 : CALL cp_cfm_release(cfm_s_ikp)
1459 74 : CALL cp_cfm_release(cfm_mos_ikp)
1460 74 : CALL cp_cfm_release(cfm_work_ikp)
1461 :
1462 74 : CALL timestop(handle)
1463 :
1464 148 : END SUBROUTINE compute_QP_energies
1465 :
1466 : ! **************************************************************************************************
1467 : !> \brief Reconstruct a k-point Hermitian operator and return its real MO diagonal.
1468 : !> \param array_ikp_n Real diagonal elements in the MO basis.
1469 : !> \param fm_Gamma Real AO operator at the Gamma point.
1470 : !> \param ikp Index of the target k-point.
1471 : !> \param qs_env Quantum simulation environment.
1472 : !> \param bs_env Band-structure environment containing the target k-points.
1473 : !> \param cfm_mos_ikp Complex MO coefficients at the target k-point.
1474 : !> \param cfm_ikp Reusable AO operator buffer; overwritten by its MO projection.
1475 : !> \param workspace Reusable projection scratch matrix with the MO coefficient layout.
1476 : ! **************************************************************************************************
1477 3532 : SUBROUTINE to_ikp_and_mo(array_ikp_n, fm_Gamma, ikp, qs_env, bs_env, cfm_mos_ikp, cfm_ikp, workspace)
1478 :
1479 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: array_ikp_n
1480 : TYPE(cp_fm_type), INTENT(IN) :: fm_Gamma
1481 : INTEGER, INTENT(IN) :: ikp
1482 : TYPE(qs_environment_type), POINTER :: qs_env
1483 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1484 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_mos_ikp
1485 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_ikp, workspace
1486 :
1487 : CHARACTER(LEN=*), PARAMETER :: routineN = 'to_ikp_and_mo'
1488 :
1489 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cdiag
1490 : INTEGER :: handle
1491 :
1492 3532 : CALL timeset(routineN, handle)
1493 :
1494 : ! get cfm_µν(k_i) from fm_µν(k=0)
1495 3532 : CALL cfm_ikp_from_fm_Gamma(cfm_ikp, fm_Gamma, ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1496 :
1497 3532 : CALL cfm_contract_ABA(cfm_mos_ikp, cfm_ikp, workspace=workspace)
1498 :
1499 : ! Taking the real MO diagonal selects the Hermitian part of the operator.
1500 10596 : ALLOCATE (cdiag(SIZE(array_ikp_n)))
1501 3532 : CALL cp_cfm_get_diag(cfm_ikp, cdiag)
1502 43776 : array_ikp_n = REAL(cdiag, KIND=dp)
1503 :
1504 3532 : CALL timestop(handle)
1505 :
1506 7064 : END SUBROUTINE to_ikp_and_mo
1507 :
1508 : ! **************************************************************************************************
1509 : !> \brief Computes bounds (AO or RI) for given atom intervals atoms_1 and atoms_2 from indices_min
1510 : !> and indices_max and returns them in bounds_out.
1511 : !> In case, atoms_3 and indices_3 are given, the bounds are computed as the intersection
1512 : !> \param bounds_out Bounds to be computed
1513 : !> \param atoms_1 First atom interval
1514 : !> \param atoms_2 Second atom interval
1515 : !> \param indices_min Minimum indices for each atom pair (typically from bs_env,
1516 : !> computed in get_i_j_atom_ranges in gw_utils.F, e.g. bs_env%min_RI_idx_from_AO_AO_atom)
1517 : !> \param indices_max Maximum indices for each atom pair (typically from bs_env,
1518 : !> computed in get_i_j_atom_ranges in gw_utils.F)
1519 : !> \param atoms_3 (Optional) Third atom interval for intersection
1520 : !> \param indices_3_start (Optional) Indices for third atom interval for intersection
1521 : !> \param indices_3_end (Optional) Indices for third atom interval for intersection
1522 : ! **************************************************************************************************
1523 31970 : SUBROUTINE get_bounds_from_atoms(bounds_out, atoms_1, atoms_2, indices_min, indices_max, &
1524 31970 : atoms_3, indices_3_start, indices_3_end)
1525 :
1526 : INTEGER, DIMENSION(2), INTENT(OUT) :: bounds_out
1527 : INTEGER, DIMENSION(2), INTENT(IN) :: atoms_1, atoms_2
1528 : INTEGER, DIMENSION(:, :) :: indices_min, indices_max
1529 : INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: atoms_3
1530 : INTEGER, DIMENSION(:), OPTIONAL :: indices_3_start, indices_3_end
1531 :
1532 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_bounds_from_atoms'
1533 :
1534 : INTEGER :: handle, i_at, j_at
1535 :
1536 31970 : CALL timeset(routineN, handle)
1537 31970 : bounds_out(1) = HUGE(0)
1538 31970 : bounds_out(2) = -1
1539 : !Loop over all atoms in the two intervals and find min/max indices
1540 100348 : DO i_at = atoms_1(1), atoms_1(2)
1541 250418 : DO j_at = atoms_2(1), atoms_2(2)
1542 150070 : bounds_out(1) = MIN(bounds_out(1), indices_min(i_at, j_at))
1543 218448 : bounds_out(2) = MAX(bounds_out(2), indices_max(i_at, j_at))
1544 : END DO
1545 : END DO
1546 :
1547 31970 : IF (PRESENT(atoms_3) .AND. PRESENT(indices_3_start) .AND. PRESENT(indices_3_end)) THEN
1548 15268 : bounds_out(1) = MAX(bounds_out(1), indices_3_start(atoms_3(1)))
1549 15268 : bounds_out(2) = MIN(bounds_out(2), indices_3_end(atoms_3(2)))
1550 : END IF
1551 :
1552 31970 : CALL timestop(handle)
1553 :
1554 31970 : END SUBROUTINE get_bounds_from_atoms
1555 :
1556 : END MODULE gw_tensor_large_cell_Gamma
|
