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