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