(git:744416f)
Loading...
Searching...
No Matches
rpa_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Utility functions for RPA calculations
10!> \par History
11!> 06.2019 Moved from rpa_ri_gpw.F [Frederick Stein]
12! **************************************************************************************************
14
15 USE cell_types, ONLY: cell_type,&
20 USE cp_cfm_types, ONLY: cp_cfm_create,&
25 USE cp_dbcsr_api, ONLY: &
46 USE dbt_api, ONLY: dbt_destroy,&
47 dbt_type
51 USE hfx_types, ONLY: block_ind_type,&
56 USE kinds, ONLY: dp
57 USE kpoint_types, ONLY: get_kpoint_info,&
60 USE mathconstants, ONLY: z_zero
68#include "./base/base_uses.f90"
69
70 IMPLICIT NONE
71
72 PRIVATE
73
74 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_util'
75
78
79CONTAINS
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 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 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 INTEGER, DIMENSION(:), POINTER :: row_blk_size
191 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 CALL timeset(routinen, handle)
196
197 ALLOCATE (cfm_mo_coeff(nspins))
198
199 DO ispin = 1, SIZE(mo_coeff)
200 CALL create_mo_coeff(cfm_mo_coeff(ispin), mo_coeff(ispin))
201 END DO
202
203 num_3c_repl = SIZE(t_3c_o, 2)
204
205 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 DO i_dim = 1, 3
209 cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
210 END DO
211 num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
212 ALLOCATE (index_to_cell_3c(3, SIZE(kpoints%index_to_cell, 2)))
213 cpassert(SIZE(kpoints%index_to_cell, 1) == 3)
214 index_to_cell_3c(:, :) = kpoints%index_to_cell(:, :)
215 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 lbound(kpoints%cell_to_index, 3):ubound(kpoints%cell_to_index, 3)))
218 cell_to_index_3c(:, :, :) = kpoints%cell_to_index(:, :, :)
219
220 ELSE
221 ALLOCATE (index_to_cell_3c(3, 1))
222 index_to_cell_3c(:, 1) = 0
223 ALLOCATE (cell_to_index_3c(0:0, 0:0, 0:0))
224 cell_to_index_3c(0, 0, 0) = 1
225 num_cells_dm = 1
226 END IF
227
228 IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
229
230 CALL get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_ri, ikp_local, first_ikp_local)
231
232 CALL cp_cfm_create(cfm_mat_q, fm_struct_sub_kp)
233 CALL cp_cfm_set_all(cfm_mat_q, z_zero)
234 ELSE
235 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 IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
241
242 NULLIFY (cell)
243 CALL get_qs_env(qs_env, cell=cell)
244 CALL get_cell(cell=cell, periodic=periodic)
245
246 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 CALL compute_wkp_w(qs_env, wkp_w, wkp_v, kpoints, cell%h_inv, periodic)
250 DEALLOCATE (wkp_v)
251
252 ELSE
253 nkp = 1
254 END IF
255
256 IF (do_kpoints_cubic_rpa) THEN
257 size_p = max(num_cells_dm/2 + 1, nkp)
258 ELSE IF (do_kpoints_from_gamma) THEN
259 size_p = max(3**(periodic(1) + periodic(2) + periodic(3)), nkp)
260 ELSE
261 size_p = 1
262 END IF
263
264 nspins_p_omega = 1
265 IF (my_open_shell) nspins_p_omega = 2
266
267 ALLOCATE (mat_p_omega(num_integ_points, size_p, nspins_p_omega))
268 DO ispin = 1, nspins_p_omega
269 DO i_kp = 1, size_p
270 DO jquad = 1, num_integ_points
271 NULLIFY (mat_p_omega(jquad, i_kp, ispin)%matrix)
272 ALLOCATE (mat_p_omega(jquad, i_kp, ispin)%matrix)
273 CALL dbcsr_create(matrix=mat_p_omega(jquad, i_kp, ispin)%matrix, &
274 template=mat_p_global%matrix)
275 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 IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
281 CALL alloc_mat_p_omega(mat_p_omega_kp, 2, size_p, mat_p_global%matrix)
282 END IF
283
284 IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
285 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 ALLOCATE (has_mat_p_blocks(num_cells_dm/2 + 1, cut_memory, cut_memory, num_3c_repl, num_3c_repl))
289 has_mat_p_blocks = .true.
290
291 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 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 mat_p_global%matrix, dimen_ri, dimen_ri_red, first_ikp_local, ikp_local, fm_struct_sub_kp)
298
299 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 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 IF (dimen_ri == dimen_ri_red) THEN
308 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 set_zero=.true.)
313
314 END IF
315
316 ! Then its DBCSR counter part
317 IF (.NOT. (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma)) THEN
318 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 ALLOCATE (mat_work)
323 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 IF (do_ri_sigma_x .OR. do_ic_model) THEN
327
328 NULLIFY (mat_minvvminv%matrix)
329 ALLOCATE (mat_minvvminv%matrix)
330 CALL dbcsr_create(mat_minvvminv%matrix, template=mat_p_global%matrix)
331 CALL dbcsr_set(mat_minvvminv%matrix, 0.0_dp)
332
333 ! for kpoints we compute SinvVSinv later with kpoints
334 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 0.0_dp, mat_minvvminv%matrix, filter_eps=eps_filter_im_time)
339
340 END IF
341
342 END IF
343
344 IF (do_ri_sigma_x) THEN
345
346 NULLIFY (mat_dm%matrix)
347 ALLOCATE (mat_dm%matrix)
348 CALL dbcsr_create(mat_dm%matrix, template=matrix_s(1)%matrix)
349
350 END IF
351
352 CALL timestop(handle)
353
354 END SUBROUTINE alloc_im_time
355
356! **************************************************************************************************
357!> \brief ...
358!> \param cfm_mo_coeff ...
359!> \param mo_coeff ...
360! **************************************************************************************************
361 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 CALL timeset(routinen, handle)
371
372 CALL cp_cfm_create(cfm_mo_coeff, mo_coeff%matrix_struct)
373 CALL cp_fm_to_cfm(msourcer=mo_coeff, mtarget=cfm_mo_coeff)
374
375 CALL timestop(handle)
376
377 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 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 CALL timeset(routinen, handle)
396
397 NULLIFY (mat_p_omega)
398 CALL dbcsr_allocate_matrix_set(mat_p_omega, num_integ_points, size_p)
399 DO i_kp = 1, size_p
400 DO jquad = 1, num_integ_points
401 ALLOCATE (mat_p_omega(jquad, i_kp)%matrix)
402 CALL dbcsr_create(matrix=mat_p_omega(jquad, i_kp)%matrix, &
403 template=template)
404 CALL dbcsr_set(mat_p_omega(jquad, i_kp)%matrix, 0.0_dp)
405 END DO
406 END DO
407
408 CALL timestop(handle)
409
410 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 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 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 CALL timeset(routinen, handle)
449
450 do_kpoints = .false.
451 IF (PRESENT(ikp_local) .AND. PRESENT(fm_struct_sub_kp)) THEN
452 do_kpoints = .true.
453 END IF
454
455 ! Get the fm_struct for fm_mat_L
456 NULLIFY (fm_struct)
457 IF (dimen_ri == dimen_ri_red) THEN
458 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 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 ALLOCATE (fm_mat_l(SIZE(fm_matrix_minv_l_kpoints, 1), SIZE(fm_matrix_minv_l_kpoints, 2)))
466 DO j_size = 1, SIZE(fm_matrix_minv_l_kpoints, 2)
467 DO ikp = 1, SIZE(fm_matrix_minv_l_kpoints, 1)
468 IF (do_kpoints) THEN
469 IF (ikp == first_ikp_local .OR. ikp_local == -1) THEN
470 CALL cp_fm_create(fm_mat_l(ikp, j_size), fm_struct_sub_kp)
471 CALL cp_fm_set_all(fm_mat_l(ikp, j_size), 0.0_dp)
472 END IF
473 ELSE
474 CALL cp_fm_create(fm_mat_l(ikp, j_size), fm_struct)
475 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 IF (dimen_ri == dimen_ri_red) THEN
482 fm_struct => fm_mat_l(first_ikp_local, 1)%matrix_struct
483 ELSE
484 CALL cp_fm_struct_release(fm_struct)
485
486 ! Create a fm_struct with transposed sizes
487 NULLIFY (fm_struct)
488 CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_ri, ncol_global=dimen_ri_red, &
489 template_fmstruct=fm_mat_l(first_ikp_local, 1)%matrix_struct) !, force_block=.TRUE.)
490 END IF
491
492 ! Allocate buffer matrix
493 CALL cp_fm_create(fm_mat_l_transposed, fm_struct)
494 CALL cp_fm_set_all(matrix=fm_mat_l_transposed, alpha=0.0_dp)
495
496 IF (dimen_ri /= dimen_ri_red) CALL cp_fm_struct_release(fm_struct)
497
498 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 DO j_size = 1, SIZE(fm_matrix_minv_l_kpoints, 2)
504 DO ikp = 1, SIZE(fm_matrix_minv_l_kpoints, 1)
505 IF (do_kpoints) THEN
506 IF (ikp_local == ikp .OR. ikp_local == -1) THEN
507 CALL cp_fm_copy_general(fm_matrix_minv_l_kpoints(ikp, j_size), fm_mat_l_transposed, para_env)
508 CALL cp_fm_to_fm(fm_mat_l_transposed, fm_mat_l(ikp, j_size))
509 ELSE
510 CALL cp_fm_copy_general(fm_matrix_minv_l_kpoints(ikp, j_size), fmdummy, para_env)
511 END IF
512 ELSE
513 CALL cp_fm_copy_general(fm_matrix_minv_l_kpoints(ikp, j_size), fm_mat_l_transposed, blacs_env%para_env)
514 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 CALL cp_fm_release(fm_matrix_minv_l_kpoints)
521 ! Release buffer
522 CALL cp_fm_release(fm_mat_l_transposed)
523
524 my_allocate_mat_l = .true.
525 IF (PRESENT(allocate_mat_l)) my_allocate_mat_l = allocate_mat_l
526
527 IF (my_allocate_mat_l) THEN
528 ! Create sparse variant of L
529 NULLIFY (mat_l%matrix)
530 ALLOCATE (mat_l%matrix)
531 IF (dimen_ri == dimen_ri_red) THEN
532 CALL dbcsr_create(mat_l%matrix, template=mat_template)
533 ELSE
534 CALL dbcsr_get_info(mat_template, nblkrows_total=nblk, col_blk_size=col_blk_size)
535
536 CALL calculate_equal_blk_size(row_blk_size, dimen_ri_red, nblk)
537
538 CALL dbcsr_create(mat_l%matrix, template=mat_template, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
539
540 DEALLOCATE (row_blk_size)
541 END IF
542
543 IF (.NOT. (do_kpoints)) THEN
544 CALL copy_fm_to_dbcsr(fm_mat_l(1, 1), mat_l%matrix)
545 END IF
546
547 END IF
548
549 CALL timestop(handle)
550
551 END SUBROUTINE reorder_mat_l
552
553! **************************************************************************************************
554!> \brief ...
555!> \param blk_size_new ...
556!> \param dimen_RI_red ...
557!> \param nblk ...
558! **************************************************************************************************
559 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 NULLIFY (blk_size_new)
566 ALLOCATE (blk_size_new(nblk))
567
568 remainder = mod(dimen_ri_red, nblk)
569 col_per_blk = dimen_ri_red/nblk
570
571 ! Determine a new distribution for the columns (corresponding to the number of columns)
572 IF (remainder > 0) blk_size_new(1:remainder) = col_per_blk + 1
573 blk_size_new(remainder + 1:nblk) = col_per_blk
574
575 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 SUBROUTINE calc_mat_q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, virtual, &
601 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 CALL timeset(routinen, handle)
623
624 IF (do_ri_sos_laplace_mp2) THEN
625 ! the first index of tau_tj starts with 0 (see mp2_weights)
626 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 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 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 IF (do_bse .AND. jquad == num_integ_points .AND. count_ev_sc_gw == 1) THEN
639 CALL cp_fm_to_fm(fm_mat_q_gemm, fm_mat_q_static_bse_gemm)
640 END IF
641 CALL timestop(handle)
642
643 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices
665 REAL(kind=dp) :: eigen_diff
666
667 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 col_indices=col_indices)
673
674!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
675!$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 CALL timestop(handle)
689
690 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices
716 REAL(kind=dp) :: eigen_diff
717
718 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 col_indices=col_indices)
724
725 ! update G matrix with the new value of omega
726 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 !$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 !$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 CALL timestop(handle)
761
762 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 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 CALL timeset(routinen, handle)
788
789 CALL dgemm_counter_start(dgemm_counter)
790 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 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 ia=1, ja=1, beta=0.0_dp, matrix_c=fm_mat_q_gemm)
800 CASE DEFAULT
801 cpabort("Unknown mm_style for contract_S_to_Q")
802 END SELECT
803 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 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 fm_mat_q_gemm%matrix_struct%context)
809
810 CALL timestop(handle)
811
812 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
831 TYPE(mp_para_env_type), POINTER :: para_env
832
833 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 para_env=para_env)
841
842 ! calculate the trace of Q and add 1 on the diagonal
843 trace_qomega = 0.0_dp
844!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
845!$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 CALL para_env%sum(trace_qomega)
857
858 CALL timestop(handle)
859
860 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
885 REAL(kind=dp) :: fcomega
886 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: q_log
887
888 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 col_indices=col_indices)
895
896 ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition
897 CALL cp_fm_cholesky_decompose(matrix=fm_mat_q, n=dimen_ri, info_out=info_chol)
898 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 "Code will abort.")
906 END IF
907
908 cpassert(info_chol == 0)
909
910 ALLOCATE (q_log(dimen_ri))
911 q_log = 0.0_dp
912!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
913!$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 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 fcomega = 0.0_dp
928 DO iib = 1, dimen_ri
929 IF (modulo(iib, para_env_rpa%num_pe) /= para_env_rpa%mepos) cycle
930 fcomega = fcomega + (q_log(iib) - trace_qomega(iib))/2.0_dp
931 END DO
932 erpa = erpa + fcomega*wjquad
933
934 DEALLOCATE (q_log)
935
936 CALL timestop(handle)
937
938 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 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 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 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 color_sub_kp = para_env%mepos/num_proc_per_kp
975 ALLOCATE (para_env_sub_kp)
976 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 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 CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp)
984
985 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 ncol_global=dimen_ri, para_env=para_env_sub_kp)
988
989 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 ikp_local = -1
994 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 CALL mp_para_env_release(para_env_sub_kp)
1007
1008 CALL timestop(handle)
1009
1010 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 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 CALL timeset(routinen, handle)
1089
1090 nspins = SIZE(cfm_mo_coeff)
1091 my_open_shell = (nspins == 2)
1092
1093 DO ispin = 1, SIZE(cfm_mo_coeff)
1094 CALL cp_cfm_release(cfm_mo_coeff(ispin))
1095 END DO
1096 CALL cp_fm_release(fm_mat_minv_l_kpoints)
1097 CALL cp_fm_release(fm_mat_l_kpoints)
1098
1099 IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
1100 CALL cp_fm_release(fm_matrix_minv_vtrunc_minv)
1101 CALL cp_fm_release(fm_matrix_minv)
1102 END IF
1103
1104 CALL cp_fm_release(fm_mat_work)
1105
1106 IF (.NOT. (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma)) THEN
1107 CALL dbcsr_release(mat_work)
1108 DEALLOCATE (mat_work)
1109 END IF
1110
1111 CALL dbcsr_release(mat_l%matrix)
1112 DEALLOCATE (mat_l%matrix)
1113
1114 IF (do_ri_sigma_x .OR. do_ic_model) THEN
1115 CALL dbcsr_release(mat_minvvminv%matrix)
1116 DEALLOCATE (mat_minvvminv%matrix)
1117 END IF
1118 IF (do_ri_sigma_x) THEN
1119 CALL dbcsr_release(mat_dm%matrix)
1120 DEALLOCATE (mat_dm%matrix)
1121 END IF
1122
1123 DEALLOCATE (index_to_cell_3c, cell_to_index_3c)
1124
1125 IF (ALLOCATED(mat_p_omega)) THEN
1126 DO ispin = 1, SIZE(mat_p_omega, 3)
1127 DO i_kp = 1, SIZE(mat_p_omega, 2)
1128 DO jquad = 1, SIZE(mat_p_omega, 1)
1129 CALL dbcsr_deallocate_matrix(mat_p_omega(jquad, i_kp, ispin)%matrix)
1130 END DO
1131 END DO
1132 END DO
1133 DEALLOCATE (mat_p_omega)
1134 END IF
1135
1136 DO j_size = 1, SIZE(t_3c_o, 2)
1137 DO i_size = 1, SIZE(t_3c_o, 1)
1138 CALL dbt_destroy(t_3c_o(i_size, j_size))
1139 END DO
1140 END DO
1141
1142 DEALLOCATE (t_3c_o)
1143 CALL dbt_destroy(t_3c_m)
1144
1145 DEALLOCATE (has_mat_p_blocks)
1146
1147 IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
1148 CALL cp_cfm_release(cfm_mat_q)
1149 CALL cp_fm_release(fm_mat_ri_global_work)
1150 CALL dbcsr_deallocate_matrix_set(mat_p_omega_kp)
1151 DEALLOCATE (wkp_w)
1152 END IF
1153
1154 cut_memory = SIZE(t_3c_o_compressed, 3)
1155
1156 DEALLOCATE (t_3c_o_ind)
1157 DO i_mem = 1, cut_memory
1158 DO j_size = 1, SIZE(t_3c_o_compressed, 2)
1159 DO i_size = 1, SIZE(t_3c_o_compressed, 1)
1160 CALL dealloc_containers(t_3c_o_compressed(i_size, j_size, i_mem), unused)
1161 END DO
1162 END DO
1163 END DO
1164 DEALLOCATE (t_3c_o_compressed)
1165
1166 IF (do_kpoints_from_gamma) THEN
1167 CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_G)
1168 IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
1169 CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma)
1170 CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc)
1171 END IF
1172 END IF
1173
1174 CALL timestop(handle)
1175
1176 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 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 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 0.0_dp, mat_work, filter_eps=eps_filter_im_time)
1209
1210 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 0.0_dp, fm_mat_q)
1214
1215 ! Reset mat_work to save memory
1216 CALL dbcsr_set(mat_work, 0.0_dp)
1217 CALL dbcsr_filter(mat_work, 1.0_dp)
1218
1219 CALL timestop(handle)
1220
1221 END SUBROUTINE contract_p_omega_with_mat_l
1222
1223END MODULE rpa_util
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
subroutine, public dbcsr_deallocate_matrix(matrix)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
subroutine, public cp_fm_syrk(uplo, trans, k, alpha, matrix_a, ia, ja, beta, matrix_c)
performs a rank-k update of a symmetric matrix_c matrix_c = beta * matrix_c + alpha * matrix_a * tran...
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_to_fm_submat_general(source, destination, nrows, ncols, s_firstrow, s_firstcol, d_firstrow, d_firstcol, global_context)
General copy of a submatrix of fm matrix to a submatrix of another fm matrix. The two matrices can ha...
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
Counters to determine the performance of parallel DGEMMs.
subroutine, public dgemm_counter_start(dgemm_counter)
start timer of the counter
subroutine, public dgemm_counter_stop(dgemm_counter, size1, size2, size3)
stop timer of the counter and provide matrix sizes
Types and set/get functions for HFX.
Definition hfx_types.F:16
subroutine, public dealloc_containers(data, memory_usage)
...
Definition hfx_types.F:2957
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public wfc_mm_style_syrk
integer, parameter, public wfc_mm_style_gemm
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Types and basic routines needed for a kpoint calculation.
subroutine, public kpoint_release(kpoint)
Release a kpoint environment, deallocate all data.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_zero
Interface to the message passing library MPI.
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Routines to calculate MP2 energy with laplace approach.
Definition mp2_laplace.F:13
subroutine, public calc_fm_mat_s_laplace(fm_mat_s, homo, virtual, eigenval, dajquad)
...
Definition mp2_laplace.F:39
basic linear algebra operations for full matrixes
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Routines treating GW and RPA calculations with kpoints.
subroutine, public compute_wkp_w(qs_env, wkp_w, wkp_v, kpoints, h_inv, periodic)
...
Utility functions for RPA calculations.
Definition rpa_util.F:13
subroutine, public alloc_im_time(qs_env, para_env, dimen_ri, dimen_ri_red, num_integ_points, nspins, fm_mat_q, cfm_mo_coeff, fm_matrix_minv_l_kpoints, fm_matrix_l_kpoints, mat_p_global, t_3c_o, matrix_s, kpoints, eps_filter_im_time, cut_memory, nkp, num_cells_dm, num_3c_repl, size_p, ikp_local, index_to_cell_3c, cell_to_index_3c, col_blk_size, do_ic_model, do_kpoints_cubic_rpa, do_kpoints_from_gamma, do_ri_sigma_x, my_open_shell, has_mat_p_blocks, wkp_w, cfm_mat_q, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, fm_mat_ri_global_work, fm_mat_work, mat_dm, mat_l, mat_m_p_munu_occ, mat_m_p_munu_virt, mat_minvvminv, mat_p_omega, mat_p_omega_kp, mat_work, mo_coeff)
...
Definition rpa_util.F:145
subroutine, public compute_erpa_by_freq_int(dimen_ri, trace_qomega, fm_mat_q, para_env_rpa, erpa, wjquad)
...
Definition rpa_util.F:872
subroutine, public q_trace_and_add_unit_matrix(dimen_ri, trace_qomega, fm_mat_q)
...
Definition rpa_util.F:821
subroutine, public calc_mat_q(fm_mat_s, do_ri_sos_laplace_mp2, first_cycle, virtual, eigenval, homo, omega, omega_old, jquad, mm_style, dimen_ri, dimen_ia, alpha, fm_mat_q, fm_mat_q_gemm, do_bse, fm_mat_q_static_bse_gemm, dgemm_counter, num_integ_points, count_ev_sc_gw)
...
Definition rpa_util.F:604
subroutine, public contract_p_omega_with_mat_l(mat_p_omega, mat_l, mat_work, eps_filter_im_time, fm_mat_work, dimen_ri, dimen_ri_red, fm_mat_l, fm_mat_q)
...
Definition rpa_util.F:1192
subroutine, public dealloc_im_time(cfm_mo_coeff, index_to_cell_3c, cell_to_index_3c, do_ic_model, do_kpoints_cubic_rpa, do_kpoints_from_gamma, do_ri_sigma_x, has_mat_p_blocks, wkp_w, cfm_mat_q, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, fm_matrix_minv, fm_matrix_minv_vtrunc_minv, fm_mat_ri_global_work, fm_mat_work, mat_dm, mat_l, mat_minvvminv, mat_p_omega, mat_p_omega_kp, t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, mat_work, qs_env)
...
Definition rpa_util.F:1052
subroutine, public remove_scaling_factor_rpa(fm_mat_s, virtual, eigenval_last, homo, omega_old)
...
Definition rpa_util.F:654
subroutine, public calc_fm_mat_s_rpa(fm_mat_s, first_cycle, virtual, eigenval, homo, omega, omega_old)
...
Definition rpa_util.F:704
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Represent a complex full matrix.
keeps the information about the structure of a full matrix
represent a full matrix
Contains information about kpoints.
stores all the informations relevant to an mpi environment