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