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