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