(git:374b731)
Loading...
Searching...
No Matches
gw_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE 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
100CONTAINS
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 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_sigma_x_gamma, fm_w_mic_time
119 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_sigma_c_gamma_time
120 TYPE(mp_para_env_type), POINTER :: para_env
121
122 CALL timeset(routinen, handle)
123
124 CALL get_qs_env(qs_env, para_env=para_env)
125
126 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 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 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 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 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 CALL compute_qp_energies(bs_env, qs_env, fm_sigma_x_gamma, fm_sigma_c_gamma_time)
146
147 CALL de_init_bs_env(bs_env)
148
149 CALL timestop(handle)
150
151 END SUBROUTINE gw
152
153! **************************************************************************************************
154!> \brief ...
155!> \param bs_env ...
156!> \param qs_env ...
157!> \param mat_chi_Gamma_tau ...
158! **************************************************************************************************
159 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 TYPE(dbt_type) :: t_2c_gocc, t_2c_gvir, t_3c_for_gocc, &
171 t_3c_for_gvir, t_3c_x_gocc, &
172 t_3c_x_gocc_2, t_3c_x_gvir, &
173 t_3c_x_gvir_2
174
175 CALL timeset(routinen, handle)
176
177 DO i_t = 1, bs_env%num_time_freq_points
178
179 bs_env%t1 = m_walltime()
180
181 IF (bs_env%read_chi(i_t)) THEN
182
183 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 keep_sparsity=.false.)
186
187 IF (bs_env%unit_nr > 0) THEN
188 WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F7.1,A)') &
189 χτ'Read (i,k=0) from file for time point ', i_t, ' /', &
190 bs_env%num_time_freq_points, &
191 ', Execution time', m_walltime() - bs_env%t1, ' s'
192 END IF
193
194 cycle
195
196 END IF
197
198 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 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 tau = bs_env%imag_time_points(i_t)
208
209 DO ispin = 1, bs_env%n_spin
210 CALL g_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.true., vir=.false.)
211 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 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 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 DO i_intval_idx = 1, bs_env%n_intervals_i
222 DO j_intval_idx = 1, bs_env%n_intervals_j
223 i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
224 j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
225
226 DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
227
228 il_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
229
230 CALL check_dist(i_atoms, il_atoms, qs_env, bs_env, dist_too_long_i)
231 CALL check_dist(j_atoms, il_atoms, qs_env, bs_env, dist_too_long_j)
232 IF (dist_too_long_i .OR. dist_too_long_j) cycle
233
234 ! 2. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
235 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 j_atoms, i_atoms, il_atoms)
240
241 ! 4. compute 3-center integrals (σλ|Q) ("|": truncated Coulomb operator)
242 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 i_atoms, j_atoms, il_atoms)
247
248 END DO ! IL_atoms
249
250 ! 6. reorder tensors
251 CALL dbt_copy(t_3c_x_gocc, t_3c_x_gocc_2, move_data=.true., order=[1, 3, 2])
252 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 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 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 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 t_3c_x_gocc, t_3c_x_gvir, t_3c_x_gocc_2, t_3c_x_gvir_2)
277
278 IF (bs_env%unit_nr > 0) THEN
279 WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
280 χτ'Computed (i,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
281 ', Execution time', m_walltime() - bs_env%t1, ' s'
282 END IF
283
284 END DO ! i_t
285
286 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
287
288 CALL timestop(handle)
289
290 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 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 CALL timeset(routinen, handle)
311
312 unit_nr = -1
313 IF (bs_env%para_env%is_source()) THEN
314
315 IF (idx < 10) THEN
316 WRITE (f_chi, '(3A,I1,A)') trim(bs_env%prefix), trim(mat_name), "_0", idx, ".matrix"
317 ELSE IF (idx < 100) THEN
318 WRITE (f_chi, '(3A,I2,A)') trim(bs_env%prefix), trim(mat_name), "_", idx, ".matrix"
319 ELSE
320 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 file_position="REWIND", file_status="OLD", unit_number=unit_nr)
325
326 END IF
327
328 CALL cp_fm_read_unformatted(fm, unit_nr)
329
330 IF (bs_env%para_env%is_source()) CALL close_file(unit_number=unit_nr)
331
332 CALL timestop(handle)
333
334 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 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 CALL timeset(routinen, handle)
362
363 CALL dbt_create(bs_env%t_G, t_2c_gocc)
364 CALL dbt_create(bs_env%t_G, t_2c_gvir)
365 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_gocc)
366 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_gvir)
367 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_gocc)
368 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_gvir)
369 CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_gocc_2)
370 CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_gvir_2)
371
372 CALL timestop(handle)
373
374 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 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 CALL timeset(routinen, handle)
399
400 CALL dbt_destroy(t_2c_gocc)
401 CALL dbt_destroy(t_2c_gvir)
402 CALL dbt_destroy(t_3c_for_gocc)
403 CALL dbt_destroy(t_3c_for_gvir)
404 CALL dbt_destroy(t_3c_x_gocc)
405 CALL dbt_destroy(t_3c_x_gvir)
406 CALL dbt_destroy(t_3c_x_gocc_2)
407 CALL dbt_destroy(t_3c_x_gvir_2)
408
409 CALL timestop(handle)
410
411 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices
435 REAL(kind=dp) :: tau_e
436
437 CALL timeset(routinen, handle)
438
439 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 col_indices=col_indices)
445
446 n_mo = bs_env%n_ao
447 homo = bs_env%n_occ(ispin)
448
449 CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(ispin), bs_env%fm_work_mo(1))
450
451 DO i_row_local = 1, nrow_local
452 DO j_col_local = 1, ncol_local
453
454 j_col = col_indices(j_col_local)
455
456 tau_e = abs(tau*0.5_dp*(bs_env%eigenval_scf_Gamma(j_col, ispin) - bs_env%e_fermi(ispin)))
457
458 IF (tau_e < bs_env%stabilize_exp) THEN
459 bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = &
460 bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local)*exp(-tau_e)
461 ELSE
462 bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = 0.0_dp
463 END IF
464
465 IF ((occ .AND. j_col > homo) .OR. (vir .AND. j_col <= homo)) THEN
466 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 beta=0.0_dp, matrix_c=fm_ggamma)
475
476 CALL timestop(handle)
477
478 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 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 TYPE(dbt_type) :: tensor_tmp
501
502 CALL timeset(routinen, handle)
503
504 CALL dbt_clear(tensor)
505 CALL copy_fm_to_dbcsr(fm_global, mat_global, keep_sparsity=.false.)
506 CALL dbcsr_filter(mat_global, bs_env%eps_filter)
507 IF (PRESENT(atom_ranges)) THEN
508 CALL global_matrix_to_local_matrix(mat_global, mat_local, bs_env%para_env, &
509 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 bs_env%para_env_tensor%num_pe)
513 END IF
514 CALL dbt_create(mat_local, tensor_tmp)
515 CALL dbt_copy_matrix_to_tensor(mat_local, tensor_tmp)
516 CALL dbt_copy(tensor_tmp, tensor, move_data=.true.)
517 CALL dbt_destroy(tensor_tmp)
518 CALL dbcsr_set(mat_local, 0.0_dp)
519 CALL dbcsr_filter(mat_local, 1.0_dp)
520
521 CALL timestop(handle)
522
523 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 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 CALL timeset(routinen, handle)
543
544 CALL dbt_copy_tensor_to_matrix(tensor, mat_tensor)
545 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 CALL para_env%sync()
549 CALL local_matrix_to_global_matrix(mat_tensor, mat_global, para_env)
550
551 CALL timestop(handle)
552
553 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 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 CALL timeset(routinen, handle)
575
576 CALL cp_fm_set_all(fm, 0.0_dp)
577
578 CALL copy_dbcsr_to_fm(matrix, fm)
579
580 CALL fm_write(fm, matrix_index, matrix_name, qs_env)
581
582 CALL timestop(handle)
583
584 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 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 CALL timeset(routinen, handle)
608
609 CALL get_qs_env(qs_env, input=input)
610
611 logger => cp_get_default_logger()
612
613 IF (btest(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
614
615 IF (matrix_index < 10) THEN
616 WRITE (filename, '(3A,I1)') "RESTART_", matrix_name, "_0", matrix_index
617 ELSE IF (matrix_index < 100) THEN
618 WRITE (filename, '(3A,I2)') "RESTART_", matrix_name, "_", matrix_index
619 ELSE
620 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 file_position="REWIND", file_action="WRITE")
626
627 CALL cp_fm_write_unformatted(fm, unit_nr)
628 IF (unit_nr > 0) THEN
629 CALL close_file(unit_nr)
630 END IF
631 END IF
632
633 CALL timestop(handle)
634
635 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 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 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_array
657
658 CALL timeset(routinen, handle)
659
660 ! free memory (not clear whether memory has been freed previously)
661 CALL dbt_clear(t_3c)
662
663 ALLOCATE (t_3c_array(1, 1))
664 CALL dbt_create(t_3c, t_3c_array(1, 1))
665
666 IF (PRESENT(atoms_ao_1)) THEN
667 my_atoms_ao_1 = atoms_ao_1
668 ELSE
669 my_atoms_ao_1 = [1, bs_env%n_atom]
670 END IF
671 IF (PRESENT(atoms_ao_2)) THEN
672 my_atoms_ao_2 = atoms_ao_2
673 ELSE
674 my_atoms_ao_2 = [1, bs_env%n_atom]
675 END IF
676 IF (PRESENT(atoms_ri)) THEN
677 my_atoms_ri = atoms_ri
678 ELSE
679 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 desymmetrize=.false.)
695
696 CALL dbt_copy(t_3c_array(1, 1), t_3c, move_data=.true.)
697
698 CALL dbt_destroy(t_3c_array(1, 1))
699 DEALLOCATE (t_3c_array)
700
701 CALL timestop(handle)
702
703 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 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 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 bs_env%i_ao_end_from_atom(atoms_il(2))]
734 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 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 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 filter_eps=bs_env%eps_filter)
751
752 CALL dbt_clear(t_3c_for_g)
753
754 CALL timestop(handle)
755
756 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 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
779
780 CALL timeset(routinen, handle)
781
782 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
783
784 min_dist_ao_atoms = 1.0e5_dp
785 DO atom_1 = atoms_1(1), atoms_1(2)
786 DO atom_2 = atoms_2(1), atoms_2(2)
787
788 rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
789
790 abs_rab = sqrt(rab(1)**2 + rab(2)**2 + rab(3)**2)
791
792 min_dist_ao_atoms = min(min_dist_ao_atoms, abs_rab)
793
794 END DO
795 END DO
796
797 dist_too_long = (min_dist_ao_atoms > bs_env%max_dist_AO_atoms)
798
799 CALL timestop(handle)
800
801 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 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 CALL timeset(routinen, handle)
821
822 IF (bs_env%all_W_exist) THEN
823 CALL read_w_mic_time(bs_env, mat_chi_gamma_tau, fm_w_mic_time)
824 ELSE
825 CALL compute_w_mic(bs_env, qs_env, mat_chi_gamma_tau, fm_w_mic_time)
826 END IF
827
828 CALL timestop(handle)
829
830 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 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 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
850 TYPE(cell_type), POINTER :: cell
851 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_v_kp
852 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
853 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
854
855 CALL timeset(routinen, handle)
856
857 nkp_chi_eps_w_batch = bs_env%nkp_chi_eps_W_batch
858
859 ikp_start = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + 1
860 ikp_end = min(ikp_batch*bs_env%nkp_chi_eps_W_batch, bs_env%kpoints_chi_eps_W%nkp)
861
862 NULLIFY (mat_v_kp)
863 ALLOCATE (mat_v_kp(ikp_start:ikp_end, 2))
864
865 DO ikp = ikp_start, ikp_end
866 DO re_im = 1, 2
867 NULLIFY (mat_v_kp(ikp, re_im)%matrix)
868 ALLOCATE (mat_v_kp(ikp, re_im)%matrix)
869 CALL dbcsr_create(mat_v_kp(ikp, re_im)%matrix, template=bs_env%mat_RI_RI%matrix)
870 CALL dbcsr_reserve_all_blocks(mat_v_kp(ikp, re_im)%matrix)
871 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 atomic_kind_set=atomic_kind_set)
881
882 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 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
886
887 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 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_extra
892
893 ELSE
894
895 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 ikp_end=ikp_end)
910
911 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
912
913 ALLOCATE (fm_v_kp(ikp_start:ikp_end, 2))
914 DO ikp = ikp_start, ikp_end
915 DO re_im = 1, 2
916 CALL cp_fm_create(fm_v_kp(ikp, re_im), bs_env%fm_RI_RI%matrix_struct)
917 CALL copy_dbcsr_to_fm(mat_v_kp(ikp, re_im)%matrix, fm_v_kp(ikp, re_im))
918 CALL dbcsr_deallocate_matrix(mat_v_kp(ikp, re_im)%matrix)
919 END DO
920 END DO
921 DEALLOCATE (mat_v_kp)
922
923 CALL timestop(handle)
924
925 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 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_m_ikp
949
950 CALL timeset(routinen, handle)
951
952 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 regularization_ri=bs_env%regularization_RI, ikp_ext=ikp)
959
960 IF (ikp == 1) THEN
961 CALL cp_cfm_create(cfm_v_sqrt_ikp, fm_v_kp(ikp, 1)%matrix_struct)
962 CALL cp_cfm_create(cfm_m_inv_v_sqrt_ikp, fm_v_kp(ikp, 1)%matrix_struct)
963 END IF
964 CALL cp_cfm_create(cfm_m_inv_ikp, fm_v_kp(ikp, 1)%matrix_struct)
965
966 CALL cp_fm_to_cfm(fm_m_ikp(1, 1), fm_m_ikp(1, 2), cfm_m_inv_ikp)
967 CALL cp_fm_to_cfm(fm_v_kp(ikp, 1), fm_v_kp(ikp, 2), cfm_v_sqrt_ikp)
968
969 CALL cp_fm_release(fm_m_ikp)
970
971 CALL cp_cfm_create(cfm_work, fm_v_kp(ikp, 1)%matrix_struct)
972
973 ! M(k) -> M^-1(k)
974 CALL cp_cfm_to_cfm(cfm_m_inv_ikp, cfm_work)
975 CALL cp_cfm_cholesky_decompose(matrix=cfm_m_inv_ikp, n=n_ri, info_out=info)
976 IF (info == 0) THEN
977 ! successful Cholesky decomposition
978 CALL cp_cfm_cholesky_invert(cfm_m_inv_ikp)
979 ! symmetrize the result
980 CALL cp_cfm_upper_to_full(cfm_m_inv_ikp)
981 ELSE
982 ! Cholesky decomposition not successful: use expensive diagonalization
983 CALL cp_cfm_power(cfm_work, threshold=bs_env%eps_eigval_mat_RI, exponent=-1.0_dp)
984 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 CALL cp_cfm_to_cfm(cfm_v_sqrt_ikp, cfm_work)
989 CALL cp_cfm_cholesky_decompose(matrix=cfm_v_sqrt_ikp, n=n_ri, info_out=info)
990 IF (info == 0) THEN
991 ! successful Cholesky decomposition
992 CALL clean_lower_part(cfm_v_sqrt_ikp)
993 ELSE
994 ! Cholesky decomposition not successful: use expensive diagonalization
995 CALL cp_cfm_power(cfm_work, threshold=0.0_dp, exponent=0.5_dp)
996 CALL cp_cfm_to_cfm(cfm_work, cfm_v_sqrt_ikp)
997 END IF
998 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 z_zero, cfm_m_inv_v_sqrt_ikp)
1003
1004 CALL cp_cfm_release(cfm_m_inv_ikp)
1005
1006 CALL timestop(handle)
1007
1008 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 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 CALL timeset(routinen, handle)
1026
1027 CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
1028 CALL create_fm_w_mic_time(bs_env, fm_w_mic_time)
1029
1030 DO i_t = 1, bs_env%num_time_freq_points
1031
1032 bs_env%t1 = m_walltime()
1033
1034 CALL fm_read(fm_w_mic_time(i_t), bs_env, bs_env%W_time_name, i_t)
1035
1036 IF (bs_env%unit_nr > 0) THEN
1037 WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F7.1,A)') &
1038 τ'Read W^MIC(i) from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1039 ', Execution time', m_walltime() - bs_env%t1, ' s'
1040 END IF
1041
1042 END DO
1043
1044 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1045
1046 CALL timestop(handle)
1047
1048 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 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_v_kp
1069
1070 CALL timeset(routinen, handle)
1071
1072 CALL create_fm_w_mic_time(bs_env, fm_w_mic_time)
1073
1074 DO ikp_batch = 1, bs_env%num_chi_eps_W_batches
1075
1076 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 CALL compute_v_k_by_lattice_sum(bs_env, qs_env, fm_v_kp, ikp_batch)
1080
1081 DO ikp_in_batch = 1, bs_env%nkp_chi_eps_W_batch
1082
1083 ikp = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + ikp_in_batch
1084
1085 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 cfm_v_sqrt_ikp, cfm_m_inv_v_sqrt_ikp, ikp)
1089
1090 CALL bs_env%para_env%sync()
1091
1092 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 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 cfm_v_sqrt_ikp)
1101
1102 ! Fourier trafo from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
1103 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 CALL cp_fm_release(fm_v_kp(ikp, 1))
1108 CALL cp_fm_release(fm_v_kp(ikp, 2))
1109
1110 END DO ! ikp_in_batch
1111
1112 DEALLOCATE (fm_v_kp)
1113
1114 IF (bs_env%unit_nr > 0) THEN
1115 WRITE (bs_env%unit_nr, '(T2,A,I12,A,I3,A,F7.1,A)') &
1116 τ'Computed W(i,k) for k-point batch', &
1117 ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
1118 ', Execution time', m_walltime() - bs_env%t1, ' s'
1119 END IF
1120
1121 END DO ! ikp_batch
1122
1123 IF (bs_env%approx_kp_extrapol) THEN
1124 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 CALL multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_w_mic_time)
1129
1130 DO i_t = 1, bs_env%num_time_freq_points
1131 CALL fm_write(fm_w_mic_time(i_t), i_t, bs_env%W_time_name, qs_env)
1132 END DO
1133
1134 CALL cp_cfm_release(cfm_m_inv_v_sqrt_ikp)
1135 CALL cp_cfm_release(cfm_v_sqrt_ikp)
1136 CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
1137
1138 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1139
1140 CALL timestop(handle)
1141
1142 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 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 CALL timeset(routinen, handle)
1170
1171 ! 1. Fourier transformation of χ_PQ(iτ,k=0) to χ_PQ(iω_j,k=0)
1172 CALL compute_fm_chi_gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_gamma_tau)
1173
1174 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 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 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 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 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 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 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 "RI_AUX")
1206 ! non-kpoint extrapolated
1207 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 "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 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 wkp_ext=bs_env%wkp_orig)
1220 END IF
1221 END SELECT
1222
1223 CALL cp_cfm_release(cfm_w_ikp_freq_j)
1224
1225 CALL timestop(handle)
1226
1227 END SUBROUTINE compute_fm_w_mic_freq_j
1228
1229! **************************************************************************************************
1230!> \brief ...
1231!> \param cfm_mat ...
1232! **************************************************************************************************
1233 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1241
1242 CALL timeset(routinen, handle)
1243
1244 CALL cp_cfm_get_info(matrix=cfm_mat, &
1245 nrow_local=nrow_local, ncol_local=ncol_local, &
1246 row_indices=row_indices, col_indices=col_indices)
1247
1248 DO i_row = 1, nrow_local
1249 DO j_col = 1, ncol_local
1250 i_global = row_indices(i_row)
1251 j_global = col_indices(j_col)
1252 IF (j_global < i_global) cfm_mat%local_data(i_row, j_col) = z_zero
1253 END DO
1254 END DO
1255
1256 CALL timestop(handle)
1257
1258 END SUBROUTINE clean_lower_part
1259
1260! **************************************************************************************************
1261!> \brief ...
1262!> \param bs_env ...
1263!> \param fm_W_MIC_time ...
1264! **************************************************************************************************
1265 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 CALL timeset(routinen, handle)
1275
1276 CALL cp_fm_get_info(matrix=fm_w_mic_time(1), nrow_local=nrow_local, ncol_local=ncol_local)
1277
1278 DO i_t = 1, bs_env%num_time_freq_points
1279 DO i = 1, nrow_local
1280 DO j = 1, ncol_local
1281
1282 w_extra_1 = bs_env%fm_W_MIC_freq_1_extra%local_data(i, j)
1283 w_no_extra_1 = bs_env%fm_W_MIC_freq_1_no_extra%local_data(i, j)
1284
1285 IF (abs(w_no_extra_1) > 1.0e-13) THEN
1286 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 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 *extrapol_factor
1296 END DO
1297 END DO
1298 END DO
1299
1300 CALL timestop(handle)
1301
1302 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 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 CALL timeset(routinen, handle)
1323
1324 CALL dbcsr_set(bs_env%mat_RI_RI%matrix, 0.0_dp)
1325
1326 freq_j = bs_env%imag_freq_points(j_w)
1327
1328 DO i_t = 1, bs_env%num_time_freq_points
1329
1330 time_i = bs_env%imag_time_points(i_t)
1331 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 1.0_dp, cos(time_i*freq_j)*weight_ij)
1336
1337 END DO
1338
1339 CALL copy_dbcsr_to_fm(bs_env%mat_RI_RI%matrix, fm_chi_gamma_freq)
1340
1341 CALL timestop(handle)
1342
1343 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 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: block_im, block_re, data_block
1371 TYPE(cell_type), POINTER :: cell
1372 TYPE(dbcsr_iterator_type) :: iter
1373 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1374
1375 CALL timeset(routinen, handle)
1376
1377 ! get the same blocks in mat_ikp_re and mat_ikp_im as in mat_Gamma
1378 CALL dbcsr_copy(mat_ikp_re, mat_gamma)
1379 CALL dbcsr_copy(mat_ikp_im, mat_gamma)
1380 CALL dbcsr_set(mat_ikp_re, 0.0_dp)
1381 CALL dbcsr_set(mat_ikp_im, 0.0_dp)
1382
1383 NULLIFY (cell, particle_set)
1384 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1385 CALL get_cell(cell=cell, h=hmat)
1386
1387 index_to_cell => kpoints%index_to_cell
1388
1389 num_cells = SIZE(index_to_cell, 2)
1390
1391 DO i_cell = 1, num_cells
1392
1393 CALL dbcsr_iterator_start(iter, mat_gamma)
1394 DO WHILE (dbcsr_iterator_blocks_left(iter))
1395 CALL dbcsr_iterator_next_block(iter, row, col, data_block)
1396
1397 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 (pbc(particle_set(col)%r(1:3), cell) + cell_vector(1:3))
1401 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 i_cell_is_the_minimum_image_cell = .true.
1405 DO j_cell = 1, num_cells
1406 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 (pbc(particle_set(col)%r(1:3), cell) + cell_vector_j(1:3))
1409 abs_rab_cell_j = sqrt(rab_cell_j(1)**2 + rab_cell_j(2)**2 + rab_cell_j(3)**2)
1410
1411 IF (abs_rab_cell_i > abs_rab_cell_j + 1.0e-6_dp) THEN
1412 i_cell_is_the_minimum_image_cell = .false.
1413 END IF
1414 END DO
1415
1416 IF (i_cell_is_the_minimum_image_cell) THEN
1417 NULLIFY (block_re, block_im)
1418 CALL dbcsr_get_block_p(matrix=mat_ikp_re, row=row, col=col, block=block_re, found=f)
1419 CALL dbcsr_get_block_p(matrix=mat_ikp_im, row=row, col=col, block=block_im, found=f)
1420 cpassert(all(abs(block_re) < 1.0e-10_dp))
1421 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 REAL(index_to_cell(3, i_cell), dp)*kpoints%xkp(3, ikp)
1426
1427 block_re(:, :) = cos(twopi*arg)*data_block(:, :)
1428 block_im(:, :) = sin(twopi*arg)*data_block(:, :)
1429 END IF
1430
1431 END DO
1432 CALL dbcsr_iterator_stop(iter)
1433
1434 END DO
1435
1436 CALL timestop(handle)
1437
1438 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 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 CALL timeset(routinen, handle)
1461
1462 CALL cp_cfm_create(cfm_work, cfm_chi_ikp_freq_j%matrix_struct)
1463 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 cfm_chi_ikp_freq_j, cfm_m_inv_v_sqrt_ikp, z_zero, cfm_work)
1470 CALL cp_cfm_release(cfm_chi_ikp_freq_j)
1471
1472 ! 1. b) eps_work = V^0.5(k)*M^-1(k)*work
1473 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 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 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 CALL cp_cfm_cholesky_decompose(matrix=cfm_eps_ikp_freq_j, n=n_ri, info_out=info)
1484 cpassert(info == 0)
1485
1486 ! 2. b) Inversion of ε(iω_j,k) using its Cholesky decomposition
1487 CALL cp_cfm_cholesky_invert(cfm_eps_ikp_freq_j)
1488 CALL cp_cfm_upper_to_full(cfm_eps_ikp_freq_j)
1489
1490 ! 2. c) ε^-1(iω_j,k)-Id
1491 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 z_zero, cfm_work)
1496
1497 ! 2. e) W(iw,k) = V^0.5(k)*work
1498 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 z_zero, cfm_w_ikp_freq_j)
1501
1502 CALL cp_cfm_release(cfm_work)
1503 CALL cp_cfm_release(cfm_eps_ikp_freq_j)
1504
1505 CALL timestop(handle)
1506
1507 END SUBROUTINE compute_cfm_w_ikp_freq_j
1508
1509! **************************************************************************************************
1510!> \brief ...
1511!> \param cfm ...
1512!> \param alpha ...
1513! **************************************************************************************************
1514 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1524
1525 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 col_indices=col_indices)
1532
1533 ! add 1 on the diagonal
1534 DO j_col = 1, ncol_local
1535 j_global = col_indices(j_col)
1536 DO i_row = 1, nrow_local
1537 i_global = row_indices(i_row)
1538 IF (j_global == i_global) THEN
1539 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 CALL timestop(handle)
1545
1546 END SUBROUTINE cfm_add_on_diag
1547
1548! **************************************************************************************************
1549!> \brief ...
1550!> \param bs_env ...
1551!> \param fm_W_MIC_time ...
1552! **************************************************************************************************
1553 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 CALL timeset(routinen, handle)
1562
1563 ALLOCATE (fm_w_mic_time(bs_env%num_time_freq_points))
1564 DO i_t = 1, bs_env%num_time_freq_points
1565 CALL cp_fm_create(fm_w_mic_time(i_t), bs_env%fm_RI_RI%matrix_struct)
1566 END DO
1567
1568 CALL timestop(handle)
1569
1570 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 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 CALL timeset(routinen, handle)
1591
1592 freq_j = bs_env%imag_freq_points(j_w)
1593
1594 DO i_t = 1, bs_env%num_time_freq_points
1595
1596 time_i = bs_env%imag_time_points(i_t)
1597 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 beta=weight_ij*cos(time_i*freq_j), matrix_b=fm_w_mic_freq_j)
1602
1603 END DO
1604
1605 CALL timestop(handle)
1606
1607 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 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_minv_gamma
1625
1626 CALL timeset(routinen, handle)
1627
1628 n_ri = bs_env%n_RI
1629
1630 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 regularization_ri=bs_env%regularization_RI)
1636
1637 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 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 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 fm_minv_gamma(1, 1), 0.0_dp, fm_w_mic_time(i_t))
1647
1648 END DO
1649
1650 CALL cp_fm_release(fm_work)
1651 CALL cp_fm_release(fm_minv_gamma)
1652
1653 CALL timestop(handle)
1654
1655 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 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 CALL timeset(routinen, handle)
1673
1674 ALLOCATE (fm_sigma_x_gamma(bs_env%n_spin))
1675 DO ispin = 1, bs_env%n_spin
1676 CALL cp_fm_create(fm_sigma_x_gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1677 END DO
1678
1679 IF (bs_env%Sigma_x_exists) THEN
1680 DO ispin = 1, bs_env%n_spin
1681 CALL fm_read(fm_sigma_x_gamma(ispin), bs_env, bs_env%Sigma_x_name, ispin)
1682 END DO
1683 ELSE
1684 CALL compute_sigma_x(bs_env, qs_env, fm_sigma_x_gamma)
1685 END IF
1686
1687 CALL timestop(handle)
1688
1689 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 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_vtr_gamma
1707 TYPE(dbcsr_type) :: mat_sigma_x_gamma
1708 TYPE(dbt_type) :: t_2c_d, t_2c_sigma_x, t_2c_v, t_3c_x_v
1709
1710 CALL timeset(routinen, handle)
1711
1712 bs_env%t1 = m_walltime()
1713
1714 CALL dbt_create(bs_env%t_G, t_2c_d)
1715 CALL dbt_create(bs_env%t_W, t_2c_v)
1716 CALL dbt_create(bs_env%t_G, t_2c_sigma_x)
1717 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_v)
1718 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 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 CALL multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_vtr_gamma(:, 1))
1727
1728 DO ispin = 1, bs_env%n_spin
1729
1730 ! 3. Compute density matrix D_µν
1731 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 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 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 DO i_intval_idx = 1, bs_env%n_intervals_i
1743 DO j_intval_idx = 1, bs_env%n_intervals_j
1744 i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
1745 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 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 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 mat_sigma_x_gamma, bs_env%para_env)
1761
1762 CALL write_matrix(mat_sigma_x_gamma, ispin, bs_env%Sigma_x_name, &
1763 bs_env%fm_work_mo(1), qs_env)
1764
1765 CALL copy_dbcsr_to_fm(mat_sigma_x_gamma, fm_sigma_x_gamma(ispin))
1766
1767 END DO
1768
1769 IF (bs_env%unit_nr > 0) THEN
1770 WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
1771 Σ'Computed ^x(k=0),', ' Execution time', m_walltime() - bs_env%t1, ' s'
1772 WRITE (bs_env%unit_nr, '(A)') ' '
1773 END IF
1774
1775 CALL dbcsr_release(mat_sigma_x_gamma)
1776 CALL dbt_destroy(t_2c_d)
1777 CALL dbt_destroy(t_2c_v)
1778 CALL dbt_destroy(t_2c_sigma_x)
1779 CALL dbt_destroy(t_3c_x_v)
1780 CALL cp_fm_release(fm_vtr_gamma)
1781
1782 CALL timestop(handle)
1783
1784 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 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 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_sigma_neg_tau, mat_sigma_pos_tau
1806 TYPE(dbt_type) :: t_2c_gocc, t_2c_gvir, &
1807 t_2c_sigma_neg_tau, &
1808 t_2c_sigma_pos_tau, t_2c_w, t_3c_x_w
1809
1810 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 mat_sigma_neg_tau, mat_sigma_pos_tau)
1815
1816 DO i_t = 1, bs_env%num_time_freq_points
1817
1818 DO ispin = 1, bs_env%n_spin
1819
1820 bs_env%t1 = m_walltime()
1821
1822 read_write_index = i_t + (ispin - 1)*bs_env%num_time_freq_points
1823
1824 ! read self-energy from restart
1825 IF (bs_env%Sigma_c_exists(i_t, ispin)) THEN
1826 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 keep_sparsity=.false.)
1829 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 keep_sparsity=.false.)
1832 IF (bs_env%unit_nr > 0) THEN
1833 WRITE (bs_env%unit_nr, '(T2,2A,I3,A,I3,A,F7.1,A)') Στ'Read ^c(i,k=0) ', &
1834 'from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1835 ', Execution time', m_walltime() - bs_env%t1, ' s'
1836 END IF
1837
1838 cycle
1839
1840 END IF
1841
1842 tau = bs_env%imag_time_points(i_t)
1843
1844 CALL g_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.true., vir=.false.)
1845 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 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 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 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 DO i_intval_idx = 1, bs_env%n_intervals_i
1861 DO j_intval_idx = 1, bs_env%n_intervals_j
1862 i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
1863 j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
1864
1865 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 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 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 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 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 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 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 bs_env%Sigma_n_name, bs_env%fm_work_mo(1), qs_env)
1897
1898 IF (bs_env%unit_nr > 0) THEN
1899 WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F7.1,A)') &
1900 Στ'Computed ^c(i,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1901 ', Execution time', m_walltime() - bs_env%t1, ' s'
1902 END IF
1903
1904 END DO ! ispin
1905
1906 END DO ! i_t
1907
1908 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 mat_sigma_pos_tau, mat_sigma_neg_tau)
1912
1913 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 mat_sigma_neg_tau, mat_sigma_pos_tau)
1918
1919 CALL delete_unnecessary_files(bs_env)
1920
1921 CALL timestop(handle)
1922
1923 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 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 CALL timeset(routinen, handle)
1952
1953 CALL dbt_create(bs_env%t_G, t_2c_gocc)
1954 CALL dbt_create(bs_env%t_G, t_2c_gvir)
1955 CALL dbt_create(bs_env%t_W, t_2c_w)
1956 CALL dbt_create(bs_env%t_G, t_2c_sigma_neg_tau)
1957 CALL dbt_create(bs_env%t_G, t_2c_sigma_pos_tau)
1958 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_w)
1959
1960 NULLIFY (mat_sigma_neg_tau, mat_sigma_pos_tau)
1961 ALLOCATE (mat_sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1962 ALLOCATE (mat_sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1963
1964 DO i_t = 1, bs_env%num_time_freq_points
1965 DO ispin = 1, bs_env%n_spin
1966 ALLOCATE (mat_sigma_neg_tau(i_t, ispin)%matrix)
1967 ALLOCATE (mat_sigma_pos_tau(i_t, ispin)%matrix)
1968 CALL dbcsr_create(mat_sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1969 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 CALL timestop(handle)
1974
1975 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 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 TYPE(dbt_type) :: t_3c_for_w, t_3c_x_w_tmp
1998
1999 CALL timeset(routinen, handle)
2000
2001 CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_w_tmp)
2002 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 bs_env%i_RI_end_from_atom(j_atoms(2))]
2006
2007 DO ri_intval_idx = 1, bs_env%n_intervals_inner_loop_atoms
2008 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 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 filter_eps=bs_env%eps_filter)
2024
2025 END DO ! RI_atoms
2026
2027 ! 3. reorder tensor
2028 CALL dbt_copy(t_3c_x_w_tmp, t_3c_x_w, order=[1, 2, 3], move_data=.true.)
2029
2030 CALL dbt_destroy(t_3c_x_w_tmp)
2031 CALL dbt_destroy(t_3c_for_w)
2032
2033 CALL timestop(handle)
2034
2035 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 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 TYPE(dbt_type) :: t_3c_for_g, t_3c_x_g, t_3c_x_g_2
2067
2068 CALL timeset(routinen, handle)
2069
2070 cpassert(occ .EQV. (.NOT. vir))
2071 IF (occ) sign_sigma = -1.0_dp
2072 IF (vir) sign_sigma = 1.0_dp
2073
2074 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_g)
2075 CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_g)
2076 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 bs_env%i_ao_end_from_atom(i_atoms(2))]
2080
2081 DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
2082 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 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 filter_eps=bs_env%eps_filter)
2096
2097 END DO ! IL_atoms
2098
2099 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 filter_eps=bs_env%eps_filter, move_data=clear_w, flop=flop)
2109
2110 IF (PRESENT(can_skip)) THEN
2111 IF (flop == 0_int_8) can_skip = .true.
2112 END IF
2113
2114 CALL dbt_destroy(t_3c_for_g)
2115 CALL dbt_destroy(t_3c_x_g)
2116 CALL dbt_destroy(t_3c_x_g_2)
2117
2118 CALL timestop(handle)
2119
2120 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 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 CALL timeset(routinen, handle)
2141
2142 ALLOCATE (fm_sigma_c_gamma_time(bs_env%num_time_freq_points, 2, bs_env%n_spin))
2143 DO i_t = 1, bs_env%num_time_freq_points
2144 DO ispin = 1, bs_env%n_spin
2145 DO pos_neg = 1, 2
2146 CALL cp_fm_create(fm_sigma_c_gamma_time(i_t, pos_neg, ispin), &
2147 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 fm_sigma_c_gamma_time(i_t, 1, ispin))
2151 CALL copy_dbcsr_to_fm(mat_sigma_neg_tau(i_t, ispin)%matrix, &
2152 fm_sigma_c_gamma_time(i_t, 2, ispin))
2153 END DO
2154 END DO
2155
2156 CALL timestop(handle)
2157
2158 END SUBROUTINE fill_fm_sigma_c_gamma_time
2159
2160! **************************************************************************************************
2161!> \brief ...
2162!> \param bs_env ...
2163! **************************************************************************************************
2164 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 CALL timeset(routinen, handle)
2174
2175 n_skip = 0
2176
2177 DO i_intval_idx = 1, bs_env%n_intervals_i
2178 DO j_intval_idx = 1, bs_env%n_intervals_j
2179 IF (bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx) .AND. &
2180 bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) THEN
2181 n_skip = n_skip + 1
2182 END IF
2183 END DO
2184 END DO
2185
2186 IF (bs_env%unit_nr > 0) THEN
2187 WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
2188 Στ'Sparsity of ^c(i,k=0): Percentage of skipped atom pairs:', &
2189 REAL(100*n_skip, kind=dp)/real(i_intval_idx*j_intval_idx, kind=dp), ' %'
2190 END IF
2191
2192 CALL timestop(handle)
2193
2194 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 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 CALL timeset(routinen, handle)
2223
2224 CALL dbt_destroy(t_2c_gocc)
2225 CALL dbt_destroy(t_2c_gvir)
2226 CALL dbt_destroy(t_2c_w)
2227 CALL dbt_destroy(t_2c_sigma_neg_tau)
2228 CALL dbt_destroy(t_2c_sigma_pos_tau)
2229 CALL dbt_destroy(t_3c_x_w)
2230 CALL cp_fm_release(fm_w_mic_time)
2231 CALL dbcsr_deallocate_matrix_set(mat_sigma_neg_tau)
2232 CALL dbcsr_deallocate_matrix_set(mat_sigma_pos_tau)
2233
2234 CALL timestop(handle)
2235
2236 END SUBROUTINE destroy_mat_sigma_c
2237
2238! **************************************************************************************************
2239!> \brief ...
2240!> \param bs_env ...
2241! **************************************************************************************************
2242 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 CALL timeset(routinen, handle)
2251
2252 prefix = bs_env%prefix
2253
2254 DO i_t = 1, bs_env%num_time_freq_points
2255
2256 IF (i_t < 10) THEN
2257 WRITE (f_chi, '(3A,I1,A)') trim(prefix), bs_env%chi_name, "_00", i_t, ".matrix"
2258 WRITE (f_w_t, '(3A,I1,A)') trim(prefix), bs_env%W_time_name, "_00", i_t, ".matrix"
2259 ELSE IF (i_t < 100) THEN
2260 WRITE (f_chi, '(3A,I2,A)') trim(prefix), bs_env%chi_name, "_0", i_t, ".matrix"
2261 WRITE (f_w_t, '(3A,I2,A)') trim(prefix), bs_env%W_time_name, "_0", i_t, ".matrix"
2262 ELSE
2263 cpabort('Please implement more than 99 time/frequency points.')
2264 END IF
2265
2266 CALL safe_delete(f_chi, bs_env)
2267 CALL safe_delete(f_w_t, bs_env)
2268
2269 END DO
2270
2271 CALL timestop(handle)
2272
2273 END SUBROUTINE delete_unnecessary_files
2274
2275! **************************************************************************************************
2276!> \brief ...
2277!> \param filename ...
2278!> \param bs_env ...
2279! **************************************************************************************************
2280 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 CALL timeset(routinen, handle)
2290
2291 IF (bs_env%para_env%mepos == 0) THEN
2292
2293 INQUIRE (file=trim(filename), exist=file_exists)
2294 IF (file_exists) CALL mp_file_delete(trim(filename))
2295
2296 END IF
2297
2298 CALL timestop(handle)
2299
2300 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sigma_c_n_cos_time, sigma_c_n_sin_time
2319
2320 CALL timeset(routinen, handle)
2321
2322 ALLOCATE (sigma_c_n_cos_time(bs_env%n_ao, bs_env%num_time_freq_points))
2323 ALLOCATE (sigma_c_n_sin_time(bs_env%n_ao, bs_env%num_time_freq_points))
2324
2325 sigma_c_n_cos_time(:, :) = 0.5_dp*(sigma_c_n_time(:, :, 1) + sigma_c_n_time(:, :, 2))
2326 sigma_c_n_sin_time(:, :) = 0.5_dp*(sigma_c_n_time(:, :, 1) - sigma_c_n_time(:, :, 2))
2327
2328 sigma_c_n_freq(:, :, :) = 0.0_dp
2329
2330 DO i_t = 1, bs_env%num_time_freq_points
2331
2332 DO j_w = 1, bs_env%num_time_freq_points
2333
2334 freq_j = bs_env%imag_freq_points(j_w)
2335 time_i = bs_env%imag_time_points(i_t)
2336 ! integration weights for cosine and sine transform
2337 w_cos_ij = bs_env%weights_cos_t_to_w(j_w, i_t)*cos(freq_j*time_i)
2338 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 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 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 n_occ = bs_env%n_occ(ispin)
2356 sigma_c_n_freq(1:n_occ, :, 2) = -sigma_c_n_freq(1:n_occ, :, 2)
2357
2358 CALL timestop(handle)
2359
2360 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 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 CALL timeset(routinen, handle)
2385
2386 CALL cp_cfm_create(cfm_mos_ikp, bs_env%fm_s_Gamma%matrix_struct)
2387 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 ALLOCATE (v_xc_ikp_n(bs_env%n_ao), sigma_x_ikp_n(bs_env%n_ao))
2390 ALLOCATE (sigma_c_ikp_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
2391 ALLOCATE (sigma_c_ikp_n_freq(bs_env%n_ao, bs_env%num_time_freq_points, 2))
2392
2393 DO ispin = 1, bs_env%n_spin
2394
2395 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 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 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 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 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 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 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 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 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 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 ikp, ispin)
2434
2435 END DO ! ikp_DOS
2436
2437 END DO ! ispin
2438
2439 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 CALL g0w0_hamiltonian(bs_env)
2443
2444 CALL cp_fm_release(fm_sigma_x_gamma)
2445 CALL cp_fm_release(fm_sigma_c_gamma_time)
2446 CALL cp_cfm_release(cfm_ks_ikp)
2447 CALL cp_cfm_release(cfm_s_ikp)
2448 CALL cp_cfm_release(cfm_mos_ikp)
2449 CALL cp_cfm_release(cfm_work_ikp)
2450 CALL cp_cfm_release(cfm_sigma_x_ikp)
2451
2452 CALL timestop(handle)
2453
2454 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 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 CALL timeset(routinen, handle)
2480
2481 CALL cp_fm_create(fm_ikp_mo_re, fm_gamma%matrix_struct)
2482
2483 CALL fm_gamma_ao_to_cfm_ikp_mo(fm_gamma, fm_ikp_mo_re, ikp, qs_env, bs_env, cfm_mos_ikp)
2484
2485 CALL cp_fm_get_diag(fm_ikp_mo_re, array_ikp_n)
2486
2487 CALL cp_fm_release(fm_ikp_mo_re)
2488
2489 CALL timestop(handle)
2490
2491 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 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 CALL timeset(routinen, handle)
2515
2516 CALL cp_cfm_create(cfm_ikp_ao, fm_gamma%matrix_struct)
2517 CALL cp_cfm_create(cfm_ikp_mo, fm_gamma%matrix_struct)
2518 CALL cp_cfm_create(cfm_tmp, fm_gamma%matrix_struct)
2519
2520 ! get cfm_µν(k_i) from fm_µν(k=0)
2521 CALL cfm_ikp_from_fm_gamma(cfm_ikp_ao, fm_gamma, ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2522
2523 nmo = bs_env%n_ao
2524 CALL parallel_gemm('N', 'N', nmo, nmo, nmo, z_one, cfm_ikp_ao, cfm_mos_ikp, z_zero, cfm_tmp)
2525 CALL parallel_gemm('C', 'N', nmo, nmo, nmo, z_one, cfm_mos_ikp, cfm_tmp, z_zero, cfm_ikp_mo)
2526
2527 CALL cp_cfm_to_fm(cfm_ikp_mo, fm_ikp_mo_re)
2528
2529 CALL cp_cfm_release(cfm_ikp_mo)
2530 CALL cp_cfm_release(cfm_ikp_ao)
2531 CALL cp_cfm_release(cfm_tmp)
2532
2533 CALL timestop(handle)
2534
2535 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 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 CALL timeset(routinen, handle)
2561
2562 n_mo = bs_env%n_ao
2563 ALLOCATE (dummy(n_mo), sigma_c_ikp_n_qp(n_mo))
2564 sigma_c_ikp_n_qp(:) = 0.0_dp
2565
2566 DO i_mo = 1, n_mo
2567
2568 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 0.0_dp, .true., .false., 1)
2581
2582 END DO
2583
2584 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 v_xc_ikp_n(:)
2590
2591 CALL get_fname(fname, bs_env, ikp, "SCF_and_G0W0", ispin=ispin)
2592
2593 IF (bs_env%para_env%is_source()) THEN
2594
2595 CALL open_file(trim(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE")
2596
2597 WRITE (iunit, "(A)") " "
2598 WRITE (iunit, "(A10,3F10.4)") "kpoint: ", bs_env%kpoints_DOS%xkp(:, ikp)
2599 WRITE (iunit, "(A)") " "
2600 WRITE (iunit, "(A5,A24,2A17,A16,A18)") "n", ϵ"_nk^DFT (eV)", Σ"^c_nk (eV)", &
2601 Σ"^x_nk (eV)", "v_n^xc (eV)", ϵ"_nk^G0W0 (eV)"
2602 WRITE (iunit, "(A)") " "
2603
2604 DO i_mo = 1, n_mo
2605 IF (i_mo .LE. bs_env%n_occ(ispin)) occ_vir = 'occ'
2606 IF (i_mo > bs_env%n_occ(ispin)) occ_vir = 'vir'
2607 WRITE (iunit, "(I5,3A,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', &
2608 bs_env%eigenval_scf(i_mo, ikp, ispin)*evolt, &
2609 sigma_c_ikp_n_qp(i_mo)*evolt, &
2610 sigma_x_ikp_n(i_mo)*evolt, &
2611 v_xc_ikp_n(i_mo)*evolt, &
2612 bs_env%eigenval_G0W0(i_mo, ikp, ispin)*evolt
2613 END DO
2614
2615 CALL close_file(iunit)
2616
2617 END IF
2618
2619 CALL timestop(handle)
2620
2621 END SUBROUTINE analyt_conti_and_print
2622
2623! **************************************************************************************************
2624!> \brief ...
2625!> \param bs_env ...
2626! **************************************************************************************************
2627 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 CALL timeset(routinen, handle)
2638
2639 n_mo = bs_env%n_ao
2640
2641 bs_env%band_edges_scf%DBG = 1000.0_dp
2642 bs_env%band_edges_G0W0%DBG = 1000.0_dp
2643
2644 SELECT CASE (bs_env%n_spin)
2645 CASE (1)
2646 homo = bs_env%n_occ(1)
2647 lumo = homo + 1
2648 bs_env%band_edges_scf%VBM = maxval(bs_env%eigenval_scf(1:homo, :, 1))
2649 bs_env%band_edges_scf%CBM = minval(bs_env%eigenval_scf(homo + 1:n_mo, :, 1))
2650 bs_env%band_edges_G0W0%VBM = maxval(bs_env%eigenval_G0W0(1:homo, :, 1))
2651 bs_env%band_edges_G0W0%CBM = minval(bs_env%eigenval_G0W0(homo + 1:n_mo, :, 1))
2652 CASE (2)
2653 homo_1 = bs_env%n_occ(1)
2654 lumo_1 = homo_1 + 1
2655 homo_2 = bs_env%n_occ(2)
2656 lumo_2 = homo_2 + 1
2657 bs_env%band_edges_scf%VBM = max(maxval(bs_env%eigenval_scf(1:homo_1, :, 1)), &
2658 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 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 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 minval(bs_env%eigenval_G0W0(homo_2 + 1:n_mo, :, 2)))
2665 CASE DEFAULT
2666 cpabort("Error with number of spins.")
2667 END SELECT
2668
2669 bs_env%band_edges_scf%IDBG = bs_env%band_edges_scf%CBM - bs_env%band_edges_scf%VBM
2670 bs_env%band_edges_G0W0%IDBG = bs_env%band_edges_G0W0%CBM - bs_env%band_edges_G0W0%VBM
2671
2672 DO ispin = 1, bs_env%n_spin
2673
2674 homo = bs_env%n_occ(ispin)
2675
2676 DO ikp = 1, bs_env%kpoints_DOS%nkp
2677 e_dbg_scf_at_ikp = -maxval(bs_env%eigenval_scf(1:homo, ikp, ispin)) + &
2678 minval(bs_env%eigenval_scf(homo + 1:n_mo, ikp, ispin))
2679 IF (e_dbg_scf_at_ikp < bs_env%band_edges_scf%DBG) THEN
2680 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 minval(bs_env%eigenval_G0W0(homo + 1:n_mo, ikp, ispin))
2685 IF (e_dbg_g0w0_at_ikp < bs_env%band_edges_G0W0%DBG) THEN
2686 bs_env%band_edges_G0W0%DBG = e_dbg_g0w0_at_ikp
2687 END IF
2688 END DO
2689 END DO
2690
2691 CALL timestop(handle)
2692
2693 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices
2707 REAL(kind=dp) :: e_j
2708
2709 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 col_indices=col_indices)
2717
2718 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 DO i_row_local = 1, nrow_local
2722 DO j_col_local = 1, ncol_local
2723
2724 j_col = col_indices(j_col_local)
2725
2726 ! the last k-point of eigenvalues_G0W0 is the Γ-point
2727 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 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 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 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 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 beta=0.0_dp, matrix_c=bs_env%fm_h_G0W0_Gamma)
2753
2754 CALL timestop(handle)
2755
2756 END SUBROUTINE g0w0_hamiltonian
2757
2758END MODULE gw_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:117
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:153
tensor
struct tensor_ tensor
cell_types::pbc
Definition cell_types.F:92
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:55
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:89
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
constants_operator
constants for the different operators of the 2c-integrals
Definition constants_operator.F:14
constants_operator::operator_coulomb
integer, parameter, public operator_coulomb
Definition constants_operator.F:22
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_cholesky_invert
subroutine, public cp_cfm_cholesky_invert(matrix, n, info_out)
Used to replace Cholesky decomposition by the inverse.
Definition cp_cfm_basic_linalg.F:982
cp_cfm_basic_linalg::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_basic_linalg.F:926
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:145
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:159
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:817
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:607
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:765
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:208
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:118
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:494
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:167
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:896
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:570
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1016
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition cp_fm_types.F:2147
cp_fm_types::cp_fm_read_unformatted
subroutine, public cp_fm_read_unformatted(fm, unit)
...
Definition cp_fm_types.F:2379
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
gw_communication
Definition gw_communication.F:13
gw_communication::global_matrix_to_local_matrix
subroutine, public global_matrix_to_local_matrix(mat_global, mat_local, para_env, num_pe_sub, atom_ranges)
...
Definition gw_communication.F:52
gw_communication::local_matrix_to_global_matrix
subroutine, public local_matrix_to_global_matrix(mat_local, mat_global, para_env)
...
Definition gw_communication.F:415
gw_methods
Definition gw_methods.F:13
gw_methods::gw
subroutine, public gw(qs_env, bs_env, post_scf_bandstructure_section)
Perform GW band structure calculation.
Definition gw_methods.F:111
gw_utils
Definition gw_utils.F:13
gw_utils::de_init_bs_env
subroutine, public de_init_bs_env(bs_env)
...
Definition gw_utils.F:173
gw_utils::create_and_init_bs_env_for_gw
subroutine, public create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec)
...
Definition gw_utils.F:115
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::ri_rpa_g0w0_crossing_newton
integer, parameter, public ri_rpa_g0w0_crossing_newton
Definition input_constants.F:1063
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_coulomb_2c
Routines to compute the Coulomb integral V_(alpha beta)(k) for a k-point k using lattice summation in...
Definition kpoint_coulomb_2c.F:15
kpoint_coulomb_2c::build_2c_coulomb_matrix_kp
subroutine, public build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, atomic_kind_set, size_lattice_sum, operator_type, ikp_start, ikp_end)
...
Definition kpoint_coulomb_2c.F:69
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:123
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_file_delete
subroutine, public mp_file_delete(filepath, info)
Deletes a file. Auxiliary routine to emulate 'replace' action for mp_file_open. Only the master proce...
Definition message_passing.F:3947
mp2_ri_2c
Framework for 2c-integrals for RI.
Definition mp2_ri_2c.F:14
mp2_ri_2c::ri_2c_integral_mat
subroutine, public ri_2c_integral_mat(qs_env, fm_matrix_minv_l_kpoints, fm_matrix_l, dimen_ri, ri_metric, do_kpoints, kpoints, put_mat_ks_env, regularization_ri, ikp_ext)
...
Definition mp2_ri_2c.F:554
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
post_scf_bandstructure_types
Definition post_scf_bandstructure_types.F:13
post_scf_bandstructure_utils
Definition post_scf_bandstructure_utils.F:13
post_scf_bandstructure_utils::cfm_ikp_from_fm_gamma
subroutine, public cfm_ikp_from_fm_gamma(cfm_ikp, fm_gamma, ikp, qs_env, kpoints, basis_type)
...
Definition post_scf_bandstructure_utils.F:2701
post_scf_bandstructure_utils::mic_contribution_from_ikp
subroutine, public mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, cfm_w_ikp_freq_j, ikp, kpoints, basis_type, wkp_ext)
...
Definition post_scf_bandstructure_utils.F:2833
post_scf_bandstructure_utils::get_fname
subroutine, public get_fname(fname, bs_env, ikp, scf_gw, ispin, soc)
...
Definition post_scf_bandstructure_utils.F:2320
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors.F:11
qs_tensors::build_3c_integrals
subroutine, public build_3c_integrals(t3c, filter_eps, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, int_eps, op_pos, do_kpoints, do_hfx_kpoints, desymmetrize, bounds_i, bounds_j, bounds_k, ri_range, img_to_ri_cell)
Build 3-center integral tensor.
Definition qs_tensors.F:2376
rpa_gw_kpoints_util
Routines treating GW and RPA calculations with kpoints.
Definition rpa_gw_kpoints_util.F:13
rpa_gw_kpoints_util::cp_cfm_upper_to_full
subroutine, public cp_cfm_upper_to_full(cfm_mat_q)
...
Definition rpa_gw_kpoints_util.F:1430
rpa_gw_kpoints_util::cp_cfm_power
subroutine, public cp_cfm_power(matrix, threshold, exponent, min_eigval)
...
Definition rpa_gw_kpoints_util.F:277
rpa_gw
Routines for GW, continuous development [Jan Wilhelm].
Definition rpa_gw.F:14
rpa_gw::continuation_pade
subroutine, public continuation_pade(vec_gw_energ, vec_omega_fit_gw, z_value, m_value, vec_sigma_c_gw, vec_sigma_x_minus_vxc_gw, eigenval, eigenval_scf, n_level_gw, gw_corr_lev_occ, nparam_pade, num_fit_points, crossing_search, homo, fermi_level_offset, do_gw_im_time, print_self_energy, count_ev_sc_gw, vec_gw_dos, dos_lower_bound, dos_precision, ndos, min_level_self_energy, max_level_self_energy, dos_eta, dos_min, dos_max)
perform analytic continuation with pade approximation
Definition rpa_gw.F:4299
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:138
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
post_scf_bandstructure_types::post_scf_bandstructure_type
Definition post_scf_bandstructure_types.F:52
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171