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