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