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