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