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