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 Routines treating GW and RPA calculations with kpoints
10 : !> \par History
11 : !> since 2018 continuous development [J. Wilhelm]
12 : ! **************************************************************************************************
13 : MODULE rpa_gw_kpoints_util
14 : USE cell_types, ONLY: cell_type,&
15 : get_cell,&
16 : pbc
17 : USE cp_blacs_env, ONLY: cp_blacs_env_type
18 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
19 : cp_cfm_scale_and_add_fm,&
20 : cp_cfm_uplo_to_full
21 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
22 : cp_cfm_cholesky_invert
23 : USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
24 : cp_cfm_geeig_canon,&
25 : cp_cfm_heevd
26 : USE cp_cfm_types, ONLY: cp_cfm_create,&
27 : cp_cfm_get_info,&
28 : cp_cfm_release,&
29 : cp_cfm_set_all,&
30 : cp_cfm_to_cfm,&
31 : cp_cfm_to_fm,&
32 : cp_cfm_type
33 : USE cp_control_types, ONLY: dft_control_type
34 : USE cp_dbcsr_api, ONLY: &
35 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_filter, &
36 : dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
37 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
38 : dbcsr_release, dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
39 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
40 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
41 : copy_fm_to_dbcsr,&
42 : dbcsr_allocate_matrix_set
43 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
44 : USE cp_fm_struct, ONLY: cp_fm_struct_type
45 : USE cp_fm_types, ONLY: cp_fm_copy_general,&
46 : cp_fm_create,&
47 : cp_fm_release,&
48 : cp_fm_set_all,&
49 : cp_fm_type
50 : USE hfx_types, ONLY: hfx_release
51 : USE input_constants, ONLY: cholesky_off,&
52 : kp_weights_W_auto,&
53 : kp_weights_W_tailored,&
54 : kp_weights_W_uniform
55 : USE kinds, ONLY: dp
56 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
57 : kpoint_initialize_mo_set,&
58 : kpoint_initialize_mos
59 : USE kpoint_types, ONLY: get_kpoint_info,&
60 : kpoint_env_type,&
61 : kpoint_type
62 : USE machine, ONLY: m_walltime
63 : USE mathconstants, ONLY: gaussi,&
64 : twopi,&
65 : z_one,&
66 : z_zero
67 : USE mathlib, ONLY: invmat
68 : USE message_passing, ONLY: mp_para_env_type
69 : USE parallel_gemm_api, ONLY: parallel_gemm
70 : USE particle_types, ONLY: particle_type
71 : USE qs_band_structure, ONLY: calculate_kpoints_for_bs
72 : USE qs_environment_types, ONLY: get_qs_env,&
73 : qs_environment_type
74 : USE qs_mo_types, ONLY: get_mo_set
75 : USE qs_scf_types, ONLY: qs_scf_env_type
76 : USE rpa_gw_im_time_util, ONLY: compute_weight_re_im,&
77 : get_atom_index_from_basis_function_index
78 : USE rpa_im_time, ONLY: init_cell_index_rpa
79 : USE scf_control_types, ONLY: scf_control_type
80 : USE time_frequency_grids, ONLY: time_frequency_grid_type
81 : #include "./base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 :
85 : PRIVATE
86 :
87 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_kpoints_util'
88 :
89 : PUBLIC :: invert_eps_compute_W_and_Erpa_kp, cp_cfm_power, real_space_to_kpoint_transform_rpa, &
90 : get_mat_cell_T_from_mat_gamma, get_bandstruc_and_k_dependent_MOs, &
91 : compute_wkp_W, mat_kp_from_mat_gamma
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief ...
97 : !> \param dimen_RI ...
98 : !> \param jquad ...
99 : !> \param nkp ...
100 : !> \param count_ev_sc_GW ...
101 : !> \param para_env ...
102 : !> \param Erpa ...
103 : !> \param grid ...
104 : !> \param wkp_W ...
105 : !> \param do_gw_im_time ...
106 : !> \param do_ri_Sigma_x ...
107 : !> \param do_kpoints_from_Gamma ...
108 : !> \param cfm_mat_Q ...
109 : !> \param ikp_local ...
110 : !> \param mat_P_omega ...
111 : !> \param mat_P_omega_kp ...
112 : !> \param qs_env ...
113 : !> \param eps_filter_im_time ...
114 : !> \param unit_nr ...
115 : !> \param kpoints ...
116 : !> \param fm_mat_Minv_L_kpoints ...
117 : !> \param fm_matrix_L_kpoints ...
118 : !> \param fm_mat_W ...
119 : !> \param fm_mat_RI_global_work ...
120 : !> \param mat_MinvVMinv ...
121 : !> \param fm_matrix_Minv ...
122 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
123 : ! **************************************************************************************************
124 132 : SUBROUTINE invert_eps_compute_W_and_Erpa_kp(dimen_RI, jquad, nkp, count_ev_sc_GW, para_env, &
125 264 : Erpa, grid, wkp_W, do_gw_im_time, &
126 : do_ri_Sigma_x, do_kpoints_from_Gamma, &
127 132 : cfm_mat_Q, ikp_local, mat_P_omega, mat_P_omega_kp, &
128 : qs_env, eps_filter_im_time, unit_nr, kpoints, fm_mat_Minv_L_kpoints, &
129 132 : fm_matrix_L_kpoints, fm_mat_W, &
130 : fm_mat_RI_global_work, mat_MinvVMinv, fm_matrix_Minv, &
131 : fm_matrix_Minv_Vtrunc_Minv)
132 :
133 : INTEGER, INTENT(IN) :: dimen_RI, jquad, nkp, count_ev_sc_GW
134 : TYPE(mp_para_env_type), POINTER :: para_env
135 : REAL(KIND=dp), INTENT(INOUT) :: Erpa
136 : TYPE(time_frequency_grid_type), INTENT(IN) :: grid
137 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_W
138 : LOGICAL, INTENT(IN) :: do_gw_im_time, do_ri_Sigma_x, &
139 : do_kpoints_from_Gamma
140 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_Q
141 : INTEGER, INTENT(IN) :: ikp_local
142 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega, mat_P_omega_kp
143 : TYPE(qs_environment_type), POINTER :: qs_env
144 : REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time
145 : INTEGER, INTENT(IN) :: unit_nr
146 : TYPE(kpoint_type), POINTER :: kpoints
147 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_Minv_L_kpoints, &
148 : fm_matrix_L_kpoints
149 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_W
150 : TYPE(cp_fm_type) :: fm_mat_RI_global_work
151 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_MinvVMinv
152 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv, &
153 : fm_matrix_Minv_Vtrunc_Minv
154 :
155 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invert_eps_compute_W_and_Erpa_kp'
156 :
157 : INTEGER :: handle, ikp, num_integ_points
158 : LOGICAL :: do_this_ikp
159 : REAL(KIND=dp) :: t1, t2
160 132 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: trace_Qomega
161 :
162 132 : CALL timeset(routineN, handle)
163 :
164 132 : num_integ_points = SIZE(grid%frequency)
165 :
166 132 : t1 = m_walltime()
167 :
168 132 : IF (do_kpoints_from_Gamma) THEN
169 96 : CALL get_mat_cell_T_from_mat_gamma(mat_P_omega(jquad, :), qs_env, kpoints, jquad, unit_nr)
170 : END IF
171 :
172 : CALL transform_P_from_real_space_to_kpoints(mat_P_omega, mat_P_omega_kp, &
173 132 : kpoints, eps_filter_im_time, jquad)
174 :
175 396 : ALLOCATE (trace_Qomega(dimen_RI))
176 :
177 132 : IF (unit_nr > 0) WRITE (unit_nr, '(/T3,A,1X,I3)') &
178 66 : 'GW_INFO| Computing chi and W frequency point', jquad
179 :
180 2700 : DO ikp = 1, nkp
181 :
182 : ! parallization, we either have all kpoints on all processors or a single kpoint per group
183 2568 : do_this_ikp = (ikp_local == -1) .OR. (ikp_local == 0 .AND. ikp == 1) .OR. (ikp_local == ikp)
184 : IF (.NOT. do_this_ikp) CYCLE
185 :
186 : ! 1. remove all spurious negative eigenvalues from P(iw,k), multiplication Q(iw,k) = K^H(k)P(iw,k)K(k)
187 : CALL compute_Q_kp_RPA(cfm_mat_Q, &
188 : mat_P_omega_kp, &
189 : fm_mat_Minv_L_kpoints(ikp, 1), &
190 : fm_mat_Minv_L_kpoints(ikp, 2), &
191 : fm_mat_RI_global_work, &
192 : dimen_RI, ikp, nkp, ikp_local, para_env, &
193 2568 : qs_env%mp2_env%ri_rpa_im_time%make_chi_pos_definite)
194 :
195 : ! 2. Cholesky decomposition of Id + Q(iw,k)
196 2568 : CALL cholesky_decomp_Q(cfm_mat_Q, para_env, trace_Qomega, dimen_RI)
197 :
198 : ! 3. Computing E_c^RPA = E_c^RPA + a_w/N_k*sum_k ln[det(1+Q(iw,k))-Tr(Q(iw,k))]
199 : CALL frequency_and_kpoint_integration(Erpa, cfm_mat_Q, para_env, trace_Qomega, &
200 2568 : dimen_RI, grid%frequency_weights(jquad), kpoints%wkp(ikp))
201 :
202 2700 : IF (do_gw_im_time) THEN
203 :
204 : ! compute S^-1*V*S^-1 for exchange part of the self-energy in real space as W in real space
205 2496 : IF (do_ri_Sigma_x .AND. jquad == 1 .AND. count_ev_sc_GW == 1 .AND. do_kpoints_from_Gamma) THEN
206 :
207 312 : CALL dbcsr_set(mat_MinvVMinv%matrix, 0.0_dp)
208 312 : CALL copy_fm_to_dbcsr(fm_matrix_Minv_Vtrunc_Minv(1, 1), mat_MinvVMinv%matrix, keep_sparsity=.FALSE.)
209 :
210 : END IF
211 2496 : IF (do_kpoints_from_Gamma) THEN
212 : CALL compute_Wc_real_space_tau_GW(fm_mat_W, cfm_mat_Q, &
213 : fm_matrix_L_kpoints(ikp, 1), &
214 : fm_matrix_L_kpoints(ikp, 2), &
215 : dimen_RI, jquad, &
216 : ikp, grid, &
217 2496 : ikp_local, para_env, kpoints, qs_env, wkp_W)
218 : END IF
219 :
220 : END IF
221 : END DO
222 :
223 : ! after the transform of (eps(iw)-1)^-1 from iw to it is done, multiply with V^1/2 to obtain W(it)
224 132 : IF (do_gw_im_time .AND. do_kpoints_from_Gamma .AND. jquad == num_integ_points) THEN
225 16 : CALL Wc_to_Minv_Wc_Minv(fm_mat_W, fm_matrix_Minv, para_env, dimen_RI, num_integ_points)
226 16 : CALL deallocate_kp_matrices(fm_matrix_L_kpoints, fm_mat_Minv_L_kpoints)
227 : END IF
228 :
229 132 : DEALLOCATE (trace_Qomega)
230 :
231 132 : t2 = m_walltime()
232 :
233 132 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,T56,F25.1)') 'Execution time (s):', t2 - t1
234 :
235 132 : CALL timestop(handle)
236 :
237 132 : END SUBROUTINE invert_eps_compute_W_and_Erpa_kp
238 :
239 : ! **************************************************************************************************
240 : !> \brief ...
241 : !> \param fm_matrix_L_kpoints ...
242 : !> \param fm_mat_Minv_L_kpoints ...
243 : ! **************************************************************************************************
244 16 : SUBROUTINE deallocate_kp_matrices(fm_matrix_L_kpoints, fm_mat_Minv_L_kpoints)
245 :
246 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_L_kpoints, &
247 : fm_mat_Minv_L_kpoints
248 :
249 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_kp_matrices'
250 :
251 : INTEGER :: handle
252 :
253 16 : CALL timeset(routineN, handle)
254 :
255 16 : CALL cp_fm_release(fm_mat_Minv_L_kpoints)
256 16 : CALL cp_fm_release(fm_matrix_L_kpoints)
257 :
258 16 : CALL timestop(handle)
259 :
260 16 : END SUBROUTINE deallocate_kp_matrices
261 :
262 : ! **************************************************************************************************
263 : !> \brief ...
264 : !> \param matrix ...
265 : !> \param threshold ...
266 : !> \param exponent ...
267 : !> \param min_eigval ...
268 : ! **************************************************************************************************
269 5862 : SUBROUTINE cp_cfm_power(matrix, threshold, exponent, min_eigval)
270 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix
271 : REAL(KIND=dp) :: threshold, exponent
272 : REAL(KIND=dp), OPTIONAL :: min_eigval
273 :
274 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_power'
275 :
276 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_exponent
277 : INTEGER :: handle, i, ncol_global, nrow_global
278 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
279 : TYPE(cp_cfm_type) :: cfm_work
280 :
281 5862 : CALL timeset(routineN, handle)
282 :
283 5862 : CALL cp_cfm_create(cfm_work, matrix%matrix_struct)
284 5862 : CALL cp_cfm_set_all(cfm_work, z_zero)
285 :
286 : ! Test that matrix is square
287 5862 : CALL cp_cfm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global)
288 5862 : CPASSERT(nrow_global == ncol_global)
289 17586 : ALLOCATE (eigenvalues(nrow_global), SOURCE=0.0_dp)
290 17586 : ALLOCATE (eigenvalues_exponent(nrow_global), SOURCE=z_zero)
291 :
292 : ! Diagonalize matrix: get eigenvectors and eigenvalues
293 5862 : CALL cp_cfm_heevd(matrix, cfm_work, eigenvalues)
294 :
295 299618 : DO i = 1, nrow_global
296 299618 : IF (eigenvalues(i) > threshold) THEN
297 265452 : eigenvalues_exponent(i) = CMPLX((eigenvalues(i))**(0.5_dp*exponent), threshold, KIND=dp)
298 : ELSE
299 28304 : IF (PRESENT(min_eigval)) THEN
300 0 : eigenvalues_exponent(i) = CMPLX(min_eigval, 0.0_dp, KIND=dp)
301 : ELSE
302 28304 : eigenvalues_exponent(i) = z_zero
303 : END IF
304 : END IF
305 : END DO
306 :
307 5862 : CALL cp_cfm_column_scale(cfm_work, eigenvalues_exponent)
308 :
309 : CALL parallel_gemm("N", "C", nrow_global, nrow_global, nrow_global, z_one, &
310 5862 : cfm_work, cfm_work, z_zero, matrix)
311 :
312 5862 : DEALLOCATE (eigenvalues, eigenvalues_exponent)
313 :
314 5862 : CALL cp_cfm_release(cfm_work)
315 :
316 5862 : CALL timestop(handle)
317 :
318 11724 : END SUBROUTINE cp_cfm_power
319 :
320 : ! **************************************************************************************************
321 : !> \brief ...
322 : !> \param cfm_mat_Q ...
323 : !> \param mat_P_omega_kp ...
324 : !> \param fm_mat_L_re ...
325 : !> \param fm_mat_L_im ...
326 : !> \param fm_mat_RI_global_work ...
327 : !> \param dimen_RI ...
328 : !> \param ikp ...
329 : !> \param nkp ...
330 : !> \param ikp_local ...
331 : !> \param para_env ...
332 : !> \param make_chi_pos_definite ...
333 : ! **************************************************************************************************
334 2568 : SUBROUTINE compute_Q_kp_RPA(cfm_mat_Q, mat_P_omega_kp, fm_mat_L_re, fm_mat_L_im, &
335 : fm_mat_RI_global_work, dimen_RI, ikp, nkp, ikp_local, para_env, &
336 : make_chi_pos_definite)
337 :
338 : TYPE(cp_cfm_type) :: cfm_mat_Q
339 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega_kp
340 : TYPE(cp_fm_type) :: fm_mat_L_re, fm_mat_L_im, &
341 : fm_mat_RI_global_work
342 : INTEGER, INTENT(IN) :: dimen_RI, ikp, nkp, ikp_local
343 : TYPE(mp_para_env_type), POINTER :: para_env
344 : LOGICAL, INTENT(IN) :: make_chi_pos_definite
345 :
346 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Q_kp_RPA'
347 :
348 : INTEGER :: handle
349 : TYPE(cp_cfm_type) :: cfm_mat_L, cfm_mat_work
350 : TYPE(cp_fm_type) :: fm_mat_work
351 :
352 2568 : CALL timeset(routineN, handle)
353 :
354 2568 : CALL cp_cfm_create(cfm_mat_work, fm_mat_L_re%matrix_struct)
355 2568 : CALL cp_cfm_set_all(cfm_mat_work, z_zero)
356 :
357 2568 : CALL cp_cfm_create(cfm_mat_L, fm_mat_L_re%matrix_struct)
358 2568 : CALL cp_cfm_set_all(cfm_mat_L, z_zero)
359 :
360 2568 : CALL cp_fm_create(fm_mat_work, fm_mat_L_re%matrix_struct)
361 2568 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
362 :
363 : ! 1. Convert the dbcsr matrix mat_P_omega_kp (that is chi(k,iw)) to a full matrix and
364 : ! distribute it to subgroups
365 : CALL mat_P_to_subgroup(mat_P_omega_kp, fm_mat_RI_global_work, &
366 2568 : fm_mat_work, cfm_mat_Q, ikp, nkp, ikp_local, para_env)
367 :
368 : ! 2. Remove all negative eigenvalues from chi(k,iw)
369 2568 : IF (make_chi_pos_definite) THEN
370 2568 : CALL cp_cfm_power(cfm_mat_Q, threshold=0.0_dp, exponent=1.0_dp)
371 : END IF
372 :
373 : ! 3. Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L
374 2568 : CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_L, z_one, fm_mat_L_re)
375 2568 : CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_L, gaussi, fm_mat_L_im)
376 :
377 : ! 4. work = P(iw,k)*L(k)
378 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_Q, cfm_mat_L, &
379 2568 : z_zero, cfm_mat_work)
380 :
381 : ! 5. Q(iw,k) = L^H(k)*work
382 : CALL parallel_gemm('C', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_work, &
383 2568 : z_zero, cfm_mat_Q)
384 :
385 2568 : CALL cp_cfm_release(cfm_mat_work)
386 2568 : CALL cp_cfm_release(cfm_mat_L)
387 2568 : CALL cp_fm_release(fm_mat_work)
388 :
389 2568 : CALL timestop(handle)
390 :
391 2568 : END SUBROUTINE compute_Q_kp_RPA
392 :
393 : ! **************************************************************************************************
394 : !> \brief ...
395 : !> \param mat_P_omega_kp ...
396 : !> \param fm_mat_RI_global_work ...
397 : !> \param fm_mat_work ...
398 : !> \param cfm_mat_Q ...
399 : !> \param ikp ...
400 : !> \param nkp ...
401 : !> \param ikp_local ...
402 : !> \param para_env ...
403 : ! **************************************************************************************************
404 2568 : SUBROUTINE mat_P_to_subgroup(mat_P_omega_kp, fm_mat_RI_global_work, &
405 : fm_mat_work, cfm_mat_Q, ikp, nkp, ikp_local, para_env)
406 :
407 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega_kp
408 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_RI_global_work, fm_mat_work
409 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_Q
410 : INTEGER, INTENT(IN) :: ikp, nkp, ikp_local
411 : TYPE(mp_para_env_type), POINTER :: para_env
412 :
413 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mat_P_to_subgroup'
414 :
415 : INTEGER :: handle, jkp
416 : TYPE(cp_fm_type) :: fm_dummy
417 : TYPE(dbcsr_type), POINTER :: mat_P_omega_im, mat_P_omega_re
418 :
419 2568 : CALL timeset(routineN, handle)
420 :
421 2568 : IF (ikp_local == -1) THEN
422 :
423 2568 : mat_P_omega_re => mat_P_omega_kp(1, ikp)%matrix
424 2568 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
425 2568 : CALL copy_dbcsr_to_fm(mat_P_omega_re, fm_mat_work)
426 2568 : CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_Q, z_one, fm_mat_work)
427 :
428 2568 : mat_P_omega_im => mat_P_omega_kp(2, ikp)%matrix
429 2568 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
430 2568 : CALL copy_dbcsr_to_fm(mat_P_omega_im, fm_mat_work)
431 2568 : CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_Q, gaussi, fm_mat_work)
432 :
433 : ELSE
434 :
435 0 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
436 :
437 0 : DO jkp = 1, nkp
438 :
439 0 : mat_P_omega_re => mat_P_omega_kp(1, jkp)%matrix
440 :
441 0 : CALL cp_fm_set_all(fm_mat_RI_global_work, 0.0_dp)
442 0 : CALL copy_dbcsr_to_fm(mat_P_omega_re, fm_mat_RI_global_work)
443 :
444 0 : CALL para_env%sync()
445 :
446 0 : IF (ikp_local == jkp) THEN
447 0 : CALL cp_fm_copy_general(fm_mat_RI_global_work, fm_mat_work, para_env)
448 : ELSE
449 0 : CALL cp_fm_copy_general(fm_mat_RI_global_work, fm_dummy, para_env)
450 : END IF
451 :
452 0 : CALL para_env%sync()
453 :
454 : END DO
455 :
456 0 : CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_Q, z_one, fm_mat_work)
457 :
458 0 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
459 :
460 0 : DO jkp = 1, nkp
461 :
462 0 : mat_P_omega_im => mat_P_omega_kp(2, jkp)%matrix
463 :
464 0 : CALL cp_fm_set_all(fm_mat_RI_global_work, 0.0_dp)
465 0 : CALL copy_dbcsr_to_fm(mat_P_omega_im, fm_mat_RI_global_work)
466 :
467 0 : CALL para_env%sync()
468 :
469 0 : IF (ikp_local == jkp) THEN
470 0 : CALL cp_fm_copy_general(fm_mat_RI_global_work, fm_mat_work, para_env)
471 : ELSE
472 0 : CALL cp_fm_copy_general(fm_mat_RI_global_work, fm_dummy, para_env)
473 : END IF
474 :
475 0 : CALL para_env%sync()
476 :
477 : END DO
478 :
479 0 : CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_Q, gaussi, fm_mat_work)
480 :
481 0 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
482 :
483 : END IF
484 :
485 2568 : CALL para_env%sync()
486 :
487 2568 : CALL timestop(handle)
488 :
489 2568 : END SUBROUTINE mat_P_to_subgroup
490 :
491 : ! **************************************************************************************************
492 : !> \brief ...
493 : !> \param cfm_mat_Q ...
494 : !> \param para_env ...
495 : !> \param trace_Qomega ...
496 : !> \param dimen_RI ...
497 : ! **************************************************************************************************
498 2568 : SUBROUTINE cholesky_decomp_Q(cfm_mat_Q, para_env, trace_Qomega, dimen_RI)
499 :
500 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_Q
501 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
502 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: trace_Qomega
503 : INTEGER, INTENT(IN) :: dimen_RI
504 :
505 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cholesky_decomp_Q'
506 :
507 : INTEGER :: handle, i_global, iiB, info_chol, &
508 : j_global, jjB, ncol_local, nrow_local
509 2568 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
510 : TYPE(cp_cfm_type) :: cfm_mat_Q_tmp, cfm_mat_work
511 :
512 2568 : CALL timeset(routineN, handle)
513 :
514 2568 : CALL cp_cfm_create(cfm_mat_work, cfm_mat_Q%matrix_struct)
515 2568 : CALL cp_cfm_set_all(cfm_mat_work, z_zero)
516 :
517 2568 : CALL cp_cfm_create(cfm_mat_Q_tmp, cfm_mat_Q%matrix_struct)
518 2568 : CALL cp_cfm_set_all(cfm_mat_Q_tmp, z_zero)
519 :
520 : ! get info of fm_mat_Q
521 : CALL cp_cfm_get_info(matrix=cfm_mat_Q, &
522 : nrow_local=nrow_local, &
523 : ncol_local=ncol_local, &
524 : row_indices=row_indices, &
525 2568 : col_indices=col_indices)
526 :
527 : ! calculate the trace of Q and add 1 on the diagonal
528 182640 : trace_Qomega = 0.0_dp
529 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
530 2568 : !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,trace_Qomega,cfm_mat_Q,dimen_RI)
531 : DO jjB = 1, ncol_local
532 : j_global = col_indices(jjB)
533 : DO iiB = 1, nrow_local
534 : i_global = row_indices(iiB)
535 : IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
536 : trace_Qomega(i_global) = REAL(cfm_mat_Q%local_data(iiB, jjB))
537 : cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB) + z_one
538 : END IF
539 : END DO
540 : END DO
541 362712 : CALL para_env%sum(trace_Qomega)
542 :
543 2568 : CALL cp_cfm_to_cfm(cfm_mat_Q, cfm_mat_Q_tmp)
544 :
545 2568 : CALL cp_cfm_cholesky_decompose(matrix=cfm_mat_Q, n=dimen_RI, info_out=info_chol)
546 :
547 2568 : CPASSERT(info_chol == 0)
548 :
549 2568 : CALL cp_cfm_release(cfm_mat_work)
550 2568 : CALL cp_cfm_release(cfm_mat_Q_tmp)
551 :
552 2568 : CALL timestop(handle)
553 :
554 2568 : END SUBROUTINE cholesky_decomp_Q
555 :
556 : ! **************************************************************************************************
557 : !> \brief ...
558 : !> \param Erpa ...
559 : !> \param cfm_mat_Q ...
560 : !> \param para_env ...
561 : !> \param trace_Qomega ...
562 : !> \param dimen_RI ...
563 : !> \param freq_weight ...
564 : !> \param kp_weight ...
565 : ! **************************************************************************************************
566 2568 : SUBROUTINE frequency_and_kpoint_integration(Erpa, cfm_mat_Q, para_env, trace_Qomega, &
567 : dimen_RI, freq_weight, kp_weight)
568 :
569 : REAL(KIND=dp), INTENT(INOUT) :: Erpa
570 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_Q
571 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
572 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: trace_Qomega
573 : INTEGER, INTENT(IN) :: dimen_RI
574 : REAL(KIND=dp), INTENT(IN) :: freq_weight, kp_weight
575 :
576 : CHARACTER(LEN=*), PARAMETER :: routineN = 'frequency_and_kpoint_integration'
577 :
578 : INTEGER :: handle, i_global, iiB, j_global, jjB, &
579 : ncol_local, nrow_local
580 2568 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
581 : REAL(KIND=dp) :: FComega
582 2568 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Q_log
583 :
584 2568 : CALL timeset(routineN, handle)
585 :
586 : ! get info of cholesky_decomposed(fm_mat_Q)
587 : CALL cp_cfm_get_info(matrix=cfm_mat_Q, &
588 : nrow_local=nrow_local, &
589 : ncol_local=ncol_local, &
590 : row_indices=row_indices, &
591 2568 : col_indices=col_indices)
592 :
593 7704 : ALLOCATE (Q_log(dimen_RI))
594 2568 : Q_log = 0.0_dp
595 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
596 2568 : !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Q_log,cfm_mat_Q,dimen_RI)
597 : DO jjB = 1, ncol_local
598 : j_global = col_indices(jjB)
599 : DO iiB = 1, nrow_local
600 : i_global = row_indices(iiB)
601 : IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
602 : Q_log(i_global) = 2.0_dp*LOG(REAL(cfm_mat_Q%local_data(iiB, jjB)))
603 : END IF
604 : END DO
605 : END DO
606 2568 : CALL para_env%sum(Q_log)
607 :
608 2568 : FComega = 0.0_dp
609 182640 : DO iiB = 1, dimen_RI
610 180072 : IF (MODULO(iiB, para_env%num_pe) /= para_env%mepos) CYCLE
611 : ! FComega=FComega+(LOG(Q_log(iiB))-trace_Qomega(iiB))/2.0_dp
612 182640 : FComega = FComega + (Q_log(iiB) - trace_Qomega(iiB))/2.0_dp
613 : END DO
614 :
615 2568 : Erpa = Erpa + FComega*freq_weight*kp_weight
616 :
617 2568 : DEALLOCATE (Q_log)
618 :
619 2568 : CALL timestop(handle)
620 :
621 5136 : END SUBROUTINE frequency_and_kpoint_integration
622 :
623 : ! **************************************************************************************************
624 : !> \brief ...
625 : !> \param mat_P_omega ...
626 : !> \param qs_env ...
627 : !> \param kpoints ...
628 : !> \param jquad ...
629 : !> \param unit_nr ...
630 : ! **************************************************************************************************
631 440 : SUBROUTINE get_mat_cell_T_from_mat_gamma(mat_P_omega, qs_env, kpoints, jquad, unit_nr)
632 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mat_P_omega
633 : TYPE(qs_environment_type), POINTER :: qs_env
634 : TYPE(kpoint_type), POINTER :: kpoints
635 : INTEGER, INTENT(IN) :: jquad, unit_nr
636 :
637 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_cell_T_from_mat_gamma'
638 :
639 : INTEGER :: col, handle, i_cell, i_dim, j_cell, &
640 : num_cells_P, num_integ_points, row
641 : INTEGER, DIMENSION(3) :: cell_grid_P, periodic
642 440 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_P
643 : LOGICAL :: i_cell_is_the_minimum_image_cell
644 : REAL(KIND=dp) :: abs_rab_cell_i, abs_rab_cell_j
645 : REAL(KIND=dp), DIMENSION(3) :: cell_vector, cell_vector_j, rab_cell_i, &
646 : rab_cell_j
647 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
648 440 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
649 : TYPE(cell_type), POINTER :: cell
650 : TYPE(dbcsr_iterator_type) :: iter
651 440 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
652 :
653 440 : CALL timeset(routineN, handle)
654 :
655 440 : NULLIFY (cell, particle_set)
656 : CALL get_qs_env(qs_env, cell=cell, &
657 440 : particle_set=particle_set)
658 440 : CALL get_cell(cell=cell, h=hmat, periodic=periodic)
659 :
660 1760 : DO i_dim = 1, 3
661 : ! we have at most 3 neigboring cells per dimension and at least one because
662 : ! the density response at Gamma is only divided to neighboring
663 1760 : IF (periodic(i_dim) == 1) THEN
664 880 : cell_grid_P(i_dim) = MAX(MIN((kpoints%nkp_grid(i_dim)/2)*2 - 1, 1), 3)
665 : ELSE
666 440 : cell_grid_P(i_dim) = 1
667 : END IF
668 : END DO
669 :
670 : ! overwrite the cell indices in kpoints
671 440 : CALL init_cell_index_rpa(cell_grid_P, kpoints%cell_to_index, kpoints%index_to_cell, cell)
672 :
673 440 : index_to_cell_P => kpoints%index_to_cell
674 :
675 440 : num_cells_P = SIZE(index_to_cell_P, 2)
676 :
677 440 : num_integ_points = SIZE(mat_P_omega, 1)
678 :
679 : ! first, copy the Gamma-only result from mat_P_omega(1) into all other matrices and
680 : ! remove the blocks later which do not belong to the cell index
681 3960 : DO i_cell = 2, num_cells_P
682 : CALL dbcsr_copy(mat_P_omega(i_cell)%matrix, &
683 3960 : mat_P_omega(1)%matrix)
684 : END DO
685 :
686 440 : IF (jquad == 1 .AND. unit_nr > 0) THEN
687 8 : WRITE (unit_nr, '(T3,A,T66,ES15.2)') 'GW_INFO| RI regularization parameter: ', &
688 16 : qs_env%mp2_env%ri_rpa_im_time%regularization_RI
689 8 : WRITE (unit_nr, '(T3,A,T66,ES15.2)') 'GW_INFO| eps_eigval_S: ', &
690 16 : qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S
691 8 : IF (qs_env%mp2_env%ri_rpa_im_time%make_chi_pos_definite) THEN
692 : WRITE (unit_nr, '(T3,A,T81)') &
693 8 : 'GW_INFO| Make chi(iw,k) positive definite? TRUE'
694 : ELSE
695 : WRITE (unit_nr, '(T3,A,T81)') &
696 0 : 'GW_INFO| Make chi(iw,k) positive definite? FALSE'
697 : END IF
698 :
699 : END IF
700 :
701 4400 : DO i_cell = 1, num_cells_P
702 :
703 3960 : CALL dbcsr_iterator_start(iter, mat_P_omega(i_cell)%matrix)
704 20385 : DO WHILE (dbcsr_iterator_blocks_left(iter))
705 16425 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
706 :
707 262800 : cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell_P(1:3, i_cell), dp))
708 : rab_cell_i(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
709 65700 : (pbc(particle_set(col)%r(1:3), cell) + cell_vector(1:3))
710 16425 : abs_rab_cell_i = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
711 :
712 : ! minimum image convention
713 16425 : i_cell_is_the_minimum_image_cell = .TRUE.
714 164250 : DO j_cell = 1, num_cells_P
715 2365200 : cell_vector_j(1:3) = MATMUL(hmat, REAL(index_to_cell_P(1:3, j_cell), dp))
716 : rab_cell_j(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
717 591300 : (pbc(particle_set(col)%r(1:3), cell) + cell_vector_j(1:3))
718 147825 : abs_rab_cell_j = SQRT(rab_cell_j(1)**2 + rab_cell_j(2)**2 + rab_cell_j(3)**2)
719 :
720 164250 : IF (abs_rab_cell_i > abs_rab_cell_j + 1.0E-6_dp) THEN
721 53748 : i_cell_is_the_minimum_image_cell = .FALSE.
722 : END IF
723 : END DO
724 :
725 32850 : IF (.NOT. i_cell_is_the_minimum_image_cell) THEN
726 2809552 : data_block(:, :) = data_block(:, :)*0.0_dp
727 : END IF
728 :
729 : END DO
730 8360 : CALL dbcsr_iterator_stop(iter)
731 :
732 : END DO
733 :
734 440 : CALL timestop(handle)
735 :
736 440 : END SUBROUTINE get_mat_cell_T_from_mat_gamma
737 :
738 : ! **************************************************************************************************
739 : !> \brief ...
740 : !> \param mat_P_omega ...
741 : !> \param mat_P_omega_kp ...
742 : !> \param kpoints ...
743 : !> \param eps_filter_im_time ...
744 : !> \param jquad ...
745 : ! **************************************************************************************************
746 132 : SUBROUTINE transform_P_from_real_space_to_kpoints(mat_P_omega, mat_P_omega_kp, &
747 : kpoints, eps_filter_im_time, jquad)
748 :
749 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega, mat_P_omega_kp
750 : TYPE(kpoint_type), POINTER :: kpoints
751 : REAL(kind=dp), INTENT(IN) :: eps_filter_im_time
752 : INTEGER, INTENT(IN) :: jquad
753 :
754 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_P_from_real_space_to_kpoints'
755 :
756 : INTEGER :: handle, icell, nkp, num_integ_points
757 :
758 132 : CALL timeset(routineN, handle)
759 :
760 132 : num_integ_points = SIZE(mat_P_omega, 1)
761 132 : nkp = SIZE(mat_P_omega, 2)
762 :
763 : CALL real_space_to_kpoint_transform_rpa(mat_P_omega_kp(1, :), mat_P_omega_kp(2, :), mat_P_omega(jquad, :), &
764 132 : kpoints, eps_filter_im_time)
765 :
766 2700 : DO icell = 1, SIZE(mat_P_omega, 2)
767 2568 : CALL dbcsr_set(mat_P_omega(jquad, icell)%matrix, 0.0_dp)
768 2700 : CALL dbcsr_filter(mat_P_omega(jquad, icell)%matrix, 1.0_dp)
769 : END DO
770 :
771 132 : CALL timestop(handle)
772 :
773 132 : END SUBROUTINE transform_P_from_real_space_to_kpoints
774 :
775 : ! **************************************************************************************************
776 : !> \brief ...
777 : !> \param real_mat_kp ...
778 : !> \param imag_mat_kp ...
779 : !> \param mat_real_space ...
780 : !> \param kpoints ...
781 : !> \param eps_filter_im_time ...
782 : !> \param real_mat_real_space ...
783 : ! **************************************************************************************************
784 476 : SUBROUTINE real_space_to_kpoint_transform_rpa(real_mat_kp, imag_mat_kp, mat_real_space, &
785 : kpoints, eps_filter_im_time, real_mat_real_space)
786 :
787 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: real_mat_kp, imag_mat_kp, mat_real_space
788 : TYPE(kpoint_type), POINTER :: kpoints
789 : REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time
790 : LOGICAL, INTENT(IN), OPTIONAL :: real_mat_real_space
791 :
792 : CHARACTER(LEN=*), PARAMETER :: routineN = 'real_space_to_kpoint_transform_rpa'
793 :
794 : INTEGER :: handle, i_cell, ik, nkp, num_cells
795 : INTEGER, DIMENSION(3) :: cell
796 476 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
797 : LOGICAL :: my_real_mat_real_space
798 : REAL(KIND=dp) :: arg, coskl, sinkl
799 476 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
800 : TYPE(dbcsr_type) :: mat_work
801 :
802 476 : CALL timeset(routineN, handle)
803 :
804 476 : my_real_mat_real_space = .TRUE.
805 476 : IF (PRESENT(real_mat_real_space)) my_real_mat_real_space = real_mat_real_space
806 :
807 : CALL dbcsr_create(matrix=mat_work, &
808 : template=real_mat_kp(1)%matrix, &
809 476 : matrix_type=dbcsr_type_no_symmetry)
810 476 : CALL dbcsr_reserve_all_blocks(mat_work)
811 476 : CALL dbcsr_set(mat_work, 0.0_dp)
812 :
813 : ! this kpoint environme t should be the kpoints for D(it) and X(it) created in init_cell_index_rpa
814 476 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp)
815 :
816 476 : NULLIFY (index_to_cell)
817 476 : index_to_cell => kpoints%index_to_cell
818 :
819 476 : num_cells = SIZE(index_to_cell, 2)
820 :
821 476 : CPASSERT(SIZE(mat_real_space) >= num_cells/2 + 1)
822 :
823 5374 : DO ik = 1, nkp
824 :
825 4898 : CALL dbcsr_reserve_all_blocks(real_mat_kp(ik)%matrix)
826 4898 : CALL dbcsr_reserve_all_blocks(imag_mat_kp(ik)%matrix)
827 :
828 4898 : CALL dbcsr_set(real_mat_kp(ik)%matrix, 0.0_dp)
829 4898 : CALL dbcsr_set(imag_mat_kp(ik)%matrix, 0.0_dp)
830 :
831 29172 : DO i_cell = 1, num_cells/2 + 1
832 :
833 97096 : cell(:) = index_to_cell(:, i_cell)
834 :
835 24274 : arg = REAL(cell(1), dp)*xkp(1, ik) + REAL(cell(2), dp)*xkp(2, ik) + REAL(cell(3), dp)*xkp(3, ik)
836 24274 : coskl = COS(twopi*arg)
837 24274 : sinkl = SIN(twopi*arg)
838 :
839 24274 : IF (my_real_mat_real_space) THEN
840 24034 : CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, coskl)
841 24034 : CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, sinkl)
842 : ELSE
843 240 : CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, -sinkl)
844 240 : CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, coskl)
845 : END IF
846 :
847 29172 : IF (.NOT. (cell(1) == 0 .AND. cell(2) == 0 .AND. cell(3) == 0)) THEN
848 :
849 19376 : CALL dbcsr_transposed(mat_work, mat_real_space(i_cell)%matrix)
850 :
851 19376 : IF (my_real_mat_real_space) THEN
852 19184 : CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_work, 1.0_dp, coskl)
853 19184 : CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_work, 1.0_dp, -sinkl)
854 : ELSE
855 : ! for an imaginary real-space matrix, we need to consider the imaginary unit
856 : ! and we need to take into account that the transposed gives an extra "-" sign
857 : ! because the transposed is actually Hermitian conjugate
858 192 : CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_work, 1.0_dp, -sinkl)
859 192 : CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_work, 1.0_dp, -coskl)
860 : END IF
861 :
862 19376 : CALL dbcsr_set(mat_work, 0.0_dp)
863 :
864 : END IF
865 :
866 : END DO
867 :
868 4898 : CALL dbcsr_filter(real_mat_kp(ik)%matrix, eps_filter_im_time)
869 5374 : CALL dbcsr_filter(imag_mat_kp(ik)%matrix, eps_filter_im_time)
870 :
871 : END DO
872 :
873 476 : CALL dbcsr_release(mat_work)
874 :
875 476 : CALL timestop(handle)
876 :
877 476 : END SUBROUTINE real_space_to_kpoint_transform_rpa
878 :
879 : ! **************************************************************************************************
880 : !> \brief ...
881 : !> \param mat_a ...
882 : !> \param mat_b ...
883 : !> \param alpha ...
884 : !> \param beta ...
885 : ! **************************************************************************************************
886 87300 : SUBROUTINE dbcsr_add_local(mat_a, mat_b, alpha, beta)
887 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_a, mat_b
888 : REAL(kind=dp), INTENT(IN) :: alpha, beta
889 :
890 : INTEGER :: col, row
891 : LOGICAL :: found
892 87300 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_to_compute, data_block
893 : TYPE(dbcsr_iterator_type) :: iter
894 :
895 87300 : CALL dbcsr_iterator_start(iter, mat_b)
896 456354 : DO WHILE (dbcsr_iterator_blocks_left(iter))
897 369054 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
898 :
899 369054 : NULLIFY (block_to_compute)
900 : CALL dbcsr_get_block_p(matrix=mat_a, &
901 369054 : row=row, col=col, block=block_to_compute, found=found)
902 :
903 369054 : CPASSERT(found)
904 :
905 275344560 : block_to_compute(:, :) = alpha*block_to_compute(:, :) + beta*data_block(:, :)
906 :
907 : END DO
908 87300 : CALL dbcsr_iterator_stop(iter)
909 :
910 87300 : END SUBROUTINE dbcsr_add_local
911 :
912 : ! **************************************************************************************************
913 : !> \brief ...
914 : !> \param fm_mat_W_tau ...
915 : !> \param cfm_mat_Q ...
916 : !> \param fm_mat_L_re ...
917 : !> \param fm_mat_L_im ...
918 : !> \param dimen_RI ...
919 : !> \param jquad ...
920 : !> \param ikp ...
921 : !> \param grid ...
922 : !> \param ikp_local ...
923 : !> \param para_env ...
924 : !> \param kpoints ...
925 : !> \param qs_env ...
926 : !> \param wkp_W ...
927 : ! **************************************************************************************************
928 2496 : SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm_mat_L_im, &
929 : dimen_RI, jquad, &
930 : ikp, grid, ikp_local, &
931 2496 : para_env, kpoints, qs_env, wkp_W)
932 :
933 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_W_tau
934 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_Q
935 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_L_re, fm_mat_L_im
936 : INTEGER, INTENT(IN) :: dimen_RI, jquad, ikp
937 : TYPE(time_frequency_grid_type), INTENT(IN) :: grid
938 : INTEGER, INTENT(IN) :: ikp_local
939 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
940 : TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
941 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
942 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_W
943 :
944 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Wc_real_space_tau_GW'
945 :
946 : INTEGER :: handle, handle2, i_global, iatom, iatom_old, iiB, iquad, irow, j_global, jatom, &
947 : jatom_old, jcol, jjB, jkp, ncol_local, nkp, nrow_local, num_cells, num_integ_points
948 2496 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_RI_index
949 2496 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
950 2496 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
951 : REAL(KIND=dp) :: contribution, omega, tau, weight, &
952 : weight_im, weight_re
953 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
954 2496 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
955 2496 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
956 : TYPE(cell_type), POINTER :: cell
957 : TYPE(cp_cfm_type) :: cfm_mat_L, cfm_mat_work, cfm_mat_work_2
958 : TYPE(cp_fm_type) :: fm_dummy, fm_mat_work_global, &
959 : fm_mat_work_local
960 2496 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
961 :
962 2496 : CALL timeset(routineN, handle)
963 :
964 2496 : num_integ_points = SIZE(grid%imaginary_time)
965 :
966 2496 : CALL timeset(routineN//"_1", handle2)
967 :
968 2496 : CALL cp_cfm_create(cfm_mat_work, cfm_mat_Q%matrix_struct)
969 2496 : CALL cp_cfm_set_all(cfm_mat_work, z_zero)
970 :
971 2496 : CALL cp_cfm_create(cfm_mat_work_2, cfm_mat_Q%matrix_struct)
972 2496 : CALL cp_cfm_set_all(cfm_mat_work_2, z_zero)
973 :
974 2496 : CALL cp_cfm_create(cfm_mat_L, cfm_mat_Q%matrix_struct)
975 2496 : CALL cp_cfm_set_all(cfm_mat_L, z_zero)
976 :
977 : ! Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L
978 2496 : CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_L, z_one, fm_mat_L_re)
979 2496 : CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_L, gaussi, fm_mat_L_im)
980 :
981 2496 : CALL cp_fm_create(fm_mat_work_global, fm_mat_W_tau(1)%matrix_struct)
982 2496 : CALL cp_fm_set_all(fm_mat_work_global, 0.0_dp)
983 :
984 2496 : CALL cp_fm_create(fm_mat_work_local, cfm_mat_Q%matrix_struct)
985 2496 : CALL cp_fm_set_all(fm_mat_work_local, 0.0_dp)
986 :
987 2496 : CALL timestop(handle2)
988 :
989 2496 : CALL timeset(routineN//"_2", handle2)
990 :
991 : ! calculate [1+Q(iw')]^-1
992 2496 : CALL cp_cfm_cholesky_invert(cfm_mat_Q)
993 :
994 : ! symmetrize the result
995 2496 : CALL cp_cfm_uplo_to_full(cfm_mat_Q)
996 :
997 : ! subtract exchange part by subtracing identity matrix from epsilon
998 : CALL cp_cfm_get_info(matrix=cfm_mat_Q, &
999 : nrow_local=nrow_local, &
1000 : ncol_local=ncol_local, &
1001 : row_indices=row_indices, &
1002 2496 : col_indices=col_indices)
1003 :
1004 176592 : DO jjB = 1, ncol_local
1005 174096 : j_global = col_indices(jjB)
1006 6950424 : DO iiB = 1, nrow_local
1007 6773832 : i_global = row_indices(iiB)
1008 6947928 : IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
1009 87048 : cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB) - z_one
1010 : END IF
1011 : END DO
1012 : END DO
1013 :
1014 2496 : CALL timestop(handle2)
1015 :
1016 2496 : CALL timeset(routineN//"_3", handle2)
1017 :
1018 : ! work = epsilon(iw,k)*V^1/2(k)
1019 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_Q, cfm_mat_L, &
1020 2496 : z_zero, cfm_mat_work)
1021 :
1022 : ! W(iw,k) = V^1/2(k)*work
1023 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_work, &
1024 2496 : z_zero, cfm_mat_work_2)
1025 :
1026 2496 : CALL timestop(handle2)
1027 :
1028 2496 : CALL timeset(routineN//"_4", handle2)
1029 :
1030 2496 : CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp)
1031 2496 : index_to_cell => kpoints%index_to_cell
1032 2496 : num_cells = SIZE(index_to_cell, 2)
1033 :
1034 2496 : CALL cp_cfm_set_all(cfm_mat_work, z_zero)
1035 :
1036 7488 : ALLOCATE (atom_from_RI_index(dimen_RI))
1037 :
1038 2496 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_RI_index, dimen_RI, "RI_AUX")
1039 :
1040 2496 : NULLIFY (cell, particle_set)
1041 2496 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1042 2496 : CALL get_cell(cell=cell, h=hmat)
1043 2496 : iatom_old = 0
1044 2496 : jatom_old = 0
1045 :
1046 : CALL cp_cfm_get_info(matrix=cfm_mat_Q, &
1047 : nrow_local=nrow_local, &
1048 : ncol_local=ncol_local, &
1049 : row_indices=row_indices, &
1050 2496 : col_indices=col_indices)
1051 :
1052 89544 : DO irow = 1, nrow_local
1053 6863376 : DO jcol = 1, ncol_local
1054 :
1055 6773832 : iatom = atom_from_RI_index(row_indices(irow))
1056 6773832 : jatom = atom_from_RI_index(col_indices(jcol))
1057 :
1058 6773832 : IF (iatom /= iatom_old .OR. jatom /= jatom_old) THEN
1059 :
1060 : ! symmetrize=.FALSE. necessary since we already have a symmetrized index_to_cell
1061 : CALL compute_weight_re_im(weight_re, weight_im, &
1062 : num_cells, iatom, jatom, xkp(1:3, ikp), wkp_W(ikp), &
1063 258336 : cell, index_to_cell, hmat, particle_set)
1064 :
1065 258336 : iatom_old = iatom
1066 258336 : jatom_old = jatom
1067 :
1068 : END IF
1069 :
1070 : contribution = weight_re*REAL(cfm_mat_work_2%local_data(irow, jcol)) + &
1071 6773832 : weight_im*AIMAG(cfm_mat_work_2%local_data(irow, jcol))
1072 :
1073 6860880 : fm_mat_work_local%local_data(irow, jcol) = fm_mat_work_local%local_data(irow, jcol) + contribution
1074 :
1075 : END DO
1076 : END DO
1077 :
1078 2496 : CALL timestop(handle2)
1079 :
1080 2496 : CALL timeset(routineN//"_5", handle2)
1081 :
1082 2496 : IF (ikp_local == -1) THEN
1083 :
1084 2496 : CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env)
1085 :
1086 17472 : DO iquad = 1, num_integ_points
1087 :
1088 14976 : omega = grid%frequency(jquad)
1089 14976 : tau = grid%imaginary_time(iquad)
1090 14976 : weight = grid%cosine_frequency_to_time_weights(iquad, jquad)*COS(tau*omega)
1091 :
1092 14976 : IF (jquad == 1 .AND. ikp == 1) THEN
1093 96 : CALL cp_fm_set_all(matrix=fm_mat_W_tau(iquad), alpha=0.0_dp)
1094 : END IF
1095 :
1096 17472 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_W_tau(iquad), beta=weight, matrix_b=fm_mat_work_global)
1097 :
1098 : END DO
1099 :
1100 : ELSE
1101 :
1102 0 : DO jkp = 1, nkp
1103 :
1104 0 : CALL para_env%sync()
1105 :
1106 0 : IF (ikp_local == jkp) THEN
1107 0 : CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env)
1108 : ELSE
1109 0 : CALL cp_fm_copy_general(fm_dummy, fm_mat_work_global, para_env)
1110 : END IF
1111 :
1112 0 : CALL para_env%sync()
1113 :
1114 0 : DO iquad = 1, num_integ_points
1115 :
1116 0 : omega = grid%frequency(jquad)
1117 0 : tau = grid%imaginary_time(iquad)
1118 0 : weight = grid%cosine_frequency_to_time_weights(iquad, jquad)*COS(tau*omega)
1119 :
1120 0 : IF (jquad == 1 .AND. jkp == 1) THEN
1121 0 : CALL cp_fm_set_all(matrix=fm_mat_W_tau(iquad), alpha=0.0_dp)
1122 : END IF
1123 :
1124 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_W_tau(iquad), beta=weight, &
1125 0 : matrix_b=fm_mat_work_global)
1126 :
1127 : END DO
1128 :
1129 : END DO
1130 :
1131 : END IF
1132 :
1133 2496 : CALL cp_cfm_release(cfm_mat_work)
1134 2496 : CALL cp_cfm_release(cfm_mat_work_2)
1135 2496 : CALL cp_cfm_release(cfm_mat_L)
1136 2496 : CALL cp_fm_release(fm_mat_work_global)
1137 2496 : CALL cp_fm_release(fm_mat_work_local)
1138 :
1139 2496 : DEALLOCATE (atom_from_RI_index)
1140 :
1141 2496 : CALL timestop(handle2)
1142 :
1143 2496 : CALL timestop(handle)
1144 :
1145 27456 : END SUBROUTINE compute_Wc_real_space_tau_GW
1146 :
1147 : ! **************************************************************************************************
1148 : !> \brief ...
1149 : !> \param fm_mat_W ...
1150 : !> \param fm_matrix_Minv ...
1151 : !> \param para_env ...
1152 : !> \param dimen_RI ...
1153 : !> \param num_integ_points ...
1154 : ! **************************************************************************************************
1155 16 : SUBROUTINE Wc_to_Minv_Wc_Minv(fm_mat_W, fm_matrix_Minv, para_env, dimen_RI, num_integ_points)
1156 : TYPE(cp_fm_type), DIMENSION(:) :: fm_mat_W
1157 : TYPE(cp_fm_type), DIMENSION(:, :) :: fm_matrix_Minv
1158 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
1159 : INTEGER :: dimen_RI, num_integ_points
1160 :
1161 : CHARACTER(LEN=*), PARAMETER :: routineN = 'Wc_to_Minv_Wc_Minv'
1162 :
1163 : INTEGER :: handle, jquad
1164 : TYPE(cp_fm_type) :: fm_work_Minv, fm_work_Minv_W
1165 :
1166 16 : CALL timeset(routineN, handle)
1167 :
1168 16 : CALL cp_fm_create(fm_work_Minv, fm_mat_W(1)%matrix_struct)
1169 16 : CALL cp_fm_copy_general(fm_matrix_Minv(1, 1), fm_work_Minv, para_env)
1170 :
1171 16 : CALL cp_fm_create(fm_work_Minv_W, fm_mat_W(1)%matrix_struct)
1172 :
1173 112 : DO jquad = 1, num_integ_points
1174 :
1175 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_work_Minv, fm_mat_W(jquad), &
1176 96 : 0.0_dp, fm_work_Minv_W)
1177 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_work_Minv_W, fm_work_Minv, &
1178 112 : 0.0_dp, fm_mat_W(jquad))
1179 :
1180 : END DO
1181 :
1182 16 : CALL cp_fm_release(fm_work_Minv)
1183 :
1184 16 : CALL cp_fm_release(fm_work_Minv_W)
1185 :
1186 16 : CALL timestop(handle)
1187 :
1188 16 : END SUBROUTINE Wc_to_Minv_Wc_Minv
1189 :
1190 : ! **************************************************************************************************
1191 : !> \brief ...
1192 : !> \param qs_env ...
1193 : !> \param wkp_W ...
1194 : !> \param wkp_V ...
1195 : !> \param kpoints ...
1196 : !> \param h_inv ...
1197 : !> \param periodic ...
1198 : ! **************************************************************************************************
1199 22 : SUBROUTINE compute_wkp_W(qs_env, wkp_W, wkp_V, kpoints, h_inv, periodic)
1200 :
1201 : TYPE(qs_environment_type), POINTER :: qs_env
1202 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1203 : INTENT(OUT) :: wkp_W, wkp_V
1204 : TYPE(kpoint_type), POINTER :: kpoints
1205 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv
1206 : INTEGER, DIMENSION(3) :: periodic
1207 :
1208 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_wkp_W'
1209 :
1210 : INTEGER :: handle, i_x, ikp, info, j_y, k_z, &
1211 : kpoint_weights_W_method, n_x, n_y, &
1212 : n_z, nkp, nsuperfine, num_lin_eqs
1213 : REAL(KIND=dp) :: exp_kpoints, integral, k_sq, weight
1214 : REAL(KIND=dp), DIMENSION(3) :: k_vec, x_vec
1215 22 : REAL(KIND=dp), DIMENSION(:), POINTER :: right_side, wkp, wkp_tmp
1216 22 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix_lin_eqs, xkp
1217 :
1218 22 : CALL timeset(routineN, handle)
1219 :
1220 22 : kpoint_weights_W_method = qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method
1221 :
1222 22 : CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp)
1223 :
1224 : ! we determine the kpoint weights of the Monkhors Pack mesh new
1225 : ! such that the functions 1/k^2, 1/k and const are integrated exactly
1226 : ! in the Brillouin zone
1227 : ! this is done by minimizing sum_i |w_i|^2 where w_i are the weights of
1228 : ! the i-th kpoint under the following constraints:
1229 : ! 1) 1/k^2, 1/k and const are integrated exactly
1230 : ! 2) the kpoint weights of kpoints with identical absolute value are
1231 : ! the same, of e.g. (1/8,3/8,3/8) same weight as for (3/8,1/8,3/8)
1232 : ! for 1d and 2d materials: we use ordinary Monkhorst-Pack weights, checked
1233 : ! by SUM(periodic) == 3
1234 88 : ALLOCATE (wkp_V(nkp), wkp_W(nkp))
1235 :
1236 : ! for exchange part of self-energy, we use truncated Coulomb operator that should be fine
1237 : ! with uniform weights (without k-point extrapolation)
1238 22 : IF (ALLOCATED(qs_env%mp2_env%ri_rpa_im_time%wkp_V)) THEN
1239 432 : wkp_V(:) = qs_env%mp2_env%ri_rpa_im_time%wkp_V(:)
1240 : ELSE
1241 18 : wkp_V(:) = wkp(:)
1242 : END IF
1243 :
1244 22 : IF (kpoint_weights_W_method == kp_weights_W_uniform) THEN
1245 :
1246 : ! in the k-point weights wkp, there might be k-point extrapolation included
1247 450 : wkp_W(:) = wkp(:)
1248 :
1249 0 : ELSE IF (kpoint_weights_W_method == kp_weights_W_tailored .OR. &
1250 : kpoint_weights_W_method == kp_weights_W_auto) THEN
1251 :
1252 0 : IF (kpoint_weights_W_method == kp_weights_W_tailored) THEN
1253 0 : exp_kpoints = qs_env%mp2_env%ri_rpa_im_time%exp_tailored_weights
1254 : END IF
1255 :
1256 0 : IF (kpoint_weights_W_method == kp_weights_W_auto) THEN
1257 0 : IF (SUM(periodic) == 2) exp_kpoints = -1.0_dp
1258 : END IF
1259 :
1260 : ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid
1261 0 : nsuperfine = 500
1262 0 : integral = 0.0_dp
1263 :
1264 0 : IF (periodic(1) == 1) THEN
1265 : n_x = nsuperfine
1266 : ELSE
1267 0 : n_x = 1
1268 : END IF
1269 0 : IF (periodic(2) == 1) THEN
1270 : n_y = nsuperfine
1271 : ELSE
1272 0 : n_y = 1
1273 : END IF
1274 0 : IF (periodic(3) == 1) THEN
1275 : n_z = nsuperfine
1276 : ELSE
1277 0 : n_z = 1
1278 : END IF
1279 :
1280 : ! actually, there is the factor *det_3x3(h_inv) missing to account for the
1281 : ! integration volume but for wkp det_3x3(h_inv) is needed
1282 0 : weight = 1.0_dp/(REAL(n_x, dp)*REAL(n_y, dp)*REAL(n_z, dp))
1283 0 : DO i_x = 1, n_x
1284 0 : DO j_y = 1, n_y
1285 0 : DO k_z = 1, n_z
1286 :
1287 0 : IF (periodic(1) == 1) THEN
1288 0 : x_vec(1) = (REAL(i_x - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp)
1289 : ELSE
1290 0 : x_vec(1) = 0.0_dp
1291 : END IF
1292 0 : IF (periodic(2) == 1) THEN
1293 0 : x_vec(2) = (REAL(j_y - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp)
1294 : ELSE
1295 0 : x_vec(2) = 0.0_dp
1296 : END IF
1297 0 : IF (periodic(3) == 1) THEN
1298 0 : x_vec(3) = (REAL(k_z - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp)
1299 : ELSE
1300 0 : x_vec(3) = 0.0_dp
1301 : END IF
1302 :
1303 0 : k_vec = MATMUL(h_inv(1:3, 1:3), x_vec)
1304 0 : k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2
1305 0 : integral = integral + weight*k_sq**(exp_kpoints*0.5_dp)
1306 :
1307 : END DO
1308 : END DO
1309 : END DO
1310 :
1311 0 : num_lin_eqs = nkp + 2
1312 :
1313 0 : ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs))
1314 0 : matrix_lin_eqs(:, :) = 0.0_dp
1315 :
1316 0 : DO ikp = 1, nkp
1317 :
1318 0 : k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp))
1319 0 : k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2
1320 :
1321 0 : matrix_lin_eqs(ikp, ikp) = 2.0_dp
1322 0 : matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp
1323 0 : matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp
1324 :
1325 0 : matrix_lin_eqs(ikp, nkp + 2) = k_sq**(exp_kpoints*0.5_dp)
1326 0 : matrix_lin_eqs(nkp + 2, ikp) = k_sq**(exp_kpoints*0.5_dp)
1327 :
1328 : END DO
1329 :
1330 0 : CALL invmat(matrix_lin_eqs, info)
1331 : ! check whether inversion was successful
1332 0 : CPASSERT(info == 0)
1333 :
1334 0 : ALLOCATE (right_side(num_lin_eqs))
1335 0 : right_side = 0.0_dp
1336 0 : right_side(nkp + 1) = 1.0_dp
1337 : ! divide integral by two because CP2K k-mesh already considers symmetry k <-> -k
1338 0 : right_side(nkp + 2) = integral
1339 :
1340 0 : ALLOCATE (wkp_tmp(num_lin_eqs))
1341 :
1342 0 : wkp_tmp(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side)
1343 :
1344 0 : wkp_W(1:nkp) = wkp_tmp(1:nkp)
1345 :
1346 0 : DEALLOCATE (matrix_lin_eqs, right_side, wkp_tmp)
1347 :
1348 : END IF
1349 :
1350 22 : CALL timestop(handle)
1351 :
1352 22 : END SUBROUTINE compute_wkp_W
1353 :
1354 : ! **************************************************************************************************
1355 : !> \brief ...
1356 : !> \param qs_env ...
1357 : !> \param Eigenval_kp ...
1358 : ! **************************************************************************************************
1359 16 : SUBROUTINE get_bandstruc_and_k_dependent_MOs(qs_env, Eigenval_kp)
1360 : TYPE(qs_environment_type), POINTER :: qs_env
1361 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_kp
1362 :
1363 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_bandstruc_and_k_dependent_MOs'
1364 :
1365 : INTEGER :: handle, ikp, ispin, nmo, nspins
1366 : INTEGER, DIMENSION(3) :: nkp_grid_G
1367 16 : REAL(KIND=dp), DIMENSION(:), POINTER :: ev
1368 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kpgeneral
1369 : TYPE(kpoint_type), POINTER :: kpoints_Sigma
1370 : TYPE(mp_para_env_type), POINTER :: para_env
1371 :
1372 16 : CALL timeset(routineN, handle)
1373 :
1374 : NULLIFY (qs_env%mp2_env%ri_rpa_im_time%kpoints_G, &
1375 16 : qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, &
1376 16 : qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc, &
1377 16 : para_env)
1378 :
1379 16 : nkp_grid_G(1:3) = [1, 1, 1]
1380 :
1381 16 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
1382 :
1383 : CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_G, &
1384 : "MONKHORST-PACK", para_env%num_pe, &
1385 16 : mp_grid=nkp_grid_G(1:3))
1386 :
1387 16 : IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
1388 :
1389 : ! set up k-points for GW band structure calculation, will be completed later
1390 16 : CALL get_kpgeneral_for_Sigma_kpoints(qs_env, kpgeneral)
1391 :
1392 : CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, &
1393 : "GENERAL", para_env%num_pe, &
1394 16 : kpgeneral=kpgeneral)
1395 :
1396 : CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc, &
1397 : "GENERAL", para_env%num_pe, &
1398 16 : kpgeneral=kpgeneral, with_xc_terms=.FALSE.)
1399 :
1400 16 : kpoints_Sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
1401 16 : nmo = SIZE(Eigenval_kp, 1)
1402 16 : nspins = SIZE(Eigenval_kp, 3)
1403 :
1404 48 : ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%Eigenval_Gamma(nmo))
1405 340 : qs_env%mp2_env%ri_rpa_im_time%Eigenval_Gamma(:) = Eigenval_kp(:, 1, 1)
1406 :
1407 16 : DEALLOCATE (Eigenval_kp)
1408 :
1409 80 : ALLOCATE (Eigenval_kp(nmo, kpoints_Sigma%nkp, nspins))
1410 :
1411 136 : DO ikp = 1, kpoints_Sigma%nkp
1412 :
1413 272 : DO ispin = 1, nspins
1414 :
1415 136 : ev => kpoints_Sigma%kp_env(ikp)%kpoint_env%mos(1, ispin)%eigenvalues
1416 :
1417 3016 : Eigenval_kp(:, ikp, ispin) = ev(:)
1418 :
1419 : END DO
1420 :
1421 : END DO
1422 :
1423 16 : DEALLOCATE (kpgeneral)
1424 :
1425 : END IF
1426 :
1427 16 : CALL release_hfx_stuff(qs_env)
1428 :
1429 16 : CALL timestop(handle)
1430 :
1431 16 : END SUBROUTINE get_bandstruc_and_k_dependent_MOs
1432 :
1433 : ! **************************************************************************************************
1434 : !> \brief releases part of the given qs_env in order to save memory
1435 : !> \param qs_env the object to release
1436 : ! **************************************************************************************************
1437 16 : SUBROUTINE release_hfx_stuff(qs_env)
1438 : TYPE(qs_environment_type), POINTER :: qs_env
1439 :
1440 16 : IF (ASSOCIATED(qs_env%x_data) .AND. .NOT. qs_env%mp2_env%ri_g0w0%do_ri_Sigma_x) THEN
1441 2 : CALL hfx_release(qs_env%x_data)
1442 : END IF
1443 :
1444 16 : END SUBROUTINE release_hfx_stuff
1445 :
1446 : ! **************************************************************************************************
1447 : !> \brief ...
1448 : !> \param qs_env ...
1449 : !> \param kpoints ...
1450 : !> \param scheme ...
1451 : !> \param group_size_ext ...
1452 : !> \param mp_grid ...
1453 : !> \param kpgeneral ...
1454 : !> \param with_xc_terms ...
1455 : !> \param kp_shift ...
1456 : !> \param gamma_centered ...
1457 : ! **************************************************************************************************
1458 336 : SUBROUTINE create_kp_and_calc_kp_orbitals(qs_env, kpoints, scheme, &
1459 48 : group_size_ext, mp_grid, kpgeneral, with_xc_terms, &
1460 : kp_shift, gamma_centered)
1461 :
1462 : TYPE(qs_environment_type), POINTER :: qs_env
1463 : TYPE(kpoint_type), POINTER :: kpoints
1464 : CHARACTER(LEN=*), INTENT(IN) :: scheme
1465 : INTEGER :: group_size_ext
1466 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: mp_grid
1467 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
1468 : OPTIONAL :: kpgeneral
1469 : LOGICAL, OPTIONAL :: with_xc_terms
1470 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: kp_shift
1471 : LOGICAL, INTENT(IN), OPTIONAL :: gamma_centered
1472 :
1473 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_kp_and_calc_kp_orbitals'
1474 :
1475 : INTEGER :: handle, i_dim, i_re_im, ikp, ispin, nkp, &
1476 : nspins
1477 : INTEGER, DIMENSION(3) :: cell_grid, periodic
1478 : LOGICAL :: my_with_xc_terms
1479 48 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1480 : TYPE(cell_type), POINTER :: cell
1481 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1482 : TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
1483 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1484 : TYPE(cp_fm_type) :: fm_work
1485 : TYPE(cp_fm_type), POINTER :: imos, rmos
1486 48 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_desymm
1487 48 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_ks_kp, mat_s_kp
1488 : TYPE(dft_control_type), POINTER :: dft_control
1489 : TYPE(kpoint_env_type), POINTER :: kp
1490 : TYPE(mp_para_env_type), POINTER :: para_env
1491 : TYPE(qs_scf_env_type), POINTER :: scf_env
1492 : TYPE(scf_control_type), POINTER :: scf_control
1493 :
1494 48 : CALL timeset(routineN, handle)
1495 :
1496 48 : my_with_xc_terms = .TRUE.
1497 48 : IF (PRESENT(with_xc_terms)) my_with_xc_terms = with_xc_terms
1498 :
1499 : CALL get_qs_env(qs_env, &
1500 : para_env=para_env, &
1501 : blacs_env=blacs_env, &
1502 : matrix_s=matrix_s, &
1503 : scf_env=scf_env, &
1504 : scf_control=scf_control, &
1505 48 : cell=cell)
1506 :
1507 : ! get kpoints
1508 : CALL calculate_kpoints_for_bs(kpoints, scheme, kpgeneral=kpgeneral, mp_grid=mp_grid, &
1509 : group_size_ext=group_size_ext, kp_shift=kp_shift, &
1510 64 : gamma_centered=gamma_centered)
1511 :
1512 48 : CALL kpoint_env_initialize(kpoints, para_env, blacs_env)
1513 :
1514 : ! calculate all MOs that are accessible in the given
1515 : ! Gaussian AO basis, therefore nadd=1E10
1516 48 : CALL kpoint_initialize_mos(kpoints, qs_env%mos, 2000000000)
1517 48 : CALL kpoint_initialize_mo_set(kpoints)
1518 :
1519 48 : CALL get_cell(cell=cell, periodic=periodic)
1520 :
1521 192 : DO i_dim = 1, 3
1522 : ! we have at most 3 neigboring cells per dimension and at least one because
1523 : ! the density response at Gamma is only divided to neighboring
1524 192 : IF (periodic(i_dim) == 1) THEN
1525 96 : cell_grid(i_dim) = MAX(MIN((kpoints%nkp_grid(i_dim)/2)*2 - 1, 1), 3)
1526 : ELSE
1527 48 : cell_grid(i_dim) = 1
1528 : END IF
1529 : END DO
1530 48 : CALL init_cell_index_rpa(cell_grid, kpoints%cell_to_index, kpoints%index_to_cell, cell)
1531 :
1532 : ! get S(k)
1533 48 : CALL get_qs_env(qs_env, matrix_s=matrix_s, scf_env=scf_env, scf_control=scf_control, dft_control=dft_control)
1534 :
1535 48 : NULLIFY (matrix_s_desymm)
1536 48 : CALL dbcsr_allocate_matrix_set(matrix_s_desymm, 1)
1537 48 : ALLOCATE (matrix_s_desymm(1)%matrix)
1538 : CALL dbcsr_create(matrix=matrix_s_desymm(1)%matrix, template=matrix_s(1)%matrix, &
1539 48 : matrix_type=dbcsr_type_no_symmetry)
1540 48 : CALL dbcsr_desymmetrize(matrix_s(1)%matrix, matrix_s_desymm(1)%matrix)
1541 :
1542 48 : CALL mat_kp_from_mat_gamma(qs_env, mat_s_kp, matrix_s_desymm(1)%matrix, kpoints, 1)
1543 :
1544 48 : CALL get_kpoint_info(kpoints, nkp=nkp)
1545 :
1546 48 : matrix_struct => kpoints%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct
1547 :
1548 48 : CALL cp_cfm_create(cksmat, matrix_struct)
1549 48 : CALL cp_cfm_create(csmat, matrix_struct)
1550 48 : CALL cp_cfm_create(cmos, matrix_struct)
1551 48 : CALL cp_cfm_create(cwork, matrix_struct)
1552 48 : CALL cp_fm_create(fm_work, matrix_struct)
1553 :
1554 48 : nspins = dft_control%nspins
1555 :
1556 102 : DO ispin = 1, nspins
1557 :
1558 : ! get H(k)
1559 54 : IF (my_with_xc_terms) THEN
1560 36 : CALL mat_kp_from_mat_gamma(qs_env, mat_ks_kp, qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix, kpoints, ispin)
1561 : ELSE
1562 : CALL mat_kp_from_mat_gamma(qs_env, mat_ks_kp, qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix, &
1563 18 : kpoints, ispin)
1564 : END IF
1565 :
1566 392 : DO ikp = 1, nkp
1567 :
1568 290 : CALL copy_dbcsr_to_fm(mat_ks_kp(ikp, 1)%matrix, kpoints%kp_env(ikp)%kpoint_env%wmat(1, ispin))
1569 290 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, kpoints%kp_env(ikp)%kpoint_env%wmat(1, ispin))
1570 :
1571 290 : CALL copy_dbcsr_to_fm(mat_ks_kp(ikp, 2)%matrix, kpoints%kp_env(ikp)%kpoint_env%wmat(2, ispin))
1572 290 : CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, kpoints%kp_env(ikp)%kpoint_env%wmat(2, ispin))
1573 :
1574 290 : CALL copy_dbcsr_to_fm(mat_s_kp(ikp, 1)%matrix, fm_work)
1575 290 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fm_work)
1576 :
1577 290 : CALL copy_dbcsr_to_fm(mat_s_kp(ikp, 2)%matrix, fm_work)
1578 290 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fm_work)
1579 :
1580 290 : kp => kpoints%kp_env(ikp)%kpoint_env
1581 :
1582 290 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
1583 290 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1584 :
1585 290 : IF (scf_env%cholesky_method == cholesky_off .OR. &
1586 : qs_env%mp2_env%ri_rpa_im_time%make_overlap_mat_ao_pos_definite) THEN
1587 0 : CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, scf_control%eps_eigval)
1588 : ELSE
1589 290 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
1590 : END IF
1591 :
1592 290 : CALL cp_cfm_to_fm(cmos, rmos, imos)
1593 :
1594 12382 : kp%mos(2, ispin)%eigenvalues = eigenvalues
1595 :
1596 : END DO
1597 :
1598 : END DO
1599 :
1600 304 : DO ikp = 1, nkp
1601 816 : DO i_re_im = 1, 2
1602 768 : CALL dbcsr_deallocate_matrix(mat_ks_kp(ikp, i_re_im)%matrix)
1603 : END DO
1604 : END DO
1605 48 : DEALLOCATE (mat_ks_kp)
1606 :
1607 304 : DO ikp = 1, nkp
1608 816 : DO i_re_im = 1, 2
1609 768 : CALL dbcsr_deallocate_matrix(mat_s_kp(ikp, i_re_im)%matrix)
1610 : END DO
1611 : END DO
1612 48 : DEALLOCATE (mat_s_kp)
1613 :
1614 48 : CALL dbcsr_deallocate_matrix(matrix_s_desymm(1)%matrix)
1615 48 : DEALLOCATE (matrix_s_desymm)
1616 :
1617 48 : CALL cp_cfm_release(cksmat)
1618 48 : CALL cp_cfm_release(csmat)
1619 48 : CALL cp_cfm_release(cwork)
1620 48 : CALL cp_cfm_release(cmos)
1621 48 : CALL cp_fm_release(fm_work)
1622 :
1623 48 : CALL timestop(handle)
1624 :
1625 48 : END SUBROUTINE create_kp_and_calc_kp_orbitals
1626 :
1627 : ! **************************************************************************************************
1628 : !> \brief ...
1629 : !> \param qs_env ...
1630 : !> \param mat_kp ...
1631 : !> \param mat_gamma ...
1632 : !> \param kpoints ...
1633 : !> \param ispin ...
1634 : !> \param real_mat_real_space ...
1635 : ! **************************************************************************************************
1636 114 : SUBROUTINE mat_kp_from_mat_gamma(qs_env, mat_kp, mat_gamma, kpoints, ispin, real_mat_real_space)
1637 :
1638 : TYPE(qs_environment_type), POINTER :: qs_env
1639 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_kp
1640 : TYPE(dbcsr_type) :: mat_gamma
1641 : TYPE(kpoint_type), POINTER :: kpoints
1642 : INTEGER :: ispin
1643 : LOGICAL, INTENT(IN), OPTIONAL :: real_mat_real_space
1644 :
1645 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mat_kp_from_mat_gamma'
1646 :
1647 : INTEGER :: handle, i_cell, i_re_im, ikp, nkp, &
1648 : num_cells
1649 : INTEGER, DIMENSION(3) :: periodic
1650 114 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1651 114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1652 : TYPE(cell_type), POINTER :: cell
1653 114 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_real_space
1654 :
1655 114 : CALL timeset(routineN, handle)
1656 :
1657 114 : CALL get_qs_env(qs_env, cell=cell)
1658 114 : CALL get_cell(cell=cell, periodic=periodic)
1659 114 : num_cells = 3**(periodic(1) + periodic(2) + periodic(3))
1660 :
1661 114 : NULLIFY (mat_real_space)
1662 114 : CALL dbcsr_allocate_matrix_set(mat_real_space, num_cells)
1663 1140 : DO i_cell = 1, num_cells
1664 1026 : ALLOCATE (mat_real_space(i_cell)%matrix)
1665 : CALL dbcsr_create(matrix=mat_real_space(i_cell)%matrix, &
1666 1026 : template=mat_gamma)
1667 1026 : CALL dbcsr_reserve_all_blocks(mat_real_space(i_cell)%matrix)
1668 1140 : CALL dbcsr_set(mat_real_space(i_cell)%matrix, 0.0_dp)
1669 : END DO
1670 :
1671 114 : CALL dbcsr_copy(mat_real_space(1)%matrix, mat_gamma)
1672 :
1673 114 : CALL get_mat_cell_T_from_mat_gamma(mat_real_space, qs_env, kpoints, 2, 0)
1674 :
1675 114 : NULLIFY (xkp, cell_to_index)
1676 114 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, cell_to_index=cell_to_index)
1677 :
1678 114 : IF (ispin == 1) THEN
1679 108 : NULLIFY (mat_kp)
1680 108 : CALL dbcsr_allocate_matrix_set(mat_kp, nkp, 2)
1681 668 : DO ikp = 1, nkp
1682 1788 : DO i_re_im = 1, 2
1683 1120 : ALLOCATE (mat_kp(ikp, i_re_im)%matrix)
1684 1120 : CALL dbcsr_create(matrix=mat_kp(ikp, i_re_im)%matrix, template=mat_gamma)
1685 1120 : CALL dbcsr_reserve_all_blocks(mat_kp(ikp, i_re_im)%matrix)
1686 1680 : CALL dbcsr_set(mat_kp(ikp, i_re_im)%matrix, 0.0_dp)
1687 : END DO
1688 : END DO
1689 : END IF
1690 :
1691 114 : IF (PRESENT(real_mat_real_space)) THEN
1692 : CALL real_space_to_kpoint_transform_rpa(mat_kp(:, 1), mat_kp(:, 2), mat_real_space, kpoints, 0.0_dp, &
1693 12 : real_mat_real_space)
1694 : ELSE
1695 102 : CALL real_space_to_kpoint_transform_rpa(mat_kp(:, 1), mat_kp(:, 2), mat_real_space, kpoints, 0.0_dp)
1696 : END IF
1697 :
1698 1140 : DO i_cell = 1, num_cells
1699 1140 : CALL dbcsr_deallocate_matrix(mat_real_space(i_cell)%matrix)
1700 : END DO
1701 114 : DEALLOCATE (mat_real_space)
1702 :
1703 114 : CALL timestop(handle)
1704 :
1705 114 : END SUBROUTINE mat_kp_from_mat_gamma
1706 :
1707 : ! **************************************************************************************************
1708 : !> \brief ...
1709 : !> \param qs_env ...
1710 : !> \param kpgeneral ...
1711 : ! **************************************************************************************************
1712 16 : SUBROUTINE get_kpgeneral_for_Sigma_kpoints(qs_env, kpgeneral)
1713 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1714 : REAL(kind=dp), DIMENSION(:, :), POINTER :: kpgeneral
1715 :
1716 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_kpgeneral_for_Sigma_kpoints'
1717 :
1718 : INTEGER :: handle, i_kp_in_kp_line, i_special_kp, &
1719 : i_x, ikk, j_y, k_z, n_kp_in_kp_line, &
1720 : n_special_kp
1721 16 : INTEGER, DIMENSION(:), POINTER :: nkp_grid
1722 :
1723 16 : CALL timeset(routineN, handle)
1724 :
1725 16 : n_special_kp = qs_env%mp2_env%ri_g0w0%n_special_kp
1726 16 : n_kp_in_kp_line = qs_env%mp2_env%ri_g0w0%n_kp_in_kp_line
1727 16 : IF (n_special_kp > 0) THEN
1728 14 : qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp = n_kp_in_kp_line*(n_special_kp - 1) + 1
1729 : ELSE
1730 2 : qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp = 0
1731 : END IF
1732 :
1733 : qs_env%mp2_env%ri_g0w0%nkp_self_energy_monkh_pack = qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(1)* &
1734 : qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(2)* &
1735 16 : qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(3)
1736 :
1737 : qs_env%mp2_env%ri_g0w0%nkp_self_energy = qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp + &
1738 16 : qs_env%mp2_env%ri_g0w0%nkp_self_energy_monkh_pack
1739 :
1740 48 : ALLOCATE (kpgeneral(3, qs_env%mp2_env%ri_g0w0%nkp_self_energy))
1741 :
1742 16 : IF (n_special_kp > 0) THEN
1743 :
1744 112 : kpgeneral(1:3, 1) = qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, 1)
1745 :
1746 14 : ikk = 1
1747 :
1748 28 : DO i_special_kp = 2, n_special_kp
1749 70 : DO i_kp_in_kp_line = 1, n_kp_in_kp_line
1750 :
1751 42 : ikk = ikk + 1
1752 : kpgeneral(1:3, ikk) = qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp - 1) + &
1753 : REAL(i_kp_in_kp_line, KIND=dp)/REAL(n_kp_in_kp_line, KIND=dp)* &
1754 : (qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp) - &
1755 350 : qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp - 1))
1756 :
1757 : END DO
1758 : END DO
1759 :
1760 : ELSE
1761 :
1762 : ikk = 0
1763 :
1764 : END IF
1765 :
1766 16 : nkp_grid => qs_env%mp2_env%ri_g0w0%kp_grid_Sigma
1767 :
1768 48 : DO i_x = 1, nkp_grid(1)
1769 112 : DO j_y = 1, nkp_grid(2)
1770 160 : DO k_z = 1, nkp_grid(3)
1771 64 : ikk = ikk + 1
1772 64 : kpgeneral(1, ikk) = REAL(2*i_x - nkp_grid(1) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(1), KIND=dp))
1773 64 : kpgeneral(2, ikk) = REAL(2*j_y - nkp_grid(2) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(2), KIND=dp))
1774 128 : kpgeneral(3, ikk) = REAL(2*k_z - nkp_grid(3) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(3), KIND=dp))
1775 : END DO
1776 : END DO
1777 : END DO
1778 :
1779 16 : CALL timestop(handle)
1780 :
1781 16 : END SUBROUTINE get_kpgeneral_for_Sigma_kpoints
1782 :
1783 0 : END MODULE rpa_gw_kpoints_util
|
