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