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