Line data Source code
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 : ! **************************************************************************************************
11 : MODULE 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_syrk,&
39 : cp_fm_triangular_multiply,&
40 : cp_fm_uplo_to_full
41 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
42 : USE cp_fm_diag, ONLY: cp_fm_power
43 : USE cp_fm_redistribution, ONLY: cp_fm_to_replicated_array
44 : USE cp_fm_types, ONLY: cp_fm_create,&
45 : cp_fm_get_info,&
46 : cp_fm_read_unformatted,&
47 : cp_fm_release,&
48 : cp_fm_set_all,&
49 : cp_fm_to_fm,&
50 : cp_fm_type,&
51 : cp_fm_write_unformatted
52 : USE cp_log_handling, ONLY: cp_get_default_logger,&
53 : cp_logger_type
54 : USE cp_output_handling, ONLY: cp_p_file,&
55 : cp_print_key_should_output,&
56 : cp_print_key_unit_nr
57 : USE distribution_2d_types, ONLY: distribution_2d_type
58 : USE gw_utils_fm, ONLY: cfm_contract_ABA,&
59 : fm_contract_ABA,&
60 : fm_sqrt
61 : USE input_constants, ONLY: rtp_method_bse,&
62 : rtp_method_bse_linearized
63 : USE input_section_types, ONLY: section_vals_type
64 : USE kinds, ONLY: default_path_length,&
65 : dp
66 : USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
67 : USE libint_2c_3c, ONLY: libint_potential_type
68 : USE machine, ONLY: m_walltime
69 : USE mathconstants, ONLY: z_one,&
70 : z_zero
71 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat_kp
72 : USE parallel_gemm_api, ONLY: parallel_gemm
73 : USE particle_methods, ONLY: get_particle_set
74 : USE particle_types, ONLY: particle_type
75 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
76 : USE post_scf_bandstructure_utils, ONLY: MIC_contribution_from_ikp,&
77 : cfm_ikp_from_fm_Gamma
78 : USE qs_environment_types, ONLY: get_qs_env,&
79 : qs_environment_type
80 : USE qs_kind_types, ONLY: qs_kind_type
81 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
82 : release_neighbor_list_sets
83 : USE qs_tensors, ONLY: build_2c_integrals,&
84 : build_2c_neighbor_lists
85 : USE rpa_gw_kpoints_util, ONLY: cfm_screened_interaction_from_eps_cholesky
86 : #include "./base/base_uses.f90"
87 :
88 : IMPLICIT NONE
89 : PRIVATE
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_compute_V_and_W'
91 : PUBLIC :: compute_W, &
92 : get_W_MIC, &
93 : compute_V_k_by_lattice_sum, &
94 : compute_MinvVsqrt_Vsqrt, &
95 : compute_fm_W_MIC_freq_j_from_chi, &
96 : apply_extrapol_factor, &
97 : fm_read, &
98 : fm_write, &
99 : get_V_tr_R
100 :
101 : CONTAINS
102 :
103 : ! **************************************************************************************************
104 : !> \brief Computes the screened Coulomb interaction on the imaginary-time grid, entirely in the
105 : !> RI auxiliary (PQ) basis:
106 : !> χ_PQ(iω) = Σ_τ w(ω,τ) cos(ωτ) χ_PQ(iτ) (cosine transform)
107 : !> V = B^T B, A = M^-1 B^T (Coulomb and RI metric factors)
108 : !> ε(iω) = Id + A^T χ(iω) A (dielectric function)
109 : !> W(iω) = A ( ε^-1(iω) - Id ) A^T (metric-dressed correlation part)
110 : !> W(iτ) = Σ_ω w̃(τ,ω) cos(ωτ) W(iω) (back transform)
111 : !> The RI metric is folded in before the back transform. V is the bare Coulomb matrix.
112 : !> \param bs_env ...
113 : !> \param qs_env ...
114 : !> \param mat_chi_Gamma_tau ...
115 : !> \param fm_W_time ...
116 : ! **************************************************************************************************
117 160 : SUBROUTINE compute_W(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_time)
118 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
119 : TYPE(qs_environment_type), POINTER :: qs_env
120 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
121 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time
122 :
123 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_W'
124 :
125 : INTEGER :: handle, i_t, j_w
126 : REAL(KIND=dp) :: t1
127 : TYPE(cp_fm_type) :: fm_coulomb_factor, &
128 : fm_Minv_coulomb_factor, fm_Minv_V_Minv
129 :
130 40 : CALL timeset(routineN, handle)
131 :
132 40 : t1 = m_walltime()
133 :
134 40 : CALL create_fm_W_MIC_time(bs_env, fm_W_time)
135 :
136 : ! 1. Allocate V and M matrices
137 40 : CALL cp_fm_create(fm_Minv_V_Minv, bs_env%fm_RI_RI%matrix_struct)
138 40 : CALL cp_fm_create(fm_coulomb_factor, bs_env%fm_RI_RI%matrix_struct)
139 40 : CALL cp_fm_create(fm_Minv_coulomb_factor, bs_env%fm_RI_RI%matrix_struct)
140 :
141 : ! Compute V = B^T B and A = M^-1 B^T.
142 40 : CALL compute_coulomb_factors(bs_env, qs_env, fm_Minv_V_Minv, fm_coulomb_factor, fm_Minv_coulomb_factor)
143 :
144 : ! Reuse the Coulomb buffer for M^-1 V M^-1 = A A^T.
145 40 : CALL cp_fm_syrk('U', 'N', bs_env%n_RI, 1.0_dp, fm_Minv_coulomb_factor, 1, 1, 0.0_dp, fm_Minv_V_Minv)
146 : ! The Coulomb factor is no longer needed; reuse it as workspace to complete the lower triangle.
147 40 : CALL cp_fm_uplo_to_full(fm_Minv_V_Minv, fm_coulomb_factor)
148 40 : CALL cp_fm_release(fm_coulomb_factor)
149 :
150 : ! 2. Loop over frequencies
151 590 : DO j_w = 1, bs_env%num_time_freq_points
152 : ! Fourier transformation of χ_PQ(iτ) to χ_PQ(iω_j)
153 550 : CALL compute_fm_chi_Gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
154 :
155 : ! ε(iω_j) = Id + A^T χ(iω_j) A; W(iω_j) = A (ε^-1(iω_j)-Id) A^T.
156 : CALL compute_fm_W_freq(bs_env, bs_env%fm_chi_Gamma_freq, fm_Minv_V_Minv, &
157 550 : fm_Minv_coulomb_factor, bs_env%fm_W_MIC_freq)
158 :
159 : ! Fourier transform from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
160 590 : CALL Fourier_transform_w_to_t(bs_env, fm_W_time, bs_env%fm_W_MIC_freq, j_w)
161 : END DO
162 :
163 40 : IF (bs_env%unit_nr > 0) THEN
164 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
165 20 : 'Computed W(iτ),', ' Execution time', m_walltime() - t1, ' s'
166 : END IF
167 :
168 40 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
169 :
170 : ! Cleanup
171 40 : CALL cp_fm_release(fm_Minv_V_Minv)
172 40 : CALL cp_fm_release(fm_Minv_coulomb_factor)
173 :
174 : ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0
175 40 : IF (bs_env%rtp_method == rtp_method_bse) THEN
176 14 : t1 = m_walltime()
177 14 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
178 : ! Set to zero
179 14 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
180 : ! Sum over all times
181 294 : DO i_t = 1, bs_env%num_time_freq_points
182 : ! Add the relevant structure with correct weight
183 : CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
184 294 : bs_env%time_frequency_grid%time_weights_at_zero_frequency(i_t), fm_W_time(i_t))
185 : END DO
186 : ! Done, save to file
187 14 : CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
188 : ! Report calculation
189 14 : IF (bs_env%unit_nr > 0) THEN
190 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
191 7 : 'Computed W(0),', ' Execution time', m_walltime() - t1, ' s'
192 : END IF
193 : END IF
194 :
195 40 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
196 :
197 40 : CALL timestop(handle)
198 :
199 40 : END SUBROUTINE compute_W
200 :
201 : ! **************************************************************************************************
202 : !> \brief ...
203 : !> \param bs_env ...
204 : !> \param qs_env ...
205 : !> \param mat_chi_Gamma_tau ...
206 : !> \param fm_W_MIC_time ...
207 : ! **************************************************************************************************
208 74 : SUBROUTINE get_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
209 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
210 : TYPE(qs_environment_type), POINTER :: qs_env
211 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
212 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
213 :
214 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_W_MIC'
215 :
216 : INTEGER :: handle
217 :
218 74 : CALL timeset(routineN, handle)
219 :
220 74 : IF (bs_env%all_W_exist) THEN
221 0 : CALL read_W_MIC_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
222 : ELSE
223 74 : CALL compute_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
224 : END IF
225 :
226 74 : CALL timestop(handle)
227 :
228 74 : END SUBROUTINE get_W_MIC
229 :
230 : ! **************************************************************************************************
231 : !> \brief ...
232 : !> \param bs_env ...
233 : !> \param qs_env ...
234 : !> \param mat_chi_Gamma_tau ...
235 : !> \param fm_W_MIC_time ...
236 : ! **************************************************************************************************
237 74 : SUBROUTINE compute_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
238 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
239 : TYPE(qs_environment_type), POINTER :: qs_env
240 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
241 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
242 :
243 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_W_MIC'
244 :
245 : INTEGER :: handle, i_t, ikp, ikp_batch, &
246 : ikp_in_batch, j_w
247 : REAL(KIND=dp) :: t1
248 : TYPE(cp_cfm_type) :: cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp
249 74 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_V_kp
250 :
251 74 : CALL timeset(routineN, handle)
252 :
253 74 : CALL create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
254 :
255 208 : DO ikp_batch = 1, bs_env%num_chi_eps_W_batches
256 :
257 134 : t1 = m_walltime()
258 :
259 : ! Compute V_PQ(k) = sum_R e^(ikR) <phi_P, cell 0 | 1/r | phi_Q, cell R>
260 134 : CALL compute_V_k_by_lattice_sum(bs_env, qs_env, cfm_V_kp, ikp_batch)
261 :
262 670 : DO ikp_in_batch = 1, bs_env%nkp_chi_eps_W_batch
263 :
264 536 : ikp = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + ikp_in_batch
265 :
266 536 : IF (ikp > bs_env%nkp_chi_eps_W_orig_plus_extra) THEN
267 : CYCLE
268 : END IF
269 :
270 : CALL compute_MinvVsqrt_Vsqrt(bs_env, qs_env, cfm_V_kp, &
271 324 : cfm_V_sqrt_ikp, cfm_M_inv_V_sqrt_ikp, ikp)
272 :
273 324 : CALL bs_env%para_env%sync()
274 324 : CALL cp_cfm_release(cfm_V_kp(ikp))
275 :
276 3922 : DO j_w = 1, bs_env%num_time_freq_points
277 :
278 : ! check if we need this (ikp, ω_j) combination for approximate k-point extrapolation
279 3464 : IF (bs_env%approx_kp_extrapol .AND. j_w > 1 .AND. &
280 : ikp > bs_env%nkp_chi_eps_W_orig) CYCLE
281 :
282 : CALL compute_fm_W_MIC_freq_j(bs_env, qs_env, bs_env%fm_W_MIC_freq, j_w, ikp, &
283 : mat_chi_Gamma_tau, cfm_M_inv_V_sqrt_ikp, &
284 3140 : cfm_V_sqrt_ikp)
285 :
286 : ! Fourier trafo from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
287 4000 : CALL Fourier_transform_w_to_t(bs_env, fm_W_MIC_time, bs_env%fm_W_MIC_freq, j_w)
288 :
289 : END DO ! ω_j
290 :
291 : END DO ! ikp_in_batch
292 :
293 134 : DEALLOCATE (cfm_V_kp)
294 :
295 208 : IF (bs_env%unit_nr > 0) THEN
296 : WRITE (bs_env%unit_nr, '(T2,A,I12,A,I3,A,F10.1,A)') &
297 67 : 'Computed W(iτ,k) for k-point batch', &
298 67 : ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
299 134 : ', Execution time', m_walltime() - t1, ' s'
300 : END IF
301 :
302 : END DO ! ikp_batch
303 :
304 74 : IF (bs_env%approx_kp_extrapol) THEN
305 2 : CALL apply_extrapol_factor(bs_env, fm_W_MIC_time)
306 : END IF
307 :
308 : ! M^-1(k=0) W^MIC(iτ) M^-1(k=0) -> fm_W_MIC_time
309 74 : CALL fm_contract_ABA(bs_env%fm_Minv_Gamma, fm_W_MIC_time)
310 :
311 1438 : DO i_t = 1, bs_env%num_time_freq_points
312 1438 : CALL fm_write(fm_W_MIC_time(i_t), i_t, bs_env%W_time_name, qs_env)
313 : END DO
314 :
315 74 : CALL cp_cfm_release(cfm_M_inv_V_sqrt_ikp)
316 74 : CALL cp_cfm_release(cfm_V_sqrt_ikp)
317 74 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
318 :
319 : ! W^MIC(iω=0) from W^MIC(iτ) for the post-GW consumers; the file is the RT-BSE restart's
320 74 : IF (bs_env%rtp_method == rtp_method_bse .OR. &
321 : bs_env%rtp_method == rtp_method_bse_linearized) THEN
322 58 : CALL sum_W_MIC_freq_zero(bs_env, fm_W_MIC_time)
323 58 : CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
324 16 : ELSE IF (bs_env%keep_W_freq_zero) THEN
325 6 : CALL sum_W_MIC_freq_zero(bs_env, fm_W_MIC_time)
326 : END IF
327 :
328 74 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
329 :
330 74 : CALL timestop(handle)
331 :
332 148 : END SUBROUTINE compute_W_MIC
333 :
334 : ! **************************************************************************************************
335 : !> \brief ...
336 : !> \param bs_env ...
337 : !> \param qs_env ...
338 : !> \param cfm_V_kp ...
339 : !> \param ikp_batch ...
340 : ! **************************************************************************************************
341 188 : SUBROUTINE compute_V_k_by_lattice_sum(bs_env, qs_env, cfm_V_kp, ikp_batch)
342 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
343 : TYPE(qs_environment_type), POINTER :: qs_env
344 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_V_kp
345 : INTEGER :: ikp_batch
346 :
347 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_k_by_lattice_sum'
348 :
349 : INTEGER :: handle, ikp, ikp_end, ikp_start, &
350 : nkp_chi_eps_W_batch, re_im
351 188 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
352 : TYPE(cell_type), POINTER :: cell
353 : TYPE(cp_fm_type) :: fm_V_kp_work
354 188 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_kp
355 188 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
356 188 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
357 :
358 188 : CALL timeset(routineN, handle)
359 :
360 188 : nkp_chi_eps_W_batch = bs_env%nkp_chi_eps_W_batch
361 :
362 188 : ikp_start = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + 1
363 188 : ikp_end = MIN(ikp_batch*bs_env%nkp_chi_eps_W_batch, bs_env%kpoints_chi_eps_W%nkp)
364 :
365 188 : NULLIFY (mat_V_kp)
366 2188 : ALLOCATE (mat_V_kp(ikp_start:ikp_end, 2))
367 :
368 564 : DO re_im = 1, 2
369 1624 : DO ikp = ikp_start, ikp_end
370 1060 : NULLIFY (mat_V_kp(ikp, re_im)%matrix)
371 1060 : ALLOCATE (mat_V_kp(ikp, re_im)%matrix)
372 1060 : CALL dbcsr_create(mat_V_kp(ikp, re_im)%matrix, template=bs_env%mat_RI_RI%matrix)
373 1060 : CALL dbcsr_reserve_all_blocks(mat_V_kp(ikp, re_im)%matrix)
374 1436 : CALL dbcsr_set(mat_V_kp(ikp, re_im)%matrix, 0.0_dp)
375 : END DO ! ikp
376 : END DO ! re_im
377 :
378 : CALL get_qs_env(qs_env=qs_env, &
379 : particle_set=particle_set, &
380 : cell=cell, &
381 : qs_kind_set=qs_kind_set, &
382 188 : atomic_kind_set=atomic_kind_set)
383 :
384 188 : IF (ikp_end <= bs_env%nkp_chi_eps_W_orig) THEN
385 :
386 : ! 1. 2c Coulomb integrals for the first "original" k-point grid
387 384 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
388 :
389 92 : ELSE IF (ikp_start > bs_env%nkp_chi_eps_W_orig .AND. &
390 : ikp_end <= bs_env%nkp_chi_eps_W_orig_plus_extra) THEN
391 :
392 : ! 2. 2c Coulomb integrals for the second "extrapolation" k-point grid
393 368 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_extra
394 :
395 : ELSE
396 :
397 0 : CPABORT("Error with k-point parallelization.")
398 :
399 : END IF
400 :
401 : CALL build_2c_coulomb_matrix_kp(mat_V_kp, &
402 : bs_env%kpoints_chi_eps_W, &
403 : basis_type="RI_AUX", &
404 : cell=cell, &
405 : particle_set=particle_set, &
406 : qs_kind_set=qs_kind_set, &
407 : atomic_kind_set=atomic_kind_set, &
408 : size_lattice_sum=bs_env%size_lattice_sum_V, &
409 : operator_type=operator_coulomb, &
410 : ikp_start=ikp_start, &
411 188 : ikp_end=ikp_end)
412 :
413 752 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
414 :
415 1094 : ALLOCATE (cfm_V_kp(ikp_start:ikp_end))
416 188 : CALL cp_fm_create(fm_V_kp_work, bs_env%fm_RI_RI%matrix_struct)
417 718 : DO ikp = ikp_start, ikp_end
418 530 : CALL cp_cfm_create(cfm_V_kp(ikp), bs_env%fm_RI_RI%matrix_struct)
419 530 : CALL cp_fm_set_all(fm_V_kp_work, 0.0_dp)
420 530 : CALL copy_dbcsr_to_fm(mat_V_kp(ikp, 1)%matrix, fm_V_kp_work)
421 89500 : cfm_V_kp(ikp)%local_data = CMPLX(fm_V_kp_work%local_data, 0.0_dp, KIND=dp)
422 :
423 530 : CALL cp_fm_set_all(fm_V_kp_work, 0.0_dp)
424 530 : CALL copy_dbcsr_to_fm(mat_V_kp(ikp, 2)%matrix, fm_V_kp_work)
425 : cfm_V_kp(ikp)%local_data = CMPLX(REAL(cfm_V_kp(ikp)%local_data, KIND=dp), &
426 89688 : fm_V_kp_work%local_data, KIND=dp)
427 : END DO
428 188 : CALL cp_fm_release(fm_V_kp_work)
429 188 : CALL dbcsr_deallocate_matrix_set(mat_V_kp)
430 :
431 188 : CALL timestop(handle)
432 :
433 376 : END SUBROUTINE compute_V_k_by_lattice_sum
434 :
435 : ! **************************************************************************************************
436 : !> \brief ...
437 : !> \param bs_env ...
438 : !> \param qs_env ...
439 : !> \param cfm_V_kp ...
440 : !> \param cfm_V_sqrt_ikp ...
441 : !> \param cfm_M_inv_V_sqrt_ikp ...
442 : !> \param ikp ...
443 : ! **************************************************************************************************
444 530 : SUBROUTINE compute_MinvVsqrt_Vsqrt(bs_env, qs_env, cfm_V_kp, cfm_V_sqrt_ikp, &
445 : cfm_M_inv_V_sqrt_ikp, ikp)
446 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
447 : TYPE(qs_environment_type), POINTER :: qs_env
448 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_V_kp
449 : TYPE(cp_cfm_type) :: cfm_V_sqrt_ikp, cfm_M_inv_V_sqrt_ikp
450 : INTEGER :: ikp
451 :
452 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_MinvVsqrt_Vsqrt'
453 :
454 : INTEGER :: handle, info, n_RI
455 : TYPE(cp_cfm_type) :: cfm_M_inv_ikp, cfm_work
456 530 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_M_ikp
457 :
458 530 : CALL timeset(routineN, handle)
459 :
460 530 : n_RI = bs_env%n_RI
461 :
462 : ! get here M(k) and write it to cfm_M_ikp
463 : CALL RI_2c_integral_mat_kp(qs_env, cfm_M_ikp, cfm_V_kp(ikp)%matrix_struct, n_RI, bs_env%ri_metric, &
464 : bs_env%kpoints_chi_eps_W, regularization_RI=bs_env%regularization_RI, &
465 530 : ikp_ext=ikp, do_build_cell_index=(ikp == 1))
466 :
467 530 : IF (ikp == 1) THEN
468 82 : CALL cp_cfm_create(cfm_V_sqrt_ikp, cfm_V_kp(ikp)%matrix_struct)
469 82 : CALL cp_cfm_create(cfm_M_inv_V_sqrt_ikp, cfm_V_kp(ikp)%matrix_struct)
470 : END IF
471 530 : CALL cp_cfm_create(cfm_M_inv_ikp, cfm_V_kp(ikp)%matrix_struct)
472 :
473 530 : CALL cp_cfm_to_cfm(cfm_M_ikp(1), cfm_M_inv_ikp)
474 530 : CALL cp_cfm_to_cfm(cfm_V_kp(ikp), cfm_V_sqrt_ikp)
475 :
476 530 : CALL cp_cfm_release(cfm_M_ikp(1))
477 530 : DEALLOCATE (cfm_M_ikp)
478 :
479 530 : CALL cp_cfm_create(cfm_work, cfm_V_kp(ikp)%matrix_struct)
480 :
481 : ! M(k) -> M^-1(k)
482 530 : CALL cp_cfm_to_cfm(cfm_M_inv_ikp, cfm_work)
483 530 : CALL cp_cfm_cholesky_decompose(matrix=cfm_M_inv_ikp, n=n_RI, info_out=info)
484 530 : IF (info == 0) THEN
485 : ! successful Cholesky decomposition
486 530 : CALL cp_cfm_cholesky_invert(cfm_M_inv_ikp)
487 : ! symmetrize the result
488 530 : CALL cp_cfm_uplo_to_full(cfm_M_inv_ikp)
489 : ELSE
490 : ! Cholesky decomposition not successful: use expensive diagonalization
491 0 : CALL cp_cfm_power(cfm_work, threshold=bs_env%eps_eigval_mat_RI, exponent=-1.0_dp)
492 0 : CALL cp_cfm_to_cfm(cfm_work, cfm_M_inv_ikp)
493 : END IF
494 :
495 : ! V(k) -> L(k) with L^H(k)*L(k) = V(k) [L(k) can be just considered to be V^0.5(k)]
496 530 : CALL cp_cfm_to_cfm(cfm_V_sqrt_ikp, cfm_work)
497 530 : CALL cp_cfm_cholesky_decompose(matrix=cfm_V_sqrt_ikp, n=n_RI, info_out=info)
498 530 : IF (info == 0) THEN
499 : ! successful Cholesky decomposition
500 530 : CALL clean_lower_part(cfm_V_sqrt_ikp)
501 : ELSE
502 : ! Cholesky decomposition not successful: use expensive diagonalization
503 0 : CALL cp_cfm_power(cfm_work, threshold=0.0_dp, exponent=0.5_dp)
504 0 : CALL cp_cfm_to_cfm(cfm_work, cfm_V_sqrt_ikp)
505 : END IF
506 530 : CALL cp_cfm_release(cfm_work)
507 :
508 : ! get M^-1(k)*V^0.5(k)
509 : CALL parallel_gemm("N", "C", n_RI, n_RI, n_RI, z_one, cfm_M_inv_ikp, cfm_V_sqrt_ikp, &
510 530 : z_zero, cfm_M_inv_V_sqrt_ikp)
511 :
512 530 : CALL cp_cfm_release(cfm_M_inv_ikp)
513 :
514 530 : CALL timestop(handle)
515 :
516 1060 : END SUBROUTINE compute_MinvVsqrt_Vsqrt
517 :
518 : ! **************************************************************************************************
519 : !> \brief Compute the screened interaction for one frequency and k-point.
520 : !> \param bs_env band-structure environment
521 : !> \param cfm_chi_eps_W_ikp_freq_j input chi; overwritten by W
522 : !> \param cfm_V_sqrt_ikp square root of the Coulomb matrix
523 : !> \param cfm_M_inv_V_sqrt_ikp metric-weighted Coulomb factor
524 : ! **************************************************************************************************
525 6692 : SUBROUTINE compute_cfm_W_ikp_freq_j(bs_env, cfm_chi_eps_W_ikp_freq_j, cfm_V_sqrt_ikp, &
526 : cfm_M_inv_V_sqrt_ikp)
527 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
528 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_chi_eps_W_ikp_freq_j
529 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_V_sqrt_ikp, cfm_M_inv_V_sqrt_ikp
530 :
531 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cfm_W_ikp_freq_j'
532 :
533 : INTEGER :: handle, info
534 :
535 3346 : CALL timeset(routineN, handle)
536 :
537 : ! 1. epsilon(iw,k) = Id - V^0.5(k)*M^-1(k)*chi(iw,k)*M^-1(k)*V^0.5(k)
538 :
539 : ! 1. a) Form V^0.5(k)*M^-1(k)*chi(iw,k)*M^-1(k)*V^0.5(k).
540 3346 : CALL cfm_contract_ABA(cfm_M_inv_V_sqrt_ikp, cfm_chi_eps_W_ikp_freq_j)
541 :
542 : ! 1. b) Add the identity matrix.
543 3346 : CALL cp_cfm_add_on_diag(cfm_chi_eps_W_ikp_freq_j, z_one)
544 :
545 : ! 2. W(iw,k) = V^0.5(k)*(epsilon^-1(iw,k)-Id)*V^0.5(k)
546 :
547 : ! 2. a) Cholesky-factorize epsilon before inversion.
548 3346 : CALL cp_cfm_cholesky_decompose(matrix=cfm_chi_eps_W_ikp_freq_j, n=bs_env%n_RI, info_out=info)
549 3346 : CPASSERT(info == 0)
550 :
551 : ! 2. b) Invert epsilon, subtract Id, then contract with V^0.5 on both sides.
552 3346 : CALL cfm_screened_interaction_from_eps_cholesky(cfm_chi_eps_W_ikp_freq_j, cfm_V_sqrt_ikp)
553 :
554 3346 : CALL timestop(handle)
555 :
556 3346 : END SUBROUTINE compute_cfm_W_ikp_freq_j
557 :
558 : ! **************************************************************************************************
559 : !> \brief ...
560 : !> \param bs_env ...
561 : !> \param mat_chi_Gamma_tau ...
562 : !> \param fm_W_MIC_time ...
563 : ! **************************************************************************************************
564 0 : SUBROUTINE read_W_MIC_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
565 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
566 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
567 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
568 :
569 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_W_MIC_time'
570 :
571 : INTEGER :: handle, i_t
572 : REAL(KIND=dp) :: t1
573 :
574 0 : CALL timeset(routineN, handle)
575 :
576 0 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
577 0 : CALL create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
578 :
579 0 : DO i_t = 1, bs_env%num_time_freq_points
580 :
581 0 : t1 = m_walltime()
582 :
583 0 : CALL fm_read(fm_W_MIC_time(i_t), bs_env, bs_env%W_time_name, i_t)
584 :
585 0 : IF (bs_env%unit_nr > 0) THEN
586 : WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
587 0 : 'Read W^MIC(iτ) from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
588 0 : ', Execution time', m_walltime() - t1, ' s'
589 : END IF
590 :
591 : END DO
592 :
593 0 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
594 :
595 : ! W^MIC(iω=0) for the post-GW consumers, built here because fm_W_MIC_time does not outlive
596 : ! the W stage: RT-BSE reads its own restart file, PROPERTIES%BSE sums the W(iτ) just read.
597 0 : IF (bs_env%rtp_method == rtp_method_bse .OR. &
598 : bs_env%rtp_method == rtp_method_bse_linearized) THEN
599 0 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
600 0 : t1 = m_walltime()
601 0 : CALL fm_read(bs_env%fm_W_MIC_freq_zero, bs_env, "W_freq_rtp", 0)
602 0 : IF (bs_env%unit_nr > 0) THEN
603 : WRITE (bs_env%unit_nr, '(T2,A,I3,A,I3,A,F10.1,A)') &
604 0 : 'Read W^MIC(f=0) from file for freq. point ', 1, ' /', 1, &
605 0 : ', Execution time', m_walltime() - t1, ' s'
606 : END IF
607 0 : ELSE IF (bs_env%keep_W_freq_zero) THEN
608 : ! no RTP file in a LRBSE-only run: W(iτ) was just read from the W_time files
609 0 : CALL sum_W_MIC_freq_zero(bs_env, fm_W_MIC_time)
610 : END IF
611 :
612 0 : CALL timestop(handle)
613 :
614 0 : END SUBROUTINE read_W_MIC_time
615 :
616 : ! **************************************************************************************************
617 : !> \brief ...
618 : !> \param bs_env ...
619 : !> \param qs_env ...
620 : !> \param fm_W_MIC_freq_j ...
621 : !> \param j_w ...
622 : !> \param ikp ...
623 : !> \param mat_chi_Gamma_tau ...
624 : !> \param cfm_M_inv_V_sqrt_ikp ...
625 : !> \param cfm_V_sqrt_ikp ...
626 : ! **************************************************************************************************
627 3140 : SUBROUTINE compute_fm_W_MIC_freq_j(bs_env, qs_env, fm_W_MIC_freq_j, j_w, ikp, mat_chi_Gamma_tau, &
628 : cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
629 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
630 : TYPE(qs_environment_type), POINTER :: qs_env
631 : TYPE(cp_fm_type) :: fm_W_MIC_freq_j
632 : INTEGER :: j_w, ikp
633 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
634 : TYPE(cp_cfm_type) :: cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp
635 :
636 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_MIC_freq_j'
637 :
638 : INTEGER :: handle
639 :
640 3140 : CALL timeset(routineN, handle)
641 :
642 : ! 1. Fourier transformation of χ_PQ(iτ,k=0) to χ_PQ(iω_j,k=0)
643 3140 : CALL compute_fm_chi_Gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
644 :
645 : ! 2. χ_PQ(iω_j,k=0) -> W_PQ^MIC(iω_j)
646 : CALL compute_fm_W_MIC_freq_j_from_chi(bs_env, qs_env, fm_W_MIC_freq_j, j_w, ikp, &
647 3140 : cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
648 :
649 3140 : CALL timestop(handle)
650 :
651 3140 : END SUBROUTINE compute_fm_W_MIC_freq_j
652 :
653 : ! **************************************************************************************************
654 : !> \brief W_PQ^MIC(iω_j) from χ_PQ(iω_j,k=0) in bs_env%fm_chi_Gamma_freq
655 : !> \param bs_env ...
656 : !> \param qs_env ...
657 : !> \param fm_W_MIC_freq_j ...
658 : !> \param j_w ...
659 : !> \param ikp ...
660 : !> \param cfm_M_inv_V_sqrt_ikp ...
661 : !> \param cfm_V_sqrt_ikp ...
662 : ! **************************************************************************************************
663 3346 : SUBROUTINE compute_fm_W_MIC_freq_j_from_chi(bs_env, qs_env, fm_W_MIC_freq_j, j_w, ikp, &
664 : cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
665 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
666 : TYPE(qs_environment_type), POINTER :: qs_env
667 : TYPE(cp_fm_type) :: fm_W_MIC_freq_j
668 : INTEGER :: j_w, ikp
669 : TYPE(cp_cfm_type) :: cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp
670 :
671 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_MIC_freq_j_from_chi'
672 :
673 : INTEGER :: handle
674 : TYPE(cp_cfm_type) :: cfm_chi_eps_W_ikp_freq_j
675 :
676 3346 : CALL timeset(routineN, handle)
677 :
678 3346 : CALL cp_fm_set_all(fm_W_MIC_freq_j, 0.0_dp)
679 :
680 : ! 1. Get χ_PQ(iω_j,k_i) from χ_PQ(iω_j,k=0) using the minimum image convention
681 : CALL cfm_ikp_from_fm_Gamma(cfm_chi_eps_W_ikp_freq_j, bs_env%fm_chi_Gamma_freq, &
682 3346 : ikp, qs_env, bs_env%kpoints_chi_eps_W, "RI_AUX")
683 :
684 : ! 2. Remove all negative eigenvalues from χ_PQ(iω_j,k_i)
685 3346 : CALL cp_cfm_power(cfm_chi_eps_W_ikp_freq_j, threshold=0.0_dp, exponent=1.0_dp)
686 :
687 : ! 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)
688 : CALL compute_cfm_W_ikp_freq_j(bs_env, cfm_chi_eps_W_ikp_freq_j, cfm_V_sqrt_ikp, &
689 3346 : cfm_M_inv_V_sqrt_ikp)
690 :
691 : ! 4. k-point integration W_PQ(iω_j, k_i) to W_PQ^MIC(iω_j)
692 3346 : SELECT CASE (bs_env%approx_kp_extrapol)
693 : CASE (.FALSE.)
694 : ! default: standard k-point extrapolation
695 : CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, &
696 : cfm_chi_eps_W_ikp_freq_j, ikp, &
697 3346 : bs_env%kpoints_chi_eps_W, "RI_AUX")
698 : CASE (.TRUE.)
699 : ! for approximate kpoint extrapolation: get W_PQ^MIC(iω_1) with and without k-point
700 : ! extrapolation to compute the extrapolation factor f_PQ for every PQ-matrix element,
701 : ! f_PQ = (W_PQ^MIC(iω_1) with extrapolation) / (W_PQ^MIC(iω_1) without extrapolation)
702 :
703 : ! for ω_1, we compute the k-point extrapolated result using all k-points
704 196 : IF (j_w == 1) THEN
705 :
706 : ! k-point extrapolated
707 : CALL MIC_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_extra, &
708 : cfm_chi_eps_W_ikp_freq_j, ikp, &
709 : bs_env%kpoints_chi_eps_W, &
710 52 : "RI_AUX")
711 : ! non-kpoint extrapolated
712 52 : IF (ikp <= bs_env%nkp_chi_eps_W_orig) THEN
713 : CALL MIC_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_no_extra, &
714 : cfm_chi_eps_W_ikp_freq_j, ikp, &
715 : bs_env%kpoints_chi_eps_W, &
716 16 : "RI_AUX", wkp_ext=bs_env%wkp_orig)
717 : END IF
718 :
719 : END IF
720 :
721 : ! for all ω_j, we need to compute W^MIC without k-point extrpolation
722 196 : IF (ikp <= bs_env%nkp_chi_eps_W_orig) THEN
723 : CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, &
724 : cfm_chi_eps_W_ikp_freq_j, &
725 : ikp, bs_env%kpoints_chi_eps_W, "RI_AUX", &
726 160 : wkp_ext=bs_env%wkp_orig)
727 : END IF
728 : END SELECT
729 :
730 3346 : CALL cp_cfm_release(cfm_chi_eps_W_ikp_freq_j)
731 :
732 3346 : CALL timestop(handle)
733 :
734 3346 : END SUBROUTINE compute_fm_W_MIC_freq_j_from_chi
735 :
736 : ! **************************************************************************************************
737 : !> \brief ...
738 : !> \param cfm_mat ...
739 : ! **************************************************************************************************
740 1060 : SUBROUTINE clean_lower_part(cfm_mat)
741 : TYPE(cp_cfm_type) :: cfm_mat
742 :
743 : CHARACTER(LEN=*), PARAMETER :: routineN = 'clean_lower_part'
744 :
745 : INTEGER :: handle, i_row, j_col, j_global, &
746 : ncol_local, nrow_local
747 530 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
748 :
749 530 : CALL timeset(routineN, handle)
750 :
751 : CALL cp_cfm_get_info(matrix=cfm_mat, &
752 : nrow_local=nrow_local, ncol_local=ncol_local, &
753 530 : row_indices=row_indices, col_indices=col_indices)
754 :
755 6646 : DO j_col = 1, ncol_local
756 6116 : j_global = col_indices(j_col)
757 89500 : DO i_row = 1, nrow_local
758 88970 : IF (j_global < row_indices(i_row)) cfm_mat%local_data(i_row, j_col) = z_zero
759 : END DO
760 : END DO
761 :
762 530 : CALL timestop(handle)
763 :
764 530 : END SUBROUTINE clean_lower_part
765 :
766 : ! **************************************************************************************************
767 : !> \brief ...
768 : !> \param bs_env ...
769 : !> \param fm_W_MIC_time ...
770 : ! **************************************************************************************************
771 4 : SUBROUTINE apply_extrapol_factor(bs_env, fm_W_MIC_time)
772 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
773 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
774 :
775 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_extrapol_factor'
776 : REAL(KIND=dp), PARAMETER :: eps_W_no_extra_1 = 1.0E-13_dp
777 :
778 : INTEGER :: handle, i, i_t, j, ncol_local, nrow_local
779 : REAL(KIND=dp) :: extrapol_factor, W_extra_1, W_no_extra_1
780 :
781 2 : CALL timeset(routineN, handle)
782 :
783 2 : CALL cp_fm_get_info(matrix=fm_W_MIC_time(1), nrow_local=nrow_local, ncol_local=ncol_local)
784 :
785 22 : DO i_t = 1, SIZE(fm_W_MIC_time)
786 122 : DO j = 1, ncol_local
787 370 : DO i = 1, nrow_local
788 :
789 250 : W_extra_1 = bs_env%fm_W_MIC_freq_1_extra%local_data(i, j)
790 250 : W_no_extra_1 = bs_env%fm_W_MIC_freq_1_no_extra%local_data(i, j)
791 :
792 250 : IF (ABS(W_no_extra_1) > eps_W_no_extra_1) THEN
793 190 : extrapol_factor = ABS(W_extra_1/W_no_extra_1)
794 : ELSE
795 : extrapol_factor = 1.0_dp
796 : END IF
797 :
798 : ! reset extrapolation factor if it is very large
799 190 : IF (extrapol_factor > 10.0_dp) extrapol_factor = 1.0_dp
800 :
801 : fm_W_MIC_time(i_t)%local_data(i, j) = fm_W_MIC_time(i_t)%local_data(i, j) &
802 350 : *extrapol_factor
803 : END DO
804 : END DO
805 : END DO
806 :
807 2 : CALL timestop(handle)
808 :
809 2 : END SUBROUTINE apply_extrapol_factor
810 :
811 : ! **************************************************************************************************
812 : !> \brief Computes V, its factor B with V = B^T B, and A = M^-1 B^T.
813 : !> The Coulomb matrix V is constructed by the RI-RS k-point path. The inverse metric
814 : !> M^-1(k=0) is precomputed once in gw_utils and stored in bs_env.
815 : !> \param bs_env ...
816 : !> \param qs_env ...
817 : !> \param fm_V Coulomb matrix V(k=0)
818 : !> \param fm_coulomb_factor upper Cholesky factor B, or symmetric square root if Cholesky fails
819 : !> \param fm_Minv_coulomb_factor product A = M^-1(k=0) B^T
820 : ! **************************************************************************************************
821 40 : SUBROUTINE compute_coulomb_factors(bs_env, qs_env, fm_V, fm_coulomb_factor, fm_Minv_coulomb_factor)
822 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
823 : TYPE(qs_environment_type), POINTER :: qs_env
824 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_V, fm_coulomb_factor, &
825 : fm_Minv_coulomb_factor
826 :
827 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_coulomb_factors'
828 :
829 : INTEGER :: handle
830 40 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
831 : TYPE(cell_type), POINTER :: cell
832 40 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_kp
833 40 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
834 40 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
835 :
836 40 : CALL timeset(routineN, handle)
837 :
838 40 : IF (bs_env%auto_ri%enabled) THEN
839 : ! -------------------------------------------------------------------
840 : ! 1a. The optimized AB functions span two atoms, so their previously transformed
841 : ! Coulomb matrix cannot be rebuilt by the atom-local k-point integral routine.
842 : ! -------------------------------------------------------------------
843 4 : CALL cp_fm_to_fm(bs_env%auto_ri%V_pq, fm_V)
844 : ELSE
845 : ! -------------------------------------------------------------------
846 : ! 1b. Build Coulomb Matrix V(k=0) using the kp-routine but only for ikp=1.
847 : ! -------------------------------------------------------------------
848 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, &
849 36 : qs_kind_set=qs_kind_set)
850 36 : atomic_kind_set => bs_env%ri_rs%atomic_kind_set
851 :
852 180 : ALLOCATE (mat_V_kp(1:1, 1:2))
853 36 : NULLIFY (mat_V_kp(1, 1)%matrix, mat_V_kp(1, 2)%matrix)
854 36 : ALLOCATE (mat_V_kp(1, 1)%matrix, mat_V_kp(1, 2)%matrix)
855 36 : CALL dbcsr_create(mat_V_kp(1, 1)%matrix, template=bs_env%mat_RI_RI%matrix)
856 36 : CALL dbcsr_reserve_all_blocks(mat_V_kp(1, 1)%matrix)
857 36 : CALL dbcsr_set(mat_V_kp(1, 1)%matrix, 0.0_dp)
858 36 : CALL dbcsr_create(mat_V_kp(1, 2)%matrix, template=bs_env%mat_RI_RI%matrix)
859 36 : CALL dbcsr_reserve_all_blocks(mat_V_kp(1, 2)%matrix)
860 : ! The dummy imaginary part is required only by the k-point routine interface.
861 36 : CALL dbcsr_set(mat_V_kp(1, 2)%matrix, 0.0_dp)
862 :
863 144 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
864 : CALL build_2c_coulomb_matrix_kp(mat_V_kp, bs_env%kpoints_chi_eps_W, "RI_AUX", cell, &
865 : particle_set, qs_kind_set, atomic_kind_set, &
866 36 : bs_env%size_lattice_sum_V, operator_coulomb, 1, 1)
867 36 : CALL copy_dbcsr_to_fm(mat_V_kp(1, 1)%matrix, fm_V)
868 :
869 36 : CALL dbcsr_deallocate_matrix(mat_V_kp(1, 1)%matrix)
870 36 : CALL dbcsr_deallocate_matrix(mat_V_kp(1, 2)%matrix)
871 36 : DEALLOCATE (mat_V_kp)
872 : END IF
873 :
874 : ! -----------------------------------------------------------------------
875 : ! 2. V -> B, with V = B^T B.
876 : ! -----------------------------------------------------------------------
877 40 : CALL fm_sqrt(fm_V, fm_coulomb_factor, bs_env%eps_eigval_mat_RI, bs_env%unit_nr)
878 :
879 : ! -----------------------------------------------------------------------
880 : ! 3. A = M^-1(k=0) B^T.
881 : ! -----------------------------------------------------------------------
882 : CALL parallel_gemm("N", "T", bs_env%n_RI, bs_env%n_RI, bs_env%n_RI, 1.0_dp, &
883 40 : bs_env%fm_Minv_Gamma, fm_coulomb_factor, 0.0_dp, fm_Minv_coulomb_factor)
884 :
885 40 : CALL timestop(handle)
886 :
887 40 : END SUBROUTINE compute_coulomb_factors
888 :
889 : ! **************************************************************************************************
890 : !> \brief Computes the metric-dressed screened interaction at one imaginary frequency.
891 : !> With V = B^T B and A = M^-1 B^T, form ε = Id + A^T χ A.
892 : !> For ε = U^T U, solve Z U = A, then W = Z Z^T - A A^T.
893 : !> If Cholesky fails, use the same eigenvalue cutoff as the inverse path:
894 : !> Z = A ε^(-1/2), with discarded eigenmodes set to zero.
895 : !> \param bs_env ...
896 : !> \param fm_chi_freq_j ...
897 : !> \param fm_Minv_V_Minv full symmetric matrix M^-1 V M^-1 = A A^T
898 : !> \param fm_Minv_coulomb_factor RI metric and Coulomb factor A = M^-1 B^T
899 : !> \param fm_W_freq_j full symmetric metric-dressed correlation interaction W
900 : ! **************************************************************************************************
901 2200 : SUBROUTINE compute_fm_W_freq(bs_env, fm_chi_freq_j, fm_Minv_V_Minv, fm_Minv_coulomb_factor, fm_W_freq_j)
902 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
903 : TYPE(cp_fm_type), INTENT(IN) :: fm_chi_freq_j, fm_Minv_V_Minv, &
904 : fm_Minv_coulomb_factor
905 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_W_freq_j
906 :
907 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_freq'
908 :
909 : INTEGER :: handle, info, n_dependent, n_RI
910 : TYPE(cp_fm_type) :: fm_eps_freq_j, fm_work
911 :
912 550 : CALL timeset(routineN, handle)
913 :
914 550 : n_RI = bs_env%n_RI
915 :
916 550 : CALL cp_fm_create(fm_eps_freq_j, fm_chi_freq_j%matrix_struct)
917 550 : CALL cp_fm_create(fm_work, fm_chi_freq_j%matrix_struct)
918 :
919 : ! -----------------------------------------------------------------------
920 : ! 1. ε(iω_j) = Id + A^T χ(iω_j) A, where A = M^-1 B^T and V = B^T B.
921 : ! -----------------------------------------------------------------------
922 : ! work = χ(iω_j) A
923 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
924 550 : fm_chi_freq_j, fm_Minv_coulomb_factor, 0.0_dp, fm_work)
925 :
926 : ! eps_work = A^T work
927 : CALL parallel_gemm('T', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
928 550 : fm_Minv_coulomb_factor, fm_work, 0.0_dp, fm_eps_freq_j)
929 :
930 : ! Add the identity without changing the sign convention of χ.
931 550 : CALL cp_fm_add_on_diag(fm_eps_freq_j, 1.0_dp)
932 :
933 : ! Force perfect symmetry before Cholesky to avoid info != 0 due to GEMM noise
934 550 : CALL cp_fm_uplo_to_full(fm_eps_freq_j, fm_work)
935 :
936 : ! -----------------------------------------------------------------------
937 : ! 2. W(iω_j) = A (ε^-1(iω_j) - Id) A^T = Z Z^T - A A^T.
938 : ! -----------------------------------------------------------------------
939 :
940 : ! Preserve ε because a failed Cholesky decomposition overwrites it.
941 550 : CALL cp_fm_to_fm(fm_eps_freq_j, fm_work)
942 550 : CALL cp_fm_cholesky_decompose(fm_eps_freq_j, info_out=info)
943 550 : IF (info == 0) THEN
944 : ! Z = A U^-1, so Z Z^T = A ε^-1 A^T. No explicit inverse is needed.
945 550 : CALL cp_fm_to_fm(fm_Minv_coulomb_factor, fm_work)
946 550 : CALL cp_fm_triangular_multiply(fm_eps_freq_j, fm_work, side='R', invert_tr=.TRUE.)
947 : ELSE
948 0 : CALL cp_fm_to_fm(fm_work, fm_eps_freq_j)
949 0 : CALL cp_fm_power(fm_eps_freq_j, fm_work, -0.5_dp, bs_env%eps_eigval_mat_RI, n_dependent)
950 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
951 0 : fm_Minv_coulomb_factor, fm_eps_freq_j, 0.0_dp, fm_work)
952 0 : IF (bs_env%unit_nr > 0) THEN
953 : WRITE (bs_env%unit_nr, '(T2,A)') &
954 0 : 'Cholesky decomposition failed; screening computed by diagonalization.'
955 0 : WRITE (bs_env%unit_nr, '(T2,A,T72,I9)') 'Discarded eigenmodes:', n_dependent
956 : END IF
957 : END IF
958 :
959 550 : CALL cp_fm_to_fm(fm_Minv_V_Minv, fm_W_freq_j)
960 550 : CALL cp_fm_syrk('U', 'N', n_RI, 1.0_dp, fm_work, 1, 1, -1.0_dp, fm_W_freq_j)
961 550 : CALL cp_fm_uplo_to_full(fm_W_freq_j, fm_work)
962 :
963 : ! Cleanup
964 550 : CALL cp_fm_release(fm_work)
965 550 : CALL cp_fm_release(fm_eps_freq_j)
966 :
967 550 : CALL timestop(handle)
968 :
969 550 : END SUBROUTINE compute_fm_W_freq
970 :
971 : ! **************************************************************************************************
972 : !> \brief ...
973 : !> \param fm ...
974 : !> \param bs_env ...
975 : !> \param mat_name ...
976 : !> \param idx ...
977 : ! **************************************************************************************************
978 248 : SUBROUTINE fm_read(fm, bs_env, mat_name, idx)
979 : TYPE(cp_fm_type) :: fm
980 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
981 : CHARACTER(LEN=*) :: mat_name
982 : INTEGER :: idx
983 :
984 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_read'
985 :
986 : CHARACTER(LEN=default_path_length) :: f_chi
987 : INTEGER :: handle, unit_nr
988 :
989 248 : CALL timeset(routineN, handle)
990 :
991 248 : unit_nr = -1
992 248 : IF (bs_env%para_env%is_source()) THEN
993 :
994 124 : IF (idx < 10) THEN
995 58 : WRITE (f_chi, '(3A,I1,A)') TRIM(bs_env%prefix), TRIM(mat_name), "_0", idx, ".matrix"
996 66 : ELSE IF (idx < 100) THEN
997 66 : WRITE (f_chi, '(3A,I2,A)') TRIM(bs_env%prefix), TRIM(mat_name), "_", idx, ".matrix"
998 : ELSE
999 0 : CPABORT('Please implement more than 99 time/frequency points.')
1000 : END IF
1001 :
1002 : CALL open_file(file_name=TRIM(f_chi), file_action="READ", file_form="UNFORMATTED", &
1003 124 : file_position="REWIND", file_status="OLD", unit_number=unit_nr)
1004 :
1005 : END IF
1006 :
1007 248 : CALL cp_fm_read_unformatted(fm, unit_nr)
1008 :
1009 248 : IF (bs_env%para_env%is_source()) CALL close_file(unit_number=unit_nr)
1010 :
1011 248 : CALL timestop(handle)
1012 :
1013 248 : END SUBROUTINE fm_read
1014 :
1015 : ! **************************************************************************************************
1016 : !> \brief Write a GW restart matrix when requested by the print key.
1017 : !> \param fm Matrix to write, or reusable conversion buffer when matrix is supplied.
1018 : !> \param matrix_index ...
1019 : !> \param matrix_name ...
1020 : !> \param qs_env ...
1021 : !> \param matrix Optional DBCSR source, converted to fm only when output is enabled.
1022 : ! **************************************************************************************************
1023 5834 : SUBROUTINE fm_write(fm, matrix_index, matrix_name, qs_env, matrix)
1024 : TYPE(cp_fm_type) :: fm
1025 : INTEGER :: matrix_index
1026 : CHARACTER(LEN=*) :: matrix_name
1027 : TYPE(qs_environment_type), POINTER :: qs_env
1028 : TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix
1029 :
1030 : CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
1031 : routineN = 'fm_write'
1032 :
1033 : CHARACTER(LEN=default_path_length) :: filename
1034 : INTEGER :: handle, unit_nr
1035 : TYPE(cp_logger_type), POINTER :: logger
1036 : TYPE(section_vals_type), POINTER :: input
1037 :
1038 5834 : CALL timeset(routineN, handle)
1039 :
1040 5834 : CALL get_qs_env(qs_env, input=input)
1041 :
1042 5834 : logger => cp_get_default_logger()
1043 :
1044 5834 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
1045 :
1046 1148 : IF (PRESENT(matrix)) THEN
1047 814 : CALL cp_fm_set_all(fm, 0.0_dp)
1048 814 : CALL copy_dbcsr_to_fm(matrix, fm)
1049 : END IF
1050 :
1051 1148 : IF (matrix_index < 10) THEN
1052 568 : WRITE (filename, '(3A,I1)') "RESTART_", matrix_name, "_0", matrix_index
1053 580 : ELSE IF (matrix_index < 100) THEN
1054 580 : WRITE (filename, '(3A,I2)') "RESTART_", matrix_name, "_", matrix_index
1055 : ELSE
1056 0 : CPABORT('Please implement more than 99 time/frequency points.')
1057 : END IF
1058 :
1059 : unit_nr = cp_print_key_unit_nr(logger, input, key, extension=".matrix", &
1060 : file_form="UNFORMATTED", middle_name=TRIM(filename), &
1061 1148 : file_position="REWIND", file_action="WRITE")
1062 :
1063 1148 : CALL cp_fm_write_unformatted(fm, unit_nr)
1064 1148 : IF (unit_nr > 0) THEN
1065 574 : CALL close_file(unit_nr)
1066 : END IF
1067 : END IF
1068 :
1069 5834 : CALL timestop(handle)
1070 :
1071 5834 : END SUBROUTINE fm_write
1072 :
1073 : ! **************************************************************************************************
1074 : !> \brief ...
1075 : !> \param bs_env ...
1076 : !> \param fm_chi_Gamma_freq ...
1077 : !> \param j_w ...
1078 : !> \param mat_chi_Gamma_tau ...
1079 : ! **************************************************************************************************
1080 3690 : SUBROUTINE compute_fm_chi_Gamma_freq(bs_env, fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
1081 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1082 : TYPE(cp_fm_type) :: fm_chi_Gamma_freq
1083 : INTEGER :: j_w
1084 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
1085 :
1086 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_chi_Gamma_freq'
1087 :
1088 : INTEGER :: handle, i_t
1089 : REAL(KIND=dp) :: weight_ij
1090 :
1091 3690 : CALL timeset(routineN, handle)
1092 :
1093 3690 : CALL dbcsr_set(bs_env%mat_RI_RI%matrix, 0.0_dp)
1094 :
1095 53844 : DO i_t = 1, bs_env%num_time_freq_points
1096 :
1097 50154 : weight_ij = bs_env%time_frequency_grid%cosine_time_to_frequency(j_w, i_t)
1098 :
1099 : ! Fourier transform
1100 : CALL dbcsr_add(bs_env%mat_RI_RI%matrix, mat_chi_Gamma_tau(i_t)%matrix, &
1101 53844 : 1.0_dp, weight_ij)
1102 :
1103 : END DO
1104 :
1105 3690 : CALL copy_dbcsr_to_fm(bs_env%mat_RI_RI%matrix, fm_chi_Gamma_freq)
1106 :
1107 3690 : CALL timestop(handle)
1108 :
1109 3690 : END SUBROUTINE compute_fm_chi_Gamma_freq
1110 :
1111 : ! **************************************************************************************************
1112 : !> \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
1113 : !> frequency, into bs_env%fm_W_MIC_freq_zero
1114 : !> \param bs_env time grid and its weights; receives fm_W_MIC_freq_zero
1115 : !> \param fm_W_MIC_time W^MIC_PQ(iτ_j) for every time point j
1116 : ! **************************************************************************************************
1117 64 : SUBROUTINE sum_W_MIC_freq_zero(bs_env, fm_W_MIC_time)
1118 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1119 : TYPE(cp_fm_type), DIMENSION(:) :: fm_W_MIC_time
1120 :
1121 : CHARACTER(LEN=*), PARAMETER :: routineN = 'sum_W_MIC_freq_zero'
1122 :
1123 : INTEGER :: handle, i_t
1124 : REAL(KIND=dp) :: t_start
1125 :
1126 64 : CALL timeset(routineN, handle)
1127 :
1128 64 : t_start = m_walltime()
1129 64 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
1130 64 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
1131 1344 : DO i_t = 1, bs_env%num_time_freq_points
1132 : CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
1133 : bs_env%time_frequency_grid%time_weights_at_zero_frequency(i_t), &
1134 1344 : fm_W_MIC_time(i_t))
1135 : END DO
1136 64 : IF (bs_env%unit_nr > 0) THEN
1137 : WRITE (bs_env%unit_nr, '(T2,A,I11,A,I3,A,F10.1,A)') &
1138 32 : 'Computed W(ω=0,k) for k-point batch', &
1139 32 : 1, ' /', 1, &
1140 64 : ', Execution time', m_walltime() - t_start, ' s'
1141 : END IF
1142 :
1143 64 : CALL timestop(handle)
1144 :
1145 64 : END SUBROUTINE sum_W_MIC_freq_zero
1146 :
1147 : ! **************************************************************************************************
1148 : !> \brief ...
1149 : !> \param bs_env ...
1150 : !> \param fm_W_MIC_time ...
1151 : ! **************************************************************************************************
1152 114 : SUBROUTINE create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
1153 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1154 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1155 :
1156 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_fm_W_MIC_time'
1157 :
1158 : INTEGER :: handle, i_t
1159 :
1160 114 : CALL timeset(routineN, handle)
1161 :
1162 2256 : ALLOCATE (fm_W_MIC_time(bs_env%num_time_freq_points))
1163 2028 : DO i_t = 1, bs_env%num_time_freq_points
1164 2028 : CALL cp_fm_create(fm_W_MIC_time(i_t), bs_env%fm_RI_RI%matrix_struct, set_zero=.TRUE.)
1165 : END DO
1166 :
1167 114 : CALL timestop(handle)
1168 :
1169 114 : END SUBROUTINE create_fm_W_MIC_time
1170 :
1171 : ! **************************************************************************************************
1172 : !> \brief ...
1173 : !> \param bs_env ...
1174 : !> \param fm_W_MIC_time ...
1175 : !> \param fm_W_MIC_freq_j ...
1176 : !> \param j_w ...
1177 : ! **************************************************************************************************
1178 3690 : SUBROUTINE Fourier_transform_w_to_t(bs_env, fm_W_MIC_time, fm_W_MIC_freq_j, j_w)
1179 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1180 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_time
1181 : TYPE(cp_fm_type) :: fm_W_MIC_freq_j
1182 : INTEGER :: j_w
1183 :
1184 : CHARACTER(LEN=*), PARAMETER :: routineN = 'Fourier_transform_w_to_t'
1185 :
1186 : INTEGER :: handle, i_t
1187 : REAL(KIND=dp) :: weight_ij
1188 :
1189 3690 : CALL timeset(routineN, handle)
1190 :
1191 53844 : DO i_t = 1, bs_env%num_time_freq_points
1192 :
1193 50154 : weight_ij = bs_env%time_frequency_grid%cosine_frequency_to_time(i_t, j_w)
1194 :
1195 : ! actual Fourier transform
1196 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_W_MIC_time(i_t), &
1197 53844 : beta=weight_ij, matrix_b=fm_W_MIC_freq_j)
1198 :
1199 : END DO
1200 :
1201 3690 : CALL timestop(handle)
1202 :
1203 3690 : END SUBROUTINE Fourier_transform_w_to_t
1204 :
1205 : ! **************************************************************************************************
1206 : !> \brief ...
1207 : !> \param V_tr_R ...
1208 : !> \param pot_type ...
1209 : !> \param regularization_RI ...
1210 : !> \param bs_env ...
1211 : !> \param qs_env ...
1212 : ! **************************************************************************************************
1213 168 : SUBROUTINE get_V_tr_R(V_tr_R, pot_type, regularization_RI, bs_env, qs_env)
1214 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: V_tr_R
1215 : TYPE(libint_potential_type) :: pot_type
1216 : REAL(KIND=dp) :: regularization_RI
1217 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1218 : TYPE(qs_environment_type), POINTER :: qs_env
1219 :
1220 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_V_tr_R'
1221 :
1222 : INTEGER :: handle, img, nimages_scf_desymm
1223 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_RI
1224 168 : INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
1225 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1226 168 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_V_tr_R
1227 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
1228 168 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: mat_V_tr_R
1229 : TYPE(distribution_2d_type), POINTER :: dist_2d
1230 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1231 168 : POINTER :: sab_RI
1232 168 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1233 168 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1234 :
1235 168 : CALL timeset(routineN, handle)
1236 :
1237 168 : NULLIFY (sab_RI, dist_2d)
1238 :
1239 : CALL get_qs_env(qs_env=qs_env, &
1240 : blacs_env=blacs_env, &
1241 : distribution_2d=dist_2d, &
1242 : qs_kind_set=qs_kind_set, &
1243 168 : particle_set=particle_set)
1244 :
1245 504 : ALLOCATE (sizes_RI(bs_env%n_atom))
1246 168 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=bs_env%basis_set_RI)
1247 : CALL build_2c_neighbor_lists(sab_RI, bs_env%basis_set_RI, bs_env%basis_set_RI, &
1248 : pot_type, "2c_nl_RI", qs_env, sym_ij=.FALSE., &
1249 168 : dist_2d=dist_2d)
1250 168 : CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
1251 504 : ALLOCATE (row_bsize(SIZE(sizes_RI)))
1252 336 : ALLOCATE (col_bsize(SIZE(sizes_RI)))
1253 644 : row_bsize(:) = sizes_RI
1254 644 : col_bsize(:) = sizes_RI
1255 :
1256 168 : nimages_scf_desymm = bs_env%nimages_scf_desymm
1257 2016 : ALLOCATE (mat_V_tr_R(nimages_scf_desymm))
1258 : CALL dbcsr_create(mat_V_tr_R(1), "(RI|RI)", dbcsr_dist, dbcsr_type_no_symmetry, &
1259 168 : row_bsize, col_bsize)
1260 168 : DEALLOCATE (row_bsize, col_bsize)
1261 :
1262 1512 : DO img = 2, nimages_scf_desymm
1263 1512 : CALL dbcsr_create(mat_V_tr_R(img), template=mat_V_tr_R(1))
1264 : END DO
1265 :
1266 : CALL build_2c_integrals(mat_V_tr_R, 0.0_dp, qs_env, sab_RI, bs_env%basis_set_RI, &
1267 : bs_env%basis_set_RI, pot_type, do_kpoints=.TRUE., &
1268 : ext_kpoints=bs_env%kpoints_scf_desymm, &
1269 168 : regularization_RI=regularization_RI)
1270 :
1271 2016 : ALLOCATE (fm_V_tr_R(nimages_scf_desymm))
1272 1680 : DO img = 1, nimages_scf_desymm
1273 1512 : CALL cp_fm_create(fm_V_tr_R(img), bs_env%fm_RI_RI%matrix_struct)
1274 1512 : CALL copy_dbcsr_to_fm(mat_V_tr_R(img), fm_V_tr_R(img))
1275 1680 : CALL dbcsr_release(mat_V_tr_R(img))
1276 : END DO
1277 :
1278 168 : IF (.NOT. ALLOCATED(V_tr_R)) THEN
1279 840 : ALLOCATE (V_tr_R(bs_env%n_RI, bs_env%n_RI, nimages_scf_desymm))
1280 : END IF
1281 :
1282 168 : CALL cp_fm_to_replicated_array(fm_V_tr_R, V_tr_R)
1283 :
1284 168 : CALL cp_fm_release(fm_V_tr_R)
1285 168 : CALL dbcsr_distribution_release(dbcsr_dist)
1286 168 : CALL release_neighbor_list_sets(sab_RI)
1287 :
1288 168 : CALL timestop(handle)
1289 :
1290 504 : END SUBROUTINE get_V_tr_R
1291 :
1292 : END MODULE gw_compute_V_and_W
|
