(git:6344619)
Loading...
Searching...
No Matches
gw_compute_V_and_W.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Coulomb matrices and screened interactions for imaginary-time GW.
10! **************************************************************************************************
11MODULE gw_compute_v_and_w
12 USE atomic_kind_types, ONLY: atomic_kind_type
13 USE cell_types, ONLY: cell_type
14 USE constants_operator, ONLY: operator_coulomb
15 USE cp_blacs_env, ONLY: cp_blacs_env_type
16 USE cp_cfm_basic_linalg, ONLY: cp_cfm_add_on_diag,&
17 cp_cfm_uplo_to_full
18 USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
19 cp_cfm_cholesky_invert
20 USE cp_cfm_diag, ONLY: cp_cfm_power
21 USE cp_cfm_types, ONLY: cp_cfm_create,&
22 cp_cfm_get_info,&
23 cp_cfm_release,&
24 cp_cfm_to_cfm,&
25 cp_cfm_type
26 USE cp_dbcsr_api, ONLY: &
27 dbcsr_add, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_release, &
28 dbcsr_distribution_type, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, &
29 dbcsr_type_no_symmetry
30 USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
31 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 cp_dbcsr_dist2d_to_dist,&
33 dbcsr_deallocate_matrix_set
34 USE cp_files, ONLY: close_file,&
35 open_file
36 USE cp_fm_basic_linalg, ONLY: cp_fm_add_on_diag,&
37 cp_fm_scale_and_add,&
38 cp_fm_uplo_to_full
39 USE cp_fm_redistribution, ONLY: cp_fm_to_replicated_array
40 USE cp_fm_types, ONLY: cp_fm_create,&
41 cp_fm_get_info,&
42 cp_fm_read_unformatted,&
43 cp_fm_release,&
44 cp_fm_set_all,&
45 cp_fm_to_fm,&
46 cp_fm_type,&
47 cp_fm_write_unformatted
48 USE cp_log_handling, ONLY: cp_get_default_logger,&
49 cp_logger_type
50 USE cp_output_handling, ONLY: cp_p_file,&
51 cp_print_key_should_output,&
52 cp_print_key_unit_nr
53 USE distribution_2d_types, ONLY: distribution_2d_type
54 USE gw_utils_fm, ONLY: cfm_contract_aba,&
55 fm_contract_aba,&
56 fm_invert,&
57 fm_sqrt
58 USE input_constants, ONLY: rtp_method_bse,&
59 rtp_method_bse_linearized
60 USE input_section_types, ONLY: section_vals_type
61 USE kinds, ONLY: default_path_length,&
62 dp
63 USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
64 USE libint_2c_3c, ONLY: libint_potential_type
65 USE machine, ONLY: m_walltime
66 USE mathconstants, ONLY: z_one,&
67 z_zero
68 USE mp2_ri_2c, ONLY: ri_2c_integral_mat_kp
69 USE parallel_gemm_api, ONLY: parallel_gemm
70 USE particle_methods, ONLY: get_particle_set
71 USE particle_types, ONLY: particle_type
72 USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
73 USE post_scf_bandstructure_utils, ONLY: mic_contribution_from_ikp,&
74 cfm_ikp_from_fm_gamma
75 USE qs_environment_types, ONLY: get_qs_env,&
76 qs_environment_type
77 USE qs_kind_types, ONLY: qs_kind_type
78 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
79 release_neighbor_list_sets
80 USE qs_tensors, ONLY: build_2c_integrals,&
81 build_2c_neighbor_lists
82 USE rpa_gw_kpoints_util, ONLY: cfm_screened_interaction_from_eps_cholesky
83#include "./base/base_uses.f90"
84
85 IMPLICIT NONE
86 PRIVATE
87 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_compute_V_and_W'
88 PUBLIC :: compute_w, &
89 get_w_mic, &
90 compute_v_k_by_lattice_sum, &
91 compute_minvvsqrt_vsqrt, &
92 compute_fm_w_mic_freq_j_from_chi, &
93 apply_extrapol_factor, &
94 fm_read, &
95 fm_write, &
96 get_v_tr_r
97
98CONTAINS
99
100! **************************************************************************************************
101!> \brief Computes the screened Coulomb interaction on the imaginary-time grid, entirely in the
102!> RI auxiliary (PQ) basis:
103!> χ_PQ(iω) = Σ_τ w(ω,τ) cos(ωτ) χ_PQ(iτ) (cosine transform)
104!> ε(iω) = Id - V^0.5 M^-1 χ(iω) M^-1 V^0.5 (dielectric function)
105!> W(iω) = V^0.5 ( ε^-1(iω) - Id ) V^0.5 (correlation part only)
106!> W(iτ) = Σ_ω w̃(τ,ω) cos(ωτ) W(iω) (back transform)
107!> W(iτ) <- M^-1 W(iτ) M^-1 (fold in the RI metric)
108!> where V is the bare Coulomb matrix and M the RI metric.
109!> \param bs_env ...
110!> \param qs_env ...
111!> \param mat_chi_Gamma_tau ...
112!> \param fm_W_time ...
113! **************************************************************************************************
114 SUBROUTINE compute_w(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_time)
115 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
116 TYPE(qs_environment_type), POINTER :: qs_env
117 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
118 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_time
119
120 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_W'
121
122 INTEGER :: handle, i_t, j_w
123 REAL(kind=dp) :: t1
124 TYPE(cp_fm_type) :: fm_m_inv_v_sqrt, fm_v, fm_v_sqrt
125
126 CALL timeset(routinen, handle)
127
128 t1 = m_walltime()
129
130 CALL create_fm_w_mic_time(bs_env, fm_w_time)
131
132 ! 1. Allocate V and M matrices
133 CALL cp_fm_create(fm_v, bs_env%fm_RI_RI%matrix_struct)
134 CALL cp_fm_create(fm_v_sqrt, bs_env%fm_RI_RI%matrix_struct)
135 CALL cp_fm_create(fm_m_inv_v_sqrt, bs_env%fm_RI_RI%matrix_struct)
136
137 ! Compute V and M^-1 * V^0.5
138 CALL compute_v_minvvsqrt(bs_env, qs_env, fm_v, fm_v_sqrt, fm_m_inv_v_sqrt)
139
140 ! 2. Loop over frequencies
141 DO j_w = 1, bs_env%num_time_freq_points
142 ! Fourier transformation of χ_PQ(iτ) to χ_PQ(iω_j)
143 CALL compute_fm_chi_gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_gamma_tau)
144
145 ! ε(iω_j) = Id - V^0.5*M^-1*χ(iω_j)*M^-1*V^0.5
146 ! W(iω_j) = V^0.5*(ε^-1(iω_j)-Id)*V^0.5
147 CALL compute_fm_w_freq(bs_env, bs_env%fm_chi_Gamma_freq, fm_v_sqrt, &
148 fm_m_inv_v_sqrt, bs_env%fm_W_MIC_freq)
149
150 ! Fourier transform from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
151 CALL fourier_transform_w_to_t(bs_env, fm_w_time, bs_env%fm_W_MIC_freq, j_w)
152 END DO
153
154 ! M^-1(k=0) W^MIC(iτ) M^-1(k=0) -> fm_W_time
155 CALL fm_contract_aba(bs_env%fm_Minv_Gamma, fm_w_time)
156
157 IF (bs_env%unit_nr > 0) THEN
158 WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
159 'Computed W(iτ),', ' Execution time', m_walltime() - t1, ' s'
160 END IF
161
162 CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
163
164 ! Cleanup
165 CALL cp_fm_release(fm_v)
166 CALL cp_fm_release(fm_v_sqrt)
167 CALL cp_fm_release(fm_m_inv_v_sqrt)
168
169 ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0
170 IF (bs_env%rtp_method == rtp_method_bse) THEN
171 t1 = m_walltime()
172 CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
173 ! Set to zero
174 CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
175 ! Sum over all times
176 DO i_t = 1, bs_env%num_time_freq_points
177 ! Add the relevant structure with correct weight
178 CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
179 bs_env%time_frequency_grid%time_weights_at_zero_frequency(i_t), fm_w_time(i_t))
180 END DO
181 ! Done, save to file
182 CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
183 ! Report calculation
184 IF (bs_env%unit_nr > 0) THEN
185 WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
186 'Computed W(0),', ' Execution time', m_walltime() - t1, ' s'
187 END IF
188 END IF
189
190 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
191
192 CALL timestop(handle)
193
194 END SUBROUTINE compute_w
195
196! **************************************************************************************************
197!> \brief ...
198!> \param bs_env ...
199!> \param qs_env ...
200!> \param mat_chi_Gamma_tau ...
201!> \param fm_W_MIC_time ...
202! **************************************************************************************************
203 SUBROUTINE get_w_mic(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
204 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
205 TYPE(qs_environment_type), POINTER :: qs_env
206 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
207 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
208
209 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_W_MIC'
210
211 INTEGER :: handle
212
213 CALL timeset(routinen, handle)
214
215 IF (bs_env%all_W_exist) THEN
216 CALL read_w_mic_time(bs_env, mat_chi_gamma_tau, fm_w_mic_time)
217 ELSE
218 CALL compute_w_mic(bs_env, qs_env, mat_chi_gamma_tau, fm_w_mic_time)
219 END IF
220
221 CALL timestop(handle)
222
223 END SUBROUTINE get_w_mic
224
225! **************************************************************************************************
226!> \brief ...
227!> \param bs_env ...
228!> \param qs_env ...
229!> \param mat_chi_Gamma_tau ...
230!> \param fm_W_MIC_time ...
231! **************************************************************************************************
232 SUBROUTINE compute_w_mic(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
233 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
234 TYPE(qs_environment_type), POINTER :: qs_env
235 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
236 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
237
238 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_W_MIC'
239
240 INTEGER :: handle, i_t, ikp, ikp_batch, &
241 ikp_in_batch, j_w
242 REAL(kind=dp) :: t1
243 TYPE(cp_cfm_type) :: cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp
244 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_v_kp
245
246 CALL timeset(routinen, handle)
247
248 CALL create_fm_w_mic_time(bs_env, fm_w_mic_time)
249
250 DO ikp_batch = 1, bs_env%num_chi_eps_W_batches
251
252 t1 = m_walltime()
253
254 ! Compute V_PQ(k) = sum_R e^(ikR) <phi_P, cell 0 | 1/r | phi_Q, cell R>
255 CALL compute_v_k_by_lattice_sum(bs_env, qs_env, cfm_v_kp, ikp_batch)
256
257 DO ikp_in_batch = 1, bs_env%nkp_chi_eps_W_batch
258
259 ikp = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + ikp_in_batch
260
261 IF (ikp > bs_env%nkp_chi_eps_W_orig_plus_extra) THEN
262 cycle
263 END IF
264
265 CALL compute_minvvsqrt_vsqrt(bs_env, qs_env, cfm_v_kp, &
266 cfm_v_sqrt_ikp, cfm_m_inv_v_sqrt_ikp, ikp)
267
268 CALL bs_env%para_env%sync()
269 CALL cp_cfm_release(cfm_v_kp(ikp))
270
271 DO j_w = 1, bs_env%num_time_freq_points
272
273 ! check if we need this (ikp, ω_j) combination for approximate k-point extrapolation
274 IF (bs_env%approx_kp_extrapol .AND. j_w > 1 .AND. &
275 ikp > bs_env%nkp_chi_eps_W_orig) cycle
276
277 CALL compute_fm_w_mic_freq_j(bs_env, qs_env, bs_env%fm_W_MIC_freq, j_w, ikp, &
278 mat_chi_gamma_tau, cfm_m_inv_v_sqrt_ikp, &
279 cfm_v_sqrt_ikp)
280
281 ! Fourier trafo from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
282 CALL fourier_transform_w_to_t(bs_env, fm_w_mic_time, bs_env%fm_W_MIC_freq, j_w)
283
284 END DO ! ω_j
285
286 END DO ! ikp_in_batch
287
288 DEALLOCATE (cfm_v_kp)
289
290 IF (bs_env%unit_nr > 0) THEN
291 WRITE (bs_env%unit_nr, '(T2,A,I12,A,I3,A,F10.1,A)') &
292 'Computed W(iτ,k) for k-point batch', &
293 ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
294 ', Execution time', m_walltime() - t1, ' s'
295 END IF
296
297 END DO ! ikp_batch
298
299 IF (bs_env%approx_kp_extrapol) THEN
300 CALL apply_extrapol_factor(bs_env, fm_w_mic_time)
301 END IF
302
303 ! M^-1(k=0) W^MIC(iτ) M^-1(k=0) -> fm_W_MIC_time
304 CALL fm_contract_aba(bs_env%fm_Minv_Gamma, fm_w_mic_time)
305
306 DO i_t = 1, bs_env%num_time_freq_points
307 CALL fm_write(fm_w_mic_time(i_t), i_t, bs_env%W_time_name, qs_env)
308 END DO
309
310 CALL cp_cfm_release(cfm_m_inv_v_sqrt_ikp)
311 CALL cp_cfm_release(cfm_v_sqrt_ikp)
312 CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
313
314 ! W^MIC(iω=0) from W^MIC(iτ) for the post-GW consumers; the file is the RT-BSE restart's
315 IF (bs_env%rtp_method == rtp_method_bse .OR. &
316 bs_env%rtp_method == rtp_method_bse_linearized) THEN
317 CALL sum_w_mic_freq_zero(bs_env, fm_w_mic_time)
318 CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
319 ELSE IF (bs_env%keep_W_freq_zero) THEN
320 CALL sum_w_mic_freq_zero(bs_env, fm_w_mic_time)
321 END IF
322
323 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
324
325 CALL timestop(handle)
326
327 END SUBROUTINE compute_w_mic
328
329! **************************************************************************************************
330!> \brief ...
331!> \param bs_env ...
332!> \param qs_env ...
333!> \param cfm_V_kp ...
334!> \param ikp_batch ...
335! **************************************************************************************************
336 SUBROUTINE compute_v_k_by_lattice_sum(bs_env, qs_env, cfm_V_kp, ikp_batch)
337 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
338 TYPE(qs_environment_type), POINTER :: qs_env
339 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_v_kp
340 INTEGER :: ikp_batch
341
342 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_V_k_by_lattice_sum'
343
344 INTEGER :: handle, ikp, ikp_end, ikp_start, &
345 nkp_chi_eps_w_batch, re_im
346 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
347 TYPE(cell_type), POINTER :: cell
348 TYPE(cp_fm_type) :: fm_v_kp_work
349 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_v_kp
350 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
351 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
352
353 CALL timeset(routinen, handle)
354
355 nkp_chi_eps_w_batch = bs_env%nkp_chi_eps_W_batch
356
357 ikp_start = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + 1
358 ikp_end = min(ikp_batch*bs_env%nkp_chi_eps_W_batch, bs_env%kpoints_chi_eps_W%nkp)
359
360 NULLIFY (mat_v_kp)
361 ALLOCATE (mat_v_kp(ikp_start:ikp_end, 2))
362
363 DO re_im = 1, 2
364 DO ikp = ikp_start, ikp_end
365 NULLIFY (mat_v_kp(ikp, re_im)%matrix)
366 ALLOCATE (mat_v_kp(ikp, re_im)%matrix)
367 CALL dbcsr_create(mat_v_kp(ikp, re_im)%matrix, template=bs_env%mat_RI_RI%matrix)
368 CALL dbcsr_reserve_all_blocks(mat_v_kp(ikp, re_im)%matrix)
369 CALL dbcsr_set(mat_v_kp(ikp, re_im)%matrix, 0.0_dp)
370 END DO ! ikp
371 END DO ! re_im
372
373 CALL get_qs_env(qs_env=qs_env, &
374 particle_set=particle_set, &
375 cell=cell, &
376 qs_kind_set=qs_kind_set, &
377 atomic_kind_set=atomic_kind_set)
378
379 IF (ikp_end <= bs_env%nkp_chi_eps_W_orig) THEN
380
381 ! 1. 2c Coulomb integrals for the first "original" k-point grid
382 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
383
384 ELSE IF (ikp_start > bs_env%nkp_chi_eps_W_orig .AND. &
385 ikp_end <= bs_env%nkp_chi_eps_W_orig_plus_extra) THEN
386
387 ! 2. 2c Coulomb integrals for the second "extrapolation" k-point grid
388 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_extra
389
390 ELSE
391
392 cpabort("Error with k-point parallelization.")
393
394 END IF
395
396 CALL build_2c_coulomb_matrix_kp(mat_v_kp, &
397 bs_env%kpoints_chi_eps_W, &
398 basis_type="RI_AUX", &
399 cell=cell, &
400 particle_set=particle_set, &
401 qs_kind_set=qs_kind_set, &
402 atomic_kind_set=atomic_kind_set, &
403 size_lattice_sum=bs_env%size_lattice_sum_V, &
404 operator_type=operator_coulomb, &
405 ikp_start=ikp_start, &
406 ikp_end=ikp_end)
407
408 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
409
410 ALLOCATE (cfm_v_kp(ikp_start:ikp_end))
411 CALL cp_fm_create(fm_v_kp_work, bs_env%fm_RI_RI%matrix_struct)
412 DO ikp = ikp_start, ikp_end
413 CALL cp_cfm_create(cfm_v_kp(ikp), bs_env%fm_RI_RI%matrix_struct)
414 CALL cp_fm_set_all(fm_v_kp_work, 0.0_dp)
415 CALL copy_dbcsr_to_fm(mat_v_kp(ikp, 1)%matrix, fm_v_kp_work)
416 cfm_v_kp(ikp)%local_data = cmplx(fm_v_kp_work%local_data, 0.0_dp, kind=dp)
417
418 CALL cp_fm_set_all(fm_v_kp_work, 0.0_dp)
419 CALL copy_dbcsr_to_fm(mat_v_kp(ikp, 2)%matrix, fm_v_kp_work)
420 cfm_v_kp(ikp)%local_data = cmplx(real(cfm_v_kp(ikp)%local_data, kind=dp), &
421 fm_v_kp_work%local_data, kind=dp)
422 END DO
423 CALL cp_fm_release(fm_v_kp_work)
424 CALL dbcsr_deallocate_matrix_set(mat_v_kp)
425
426 CALL timestop(handle)
427
428 END SUBROUTINE compute_v_k_by_lattice_sum
429
430! **************************************************************************************************
431!> \brief ...
432!> \param bs_env ...
433!> \param qs_env ...
434!> \param cfm_V_kp ...
435!> \param cfm_V_sqrt_ikp ...
436!> \param cfm_M_inv_V_sqrt_ikp ...
437!> \param ikp ...
438! **************************************************************************************************
439 SUBROUTINE compute_minvvsqrt_vsqrt(bs_env, qs_env, cfm_V_kp, cfm_V_sqrt_ikp, &
440 cfm_M_inv_V_sqrt_ikp, ikp)
441 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
442 TYPE(qs_environment_type), POINTER :: qs_env
443 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_v_kp
444 TYPE(cp_cfm_type) :: cfm_v_sqrt_ikp, cfm_m_inv_v_sqrt_ikp
445 INTEGER :: ikp
446
447 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_MinvVsqrt_Vsqrt'
448
449 INTEGER :: handle, info, n_ri
450 TYPE(cp_cfm_type) :: cfm_m_inv_ikp, cfm_work
451 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_m_ikp
452
453 CALL timeset(routinen, handle)
454
455 n_ri = bs_env%n_RI
456
457 ! get here M(k) and write it to cfm_M_ikp
458 CALL ri_2c_integral_mat_kp(qs_env, cfm_m_ikp, cfm_v_kp(ikp)%matrix_struct, n_ri, bs_env%ri_metric, &
459 bs_env%kpoints_chi_eps_W, regularization_ri=bs_env%regularization_RI, &
460 ikp_ext=ikp, do_build_cell_index=(ikp == 1))
461
462 IF (ikp == 1) THEN
463 CALL cp_cfm_create(cfm_v_sqrt_ikp, cfm_v_kp(ikp)%matrix_struct)
464 CALL cp_cfm_create(cfm_m_inv_v_sqrt_ikp, cfm_v_kp(ikp)%matrix_struct)
465 END IF
466 CALL cp_cfm_create(cfm_m_inv_ikp, cfm_v_kp(ikp)%matrix_struct)
467
468 CALL cp_cfm_to_cfm(cfm_m_ikp(1), cfm_m_inv_ikp)
469 CALL cp_cfm_to_cfm(cfm_v_kp(ikp), cfm_v_sqrt_ikp)
470
471 CALL cp_cfm_release(cfm_m_ikp(1))
472 DEALLOCATE (cfm_m_ikp)
473
474 CALL cp_cfm_create(cfm_work, cfm_v_kp(ikp)%matrix_struct)
475
476 ! M(k) -> M^-1(k)
477 CALL cp_cfm_to_cfm(cfm_m_inv_ikp, cfm_work)
478 CALL cp_cfm_cholesky_decompose(matrix=cfm_m_inv_ikp, n=n_ri, info_out=info)
479 IF (info == 0) THEN
480 ! successful Cholesky decomposition
481 CALL cp_cfm_cholesky_invert(cfm_m_inv_ikp)
482 ! symmetrize the result
483 CALL cp_cfm_uplo_to_full(cfm_m_inv_ikp)
484 ELSE
485 ! Cholesky decomposition not successful: use expensive diagonalization
486 CALL cp_cfm_power(cfm_work, threshold=bs_env%eps_eigval_mat_RI, exponent=-1.0_dp)
487 CALL cp_cfm_to_cfm(cfm_work, cfm_m_inv_ikp)
488 END IF
489
490 ! V(k) -> L(k) with L^H(k)*L(k) = V(k) [L(k) can be just considered to be V^0.5(k)]
491 CALL cp_cfm_to_cfm(cfm_v_sqrt_ikp, cfm_work)
492 CALL cp_cfm_cholesky_decompose(matrix=cfm_v_sqrt_ikp, n=n_ri, info_out=info)
493 IF (info == 0) THEN
494 ! successful Cholesky decomposition
495 CALL clean_lower_part(cfm_v_sqrt_ikp)
496 ELSE
497 ! Cholesky decomposition not successful: use expensive diagonalization
498 CALL cp_cfm_power(cfm_work, threshold=0.0_dp, exponent=0.5_dp)
499 CALL cp_cfm_to_cfm(cfm_work, cfm_v_sqrt_ikp)
500 END IF
501 CALL cp_cfm_release(cfm_work)
502
503 ! get M^-1(k)*V^0.5(k)
504 CALL parallel_gemm("N", "C", n_ri, n_ri, n_ri, z_one, cfm_m_inv_ikp, cfm_v_sqrt_ikp, &
505 z_zero, cfm_m_inv_v_sqrt_ikp)
506
507 CALL cp_cfm_release(cfm_m_inv_ikp)
508
509 CALL timestop(handle)
510
511 END SUBROUTINE compute_minvvsqrt_vsqrt
512
513! **************************************************************************************************
514!> \brief Compute the screened interaction for one frequency and k-point.
515!> \param bs_env band-structure environment
516!> \param cfm_chi_eps_W_ikp_freq_j input chi; overwritten by W
517!> \param cfm_V_sqrt_ikp square root of the Coulomb matrix
518!> \param cfm_M_inv_V_sqrt_ikp metric-weighted Coulomb factor
519! **************************************************************************************************
520 SUBROUTINE compute_cfm_w_ikp_freq_j(bs_env, cfm_chi_eps_W_ikp_freq_j, cfm_V_sqrt_ikp, &
521 cfm_M_inv_V_sqrt_ikp)
522 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
523 TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_chi_eps_w_ikp_freq_j
524 TYPE(cp_cfm_type), INTENT(IN) :: cfm_v_sqrt_ikp, cfm_m_inv_v_sqrt_ikp
525
526 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_cfm_W_ikp_freq_j'
527
528 INTEGER :: handle, info
529
530 CALL timeset(routinen, handle)
531
532 ! 1. epsilon(iw,k) = Id - V^0.5(k)*M^-1(k)*chi(iw,k)*M^-1(k)*V^0.5(k)
533
534 ! 1. a) Form V^0.5(k)*M^-1(k)*chi(iw,k)*M^-1(k)*V^0.5(k).
535 CALL cfm_contract_aba(cfm_m_inv_v_sqrt_ikp, cfm_chi_eps_w_ikp_freq_j)
536
537 ! 1. b) Add the identity matrix.
538 CALL cp_cfm_add_on_diag(cfm_chi_eps_w_ikp_freq_j, z_one)
539
540 ! 2. W(iw,k) = V^0.5(k)*(epsilon^-1(iw,k)-Id)*V^0.5(k)
541
542 ! 2. a) Cholesky-factorize epsilon before inversion.
543 CALL cp_cfm_cholesky_decompose(matrix=cfm_chi_eps_w_ikp_freq_j, n=bs_env%n_RI, info_out=info)
544 cpassert(info == 0)
545
546 ! 2. b) Invert epsilon, subtract Id, then contract with V^0.5 on both sides.
547 CALL cfm_screened_interaction_from_eps_cholesky(cfm_chi_eps_w_ikp_freq_j, cfm_v_sqrt_ikp)
548
549 CALL timestop(handle)
550
551 END SUBROUTINE compute_cfm_w_ikp_freq_j
552
553! **************************************************************************************************
554!> \brief ...
555!> \param bs_env ...
556!> \param mat_chi_Gamma_tau ...
557!> \param fm_W_MIC_time ...
558! **************************************************************************************************
559 SUBROUTINE read_w_mic_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
560 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
561 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
562 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
563
564 CHARACTER(LEN=*), PARAMETER :: routinen = 'read_W_MIC_time'
565
566 INTEGER :: handle, i_t
567 REAL(kind=dp) :: t1
568
569 CALL timeset(routinen, handle)
570
571 CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
572 CALL create_fm_w_mic_time(bs_env, fm_w_mic_time)
573
574 DO i_t = 1, bs_env%num_time_freq_points
575
576 t1 = m_walltime()
577
578 CALL fm_read(fm_w_mic_time(i_t), bs_env, bs_env%W_time_name, i_t)
579
580 IF (bs_env%unit_nr > 0) THEN
581 WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
582 'Read W^MIC(iτ) from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
583 ', Execution time', m_walltime() - t1, ' s'
584 END IF
585
586 END DO
587
588 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
589
590 ! W^MIC(iω=0) for the post-GW consumers, built here because fm_W_MIC_time does not outlive
591 ! the W stage: RT-BSE reads its own restart file, PROPERTIES%BSE sums the W(iτ) just read.
592 IF (bs_env%rtp_method == rtp_method_bse .OR. &
593 bs_env%rtp_method == rtp_method_bse_linearized) THEN
594 CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
595 t1 = m_walltime()
596 CALL fm_read(bs_env%fm_W_MIC_freq_zero, bs_env, "W_freq_rtp", 0)
597 IF (bs_env%unit_nr > 0) THEN
598 WRITE (bs_env%unit_nr, '(T2,A,I3,A,I3,A,F10.1,A)') &
599 'Read W^MIC(f=0) from file for freq. point ', 1, ' /', 1, &
600 ', Execution time', m_walltime() - t1, ' s'
601 END IF
602 ELSE IF (bs_env%keep_W_freq_zero) THEN
603 ! no RTP file in a LRBSE-only run: W(iτ) was just read from the W_time files
604 CALL sum_w_mic_freq_zero(bs_env, fm_w_mic_time)
605 END IF
606
607 CALL timestop(handle)
608
609 END SUBROUTINE read_w_mic_time
610
611! **************************************************************************************************
612!> \brief ...
613!> \param bs_env ...
614!> \param qs_env ...
615!> \param fm_W_MIC_freq_j ...
616!> \param j_w ...
617!> \param ikp ...
618!> \param mat_chi_Gamma_tau ...
619!> \param cfm_M_inv_V_sqrt_ikp ...
620!> \param cfm_V_sqrt_ikp ...
621! **************************************************************************************************
622 SUBROUTINE compute_fm_w_mic_freq_j(bs_env, qs_env, fm_W_MIC_freq_j, j_w, ikp, mat_chi_Gamma_tau, &
623 cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
624 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
625 TYPE(qs_environment_type), POINTER :: qs_env
626 TYPE(cp_fm_type) :: fm_w_mic_freq_j
627 INTEGER :: j_w, ikp
628 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
629 TYPE(cp_cfm_type) :: cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp
630
631 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_fm_W_MIC_freq_j'
632
633 INTEGER :: handle
634
635 CALL timeset(routinen, handle)
636
637 ! 1. Fourier transformation of χ_PQ(iτ,k=0) to χ_PQ(iω_j,k=0)
638 CALL compute_fm_chi_gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_gamma_tau)
639
640 ! 2. χ_PQ(iω_j,k=0) -> W_PQ^MIC(iω_j)
641 CALL compute_fm_w_mic_freq_j_from_chi(bs_env, qs_env, fm_w_mic_freq_j, j_w, ikp, &
642 cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp)
643
644 CALL timestop(handle)
645
646 END SUBROUTINE compute_fm_w_mic_freq_j
647
648! **************************************************************************************************
649!> \brief W_PQ^MIC(iω_j) from χ_PQ(iω_j,k=0) in bs_env%fm_chi_Gamma_freq
650!> \param bs_env ...
651!> \param qs_env ...
652!> \param fm_W_MIC_freq_j ...
653!> \param j_w ...
654!> \param ikp ...
655!> \param cfm_M_inv_V_sqrt_ikp ...
656!> \param cfm_V_sqrt_ikp ...
657! **************************************************************************************************
658 SUBROUTINE compute_fm_w_mic_freq_j_from_chi(bs_env, qs_env, fm_W_MIC_freq_j, j_w, ikp, &
659 cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
660 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
661 TYPE(qs_environment_type), POINTER :: qs_env
662 TYPE(cp_fm_type) :: fm_w_mic_freq_j
663 INTEGER :: j_w, ikp
664 TYPE(cp_cfm_type) :: cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp
665
666 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_fm_W_MIC_freq_j_from_chi'
667
668 INTEGER :: handle
669 TYPE(cp_cfm_type) :: cfm_chi_eps_w_ikp_freq_j
670
671 CALL timeset(routinen, handle)
672
673 CALL cp_fm_set_all(fm_w_mic_freq_j, 0.0_dp)
674
675 ! 1. Get χ_PQ(iω_j,k_i) from χ_PQ(iω_j,k=0) using the minimum image convention
676 CALL cfm_ikp_from_fm_gamma(cfm_chi_eps_w_ikp_freq_j, bs_env%fm_chi_Gamma_freq, &
677 ikp, qs_env, bs_env%kpoints_chi_eps_W, "RI_AUX")
678
679 ! 2. Remove all negative eigenvalues from χ_PQ(iω_j,k_i)
680 CALL cp_cfm_power(cfm_chi_eps_w_ikp_freq_j, threshold=0.0_dp, exponent=1.0_dp)
681
682 ! 3. ε(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)
683 CALL compute_cfm_w_ikp_freq_j(bs_env, cfm_chi_eps_w_ikp_freq_j, cfm_v_sqrt_ikp, &
684 cfm_m_inv_v_sqrt_ikp)
685
686 ! 4. k-point integration W_PQ(iω_j, k_i) to W_PQ^MIC(iω_j)
687 SELECT CASE (bs_env%approx_kp_extrapol)
688 CASE (.false.)
689 ! default: standard k-point extrapolation
690 CALL mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, &
691 cfm_chi_eps_w_ikp_freq_j, ikp, &
692 bs_env%kpoints_chi_eps_W, "RI_AUX")
693 CASE (.true.)
694 ! for approximate kpoint extrapolation: get W_PQ^MIC(iω_1) with and without k-point
695 ! extrapolation to compute the extrapolation factor f_PQ for every PQ-matrix element,
696 ! f_PQ = (W_PQ^MIC(iω_1) with extrapolation) / (W_PQ^MIC(iω_1) without extrapolation)
697
698 ! for ω_1, we compute the k-point extrapolated result using all k-points
699 IF (j_w == 1) THEN
700
701 ! k-point extrapolated
702 CALL mic_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_extra, &
703 cfm_chi_eps_w_ikp_freq_j, ikp, &
704 bs_env%kpoints_chi_eps_W, &
705 "RI_AUX")
706 ! non-kpoint extrapolated
707 IF (ikp <= bs_env%nkp_chi_eps_W_orig) THEN
708 CALL mic_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_no_extra, &
709 cfm_chi_eps_w_ikp_freq_j, ikp, &
710 bs_env%kpoints_chi_eps_W, &
711 "RI_AUX", wkp_ext=bs_env%wkp_orig)
712 END IF
713
714 END IF
715
716 ! for all ω_j, we need to compute W^MIC without k-point extrpolation
717 IF (ikp <= bs_env%nkp_chi_eps_W_orig) THEN
718 CALL mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, &
719 cfm_chi_eps_w_ikp_freq_j, &
720 ikp, bs_env%kpoints_chi_eps_W, "RI_AUX", &
721 wkp_ext=bs_env%wkp_orig)
722 END IF
723 END SELECT
724
725 CALL cp_cfm_release(cfm_chi_eps_w_ikp_freq_j)
726
727 CALL timestop(handle)
728
729 END SUBROUTINE compute_fm_w_mic_freq_j_from_chi
730
731! **************************************************************************************************
732!> \brief ...
733!> \param cfm_mat ...
734! **************************************************************************************************
735 SUBROUTINE clean_lower_part(cfm_mat)
736 TYPE(cp_cfm_type) :: cfm_mat
737
738 CHARACTER(LEN=*), PARAMETER :: routinen = 'clean_lower_part'
739
740 INTEGER :: handle, i_row, j_col, j_global, &
741 ncol_local, nrow_local
742 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
743
744 CALL timeset(routinen, handle)
745
746 CALL cp_cfm_get_info(matrix=cfm_mat, &
747 nrow_local=nrow_local, ncol_local=ncol_local, &
748 row_indices=row_indices, col_indices=col_indices)
749
750 DO j_col = 1, ncol_local
751 j_global = col_indices(j_col)
752 DO i_row = 1, nrow_local
753 IF (j_global < row_indices(i_row)) cfm_mat%local_data(i_row, j_col) = z_zero
754 END DO
755 END DO
756
757 CALL timestop(handle)
758
759 END SUBROUTINE clean_lower_part
760
761! **************************************************************************************************
762!> \brief ...
763!> \param bs_env ...
764!> \param fm_W_MIC_time ...
765! **************************************************************************************************
766 SUBROUTINE apply_extrapol_factor(bs_env, fm_W_MIC_time)
767 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
768 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
769
770 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_extrapol_factor'
771 REAL(kind=dp), PARAMETER :: eps_w_no_extra_1 = 1.0e-13_dp
772
773 INTEGER :: handle, i, i_t, j, ncol_local, nrow_local
774 REAL(kind=dp) :: extrapol_factor, w_extra_1, w_no_extra_1
775
776 CALL timeset(routinen, handle)
777
778 CALL cp_fm_get_info(matrix=fm_w_mic_time(1), nrow_local=nrow_local, ncol_local=ncol_local)
779
780 DO i_t = 1, SIZE(fm_w_mic_time)
781 DO j = 1, ncol_local
782 DO i = 1, nrow_local
783
784 w_extra_1 = bs_env%fm_W_MIC_freq_1_extra%local_data(i, j)
785 w_no_extra_1 = bs_env%fm_W_MIC_freq_1_no_extra%local_data(i, j)
786
787 IF (abs(w_no_extra_1) > eps_w_no_extra_1) THEN
788 extrapol_factor = abs(w_extra_1/w_no_extra_1)
789 ELSE
790 extrapol_factor = 1.0_dp
791 END IF
792
793 ! reset extrapolation factor if it is very large
794 IF (extrapol_factor > 10.0_dp) extrapol_factor = 1.0_dp
795
796 fm_w_mic_time(i_t)%local_data(i, j) = fm_w_mic_time(i_t)%local_data(i, j) &
797 *extrapol_factor
798 END DO
799 END DO
800 END DO
801
802 CALL timestop(handle)
803
804 END SUBROUTINE apply_extrapol_factor
805
806! **************************************************************************************************
807!> \brief Computes V, V^0.5, and M^-1 V^0.5 for the RI-RS dielectric function.
808!> The Coulomb matrix V is constructed by the RI-RS k-point path. The inverse metric
809!> M^-1(k=0) is precomputed once in gw_utils and stored in bs_env.
810!> \param bs_env ...
811!> \param qs_env ...
812!> \param fm_V Coulomb matrix V(k=0)
813!> \param fm_V_sqrt symmetric factor V^0.5
814!> \param fm_Minv_Vsqrt product M^-1(k=0) V^0.5
815! **************************************************************************************************
816 SUBROUTINE compute_v_minvvsqrt(bs_env, qs_env, fm_V, fm_V_sqrt, fm_Minv_Vsqrt)
817 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
818 TYPE(qs_environment_type), POINTER :: qs_env
819 TYPE(cp_fm_type), INTENT(INOUT) :: fm_v, fm_v_sqrt, fm_minv_vsqrt
820
821 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_V_MinvVsqrt'
822
823 INTEGER :: handle
824 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
825 TYPE(cell_type), POINTER :: cell
826 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_v_kp
827 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
828 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
829
830 CALL timeset(routinen, handle)
831
832 IF (bs_env%auto_ri%enabled) THEN
833 ! -------------------------------------------------------------------
834 ! 1a. The optimized AB functions span two atoms, so their previously transformed
835 ! Coulomb matrix cannot be rebuilt by the atom-local k-point integral routine.
836 ! -------------------------------------------------------------------
837 CALL cp_fm_to_fm(bs_env%auto_ri%V_pq, fm_v)
838 ELSE
839 ! -------------------------------------------------------------------
840 ! 1b. Build Coulomb Matrix V(k=0) using the kp-routine but only for ikp=1.
841 ! -------------------------------------------------------------------
842 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, &
843 qs_kind_set=qs_kind_set)
844 atomic_kind_set => bs_env%ri_rs%atomic_kind_set
845
846 ALLOCATE (mat_v_kp(1:1, 1:2))
847 NULLIFY (mat_v_kp(1, 1)%matrix, mat_v_kp(1, 2)%matrix)
848 ALLOCATE (mat_v_kp(1, 1)%matrix, mat_v_kp(1, 2)%matrix)
849 CALL dbcsr_create(mat_v_kp(1, 1)%matrix, template=bs_env%mat_RI_RI%matrix)
850 CALL dbcsr_reserve_all_blocks(mat_v_kp(1, 1)%matrix)
851 CALL dbcsr_set(mat_v_kp(1, 1)%matrix, 0.0_dp)
852 CALL dbcsr_create(mat_v_kp(1, 2)%matrix, template=bs_env%mat_RI_RI%matrix)
853 CALL dbcsr_reserve_all_blocks(mat_v_kp(1, 2)%matrix)
854 ! The dummy imaginary part is required only by the k-point routine interface.
855 CALL dbcsr_set(mat_v_kp(1, 2)%matrix, 0.0_dp)
856
857 bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
858 CALL build_2c_coulomb_matrix_kp(mat_v_kp, bs_env%kpoints_chi_eps_W, "RI_AUX", cell, &
859 particle_set, qs_kind_set, atomic_kind_set, &
860 bs_env%size_lattice_sum_V, operator_coulomb, 1, 1)
861 CALL copy_dbcsr_to_fm(mat_v_kp(1, 1)%matrix, fm_v)
862
863 CALL dbcsr_deallocate_matrix(mat_v_kp(1, 1)%matrix)
864 CALL dbcsr_deallocate_matrix(mat_v_kp(1, 2)%matrix)
865 DEALLOCATE (mat_v_kp)
866 END IF
867
868 ! -----------------------------------------------------------------------
869 ! 2. V -> V^0.5.
870 ! -----------------------------------------------------------------------
871 CALL fm_sqrt(fm_v, fm_v_sqrt, bs_env%eps_eigval_mat_RI, bs_env%unit_nr)
872
873 ! -----------------------------------------------------------------------
874 ! 3. M^-1(k=0) V^0.5.
875 ! -----------------------------------------------------------------------
876 CALL parallel_gemm("N", "T", bs_env%n_RI, bs_env%n_RI, bs_env%n_RI, 1.0_dp, &
877 bs_env%fm_Minv_Gamma, fm_v_sqrt, 0.0_dp, fm_minv_vsqrt)
878
879 CALL timestop(handle)
880
881 END SUBROUTINE compute_v_minvvsqrt
882
883! **************************************************************************************************
884!> \brief Computes the screened interaction at one imaginary frequency:
885!> ε(iω_j) = Id - (M^-1 V^0.5)^T χ(iω_j) (M^-1 V^0.5)
886!> W(iω_j) = V^0.5^T ( ε^-1(iω_j) - Id ) V^0.5
887!> ε is inverted via Cholesky; if that fails due to conditioning, via
888!> eigendecomposition (cp_fm_power) with eigenvalue filtering.
889!> \param bs_env ...
890!> \param fm_chi_freq_j ...
891!> \param fm_V_sqrt ...
892!> \param fm_Minv_Vsqrt ...
893!> \param fm_W_freq_j ...
894! **************************************************************************************************
895 SUBROUTINE compute_fm_w_freq(bs_env, fm_chi_freq_j, fm_V_sqrt, fm_Minv_Vsqrt, fm_W_freq_j)
896 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
897 TYPE(cp_fm_type), INTENT(IN) :: fm_chi_freq_j, fm_v_sqrt, fm_minv_vsqrt
898 TYPE(cp_fm_type), INTENT(INOUT) :: fm_w_freq_j
899
900 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_fm_W_freq'
901
902 INTEGER :: handle, n_ri
903 TYPE(cp_fm_type) :: fm_eps_freq_j, fm_work
904
905 CALL timeset(routinen, handle)
906
907 n_ri = bs_env%n_RI
908
909 CALL cp_fm_create(fm_eps_freq_j, fm_chi_freq_j%matrix_struct)
910 CALL cp_fm_create(fm_work, fm_chi_freq_j%matrix_struct)
911
912 ! -----------------------------------------------------------------------
913 ! 1. ε(iω_j) = Id - (M^-1 * V^0.5)^T * χ(iω_j) * (M^-1 * V^0.5)
914 ! -----------------------------------------------------------------------
915 ! work = χ(iω_j) * (M^-1 * V^0.5)
916 CALL parallel_gemm('N', 'N', n_ri, n_ri, n_ri, 1.0_dp, &
917 fm_chi_freq_j, fm_minv_vsqrt, 0.0_dp, fm_work)
918
919 ! eps_work = (M^-1 * V^0.5)^T * work
920 CALL parallel_gemm('T', 'N', n_ri, n_ri, n_ri, 1.0_dp, &
921 fm_minv_vsqrt, fm_work, 0.0_dp, fm_eps_freq_j)
922
923 ! ε(iω_j) = Id - eps_work --> -eps_work + Id
924 CALL cp_fm_add_on_diag(fm_eps_freq_j, 1.0_dp)
925
926 ! Force perfect symmetry before Cholesky to avoid info != 0 due to GEMM noise
927 CALL cp_fm_uplo_to_full(fm_eps_freq_j, fm_work)
928
929 ! -----------------------------------------------------------------------
930 ! 2. W(iω_j) = V^0.5^T * (ε^-1(iω_j) - Id) * V^0.5
931 ! -----------------------------------------------------------------------
932
933 ! a) Invert ε by Cholesky decomposition or, if that fails, by diagonalization.
934 CALL fm_invert(fm_eps_freq_j, bs_env%eps_eigval_mat_RI, bs_env%unit_nr)
935
936 ! b) ε^-1(iω_j) - Id
937 CALL cp_fm_add_on_diag(fm_eps_freq_j, -1.0_dp)
938
939 ! c) W(iω_j) = V^0.5^T * (ε^-1(iω_j) - Id) * V^0.5
940 CALL fm_contract_aba(fm_v_sqrt, fm_eps_freq_j, fm_w_freq_j)
941
942 ! Cleanup
943 CALL cp_fm_release(fm_work)
944 CALL cp_fm_release(fm_eps_freq_j)
945
946 CALL timestop(handle)
947
948 END SUBROUTINE compute_fm_w_freq
949
950! **************************************************************************************************
951!> \brief ...
952!> \param fm ...
953!> \param bs_env ...
954!> \param mat_name ...
955!> \param idx ...
956! **************************************************************************************************
957 SUBROUTINE fm_read(fm, bs_env, mat_name, idx)
958 TYPE(cp_fm_type) :: fm
959 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
960 CHARACTER(LEN=*) :: mat_name
961 INTEGER :: idx
962
963 CHARACTER(LEN=*), PARAMETER :: routinen = 'fm_read'
964
965 CHARACTER(LEN=default_path_length) :: f_chi
966 INTEGER :: handle, unit_nr
967
968 CALL timeset(routinen, handle)
969
970 unit_nr = -1
971 IF (bs_env%para_env%is_source()) THEN
972
973 IF (idx < 10) THEN
974 WRITE (f_chi, '(3A,I1,A)') trim(bs_env%prefix), trim(mat_name), "_0", idx, ".matrix"
975 ELSE IF (idx < 100) THEN
976 WRITE (f_chi, '(3A,I2,A)') trim(bs_env%prefix), trim(mat_name), "_", idx, ".matrix"
977 ELSE
978 cpabort('Please implement more than 99 time/frequency points.')
979 END IF
980
981 CALL open_file(file_name=trim(f_chi), file_action="READ", file_form="UNFORMATTED", &
982 file_position="REWIND", file_status="OLD", unit_number=unit_nr)
983
984 END IF
985
986 CALL cp_fm_read_unformatted(fm, unit_nr)
987
988 IF (bs_env%para_env%is_source()) CALL close_file(unit_number=unit_nr)
989
990 CALL timestop(handle)
991
992 END SUBROUTINE fm_read
993
994! **************************************************************************************************
995!> \brief Write a GW restart matrix when requested by the print key.
996!> \param fm Matrix to write, or reusable conversion buffer when matrix is supplied.
997!> \param matrix_index ...
998!> \param matrix_name ...
999!> \param qs_env ...
1000!> \param matrix Optional DBCSR source, converted to fm only when output is enabled.
1001! **************************************************************************************************
1002 SUBROUTINE fm_write(fm, matrix_index, matrix_name, qs_env, matrix)
1003 TYPE(cp_fm_type) :: fm
1004 INTEGER :: matrix_index
1005 CHARACTER(LEN=*) :: matrix_name
1006 TYPE(qs_environment_type), POINTER :: qs_env
1007 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix
1008
1009 CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
1010 routinen = 'fm_write'
1011
1012 CHARACTER(LEN=default_path_length) :: filename
1013 INTEGER :: handle, unit_nr
1014 TYPE(cp_logger_type), POINTER :: logger
1015 TYPE(section_vals_type), POINTER :: input
1016
1017 CALL timeset(routinen, handle)
1018
1019 CALL get_qs_env(qs_env, input=input)
1020
1021 logger => cp_get_default_logger()
1022
1023 IF (btest(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
1024
1025 IF (PRESENT(matrix)) THEN
1026 CALL cp_fm_set_all(fm, 0.0_dp)
1027 CALL copy_dbcsr_to_fm(matrix, fm)
1028 END IF
1029
1030 IF (matrix_index < 10) THEN
1031 WRITE (filename, '(3A,I1)') "RESTART_", matrix_name, "_0", matrix_index
1032 ELSE IF (matrix_index < 100) THEN
1033 WRITE (filename, '(3A,I2)') "RESTART_", matrix_name, "_", matrix_index
1034 ELSE
1035 cpabort('Please implement more than 99 time/frequency points.')
1036 END IF
1037
1038 unit_nr = cp_print_key_unit_nr(logger, input, key, extension=".matrix", &
1039 file_form="UNFORMATTED", middle_name=trim(filename), &
1040 file_position="REWIND", file_action="WRITE")
1041
1042 CALL cp_fm_write_unformatted(fm, unit_nr)
1043 IF (unit_nr > 0) THEN
1044 CALL close_file(unit_nr)
1045 END IF
1046 END IF
1047
1048 CALL timestop(handle)
1049
1050 END SUBROUTINE fm_write
1051
1052! **************************************************************************************************
1053!> \brief ...
1054!> \param bs_env ...
1055!> \param fm_chi_Gamma_freq ...
1056!> \param j_w ...
1057!> \param mat_chi_Gamma_tau ...
1058! **************************************************************************************************
1059 SUBROUTINE compute_fm_chi_gamma_freq(bs_env, fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
1060 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1061 TYPE(cp_fm_type) :: fm_chi_gamma_freq
1062 INTEGER :: j_w
1063 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_gamma_tau
1064
1065 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_fm_chi_Gamma_freq'
1066
1067 INTEGER :: handle, i_t
1068 REAL(kind=dp) :: weight_ij
1069
1070 CALL timeset(routinen, handle)
1071
1072 CALL dbcsr_set(bs_env%mat_RI_RI%matrix, 0.0_dp)
1073
1074 DO i_t = 1, bs_env%num_time_freq_points
1075
1076 weight_ij = bs_env%time_frequency_grid%cosine_time_to_frequency(j_w, i_t)
1077
1078 ! Fourier transform
1079 CALL dbcsr_add(bs_env%mat_RI_RI%matrix, mat_chi_gamma_tau(i_t)%matrix, &
1080 1.0_dp, weight_ij)
1081
1082 END DO
1083
1084 CALL copy_dbcsr_to_fm(bs_env%mat_RI_RI%matrix, fm_chi_gamma_freq)
1085
1086 CALL timestop(handle)
1087
1088 END SUBROUTINE compute_fm_chi_gamma_freq
1089
1090! **************************************************************************************************
1091!> \brief W^MIC_PQ(iω=0) = sum_j w_j W^MIC_PQ(iτ_j) with the time weights w_j of the grid at zero
1092!> frequency, into bs_env%fm_W_MIC_freq_zero
1093!> \param bs_env time grid and its weights; receives fm_W_MIC_freq_zero
1094!> \param fm_W_MIC_time W^MIC_PQ(iτ_j) for every time point j
1095! **************************************************************************************************
1096 SUBROUTINE sum_w_mic_freq_zero(bs_env, fm_W_MIC_time)
1097 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1098 TYPE(cp_fm_type), DIMENSION(:) :: fm_w_mic_time
1099
1100 CHARACTER(LEN=*), PARAMETER :: routinen = 'sum_W_MIC_freq_zero'
1101
1102 INTEGER :: handle, i_t
1103 REAL(kind=dp) :: t_start
1104
1105 CALL timeset(routinen, handle)
1106
1107 t_start = m_walltime()
1108 CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
1109 CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
1110 DO i_t = 1, bs_env%num_time_freq_points
1111 CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
1112 bs_env%time_frequency_grid%time_weights_at_zero_frequency(i_t), &
1113 fm_w_mic_time(i_t))
1114 END DO
1115 IF (bs_env%unit_nr > 0) THEN
1116 WRITE (bs_env%unit_nr, '(T2,A,I11,A,I3,A,F10.1,A)') &
1117 'Computed W(ω=0,k) for k-point batch', &
1118 1, ' /', 1, &
1119 ', Execution time', m_walltime() - t_start, ' s'
1120 END IF
1121
1122 CALL timestop(handle)
1123
1124 END SUBROUTINE sum_w_mic_freq_zero
1125
1126! **************************************************************************************************
1127!> \brief ...
1128!> \param bs_env ...
1129!> \param fm_W_MIC_time ...
1130! **************************************************************************************************
1131 SUBROUTINE create_fm_w_mic_time(bs_env, fm_W_MIC_time)
1132 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1133 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
1134
1135 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_fm_W_MIC_time'
1136
1137 INTEGER :: handle, i_t
1138
1139 CALL timeset(routinen, handle)
1140
1141 ALLOCATE (fm_w_mic_time(bs_env%num_time_freq_points))
1142 DO i_t = 1, bs_env%num_time_freq_points
1143 CALL cp_fm_create(fm_w_mic_time(i_t), bs_env%fm_RI_RI%matrix_struct, set_zero=.true.)
1144 END DO
1145
1146 CALL timestop(handle)
1147
1148 END SUBROUTINE create_fm_w_mic_time
1149
1150! **************************************************************************************************
1151!> \brief ...
1152!> \param bs_env ...
1153!> \param fm_W_MIC_time ...
1154!> \param fm_W_MIC_freq_j ...
1155!> \param j_w ...
1156! **************************************************************************************************
1157 SUBROUTINE fourier_transform_w_to_t(bs_env, fm_W_MIC_time, fm_W_MIC_freq_j, j_w)
1158 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1159 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_w_mic_time
1160 TYPE(cp_fm_type) :: fm_w_mic_freq_j
1161 INTEGER :: j_w
1162
1163 CHARACTER(LEN=*), PARAMETER :: routinen = 'Fourier_transform_w_to_t'
1164
1165 INTEGER :: handle, i_t
1166 REAL(kind=dp) :: weight_ij
1167
1168 CALL timeset(routinen, handle)
1169
1170 DO i_t = 1, bs_env%num_time_freq_points
1171
1172 weight_ij = bs_env%time_frequency_grid%cosine_frequency_to_time(i_t, j_w)
1173
1174 ! actual Fourier transform
1175 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_w_mic_time(i_t), &
1176 beta=weight_ij, matrix_b=fm_w_mic_freq_j)
1177
1178 END DO
1179
1180 CALL timestop(handle)
1181
1182 END SUBROUTINE fourier_transform_w_to_t
1183
1184! **************************************************************************************************
1185!> \brief ...
1186!> \param V_tr_R ...
1187!> \param pot_type ...
1188!> \param regularization_RI ...
1189!> \param bs_env ...
1190!> \param qs_env ...
1191! **************************************************************************************************
1192 SUBROUTINE get_v_tr_r(V_tr_R, pot_type, regularization_RI, bs_env, qs_env)
1193 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: v_tr_r
1194 TYPE(libint_potential_type) :: pot_type
1195 REAL(kind=dp) :: regularization_ri
1196 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1197 TYPE(qs_environment_type), POINTER :: qs_env
1198
1199 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_V_tr_R'
1200
1201 INTEGER :: handle, img, nimages_scf_desymm
1202 INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_ri
1203 INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
1204 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1205 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_v_tr_r
1206 TYPE(dbcsr_distribution_type) :: dbcsr_dist
1207 TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: mat_v_tr_r
1208 TYPE(distribution_2d_type), POINTER :: dist_2d
1209 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1210 POINTER :: sab_ri
1211 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1212 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1213
1214 CALL timeset(routinen, handle)
1215
1216 NULLIFY (sab_ri, dist_2d)
1217
1218 CALL get_qs_env(qs_env=qs_env, &
1219 blacs_env=blacs_env, &
1220 distribution_2d=dist_2d, &
1221 qs_kind_set=qs_kind_set, &
1222 particle_set=particle_set)
1223
1224 ALLOCATE (sizes_ri(bs_env%n_atom))
1225 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=bs_env%basis_set_RI)
1226 CALL build_2c_neighbor_lists(sab_ri, bs_env%basis_set_RI, bs_env%basis_set_RI, &
1227 pot_type, "2c_nl_RI", qs_env, sym_ij=.false., &
1228 dist_2d=dist_2d)
1229 CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
1230 ALLOCATE (row_bsize(SIZE(sizes_ri)))
1231 ALLOCATE (col_bsize(SIZE(sizes_ri)))
1232 row_bsize(:) = sizes_ri
1233 col_bsize(:) = sizes_ri
1234
1235 nimages_scf_desymm = bs_env%nimages_scf_desymm
1236 ALLOCATE (mat_v_tr_r(nimages_scf_desymm))
1237 CALL dbcsr_create(mat_v_tr_r(1), "(RI|RI)", dbcsr_dist, dbcsr_type_no_symmetry, &
1238 row_bsize, col_bsize)
1239 DEALLOCATE (row_bsize, col_bsize)
1240
1241 DO img = 2, nimages_scf_desymm
1242 CALL dbcsr_create(mat_v_tr_r(img), template=mat_v_tr_r(1))
1243 END DO
1244
1245 CALL build_2c_integrals(mat_v_tr_r, 0.0_dp, qs_env, sab_ri, bs_env%basis_set_RI, &
1246 bs_env%basis_set_RI, pot_type, do_kpoints=.true., &
1247 ext_kpoints=bs_env%kpoints_scf_desymm, &
1248 regularization_ri=regularization_ri)
1249
1250 ALLOCATE (fm_v_tr_r(nimages_scf_desymm))
1251 DO img = 1, nimages_scf_desymm
1252 CALL cp_fm_create(fm_v_tr_r(img), bs_env%fm_RI_RI%matrix_struct)
1253 CALL copy_dbcsr_to_fm(mat_v_tr_r(img), fm_v_tr_r(img))
1254 CALL dbcsr_release(mat_v_tr_r(img))
1255 END DO
1256
1257 IF (.NOT. ALLOCATED(v_tr_r)) THEN
1258 ALLOCATE (v_tr_r(bs_env%n_RI, bs_env%n_RI, nimages_scf_desymm))
1259 END IF
1260
1261 CALL cp_fm_to_replicated_array(fm_v_tr_r, v_tr_r)
1262
1263 CALL cp_fm_release(fm_v_tr_r)
1264 CALL dbcsr_distribution_release(dbcsr_dist)
1265 CALL release_neighbor_list_sets(sab_ri)
1266
1267 CALL timestop(handle)
1268
1269 END SUBROUTINE get_v_tr_r
1270
1271END MODULE gw_compute_v_and_w
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:162
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:103
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
gw_utils_fm::fm_contract_aba
Definition gw_utils_fm.F:44
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
constants_operator
constants for the different operators of the 2c-integrals
Definition constants_operator.F:14
constants_operator::operator_coulomb
integer, parameter, public operator_coulomb
Definition constants_operator.F:22
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_uplo_to_full
subroutine, public cp_cfm_uplo_to_full(matrix, workspace, uplo)
...
Definition cp_cfm_basic_linalg.F:1351
cp_cfm_basic_linalg::cp_cfm_add_on_diag
subroutine, public cp_cfm_add_on_diag(matrix, alpha, n_active)
Adds a scalar to the diagonal of a distributed complex full matrix.
Definition cp_cfm_basic_linalg.F:1448
cp_cfm_cholesky
various cholesky decomposition related routines
Definition cp_cfm_cholesky.F:13
cp_cfm_cholesky::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_cholesky.F:53
cp_cfm_cholesky::cp_cfm_cholesky_invert
subroutine, public cp_cfm_cholesky_invert(matrix, n, info_out)
Used to replace Cholesky decomposition by the inverse.
Definition cp_cfm_cholesky.F:124
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_power
subroutine, public cp_cfm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
Raise the positive spectral part of a Hermitian matrix to a real power. Eigenmodes with eigenvalues b...
Definition cp_cfm_diag.F:239
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:643
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_reserve_all_blocks
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
Definition cp_dbcsr_contrib.F:300
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:233
cp_dbcsr_operations::cp_dbcsr_dist2d_to_dist
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
Definition cp_dbcsr_operations.F:696
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_add_on_diag
subroutine, public cp_fm_add_on_diag(matrix, alpha, n_active)
Adds a scalar to the diagonal of a distributed real full matrix.
Definition cp_fm_basic_linalg.F:2941
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:171
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1585
cp_fm_redistribution
Redistribution of full matrices and distributed arrays.
Definition cp_fm_redistribution.F:11
cp_fm_redistribution::cp_fm_to_replicated_array
subroutine, public cp_fm_to_replicated_array(fm_s, array_s, weight, add)
...
Definition cp_fm_redistribution.F:282
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1045
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition cp_fm_types.F:2177
cp_fm_types::cp_fm_read_unformatted
subroutine, public cp_fm_read_unformatted(fm, unit)
...
Definition cp_fm_types.F:2408
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:166
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:838
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:104
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:330
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition distribution_2d_types.F:15
gw_compute_v_and_w
Coulomb matrices and screened interactions for imaginary-time GW.
Definition gw_compute_V_and_W.F:11
gw_compute_v_and_w::fm_write
subroutine, public fm_write(fm, matrix_index, matrix_name, qs_env, matrix)
Write a GW restart matrix when requested by the print key.
Definition gw_compute_V_and_W.F:1003
gw_compute_v_and_w::compute_minvvsqrt_vsqrt
subroutine, public compute_minvvsqrt_vsqrt(bs_env, qs_env, cfm_v_kp, cfm_v_sqrt_ikp, cfm_m_inv_v_sqrt_ikp, ikp)
...
Definition gw_compute_V_and_W.F:441
gw_compute_v_and_w::get_w_mic
subroutine, public get_w_mic(bs_env, qs_env, mat_chi_gamma_tau, fm_w_mic_time)
...
Definition gw_compute_V_and_W.F:204
gw_compute_v_and_w::compute_fm_w_mic_freq_j_from_chi
subroutine, public compute_fm_w_mic_freq_j_from_chi(bs_env, qs_env, fm_w_mic_freq_j, j_w, ikp, cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp)
W_PQ^MIC(iω_j) from χ_PQ(iω_j,k=0) in bs_envfm_chi_Gamma_freq.
Definition gw_compute_V_and_W.F:660
gw_compute_v_and_w::fm_read
subroutine, public fm_read(fm, bs_env, mat_name, idx)
...
Definition gw_compute_V_and_W.F:958
gw_compute_v_and_w::get_v_tr_r
subroutine, public get_v_tr_r(v_tr_r, pot_type, regularization_ri, bs_env, qs_env)
...
Definition gw_compute_V_and_W.F:1193
gw_compute_v_and_w::compute_v_k_by_lattice_sum
subroutine, public compute_v_k_by_lattice_sum(bs_env, qs_env, cfm_v_kp, ikp_batch)
...
Definition gw_compute_V_and_W.F:337
gw_compute_v_and_w::compute_w
subroutine, public compute_w(bs_env, qs_env, mat_chi_gamma_tau, fm_w_time)
Computes the screened Coulomb interaction on the imaginary-time grid, entirely in the RI auxiliary (P...
Definition gw_compute_V_and_W.F:115
gw_compute_v_and_w::apply_extrapol_factor
subroutine, public apply_extrapol_factor(bs_env, fm_w_mic_time)
...
Definition gw_compute_V_and_W.F:767
gw_utils_fm
Full-matrix operations not provided by the CP2K FM packages.
Definition gw_utils_fm.F:13
gw_utils_fm::fm_invert
subroutine, public fm_invert(matrix_a, eigenvalue_threshold, unit_nr)
Inverts a symmetric matrix. First, Cholesky decomposition is tried. If it fails, the matrix is diagon...
Definition gw_utils_fm.F:186
gw_utils_fm::fm_sqrt
subroutine, public fm_sqrt(matrix_a, matrix_b, eigenvalue_threshold, unit_nr)
For input A, computes B such that B^T B=A. First, Cholesky decomposition is tried....
Definition gw_utils_fm.F:232
gw_utils_fm::cfm_contract_aba
subroutine, public cfm_contract_aba(matrix_a, matrix_b, matrix_c, workspace)
Computes A^H B A for complex full matrices.
Definition gw_utils_fm.F:59
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::rtp_method_bse_linearized
integer, parameter, public rtp_method_bse_linearized
Definition input_constants.F:1013
input_constants::rtp_method_bse
integer, parameter, public rtp_method_bse
Definition input_constants.F:1013
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kpoint_coulomb_2c
Routines to compute the Coulomb integral V_(alpha beta)(k) for a k-point k using lattice summation in...
Definition kpoint_coulomb_2c.F:15
kpoint_coulomb_2c::build_2c_coulomb_matrix_kp
subroutine, public build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, atomic_kind_set, size_lattice_sum, operator_type, ikp_start, ikp_end)
...
Definition kpoint_coulomb_2c.F:76
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:144
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
mp2_ri_2c
Framework for 2c-integrals for RI.
Definition mp2_ri_2c.F:14
mp2_ri_2c::ri_2c_integral_mat_kp
subroutine, public ri_2c_integral_mat_kp(qs_env, cfm_matrix_m_kpoints, fm_matrix_struct, dimen_ri, ri_metric, kpoints, put_mat_ks_env, regularization_ri, ikp_ext, do_build_cell_index)
Computes the complex k-point RI metric as native CFM matrices.
Definition mp2_ri_2c.F:607
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:123
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
post_scf_bandstructure_types
Definition post_scf_bandstructure_types.F:13
post_scf_bandstructure_utils
Definition post_scf_bandstructure_utils.F:13
post_scf_bandstructure_utils::cfm_ikp_from_fm_gamma
subroutine, public cfm_ikp_from_fm_gamma(cfm_ikp, fm_gamma, ikp, qs_env, kpoints, basis_type)
Reconstruct a k-point operator from a real Gamma matrix using one MIC phase.
Definition post_scf_bandstructure_utils.F:3028
post_scf_bandstructure_utils::mic_contribution_from_ikp
subroutine, public mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, cfm_w_ikp_freq_j, ikp, kpoints, basis_type, wkp_ext)
...
Definition post_scf_bandstructure_utils.F:3069
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, 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, matrix_vhxc, 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, xcint_weights, 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.
Definition qs_environment_types.F:544
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition qs_neighbor_list_types.F:1096
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors.F:11
qs_tensors::build_2c_integrals
subroutine, public build_2c_integrals(t2c, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints, do_hfx_kpoints, ext_kpoints, regularization_ri)
...
Definition qs_tensors.F:3156
qs_tensors::build_2c_neighbor_lists
subroutine, public build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, sym_ij, molecular, dist_2d, pot_to_rad)
Build 2-center neighborlists adapted to different operators This mainly wraps build_neighbor_lists fo...
Definition qs_tensors.F:144
rpa_gw_kpoints_util
Routines treating GW and RPA calculations with kpoints.
Definition rpa_gw_kpoints_util.F:13
rpa_gw_kpoints_util::cfm_screened_interaction_from_eps_cholesky
subroutine, public cfm_screened_interaction_from_eps_cholesky(eps_chol, b, wc, n_active)
Form the screened interaction from a Cholesky-factorized dielectric matrix.
Definition rpa_gw_kpoints_util.F:106
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_types.F:62
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
particle_types::particle_type
Definition particle_types.F:35
post_scf_bandstructure_types::post_scf_bandstructure_type
Definition post_scf_bandstructure_types.F:65
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67