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_large_cell_gamma
16 : USE atomic_kind_types, ONLY: atomic_kind_type
17 : USE bibliography, ONLY: Graml2024,&
18 : cite_reference
19 : USE cell_types, ONLY: cell_type,&
20 : get_cell,&
21 : pbc
22 : USE constants_operator, ONLY: operator_coulomb
23 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_uplo_to_full
24 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
25 : cp_cfm_cholesky_invert
26 : USE cp_cfm_diag, ONLY: cp_cfm_geeig
27 : USE cp_cfm_types, ONLY: cp_cfm_create,&
28 : cp_cfm_get_info,&
29 : cp_cfm_release,&
30 : cp_cfm_to_cfm,&
31 : cp_cfm_to_fm,&
32 : cp_cfm_type,&
33 : cp_fm_to_cfm
34 : USE cp_dbcsr_api, ONLY: &
35 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_block_p, &
36 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
37 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_release, dbcsr_set, &
38 : dbcsr_type
39 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
40 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
41 : copy_fm_to_dbcsr,&
42 : dbcsr_deallocate_matrix_set
43 : USE cp_files, ONLY: close_file,&
44 : open_file
45 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
46 : USE cp_fm_diag, ONLY: cp_fm_power
47 : USE cp_fm_types, ONLY: &
48 : cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_read_unformatted, cp_fm_release, &
49 : cp_fm_set_all, cp_fm_to_fm, cp_fm_type, cp_fm_write_unformatted
50 : USE cp_log_handling, ONLY: cp_get_default_logger,&
51 : cp_logger_type
52 : USE cp_output_handling, ONLY: cp_p_file,&
53 : cp_print_key_should_output,&
54 : cp_print_key_unit_nr
55 : USE dbt_api, ONLY: dbt_clear,&
56 : dbt_contract,&
57 : dbt_copy,&
58 : dbt_create,&
59 : dbt_destroy,&
60 : dbt_filter,&
61 : dbt_type
62 : USE gw_communication, ONLY: fm_to_local_tensor,&
63 : local_dbt_to_global_mat
64 : USE gw_utils, ONLY: analyt_conti_and_print,&
65 : de_init_bs_env,&
66 : time_to_freq
67 : USE input_constants, ONLY: G0W0,&
68 : rtp_method_bse,&
69 : rtp_method_bse_linearized
70 : USE input_section_types, ONLY: section_vals_type
71 : USE kinds, ONLY: default_path_length,&
72 : dp,&
73 : int_8
74 : USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
75 : USE kpoint_types, ONLY: kpoint_type
76 : USE machine, ONLY: m_walltime
77 : USE mathconstants, ONLY: twopi,&
78 : z_one,&
79 : z_zero
80 : USE message_passing, ONLY: mp_file_delete
81 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat
82 : USE parallel_gemm_api, ONLY: parallel_gemm
83 : USE particle_types, ONLY: particle_type
84 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
85 : USE post_scf_bandstructure_utils, ONLY: MIC_contribution_from_ikp,&
86 : cfm_ikp_from_fm_Gamma,&
87 : get_all_VBM_CBM_bandgaps
88 : USE qs_environment_types, ONLY: get_qs_env,&
89 : qs_environment_type
90 : USE qs_kind_types, ONLY: qs_kind_type
91 : USE qs_tensors, ONLY: build_3c_integrals
92 : USE rpa_gw_kpoints_util, ONLY: cp_cfm_power
93 : #include "./base/base_uses.f90"
94 :
95 : IMPLICIT NONE
96 :
97 : PRIVATE
98 :
99 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_large_cell_gamma'
100 :
101 : PUBLIC :: gw_calc_large_cell_Gamma, &
102 : compute_3c_integrals, G_occ_vir, fm_read, write_matrix, &
103 : fill_fm_Sigma_c_Gamma_time, delete_unnecessary_files, &
104 : multiply_fm_W_MIC_time_with_Minv_Gamma, get_W_MIC, &
105 : create_fm_W_MIC_time, Fourier_transform_w_to_t, &
106 : compute_fm_chi_Gamma_freq, compute_QP_energies, fm_write
107 :
108 : CONTAINS
109 :
110 : ! **************************************************************************************************
111 : !> \brief Perform GW band structure calculation
112 : !> \param qs_env ...
113 : !> \param bs_env ...
114 : !> \par History
115 : !> * 07.2023 created [Jan Wilhelm]
116 : ! **************************************************************************************************
117 68 : SUBROUTINE gw_calc_large_cell_Gamma(qs_env, bs_env)
118 : TYPE(qs_environment_type), POINTER :: qs_env
119 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
120 :
121 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_large_cell_Gamma'
122 :
123 : INTEGER :: handle
124 68 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma, fm_W_MIC_time
125 68 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
126 :
127 68 : CALL timeset(routineN, handle)
128 :
129 68 : CALL cite_reference(Graml2024)
130 :
131 : ! G^occ_µλ(i|τ|,k=0) = sum_n^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
132 : ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
133 : ! χ_PQ(iτ,k=0) = sum_λν [sum_µ (µν|P) G^occ_µλ(i|τ|)] [sum_σ (σλ|Q) G^vir_σν(i|τ|)]
134 68 : CALL get_mat_chi_Gamma_tau(bs_env, qs_env, bs_env%mat_chi_Gamma_tau)
135 :
136 : ! χ_PQ(iτ,k=0) -> χ_PQ(iω,k) -> ε_PQ(iω,k) -> W_PQ(iω,k) -> W^MIC_PQ(iτ) -> M^-1*W^MIC*M^-1
137 68 : CALL get_W_MIC(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_MIC_time)
138 :
139 : ! 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>
140 : ! Σ^x_λσ(k=0) = sum_νQ [sum_P (νσ|P) V^trunc_PQ] [sum_µ (λµ|Q) D_µν)]
141 68 : CALL get_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
142 :
143 : ! Σ^c_λσ(iτ,k=0) = sum_νQ [sum_P (νσ|P) W^MIC_PQ(iτ)] [sum_µ (λµ|Q) G^occ_µν(i|τ|)], τ < 0
144 : ! Σ^c_λσ(iτ,k=0) = sum_νQ [sum_P (νσ|P) W^MIC_PQ(iτ)] [sum_µ (λµ|Q) G^vir_µν(i|τ|)], τ > 0
145 68 : CALL get_Sigma_c(bs_env, qs_env, fm_W_MIC_time, fm_Sigma_c_Gamma_time)
146 :
147 : ! Σ^c_λσ(iτ,k=0) -> Σ^c_nn(ϵ,k); ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ,k) + Σ^x_nn(k) - v^xc_nn(k)
148 68 : CALL compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
149 :
150 68 : CALL de_init_bs_env(qs_env, bs_env)
151 :
152 68 : CALL timestop(handle)
153 :
154 68 : END SUBROUTINE gw_calc_large_cell_Gamma
155 :
156 : ! **************************************************************************************************
157 : !> \brief ...
158 : !> \param bs_env ...
159 : !> \param qs_env ...
160 : !> \param mat_chi_Gamma_tau ...
161 : ! **************************************************************************************************
162 68 : SUBROUTINE get_mat_chi_Gamma_tau(bs_env, qs_env, mat_chi_Gamma_tau)
163 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
164 : TYPE(qs_environment_type), POINTER :: qs_env
165 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
166 :
167 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_chi_Gamma_tau'
168 :
169 : INTEGER :: handle, i_intval_idx, i_t, inner_loop_atoms_interval_index, ispin, j_intval_idx
170 : INTEGER(KIND=int_8) :: flop
171 : INTEGER, DIMENSION(2) :: bounds_P, bounds_Q, i_atoms, IL_atoms, &
172 : j_atoms
173 : INTEGER, DIMENSION(2, 2) :: bounds_comb
174 : LOGICAL :: dist_too_long_i, dist_too_long_j
175 : REAL(KIND=dp) :: t1, tau
176 1700 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
177 1156 : t_3c_for_Gvir, t_3c_x_Gocc, &
178 1156 : t_3c_x_Gocc_2, t_3c_x_Gvir, &
179 612 : t_3c_x_Gvir_2
180 :
181 68 : CALL timeset(routineN, handle)
182 :
183 1312 : DO i_t = 1, bs_env%num_time_freq_points
184 :
185 1244 : t1 = m_walltime()
186 :
187 1244 : IF (bs_env%read_chi(i_t)) THEN
188 :
189 0 : CALL fm_read(bs_env%fm_RI_RI, bs_env, bs_env%chi_name, i_t)
190 :
191 : CALL copy_fm_to_dbcsr(bs_env%fm_RI_RI, mat_chi_Gamma_tau(i_t)%matrix, &
192 0 : keep_sparsity=.FALSE.)
193 :
194 0 : IF (bs_env%unit_nr > 0) THEN
195 : WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
196 0 : 'Read χ(iτ,k=0) from file for time point ', i_t, ' /', &
197 0 : bs_env%num_time_freq_points, &
198 0 : ', Execution time', m_walltime() - t1, ' s'
199 : END IF
200 :
201 : CYCLE
202 :
203 : END IF
204 :
205 1244 : IF (.NOT. bs_env%calc_chi(i_t)) CYCLE
206 :
207 : CALL create_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
208 1144 : t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2, bs_env)
209 :
210 : ! 1. compute G^occ and G^vir
211 : ! Background: G^σ(iτ) = G^occ,σ(iτ) * Θ(-τ) + G^vir,σ(iτ) * Θ(τ), σ ∈ {↑,↓}
212 : ! G^occ,σ_µλ(i|τ|,k=0) = sum_n^occ C^σ_µn(k=0) e^(-|(ϵ^σ_nk=0-ϵ_F)τ|) C^σ_λn(k=0)
213 : ! G^vir,σ_µλ(i|τ|,k=0) = sum_n^vir C^σ_µn(k=0) e^(-|(ϵ^σ_nk=0-ϵ_F)τ|) C^σ_λn(k=0)
214 1144 : tau = bs_env%imag_time_points(i_t)
215 :
216 2428 : DO ispin = 1, bs_env%n_spin
217 1284 : CALL G_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.TRUE., vir=.FALSE.)
218 1284 : CALL G_occ_vir(bs_env, tau, bs_env%fm_Gvir, ispin, occ=.FALSE., vir=.TRUE.)
219 :
220 : CALL fm_to_local_tensor(bs_env%fm_Gocc, bs_env%mat_ao_ao%matrix, &
221 : bs_env%mat_ao_ao_tensor%matrix, t_2c_Gocc, bs_env, &
222 1284 : bs_env%atoms_j_t_group)
223 : CALL fm_to_local_tensor(bs_env%fm_Gvir, bs_env%mat_ao_ao%matrix, &
224 : bs_env%mat_ao_ao_tensor%matrix, t_2c_Gvir, bs_env, &
225 1284 : bs_env%atoms_i_t_group)
226 :
227 : ! every group has its own range of i_atoms and j_atoms; only deal with a
228 : ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
229 3712 : DO i_intval_idx = 1, bs_env%n_intervals_i
230 3852 : DO j_intval_idx = 1, bs_env%n_intervals_j
231 3852 : i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
232 3852 : j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
233 :
234 1284 : IF (bs_env%skip_chi(i_intval_idx, j_intval_idx)) THEN
235 : ! Do that only after first timestep to avoid skips due to vanishing G
236 : ! caused by gaps
237 14 : IF (i_t == 2) THEN
238 0 : bs_env%n_skip_chi = bs_env%n_skip_chi + 1
239 : END IF
240 : CYCLE
241 : END IF
242 :
243 2540 : DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
244 :
245 3810 : IL_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
246 : ! Idea: Use sparsity in 3c integrals behind χ_PQ(iτ,k=0)
247 : ! -> λ bounds from j_atoms -> sparse in IL_atoms through σ in
248 : ! N_Qλν(iτ) = sum_σ (Qλ|σ) G^vir_νσ(i|τ|,k=0)
249 : ! -> ν bounds from i_atoms -> sparse in IL_atoms through µ in
250 : ! M_Pνλ(iτ) = sum_µ (Pν|µ) G^occ_λµ(i|τ|,k=0)
251 1270 : CALL check_dist(i_atoms, IL_atoms, qs_env, bs_env, dist_too_long_i)
252 1270 : CALL check_dist(j_atoms, IL_atoms, qs_env, bs_env, dist_too_long_j)
253 1270 : IF (.NOT. dist_too_long_i) THEN
254 : ! 2. compute 3-center integrals (Pν|µ) ("|": truncated Coulomb operator)
255 : CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gocc, &
256 1270 : atoms_AO_1=i_atoms, atoms_AO_2=IL_atoms)
257 : ! 3. tensor operation M_Pνλ(iτ) = sum_µ (Pν|µ) G^occ_λµ(i|τ|,k=0)
258 : CALL G_times_3c(t_3c_for_Gocc, t_2c_Gocc, t_3c_x_Gocc, bs_env, &
259 1270 : j_atoms, i_atoms, IL_atoms)
260 : END IF
261 2540 : IF (.NOT. dist_too_long_j) THEN
262 : ! 4. compute 3-center integrals (Qλ|σ) ("|": truncated Coulomb operator)
263 : CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gvir, &
264 1270 : atoms_AO_1=j_atoms, atoms_AO_2=IL_atoms)
265 : ! 5. tensor operation N_Qλν(iτ) = sum_σ (Qλ|σ) G^vir_νσ(i|τ|,k=0)
266 : CALL G_times_3c(t_3c_for_Gvir, t_2c_Gvir, t_3c_x_Gvir, bs_env, &
267 1270 : i_atoms, j_atoms, IL_atoms)
268 : END IF
269 : END DO ! IL_atoms
270 :
271 : ! 6. reorder tensors: M_Pνλ -> M_Pλν
272 1270 : CALL dbt_copy(t_3c_x_Gocc, t_3c_x_Gocc_2, move_data=.TRUE., order=[1, 3, 2])
273 1270 : CALL dbt_copy(t_3c_x_Gvir, t_3c_x_Gvir_2, move_data=.TRUE.)
274 :
275 : ! 7. tensor operation χ_PQ(iτ,k=0) = sum_λν M_Pλν(iτ) N_Qλν(iτ),
276 : ! Bounds:
277 : ! "comb" (combined index)
278 : ! -> λ bounds from j_atoms
279 : ! -> ν bounds from i_atoms
280 : ! P -> sparse in ν (see 3.)
281 : ! Q -> sparse in λ (see 5.)
282 : bounds_comb(1:2, 1) = [bs_env%i_ao_start_from_atom(j_atoms(1)), &
283 3810 : bs_env%i_ao_end_from_atom(j_atoms(2))]
284 : bounds_comb(1:2, 2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
285 3810 : bs_env%i_ao_end_from_atom(i_atoms(2))]
286 :
287 : CALL get_bounds_from_atoms(bounds_P, i_atoms, [1, bs_env%n_atom], &
288 : bs_env%min_RI_idx_from_AO_AO_atom, &
289 3810 : bs_env%max_RI_idx_from_AO_AO_atom)
290 : CALL get_bounds_from_atoms(bounds_Q, [1, bs_env%n_atom], j_atoms, &
291 : bs_env%min_RI_idx_from_AO_AO_atom, &
292 3810 : bs_env%max_RI_idx_from_AO_AO_atom)
293 :
294 1270 : IF (bounds_Q(1) > bounds_Q(2) .OR. bounds_P(1) > bounds_P(2)) THEN
295 0 : flop = 0_int_8
296 : ELSE
297 : CALL dbt_contract(alpha=bs_env%spin_degeneracy, &
298 : tensor_1=t_3c_x_Gocc_2, tensor_2=t_3c_x_Gvir_2, &
299 : beta=1.0_dp, tensor_3=bs_env%t_chi, &
300 : contract_1=[2, 3], notcontract_1=[1], map_1=[1], &
301 : contract_2=[2, 3], notcontract_2=[1], map_2=[2], &
302 : bounds_1=bounds_comb, &
303 : bounds_2=bounds_P, &
304 : bounds_3=bounds_Q, &
305 1270 : filter_eps=bs_env%eps_filter, move_data=.FALSE., flop=flop)
306 : END IF
307 2554 : IF (flop == 0_int_8) bs_env%skip_chi(i_intval_idx, j_intval_idx) = .TRUE.
308 :
309 : END DO ! j_atoms
310 : END DO ! i_atoms
311 : END DO ! ispin
312 :
313 : ! 8. communicate data of χ_PQ(iτ,k=0) in tensor bs_env%t_chi (which local in the
314 : ! subgroup) to the global dbcsr matrix mat_chi_Gamma_tau (which stores
315 : ! χ_PQ(iτ,k=0) for all time points)
316 : CALL local_dbt_to_global_mat(bs_env%t_chi, bs_env%mat_RI_RI_tensor%matrix, &
317 1144 : mat_chi_Gamma_tau(i_t)%matrix, bs_env%para_env)
318 :
319 : CALL write_matrix(mat_chi_Gamma_tau(i_t)%matrix, i_t, bs_env%chi_name, &
320 1144 : bs_env%fm_RI_RI, qs_env)
321 :
322 : CALL destroy_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
323 1144 : t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2)
324 :
325 1212 : IF (bs_env%unit_nr > 0) THEN
326 : WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F10.1,A)') &
327 572 : 'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
328 1144 : ', Execution time', m_walltime() - t1, ' s'
329 : END IF
330 :
331 : END DO ! i_t
332 :
333 68 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
334 :
335 68 : CALL timestop(handle)
336 :
337 68 : END SUBROUTINE get_mat_chi_Gamma_tau
338 :
339 : ! **************************************************************************************************
340 : !> \brief ...
341 : !> \param fm ...
342 : !> \param bs_env ...
343 : !> \param mat_name ...
344 : !> \param idx ...
345 : ! **************************************************************************************************
346 352 : SUBROUTINE fm_read(fm, bs_env, mat_name, idx)
347 : TYPE(cp_fm_type) :: fm
348 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
349 : CHARACTER(LEN=*) :: mat_name
350 : INTEGER :: idx
351 :
352 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_read'
353 :
354 : CHARACTER(LEN=default_path_length) :: f_chi
355 : INTEGER :: handle, unit_nr
356 :
357 352 : CALL timeset(routineN, handle)
358 :
359 352 : unit_nr = -1
360 352 : IF (bs_env%para_env%is_source()) THEN
361 :
362 176 : IF (idx < 10) THEN
363 87 : WRITE (f_chi, '(3A,I1,A)') TRIM(bs_env%prefix), TRIM(mat_name), "_0", idx, ".matrix"
364 89 : ELSE IF (idx < 100) THEN
365 89 : WRITE (f_chi, '(3A,I2,A)') TRIM(bs_env%prefix), TRIM(mat_name), "_", idx, ".matrix"
366 : ELSE
367 0 : CPABORT('Please implement more than 99 time/frequency points.')
368 : END IF
369 :
370 : CALL open_file(file_name=TRIM(f_chi), file_action="READ", file_form="UNFORMATTED", &
371 176 : file_position="REWIND", file_status="OLD", unit_number=unit_nr)
372 :
373 : END IF
374 :
375 352 : CALL cp_fm_read_unformatted(fm, unit_nr)
376 :
377 352 : IF (bs_env%para_env%is_source()) CALL close_file(unit_number=unit_nr)
378 :
379 352 : CALL timestop(handle)
380 :
381 352 : END SUBROUTINE fm_read
382 :
383 : ! **************************************************************************************************
384 : !> \brief ...
385 : !> \param t_2c_Gocc ...
386 : !> \param t_2c_Gvir ...
387 : !> \param t_3c_for_Gocc ...
388 : !> \param t_3c_for_Gvir ...
389 : !> \param t_3c_x_Gocc ...
390 : !> \param t_3c_x_Gvir ...
391 : !> \param t_3c_x_Gocc_2 ...
392 : !> \param t_3c_x_Gvir_2 ...
393 : !> \param bs_env ...
394 : ! **************************************************************************************************
395 1144 : SUBROUTINE create_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
396 : t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2, bs_env)
397 :
398 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
399 : t_3c_for_Gvir, t_3c_x_Gocc, &
400 : t_3c_x_Gvir, t_3c_x_Gocc_2, &
401 : t_3c_x_Gvir_2
402 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
403 :
404 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_tensors_chi'
405 :
406 : INTEGER :: handle
407 :
408 1144 : CALL timeset(routineN, handle)
409 :
410 1144 : CALL dbt_create(bs_env%t_G, t_2c_Gocc, name="Gocc 2c (AO|AO)")
411 1144 : CALL dbt_create(bs_env%t_G, t_2c_Gvir, name="Gvir 2c (AO|AO)")
412 1144 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gocc, name="Gocc 3c (RI AO|AO)")
413 1144 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gvir, name="Gvir 3c (RI AO|AO)")
414 1144 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gocc, name="xGocc 3c (RI AO|AO)")
415 1144 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gvir, name="xGvir 3c (RI AO|AO)")
416 1144 : CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gocc_2, name="x2Gocc 3c (RI AO|AO)")
417 1144 : CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gvir_2, name="x2Gvir 3c (RI AO|AO)")
418 :
419 1144 : CALL timestop(handle)
420 :
421 1144 : END SUBROUTINE create_tensors_chi
422 :
423 : ! **************************************************************************************************
424 : !> \brief ...
425 : !> \param t_2c_Gocc ...
426 : !> \param t_2c_Gvir ...
427 : !> \param t_3c_for_Gocc ...
428 : !> \param t_3c_for_Gvir ...
429 : !> \param t_3c_x_Gocc ...
430 : !> \param t_3c_x_Gvir ...
431 : !> \param t_3c_x_Gocc_2 ...
432 : !> \param t_3c_x_Gvir_2 ...
433 : ! **************************************************************************************************
434 1144 : SUBROUTINE destroy_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
435 : t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2)
436 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
437 : t_3c_for_Gvir, t_3c_x_Gocc, &
438 : t_3c_x_Gvir, t_3c_x_Gocc_2, &
439 : t_3c_x_Gvir_2
440 :
441 : CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_tensors_chi'
442 :
443 : INTEGER :: handle
444 :
445 1144 : CALL timeset(routineN, handle)
446 :
447 1144 : CALL dbt_destroy(t_2c_Gocc)
448 1144 : CALL dbt_destroy(t_2c_Gvir)
449 1144 : CALL dbt_destroy(t_3c_for_Gocc)
450 1144 : CALL dbt_destroy(t_3c_for_Gvir)
451 1144 : CALL dbt_destroy(t_3c_x_Gocc)
452 1144 : CALL dbt_destroy(t_3c_x_Gvir)
453 1144 : CALL dbt_destroy(t_3c_x_Gocc_2)
454 1144 : CALL dbt_destroy(t_3c_x_Gvir_2)
455 :
456 1144 : CALL timestop(handle)
457 :
458 1144 : END SUBROUTINE destroy_tensors_chi
459 :
460 : ! **************************************************************************************************
461 : !> \brief ...
462 : !> \param matrix ...
463 : !> \param matrix_index ...
464 : !> \param matrix_name ...
465 : !> \param fm ...
466 : !> \param qs_env ...
467 : ! **************************************************************************************************
468 3782 : SUBROUTINE write_matrix(matrix, matrix_index, matrix_name, fm, qs_env)
469 : TYPE(dbcsr_type) :: matrix
470 : INTEGER :: matrix_index
471 : CHARACTER(LEN=*) :: matrix_name
472 : TYPE(cp_fm_type), INTENT(IN), POINTER :: fm
473 : TYPE(qs_environment_type), POINTER :: qs_env
474 :
475 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_matrix'
476 :
477 : INTEGER :: handle
478 :
479 3782 : CALL timeset(routineN, handle)
480 :
481 3782 : CALL cp_fm_set_all(fm, 0.0_dp)
482 :
483 3782 : CALL copy_dbcsr_to_fm(matrix, fm)
484 :
485 3782 : CALL fm_write(fm, matrix_index, matrix_name, qs_env)
486 :
487 3782 : CALL timestop(handle)
488 :
489 3782 : END SUBROUTINE write_matrix
490 :
491 : ! **************************************************************************************************
492 : !> \brief ...
493 : !> \param fm ...
494 : !> \param matrix_index ...
495 : !> \param matrix_name ...
496 : !> \param qs_env ...
497 : ! **************************************************************************************************
498 4994 : SUBROUTINE fm_write(fm, matrix_index, matrix_name, qs_env)
499 : TYPE(cp_fm_type) :: fm
500 : INTEGER :: matrix_index
501 : CHARACTER(LEN=*) :: matrix_name
502 : TYPE(qs_environment_type), POINTER :: qs_env
503 :
504 : CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
505 : routineN = 'fm_write'
506 :
507 : CHARACTER(LEN=default_path_length) :: filename
508 : INTEGER :: handle, unit_nr
509 : TYPE(cp_logger_type), POINTER :: logger
510 : TYPE(section_vals_type), POINTER :: input
511 :
512 4994 : CALL timeset(routineN, handle)
513 :
514 4994 : CALL get_qs_env(qs_env, input=input)
515 :
516 4994 : logger => cp_get_default_logger()
517 :
518 4994 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
519 :
520 944 : IF (matrix_index < 10) THEN
521 456 : WRITE (filename, '(3A,I1)') "RESTART_", matrix_name, "_0", matrix_index
522 488 : ELSE IF (matrix_index < 100) THEN
523 488 : WRITE (filename, '(3A,I2)') "RESTART_", matrix_name, "_", matrix_index
524 : ELSE
525 0 : CPABORT('Please implement more than 99 time/frequency points.')
526 : END IF
527 :
528 : unit_nr = cp_print_key_unit_nr(logger, input, key, extension=".matrix", &
529 : file_form="UNFORMATTED", middle_name=TRIM(filename), &
530 944 : file_position="REWIND", file_action="WRITE")
531 :
532 944 : CALL cp_fm_write_unformatted(fm, unit_nr)
533 944 : IF (unit_nr > 0) THEN
534 472 : CALL close_file(unit_nr)
535 : END IF
536 : END IF
537 :
538 4994 : CALL timestop(handle)
539 :
540 4994 : END SUBROUTINE fm_write
541 :
542 : ! **************************************************************************************************
543 : !> \brief ...
544 : !> \param bs_env ...
545 : !> \param tau ...
546 : !> \param fm_G_Gamma ...
547 : !> \param ispin ...
548 : !> \param occ ...
549 : !> \param vir ...
550 : ! **************************************************************************************************
551 16016 : SUBROUTINE G_occ_vir(bs_env, tau, fm_G_Gamma, ispin, occ, vir)
552 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
553 : REAL(KIND=dp) :: tau
554 : TYPE(cp_fm_type) :: fm_G_Gamma
555 : INTEGER :: ispin
556 : LOGICAL :: occ, vir
557 :
558 : CHARACTER(LEN=*), PARAMETER :: routineN = 'G_occ_vir'
559 :
560 : INTEGER :: handle, homo, i_row_local, j_col, &
561 : j_col_local, n_mo, ncol_local, &
562 : nrow_local
563 8008 : INTEGER, DIMENSION(:), POINTER :: col_indices
564 : REAL(KIND=dp) :: tau_E
565 :
566 8008 : CALL timeset(routineN, handle)
567 :
568 8008 : CPASSERT(occ .NEQV. vir)
569 :
570 : CALL cp_fm_get_info(matrix=bs_env%fm_work_mo(1), &
571 : nrow_local=nrow_local, &
572 : ncol_local=ncol_local, &
573 8008 : col_indices=col_indices)
574 :
575 8008 : n_mo = bs_env%n_ao
576 8008 : homo = bs_env%n_occ(ispin)
577 :
578 8008 : CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(ispin), bs_env%fm_work_mo(1))
579 :
580 49274 : DO i_row_local = 1, nrow_local
581 523600 : DO j_col_local = 1, ncol_local
582 :
583 474326 : j_col = col_indices(j_col_local)
584 :
585 474326 : tau_E = ABS(tau*0.5_dp*(bs_env%eigenval_scf_Gamma(j_col, ispin) - bs_env%e_fermi(ispin)))
586 :
587 474326 : IF (tau_E < bs_env%stabilize_exp) THEN
588 : bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = &
589 472638 : bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local)*EXP(-tau_E)
590 : ELSE
591 1688 : bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = 0.0_dp
592 : END IF
593 :
594 515592 : IF ((occ .AND. j_col > homo) .OR. (vir .AND. j_col <= homo)) THEN
595 239310 : bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = 0.0_dp
596 : END IF
597 :
598 : END DO
599 : END DO
600 :
601 : CALL parallel_gemm(transa="N", transb="T", m=n_mo, n=n_mo, k=n_mo, alpha=1.0_dp, &
602 : matrix_a=bs_env%fm_work_mo(1), matrix_b=bs_env%fm_work_mo(1), &
603 8008 : beta=0.0_dp, matrix_c=fm_G_Gamma)
604 :
605 8008 : CALL timestop(handle)
606 :
607 8008 : END SUBROUTINE G_occ_vir
608 :
609 : ! **************************************************************************************************
610 : !> \brief ...
611 : !> \param qs_env ...
612 : !> \param bs_env ...
613 : !> \param t_3c ...
614 : !> \param atoms_AO_1 ...
615 : !> \param atoms_AO_2 ...
616 : !> \param atoms_RI ...
617 : ! **************************************************************************************************
618 6528 : SUBROUTINE compute_3c_integrals(qs_env, bs_env, t_3c, atoms_AO_1, atoms_AO_2, atoms_RI)
619 : TYPE(qs_environment_type), POINTER :: qs_env
620 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
621 : TYPE(dbt_type) :: t_3c
622 : INTEGER, DIMENSION(2), OPTIONAL :: atoms_AO_1, atoms_AO_2, atoms_RI
623 :
624 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_integrals'
625 :
626 : INTEGER :: handle
627 6528 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_array
628 :
629 6528 : CALL timeset(routineN, handle)
630 :
631 : ! free memory (not clear whether memory has been freed previously)
632 6528 : CALL dbt_clear(t_3c)
633 :
634 71808 : ALLOCATE (t_3c_array(1, 1))
635 6528 : CALL dbt_create(t_3c, t_3c_array(1, 1))
636 :
637 : CALL build_3c_integrals(t_3c_array, &
638 : bs_env%eps_filter, &
639 : qs_env, &
640 : bs_env%nl_3c, &
641 : int_eps=bs_env%eps_filter, &
642 : basis_i=bs_env%basis_set_RI, &
643 : basis_j=bs_env%basis_set_AO, &
644 : basis_k=bs_env%basis_set_AO, &
645 : potential_parameter=bs_env%ri_metric, &
646 : bounds_i=atoms_RI, &
647 : bounds_j=atoms_AO_1, &
648 : bounds_k=atoms_AO_2, &
649 6528 : desymmetrize=.FALSE.)
650 :
651 6528 : CALL dbt_filter(t_3c_array(1, 1), bs_env%eps_filter)
652 :
653 6528 : CALL dbt_copy(t_3c_array(1, 1), t_3c, move_data=.TRUE.)
654 :
655 6528 : CALL dbt_destroy(t_3c_array(1, 1))
656 13056 : DEALLOCATE (t_3c_array)
657 :
658 6528 : CALL timestop(handle)
659 :
660 13056 : END SUBROUTINE compute_3c_integrals
661 :
662 : ! **************************************************************************************************
663 : !> \brief ...
664 : !> \param t_3c_for_G ...
665 : !> \param t_G ...
666 : !> \param t_M ...
667 : !> \param bs_env ...
668 : !> \param atoms_AO_1 ...
669 : !> \param atoms_AO_2 ...
670 : !> \param atoms_IL ...
671 : ! **************************************************************************************************
672 2540 : SUBROUTINE G_times_3c(t_3c_for_G, t_G, t_M, bs_env, atoms_AO_1, atoms_AO_2, atoms_IL)
673 : TYPE(dbt_type) :: t_3c_for_G, t_G, t_M
674 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
675 : INTEGER, DIMENSION(2) :: atoms_AO_1, atoms_AO_2, atoms_IL
676 :
677 : CHARACTER(LEN=*), PARAMETER :: routineN = 'G_times_3c'
678 :
679 : INTEGER :: handle
680 : INTEGER(KIND=int_8) :: flop
681 : INTEGER, DIMENSION(2) :: bounds_ao_1, bounds_IL
682 : INTEGER, DIMENSION(2, 2) :: bounds_comb
683 :
684 2540 : CALL timeset(routineN, handle)
685 :
686 : ! Bounds reduce needed memory and therefore scaling behavior
687 : ! Operations are of the form, e.g, M_Pνλ = sum_µ (Pν|µ) G_λµ
688 : ! "comb" (combined index)
689 : ! -> P sparse in ν and µ
690 : ! -> λ bounds from j_atoms (via atoms_AO_1)
691 : ! µ bounds from inner loop "IL" indices and sparse in P and ν
692 : ! ν bounds from i_atoms (via atoms_AO_2) and sparse in P and µ
693 :
694 : ! µ index
695 : CALL get_bounds_from_atoms(bounds_IL, [1, bs_env%n_atom], atoms_AO_2, &
696 : bs_env%min_AO_idx_from_RI_AO_atom, &
697 : bs_env%max_AO_idx_from_RI_AO_atom, &
698 : atoms_3=atoms_IL, &
699 : indices_3_start=bs_env%i_ao_start_from_atom, &
700 7620 : indices_3_end=bs_env%i_ao_end_from_atom)
701 :
702 : ! P index
703 : CALL get_bounds_from_atoms(bounds_comb(:, 1), atoms_IL, atoms_AO_2, &
704 : bs_env%min_RI_idx_from_AO_AO_atom, &
705 2540 : bs_env%max_RI_idx_from_AO_AO_atom)
706 :
707 : ! ν index
708 : CALL get_bounds_from_atoms(bounds_comb(:, 2), [1, bs_env%n_atom], atoms_IL, &
709 : bs_env%min_AO_idx_from_RI_AO_atom, &
710 : bs_env%max_AO_idx_from_RI_AO_atom, &
711 : atoms_3=atoms_AO_2, &
712 : indices_3_start=bs_env%i_ao_start_from_atom, &
713 7620 : indices_3_end=bs_env%i_ao_end_from_atom)
714 :
715 : ! λ index
716 : bounds_ao_1(1:2) = [bs_env%i_ao_start_from_atom(atoms_AO_1(1)), &
717 7620 : bs_env%i_ao_end_from_atom(atoms_AO_1(2))]
718 :
719 2540 : IF (bounds_IL(1) > bounds_IL(2) .OR. bounds_comb(1, 2) > bounds_comb(2, 2)) THEN
720 : flop = 0_int_8
721 : ELSE
722 : CALL dbt_contract(alpha=1.0_dp, &
723 : tensor_1=t_3c_for_G, &
724 : tensor_2=t_G, &
725 : beta=1.0_dp, &
726 : tensor_3=t_M, &
727 : contract_1=[3], notcontract_1=[1, 2], map_1=[1, 2], &
728 : contract_2=[2], notcontract_2=[1], map_2=[3], &
729 : bounds_1=bounds_IL, &
730 : bounds_2=bounds_comb, &
731 : bounds_3=bounds_ao_1, &
732 : flop=flop, &
733 2540 : filter_eps=bs_env%eps_filter)
734 : END IF
735 :
736 2540 : CALL dbt_clear(t_3c_for_G)
737 :
738 2540 : CALL timestop(handle)
739 :
740 2540 : END SUBROUTINE G_times_3c
741 :
742 : ! **************************************************************************************************
743 : !> \brief ...
744 : !> \param atoms_1 ...
745 : !> \param atoms_2 ...
746 : !> \param qs_env ...
747 : !> \param bs_env ...
748 : !> \param dist_too_long ...
749 : ! **************************************************************************************************
750 2540 : SUBROUTINE check_dist(atoms_1, atoms_2, qs_env, bs_env, dist_too_long)
751 : INTEGER, DIMENSION(2) :: atoms_1, atoms_2
752 : TYPE(qs_environment_type), POINTER :: qs_env
753 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
754 : LOGICAL :: dist_too_long
755 :
756 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_dist'
757 :
758 : INTEGER :: atom_1, atom_2, handle
759 : REAL(dp) :: abs_rab, min_dist_AO_atoms
760 : REAL(KIND=dp), DIMENSION(3) :: rab
761 : TYPE(cell_type), POINTER :: cell
762 2540 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
763 :
764 2540 : CALL timeset(routineN, handle)
765 :
766 2540 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
767 :
768 2540 : min_dist_AO_atoms = HUGE(1.0_dp)
769 7668 : DO atom_1 = atoms_1(1), atoms_1(2)
770 18068 : DO atom_2 = atoms_2(1), atoms_2(2)
771 10400 : rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
772 :
773 10400 : abs_rab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
774 :
775 15528 : min_dist_AO_atoms = MIN(min_dist_AO_atoms, abs_rab)
776 : END DO
777 : END DO
778 :
779 2540 : dist_too_long = (min_dist_AO_atoms > bs_env%max_dist_AO_atoms)
780 :
781 2540 : CALL timestop(handle)
782 :
783 2540 : END SUBROUTINE check_dist
784 :
785 : ! **************************************************************************************************
786 : !> \brief ...
787 : !> \param bs_env ...
788 : !> \param qs_env ...
789 : !> \param mat_chi_Gamma_tau ...
790 : !> \param fm_W_MIC_time ...
791 : ! **************************************************************************************************
792 68 : SUBROUTINE get_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
793 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
794 : TYPE(qs_environment_type), POINTER :: qs_env
795 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
796 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
797 :
798 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_W_MIC'
799 :
800 : INTEGER :: handle
801 :
802 68 : CALL timeset(routineN, handle)
803 :
804 68 : IF (bs_env%all_W_exist) THEN
805 6 : CALL read_W_MIC_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
806 : ELSE
807 62 : CALL compute_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
808 : END IF
809 :
810 68 : CALL timestop(handle)
811 :
812 68 : END SUBROUTINE get_W_MIC
813 :
814 : ! **************************************************************************************************
815 : !> \brief ...
816 : !> \param bs_env ...
817 : !> \param qs_env ...
818 : !> \param fm_V_kp ...
819 : !> \param ikp_batch ...
820 : ! **************************************************************************************************
821 110 : SUBROUTINE compute_V_k_by_lattice_sum(bs_env, qs_env, fm_V_kp, ikp_batch)
822 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
823 : TYPE(qs_environment_type), POINTER :: qs_env
824 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_V_kp
825 : INTEGER :: ikp_batch
826 :
827 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_k_by_lattice_sum'
828 :
829 : INTEGER :: handle, ikp, ikp_end, ikp_start, &
830 : nkp_chi_eps_W_batch, re_im
831 110 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
832 : TYPE(cell_type), POINTER :: cell
833 110 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_kp
834 110 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
835 110 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
836 :
837 110 : CALL timeset(routineN, handle)
838 :
839 110 : nkp_chi_eps_W_batch = bs_env%nkp_chi_eps_W_batch
840 :
841 110 : ikp_start = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + 1
842 110 : ikp_end = MIN(ikp_batch*bs_env%nkp_chi_eps_W_batch, bs_env%kpoints_chi_eps_W%nkp)
843 :
844 110 : NULLIFY (mat_V_kp)
845 1184 : ALLOCATE (mat_V_kp(ikp_start:ikp_end, 2))
846 :
847 330 : DO re_im = 1, 2
848 854 : DO ikp = ikp_start, ikp_end
849 524 : NULLIFY (mat_V_kp(ikp, re_im)%matrix)
850 524 : ALLOCATE (mat_V_kp(ikp, re_im)%matrix)
851 524 : CALL dbcsr_create(mat_V_kp(ikp, re_im)%matrix, template=bs_env%mat_RI_RI%matrix)
852 524 : CALL dbcsr_reserve_all_blocks(mat_V_kp(ikp, re_im)%matrix)
853 744 : CALL dbcsr_set(mat_V_kp(ikp, re_im)%matrix, 0.0_dp)
854 : END DO ! ikp
855 : END DO ! re_im
856 :
857 : CALL get_qs_env(qs_env=qs_env, &
858 : particle_set=particle_set, &
859 : cell=cell, &
860 : qs_kind_set=qs_kind_set, &
861 110 : atomic_kind_set=atomic_kind_set)
862 :
863 110 : IF (ikp_end <= bs_env%nkp_chi_eps_W_orig) THEN
864 :
865 : ! 1. 2c Coulomb integrals for the first "original" k-point grid
866 280 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
867 :
868 40 : ELSE IF (ikp_start > bs_env%nkp_chi_eps_W_orig .AND. &
869 : ikp_end <= bs_env%nkp_chi_eps_W_orig_plus_extra) THEN
870 :
871 : ! 2. 2c Coulomb integrals for the second "extrapolation" k-point grid
872 160 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_extra
873 :
874 : ELSE
875 :
876 0 : CPABORT("Error with k-point parallelization.")
877 :
878 : END IF
879 :
880 : CALL build_2c_coulomb_matrix_kp(mat_V_kp, &
881 : bs_env%kpoints_chi_eps_W, &
882 : basis_type="RI_AUX", &
883 : cell=cell, &
884 : particle_set=particle_set, &
885 : qs_kind_set=qs_kind_set, &
886 : atomic_kind_set=atomic_kind_set, &
887 : size_lattice_sum=bs_env%size_lattice_sum_V, &
888 : operator_type=operator_coulomb, &
889 : ikp_start=ikp_start, &
890 110 : ikp_end=ikp_end)
891 :
892 440 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
893 :
894 1184 : ALLOCATE (fm_V_kp(ikp_start:ikp_end, 2))
895 330 : DO re_im = 1, 2
896 854 : DO ikp = ikp_start, ikp_end
897 524 : CALL cp_fm_create(fm_V_kp(ikp, re_im), bs_env%fm_RI_RI%matrix_struct)
898 524 : CALL copy_dbcsr_to_fm(mat_V_kp(ikp, re_im)%matrix, fm_V_kp(ikp, re_im))
899 744 : CALL dbcsr_deallocate_matrix(mat_V_kp(ikp, re_im)%matrix)
900 : END DO
901 : END DO
902 110 : DEALLOCATE (mat_V_kp)
903 :
904 110 : CALL timestop(handle)
905 :
906 110 : END SUBROUTINE compute_V_k_by_lattice_sum
907 :
908 : ! **************************************************************************************************
909 : !> \brief ...
910 : !> \param bs_env ...
911 : !> \param qs_env ...
912 : !> \param fm_V_kp ...
913 : !> \param cfm_V_sqrt_ikp ...
914 : !> \param cfm_M_inv_V_sqrt_ikp ...
915 : !> \param ikp ...
916 : ! **************************************************************************************************
917 262 : SUBROUTINE compute_MinvVsqrt_Vsqrt(bs_env, qs_env, fm_V_kp, cfm_V_sqrt_ikp, &
918 : cfm_M_inv_V_sqrt_ikp, ikp)
919 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
920 : TYPE(qs_environment_type), POINTER :: qs_env
921 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_V_kp
922 : TYPE(cp_cfm_type) :: cfm_V_sqrt_ikp, cfm_M_inv_V_sqrt_ikp
923 : INTEGER :: ikp
924 :
925 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_MinvVsqrt_Vsqrt'
926 :
927 : INTEGER :: handle, info, n_RI
928 : TYPE(cp_cfm_type) :: cfm_M_inv_ikp, cfm_work
929 262 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_M_ikp
930 :
931 262 : CALL timeset(routineN, handle)
932 :
933 262 : n_RI = bs_env%n_RI
934 :
935 : ! get here M(k) and write it to fm_M_ikp
936 : CALL RI_2c_integral_mat(qs_env, fm_M_ikp, fm_V_kp(ikp, 1), &
937 : n_RI, bs_env%ri_metric, do_kpoints=.TRUE., &
938 : kpoints=bs_env%kpoints_chi_eps_W, &
939 : regularization_RI=bs_env%regularization_RI, ikp_ext=ikp, &
940 262 : do_build_cell_index=(ikp == 1))
941 :
942 262 : IF (ikp == 1) THEN
943 62 : CALL cp_cfm_create(cfm_V_sqrt_ikp, fm_V_kp(ikp, 1)%matrix_struct)
944 62 : CALL cp_cfm_create(cfm_M_inv_V_sqrt_ikp, fm_V_kp(ikp, 1)%matrix_struct)
945 : END IF
946 262 : CALL cp_cfm_create(cfm_M_inv_ikp, fm_V_kp(ikp, 1)%matrix_struct)
947 :
948 262 : CALL cp_fm_to_cfm(fm_M_ikp(1, 1), fm_M_ikp(1, 2), cfm_M_inv_ikp)
949 262 : CALL cp_fm_to_cfm(fm_V_kp(ikp, 1), fm_V_kp(ikp, 2), cfm_V_sqrt_ikp)
950 :
951 262 : CALL cp_fm_release(fm_M_ikp)
952 :
953 262 : CALL cp_cfm_create(cfm_work, fm_V_kp(ikp, 1)%matrix_struct)
954 :
955 : ! M(k) -> M^-1(k)
956 262 : CALL cp_cfm_to_cfm(cfm_M_inv_ikp, cfm_work)
957 262 : CALL cp_cfm_cholesky_decompose(matrix=cfm_M_inv_ikp, n=n_RI, info_out=info)
958 262 : IF (info == 0) THEN
959 : ! successful Cholesky decomposition
960 262 : CALL cp_cfm_cholesky_invert(cfm_M_inv_ikp)
961 : ! symmetrize the result
962 262 : CALL cp_cfm_uplo_to_full(cfm_M_inv_ikp)
963 : ELSE
964 : ! Cholesky decomposition not successful: use expensive diagonalization
965 0 : CALL cp_cfm_power(cfm_work, threshold=bs_env%eps_eigval_mat_RI, exponent=-1.0_dp)
966 0 : CALL cp_cfm_to_cfm(cfm_work, cfm_M_inv_ikp)
967 : END IF
968 :
969 : ! V(k) -> L(k) with L^H(k)*L(k) = V(k) [L(k) can be just considered to be V^0.5(k)]
970 262 : CALL cp_cfm_to_cfm(cfm_V_sqrt_ikp, cfm_work)
971 262 : CALL cp_cfm_cholesky_decompose(matrix=cfm_V_sqrt_ikp, n=n_RI, info_out=info)
972 262 : IF (info == 0) THEN
973 : ! successful Cholesky decomposition
974 262 : CALL clean_lower_part(cfm_V_sqrt_ikp)
975 : ELSE
976 : ! Cholesky decomposition not successful: use expensive diagonalization
977 0 : CALL cp_cfm_power(cfm_work, threshold=0.0_dp, exponent=0.5_dp)
978 0 : CALL cp_cfm_to_cfm(cfm_work, cfm_V_sqrt_ikp)
979 : END IF
980 262 : CALL cp_cfm_release(cfm_work)
981 :
982 : ! get M^-1(k)*V^0.5(k)
983 : CALL parallel_gemm("N", "C", n_RI, n_RI, n_RI, z_one, cfm_M_inv_ikp, cfm_V_sqrt_ikp, &
984 262 : z_zero, cfm_M_inv_V_sqrt_ikp)
985 :
986 262 : CALL cp_cfm_release(cfm_M_inv_ikp)
987 :
988 262 : CALL timestop(handle)
989 :
990 524 : END SUBROUTINE compute_MinvVsqrt_Vsqrt
991 :
992 : ! **************************************************************************************************
993 : !> \brief ...
994 : !> \param bs_env ...
995 : !> \param mat_chi_Gamma_tau ...
996 : !> \param fm_W_MIC_time ...
997 : ! **************************************************************************************************
998 6 : SUBROUTINE read_W_MIC_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
999 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1000 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
1001 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1002 :
1003 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_W_MIC_time'
1004 :
1005 : INTEGER :: handle, i_t
1006 : REAL(KIND=dp) :: t1
1007 :
1008 6 : CALL timeset(routineN, handle)
1009 :
1010 6 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
1011 6 : CALL create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
1012 :
1013 106 : DO i_t = 1, bs_env%num_time_freq_points
1014 :
1015 100 : t1 = m_walltime()
1016 :
1017 100 : CALL fm_read(fm_W_MIC_time(i_t), bs_env, bs_env%W_time_name, i_t)
1018 :
1019 106 : IF (bs_env%unit_nr > 0) THEN
1020 : WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
1021 50 : 'Read W^MIC(iτ) from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1022 100 : ', Execution time', m_walltime() - t1, ' s'
1023 : END IF
1024 :
1025 : END DO
1026 :
1027 6 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1028 :
1029 : ! Marek : Reading of the W(w=0) potential for RTP
1030 : ! TODO : is the condition bs_env%all_W_exist sufficient for reading?
1031 : ! This block builds
1032 : ! bs_env%fm_W_MIC_freq_zero specifically for RT-BSE consumption (read by
1033 : ! rt_bse_linearized.F initialize_cohsex_selfenergy and by
1034 : ! rt_bse_ri_rs.F rt_bse_ri_rs_ensure_W0_grid). RT-BSE-specific compute
1035 : ! embedded in GW; left here because moving it would require keeping
1036 : ! fm_W_MIC_time alive past compute_W_MIC.
1037 6 : IF (bs_env%rtp_method == rtp_method_bse .OR. &
1038 : bs_env%rtp_method == rtp_method_bse_linearized) THEN
1039 4 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
1040 4 : t1 = m_walltime()
1041 4 : CALL fm_read(bs_env%fm_W_MIC_freq_zero, bs_env, "W_freq_rtp", 0)
1042 4 : IF (bs_env%unit_nr > 0) THEN
1043 : WRITE (bs_env%unit_nr, '(T2,A,I3,A,I3,A,F10.1,A)') &
1044 2 : 'Read W^MIC(f=0) from file for freq. point ', 1, ' /', 1, &
1045 4 : ', Execution time', m_walltime() - t1, ' s'
1046 : END IF
1047 : END IF
1048 :
1049 6 : CALL timestop(handle)
1050 :
1051 6 : END SUBROUTINE read_W_MIC_time
1052 :
1053 : ! **************************************************************************************************
1054 : !> \brief ...
1055 : !> \param bs_env ...
1056 : !> \param qs_env ...
1057 : !> \param mat_chi_Gamma_tau ...
1058 : !> \param fm_W_MIC_time ...
1059 : ! **************************************************************************************************
1060 62 : SUBROUTINE compute_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
1061 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1062 : TYPE(qs_environment_type), POINTER :: qs_env
1063 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
1064 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1065 :
1066 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_W_MIC'
1067 :
1068 : INTEGER :: handle, i_t, ikp, ikp_batch, &
1069 : ikp_in_batch, j_w
1070 : REAL(KIND=dp) :: t1
1071 : TYPE(cp_cfm_type) :: cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp
1072 62 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_V_kp
1073 :
1074 62 : CALL timeset(routineN, handle)
1075 :
1076 62 : CALL create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
1077 :
1078 172 : DO ikp_batch = 1, bs_env%num_chi_eps_W_batches
1079 :
1080 110 : t1 = m_walltime()
1081 :
1082 : ! Compute V_PQ(k) = sum_R e^(ikR) <phi_P, cell 0 | 1/r | phi_Q, cell R>
1083 110 : CALL compute_V_k_by_lattice_sum(bs_env, qs_env, fm_V_kp, ikp_batch)
1084 :
1085 550 : DO ikp_in_batch = 1, bs_env%nkp_chi_eps_W_batch
1086 :
1087 440 : ikp = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + ikp_in_batch
1088 :
1089 440 : IF (ikp > bs_env%nkp_chi_eps_W_orig_plus_extra) CYCLE
1090 :
1091 : CALL compute_MinvVsqrt_Vsqrt(bs_env, qs_env, fm_V_kp, &
1092 262 : cfm_V_sqrt_ikp, cfm_M_inv_V_sqrt_ikp, ikp)
1093 :
1094 262 : CALL bs_env%para_env%sync()
1095 262 : CALL cp_fm_release(fm_V_kp(ikp, 1))
1096 262 : CALL cp_fm_release(fm_V_kp(ikp, 2))
1097 :
1098 3116 : DO j_w = 1, bs_env%num_time_freq_points
1099 :
1100 : ! check if we need this (ikp, ω_j) combination for approximate k-point extrapolation
1101 2744 : IF (bs_env%approx_kp_extrapol .AND. j_w > 1 .AND. &
1102 : ikp > bs_env%nkp_chi_eps_W_orig) CYCLE
1103 :
1104 : CALL compute_fm_W_MIC_freq_j(bs_env, qs_env, bs_env%fm_W_MIC_freq, j_w, ikp, &
1105 : mat_chi_Gamma_tau, cfm_M_inv_V_sqrt_ikp, &
1106 2420 : cfm_V_sqrt_ikp)
1107 :
1108 : ! Fourier trafo from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
1109 3184 : CALL Fourier_transform_w_to_t(bs_env, fm_W_MIC_time, bs_env%fm_W_MIC_freq, j_w)
1110 :
1111 : END DO ! ω_j
1112 :
1113 : END DO ! ikp_in_batch
1114 :
1115 110 : DEALLOCATE (fm_V_kp)
1116 :
1117 172 : IF (bs_env%unit_nr > 0) THEN
1118 : WRITE (bs_env%unit_nr, '(T2,A,I12,A,I3,A,F10.1,A)') &
1119 55 : 'Computed W(iτ,k) for k-point batch', &
1120 55 : ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
1121 110 : ', Execution time', m_walltime() - t1, ' s'
1122 : END IF
1123 :
1124 : END DO ! ikp_batch
1125 :
1126 62 : IF (bs_env%approx_kp_extrapol) THEN
1127 2 : CALL apply_extrapol_factor(bs_env, fm_W_MIC_time)
1128 : END IF
1129 :
1130 : ! M^-1(k=0)*W^MIC(iτ)*M^-1(k=0)
1131 62 : CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_W_MIC_time)
1132 :
1133 1206 : DO i_t = 1, bs_env%num_time_freq_points
1134 1206 : CALL fm_write(fm_W_MIC_time(i_t), i_t, bs_env%W_time_name, qs_env)
1135 : END DO
1136 :
1137 62 : CALL cp_cfm_release(cfm_M_inv_V_sqrt_ikp)
1138 62 : CALL cp_cfm_release(cfm_V_sqrt_ikp)
1139 62 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
1140 :
1141 : ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0
1142 : ! Same RT-BSE coupling as read_W_MIC_time.
1143 62 : IF (bs_env%rtp_method == rtp_method_bse .OR. &
1144 : bs_env%rtp_method == rtp_method_bse_linearized) THEN
1145 54 : t1 = m_walltime()
1146 54 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
1147 : ! Set to zero
1148 54 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
1149 : ! Sum over all times
1150 1134 : DO i_t = 1, bs_env%num_time_freq_points
1151 : ! Add the relevant structure with correct weight
1152 : CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
1153 1134 : bs_env%imag_time_weights_freq_zero(i_t), fm_W_MIC_time(i_t))
1154 : END DO
1155 : ! Done, save to file
1156 54 : CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
1157 : ! Report calculation
1158 54 : IF (bs_env%unit_nr > 0) THEN
1159 : WRITE (bs_env%unit_nr, '(T2,A,I11,A,I3,A,F10.1,A)') &
1160 27 : 'Computed W(f=0,k) for k-point batch', &
1161 27 : 1, ' /', 1, &
1162 54 : ', Execution time', m_walltime() - t1, ' s'
1163 : END IF
1164 : END IF
1165 :
1166 62 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1167 :
1168 62 : CALL timestop(handle)
1169 :
1170 124 : END SUBROUTINE compute_W_MIC
1171 :
1172 : ! **************************************************************************************************
1173 : !> \brief ...
1174 : !> \param bs_env ...
1175 : !> \param qs_env ...
1176 : !> \param fm_W_MIC_freq_j ...
1177 : !> \param j_w ...
1178 : !> \param ikp ...
1179 : !> \param mat_chi_Gamma_tau ...
1180 : !> \param cfm_M_inv_V_sqrt_ikp ...
1181 : !> \param cfm_V_sqrt_ikp ...
1182 : ! **************************************************************************************************
1183 2420 : SUBROUTINE compute_fm_W_MIC_freq_j(bs_env, qs_env, fm_W_MIC_freq_j, j_w, ikp, mat_chi_Gamma_tau, &
1184 : cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
1185 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1186 : TYPE(qs_environment_type), POINTER :: qs_env
1187 : TYPE(cp_fm_type) :: fm_W_MIC_freq_j
1188 : INTEGER :: j_w, ikp
1189 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
1190 : TYPE(cp_cfm_type) :: cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp
1191 :
1192 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_MIC_freq_j'
1193 :
1194 : INTEGER :: handle
1195 : TYPE(cp_cfm_type) :: cfm_chi_ikp_freq_j, cfm_W_ikp_freq_j
1196 :
1197 2420 : CALL timeset(routineN, handle)
1198 :
1199 : ! 1. Fourier transformation of χ_PQ(iτ,k=0) to χ_PQ(iω_j,k=0)
1200 2420 : CALL compute_fm_chi_Gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
1201 :
1202 2420 : CALL cp_fm_set_all(fm_W_MIC_freq_j, 0.0_dp)
1203 :
1204 : ! 2. Get χ_PQ(iω_j,k_i) from χ_PQ(iω_j,k=0) using the minimum image convention
1205 : CALL cfm_ikp_from_fm_Gamma(cfm_chi_ikp_freq_j, bs_env%fm_chi_Gamma_freq, &
1206 2420 : ikp, qs_env, bs_env%kpoints_chi_eps_W, "RI_AUX")
1207 :
1208 : ! 3. Remove all negative eigenvalues from χ_PQ(iω_j,k_i)
1209 2420 : CALL cp_cfm_power(cfm_chi_ikp_freq_j, threshold=0.0_dp, exponent=1.0_dp)
1210 :
1211 : ! 4. ε(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)
1212 : ! W(iω_j,k_i) = V^0.5(k_i)*(ε^-1(iω_j,k_i)-Id)*V^0.5(k_i)
1213 : CALL compute_cfm_W_ikp_freq_j(bs_env, cfm_chi_ikp_freq_j, cfm_V_sqrt_ikp, &
1214 2420 : cfm_M_inv_V_sqrt_ikp, cfm_W_ikp_freq_j)
1215 :
1216 : ! 5. k-point integration W_PQ(iω_j, k_i) to W_PQ^MIC(iω_j)
1217 2420 : SELECT CASE (bs_env%approx_kp_extrapol)
1218 : CASE (.FALSE.)
1219 : ! default: standard k-point extrapolation
1220 : CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, cfm_W_ikp_freq_j, ikp, &
1221 2420 : bs_env%kpoints_chi_eps_W, "RI_AUX")
1222 : CASE (.TRUE.)
1223 : ! for approximate kpoint extrapolation: get W_PQ^MIC(iω_1) with and without k-point
1224 : ! extrapolation to compute the extrapolation factor f_PQ for every PQ-matrix element,
1225 : ! f_PQ = (W_PQ^MIC(iω_1) with extrapolation) / (W_PQ^MIC(iω_1) without extrapolation)
1226 :
1227 : ! for ω_1, we compute the k-point extrapolated result using all k-points
1228 196 : IF (j_w == 1) THEN
1229 :
1230 : ! k-point extrapolated
1231 : CALL MIC_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_extra, &
1232 : cfm_W_ikp_freq_j, ikp, bs_env%kpoints_chi_eps_W, &
1233 52 : "RI_AUX")
1234 : ! non-kpoint extrapolated
1235 52 : IF (ikp <= bs_env%nkp_chi_eps_W_orig) THEN
1236 : CALL MIC_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_no_extra, &
1237 : cfm_W_ikp_freq_j, ikp, bs_env%kpoints_chi_eps_W, &
1238 16 : "RI_AUX", wkp_ext=bs_env%wkp_orig)
1239 : END IF
1240 :
1241 : END IF
1242 :
1243 : ! for all ω_j, we need to compute W^MIC without k-point extrpolation
1244 196 : IF (ikp <= bs_env%nkp_chi_eps_W_orig) THEN
1245 : CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, cfm_W_ikp_freq_j, &
1246 : ikp, bs_env%kpoints_chi_eps_W, "RI_AUX", &
1247 160 : wkp_ext=bs_env%wkp_orig)
1248 : END IF
1249 : END SELECT
1250 :
1251 2420 : CALL cp_cfm_release(cfm_W_ikp_freq_j)
1252 :
1253 2420 : CALL timestop(handle)
1254 :
1255 2420 : END SUBROUTINE compute_fm_W_MIC_freq_j
1256 :
1257 : ! **************************************************************************************************
1258 : !> \brief ...
1259 : !> \param cfm_mat ...
1260 : ! **************************************************************************************************
1261 524 : SUBROUTINE clean_lower_part(cfm_mat)
1262 : TYPE(cp_cfm_type) :: cfm_mat
1263 :
1264 : CHARACTER(LEN=*), PARAMETER :: routineN = 'clean_lower_part'
1265 :
1266 : INTEGER :: handle, i_row, j_col, j_global, &
1267 : ncol_local, nrow_local
1268 262 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1269 :
1270 262 : CALL timeset(routineN, handle)
1271 :
1272 : CALL cp_cfm_get_info(matrix=cfm_mat, &
1273 : nrow_local=nrow_local, ncol_local=ncol_local, &
1274 262 : row_indices=row_indices, col_indices=col_indices)
1275 :
1276 3858 : DO j_col = 1, ncol_local
1277 3596 : j_global = col_indices(j_col)
1278 57458 : DO i_row = 1, nrow_local
1279 57196 : IF (j_global < row_indices(i_row)) cfm_mat%local_data(i_row, j_col) = z_zero
1280 : END DO
1281 : END DO
1282 :
1283 262 : CALL timestop(handle)
1284 :
1285 262 : END SUBROUTINE clean_lower_part
1286 :
1287 : ! **************************************************************************************************
1288 : !> \brief ...
1289 : !> \param bs_env ...
1290 : !> \param fm_W_MIC_time ...
1291 : ! **************************************************************************************************
1292 4 : SUBROUTINE apply_extrapol_factor(bs_env, fm_W_MIC_time)
1293 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1294 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1295 :
1296 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_extrapol_factor'
1297 :
1298 : INTEGER :: handle, i, i_t, j, ncol_local, nrow_local
1299 : REAL(KIND=dp) :: extrapol_factor, W_extra_1, W_no_extra_1
1300 :
1301 2 : CALL timeset(routineN, handle)
1302 :
1303 2 : CALL cp_fm_get_info(matrix=fm_W_MIC_time(1), nrow_local=nrow_local, ncol_local=ncol_local)
1304 :
1305 22 : DO i_t = 1, bs_env%num_time_freq_points
1306 122 : DO j = 1, ncol_local
1307 370 : DO i = 1, nrow_local
1308 :
1309 250 : W_extra_1 = bs_env%fm_W_MIC_freq_1_extra%local_data(i, j)
1310 250 : W_no_extra_1 = bs_env%fm_W_MIC_freq_1_no_extra%local_data(i, j)
1311 :
1312 250 : IF (ABS(W_no_extra_1) > 1.0E-13) THEN
1313 190 : extrapol_factor = ABS(W_extra_1/W_no_extra_1)
1314 : ELSE
1315 : extrapol_factor = 1.0_dp
1316 : END IF
1317 :
1318 : ! reset extrapolation factor if it is very large
1319 190 : IF (extrapol_factor > 10.0_dp) extrapol_factor = 1.0_dp
1320 :
1321 : fm_W_MIC_time(i_t)%local_data(i, j) = fm_W_MIC_time(i_t)%local_data(i, j) &
1322 350 : *extrapol_factor
1323 : END DO
1324 : END DO
1325 : END DO
1326 :
1327 2 : CALL timestop(handle)
1328 :
1329 2 : END SUBROUTINE apply_extrapol_factor
1330 :
1331 : ! **************************************************************************************************
1332 : !> \brief ...
1333 : !> \param bs_env ...
1334 : !> \param fm_chi_Gamma_freq ...
1335 : !> \param j_w ...
1336 : !> \param mat_chi_Gamma_tau ...
1337 : ! **************************************************************************************************
1338 2930 : SUBROUTINE compute_fm_chi_Gamma_freq(bs_env, fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
1339 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1340 : TYPE(cp_fm_type) :: fm_chi_Gamma_freq
1341 : INTEGER :: j_w
1342 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
1343 :
1344 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_chi_Gamma_freq'
1345 :
1346 : INTEGER :: handle, i_t
1347 : REAL(KIND=dp) :: freq_j, time_i, weight_ij
1348 :
1349 2930 : CALL timeset(routineN, handle)
1350 :
1351 2930 : CALL dbcsr_set(bs_env%mat_RI_RI%matrix, 0.0_dp)
1352 :
1353 2930 : freq_j = bs_env%imag_freq_points(j_w)
1354 :
1355 43484 : DO i_t = 1, bs_env%num_time_freq_points
1356 :
1357 40554 : time_i = bs_env%imag_time_points(i_t)
1358 40554 : weight_ij = bs_env%weights_cos_t_to_w(j_w, i_t)
1359 :
1360 : ! actual Fourier transform
1361 : CALL dbcsr_add(bs_env%mat_RI_RI%matrix, mat_chi_Gamma_tau(i_t)%matrix, &
1362 43484 : 1.0_dp, COS(time_i*freq_j)*weight_ij)
1363 :
1364 : END DO
1365 :
1366 2930 : CALL copy_dbcsr_to_fm(bs_env%mat_RI_RI%matrix, fm_chi_Gamma_freq)
1367 :
1368 2930 : CALL timestop(handle)
1369 :
1370 2930 : END SUBROUTINE compute_fm_chi_Gamma_freq
1371 :
1372 : ! **************************************************************************************************
1373 : !> \brief ...
1374 : !> \param mat_ikp_re ...
1375 : !> \param mat_ikp_im ...
1376 : !> \param mat_Gamma ...
1377 : !> \param kpoints ...
1378 : !> \param ikp ...
1379 : !> \param qs_env ...
1380 : ! **************************************************************************************************
1381 0 : SUBROUTINE mat_ikp_from_mat_Gamma(mat_ikp_re, mat_ikp_im, mat_Gamma, kpoints, ikp, qs_env)
1382 : TYPE(dbcsr_type) :: mat_ikp_re, mat_ikp_im, mat_Gamma
1383 : TYPE(kpoint_type), POINTER :: kpoints
1384 : INTEGER :: ikp
1385 : TYPE(qs_environment_type), POINTER :: qs_env
1386 :
1387 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mat_ikp_from_mat_Gamma'
1388 :
1389 : INTEGER :: col, handle, i_cell, j_cell, num_cells, &
1390 : row
1391 0 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1392 : LOGICAL :: f, i_cell_is_the_minimum_image_cell
1393 : REAL(KIND=dp) :: abs_rab_cell_i, abs_rab_cell_j, arg
1394 : REAL(KIND=dp), DIMENSION(3) :: cell_vector, cell_vector_j, rab_cell_i, &
1395 : rab_cell_j
1396 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1397 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_im, block_re, data_block
1398 : TYPE(cell_type), POINTER :: cell
1399 : TYPE(dbcsr_iterator_type) :: iter
1400 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1401 :
1402 0 : CALL timeset(routineN, handle)
1403 :
1404 : ! get the same blocks in mat_ikp_re and mat_ikp_im as in mat_Gamma
1405 0 : CALL dbcsr_copy(mat_ikp_re, mat_Gamma)
1406 0 : CALL dbcsr_copy(mat_ikp_im, mat_Gamma)
1407 0 : CALL dbcsr_set(mat_ikp_re, 0.0_dp)
1408 0 : CALL dbcsr_set(mat_ikp_im, 0.0_dp)
1409 :
1410 0 : NULLIFY (cell, particle_set)
1411 0 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1412 0 : CALL get_cell(cell=cell, h=hmat)
1413 :
1414 0 : index_to_cell => kpoints%index_to_cell
1415 :
1416 0 : num_cells = SIZE(index_to_cell, 2)
1417 :
1418 0 : DO i_cell = 1, num_cells
1419 :
1420 0 : CALL dbcsr_iterator_start(iter, mat_Gamma)
1421 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1422 0 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
1423 :
1424 0 : cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, i_cell), dp))
1425 :
1426 : rab_cell_i(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
1427 0 : (pbc(particle_set(col)%r(1:3), cell) + cell_vector(1:3))
1428 0 : abs_rab_cell_i = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
1429 :
1430 : ! minimum image convention
1431 0 : i_cell_is_the_minimum_image_cell = .TRUE.
1432 0 : DO j_cell = 1, num_cells
1433 0 : cell_vector_j(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, j_cell), dp))
1434 : rab_cell_j(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
1435 0 : (pbc(particle_set(col)%r(1:3), cell) + cell_vector_j(1:3))
1436 0 : abs_rab_cell_j = SQRT(rab_cell_j(1)**2 + rab_cell_j(2)**2 + rab_cell_j(3)**2)
1437 :
1438 0 : IF (abs_rab_cell_i > abs_rab_cell_j + 1.0E-6_dp) THEN
1439 0 : i_cell_is_the_minimum_image_cell = .FALSE.
1440 : END IF
1441 : END DO
1442 :
1443 0 : IF (i_cell_is_the_minimum_image_cell) THEN
1444 0 : NULLIFY (block_re, block_im)
1445 0 : CALL dbcsr_get_block_p(matrix=mat_ikp_re, row=row, col=col, block=block_re, found=f)
1446 0 : CALL dbcsr_get_block_p(matrix=mat_ikp_im, row=row, col=col, block=block_im, found=f)
1447 0 : CPASSERT(ALL(ABS(block_re) < 1.0E-10_dp))
1448 0 : CPASSERT(ALL(ABS(block_im) < 1.0E-10_dp))
1449 :
1450 : arg = REAL(index_to_cell(1, i_cell), dp)*kpoints%xkp(1, ikp) + &
1451 : REAL(index_to_cell(2, i_cell), dp)*kpoints%xkp(2, ikp) + &
1452 0 : REAL(index_to_cell(3, i_cell), dp)*kpoints%xkp(3, ikp)
1453 :
1454 0 : block_re(:, :) = COS(twopi*arg)*data_block(:, :)
1455 0 : block_im(:, :) = SIN(twopi*arg)*data_block(:, :)
1456 : END IF
1457 :
1458 : END DO
1459 0 : CALL dbcsr_iterator_stop(iter)
1460 :
1461 : END DO
1462 :
1463 0 : CALL timestop(handle)
1464 :
1465 0 : END SUBROUTINE mat_ikp_from_mat_Gamma
1466 :
1467 : ! **************************************************************************************************
1468 : !> \brief ...
1469 : !> \param bs_env ...
1470 : !> \param cfm_chi_ikp_freq_j ...
1471 : !> \param cfm_V_sqrt_ikp ...
1472 : !> \param cfm_M_inv_V_sqrt_ikp ...
1473 : !> \param cfm_W_ikp_freq_j ...
1474 : ! **************************************************************************************************
1475 12100 : SUBROUTINE compute_cfm_W_ikp_freq_j(bs_env, cfm_chi_ikp_freq_j, cfm_V_sqrt_ikp, &
1476 : cfm_M_inv_V_sqrt_ikp, cfm_W_ikp_freq_j)
1477 :
1478 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1479 : TYPE(cp_cfm_type) :: cfm_chi_ikp_freq_j, cfm_V_sqrt_ikp, &
1480 : cfm_M_inv_V_sqrt_ikp, cfm_W_ikp_freq_j
1481 :
1482 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cfm_W_ikp_freq_j'
1483 :
1484 : INTEGER :: handle, info, n_RI
1485 : TYPE(cp_cfm_type) :: cfm_eps_ikp_freq_j, cfm_work
1486 :
1487 2420 : CALL timeset(routineN, handle)
1488 :
1489 2420 : CALL cp_cfm_create(cfm_work, cfm_chi_ikp_freq_j%matrix_struct)
1490 2420 : n_RI = bs_env%n_RI
1491 :
1492 : ! 1. ε(iω_j,k) = Id - V^0.5(k)*M^-1(k)*χ(iω_j,k)*M^-1(k)*V^0.5(k)
1493 :
1494 : ! 1. a) work = χ(iω_j,k)*M^-1(k)*V^0.5(k)
1495 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, z_one, &
1496 2420 : cfm_chi_ikp_freq_j, cfm_M_inv_V_sqrt_ikp, z_zero, cfm_work)
1497 2420 : CALL cp_cfm_release(cfm_chi_ikp_freq_j)
1498 :
1499 : ! 1. b) eps_work = V^0.5(k)*M^-1(k)*work
1500 2420 : CALL cp_cfm_create(cfm_eps_ikp_freq_j, cfm_work%matrix_struct)
1501 : CALL parallel_gemm('C', 'N', n_RI, n_RI, n_RI, z_one, &
1502 2420 : cfm_M_inv_V_sqrt_ikp, cfm_work, z_zero, cfm_eps_ikp_freq_j)
1503 :
1504 : ! 1. c) ε(iω_j,k) = eps_work - Id
1505 2420 : CALL cfm_add_on_diag(cfm_eps_ikp_freq_j, z_one)
1506 :
1507 : ! 2. W(iω_j,k) = V^0.5(k)*(ε^-1(iω_j,k)-Id)*V^0.5(k)
1508 :
1509 : ! 2. a) Cholesky decomposition of ε(iω_j,k) as preparation for inversion
1510 2420 : CALL cp_cfm_cholesky_decompose(matrix=cfm_eps_ikp_freq_j, n=n_RI, info_out=info)
1511 2420 : CPASSERT(info == 0)
1512 :
1513 : ! 2. b) Inversion of ε(iω_j,k) using its Cholesky decomposition
1514 2420 : CALL cp_cfm_cholesky_invert(cfm_eps_ikp_freq_j)
1515 2420 : CALL cp_cfm_uplo_to_full(cfm_eps_ikp_freq_j)
1516 :
1517 : ! 2. c) ε^-1(iω_j,k)-Id
1518 2420 : CALL cfm_add_on_diag(cfm_eps_ikp_freq_j, -z_one)
1519 :
1520 : ! 2. d) work = (ε^-1(iω_j,k)-Id)*V^0.5(k)
1521 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, z_one, cfm_eps_ikp_freq_j, cfm_V_sqrt_ikp, &
1522 2420 : z_zero, cfm_work)
1523 :
1524 : ! 2. e) W(iw,k) = V^0.5(k)*work
1525 2420 : CALL cp_cfm_create(cfm_W_ikp_freq_j, cfm_work%matrix_struct)
1526 : CALL parallel_gemm('C', 'N', n_RI, n_RI, n_RI, z_one, cfm_V_sqrt_ikp, cfm_work, &
1527 2420 : z_zero, cfm_W_ikp_freq_j)
1528 :
1529 2420 : CALL cp_cfm_release(cfm_work)
1530 2420 : CALL cp_cfm_release(cfm_eps_ikp_freq_j)
1531 :
1532 2420 : CALL timestop(handle)
1533 :
1534 2420 : END SUBROUTINE compute_cfm_W_ikp_freq_j
1535 :
1536 : ! **************************************************************************************************
1537 : !> \brief ...
1538 : !> \param cfm ...
1539 : !> \param alpha ...
1540 : ! **************************************************************************************************
1541 9680 : SUBROUTINE cfm_add_on_diag(cfm, alpha)
1542 :
1543 : TYPE(cp_cfm_type) :: cfm
1544 : COMPLEX(KIND=dp) :: alpha
1545 :
1546 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cfm_add_on_diag'
1547 :
1548 : INTEGER :: handle, i_row, j_col, j_global, &
1549 : ncol_local, nrow_local
1550 4840 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1551 :
1552 4840 : CALL timeset(routineN, handle)
1553 :
1554 : CALL cp_cfm_get_info(matrix=cfm, &
1555 : nrow_local=nrow_local, &
1556 : ncol_local=ncol_local, &
1557 : row_indices=row_indices, &
1558 4840 : col_indices=col_indices)
1559 :
1560 : ! add 1 on the diagonal
1561 111744 : DO j_col = 1, ncol_local
1562 106904 : j_global = col_indices(j_col)
1563 2128460 : DO i_row = 1, nrow_local
1564 2123620 : IF (j_global == row_indices(i_row)) THEN
1565 53452 : cfm%local_data(i_row, j_col) = cfm%local_data(i_row, j_col) + alpha
1566 : END IF
1567 : END DO
1568 : END DO
1569 :
1570 4840 : CALL timestop(handle)
1571 :
1572 4840 : END SUBROUTINE cfm_add_on_diag
1573 :
1574 : ! **************************************************************************************************
1575 : !> \brief ...
1576 : !> \param bs_env ...
1577 : !> \param fm_W_MIC_time ...
1578 : ! **************************************************************************************************
1579 104 : SUBROUTINE create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
1580 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1581 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1582 :
1583 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_fm_W_MIC_time'
1584 :
1585 : INTEGER :: handle, i_t
1586 :
1587 104 : CALL timeset(routineN, handle)
1588 :
1589 2066 : ALLOCATE (fm_W_MIC_time(bs_env%num_time_freq_points))
1590 1858 : DO i_t = 1, bs_env%num_time_freq_points
1591 1858 : CALL cp_fm_create(fm_W_MIC_time(i_t), bs_env%fm_RI_RI%matrix_struct, set_zero=.TRUE.)
1592 : END DO
1593 :
1594 104 : CALL timestop(handle)
1595 :
1596 104 : END SUBROUTINE create_fm_W_MIC_time
1597 :
1598 : ! **************************************************************************************************
1599 : !> \brief ...
1600 : !> \param bs_env ...
1601 : !> \param fm_W_MIC_time ...
1602 : !> \param fm_W_MIC_freq_j ...
1603 : !> \param j_w ...
1604 : ! **************************************************************************************************
1605 2930 : SUBROUTINE Fourier_transform_w_to_t(bs_env, fm_W_MIC_time, fm_W_MIC_freq_j, j_w)
1606 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1607 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1608 : TYPE(cp_fm_type) :: fm_W_MIC_freq_j
1609 : INTEGER :: j_w
1610 :
1611 : CHARACTER(LEN=*), PARAMETER :: routineN = 'Fourier_transform_w_to_t'
1612 :
1613 : INTEGER :: handle, i_t
1614 : REAL(KIND=dp) :: freq_j, time_i, weight_ij
1615 :
1616 2930 : CALL timeset(routineN, handle)
1617 :
1618 2930 : freq_j = bs_env%imag_freq_points(j_w)
1619 :
1620 43484 : DO i_t = 1, bs_env%num_time_freq_points
1621 :
1622 40554 : time_i = bs_env%imag_time_points(i_t)
1623 40554 : weight_ij = bs_env%weights_cos_w_to_t(i_t, j_w)
1624 :
1625 : ! actual Fourier transform
1626 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_W_MIC_time(i_t), &
1627 43484 : beta=weight_ij*COS(time_i*freq_j), matrix_b=fm_W_MIC_freq_j)
1628 :
1629 : END DO
1630 :
1631 2930 : CALL timestop(handle)
1632 :
1633 2930 : END SUBROUTINE Fourier_transform_w_to_t
1634 :
1635 : ! **************************************************************************************************
1636 : !> \brief ...
1637 : !> \param bs_env ...
1638 : !> \param qs_env ...
1639 : !> \param fm_W_MIC_time ...
1640 : ! **************************************************************************************************
1641 240 : SUBROUTINE multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_W_MIC_time)
1642 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1643 : TYPE(qs_environment_type), POINTER :: qs_env
1644 : TYPE(cp_fm_type), DIMENSION(:) :: fm_W_MIC_time
1645 :
1646 : CHARACTER(LEN=*), PARAMETER :: routineN = 'multiply_fm_W_MIC_time_with_Minv_Gamma'
1647 :
1648 : INTEGER :: handle, i_t, n_RI, ndep
1649 : TYPE(cp_fm_type) :: fm_work
1650 240 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Minv_Gamma
1651 :
1652 240 : CALL timeset(routineN, handle)
1653 :
1654 240 : n_RI = bs_env%n_RI
1655 :
1656 240 : CALL cp_fm_create(fm_work, fm_W_MIC_time(1)%matrix_struct)
1657 :
1658 : ! compute Gamma-only RI-metric matrix M(k=0); no regularization
1659 : CALL RI_2c_integral_mat(qs_env, fm_Minv_Gamma, fm_W_MIC_time(1), n_RI, &
1660 240 : bs_env%ri_metric, do_kpoints=.FALSE.)
1661 :
1662 240 : CALL cp_fm_power(fm_Minv_Gamma(1, 1), fm_work, -1.0_dp, 0.0_dp, ndep)
1663 :
1664 : ! M^-1(k=0)*W^MIC(iτ)*M^-1(k=0)
1665 2036 : DO i_t = 1, SIZE(fm_W_MIC_time)
1666 :
1667 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, fm_Minv_Gamma(1, 1), &
1668 1796 : fm_W_MIC_time(i_t), 0.0_dp, fm_work)
1669 :
1670 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, fm_work, &
1671 2036 : fm_Minv_Gamma(1, 1), 0.0_dp, fm_W_MIC_time(i_t))
1672 :
1673 : END DO
1674 :
1675 240 : CALL cp_fm_release(fm_work)
1676 240 : CALL cp_fm_release(fm_Minv_Gamma)
1677 :
1678 240 : CALL timestop(handle)
1679 :
1680 480 : END SUBROUTINE multiply_fm_W_MIC_time_with_Minv_Gamma
1681 :
1682 : ! **************************************************************************************************
1683 : !> \brief ...
1684 : !> \param bs_env ...
1685 : !> \param qs_env ...
1686 : !> \param fm_Sigma_x_Gamma ...
1687 : ! **************************************************************************************************
1688 68 : SUBROUTINE get_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
1689 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1690 : TYPE(qs_environment_type), POINTER :: qs_env
1691 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
1692 :
1693 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_Sigma_x'
1694 :
1695 : INTEGER :: handle, ispin
1696 :
1697 68 : CALL timeset(routineN, handle)
1698 :
1699 282 : ALLOCATE (fm_Sigma_x_Gamma(bs_env%n_spin))
1700 146 : DO ispin = 1, bs_env%n_spin
1701 146 : CALL cp_fm_create(fm_Sigma_x_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1702 : END DO
1703 :
1704 68 : IF (bs_env%Sigma_x_exists) THEN
1705 14 : DO ispin = 1, bs_env%n_spin
1706 14 : CALL fm_read(fm_Sigma_x_Gamma(ispin), bs_env, bs_env%Sigma_x_name, ispin)
1707 : END DO
1708 : ELSE
1709 62 : CALL compute_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
1710 : END IF
1711 :
1712 68 : CALL timestop(handle)
1713 :
1714 68 : END SUBROUTINE get_Sigma_x
1715 :
1716 : ! **************************************************************************************************
1717 : !> \brief ...
1718 : !> \param bs_env ...
1719 : !> \param qs_env ...
1720 : !> \param fm_Sigma_x_Gamma ...
1721 : ! **************************************************************************************************
1722 62 : SUBROUTINE compute_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
1723 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1724 : TYPE(qs_environment_type), POINTER :: qs_env
1725 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
1726 :
1727 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_x'
1728 :
1729 : INTEGER :: handle, i_intval_idx, ispin, j_intval_idx
1730 : INTEGER, DIMENSION(2) :: i_atoms, j_atoms
1731 : REAL(KIND=dp) :: t1
1732 62 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Vtr_Gamma
1733 : TYPE(dbcsr_type) :: mat_Sigma_x_Gamma
1734 2046 : TYPE(dbt_type) :: t_2c_D, t_2c_Sigma_x, t_2c_V, t_3c_x_V
1735 :
1736 62 : CALL timeset(routineN, handle)
1737 :
1738 62 : t1 = m_walltime()
1739 :
1740 62 : CALL dbt_create(bs_env%t_G, t_2c_D)
1741 62 : CALL dbt_create(bs_env%t_W, t_2c_V)
1742 62 : CALL dbt_create(bs_env%t_G, t_2c_Sigma_x)
1743 62 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_V)
1744 62 : CALL dbcsr_create(mat_Sigma_x_Gamma, template=bs_env%mat_ao_ao%matrix)
1745 :
1746 : ! 1. Compute truncated Coulomb operator matrix V^tr(k=0) (cutoff rad: cellsize/2)
1747 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI, bs_env%n_RI, &
1748 62 : bs_env%trunc_coulomb, do_kpoints=.FALSE.)
1749 :
1750 : ! 2. Compute M^-1(k=0) and get M^-1(k=0)*V^tr(k=0)*M^-1(k=0)
1751 62 : CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_Vtr_Gamma(:, 1))
1752 :
1753 132 : DO ispin = 1, bs_env%n_spin
1754 :
1755 : ! 3. Compute density matrix D_µν
1756 70 : CALL G_occ_vir(bs_env, 0.0_dp, bs_env%fm_work_mo(2), ispin, occ=.TRUE., vir=.FALSE.)
1757 :
1758 : CALL fm_to_local_tensor(bs_env%fm_work_mo(2), bs_env%mat_ao_ao%matrix, &
1759 : bs_env%mat_ao_ao_tensor%matrix, t_2c_D, bs_env, &
1760 70 : bs_env%atoms_i_t_group)
1761 :
1762 : CALL fm_to_local_tensor(fm_Vtr_Gamma(1, 1), bs_env%mat_RI_RI%matrix, &
1763 : bs_env%mat_RI_RI_tensor%matrix, t_2c_V, bs_env, &
1764 70 : bs_env%atoms_j_t_group)
1765 :
1766 : ! every group has its own range of i_atoms and j_atoms; only deal with a
1767 : ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
1768 140 : DO i_intval_idx = 1, bs_env%n_intervals_i
1769 210 : DO j_intval_idx = 1, bs_env%n_intervals_j
1770 210 : i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
1771 210 : j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
1772 :
1773 : ! 4. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
1774 : ! 5. M_Qνσ(iτ) = sum_P (νσ|P) (M^-1(k=0)*V^tr(k=0)*M^-1(k=0))_QP(iτ)
1775 70 : CALL compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_V, t_2c_V)
1776 :
1777 : ! 6. tensor operations with D and computation of Σ^x
1778 : ! Σ^x_λσ(k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) D_νµ
1779 : CALL contract_to_Sigma(t_2c_D, t_3c_x_V, t_2c_Sigma_x, i_atoms, j_atoms, &
1780 140 : qs_env, bs_env, occ=.TRUE., vir=.FALSE.)
1781 :
1782 : END DO ! j_atoms
1783 : END DO ! i_atoms
1784 :
1785 : CALL local_dbt_to_global_mat(t_2c_Sigma_x, bs_env%mat_ao_ao_tensor%matrix, &
1786 70 : mat_Sigma_x_Gamma, bs_env%para_env)
1787 :
1788 : CALL write_matrix(mat_Sigma_x_Gamma, ispin, bs_env%Sigma_x_name, &
1789 70 : bs_env%fm_work_mo(1), qs_env)
1790 :
1791 132 : CALL copy_dbcsr_to_fm(mat_Sigma_x_Gamma, fm_Sigma_x_Gamma(ispin))
1792 :
1793 : END DO ! ispin
1794 :
1795 62 : IF (bs_env%unit_nr > 0) THEN
1796 : WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
1797 31 : 'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
1798 31 : WRITE (bs_env%unit_nr, '(A)') ' '
1799 : END IF
1800 :
1801 62 : CALL dbcsr_release(mat_Sigma_x_Gamma)
1802 62 : CALL dbt_destroy(t_2c_D)
1803 62 : CALL dbt_destroy(t_2c_V)
1804 62 : CALL dbt_destroy(t_2c_Sigma_x)
1805 62 : CALL dbt_destroy(t_3c_x_V)
1806 62 : CALL cp_fm_release(fm_Vtr_Gamma)
1807 :
1808 62 : CALL timestop(handle)
1809 :
1810 124 : END SUBROUTINE compute_Sigma_x
1811 :
1812 : ! **************************************************************************************************
1813 : !> \brief ...
1814 : !> \param bs_env ...
1815 : !> \param qs_env ...
1816 : !> \param fm_W_MIC_time ...
1817 : !> \param fm_Sigma_c_Gamma_time ...
1818 : ! **************************************************************************************************
1819 68 : SUBROUTINE get_Sigma_c(bs_env, qs_env, fm_W_MIC_time, fm_Sigma_c_Gamma_time)
1820 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1821 : TYPE(qs_environment_type), POINTER :: qs_env
1822 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1823 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
1824 :
1825 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_Sigma_c'
1826 :
1827 : INTEGER :: handle, i_intval_idx, i_t, ispin, &
1828 : j_intval_idx, read_write_index
1829 : INTEGER, DIMENSION(2) :: i_atoms, j_atoms
1830 : REAL(KIND=dp) :: t1, tau
1831 68 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
1832 1156 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, &
1833 612 : t_2c_Sigma_neg_tau, &
1834 1700 : t_2c_Sigma_pos_tau, t_2c_W, t_3c_x_W
1835 :
1836 68 : CALL timeset(routineN, handle)
1837 :
1838 : CALL create_mat_for_Sigma_c(bs_env, t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
1839 : t_2c_Sigma_pos_tau, t_3c_x_W, &
1840 68 : mat_Sigma_neg_tau, mat_Sigma_pos_tau)
1841 :
1842 1312 : DO i_t = 1, bs_env%num_time_freq_points
1843 :
1844 2716 : DO ispin = 1, bs_env%n_spin
1845 :
1846 1404 : t1 = m_walltime()
1847 :
1848 1404 : read_write_index = i_t + (ispin - 1)*bs_env%num_time_freq_points
1849 :
1850 : ! read self-energy from restart
1851 1404 : IF (bs_env%Sigma_c_exists(i_t, ispin)) THEN
1852 120 : CALL fm_read(bs_env%fm_work_mo(1), bs_env, bs_env%Sigma_p_name, read_write_index)
1853 : CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_Sigma_pos_tau(i_t, ispin)%matrix, &
1854 120 : keep_sparsity=.FALSE.)
1855 120 : CALL fm_read(bs_env%fm_work_mo(1), bs_env, bs_env%Sigma_n_name, read_write_index)
1856 : CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_Sigma_neg_tau(i_t, ispin)%matrix, &
1857 120 : keep_sparsity=.FALSE.)
1858 120 : IF (bs_env%unit_nr > 0) THEN
1859 60 : WRITE (bs_env%unit_nr, '(T2,2A,I3,A,I3,A,F10.1,A)') 'Read Σ^c(iτ,k=0) ', &
1860 60 : 'from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1861 120 : ', Execution time', m_walltime() - t1, ' s'
1862 : END IF
1863 :
1864 : CYCLE
1865 :
1866 : END IF
1867 :
1868 1284 : tau = bs_env%imag_time_points(i_t)
1869 :
1870 1284 : CALL G_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.TRUE., vir=.FALSE.)
1871 1284 : CALL G_occ_vir(bs_env, tau, bs_env%fm_Gvir, ispin, occ=.FALSE., vir=.TRUE.)
1872 :
1873 : ! fm G^occ, G^vir and W to local tensor
1874 : CALL fm_to_local_tensor(bs_env%fm_Gocc, bs_env%mat_ao_ao%matrix, &
1875 : bs_env%mat_ao_ao_tensor%matrix, t_2c_Gocc, bs_env, &
1876 1284 : bs_env%atoms_i_t_group)
1877 : CALL fm_to_local_tensor(bs_env%fm_Gvir, bs_env%mat_ao_ao%matrix, &
1878 : bs_env%mat_ao_ao_tensor%matrix, t_2c_Gvir, bs_env, &
1879 1284 : bs_env%atoms_i_t_group)
1880 : CALL fm_to_local_tensor(fm_W_MIC_time(i_t), bs_env%mat_RI_RI%matrix, &
1881 : bs_env%mat_RI_RI_tensor%matrix, t_2c_W, bs_env, &
1882 1284 : bs_env%atoms_j_t_group)
1883 :
1884 : ! every group has its own range of i_atoms and j_atoms; only deal with a
1885 : ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
1886 2568 : DO i_intval_idx = 1, bs_env%n_intervals_i
1887 3852 : DO j_intval_idx = 1, bs_env%n_intervals_j
1888 3852 : i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
1889 3852 : j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
1890 :
1891 1284 : IF (bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx) .AND. &
1892 : bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) THEN
1893 : ! Do that only after first timestep to avoid skips due to vanishing G
1894 : ! caused by gaps
1895 18 : IF (i_t == 2) THEN
1896 0 : bs_env%n_skip_sigma = bs_env%n_skip_sigma + 1
1897 : END IF
1898 : CYCLE
1899 : END IF
1900 :
1901 : ! 1. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
1902 : ! 2. tensor operation M_Qνσ(iτ) = sum_P (νσ|P) W^MIC_QP(iτ)
1903 1266 : CALL compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_W, t_2c_W)
1904 :
1905 : ! 3. Σ_λσ(iτ,k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) G^occ_νµ(i|τ|) for τ < 0
1906 : ! (recall M_Qνσ(iτ) = M_Qνσ(-iτ) because W^MIC_PQ(iτ) = W^MIC_PQ(-iτ) )
1907 : CALL contract_to_Sigma(t_2c_Gocc, t_3c_x_W, t_2c_Sigma_neg_tau, i_atoms, j_atoms, &
1908 : qs_env, bs_env, occ=.TRUE., vir=.FALSE., &
1909 1266 : can_skip=bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx))
1910 :
1911 : ! Σ_λσ(iτ,k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) G^vir_νµ(i|τ|) for τ > 0
1912 : CALL contract_to_Sigma(t_2c_Gvir, t_3c_x_W, t_2c_Sigma_pos_tau, i_atoms, j_atoms, &
1913 : qs_env, bs_env, occ=.FALSE., vir=.TRUE., &
1914 2550 : can_skip=bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx))
1915 :
1916 : END DO ! j_atoms
1917 : END DO ! i_atoms
1918 :
1919 : ! 4. communicate data tensor t_2c_Sigma (which is local in the subgroup)
1920 : ! to the global dbcsr matrix mat_Sigma_pos/neg_tau (which stores Σ for all iτ)
1921 : CALL local_dbt_to_global_mat(t_2c_Sigma_neg_tau, bs_env%mat_ao_ao_tensor%matrix, &
1922 1284 : mat_Sigma_neg_tau(i_t, ispin)%matrix, bs_env%para_env)
1923 : CALL local_dbt_to_global_mat(t_2c_Sigma_pos_tau, bs_env%mat_ao_ao_tensor%matrix, &
1924 1284 : mat_Sigma_pos_tau(i_t, ispin)%matrix, bs_env%para_env)
1925 :
1926 : CALL write_matrix(mat_Sigma_pos_tau(i_t, ispin)%matrix, read_write_index, &
1927 1284 : bs_env%Sigma_p_name, bs_env%fm_work_mo(1), qs_env)
1928 : CALL write_matrix(mat_Sigma_neg_tau(i_t, ispin)%matrix, read_write_index, &
1929 1284 : bs_env%Sigma_n_name, bs_env%fm_work_mo(1), qs_env)
1930 :
1931 2528 : IF (bs_env%unit_nr > 0) THEN
1932 : WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F10.1,A)') &
1933 642 : 'Computed Σ^c(iτ,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1934 1284 : ', Execution time', m_walltime() - t1, ' s'
1935 : END IF
1936 :
1937 : END DO ! ispin
1938 :
1939 : END DO ! i_t
1940 :
1941 68 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1942 :
1943 : CALL fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
1944 68 : mat_Sigma_pos_tau, mat_Sigma_neg_tau)
1945 :
1946 68 : CALL print_skipping(bs_env)
1947 :
1948 : CALL destroy_mat_Sigma_c(t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
1949 : t_2c_Sigma_pos_tau, t_3c_x_W, fm_W_MIC_time, &
1950 68 : mat_Sigma_neg_tau, mat_Sigma_pos_tau)
1951 :
1952 68 : CALL delete_unnecessary_files(bs_env)
1953 :
1954 68 : CALL timestop(handle)
1955 :
1956 136 : END SUBROUTINE get_Sigma_c
1957 :
1958 : ! **************************************************************************************************
1959 : !> \brief ...
1960 : !> \param bs_env ...
1961 : !> \param t_2c_Gocc ...
1962 : !> \param t_2c_Gvir ...
1963 : !> \param t_2c_W ...
1964 : !> \param t_2c_Sigma_neg_tau ...
1965 : !> \param t_2c_Sigma_pos_tau ...
1966 : !> \param t_3c_x_W ...
1967 : !> \param mat_Sigma_neg_tau ...
1968 : !> \param mat_Sigma_pos_tau ...
1969 : ! **************************************************************************************************
1970 68 : SUBROUTINE create_mat_for_Sigma_c(bs_env, t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
1971 : t_2c_Sigma_pos_tau, t_3c_x_W, &
1972 : mat_Sigma_neg_tau, mat_Sigma_pos_tau)
1973 :
1974 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1975 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_2c_W, &
1976 : t_2c_Sigma_neg_tau, &
1977 : t_2c_Sigma_pos_tau, t_3c_x_W
1978 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
1979 :
1980 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_mat_for_Sigma_c'
1981 :
1982 : INTEGER :: handle, i_t, ispin
1983 :
1984 68 : CALL timeset(routineN, handle)
1985 :
1986 68 : CALL dbt_create(bs_env%t_G, t_2c_Gocc)
1987 68 : CALL dbt_create(bs_env%t_G, t_2c_Gvir)
1988 68 : CALL dbt_create(bs_env%t_W, t_2c_W)
1989 68 : CALL dbt_create(bs_env%t_G, t_2c_Sigma_neg_tau)
1990 68 : CALL dbt_create(bs_env%t_G, t_2c_Sigma_pos_tau)
1991 68 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_W)
1992 :
1993 68 : NULLIFY (mat_Sigma_neg_tau, mat_Sigma_pos_tau)
1994 1754 : ALLOCATE (mat_Sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1995 1754 : ALLOCATE (mat_Sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1996 :
1997 146 : DO ispin = 1, bs_env%n_spin
1998 1550 : DO i_t = 1, bs_env%num_time_freq_points
1999 1404 : ALLOCATE (mat_Sigma_neg_tau(i_t, ispin)%matrix)
2000 1404 : ALLOCATE (mat_Sigma_pos_tau(i_t, ispin)%matrix)
2001 1404 : CALL dbcsr_create(mat_Sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
2002 1482 : CALL dbcsr_create(mat_Sigma_pos_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
2003 : END DO
2004 : END DO
2005 :
2006 68 : CALL timestop(handle)
2007 :
2008 68 : END SUBROUTINE create_mat_for_Sigma_c
2009 :
2010 : ! **************************************************************************************************
2011 : !> \brief ...
2012 : !> \param qs_env ...
2013 : !> \param bs_env ...
2014 : !> \param i_atoms ...
2015 : !> \param j_atoms ...
2016 : !> \param t_3c_x_W ...
2017 : !> \param t_2c_W ...
2018 : ! **************************************************************************************************
2019 1336 : SUBROUTINE compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_W, t_2c_W)
2020 :
2021 : TYPE(qs_environment_type), POINTER :: qs_env
2022 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2023 : INTEGER, DIMENSION(2) :: i_atoms, j_atoms
2024 : TYPE(dbt_type) :: t_3c_x_W, t_2c_W
2025 :
2026 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_and_contract_W'
2027 :
2028 : INTEGER :: handle, RI_intval_idx
2029 : INTEGER(KIND=int_8) :: flop
2030 : INTEGER, DIMENSION(2) :: bounds_P, bounds_Q, RI_atoms
2031 : INTEGER, DIMENSION(2, 2) :: bounds_ao
2032 22712 : TYPE(dbt_type) :: t_3c_for_W, t_3c_x_W_tmp
2033 :
2034 1336 : CALL timeset(routineN, handle)
2035 :
2036 1336 : CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_W_tmp)
2037 1336 : CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_for_W)
2038 :
2039 : ! final layout will be: M_Qνσ(iτ) = sum_P (P|νσ) W^MIC_QP(iτ)
2040 : ! Bounds:
2041 : ! "AO"
2042 : ! -> ν (AO_1 in compute_3c_integrals) bounds from i_atoms and sparse in σ and P
2043 : ! -> σ (AO_2 in compute_3c_integrals) sparse in ν and P
2044 : ! Q bounds from j_atoms
2045 : ! P bounds from inner loop indices and sparse in ν and σ
2046 :
2047 : bounds_Q(1:2) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
2048 4008 : bs_env%i_RI_end_from_atom(j_atoms(2))]
2049 :
2050 2672 : DO RI_intval_idx = 1, bs_env%n_intervals_inner_loop_atoms
2051 4008 : RI_atoms = bs_env%inner_loop_atom_intervals(1:2, RI_intval_idx)
2052 :
2053 : CALL get_bounds_from_atoms(bounds_P, i_atoms, [1, bs_env%n_atom], &
2054 : bs_env%min_RI_idx_from_AO_AO_atom, &
2055 : bs_env%max_RI_idx_from_AO_AO_atom, &
2056 : atoms_3=RI_atoms, &
2057 : indices_3_start=bs_env%i_RI_start_from_atom, &
2058 4008 : indices_3_end=bs_env%i_RI_end_from_atom)
2059 :
2060 : ! σ
2061 : CALL get_bounds_from_atoms(bounds_ao(:, 2), RI_atoms, i_atoms, &
2062 : bs_env%min_AO_idx_from_RI_AO_atom, &
2063 1336 : bs_env%max_AO_idx_from_RI_AO_atom)
2064 : ! ν
2065 : CALL get_bounds_from_atoms(bounds_ao(:, 1), RI_atoms, [1, bs_env%n_atom], &
2066 : bs_env%min_AO_idx_from_RI_AO_atom, &
2067 : bs_env%max_AO_idx_from_RI_AO_atom, &
2068 : atoms_3=i_atoms, &
2069 : indices_3_start=bs_env%i_ao_start_from_atom, &
2070 4008 : indices_3_end=bs_env%i_ao_end_from_atom)
2071 :
2072 1336 : IF (bounds_P(1) > bounds_P(2) .OR. bounds_ao(1, 2) > bounds_ao(2, 2)) THEN
2073 : CYCLE
2074 : END IF
2075 :
2076 : ! 1. compute 3-center integrals (P|µν) ("|": truncated Coulomb operator)
2077 : CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_W, &
2078 1336 : atoms_AO_1=i_atoms, atoms_RI=RI_atoms)
2079 :
2080 : ! 2. tensor operation M_Qνσ(iτ) = sum_P W^MIC_QP(iτ) (P|νσ)
2081 : CALL dbt_contract(alpha=1.0_dp, &
2082 : tensor_1=t_2c_W, &
2083 : tensor_2=t_3c_for_W, &
2084 : beta=1.0_dp, &
2085 : tensor_3=t_3c_x_W_tmp, &
2086 : contract_1=[2], notcontract_1=[1], map_1=[1], &
2087 : contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3], &
2088 : bounds_1=bounds_P, &
2089 : bounds_2=bounds_Q, &
2090 : bounds_3=bounds_ao, &
2091 : flop=flop, &
2092 : move_data=.FALSE., &
2093 2672 : filter_eps=bs_env%eps_filter)
2094 :
2095 : END DO ! RI_atoms
2096 :
2097 : ! 3. reorder tensor
2098 1336 : CALL dbt_copy(t_3c_x_W_tmp, t_3c_x_W, order=[1, 2, 3], move_data=.TRUE.)
2099 :
2100 1336 : CALL dbt_destroy(t_3c_x_W_tmp)
2101 1336 : CALL dbt_destroy(t_3c_for_W)
2102 :
2103 1336 : CALL timestop(handle)
2104 :
2105 1336 : END SUBROUTINE compute_3c_and_contract_W
2106 :
2107 : ! **************************************************************************************************
2108 : !> \brief ...
2109 : !> \param t_2c_G ...
2110 : !> \param t_3c_x_W ...
2111 : !> \param t_2c_Sigma ...
2112 : !> \param i_atoms ...
2113 : !> \param j_atoms ...
2114 : !> \param qs_env ...
2115 : !> \param bs_env ...
2116 : !> \param occ ...
2117 : !> \param vir ...
2118 : !> \param can_skip ...
2119 : ! **************************************************************************************************
2120 2602 : SUBROUTINE contract_to_Sigma(t_2c_G, t_3c_x_W, t_2c_Sigma, i_atoms, j_atoms, qs_env, bs_env, &
2121 : occ, vir, can_skip)
2122 : TYPE(dbt_type) :: t_2c_G, t_3c_x_W, t_2c_Sigma
2123 : INTEGER, DIMENSION(2) :: i_atoms, j_atoms
2124 : TYPE(qs_environment_type), POINTER :: qs_env
2125 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2126 : LOGICAL :: occ, vir
2127 : LOGICAL, OPTIONAL :: can_skip
2128 :
2129 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_to_Sigma'
2130 :
2131 : INTEGER :: handle, inner_loop_atoms_interval_index
2132 : INTEGER(KIND=int_8) :: flop
2133 : INTEGER, DIMENSION(2) :: bounds_lambda, bounds_mu, bounds_nu, &
2134 : bounds_sigma, IL_atoms
2135 : INTEGER, DIMENSION(2, 2) :: bounds_comb
2136 : REAL(KIND=dp) :: sign_Sigma
2137 65050 : TYPE(dbt_type) :: t_3c_for_G, t_3c_x_G, t_3c_x_G_2
2138 :
2139 2602 : CALL timeset(routineN, handle)
2140 :
2141 2602 : CPASSERT(occ .EQV. (.NOT. vir))
2142 2602 : IF (occ) sign_Sigma = -1.0_dp
2143 2602 : IF (vir) sign_Sigma = 1.0_dp
2144 :
2145 2602 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_G)
2146 2602 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_G)
2147 2602 : CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_G_2)
2148 :
2149 : ! Here, in the first step e.g., is computed: N_Qλν = sum_µ (Qλ|µ) G_νµ
2150 : ! Afterwards e.g., is computed: Σ_λσ = sum_νQ M_Qνσ N_Qνλ (after reordering)
2151 : ! Bounds:
2152 : ! "comb" (combined index)
2153 : ! -> Q bounds from j_atoms and sparse in λ
2154 : ! -> λ (AO_1 in compute_3c_integrals) sparse in Q and µ
2155 : ! µ (AO_2 in compute_3c_integrals) bounds from inner loop "IL" indices and sparse in Q and λ
2156 : ! ν bounds from i_atoms
2157 : ! σ sparse in ν
2158 :
2159 : ! ν
2160 : bounds_nu(1:2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
2161 7806 : bs_env%i_ao_end_from_atom(i_atoms(2))]
2162 :
2163 5204 : DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
2164 7806 : IL_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
2165 :
2166 : ! µ
2167 : CALL get_bounds_from_atoms(bounds_mu, j_atoms, [1, bs_env%n_atom], &
2168 : bs_env%min_AO_idx_from_RI_AO_atom, &
2169 : bs_env%max_AO_idx_from_RI_AO_atom, &
2170 : atoms_3=IL_atoms, &
2171 : indices_3_start=bs_env%i_ao_start_from_atom, &
2172 7806 : indices_3_end=bs_env%i_ao_end_from_atom)
2173 :
2174 : ! Q
2175 : CALL get_bounds_from_atoms(bounds_comb(:, 1), IL_atoms, [1, bs_env%n_atom], &
2176 : bs_env%min_RI_idx_from_AO_AO_atom, &
2177 : bs_env%max_RI_idx_from_AO_AO_atom, &
2178 : atoms_3=j_atoms, &
2179 : indices_3_start=bs_env%i_RI_start_from_atom, &
2180 7806 : indices_3_end=bs_env%i_RI_end_from_atom)
2181 :
2182 : ! λ
2183 : CALL get_bounds_from_atoms(bounds_comb(:, 2), j_atoms, IL_atoms, &
2184 : bs_env%min_AO_idx_from_RI_AO_atom, &
2185 2602 : bs_env%max_AO_idx_from_RI_AO_atom)
2186 :
2187 2602 : IF (bounds_mu(1) > bounds_mu(2) .OR. bounds_comb(1, 1) > bounds_comb(2, 1) .OR. &
2188 : bounds_comb(1, 2) > bounds_comb(2, 2)) THEN
2189 : CYCLE
2190 : END IF
2191 :
2192 : CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_G, &
2193 2602 : atoms_RI=j_atoms, atoms_AO_2=IL_atoms)
2194 :
2195 : CALL dbt_contract(alpha=1.0_dp, &
2196 : tensor_1=t_2c_G, &
2197 : tensor_2=t_3c_for_G, &
2198 : beta=1.0_dp, &
2199 : tensor_3=t_3c_x_G, &
2200 : contract_1=[2], notcontract_1=[1], map_1=[3], &
2201 : contract_2=[3], notcontract_2=[1, 2], map_2=[1, 2], &
2202 : bounds_1=bounds_mu, &
2203 : bounds_2=bounds_nu, &
2204 : bounds_3=bounds_comb, &
2205 : flop=flop, &
2206 : move_data=.FALSE., &
2207 5204 : filter_eps=bs_env%eps_filter)
2208 : END DO ! IL_atoms
2209 :
2210 : ! Reordering: N_Qλν -> N_Qνλ
2211 2602 : CALL dbt_copy(t_3c_x_G, t_3c_x_G_2, order=[1, 3, 2], move_data=.TRUE.)
2212 :
2213 : ! Here, the last contraction is done, e.g., Σ_λσ = sum_νQ M_Qνσ N_Qνλ
2214 : ! Bounds as above, new "comb" with upper ingredients
2215 : bounds_comb(1:2, 1) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
2216 7806 : bs_env%i_RI_end_from_atom(j_atoms(2))]
2217 7806 : bounds_comb(1:2, 2) = bounds_nu(1:2)
2218 :
2219 : CALL get_bounds_from_atoms(bounds_lambda, j_atoms, [1, bs_env%n_atom], &
2220 : bs_env%min_AO_idx_from_RI_AO_atom, &
2221 7806 : bs_env%max_AO_idx_from_RI_AO_atom)
2222 : CALL get_bounds_from_atoms(bounds_sigma, [1, bs_env%n_atom], i_atoms, &
2223 : bs_env%min_AO_idx_from_RI_AO_atom, &
2224 7806 : bs_env%max_AO_idx_from_RI_AO_atom)
2225 :
2226 2602 : IF (bounds_sigma(1) > bounds_sigma(2) .OR. bounds_lambda(1) > bounds_lambda(2)) THEN
2227 0 : flop = 0_int_8
2228 : ELSE
2229 : CALL dbt_contract(alpha=sign_Sigma, &
2230 : tensor_1=t_3c_x_W, &
2231 : tensor_2=t_3c_x_G_2, &
2232 : beta=1.0_dp, &
2233 : tensor_3=t_2c_Sigma, &
2234 : contract_1=[1, 2], notcontract_1=[3], map_1=[1], &
2235 : contract_2=[1, 2], notcontract_2=[3], map_2=[2], &
2236 : bounds_1=bounds_comb, &
2237 : bounds_2=bounds_sigma, &
2238 : bounds_3=bounds_lambda, &
2239 2602 : filter_eps=bs_env%eps_filter, move_data=.FALSE., flop=flop)
2240 : END IF
2241 :
2242 2602 : IF (PRESENT(can_skip)) THEN
2243 2532 : IF (flop == 0_int_8) can_skip = .TRUE.
2244 : END IF
2245 :
2246 2602 : CALL dbt_destroy(t_3c_for_G)
2247 2602 : CALL dbt_destroy(t_3c_x_G)
2248 2602 : CALL dbt_destroy(t_3c_x_G_2)
2249 :
2250 2602 : CALL timestop(handle)
2251 :
2252 2602 : END SUBROUTINE contract_to_Sigma
2253 :
2254 : ! **************************************************************************************************
2255 : !> \brief ...
2256 : !> \param fm_Sigma_c_Gamma_time ...
2257 : !> \param bs_env ...
2258 : !> \param mat_Sigma_pos_tau ...
2259 : !> \param mat_Sigma_neg_tau ...
2260 : ! **************************************************************************************************
2261 118 : SUBROUTINE fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
2262 : mat_Sigma_pos_tau, mat_Sigma_neg_tau)
2263 :
2264 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
2265 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2266 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_pos_tau, mat_Sigma_neg_tau
2267 :
2268 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_fm_Sigma_c_Gamma_time'
2269 :
2270 : INTEGER :: handle, i_t, ispin, pos_neg
2271 :
2272 118 : CALL timeset(routineN, handle)
2273 :
2274 5280 : ALLOCATE (fm_Sigma_c_Gamma_time(bs_env%num_time_freq_points, 2, bs_env%n_spin))
2275 258 : DO ispin = 1, bs_env%n_spin
2276 2452 : DO i_t = 1, bs_env%num_time_freq_points
2277 6582 : DO pos_neg = 1, 2
2278 : CALL cp_fm_create(fm_Sigma_c_Gamma_time(i_t, pos_neg, ispin), &
2279 6582 : bs_env%fm_s_Gamma%matrix_struct)
2280 : END DO
2281 : CALL copy_dbcsr_to_fm(mat_Sigma_pos_tau(i_t, ispin)%matrix, &
2282 2194 : fm_Sigma_c_Gamma_time(i_t, 1, ispin))
2283 : CALL copy_dbcsr_to_fm(mat_Sigma_neg_tau(i_t, ispin)%matrix, &
2284 2334 : fm_Sigma_c_Gamma_time(i_t, 2, ispin))
2285 : END DO
2286 : END DO
2287 :
2288 118 : CALL timestop(handle)
2289 :
2290 118 : END SUBROUTINE fill_fm_Sigma_c_Gamma_time
2291 :
2292 : ! **************************************************************************************************
2293 : !> \brief ...
2294 : !> \param bs_env ...
2295 : ! **************************************************************************************************
2296 68 : SUBROUTINE print_skipping(bs_env)
2297 :
2298 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2299 :
2300 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_skipping'
2301 :
2302 : INTEGER :: handle, n_pairs
2303 :
2304 68 : CALL timeset(routineN, handle)
2305 :
2306 68 : n_pairs = bs_env%n_intervals_i*bs_env%n_intervals_j*bs_env%n_spin
2307 :
2308 68 : CALL bs_env%para_env_tensor%sum(bs_env%n_skip_sigma)
2309 68 : CALL bs_env%para_env_tensor%sum(bs_env%n_skip_chi)
2310 68 : CALL bs_env%para_env_tensor%sum(n_pairs)
2311 :
2312 68 : IF (bs_env%unit_nr > 0) THEN
2313 : WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
2314 34 : 'Sparsity of Σ^c(iτ,k=0): Percentage of skipped atom pairs:', &
2315 68 : REAL(100*bs_env%n_skip_sigma, KIND=dp)/REAL(n_pairs, KIND=dp), ' %'
2316 : WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
2317 34 : 'Sparsity of χ(iτ,k=0): Percentage of skipped atom pairs:', &
2318 68 : REAL(100*bs_env%n_skip_chi, KIND=dp)/REAL(n_pairs, KIND=dp), ' %'
2319 : END IF
2320 :
2321 68 : CALL timestop(handle)
2322 :
2323 68 : END SUBROUTINE print_skipping
2324 :
2325 : ! **************************************************************************************************
2326 : !> \brief ...
2327 : !> \param t_2c_Gocc ...
2328 : !> \param t_2c_Gvir ...
2329 : !> \param t_2c_W ...
2330 : !> \param t_2c_Sigma_neg_tau ...
2331 : !> \param t_2c_Sigma_pos_tau ...
2332 : !> \param t_3c_x_W ...
2333 : !> \param fm_W_MIC_time ...
2334 : !> \param mat_Sigma_neg_tau ...
2335 : !> \param mat_Sigma_pos_tau ...
2336 : ! **************************************************************************************************
2337 68 : SUBROUTINE destroy_mat_Sigma_c(t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
2338 : t_2c_Sigma_pos_tau, t_3c_x_W, fm_W_MIC_time, &
2339 : mat_Sigma_neg_tau, mat_Sigma_pos_tau)
2340 :
2341 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_2c_W, &
2342 : t_2c_Sigma_neg_tau, &
2343 : t_2c_Sigma_pos_tau, t_3c_x_W
2344 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
2345 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
2346 :
2347 : CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_mat_Sigma_c'
2348 :
2349 : INTEGER :: handle
2350 :
2351 68 : CALL timeset(routineN, handle)
2352 :
2353 68 : CALL dbt_destroy(t_2c_Gocc)
2354 68 : CALL dbt_destroy(t_2c_Gvir)
2355 68 : CALL dbt_destroy(t_2c_W)
2356 68 : CALL dbt_destroy(t_2c_Sigma_neg_tau)
2357 68 : CALL dbt_destroy(t_2c_Sigma_pos_tau)
2358 68 : CALL dbt_destroy(t_3c_x_W)
2359 68 : CALL cp_fm_release(fm_W_MIC_time)
2360 68 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_neg_tau)
2361 68 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_pos_tau)
2362 :
2363 68 : CALL timestop(handle)
2364 :
2365 68 : END SUBROUTINE destroy_mat_Sigma_c
2366 :
2367 : ! **************************************************************************************************
2368 : !> \brief ...
2369 : !> \param bs_env ...
2370 : ! **************************************************************************************************
2371 104 : SUBROUTINE delete_unnecessary_files(bs_env)
2372 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2373 :
2374 : CHARACTER(LEN=*), PARAMETER :: routineN = 'delete_unnecessary_files'
2375 :
2376 : CHARACTER(LEN=default_path_length) :: f_chi, f_W_t, prefix
2377 : INTEGER :: handle, i_t
2378 :
2379 104 : CALL timeset(routineN, handle)
2380 :
2381 104 : prefix = bs_env%prefix
2382 :
2383 1858 : DO i_t = 1, bs_env%num_time_freq_points
2384 :
2385 1754 : IF (i_t < 10) THEN
2386 924 : WRITE (f_chi, '(3A,I1,A)') TRIM(prefix), bs_env%chi_name, "_00", i_t, ".matrix"
2387 924 : WRITE (f_W_t, '(3A,I1,A)') TRIM(prefix), bs_env%W_time_name, "_00", i_t, ".matrix"
2388 830 : ELSE IF (i_t < 100) THEN
2389 830 : WRITE (f_chi, '(3A,I2,A)') TRIM(prefix), bs_env%chi_name, "_0", i_t, ".matrix"
2390 830 : WRITE (f_W_t, '(3A,I2,A)') TRIM(prefix), bs_env%W_time_name, "_0", i_t, ".matrix"
2391 : ELSE
2392 0 : CPABORT('Please implement more than 99 time/frequency points.')
2393 : END IF
2394 :
2395 1754 : CALL safe_delete(f_chi, bs_env)
2396 1858 : CALL safe_delete(f_W_t, bs_env)
2397 :
2398 : END DO
2399 :
2400 104 : CALL timestop(handle)
2401 :
2402 104 : END SUBROUTINE delete_unnecessary_files
2403 :
2404 : ! **************************************************************************************************
2405 : !> \brief ...
2406 : !> \param filename ...
2407 : !> \param bs_env ...
2408 : ! **************************************************************************************************
2409 3508 : SUBROUTINE safe_delete(filename, bs_env)
2410 : CHARACTER(LEN=*) :: filename
2411 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2412 :
2413 : CHARACTER(LEN=*), PARAMETER :: routineN = 'safe_delete'
2414 :
2415 : INTEGER :: handle
2416 : LOGICAL :: file_exists
2417 :
2418 3508 : CALL timeset(routineN, handle)
2419 :
2420 3508 : IF (bs_env%para_env%mepos == 0) THEN
2421 :
2422 1754 : INQUIRE (file=TRIM(filename), exist=file_exists)
2423 1754 : IF (file_exists) CALL mp_file_delete(TRIM(filename))
2424 :
2425 : END IF
2426 :
2427 3508 : CALL timestop(handle)
2428 :
2429 3508 : END SUBROUTINE safe_delete
2430 :
2431 : ! **************************************************************************************************
2432 : !> \brief ...
2433 : !> \param bs_env ...
2434 : !> \param qs_env ...
2435 : !> \param fm_Sigma_x_Gamma ...
2436 : !> \param fm_Sigma_c_Gamma_time ...
2437 : ! **************************************************************************************************
2438 68 : SUBROUTINE compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
2439 :
2440 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2441 : TYPE(qs_environment_type), POINTER :: qs_env
2442 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
2443 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
2444 :
2445 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies'
2446 :
2447 : INTEGER :: handle, ikp, ispin, j_t
2448 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Sigma_x_ikp_n, V_xc_ikp_n
2449 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Sigma_c_ikp_n_freq, Sigma_c_ikp_n_time
2450 : TYPE(cp_cfm_type) :: cfm_ks_ikp, cfm_mos_ikp, cfm_s_ikp, &
2451 : cfm_Sigma_x_ikp, cfm_work_ikp
2452 :
2453 68 : CALL timeset(routineN, handle)
2454 :
2455 68 : CALL cp_cfm_create(cfm_mos_ikp, bs_env%fm_s_Gamma%matrix_struct)
2456 68 : CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_s_Gamma%matrix_struct)
2457 : ! JW TODO: fully distribute these arrays at given time; also eigenvalues in bs_env
2458 272 : ALLOCATE (V_xc_ikp_n(bs_env%n_ao), Sigma_x_ikp_n(bs_env%n_ao))
2459 340 : ALLOCATE (Sigma_c_ikp_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
2460 204 : ALLOCATE (Sigma_c_ikp_n_freq(bs_env%n_ao, bs_env%num_time_freq_points, 2))
2461 :
2462 146 : DO ispin = 1, bs_env%n_spin
2463 :
2464 236 : DO ikp = 1, bs_env%nkp_bs_and_DOS
2465 :
2466 : ! 1. get H^KS_µν(k_i) from H^KS_µν(k=0)
2467 : CALL cfm_ikp_from_fm_Gamma(cfm_ks_ikp, bs_env%fm_ks_Gamma(ispin), &
2468 90 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2469 :
2470 : ! 2. get S_µν(k_i) from S_µν(k=0)
2471 : CALL cfm_ikp_from_fm_Gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
2472 90 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2473 :
2474 : ! 3. Diagonalize (Roothaan-Hall): H_KS(k_i)*C(k_i) = S(k_i)*C(k_i)*ϵ(k_i)
2475 : CALL cp_cfm_geeig(cfm_ks_ikp, cfm_s_ikp, cfm_mos_ikp, &
2476 90 : bs_env%eigenval_scf(:, ikp, ispin), cfm_work_ikp)
2477 :
2478 : ! 4. V^xc_µν(k=0) -> V^xc_µν(k_i) -> V^xc_nn(k_i)
2479 : CALL to_ikp_and_mo(V_xc_ikp_n, bs_env%fm_V_xc_Gamma(ispin), &
2480 90 : ikp, qs_env, bs_env, cfm_mos_ikp)
2481 :
2482 : ! 5. Σ^x_µν(k=0) -> Σ^x_µν(k_i) -> Σ^x_nn(k_i)
2483 : CALL to_ikp_and_mo(Sigma_x_ikp_n, fm_Sigma_x_Gamma(ispin), &
2484 90 : ikp, qs_env, bs_env, cfm_mos_ikp)
2485 :
2486 : ! 6. Σ^c_µν(k=0,+/-i|τ_j|) -> Σ^c_µν(k_i,+/-i|τ_j|) -> Σ^c_nn(k_i,+/-i|τ_j|)
2487 1598 : DO j_t = 1, bs_env%num_time_freq_points
2488 : CALL to_ikp_and_mo(Sigma_c_ikp_n_time(:, j_t, 1), &
2489 : fm_Sigma_c_Gamma_time(j_t, 1, ispin), &
2490 1508 : ikp, qs_env, bs_env, cfm_mos_ikp)
2491 : CALL to_ikp_and_mo(Sigma_c_ikp_n_time(:, j_t, 2), &
2492 : fm_Sigma_c_Gamma_time(j_t, 2, ispin), &
2493 1598 : ikp, qs_env, bs_env, cfm_mos_ikp)
2494 : END DO
2495 :
2496 : ! 7. Σ^c_nn(k_i,iτ) -> Σ^c_nn(k_i,iω)
2497 90 : CALL time_to_freq(bs_env, Sigma_c_ikp_n_time, Sigma_c_ikp_n_freq, ispin)
2498 :
2499 : ! 8. Analytic continuation Σ^c_nn(k_i,iω) -> Σ^c_nn(k_i,ϵ) and
2500 : ! ϵ_nk_i^GW = ϵ_nk_i^DFT + Σ^c_nn(k_i,ϵ) + Σ^x_nn(k_i) - v^xc_nn(k_i)
2501 : CALL analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_xc_ikp_n, &
2502 168 : bs_env%eigenval_scf(:, ikp, ispin), ikp, ispin)
2503 :
2504 : END DO ! ikp_DOS
2505 :
2506 : END DO ! ispin
2507 :
2508 68 : CALL get_all_VBM_CBM_bandgaps(bs_env)
2509 :
2510 : ! Σ^x is releases here in case of G0W0
2511 68 : IF (bs_env%gw_flavour == G0W0) CALL cp_fm_release(fm_Sigma_x_Gamma)
2512 68 : CALL cp_fm_release(fm_Sigma_c_Gamma_time)
2513 68 : CALL cp_cfm_release(cfm_ks_ikp)
2514 68 : CALL cp_cfm_release(cfm_s_ikp)
2515 68 : CALL cp_cfm_release(cfm_mos_ikp)
2516 68 : CALL cp_cfm_release(cfm_work_ikp)
2517 68 : CALL cp_cfm_release(cfm_Sigma_x_ikp)
2518 :
2519 68 : CALL timestop(handle)
2520 :
2521 136 : END SUBROUTINE compute_QP_energies
2522 :
2523 : ! **************************************************************************************************
2524 : !> \brief ...
2525 : !> \param array_ikp_n ...
2526 : !> \param fm_Gamma ...
2527 : !> \param ikp ...
2528 : !> \param qs_env ...
2529 : !> \param bs_env ...
2530 : !> \param cfm_mos_ikp ...
2531 : ! **************************************************************************************************
2532 3196 : SUBROUTINE to_ikp_and_mo(array_ikp_n, fm_Gamma, ikp, qs_env, bs_env, cfm_mos_ikp)
2533 :
2534 : REAL(KIND=dp), DIMENSION(:) :: array_ikp_n
2535 : TYPE(cp_fm_type) :: fm_Gamma
2536 : INTEGER :: ikp
2537 : TYPE(qs_environment_type), POINTER :: qs_env
2538 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2539 : TYPE(cp_cfm_type) :: cfm_mos_ikp
2540 :
2541 : CHARACTER(LEN=*), PARAMETER :: routineN = 'to_ikp_and_mo'
2542 :
2543 : INTEGER :: handle
2544 : TYPE(cp_fm_type) :: fm_ikp_mo_re
2545 :
2546 3196 : CALL timeset(routineN, handle)
2547 :
2548 3196 : CALL cp_fm_create(fm_ikp_mo_re, fm_Gamma%matrix_struct)
2549 :
2550 3196 : CALL fm_Gamma_ao_to_cfm_ikp_mo(fm_Gamma, fm_ikp_mo_re, ikp, qs_env, bs_env, cfm_mos_ikp)
2551 :
2552 3196 : CALL cp_fm_get_diag(fm_ikp_mo_re, array_ikp_n)
2553 :
2554 3196 : CALL cp_fm_release(fm_ikp_mo_re)
2555 :
2556 3196 : CALL timestop(handle)
2557 :
2558 3196 : END SUBROUTINE to_ikp_and_mo
2559 :
2560 : ! **************************************************************************************************
2561 : !> \brief ...
2562 : !> \param fm_Gamma ...
2563 : !> \param fm_ikp_mo_re ...
2564 : !> \param ikp ...
2565 : !> \param qs_env ...
2566 : !> \param bs_env ...
2567 : !> \param cfm_mos_ikp ...
2568 : ! **************************************************************************************************
2569 12784 : SUBROUTINE fm_Gamma_ao_to_cfm_ikp_mo(fm_Gamma, fm_ikp_mo_re, ikp, qs_env, bs_env, cfm_mos_ikp)
2570 : TYPE(cp_fm_type) :: fm_Gamma, fm_ikp_mo_re
2571 : INTEGER :: ikp
2572 : TYPE(qs_environment_type), POINTER :: qs_env
2573 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2574 : TYPE(cp_cfm_type) :: cfm_mos_ikp
2575 :
2576 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_Gamma_ao_to_cfm_ikp_mo'
2577 :
2578 : INTEGER :: handle, nmo
2579 : TYPE(cp_cfm_type) :: cfm_ikp_ao, cfm_ikp_mo, cfm_tmp
2580 :
2581 3196 : CALL timeset(routineN, handle)
2582 :
2583 3196 : CALL cp_cfm_create(cfm_ikp_ao, fm_Gamma%matrix_struct)
2584 3196 : CALL cp_cfm_create(cfm_ikp_mo, fm_Gamma%matrix_struct)
2585 3196 : CALL cp_cfm_create(cfm_tmp, fm_Gamma%matrix_struct)
2586 :
2587 : ! get cfm_µν(k_i) from fm_µν(k=0)
2588 3196 : CALL cfm_ikp_from_fm_Gamma(cfm_ikp_ao, fm_Gamma, ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2589 :
2590 3196 : nmo = bs_env%n_ao
2591 3196 : CALL parallel_gemm('N', 'N', nmo, nmo, nmo, z_one, cfm_ikp_ao, cfm_mos_ikp, z_zero, cfm_tmp)
2592 3196 : CALL parallel_gemm('C', 'N', nmo, nmo, nmo, z_one, cfm_mos_ikp, cfm_tmp, z_zero, cfm_ikp_mo)
2593 :
2594 3196 : CALL cp_cfm_to_fm(cfm_ikp_mo, fm_ikp_mo_re)
2595 :
2596 3196 : CALL cp_cfm_release(cfm_ikp_mo)
2597 3196 : CALL cp_cfm_release(cfm_ikp_ao)
2598 3196 : CALL cp_cfm_release(cfm_tmp)
2599 :
2600 3196 : CALL timestop(handle)
2601 :
2602 3196 : END SUBROUTINE fm_Gamma_ao_to_cfm_ikp_mo
2603 :
2604 : ! **************************************************************************************************
2605 : !> \brief Computes bounds (AO or RI) for given atom intervals atoms_1 and atoms_2 from indices_min
2606 : !> and indices_max and returns them in bounds_out.
2607 : !> In case, atoms_3 and indices_3 are given, the bounds are computed as the intersection
2608 : !> \param bounds_out Bounds to be computed
2609 : !> \param atoms_1 First atom interval
2610 : !> \param atoms_2 Second atom interval
2611 : !> \param indices_min Minimum indices for each atom pair (typically from bs_env,
2612 : !> computed in get_i_j_atom_ranges in gw_utils.F, e.g. bs_env%min_RI_idx_from_AO_AO_atom)
2613 : !> \param indices_max Maximum indices for each atom pair (typically from bs_env,
2614 : !> computed in get_i_j_atom_ranges in gw_utils.F)
2615 : !> \param atoms_3 (Optional) Third atom interval for intersection
2616 : !> \param indices_3_start (Optional) Indices for third atom interval for intersection
2617 : !> \param indices_3_end (Optional) Indices for third atom interval for intersection
2618 : ! **************************************************************************************************
2619 27178 : SUBROUTINE get_bounds_from_atoms(bounds_out, atoms_1, atoms_2, indices_min, indices_max, &
2620 27178 : atoms_3, indices_3_start, indices_3_end)
2621 :
2622 : INTEGER, DIMENSION(2), INTENT(OUT) :: bounds_out
2623 : INTEGER, DIMENSION(2), INTENT(IN) :: atoms_1, atoms_2
2624 : INTEGER, DIMENSION(:, :) :: indices_min, indices_max
2625 : INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: atoms_3
2626 : INTEGER, DIMENSION(:), OPTIONAL :: indices_3_start, indices_3_end
2627 :
2628 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_bounds_from_atoms'
2629 :
2630 : INTEGER :: handle, i_at, j_at
2631 :
2632 27178 : CALL timeset(routineN, handle)
2633 27178 : bounds_out(1) = HUGE(0)
2634 27178 : bounds_out(2) = -1
2635 : !Loop over all atoms in the two intervals and find min/max indices
2636 82070 : DO i_at = atoms_1(1), atoms_1(2)
2637 193462 : DO j_at = atoms_2(1), atoms_2(2)
2638 111392 : bounds_out(1) = MIN(bounds_out(1), indices_min(i_at, j_at))
2639 166284 : bounds_out(2) = MAX(bounds_out(2), indices_max(i_at, j_at))
2640 : END DO
2641 : END DO
2642 :
2643 27178 : IF (PRESENT(atoms_3) .AND. PRESENT(indices_3_start) .AND. PRESENT(indices_3_end)) THEN
2644 12956 : bounds_out(1) = MAX(bounds_out(1), indices_3_start(atoms_3(1)))
2645 12956 : bounds_out(2) = MIN(bounds_out(2), indices_3_end(atoms_3(2)))
2646 : END IF
2647 :
2648 27178 : CALL timestop(handle)
2649 :
2650 27178 : END SUBROUTINE get_bounds_from_atoms
2651 :
2652 : END MODULE gw_large_cell_gamma
|
