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 Utility functions for RPA calculations
10 : !> \par History
11 : !> 06.2019 Moved from rpa_ri_gpw.F [Frederick Stein]
12 : ! **************************************************************************************************
13 : MODULE rpa_util
14 :
15 : USE cell_types, ONLY: cell_type,&
16 : get_cell
17 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
18 : cp_blacs_env_release,&
19 : cp_blacs_env_type
20 : USE cp_cfm_types, ONLY: cp_cfm_copy_general,&
21 : cp_cfm_create,&
22 : cp_cfm_release,&
23 : cp_cfm_set_all,&
24 : cp_cfm_type,&
25 : cp_fm_to_cfm
26 : USE cp_dbcsr_api, ONLY: &
27 : dbcsr_create, dbcsr_deallocate_matrix, dbcsr_filter, dbcsr_get_info, dbcsr_multiply, &
28 : dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
29 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
30 : copy_fm_to_dbcsr,&
31 : dbcsr_allocate_matrix_set,&
32 : dbcsr_deallocate_matrix_set
33 : USE cp_fm_basic_linalg, ONLY: cp_fm_syrk,&
34 : cp_fm_transpose
35 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
36 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37 : cp_fm_struct_release,&
38 : cp_fm_struct_type
39 : USE cp_fm_types, ONLY: &
40 : cp_fm_copy_general, cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_release, &
41 : cp_fm_set_all, cp_fm_to_fm, cp_fm_to_fm_submat_general, cp_fm_type
42 : USE dbt_api, ONLY: dbt_destroy,&
43 : dbt_type
44 : USE dgemm_counter_types, ONLY: dgemm_counter_start,&
45 : dgemm_counter_stop,&
46 : dgemm_counter_type
47 : USE hfx_types, ONLY: block_ind_type,&
48 : dealloc_containers,&
49 : hfx_compression_type
50 : USE input_constants, ONLY: wfc_mm_style_gemm,&
51 : wfc_mm_style_syrk
52 : USE kinds, ONLY: dp
53 : USE kpoint_types, ONLY: get_kpoint_info,&
54 : kpoint_release,&
55 : kpoint_type
56 : USE mathconstants, ONLY: z_zero
57 : USE message_passing, ONLY: mp_para_env_release,&
58 : mp_para_env_type
59 : USE mp2_laplace, ONLY: calc_fm_mat_S_laplace
60 : USE parallel_gemm_api, ONLY: parallel_gemm
61 : USE qs_environment_types, ONLY: get_qs_env,&
62 : qs_environment_type
63 : USE rpa_gw_kpoints_util, ONLY: compute_wkp_W
64 : #include "./base/base_uses.f90"
65 :
66 : IMPLICIT NONE
67 :
68 : PRIVATE
69 :
70 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_util'
71 :
72 : PUBLIC :: compute_Erpa_by_freq_int, alloc_im_time, calc_mat_Q, &
73 : dealloc_im_time, contract_P_omega_with_mat_L, calc_fm_mat_S_rpa, remove_scaling_factor_rpa
74 :
75 : CONTAINS
76 :
77 : ! **************************************************************************************************
78 : !> \brief ...
79 : !> \param qs_env ...
80 : !> \param para_env ...
81 : !> \param dimen_RI ...
82 : !> \param dimen_RI_red ...
83 : !> \param num_integ_points ...
84 : !> \param nspins ...
85 : !> \param create_cfm_mo_coeff whether to promote the real Gamma-point coefficients to CFM
86 : !> \param fm_mat_Q ...
87 : !> \param cfm_mo_coeff ...
88 : !> \param cfm_matrix_Minv_L_kpoints ...
89 : !> \param cfm_matrix_L_kpoints ...
90 : !> \param fm_matrix_Minv_L_gamma ...
91 : !> \param mat_P_global ...
92 : !> \param t_3c_O ...
93 : !> \param matrix_s ...
94 : !> \param kpoints ...
95 : !> \param eps_filter_im_time ...
96 : !> \param cut_memory ...
97 : !> \param nkp ...
98 : !> \param num_cells_dm ...
99 : !> \param num_3c_repl ...
100 : !> \param size_P ...
101 : !> \param ikp_local ...
102 : !> \param index_to_cell_3c ...
103 : !> \param cell_to_index_3c ...
104 : !> \param col_blk_size ...
105 : !> \param do_ic_model ...
106 : !> \param do_kpoints_cubic_RPA ...
107 : !> \param do_kpoints_from_Gamma ...
108 : !> \param do_ri_Sigma_x ...
109 : !> \param my_open_shell ...
110 : !> \param has_mat_P_blocks ...
111 : !> \param wkp_W ...
112 : !> \param cfm_mat_Q ...
113 : !> \param cfm_mat_Minv_L_kpoints ...
114 : !> \param cfm_mat_L_kpoints ...
115 : !> \param fm_mat_Minv_L_gamma ...
116 : !> \param cfm_mat_RI_global_work ...
117 : !> \param fm_mat_work ...
118 : !> \param mat_dm ...
119 : !> \param mat_L ...
120 : !> \param mat_M_P_munu_occ ...
121 : !> \param mat_M_P_munu_virt ...
122 : !> \param mat_MinvVMinv ...
123 : !> \param mat_P_omega ...
124 : !> \param mat_P_omega_kp ...
125 : !> \param mat_work ...
126 : !> \param mo_coeff ...
127 : ! **************************************************************************************************
128 288 : SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_points, nspins, &
129 : create_cfm_mo_coeff, &
130 : fm_mat_Q, cfm_mo_coeff, &
131 : cfm_matrix_Minv_L_kpoints, cfm_matrix_L_kpoints, fm_matrix_Minv_L_gamma, mat_P_global, &
132 : t_3c_O, matrix_s, kpoints, eps_filter_im_time, &
133 : cut_memory, nkp, num_cells_dm, num_3c_repl, &
134 : size_P, ikp_local, &
135 : index_to_cell_3c, &
136 : cell_to_index_3c, &
137 : col_blk_size, &
138 : do_ic_model, do_kpoints_cubic_RPA, &
139 : do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, &
140 : has_mat_P_blocks, wkp_W, &
141 : cfm_mat_Q, cfm_mat_Minv_L_kpoints, cfm_mat_L_kpoints, fm_mat_Minv_L_gamma, &
142 : cfm_mat_RI_global_work, fm_mat_work, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, &
143 : mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
144 144 : mat_work, mo_coeff)
145 :
146 : TYPE(qs_environment_type), POINTER :: qs_env
147 : TYPE(mp_para_env_type), POINTER :: para_env
148 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, &
149 : num_integ_points, nspins
150 : LOGICAL, INTENT(IN) :: create_cfm_mo_coeff
151 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
152 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_mo_coeff, cfm_matrix_Minv_L_kpoints, &
153 : cfm_matrix_L_kpoints
154 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv_L_gamma
155 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_P_global
156 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
157 : INTENT(INOUT) :: t_3c_O
158 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
159 : TYPE(kpoint_type), POINTER :: kpoints
160 : REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time
161 : INTEGER, INTENT(IN) :: cut_memory
162 : INTEGER, INTENT(OUT) :: nkp, num_cells_dm, num_3c_repl, size_P, &
163 : ikp_local
164 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell_3c
165 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
166 : INTENT(OUT) :: cell_to_index_3c
167 : INTEGER, DIMENSION(:), POINTER :: col_blk_size
168 : LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA, &
169 : do_kpoints_from_Gamma, do_ri_Sigma_x, &
170 : my_open_shell
171 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), &
172 : INTENT(OUT) :: has_mat_P_blocks
173 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
174 : INTENT(OUT) :: wkp_W
175 : TYPE(cp_cfm_type), INTENT(OUT) :: cfm_mat_Q
176 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_mat_Minv_L_kpoints, cfm_mat_L_kpoints
177 : TYPE(cp_fm_type), INTENT(OUT) :: fm_mat_Minv_L_gamma
178 : TYPE(cp_cfm_type), INTENT(OUT) :: cfm_mat_RI_global_work
179 : TYPE(cp_fm_type), INTENT(OUT) :: fm_mat_work
180 : TYPE(dbcsr_p_type), INTENT(OUT) :: mat_dm, mat_L, mat_M_P_munu_occ, &
181 : mat_M_P_munu_virt, mat_MinvVMinv
182 : TYPE(dbcsr_p_type), ALLOCATABLE, &
183 : DIMENSION(:, :, :) :: mat_P_omega
184 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_kp
185 : TYPE(dbcsr_type), POINTER :: mat_work
186 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
187 :
188 : CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_im_time'
189 :
190 : INTEGER :: cell_grid_dm(3), first_ikp_local, &
191 : handle, i_dim, i_kp, ispin, jquad, &
192 : nspins_P_omega, periodic(3)
193 144 : INTEGER, DIMENSION(:), POINTER :: row_blk_size
194 144 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wkp_V
195 : TYPE(cell_type), POINTER :: cell
196 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_sub_kp
197 :
198 144 : CALL timeset(routineN, handle)
199 :
200 144 : NULLIFY (mat_L%matrix)
201 :
202 144 : IF (create_cfm_mo_coeff) THEN
203 212 : ALLOCATE (cfm_mo_coeff(nspins))
204 112 : DO ispin = 1, SIZE(mo_coeff)
205 112 : CALL create_mo_coeff(cfm_mo_coeff(ispin), mo_coeff(ispin))
206 : END DO
207 : END IF
208 :
209 144 : num_3c_repl = SIZE(t_3c_O, 2)
210 :
211 144 : IF (do_kpoints_cubic_RPA) THEN
212 : ! we always use an odd number of image cells
213 : ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical
214 24 : DO i_dim = 1, 3
215 24 : cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
216 : END DO
217 6 : num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
218 18 : ALLOCATE (index_to_cell_3c(3, SIZE(kpoints%index_to_cell, 2)))
219 6 : CPASSERT(SIZE(kpoints%index_to_cell, 1) == 3)
220 126 : index_to_cell_3c(:, :) = kpoints%index_to_cell(:, :)
221 0 : ALLOCATE (cell_to_index_3c(LBOUND(kpoints%cell_to_index, 1):UBOUND(kpoints%cell_to_index, 1), &
222 : LBOUND(kpoints%cell_to_index, 2):UBOUND(kpoints%cell_to_index, 2), &
223 30 : LBOUND(kpoints%cell_to_index, 3):UBOUND(kpoints%cell_to_index, 3)))
224 96 : cell_to_index_3c(:, :, :) = kpoints%cell_to_index(:, :, :)
225 :
226 : ELSE
227 138 : ALLOCATE (index_to_cell_3c(3, 1))
228 552 : index_to_cell_3c(:, 1) = 0
229 138 : ALLOCATE (cell_to_index_3c(0:0, 0:0, 0:0))
230 138 : cell_to_index_3c(0, 0, 0) = 1
231 138 : num_cells_dm = 1
232 : END IF
233 :
234 144 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
235 :
236 22 : CALL get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_RI, ikp_local, first_ikp_local)
237 :
238 22 : CALL cp_cfm_create(cfm_mat_Q, fm_struct_sub_kp)
239 22 : CALL cp_cfm_set_all(cfm_mat_Q, z_zero)
240 : ELSE
241 122 : first_ikp_local = 1
242 : END IF
243 :
244 : ! if we do kpoints, mat_P has a kpoint and mat_P_omega has the inted
245 : ! mat_P(tau, kpoint)
246 144 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
247 :
248 22 : NULLIFY (cell)
249 22 : CALL get_qs_env(qs_env, cell=cell)
250 22 : CALL get_cell(cell=cell, periodic=periodic)
251 :
252 22 : CALL get_kpoint_info(kpoints, nkp=nkp)
253 : ! compute k-point weights such that functions 1/k^2, 1/k and const function are
254 : ! integrated correctly
255 22 : CALL compute_wkp_W(qs_env, wkp_W, wkp_V, kpoints, cell%h_inv, periodic)
256 22 : DEALLOCATE (wkp_V)
257 :
258 : ELSE
259 122 : nkp = 1
260 : END IF
261 :
262 144 : IF (do_kpoints_cubic_RPA) THEN
263 6 : size_P = MAX(num_cells_dm/2 + 1, nkp)
264 138 : ELSE IF (do_kpoints_from_Gamma) THEN
265 16 : size_P = MAX(3**(periodic(1) + periodic(2) + periodic(3)), nkp)
266 : ELSE
267 122 : size_P = 1
268 : END IF
269 :
270 144 : nspins_P_omega = 1
271 144 : IF (my_open_shell) nspins_P_omega = 2
272 :
273 5518 : ALLOCATE (mat_P_omega(num_integ_points, size_P, nspins_P_omega))
274 320 : DO ispin = 1, nspins_P_omega
275 952 : DO i_kp = 1, size_P
276 4798 : DO jquad = 1, num_integ_points
277 3990 : NULLIFY (mat_P_omega(jquad, i_kp, ispin)%matrix)
278 3990 : ALLOCATE (mat_P_omega(jquad, i_kp, ispin)%matrix)
279 : CALL dbcsr_create(matrix=mat_P_omega(jquad, i_kp, ispin)%matrix, &
280 3990 : template=mat_P_global%matrix)
281 4622 : CALL dbcsr_set(mat_P_omega(jquad, i_kp, ispin)%matrix, 0.0_dp)
282 : END DO
283 : END DO
284 : END DO
285 :
286 144 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
287 22 : CALL alloc_mat_P_omega(mat_P_omega_kp, 2, size_P, mat_P_global%matrix)
288 : END IF
289 :
290 144 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
291 22 : CALL cp_cfm_create(cfm_mat_RI_global_work, cfm_matrix_Minv_L_kpoints(1)%matrix_struct, set_zero=.TRUE.)
292 : END IF
293 :
294 1008 : ALLOCATE (has_mat_P_blocks(num_cells_dm/2 + 1, cut_memory, cut_memory, num_3c_repl, num_3c_repl))
295 3352 : has_mat_P_blocks = .TRUE.
296 :
297 144 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
298 : CALL reorder_cfm_mat_L(cfm_mat_Minv_L_kpoints, cfm_matrix_Minv_L_kpoints, fm_struct_sub_kp, para_env, &
299 22 : first_ikp_local, ikp_local)
300 :
301 : CALL reorder_cfm_mat_L(cfm_mat_L_kpoints, cfm_matrix_L_kpoints, fm_struct_sub_kp, para_env, &
302 22 : first_ikp_local, ikp_local)
303 :
304 22 : CALL cp_fm_struct_release(fm_struct_sub_kp)
305 :
306 : ELSE
307 : CALL reorder_mat_L(fm_mat_Minv_L_gamma, fm_matrix_Minv_L_gamma, fm_mat_Q%matrix_struct, mat_L, &
308 122 : mat_P_global%matrix, dimen_RI, dimen_RI_red)
309 : END IF
310 :
311 : ! Create Scalapack working matrix for the contraction with the metric
312 144 : IF (dimen_RI == dimen_RI_red) THEN
313 140 : CALL cp_fm_create(fm_mat_work, fm_mat_Q%matrix_struct, set_zero=.TRUE.)
314 :
315 : ELSE
316 : CALL cp_fm_create(fm_mat_work, fm_mat_Q%matrix_struct, nrow=dimen_RI, ncol=dimen_RI_red, &
317 4 : set_zero=.TRUE.)
318 :
319 : END IF
320 :
321 : ! Then its DBCSR counter part
322 144 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
323 122 : CALL dbcsr_get_info(mat_L%matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size)
324 :
325 : ! Create mat_work having the shape of the transposed of mat_L (compare with contract_P_omega_with_mat_L)
326 : NULLIFY (mat_work)
327 122 : ALLOCATE (mat_work)
328 122 : CALL dbcsr_create(mat_work, template=mat_L%matrix, row_blk_size=col_blk_size, col_blk_size=row_blk_size)
329 : END IF
330 :
331 144 : IF (do_ri_Sigma_x .OR. do_ic_model) THEN
332 :
333 : NULLIFY (mat_MinvVMinv%matrix)
334 114 : ALLOCATE (mat_MinvVMinv%matrix)
335 114 : CALL dbcsr_create(mat_MinvVMinv%matrix, template=mat_P_global%matrix)
336 114 : CALL dbcsr_set(mat_MinvVMinv%matrix, 0.0_dp)
337 :
338 : ! for kpoints we compute SinvVSinv later with kpoints
339 114 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
340 :
341 : ! get the Coulomb matrix for Sigma_x = G*V
342 : CALL dbcsr_multiply("T", "N", 1.0_dp, mat_L%matrix, mat_L%matrix, &
343 96 : 0.0_dp, mat_MinvVMinv%matrix, filter_eps=eps_filter_im_time)
344 :
345 : END IF
346 :
347 : END IF
348 :
349 144 : IF (do_ri_Sigma_x) THEN
350 :
351 : NULLIFY (mat_dm%matrix)
352 114 : ALLOCATE (mat_dm%matrix)
353 114 : CALL dbcsr_create(mat_dm%matrix, template=matrix_s(1)%matrix)
354 :
355 : END IF
356 :
357 144 : CALL timestop(handle)
358 :
359 288 : END SUBROUTINE alloc_im_time
360 :
361 : ! **************************************************************************************************
362 : !> \brief ...
363 : !> \param cfm_mo_coeff ...
364 : !> \param mo_coeff ...
365 : ! **************************************************************************************************
366 124 : SUBROUTINE create_mo_coeff(cfm_mo_coeff, mo_coeff)
367 :
368 : TYPE(cp_cfm_type), INTENT(OUT) :: cfm_mo_coeff
369 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
370 :
371 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_mo_coeff'
372 :
373 : INTEGER :: handle
374 :
375 62 : CALL timeset(routineN, handle)
376 :
377 62 : CALL cp_cfm_create(cfm_mo_coeff, mo_coeff%matrix_struct)
378 62 : CALL cp_fm_to_cfm(msourcer=mo_coeff, mtarget=cfm_mo_coeff)
379 :
380 62 : CALL timestop(handle)
381 :
382 62 : END SUBROUTINE create_mo_coeff
383 :
384 : ! **************************************************************************************************
385 : !> \brief Reorders complex k-point matrices to the subgroup distribution.
386 : !> \param cfm_mat_L ...
387 : !> \param cfm_matrix_L_kpoints ...
388 : !> \param fm_struct_sub_kp ...
389 : !> \param para_env ...
390 : !> \param first_ikp_local ...
391 : !> \param ikp_local ...
392 : ! **************************************************************************************************
393 44 : SUBROUTINE reorder_cfm_mat_L(cfm_mat_L, cfm_matrix_L_kpoints, fm_struct_sub_kp, para_env, first_ikp_local, ikp_local)
394 :
395 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_mat_L, cfm_matrix_L_kpoints
396 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_sub_kp
397 : TYPE(mp_para_env_type), POINTER :: para_env
398 : INTEGER, INTENT(IN) :: first_ikp_local, ikp_local
399 :
400 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reorder_cfm_mat_L'
401 :
402 : INTEGER :: handle, ikp
403 :
404 44 : CALL timeset(routineN, handle)
405 :
406 988 : ALLOCATE (cfm_mat_L(SIZE(cfm_matrix_L_kpoints)))
407 900 : DO ikp = 1, SIZE(cfm_matrix_L_kpoints)
408 856 : IF (ikp == first_ikp_local .OR. ikp_local == -1) THEN
409 856 : CALL cp_cfm_create(cfm_mat_L(ikp), fm_struct_sub_kp, set_zero=.TRUE.)
410 : END IF
411 856 : CALL cp_cfm_copy_general(cfm_matrix_L_kpoints(ikp), cfm_mat_L(ikp), para_env)
412 900 : IF (ASSOCIATED(cfm_matrix_L_kpoints(ikp)%matrix_struct)) CALL cp_cfm_release(cfm_matrix_L_kpoints(ikp))
413 : END DO
414 44 : DEALLOCATE (cfm_matrix_L_kpoints)
415 :
416 44 : CALL timestop(handle)
417 :
418 44 : END SUBROUTINE reorder_cfm_mat_L
419 :
420 : ! **************************************************************************************************
421 : !> \brief ...
422 : !> \param mat_P_omega ...
423 : !> \param num_integ_points ...
424 : !> \param size_P ...
425 : !> \param template ...
426 : ! **************************************************************************************************
427 22 : SUBROUTINE alloc_mat_P_omega(mat_P_omega, num_integ_points, size_P, template)
428 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega
429 : INTEGER, INTENT(IN) :: num_integ_points, size_P
430 : TYPE(dbcsr_type), POINTER :: template
431 :
432 : CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_mat_P_omega'
433 :
434 : INTEGER :: handle, i_kp, jquad
435 :
436 22 : CALL timeset(routineN, handle)
437 :
438 22 : NULLIFY (mat_P_omega)
439 22 : CALL dbcsr_allocate_matrix_set(mat_P_omega, num_integ_points, size_P)
440 450 : DO i_kp = 1, size_P
441 1306 : DO jquad = 1, num_integ_points
442 856 : ALLOCATE (mat_P_omega(jquad, i_kp)%matrix)
443 : CALL dbcsr_create(matrix=mat_P_omega(jquad, i_kp)%matrix, &
444 856 : template=template)
445 1284 : CALL dbcsr_set(mat_P_omega(jquad, i_kp)%matrix, 0.0_dp)
446 : END DO
447 : END DO
448 :
449 22 : CALL timestop(handle)
450 :
451 22 : END SUBROUTINE alloc_mat_P_omega
452 :
453 : ! **************************************************************************************************
454 : !> \brief ...
455 : !> \param fm_mat_L ...
456 : !> \param fm_matrix_Minv_L_gamma ...
457 : !> \param fm_struct_template ...
458 : !> \param mat_L ...
459 : !> \param mat_template ...
460 : !> \param dimen_RI ...
461 : !> \param dimen_RI_red ...
462 : ! **************************************************************************************************
463 366 : SUBROUTINE reorder_mat_L(fm_mat_L, fm_matrix_Minv_L_gamma, fm_struct_template, mat_L, mat_template, &
464 : dimen_RI, dimen_RI_red)
465 : TYPE(cp_fm_type), INTENT(OUT) :: fm_mat_L
466 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv_L_gamma
467 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_template
468 : TYPE(dbcsr_p_type), INTENT(OUT) :: mat_L
469 : TYPE(dbcsr_type), INTENT(IN) :: mat_template
470 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
471 :
472 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reorder_mat_L'
473 :
474 : INTEGER :: handle, nblk
475 122 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
476 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
477 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
478 : TYPE(cp_fm_type) :: fm_mat_L_transposed
479 :
480 122 : CALL timeset(routineN, handle)
481 :
482 : ! Get the fm_struct for fm_mat_L
483 122 : NULLIFY (fm_struct)
484 122 : IF (dimen_RI == dimen_RI_red) THEN
485 118 : fm_struct => fm_struct_template
486 : ELSE
487 : ! The template is assumed to be square such that we need a new fm_struct if dimensions are not equal
488 4 : CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI_red, ncol_global=dimen_RI, template_fmstruct=fm_struct_template)
489 : END IF
490 :
491 122 : CALL cp_fm_create(fm_mat_L, fm_struct)
492 122 : CALL cp_fm_set_all(fm_mat_L, 0.0_dp)
493 :
494 : ! For the transposed matric we need a different fm_struct
495 122 : IF (dimen_RI == dimen_RI_red) THEN
496 118 : fm_struct => fm_mat_L%matrix_struct
497 : ELSE
498 4 : CALL cp_fm_struct_release(fm_struct)
499 :
500 : ! Create a fm_struct with transposed sizes
501 4 : NULLIFY (fm_struct)
502 : CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI, ncol_global=dimen_RI_red, &
503 4 : template_fmstruct=fm_mat_L%matrix_struct)
504 : END IF
505 :
506 122 : CALL cp_fm_create(fm_mat_L_transposed, fm_struct)
507 122 : CALL cp_fm_set_all(matrix=fm_mat_L_transposed, alpha=0.0_dp)
508 :
509 122 : IF (dimen_RI /= dimen_RI_red) CALL cp_fm_struct_release(fm_struct)
510 :
511 122 : CALL cp_fm_get_info(fm_mat_L_transposed, context=blacs_env)
512 :
513 122 : CALL cp_fm_copy_general(fm_matrix_Minv_L_gamma(1, 1), fm_mat_L_transposed, blacs_env%para_env)
514 122 : CALL cp_fm_transpose(fm_mat_L_transposed, fm_mat_L)
515 :
516 : ! Release old matrix
517 122 : CALL cp_fm_release(fm_matrix_Minv_L_gamma)
518 : ! Release buffer
519 122 : CALL cp_fm_release(fm_mat_L_transposed)
520 :
521 : NULLIFY (mat_L%matrix)
522 122 : ALLOCATE (mat_L%matrix)
523 122 : IF (dimen_RI == dimen_RI_red) THEN
524 118 : CALL dbcsr_create(mat_L%matrix, template=mat_template)
525 : ELSE
526 4 : CALL dbcsr_get_info(mat_template, nblkrows_total=nblk, col_blk_size=col_blk_size)
527 :
528 4 : CALL calculate_equal_blk_size(row_blk_size, dimen_RI_red, nblk)
529 :
530 4 : CALL dbcsr_create(mat_L%matrix, template=mat_template, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
531 :
532 4 : DEALLOCATE (row_blk_size)
533 : END IF
534 :
535 122 : CALL copy_fm_to_dbcsr(fm_mat_L, mat_L%matrix)
536 :
537 122 : CALL timestop(handle)
538 :
539 122 : END SUBROUTINE reorder_mat_L
540 :
541 : ! **************************************************************************************************
542 : !> \brief ...
543 : !> \param blk_size_new ...
544 : !> \param dimen_RI_red ...
545 : !> \param nblk ...
546 : ! **************************************************************************************************
547 4 : SUBROUTINE calculate_equal_blk_size(blk_size_new, dimen_RI_red, nblk)
548 : INTEGER, DIMENSION(:), POINTER :: blk_size_new
549 : INTEGER, INTENT(IN) :: dimen_RI_red, nblk
550 :
551 : INTEGER :: col_per_blk, remainder
552 :
553 4 : NULLIFY (blk_size_new)
554 12 : ALLOCATE (blk_size_new(nblk))
555 :
556 4 : remainder = MOD(dimen_RI_red, nblk)
557 4 : col_per_blk = dimen_RI_red/nblk
558 :
559 : ! Determine a new distribution for the columns (corresponding to the number of columns)
560 10 : IF (remainder > 0) blk_size_new(1:remainder) = col_per_blk + 1
561 10 : blk_size_new(remainder + 1:nblk) = col_per_blk
562 :
563 4 : END SUBROUTINE calculate_equal_blk_size
564 :
565 : ! **************************************************************************************************
566 : !> \brief ...
567 : !> \param fm_mat_S ...
568 : !> \param do_ri_sos_laplace_mp2 ...
569 : !> \param first_cycle ...
570 : !> \param virtual ...
571 : !> \param Eigenval ...
572 : !> \param homo ...
573 : !> \param omega ...
574 : !> \param omega_old ...
575 : !> \param jquad ...
576 : !> \param mm_style ...
577 : !> \param dimen_RI ...
578 : !> \param dimen_ia ...
579 : !> \param alpha ...
580 : !> \param fm_mat_Q ...
581 : !> \param fm_mat_Q_gemm ...
582 : !> \param do_bse ...
583 : !> \param fm_mat_Q_static_bse_gemm ...
584 : !> \param dgemm_counter ...
585 : !> \param num_integ_points ...
586 : !> \param count_ev_sc_GW ...
587 : ! **************************************************************************************************
588 25038 : SUBROUTINE calc_mat_Q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, virtual, &
589 12519 : Eigenval, homo, omega, omega_old, jquad, mm_style, dimen_RI, dimen_ia, alpha, fm_mat_Q, fm_mat_Q_gemm, &
590 : do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
591 : num_integ_points, count_ev_sc_GW)
592 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
593 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, first_cycle
594 : INTEGER, INTENT(IN) :: virtual
595 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
596 : INTEGER, INTENT(IN) :: homo
597 : REAL(KIND=dp), INTENT(IN) :: omega, omega_old
598 : INTEGER, INTENT(IN) :: jquad, mm_style, dimen_RI, dimen_ia
599 : REAL(KIND=dp), INTENT(IN) :: alpha
600 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q, fm_mat_Q_gemm
601 : LOGICAL, INTENT(IN) :: do_bse
602 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q_static_bse_gemm
603 : TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter
604 : INTEGER, INTENT(IN) :: num_integ_points, count_ev_sc_GW
605 :
606 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mat_Q'
607 :
608 : INTEGER :: handle
609 :
610 12519 : CALL timeset(routineN, handle)
611 :
612 12519 : IF (do_ri_sos_laplace_mp2) THEN
613 : ! the first index of tau_tj starts with 0 (see mp2_weights)
614 88 : CALL calc_fm_mat_S_laplace(fm_mat_S, homo, virtual, Eigenval, omega - omega_old)
615 : ELSE
616 : CALL calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval, &
617 12431 : homo, omega, omega_old)
618 : END IF
619 :
620 : CALL contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, &
621 12519 : fm_mat_Q, dgemm_counter)
622 : ! fm_mat_Q_static_bse_gemm does not enter W_ijab (A matrix in TDA), but only full ABBA
623 : ! (since only B_ij_bar enters W_ijab)
624 : ! Changing jquad, since omega=0 is at last idx
625 : ! We enforce W0 for BSE in case of evGW
626 12519 : IF (do_bse .AND. jquad == num_integ_points .AND. count_ev_sc_GW == 1) THEN
627 64 : CALL cp_fm_to_fm(fm_mat_Q_gemm, fm_mat_Q_static_bse_gemm)
628 : END IF
629 12519 : CALL timestop(handle)
630 :
631 12519 : END SUBROUTINE calc_mat_Q
632 :
633 : ! **************************************************************************************************
634 : !> \brief ...
635 : !> \param fm_mat_S ...
636 : !> \param virtual ...
637 : !> \param Eigenval_last ...
638 : !> \param homo ...
639 : !> \param omega_old ...
640 : ! **************************************************************************************************
641 282 : SUBROUTINE remove_scaling_factor_rpa(fm_mat_S, virtual, Eigenval_last, homo, omega_old)
642 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
643 : INTEGER, INTENT(IN) :: virtual
644 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval_last
645 : INTEGER, INTENT(IN) :: homo
646 : REAL(KIND=dp), INTENT(IN) :: omega_old
647 :
648 : CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_scaling_factor_rpa'
649 :
650 : INTEGER :: avirt, handle, i_global, iiB, iocc, &
651 : ncol_local
652 282 : INTEGER, DIMENSION(:), POINTER :: col_indices
653 : REAL(KIND=dp) :: eigen_diff
654 :
655 282 : CALL timeset(routineN, handle)
656 :
657 : ! get info of fm_mat_S
658 : CALL cp_fm_get_info(matrix=fm_mat_S, &
659 : ncol_local=ncol_local, &
660 282 : col_indices=col_indices)
661 :
662 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
663 282 : !$OMP SHARED(ncol_local,col_indices,Eigenval_last,fm_mat_S,virtual,homo,omega_old)
664 : DO iiB = 1, ncol_local
665 : i_global = col_indices(iiB)
666 :
667 : iocc = MAX(1, i_global - 1)/virtual + 1
668 : avirt = i_global - (iocc - 1)*virtual
669 : eigen_diff = Eigenval_last(avirt + homo) - Eigenval_last(iocc)
670 :
671 : fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)/ &
672 : SQRT(eigen_diff/(eigen_diff**2 + omega_old**2))
673 :
674 : END DO
675 :
676 282 : CALL timestop(handle)
677 :
678 282 : END SUBROUTINE remove_scaling_factor_rpa
679 :
680 : ! **************************************************************************************************
681 : !> \brief ...
682 : !> \param fm_mat_S ...
683 : !> \param first_cycle ...
684 : !> \param virtual ...
685 : !> \param Eigenval ...
686 : !> \param homo ...
687 : !> \param omega ...
688 : !> \param omega_old ...
689 : ! **************************************************************************************************
690 12521 : SUBROUTINE calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval, homo, &
691 : omega, omega_old)
692 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
693 : LOGICAL, INTENT(IN) :: first_cycle
694 : INTEGER, INTENT(IN) :: virtual
695 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
696 : INTEGER, INTENT(IN) :: homo
697 : REAL(KIND=dp), INTENT(IN) :: omega, omega_old
698 :
699 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_fm_mat_S_rpa'
700 :
701 : INTEGER :: avirt, handle, i_global, iiB, iocc, &
702 : ncol_local
703 12521 : INTEGER, DIMENSION(:), POINTER :: col_indices
704 : REAL(KIND=dp) :: eigen_diff
705 :
706 12521 : CALL timeset(routineN, handle)
707 :
708 : ! get info of fm_mat_S
709 : CALL cp_fm_get_info(matrix=fm_mat_S, &
710 : ncol_local=ncol_local, &
711 12521 : col_indices=col_indices)
712 :
713 : ! update G matrix with the new value of omega
714 12521 : IF (first_cycle) THEN
715 : ! In this case just update the matrix (symmetric form) with
716 : ! SQRT((epsi_a-epsi_i)/((epsi_a-epsi_i)**2+omega**2))
717 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
718 418 : !$OMP SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega)
719 : DO iiB = 1, ncol_local
720 : i_global = col_indices(iiB)
721 :
722 : iocc = MAX(1, i_global - 1)/virtual + 1
723 : avirt = i_global - (iocc - 1)*virtual
724 : eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
725 :
726 : fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)* &
727 : SQRT(eigen_diff/(eigen_diff**2 + omega**2))
728 :
729 : END DO
730 : ELSE
731 : ! In this case the update has to remove the old omega component thus
732 : ! SQRT(((epsi_a-epsi_i)**2+omega_old**2)/((epsi_a-epsi_i)**2+omega**2))
733 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
734 12103 : !$OMP SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega,omega_old)
735 : DO iiB = 1, ncol_local
736 : i_global = col_indices(iiB)
737 :
738 : iocc = MAX(1, i_global - 1)/virtual + 1
739 : avirt = i_global - (iocc - 1)*virtual
740 : eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
741 :
742 : fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)* &
743 : SQRT((eigen_diff**2 + omega_old**2)/(eigen_diff**2 + omega**2))
744 :
745 : END DO
746 : END IF
747 :
748 12521 : CALL timestop(handle)
749 :
750 12521 : END SUBROUTINE calc_fm_mat_S_rpa
751 :
752 : ! **************************************************************************************************
753 : !> \brief ...
754 : !> \param mm_style ...
755 : !> \param dimen_RI ...
756 : !> \param dimen_ia ...
757 : !> \param alpha ...
758 : !> \param fm_mat_S ...
759 : !> \param fm_mat_Q_gemm ...
760 : !> \param fm_mat_Q ...
761 : !> \param dgemm_counter ...
762 : ! **************************************************************************************************
763 12519 : SUBROUTINE contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, &
764 : fm_mat_Q, dgemm_counter)
765 :
766 : INTEGER, INTENT(IN) :: mm_style, dimen_RI, dimen_ia
767 : REAL(KIND=dp), INTENT(IN) :: alpha
768 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S, fm_mat_Q_gemm, fm_mat_Q
769 : TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter
770 :
771 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_S_to_Q'
772 :
773 : INTEGER :: handle
774 :
775 12519 : CALL timeset(routineN, handle)
776 :
777 12519 : CALL dgemm_counter_start(dgemm_counter)
778 25028 : SELECT CASE (mm_style)
779 : CASE (wfc_mm_style_gemm)
780 : ! waste-fully computes the full symmetrix matrix, but maybe faster than cp_fm_syrk for optimized cp_fm_gemm !!!
781 : CALL parallel_gemm(transa="N", transb="T", m=dimen_RI, n=dimen_RI, k=dimen_ia, alpha=alpha, &
782 : matrix_a=fm_mat_S, matrix_b=fm_mat_S, beta=0.0_dp, &
783 12509 : matrix_c=fm_mat_Q_gemm)
784 : CASE (wfc_mm_style_syrk)
785 : ! will only compute the upper half of the matrix, which is fine, since we only use it for cholesky later
786 : CALL cp_fm_syrk(uplo='U', trans='N', k=dimen_ia, alpha=alpha, matrix_a=fm_mat_S, &
787 10 : ia=1, ja=1, beta=0.0_dp, matrix_c=fm_mat_Q_gemm)
788 : CASE DEFAULT
789 12519 : CPABORT("Unknown mm_style for contract_S_to_Q")
790 : END SELECT
791 12519 : CALL dgemm_counter_stop(dgemm_counter, dimen_RI, dimen_RI, dimen_ia)
792 :
793 : ! copy/redistribute fm_mat_Q_gemm to fm_mat_Q
794 12519 : CALL cp_fm_set_all(matrix=fm_mat_Q, alpha=0.0_dp)
795 : CALL cp_fm_to_fm_submat_general(fm_mat_Q_gemm, fm_mat_Q, dimen_RI, dimen_RI, 1, 1, 1, 1, &
796 12519 : fm_mat_Q_gemm%matrix_struct%context)
797 :
798 12519 : CALL timestop(handle)
799 :
800 12519 : END SUBROUTINE contract_S_to_Q
801 :
802 : ! **************************************************************************************************
803 : !> \brief Factorize epsilon in place and accumulate the frequency-weighted RPA energy.
804 : !> \param dimen_RI dimension of the active RI block
805 : !> \param diag_Qomega global diagonal elements of Q
806 : !> \param fm_mat_Q epsilon on entry, its upper Cholesky factor on exit
807 : !> \param para_env_RPA ...
808 : !> \param Erpa ...
809 : !> \param wjquad ...
810 : ! **************************************************************************************************
811 12721 : SUBROUTINE compute_Erpa_by_freq_int(dimen_RI, diag_Qomega, fm_mat_Q, para_env_RPA, Erpa, wjquad)
812 :
813 : INTEGER, INTENT(IN) :: dimen_RI
814 : REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(IN) :: diag_Qomega
815 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
816 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA
817 : REAL(KIND=dp), INTENT(INOUT) :: Erpa
818 : REAL(KIND=dp), INTENT(IN) :: wjquad
819 :
820 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Erpa_by_freq_int'
821 :
822 : INTEGER :: handle, iiB, info_chol
823 : REAL(KIND=dp) :: FComega
824 12721 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Q_log
825 :
826 12721 : CALL timeset(routineN, handle)
827 :
828 : ! epsilon = U^T U; log(det(epsilon)) = 2*sum(log(U_ii)).
829 12721 : CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q, n=dimen_RI, info_out=info_chol)
830 12721 : IF (info_chol /= 0) THEN
831 : CALL cp_warn(__LOCATION__, &
832 : "The Cholesky decomposition before inverting the RPA matrix / dielectric "// &
833 : "function failed. "// &
834 : "In case of low-scaling RPA/GW, decreasing EPS_FILTER in the &LOW_SCALING "// &
835 : "section might "// &
836 : "increase the overall accuracy making the matrix positive definite. "// &
837 0 : "Code will abort.")
838 : END IF
839 :
840 12721 : CPASSERT(info_chol == 0)
841 :
842 38163 : ALLOCATE (Q_log(dimen_RI))
843 12721 : CALL cp_fm_get_diag(fm_mat_Q, Q_log, n_active=dimen_RI)
844 1072186 : Q_log = 2.0_dp*LOG(Q_log)
845 :
846 : ! the following frequency integration is Eq. (27) in M. Del Ben et al., JCTC 9, 2654 (2013)
847 : ! (https://doi.org/10.1021/ct4002202)
848 : FComega = 0.0_dp
849 1072186 : DO iiB = 1, dimen_RI
850 1059465 : IF (MODULO(iiB, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE
851 1072186 : FComega = FComega + (Q_log(iiB) - diag_Qomega(iiB))/2.0_dp
852 : END DO
853 12721 : Erpa = Erpa + FComega*wjquad
854 :
855 12721 : DEALLOCATE (Q_log)
856 :
857 12721 : CALL timestop(handle)
858 :
859 12721 : END SUBROUTINE compute_Erpa_by_freq_int
860 :
861 : ! **************************************************************************************************
862 : !> \brief ...
863 : !> \param fm_struct_sub_kp ...
864 : !> \param para_env ...
865 : !> \param dimen_RI ...
866 : !> \param ikp_local ...
867 : !> \param first_ikp_local ...
868 : ! **************************************************************************************************
869 22 : SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_RI, &
870 : ikp_local, first_ikp_local)
871 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_sub_kp
872 : TYPE(mp_para_env_type), POINTER :: para_env
873 : INTEGER, INTENT(IN) :: dimen_RI
874 : INTEGER, INTENT(OUT) :: ikp_local, first_ikp_local
875 :
876 : CHARACTER(len=*), PARAMETER :: routineN = 'get_sub_para_kp'
877 :
878 : INTEGER :: color_sub_kp, handle, num_proc_per_kp
879 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub_kp
880 : TYPE(mp_para_env_type), POINTER :: para_env_sub_kp
881 :
882 22 : CALL timeset(routineN, handle)
883 :
884 : ! we use all processors for every k-point, subgroups for cp_cfm_heevd only seems to work for
885 : ! very small subgroups with 1, 2, or 3 MPI ranks. For more MPI-ranks, eigenvalues and
886 : ! eigenvectors coming out of cp_cfm_heevd are totally wrong unfortunately.
887 22 : num_proc_per_kp = para_env%num_pe
888 :
889 : ! IF(nkp > para_env%num_pe) THEN
890 : ! num_proc_per_kp = para_env%num_pe
891 : ! ELSE
892 : ! num_proc_per_kp = para_env%num_pe/nkp
893 : ! END IF
894 :
895 22 : color_sub_kp = para_env%mepos/num_proc_per_kp
896 22 : ALLOCATE (para_env_sub_kp)
897 22 : CALL para_env_sub_kp%from_split(para_env, color_sub_kp)
898 :
899 : ! grid_2d(1) = 1
900 : ! grid_2d(2) = para_env_sub_kp%num_pe
901 :
902 22 : NULLIFY (blacs_env_sub_kp)
903 : ! CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp, grid_2d=grid_2d)
904 22 : CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp)
905 :
906 22 : NULLIFY (fm_struct_sub_kp)
907 : CALL cp_fm_struct_create(fm_struct_sub_kp, context=blacs_env_sub_kp, nrow_global=dimen_RI, &
908 22 : ncol_global=dimen_RI, para_env=para_env_sub_kp)
909 :
910 22 : CALL cp_blacs_env_release(blacs_env_sub_kp)
911 :
912 : ! IF(nkp > para_env%num_pe) THEN
913 : ! every processor has all ikp's
914 22 : ikp_local = -1
915 22 : first_ikp_local = 1
916 : ! ELSE
917 : ! ikp_local = 0
918 : ! first_ikp_local = 1
919 : ! DO ikp = 1, nkp
920 : ! IF(MOD(ikp-1, para_env%num_pe/num_proc_per_kp) == color_sub_kp) THEN
921 : ! ikp_local = ikp
922 : ! first_ikp_local = ikp
923 : ! END IF
924 : ! END DO
925 : ! END IF
926 :
927 22 : CALL mp_para_env_release(para_env_sub_kp)
928 :
929 22 : CALL timestop(handle)
930 :
931 22 : END SUBROUTINE get_sub_para_kp
932 :
933 : ! **************************************************************************************************
934 : !> \brief ...
935 : !> \param cfm_mo_coeff ...
936 : !> \param index_to_cell_3c ...
937 : !> \param cell_to_index_3c ...
938 : !> \param do_ic_model ...
939 : !> \param do_kpoints_cubic_RPA ...
940 : !> \param do_kpoints_from_Gamma ...
941 : !> \param do_ri_Sigma_x ...
942 : !> \param has_mat_P_blocks ...
943 : !> \param wkp_W ...
944 : !> \param cfm_mat_Q ...
945 : !> \param cfm_mat_Minv_L_kpoints ...
946 : !> \param cfm_mat_L_kpoints ...
947 : !> \param fm_mat_Minv_L_gamma ...
948 : !> \param fm_matrix_Minv ...
949 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
950 : !> \param cfm_mat_RI_global_work ...
951 : !> \param fm_mat_work ...
952 : !> \param mat_dm ...
953 : !> \param mat_L ...
954 : !> \param mat_MinvVMinv ...
955 : !> \param mat_P_omega ...
956 : !> \param mat_P_omega_kp ...
957 : !> \param t_3c_M ...
958 : !> \param t_3c_O ...
959 : !> \param t_3c_O_compressed ...
960 : !> \param t_3c_O_ind ...
961 : !> \param mat_work ...
962 : !> \param qs_env ...
963 : ! **************************************************************************************************
964 144 : SUBROUTINE dealloc_im_time(cfm_mo_coeff, index_to_cell_3c, &
965 : cell_to_index_3c, do_ic_model, &
966 : do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, &
967 : has_mat_P_blocks, &
968 : wkp_W, cfm_mat_Q, cfm_mat_Minv_L_kpoints, cfm_mat_L_kpoints, &
969 : fm_mat_Minv_L_gamma, &
970 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
971 : cfm_mat_RI_global_work, fm_mat_work, mat_dm, mat_L, &
972 : mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
973 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
974 : mat_work, qs_env)
975 :
976 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:), &
977 : INTENT(INOUT) :: cfm_mo_coeff
978 : INTEGER, ALLOCATABLE, DIMENSION(:, :), &
979 : INTENT(INOUT) :: index_to_cell_3c
980 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
981 : INTENT(INOUT) :: cell_to_index_3c
982 : LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA, &
983 : do_kpoints_from_Gamma, do_ri_Sigma_x
984 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), &
985 : INTENT(INOUT) :: has_mat_P_blocks
986 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
987 : INTENT(INOUT) :: wkp_W
988 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_mat_Q
989 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_mat_Minv_L_kpoints, cfm_mat_L_kpoints
990 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_Minv_L_gamma
991 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv, &
992 : fm_matrix_Minv_Vtrunc_Minv
993 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_mat_RI_global_work
994 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_work
995 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_dm, mat_L, mat_MinvVMinv
996 : TYPE(dbcsr_p_type), ALLOCATABLE, &
997 : DIMENSION(:, :, :), INTENT(INOUT) :: mat_P_omega
998 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_kp
999 : TYPE(dbt_type) :: t_3c_M
1000 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O
1001 : TYPE(hfx_compression_type), ALLOCATABLE, &
1002 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_compressed
1003 : TYPE(block_ind_type), ALLOCATABLE, &
1004 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_ind
1005 : TYPE(dbcsr_type), POINTER :: mat_work
1006 : TYPE(qs_environment_type), POINTER :: qs_env
1007 :
1008 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dealloc_im_time'
1009 :
1010 : INTEGER :: cut_memory, handle, i_kp, i_mem, i_size, &
1011 : ispin, j_size, jquad, unused
1012 :
1013 144 : CALL timeset(routineN, handle)
1014 :
1015 144 : IF (ALLOCATED(cfm_mo_coeff)) THEN
1016 112 : DO ispin = 1, SIZE(cfm_mo_coeff)
1017 112 : CALL cp_cfm_release(cfm_mo_coeff(ispin))
1018 : END DO
1019 50 : DEALLOCATE (cfm_mo_coeff)
1020 : END IF
1021 144 : IF (ALLOCATED(cfm_mat_Minv_L_kpoints)) THEN
1022 18 : DO i_kp = 1, SIZE(cfm_mat_Minv_L_kpoints)
1023 18 : IF (ASSOCIATED(cfm_mat_Minv_L_kpoints(i_kp)%matrix_struct)) THEN
1024 12 : CALL cp_cfm_release(cfm_mat_Minv_L_kpoints(i_kp))
1025 : END IF
1026 : END DO
1027 6 : DEALLOCATE (cfm_mat_Minv_L_kpoints)
1028 : END IF
1029 144 : IF (ALLOCATED(cfm_mat_L_kpoints)) THEN
1030 18 : DO i_kp = 1, SIZE(cfm_mat_L_kpoints)
1031 18 : IF (ASSOCIATED(cfm_mat_L_kpoints(i_kp)%matrix_struct)) CALL cp_cfm_release(cfm_mat_L_kpoints(i_kp))
1032 : END DO
1033 6 : DEALLOCATE (cfm_mat_L_kpoints)
1034 : END IF
1035 144 : CALL cp_fm_release(fm_mat_Minv_L_gamma)
1036 :
1037 144 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
1038 22 : CALL cp_fm_release(fm_matrix_Minv_Vtrunc_Minv)
1039 22 : CALL cp_fm_release(fm_matrix_Minv)
1040 : END IF
1041 :
1042 144 : CALL cp_fm_release(fm_mat_work)
1043 :
1044 144 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
1045 122 : CALL dbcsr_release(mat_work)
1046 122 : DEALLOCATE (mat_work)
1047 : END IF
1048 :
1049 144 : IF (ASSOCIATED(mat_L%matrix)) THEN
1050 122 : CALL dbcsr_release(mat_L%matrix)
1051 122 : DEALLOCATE (mat_L%matrix)
1052 : END IF
1053 :
1054 144 : IF (do_ri_Sigma_x .OR. do_ic_model) THEN
1055 114 : CALL dbcsr_release(mat_MinvVMinv%matrix)
1056 114 : DEALLOCATE (mat_MinvVMinv%matrix)
1057 : END IF
1058 144 : IF (do_ri_Sigma_x) THEN
1059 114 : CALL dbcsr_release(mat_dm%matrix)
1060 114 : DEALLOCATE (mat_dm%matrix)
1061 : END IF
1062 :
1063 144 : DEALLOCATE (index_to_cell_3c, cell_to_index_3c)
1064 :
1065 144 : IF (ALLOCATED(mat_P_omega)) THEN
1066 320 : DO ispin = 1, SIZE(mat_P_omega, 3)
1067 952 : DO i_kp = 1, SIZE(mat_P_omega, 2)
1068 4798 : DO jquad = 1, SIZE(mat_P_omega, 1)
1069 4622 : CALL dbcsr_deallocate_matrix(mat_P_omega(jquad, i_kp, ispin)%matrix)
1070 : END DO
1071 : END DO
1072 : END DO
1073 144 : DEALLOCATE (mat_P_omega)
1074 : END IF
1075 :
1076 312 : DO j_size = 1, SIZE(t_3c_O, 2)
1077 600 : DO i_size = 1, SIZE(t_3c_O, 1)
1078 456 : CALL dbt_destroy(t_3c_O(i_size, j_size))
1079 : END DO
1080 : END DO
1081 :
1082 432 : DEALLOCATE (t_3c_O)
1083 144 : CALL dbt_destroy(t_3c_M)
1084 :
1085 144 : DEALLOCATE (has_mat_P_blocks)
1086 :
1087 144 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
1088 22 : CALL cp_cfm_release(cfm_mat_Q)
1089 22 : CALL cp_cfm_release(cfm_mat_RI_global_work)
1090 22 : CALL dbcsr_deallocate_matrix_set(mat_P_omega_kp)
1091 22 : DEALLOCATE (wkp_W)
1092 : END IF
1093 :
1094 144 : cut_memory = SIZE(t_3c_O_compressed, 3)
1095 :
1096 620 : DEALLOCATE (t_3c_O_ind)
1097 428 : DO i_mem = 1, cut_memory
1098 744 : DO j_size = 1, SIZE(t_3c_O_compressed, 2)
1099 1076 : DO i_size = 1, SIZE(t_3c_O_compressed, 1)
1100 792 : CALL dealloc_containers(t_3c_O_compressed(i_size, j_size, i_mem), unused)
1101 : END DO
1102 : END DO
1103 : END DO
1104 144 : DEALLOCATE (t_3c_O_compressed)
1105 :
1106 144 : IF (do_kpoints_from_Gamma) THEN
1107 16 : CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_G)
1108 16 : IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
1109 16 : CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma)
1110 16 : CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc)
1111 : END IF
1112 : END IF
1113 :
1114 144 : CALL timestop(handle)
1115 :
1116 144 : END SUBROUTINE dealloc_im_time
1117 :
1118 : ! **************************************************************************************************
1119 : !> \brief ...
1120 : !> \param mat_P_omega ...
1121 : !> \param mat_L ...
1122 : !> \param mat_work ...
1123 : !> \param eps_filter_im_time ...
1124 : !> \param fm_mat_work ...
1125 : !> \param dimen_RI ...
1126 : !> \param dimen_RI_red ...
1127 : !> \param fm_mat_L ...
1128 : !> \param fm_mat_Q ...
1129 : ! **************************************************************************************************
1130 1338 : SUBROUTINE contract_P_omega_with_mat_L(mat_P_omega, mat_L, mat_work, eps_filter_im_time, fm_mat_work, dimen_RI, &
1131 : dimen_RI_red, fm_mat_L, fm_mat_Q)
1132 :
1133 : TYPE(dbcsr_type), INTENT(IN) :: mat_P_omega, mat_L
1134 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_work
1135 : REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time
1136 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_work
1137 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
1138 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_L, fm_mat_Q
1139 :
1140 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_P_omega_with_mat_L'
1141 :
1142 : INTEGER :: handle
1143 :
1144 1338 : CALL timeset(routineN, handle)
1145 :
1146 : ! multiplication with RI metric/Coulomb operator
1147 : CALL dbcsr_multiply("N", "T", 1.0_dp, mat_P_omega, mat_L, &
1148 1338 : 0.0_dp, mat_work, filter_eps=eps_filter_im_time)
1149 :
1150 1338 : CALL copy_dbcsr_to_fm(mat_work, fm_mat_work)
1151 :
1152 : CALL parallel_gemm('N', 'N', dimen_RI_red, dimen_RI_red, dimen_RI, 1.0_dp, fm_mat_L, fm_mat_work, &
1153 1338 : 0.0_dp, fm_mat_Q)
1154 :
1155 : ! Reset mat_work to save memory
1156 1338 : CALL dbcsr_set(mat_work, 0.0_dp)
1157 1338 : CALL dbcsr_filter(mat_work, 1.0_dp)
1158 :
1159 1338 : CALL timestop(handle)
1160 :
1161 1338 : END SUBROUTINE contract_P_omega_with_mat_L
1162 :
1163 : END MODULE rpa_util
|
