(git:b977e33)
rpa_main.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 Routines to calculate RI-RPA energy
10 !> \par History
11 !> 06.2012 created [Mauro Del Ben]
12 !> 04.2015 GW routines added [Jan Wilhelm]
13 !> 10.2015 Cubic-scaling RPA routines added [Jan Wilhelm]
14 !> 10.2018 Cubic-scaling SOS-MP2 added [Frederick Stein]
15 !> 03.2019 Refactoring [Frederick Stein]
16 ! **************************************************************************************************
17 MODULE rpa_main
18  USE bibliography, ONLY: &
19  bates2013, delben2013, delben2015, ren2011, ren2013, wilhelm2016a, wilhelm2016b, &
20  wilhelm2017, wilhelm2018, cite_reference
21  USE bse_main, ONLY: start_bse_calculation
22  USE cp_blacs_env, ONLY: cp_blacs_env_create,&
23  cp_blacs_env_release,&
24  cp_blacs_env_type
25  USE cp_cfm_types, ONLY: cp_cfm_type
26  USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
27  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
28  cp_fm_struct_release,&
29  cp_fm_struct_type
30  USE cp_fm_types, ONLY: cp_fm_create,&
31  cp_fm_release,&
32  cp_fm_set_all,&
33  cp_fm_to_fm,&
34  cp_fm_type
35  USE dbcsr_api, ONLY: dbcsr_add,&
36  dbcsr_clear,&
37  dbcsr_get_info,&
38  dbcsr_p_type,&
39  dbcsr_release,&
40  dbcsr_type
41  USE dbt_api, ONLY: dbt_type
42  USE dgemm_counter_types, ONLY: dgemm_counter_init,&
43  dgemm_counter_type,&
44  dgemm_counter_write
45  USE group_dist_types, ONLY: create_group_dist,&
46  get_group_dist,&
47  group_dist_d1_type,&
48  maxsize,&
49  release_group_dist
50  USE hfx_types, ONLY: block_ind_type,&
51  hfx_compression_type
52  USE input_constants, ONLY: wfc_mm_style_gemm
53  USE kinds, ONLY: dp,&
54  int_8
55  USE kpoint_types, ONLY: kpoint_type
56  USE machine, ONLY: m_flush,&
57  m_memory
58  USE mathconstants, ONLY: pi,&
59  z_zero
60  USE message_passing, ONLY: mp_comm_type,&
61  mp_para_env_release,&
62  mp_para_env_type
63  USE minimax_exp, ONLY: check_exp_minimax_range
64  USE mp2_grids, ONLY: get_clenshaw_grid,&
65  get_minimax_grid
66  USE mp2_laplace, ONLY: sos_mp2_postprocessing
67  USE mp2_ri_grad_util, ONLY: array2fm
68  USE mp2_types, ONLY: mp2_type,&
69  three_dim_real_array,&
70  two_dim_int_array,&
71  two_dim_real_array
72  USE qs_environment_types, ONLY: get_qs_env,&
73  qs_environment_type
74  USE rpa_axk, ONLY: compute_axk_ener
75  USE rpa_grad, ONLY: rpa_grad_copy_q,&
76  rpa_grad_create,&
77  rpa_grad_finalize,&
78  rpa_grad_matrix_operations,&
79  rpa_grad_needed_mem,&
80  rpa_grad_type
81  USE rpa_gw, ONLY: allocate_matrices_gw,&
82  allocate_matrices_gw_im_time,&
83  compute_gw_self_energy,&
84  compute_qp_energies,&
85  compute_w_cubic_gw,&
86  deallocate_matrices_gw,&
87  deallocate_matrices_gw_im_time,&
88  get_fermi_level_offset
89  USE rpa_gw_ic, ONLY: calculate_ic_correction
90  USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_mos,&
91  invert_eps_compute_w_and_erpa_kp
92  USE rpa_im_time, ONLY: compute_mat_p_omega,&
93  zero_mat_p_omega
94  USE rpa_im_time_force_methods, ONLY: calc_laplace_loop_forces,&
95  calc_post_loop_forces,&
96  calc_rpa_loop_forces,&
97  init_im_time_forces,&
98  keep_initial_quad
99  USE rpa_im_time_force_types, ONLY: im_time_force_release,&
100  im_time_force_type
101  USE rpa_util, ONLY: q_trace_and_add_unit_matrix,&
102  alloc_im_time,&
103  calc_mat_q,&
104  compute_erpa_by_freq_int,&
105  contract_p_omega_with_mat_l,&
106  dealloc_im_time,&
107  remove_scaling_factor_rpa
108  USE util, ONLY: get_limit
109 #include "./base/base_uses.f90"
110 
111  IMPLICIT NONE
112 
113  PRIVATE
114 
115  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_main'
116 
117  PUBLIC :: rpa_ri_compute_en
118 
119 CONTAINS
120 
121 ! **************************************************************************************************
122 !> \brief ...
123 !> \param qs_env ...
124 !> \param Erpa ...
125 !> \param mp2_env ...
126 !> \param BIb_C ...
127 !> \param BIb_C_gw ...
128 !> \param BIb_C_bse_ij ...
129 !> \param BIb_C_bse_ab ...
130 !> \param para_env ...
131 !> \param para_env_sub ...
132 !> \param color_sub ...
133 !> \param gd_array ...
134 !> \param gd_B_virtual ...
135 !> \param gd_B_all ...
136 !> \param gd_B_occ_bse ...
137 !> \param gd_B_virt_bse ...
138 !> \param mo_coeff ...
139 !> \param fm_matrix_PQ ...
140 !> \param fm_matrix_L_kpoints ...
141 !> \param fm_matrix_Minv_L_kpoints ...
142 !> \param fm_matrix_Minv ...
143 !> \param fm_matrix_Minv_Vtrunc_Minv ...
144 !> \param kpoints ...
145 !> \param Eigenval ...
146 !> \param nmo ...
147 !> \param homo ...
148 !> \param dimen_RI ...
149 !> \param dimen_RI_red ...
150 !> \param gw_corr_lev_occ ...
151 !> \param gw_corr_lev_virt ...
152 !> \param unit_nr ...
153 !> \param do_ri_sos_laplace_mp2 ...
154 !> \param my_do_gw ...
155 !> \param do_im_time ...
156 !> \param do_bse ...
157 !> \param matrix_s ...
158 !> \param mat_munu ...
159 !> \param mat_P_global ...
160 !> \param t_3c_M ...
161 !> \param t_3c_O ...
162 !> \param t_3c_O_compressed ...
163 !> \param t_3c_O_ind ...
164 !> \param starts_array_mc ...
165 !> \param ends_array_mc ...
166 !> \param starts_array_mc_block ...
167 !> \param ends_array_mc_block ...
168 !> \param calc_forces ...
169 ! **************************************************************************************************
170  SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, &
171  para_env, para_env_sub, color_sub, &
172  gd_array, gd_B_virtual, gd_B_all, gd_B_occ_bse, gd_B_virt_bse, &
173  mo_coeff, fm_matrix_PQ, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
174  fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, kpoints, &
175  Eigenval, nmo, homo, dimen_RI, dimen_RI_red, gw_corr_lev_occ, gw_corr_lev_virt, &
176  unit_nr, do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_bse, matrix_s, &
177  mat_munu, mat_P_global, &
178  t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
179  starts_array_mc, ends_array_mc, &
180  starts_array_mc_block, ends_array_mc_block, calc_forces)
181 
182  TYPE(qs_environment_type), POINTER :: qs_env
183  REAL(kind=dp), INTENT(OUT) :: erpa
184  TYPE(mp2_type), INTENT(INOUT) :: mp2_env
185  TYPE(three_dim_real_array), DIMENSION(:), &
186  INTENT(INOUT) :: bib_c, bib_c_gw
187  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
188  INTENT(INOUT) :: bib_c_bse_ij, bib_c_bse_ab
189  TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
190  INTEGER, INTENT(INOUT) :: color_sub
191  TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_array
192  TYPE(group_dist_d1_type), DIMENSION(:), &
193  INTENT(INOUT) :: gd_b_virtual
194  TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_b_all, gd_b_occ_bse, gd_b_virt_bse
195  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
196  TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_pq
197  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_l_kpoints, &
198  fm_matrix_minv_l_kpoints, &
199  fm_matrix_minv, &
200  fm_matrix_minv_vtrunc_minv
201  TYPE(kpoint_type), POINTER :: kpoints
202  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
203  INTENT(INOUT) :: eigenval
204  INTEGER, INTENT(IN) :: nmo
205  INTEGER, DIMENSION(:), INTENT(IN) :: homo
206  INTEGER, INTENT(IN) :: dimen_ri, dimen_ri_red
207  INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt
208  INTEGER, INTENT(IN) :: unit_nr
209  LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, my_do_gw, &
210  do_im_time, do_bse
211  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
212  TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
213  TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_p_global
214  TYPE(dbt_type) :: t_3c_m
215  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_o
216  TYPE(hfx_compression_type), ALLOCATABLE, &
217  DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_o_compressed
218  TYPE(block_ind_type), ALLOCATABLE, &
219  DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_o_ind
220  INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
221  starts_array_mc_block, &
222  ends_array_mc_block
223  LOGICAL, INTENT(IN) :: calc_forces
224 
225  CHARACTER(LEN=*), PARAMETER :: routinen = 'rpa_ri_compute_en'
226 
227  INTEGER :: best_integ_group_size, best_num_integ_point, color_rpa_group, dimen_homo_square, &
228  dimen_nm_gw, dimen_virt_square, handle, handle2, handle3, ierr, iib, &
229  input_num_integ_groups, integ_group_size, ispin, jjb, min_integ_group_size, &
230  my_ab_comb_bse_end, my_ab_comb_bse_size, my_ab_comb_bse_start, my_group_l_end, &
231  my_group_l_size, my_group_l_start, my_ij_comb_bse_end, my_ij_comb_bse_size, &
232  my_ij_comb_bse_start, my_nm_gw_end, my_nm_gw_size, my_nm_gw_start, ncol_block_mat, &
233  ngroup, nrow_block_mat, nspins, num_integ_group, num_integ_points, pos_integ_group
234  INTEGER(KIND=int_8) :: mem
235  INTEGER, ALLOCATABLE, DIMENSION(:) :: dimen_ia, my_ia_end, my_ia_size, &
236  my_ia_start, virtual
237  LOGICAL :: do_kpoints_from_gamma, do_minimax_quad, &
238  my_open_shell, skip_integ_group_opt
239  REAL(kind=dp) :: allowed_memory, avail_mem, e_range, emax, emin, mem_for_iak, mem_for_qk, &
240  mem_min, mem_per_group, mem_per_rank, mem_per_repl, mem_real
241  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: eigenval_kp
242  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_q, fm_mat_q_gemm, fm_mat_s, &
243  fm_mat_s_gw
244  TYPE(cp_fm_type), DIMENSION(1) :: fm_mat_r_gw, fm_mat_s_ab_bse, &
245  fm_mat_s_ij_bse
246  TYPE(mp_para_env_type), POINTER :: para_env_rpa
247  TYPE(two_dim_real_array), ALLOCATABLE, &
248  DIMENSION(:) :: bib_c_2d, bib_c_2d_gw
249  TYPE(two_dim_real_array), DIMENSION(1) :: bib_c_2d_bse_ab, bib_c_2d_bse_ij
250 
251  CALL timeset(routinen, handle)
252 
253  CALL cite_reference(delben2013)
254  CALL cite_reference(delben2015)
255 
256  IF (mp2_env%ri_rpa%do_ri_axk) THEN
257  CALL cite_reference(bates2013)
258  END IF
259  IF (mp2_env%ri_rpa%do_rse) THEN
260  CALL cite_reference(ren2011)
261  CALL cite_reference(ren2013)
262  END IF
263 
264  IF (my_do_gw) THEN
265  CALL cite_reference(wilhelm2016a)
266  CALL cite_reference(wilhelm2017)
267  CALL cite_reference(wilhelm2018)
268  END IF
269 
270  IF (do_im_time) THEN
271  CALL cite_reference(wilhelm2016b)
272  END IF
273 
274  nspins = SIZE(homo)
275  my_open_shell = (nspins == 2)
276  ALLOCATE (virtual(nspins), dimen_ia(nspins), my_ia_end(nspins), my_ia_start(nspins), my_ia_size(nspins))
277  virtual(:) = nmo - homo(:)
278  dimen_ia(:) = virtual(:)*homo(:)
279 
280  ALLOCATE (eigenval_kp(nmo, 1, nspins))
281  eigenval_kp(:, 1, :) = eigenval(:, :)
282 
283  IF (do_im_time) mp2_env%ri_rpa%minimax_quad = .true.
284  do_minimax_quad = mp2_env%ri_rpa%minimax_quad
285 
286  IF (do_ri_sos_laplace_mp2) THEN
287  num_integ_points = mp2_env%ri_laplace%n_quadrature
288  input_num_integ_groups = mp2_env%ri_laplace%num_integ_groups
289 
290  ! check the range for the minimax approximation
291  e_range = mp2_env%e_range
292  IF (mp2_env%e_range <= 1.0_dp .OR. mp2_env%e_gap <= 0.0_dp) THEN
293  emin = huge(dp)
294  emax = 0.0_dp
295  DO ispin = 1, nspins
296  IF (homo(ispin) > 0) THEN
297  emin = min(emin, 2.0_dp*(eigenval(homo(ispin) + 1, ispin) - eigenval(homo(ispin), ispin)))
298  emax = max(emax, 2.0_dp*(maxval(eigenval(:, ispin)) - minval(eigenval(:, ispin))))
299  END IF
300  END DO
301  e_range = emax/emin
302  END IF
303  IF (e_range < 2.0_dp) e_range = 2.0_dp
304  ierr = 0
305  CALL check_exp_minimax_range(num_integ_points, e_range, ierr)
306  IF (ierr /= 0) THEN
307  jjb = num_integ_points - 1
308  DO iib = 1, jjb
309  num_integ_points = num_integ_points - 1
310  ierr = 0
311  CALL check_exp_minimax_range(num_integ_points, e_range, ierr)
312  IF (ierr == 0) EXIT
313  END DO
314  END IF
315  cpassert(num_integ_points >= 1)
316  ELSE
317  num_integ_points = mp2_env%ri_rpa%rpa_num_quad_points
318  input_num_integ_groups = mp2_env%ri_rpa%rpa_num_integ_groups
319 
320  IF (my_do_gw .AND. do_minimax_quad) THEN
321  IF (num_integ_points > 34) THEN
322  IF (unit_nr > 0) &
323  CALL cp_warn(__location__, &
324  "The required number of quadrature point exceeds the maximum possible in the "// &
325  "Minimax quadrature scheme. The number of quadrature point has been reset to 30.")
326  num_integ_points = 30
327  END IF
328  ELSE
329  IF (do_minimax_quad .AND. num_integ_points > 20) THEN
330  IF (unit_nr > 0) &
331  CALL cp_warn(__location__, &
332  "The required number of quadrature point exceeds the maximum possible in the "// &
333  "Minimax quadrature scheme. The number of quadrature point has been reset to 20.")
334  num_integ_points = 20
335  END IF
336  END IF
337  END IF
338  allowed_memory = mp2_env%mp2_memory
339 
340  CALL get_group_dist(gd_array, color_sub, my_group_l_start, my_group_l_end, my_group_l_size)
341 
342  ngroup = para_env%num_pe/para_env_sub%num_pe
343 
344  ! for imaginary time or periodic GW or BSE, we use all processors for a single frequency/time point
345  IF (do_im_time .OR. mp2_env%ri_g0w0%do_periodic .OR. do_bse) THEN
346 
347  integ_group_size = ngroup
348  best_num_integ_point = num_integ_points
349 
350  ELSE
351 
352  ! Calculate available memory and create integral group according to that
353  ! mem_for_iaK is the memory needed for storing the 3 centre integrals
354  mem_for_iak = real(sum(dimen_ia), kind=dp)*dimen_ri_red*8.0_dp/(1024_dp**2)
355  mem_for_qk = real(dimen_ri_red, kind=dp)*nspins*dimen_ri_red*8.0_dp/(1024_dp**2)
356 
357  CALL m_memory(mem)
358  mem_real = (mem + 1024*1024 - 1)/(1024*1024)
359  CALL para_env%min(mem_real)
360 
361  mem_per_rank = 0.0_dp
362 
363  ! B_ia_P
364  mem_per_repl = mem_for_iak
365  ! Q (regular and for dgemm)
366  mem_per_repl = mem_per_repl + 2.0_dp*mem_for_qk
367 
368  IF (calc_forces) CALL rpa_grad_needed_mem(homo, virtual, dimen_ri_red, mem_per_rank, mem_per_repl, do_ri_sos_laplace_mp2)
369 
370  mem_min = mem_per_repl/para_env%num_pe + mem_per_rank
371 
372  IF (unit_nr > 0) THEN
373  WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum required memory per MPI process:', mem_min, ' MiB'
374  WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Available memory per MPI process:', mem_real, ' MiB'
375  END IF
376 
377  ! Use only the allowed amount of memory
378  mem_real = min(mem_real, allowed_memory)
379  ! For the memory estimate, we require the amount of required memory per replication group and the available memory
380  mem_real = mem_real - mem_per_rank
381 
382  mem_per_group = mem_real*para_env_sub%num_pe
383 
384  ! here we try to find the best rpa/laplace group size
385  skip_integ_group_opt = .false.
386 
387  ! Check the input number of integration groups
388  IF (input_num_integ_groups > 0) THEN
389  IF (mod(num_integ_points, input_num_integ_groups) == 0) THEN
390  IF (mod(ngroup, input_num_integ_groups) == 0) THEN
391  best_integ_group_size = ngroup/input_num_integ_groups
392  best_num_integ_point = num_integ_points/input_num_integ_groups
393  skip_integ_group_opt = .true.
394  ELSE
395  IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Total number of groups not multiple of NUM_INTEG_GROUPS'
396  END IF
397  ELSE
398  IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'NUM_QUAD_POINTS not multiple of NUM_INTEG_GROUPS'
399  END IF
400  END IF
401 
402  IF (.NOT. skip_integ_group_opt) THEN
403  best_integ_group_size = ngroup
404  best_num_integ_point = num_integ_points
405 
406  min_integ_group_size = max(1, ngroup/num_integ_points)
407 
408  integ_group_size = min_integ_group_size - 1
409  DO iib = min_integ_group_size + 1, ngroup
410  integ_group_size = integ_group_size + 1
411 
412  ! check that the ngroup is a multiple of integ_group_size
413  IF (mod(ngroup, integ_group_size) /= 0) cycle
414 
415  ! check for memory
416  avail_mem = integ_group_size*mem_per_group
417  IF (avail_mem < mem_per_repl) cycle
418 
419  ! check that the integration groups have the same size
420  num_integ_group = ngroup/integ_group_size
421  IF (mod(num_integ_points, num_integ_group) /= 0) cycle
422 
423  best_num_integ_point = num_integ_points/num_integ_group
424  best_integ_group_size = integ_group_size
425 
426  EXIT
427 
428  END DO
429  END IF
430 
431  integ_group_size = best_integ_group_size
432 
433  END IF
434 
435  IF (unit_nr > 0 .AND. .NOT. do_im_time) THEN
436  IF (do_ri_sos_laplace_mp2) THEN
437  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
438  "RI_INFO| Group size for laplace numerical integration:", integ_group_size*para_env_sub%num_pe
439  WRITE (unit=unit_nr, fmt="(T3,A)") &
440  "INTEG_INFO| MINIMAX approximation"
441  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
442  "INTEG_INFO| Number of integration points:", num_integ_points
443  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
444  "INTEG_INFO| Number of integration points per Laplace group:", best_num_integ_point
445  ELSE
446  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
447  "RI_INFO| Group size for frequency integration:", integ_group_size*para_env_sub%num_pe
448  IF (do_minimax_quad) THEN
449  WRITE (unit=unit_nr, fmt="(T3,A)") &
450  "INTEG_INFO| MINIMAX quadrature"
451  ELSE
452  WRITE (unit=unit_nr, fmt="(T3,A)") &
453  "INTEG_INFO| Clenshaw-Curtius quadrature"
454  END IF
455  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
456  "INTEG_INFO| Number of integration points:", num_integ_points
457  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
458  "INTEG_INFO| Number of integration points per RPA group:", best_num_integ_point
459  END IF
460  CALL m_flush(unit_nr)
461  END IF
462 
463  num_integ_group = ngroup/integ_group_size
464 
465  pos_integ_group = mod(color_sub, integ_group_size)
466  color_rpa_group = color_sub/integ_group_size
467 
468  CALL timeset(routinen//"_reorder", handle2)
469 
470  ! not necessary for imaginary time
471 
472  ALLOCATE (bib_c_2d(nspins))
473  IF (.NOT. do_im_time) THEN
474 
475  ! reorder the local data in such a way to help the next stage of matrix creation
476  ! now the data inside the group are divided into a ia x K matrix
477  DO ispin = 1, nspins
478  CALL calculate_bib_c_2d(bib_c_2d(ispin)%array, bib_c(ispin)%array, para_env_sub, dimen_ia(ispin), &
479  homo(ispin), virtual(ispin), gd_b_virtual(ispin), &
480  my_ia_size(ispin), my_ia_start(ispin), my_ia_end(ispin), my_group_l_size)
481 
482  DEALLOCATE (bib_c(ispin)%array)
483  CALL release_group_dist(gd_b_virtual(ispin))
484 
485  END DO
486 
487  ! in the GW case, BIb_C_2D_gw is an nm x K matrix, with n: number of corr GW levels, m=nmo
488  IF (my_do_gw) THEN
489  ALLOCATE (bib_c_2d_gw(nspins))
490 
491  CALL timeset(routinen//"_reorder_gw", handle3)
492 
493  dimen_nm_gw = nmo*(gw_corr_lev_occ(1) + gw_corr_lev_virt(1))
494 
495  ! The same for open shell
496  DO ispin = 1, nspins
497  CALL calculate_bib_c_2d(bib_c_2d_gw(ispin)%array, bib_c_gw(ispin)%array, para_env_sub, dimen_nm_gw, &
498  gw_corr_lev_occ(ispin) + gw_corr_lev_virt(ispin), nmo, gd_b_all, &
499  my_nm_gw_size, my_nm_gw_start, my_nm_gw_end, my_group_l_size)
500  DEALLOCATE (bib_c_gw(ispin)%array)
501  END DO
502 
503  CALL release_group_dist(gd_b_all)
504 
505  CALL timestop(handle3)
506 
507  END IF
508  END IF
509 
510  IF (do_bse) THEN
511 
512  CALL timeset(routinen//"_reorder_bse1", handle3)
513 
514  dimen_homo_square = homo(1)**2
515 
516  CALL calculate_bib_c_2d(bib_c_2d_bse_ij(1)%array, bib_c_bse_ij, para_env_sub, dimen_homo_square, &
517  homo(1), homo(1), gd_b_occ_bse, &
518  my_ij_comb_bse_size, my_ij_comb_bse_start, my_ij_comb_bse_end, my_group_l_size)
519 
520  DEALLOCATE (bib_c_bse_ij)
521  CALL release_group_dist(gd_b_occ_bse)
522 
523  CALL timestop(handle3)
524 
525  CALL timeset(routinen//"_reorder_bse2", handle3)
526 
527  dimen_virt_square = virtual(1)**2
528 
529  CALL calculate_bib_c_2d(bib_c_2d_bse_ab(1)%array, bib_c_bse_ab, para_env_sub, dimen_virt_square, &
530  virtual(1), virtual(1), gd_b_virt_bse, &
531  my_ab_comb_bse_size, my_ab_comb_bse_start, my_ab_comb_bse_end, my_group_l_size)
532 
533  DEALLOCATE (bib_c_bse_ab)
534  CALL release_group_dist(gd_b_virt_bse)
535 
536  CALL timestop(handle3)
537 
538  END IF
539 
540  CALL timestop(handle2)
541 
542  IF (num_integ_group > 1) THEN
543  ALLOCATE (para_env_rpa)
544  CALL para_env_rpa%from_split(para_env, color_rpa_group)
545  ELSE
546  para_env_rpa => para_env
547  END IF
548 
549  ! now create the matrices needed for the calculation, Q, S and G
550  ! Q and G will have omega dependence
551 
552  IF (do_im_time) THEN
553  ALLOCATE (fm_mat_q(nspins), fm_mat_q_gemm(1), fm_mat_s(1))
554  ELSE
555  ALLOCATE (fm_mat_q(nspins), fm_mat_q_gemm(nspins), fm_mat_s(nspins))
556  END IF
557 
558  CALL create_integ_mat(bib_c_2d, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
559  dimen_ri_red, dimen_ia, color_rpa_group, &
560  mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
561  my_ia_size, my_ia_start, my_ia_end, &
562  my_group_l_size, my_group_l_start, my_group_l_end, &
563  para_env_rpa, fm_mat_s, nrow_block_mat, ncol_block_mat, &
564  dimen_ia_for_block_size=dimen_ia(1), &
565  do_im_time=do_im_time, fm_mat_q_gemm=fm_mat_q_gemm, fm_mat_q=fm_mat_q, qs_env=qs_env)
566 
567  DEALLOCATE (bib_c_2d, my_ia_end, my_ia_size, my_ia_start)
568 
569  ! for GW, we need other matrix fm_mat_S, always allocate the container to prevent crying compilers
570  ALLOCATE (fm_mat_s_gw(nspins))
571  IF (my_do_gw .AND. .NOT. do_im_time) THEN
572 
573  CALL create_integ_mat(bib_c_2d_gw, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
574  dimen_ri_red, [dimen_nm_gw, dimen_nm_gw], color_rpa_group, &
575  mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
576  [my_nm_gw_size, my_nm_gw_size], [my_nm_gw_start, my_nm_gw_start], [my_nm_gw_end, my_nm_gw_end], &
577  my_group_l_size, my_group_l_start, my_group_l_end, &
578  para_env_rpa, fm_mat_s_gw, nrow_block_mat, ncol_block_mat, &
579  fm_mat_q(1)%matrix_struct%context, fm_mat_q(1)%matrix_struct%context, &
580  fm_mat_q=fm_mat_r_gw)
581  DEALLOCATE (bib_c_2d_gw)
582 
583  END IF
584 
585  ! for Bethe-Salpeter, we need other matrix fm_mat_S
586  IF (do_bse) THEN
587  CALL create_integ_mat(bib_c_2d_bse_ij, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
588  dimen_ri_red, [dimen_homo_square], color_rpa_group, &
589  mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
590  [my_ij_comb_bse_size], [my_ij_comb_bse_start], [my_ij_comb_bse_end], &
591  my_group_l_size, my_group_l_start, my_group_l_end, &
592  para_env_rpa, fm_mat_s_ij_bse, nrow_block_mat, ncol_block_mat, &
593  fm_mat_q(1)%matrix_struct%context, fm_mat_q(1)%matrix_struct%context)
594 
595  CALL create_integ_mat(bib_c_2d_bse_ab, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
596  dimen_ri_red, [dimen_virt_square], color_rpa_group, &
597  mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
598  [my_ab_comb_bse_size], [my_ab_comb_bse_start], [my_ab_comb_bse_end], &
599  my_group_l_size, my_group_l_start, my_group_l_end, &
600  para_env_rpa, fm_mat_s_ab_bse, nrow_block_mat, ncol_block_mat, &
601  fm_mat_q(1)%matrix_struct%context, fm_mat_q(1)%matrix_struct%context)
602 
603  END IF
604 
605  do_kpoints_from_gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
606  IF (do_kpoints_from_gamma) THEN
607  CALL get_bandstruc_and_k_dependent_mos(qs_env, eigenval_kp)
608  END IF
609 
610  ! Now start the RPA calculation
611  ! fm_mo_coeff_occ, fm_mo_coeff_virt will be deallocated here
612  CALL rpa_num_int(qs_env, erpa, mp2_env, para_env, para_env_rpa, para_env_sub, unit_nr, &
613  homo, virtual, dimen_ri, dimen_ri_red, dimen_ia, dimen_nm_gw, &
614  eigenval_kp, num_integ_points, num_integ_group, color_rpa_group, &
615  fm_matrix_pq, fm_mat_s, fm_mat_q_gemm, fm_mat_q, fm_mat_s_gw, fm_mat_r_gw(1), &
616  fm_mat_s_ij_bse(1), fm_mat_s_ab_bse(1), &
617  my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
618  do_minimax_quad, &
619  do_im_time, mo_coeff, &
620  fm_matrix_l_kpoints, fm_matrix_minv_l_kpoints, &
621  fm_matrix_minv, fm_matrix_minv_vtrunc_minv, &
622  mat_munu, mat_p_global, &
623  t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, &
624  starts_array_mc, ends_array_mc, &
625  starts_array_mc_block, ends_array_mc_block, &
626  matrix_s, do_kpoints_from_gamma, kpoints, gd_array, color_sub, &
627  do_ri_sos_laplace_mp2=do_ri_sos_laplace_mp2, calc_forces=calc_forces)
628 
629  CALL release_group_dist(gd_array)
630 
631  IF (num_integ_group > 1) CALL mp_para_env_release(para_env_rpa)
632 
633  IF (.NOT. do_im_time) THEN
634  CALL cp_fm_release(fm_mat_q_gemm)
635  CALL cp_fm_release(fm_mat_s)
636  END IF
637  CALL cp_fm_release(fm_mat_q)
638 
639  IF (my_do_gw .AND. .NOT. do_im_time) THEN
640  CALL cp_fm_release(fm_mat_s_gw)
641  CALL cp_fm_release(fm_mat_r_gw(1))
642  END IF
643 
644  IF (do_bse) THEN
645  CALL cp_fm_release(fm_mat_s_ij_bse(1))
646  CALL cp_fm_release(fm_mat_s_ab_bse(1))
647  END IF
648 
649  IF (mp2_env%ri_rpa%do_ri_axk) THEN
650  CALL dbcsr_release(mp2_env%ri_rpa%mo_coeff_o)
651  DEALLOCATE (mp2_env%ri_rpa%mo_coeff_o)
652  CALL dbcsr_release(mp2_env%ri_rpa%mo_coeff_v)
653  DEALLOCATE (mp2_env%ri_rpa%mo_coeff_v)
654  END IF
655 
656  CALL timestop(handle)
657 
658  END SUBROUTINE rpa_ri_compute_en
659 
660 ! **************************************************************************************************
661 !> \brief reorder the local data in such a way to help the next stage of matrix creation;
662 !> now the data inside the group are divided into a ia x K matrix (BIb_C_2D);
663 !> Subroutine created to avoid massive double coding
664 !> \param BIb_C_2D ...
665 !> \param BIb_C ...
666 !> \param para_env_sub ...
667 !> \param dimen_ia ...
668 !> \param homo ...
669 !> \param virtual ...
670 !> \param gd_B_virtual ...
671 !> \param my_ia_size ...
672 !> \param my_ia_start ...
673 !> \param my_ia_end ...
674 !> \param my_group_L_size ...
675 !> \author Jan Wilhelm, 03/2015
676 ! **************************************************************************************************
677  SUBROUTINE calculate_bib_c_2d(BIb_C_2D, BIb_C, para_env_sub, dimen_ia, homo, virtual, &
678  gd_B_virtual, &
679  my_ia_size, my_ia_start, my_ia_end, my_group_L_size)
680 
681  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
682  INTENT(OUT) :: bib_c_2d
683  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
684  INTENT(IN) :: bib_c
685  TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
686  INTEGER, INTENT(IN) :: dimen_ia, homo, virtual
687  TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_b_virtual
688  INTEGER :: my_ia_size, my_ia_start, my_ia_end, &
689  my_group_l_size
690 
691  INTEGER, PARAMETER :: occ_chunk = 128
692 
693  INTEGER :: ia_global, iib, itmp(2), jjb, my_b_size, my_b_virtual_start, occ_high, occ_low, &
694  proc_receive, proc_send, proc_shift, rec_b_size, rec_b_virtual_end, rec_b_virtual_start
695  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), TARGET :: bib_c_rec_1d
696  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: bib_c_rec
697 
698  itmp = get_limit(dimen_ia, para_env_sub%num_pe, para_env_sub%mepos)
699  my_ia_start = itmp(1)
700  my_ia_end = itmp(2)
701  my_ia_size = my_ia_end - my_ia_start + 1
702 
703  CALL get_group_dist(gd_b_virtual, para_env_sub%mepos, sizes=my_b_size, starts=my_b_virtual_start)
704 
705  ! reorder data
706  ALLOCATE (bib_c_2d(my_group_l_size, my_ia_size))
707 
708 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,ia_global) &
709 !$OMP SHARED(homo,my_B_size,virtual,my_B_virtual_start,my_ia_start,my_ia_end,BIb_C,BIb_C_2D,&
710 !$OMP my_group_L_size)
711  DO iib = 1, homo
712  DO jjb = 1, my_b_size
713  ia_global = (iib - 1)*virtual + my_b_virtual_start + jjb - 1
714  IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
715  bib_c_2d(1:my_group_l_size, ia_global - my_ia_start + 1) = bib_c(1:my_group_l_size, jjb, iib)
716  END IF
717  END DO
718  END DO
719 
720  IF (para_env_sub%num_pe > 1) THEN
721  ALLOCATE (bib_c_rec_1d(int(my_group_l_size, int_8)*maxsize(gd_b_virtual)*min(homo, occ_chunk)))
722  DO proc_shift = 1, para_env_sub%num_pe - 1
723  proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
724  proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
725 
726  CALL get_group_dist(gd_b_virtual, proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
727 
728  ! do this in chunks to avoid high memory overhead
729  DO occ_low = 1, homo, occ_chunk
730  occ_high = min(homo, occ_low + occ_chunk - 1)
731  bib_c_rec(1:my_group_l_size, 1:rec_b_size, 1:occ_high - occ_low + 1) => &
732  bib_c_rec_1d(1:int(my_group_l_size, int_8)*rec_b_size*(occ_high - occ_low + 1))
733  CALL para_env_sub%sendrecv(bib_c(:, :, occ_low:occ_high), proc_send, &
734  bib_c_rec(:, :, 1:occ_high - occ_low + 1), proc_receive)
735 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,ia_global) &
736 !$OMP SHARED(occ_low,occ_high,rec_B_size,virtual,rec_B_virtual_start,my_ia_start,my_ia_end,BIb_C_rec,BIb_C_2D,&
737 !$OMP my_group_L_size)
738  DO iib = occ_low, occ_high
739  DO jjb = 1, rec_b_size
740  ia_global = (iib - 1)*virtual + rec_b_virtual_start + jjb - 1
741  IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
742  bib_c_2d(1:my_group_l_size, ia_global - my_ia_start + 1) = bib_c_rec(1:my_group_l_size, jjb, iib - occ_low + 1)
743  END IF
744  END DO
745  END DO
746  END DO
747 
748  END DO
749  DEALLOCATE (bib_c_rec_1d)
750  END IF
751 
752  END SUBROUTINE calculate_bib_c_2d
753 
754 ! **************************************************************************************************
755 !> \brief ...
756 !> \param BIb_C_2D ...
757 !> \param para_env ...
758 !> \param para_env_sub ...
759 !> \param color_sub ...
760 !> \param ngroup ...
761 !> \param integ_group_size ...
762 !> \param dimen_RI ...
763 !> \param dimen_ia ...
764 !> \param color_rpa_group ...
765 !> \param ext_row_block_size ...
766 !> \param ext_col_block_size ...
767 !> \param unit_nr ...
768 !> \param my_ia_size ...
769 !> \param my_ia_start ...
770 !> \param my_ia_end ...
771 !> \param my_group_L_size ...
772 !> \param my_group_L_start ...
773 !> \param my_group_L_end ...
774 !> \param para_env_RPA ...
775 !> \param fm_mat_S ...
776 !> \param nrow_block_mat ...
777 !> \param ncol_block_mat ...
778 !> \param blacs_env_ext ...
779 !> \param blacs_env_ext_S ...
780 !> \param dimen_ia_for_block_size ...
781 !> \param do_im_time ...
782 !> \param fm_mat_Q_gemm ...
783 !> \param fm_mat_Q ...
784 !> \param qs_env ...
785 ! **************************************************************************************************
786  SUBROUTINE create_integ_mat(BIb_C_2D, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
787  dimen_RI, dimen_ia, color_rpa_group, &
788  ext_row_block_size, ext_col_block_size, unit_nr, &
789  my_ia_size, my_ia_start, my_ia_end, &
790  my_group_L_size, my_group_L_start, my_group_L_end, &
791  para_env_RPA, fm_mat_S, nrow_block_mat, ncol_block_mat, &
792  blacs_env_ext, blacs_env_ext_S, dimen_ia_for_block_size, &
793  do_im_time, fm_mat_Q_gemm, fm_mat_Q, qs_env)
794 
795  TYPE(two_dim_real_array), DIMENSION(:), &
796  INTENT(INOUT) :: bib_c_2d
797  TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_sub
798  INTEGER, INTENT(IN) :: color_sub, ngroup, integ_group_size, &
799  dimen_ri
800  INTEGER, DIMENSION(:), INTENT(IN) :: dimen_ia
801  INTEGER, INTENT(IN) :: color_rpa_group, ext_row_block_size, &
802  ext_col_block_size, unit_nr
803  INTEGER, DIMENSION(:), INTENT(IN) :: my_ia_size, my_ia_start, my_ia_end
804  INTEGER, INTENT(IN) :: my_group_l_size, my_group_l_start, &
805  my_group_l_end
806  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_rpa
807  TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_s
808  INTEGER, INTENT(INOUT) :: nrow_block_mat, ncol_block_mat
809  TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: blacs_env_ext, blacs_env_ext_s
810  INTEGER, INTENT(IN), OPTIONAL :: dimen_ia_for_block_size
811  LOGICAL, INTENT(IN), OPTIONAL :: do_im_time
812  TYPE(cp_fm_type), DIMENSION(:), OPTIONAL :: fm_mat_q_gemm, fm_mat_q
813  TYPE(qs_environment_type), INTENT(IN), OPTIONAL, &
814  POINTER :: qs_env
815 
816  CHARACTER(LEN=*), PARAMETER :: routinen = 'create_integ_mat'
817 
818  INTEGER :: col_row_proc_ratio, grid_2d(2), handle, &
819  iproc, iproc_col, iproc_row, ispin, &
820  mepos_in_rpa_group
821  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: group_grid_2_mepos
822  LOGICAL :: my_blacs_ext, my_blacs_s_ext, &
823  my_do_im_time
824  TYPE(cp_blacs_env_type), POINTER :: blacs_env, blacs_env_q
825  TYPE(cp_fm_struct_type), POINTER :: fm_struct
826  TYPE(group_dist_d1_type) :: gd_ia, gd_l
827 
828  CALL timeset(routinen, handle)
829 
830  cpassert(PRESENT(blacs_env_ext) .OR. PRESENT(dimen_ia_for_block_size))
831 
832  my_blacs_ext = .false.
833  IF (PRESENT(blacs_env_ext)) my_blacs_ext = .true.
834 
835  my_blacs_s_ext = .false.
836  IF (PRESENT(blacs_env_ext_s)) my_blacs_s_ext = .true.
837 
838  my_do_im_time = .false.
839  IF (PRESENT(do_im_time)) my_do_im_time = do_im_time
840 
841  NULLIFY (blacs_env)
842  ! create the RPA blacs env
843  IF (my_blacs_s_ext) THEN
844  blacs_env => blacs_env_ext_s
845  ELSE
846  IF (para_env_rpa%num_pe > 1) THEN
847  col_row_proc_ratio = max(1, dimen_ia_for_block_size/dimen_ri)
848 
849  iproc_col = min(max(int(sqrt(real(para_env_rpa%num_pe*col_row_proc_ratio, kind=dp))), 1), para_env_rpa%num_pe) + 1
850  DO iproc = 1, para_env_rpa%num_pe
851  iproc_col = iproc_col - 1
852  IF (mod(para_env_rpa%num_pe, iproc_col) == 0) EXIT
853  END DO
854 
855  iproc_row = para_env_rpa%num_pe/iproc_col
856  grid_2d(1) = iproc_row
857  grid_2d(2) = iproc_col
858  ELSE
859  grid_2d = 1
860  END IF
861  CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_rpa, grid_2d=grid_2d)
862 
863  IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
864  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
865  "MATRIX_INFO| Number row processes:", grid_2d(1)
866  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
867  "MATRIX_INFO| Number column processes:", grid_2d(2)
868  END IF
869 
870  ! define the block_size for the row
871  IF (ext_row_block_size > 0) THEN
872  nrow_block_mat = ext_row_block_size
873  ELSE
874  nrow_block_mat = max(1, dimen_ri/grid_2d(1)/2)
875  END IF
876 
877  ! define the block_size for the column
878  IF (ext_col_block_size > 0) THEN
879  ncol_block_mat = ext_col_block_size
880  ELSE
881  ncol_block_mat = max(1, dimen_ia_for_block_size/grid_2d(2)/2)
882  END IF
883 
884  IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
885  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
886  "MATRIX_INFO| Row block size:", nrow_block_mat
887  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
888  "MATRIX_INFO| Column block size:", ncol_block_mat
889  END IF
890  END IF
891 
892  IF (.NOT. my_do_im_time) THEN
893  DO ispin = 1, SIZE(bib_c_2d)
894  NULLIFY (fm_struct)
895  IF (my_blacs_ext) THEN
896  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_ri, &
897  ncol_global=dimen_ia(ispin), para_env=para_env_rpa)
898  ELSE
899  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_ri, &
900  ncol_global=dimen_ia(ispin), para_env=para_env_rpa, &
901  nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.true.)
902 
903  END IF ! external blacs_env
904 
905  CALL create_group_dist(gd_ia, my_ia_start(ispin), my_ia_end(ispin), my_ia_size(ispin), para_env_rpa)
906  CALL create_group_dist(gd_l, my_group_l_start, my_group_l_end, my_group_l_size, para_env_rpa)
907 
908  ! create the info array
909 
910  mepos_in_rpa_group = mod(color_sub, integ_group_size)
911  ALLOCATE (group_grid_2_mepos(0:integ_group_size - 1, 0:para_env_sub%num_pe - 1))
912  group_grid_2_mepos = 0
913  group_grid_2_mepos(mepos_in_rpa_group, para_env_sub%mepos) = para_env_rpa%mepos
914  CALL para_env_rpa%sum(group_grid_2_mepos)
915 
916  CALL array2fm(bib_c_2d(ispin)%array, fm_struct, my_group_l_start, my_group_l_end, &
917  my_ia_start(ispin), my_ia_end(ispin), gd_l, gd_ia, &
918  group_grid_2_mepos, ngroup, para_env_sub%num_pe, fm_mat_s(ispin), &
919  integ_group_size, color_rpa_group)
920 
921  DEALLOCATE (group_grid_2_mepos)
922  CALL cp_fm_struct_release(fm_struct)
923 
924  ! deallocate the info array
925  CALL release_group_dist(gd_l)
926  CALL release_group_dist(gd_ia)
927 
928  ! sum the local data across processes belonging to different RPA group.
929  IF (para_env_rpa%num_pe /= para_env%num_pe) THEN
930  block
931  TYPE(mp_comm_type) :: comm_exchange
932  comm_exchange = fm_mat_s(ispin)%matrix_struct%context%interconnect(para_env)
933  CALL comm_exchange%sum(fm_mat_s(ispin)%local_data)
934  CALL comm_exchange%free()
935  END block
936  END IF
937  END DO
938  END IF
939 
940  IF (PRESENT(fm_mat_q_gemm) .AND. .NOT. my_do_im_time) THEN
941  ! create the Q matrix dimen_RIxdimen_RI where the result of the mat-mat-mult will be stored
942  NULLIFY (fm_struct)
943  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_ri, &
944  ncol_global=dimen_ri, para_env=para_env_rpa, &
945  nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.true.)
946  DO ispin = 1, SIZE(fm_mat_q_gemm)
947  CALL cp_fm_create(fm_mat_q_gemm(ispin), fm_struct, name="fm_mat_Q_gemm")
948  END DO
949  CALL cp_fm_struct_release(fm_struct)
950  END IF
951 
952  IF (PRESENT(fm_mat_q)) THEN
953  NULLIFY (blacs_env_q)
954  IF (my_blacs_ext) THEN
955  blacs_env_q => blacs_env_ext
956  ELSE IF (para_env_rpa%num_pe == para_env%num_pe .AND. PRESENT(qs_env)) THEN
957  CALL get_qs_env(qs_env, blacs_env=blacs_env_q)
958  ELSE
959  CALL cp_blacs_env_create(blacs_env=blacs_env_q, para_env=para_env_rpa)
960  END IF
961  NULLIFY (fm_struct)
962  CALL cp_fm_struct_create(fm_struct, context=blacs_env_q, nrow_global=dimen_ri, &
963  ncol_global=dimen_ri, para_env=para_env_rpa)
964  DO ispin = 1, SIZE(fm_mat_q)
965  CALL cp_fm_create(fm_mat_q(ispin), fm_struct, name="fm_mat_Q")
966  END DO
967 
968  CALL cp_fm_struct_release(fm_struct)
969 
970  IF (.NOT. (my_blacs_ext .OR. (para_env_rpa%num_pe == para_env%num_pe .AND. PRESENT(qs_env)))) &
971  CALL cp_blacs_env_release(blacs_env_q)
972  END IF
973 
974  ! release blacs_env
975  IF (.NOT. my_blacs_s_ext) CALL cp_blacs_env_release(blacs_env)
976 
977  CALL timestop(handle)
978 
979  END SUBROUTINE create_integ_mat
980 
981 ! **************************************************************************************************
982 !> \brief ...
983 !> \param qs_env ...
984 !> \param Erpa ...
985 !> \param mp2_env ...
986 !> \param para_env ...
987 !> \param para_env_RPA ...
988 !> \param para_env_sub ...
989 !> \param unit_nr ...
990 !> \param homo ...
991 !> \param virtual ...
992 !> \param dimen_RI ...
993 !> \param dimen_RI_red ...
994 !> \param dimen_ia ...
995 !> \param dimen_nm_gw ...
996 !> \param Eigenval ...
997 !> \param num_integ_points ...
998 !> \param num_integ_group ...
999 !> \param color_rpa_group ...
1000 !> \param fm_matrix_PQ ...
1001 !> \param fm_mat_S ...
1002 !> \param fm_mat_Q_gemm ...
1003 !> \param fm_mat_Q ...
1004 !> \param fm_mat_S_gw ...
1005 !> \param fm_mat_R_gw ...
1006 !> \param fm_mat_S_ij_bse ...
1007 !> \param fm_mat_S_ab_bse ...
1008 !> \param my_do_gw ...
1009 !> \param do_bse ...
1010 !> \param gw_corr_lev_occ ...
1011 !> \param gw_corr_lev_virt ...
1012 !> \param do_minimax_quad ...
1013 !> \param do_im_time ...
1014 !> \param mo_coeff ...
1015 !> \param fm_matrix_L_kpoints ...
1016 !> \param fm_matrix_Minv_L_kpoints ...
1017 !> \param fm_matrix_Minv ...
1018 !> \param fm_matrix_Minv_Vtrunc_Minv ...
1019 !> \param mat_munu ...
1020 !> \param mat_P_global ...
1021 !> \param t_3c_M ...
1022 !> \param t_3c_O ...
1023 !> \param t_3c_O_compressed ...
1024 !> \param t_3c_O_ind ...
1025 !> \param starts_array_mc ...
1026 !> \param ends_array_mc ...
1027 !> \param starts_array_mc_block ...
1028 !> \param ends_array_mc_block ...
1029 !> \param matrix_s ...
1030 !> \param do_kpoints_from_Gamma ...
1031 !> \param kpoints ...
1032 !> \param gd_array ...
1033 !> \param color_sub ...
1034 !> \param do_ri_sos_laplace_mp2 ...
1035 !> \param calc_forces ...
1036 ! **************************************************************************************************
1037  SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_sub, unit_nr, &
1038  homo, virtual, dimen_RI, dimen_RI_red, dimen_ia, dimen_nm_gw, &
1039  Eigenval, num_integ_points, num_integ_group, color_rpa_group, &
1040  fm_matrix_PQ, fm_mat_S, fm_mat_Q_gemm, fm_mat_Q, fm_mat_S_gw, fm_mat_R_gw, &
1041  fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
1042  my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
1043  do_minimax_quad, do_im_time, mo_coeff, &
1044  fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
1045  fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
1046  mat_munu, mat_P_global, &
1047  t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1048  starts_array_mc, ends_array_mc, &
1049  starts_array_mc_block, ends_array_mc_block, &
1050  matrix_s, do_kpoints_from_Gamma, kpoints, gd_array, color_sub, &
1051  do_ri_sos_laplace_mp2, calc_forces)
1052 
1053  TYPE(qs_environment_type), POINTER :: qs_env
1054  REAL(kind=dp), INTENT(OUT) :: erpa
1055  TYPE(mp2_type) :: mp2_env
1056  TYPE(mp_para_env_type), POINTER :: para_env, para_env_rpa, para_env_sub
1057  INTEGER, INTENT(IN) :: unit_nr
1058  INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
1059  INTEGER, INTENT(IN) :: dimen_ri, dimen_ri_red
1060  INTEGER, DIMENSION(:), INTENT(IN) :: dimen_ia
1061  INTEGER, INTENT(IN) :: dimen_nm_gw
1062  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1063  INTENT(INOUT) :: eigenval
1064  INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, &
1065  color_rpa_group
1066  TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_pq
1067  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_s, fm_mat_q_gemm, fm_mat_q, &
1068  fm_mat_s_gw
1069  TYPE(cp_fm_type), INTENT(IN) :: fm_mat_r_gw, fm_mat_s_ij_bse, &
1070  fm_mat_s_ab_bse
1071  LOGICAL, INTENT(IN) :: my_do_gw, do_bse
1072  INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt
1073  LOGICAL, INTENT(IN) :: do_minimax_quad, do_im_time
1074  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1075  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_l_kpoints, &
1076  fm_matrix_minv_l_kpoints, &
1077  fm_matrix_minv, &
1078  fm_matrix_minv_vtrunc_minv
1079  TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
1080  TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_p_global
1081  TYPE(dbt_type), INTENT(INOUT) :: t_3c_m
1082  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
1083  INTENT(INOUT) :: t_3c_o
1084  TYPE(hfx_compression_type), ALLOCATABLE, &
1085  DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_o_compressed
1086  TYPE(block_ind_type), ALLOCATABLE, &
1087  DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_o_ind
1088  INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
1089  starts_array_mc_block, &
1090  ends_array_mc_block
1091  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1092  LOGICAL :: do_kpoints_from_gamma
1093  TYPE(kpoint_type), POINTER :: kpoints
1094  TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1095  INTEGER, INTENT(IN) :: color_sub
1096  LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, calc_forces
1097 
1098  CHARACTER(LEN=*), PARAMETER :: routinen = 'rpa_num_int'
1099 
1100  COMPLEX(KIND=dp), ALLOCATABLE, &
1101  DIMENSION(:, :, :, :) :: vec_sigma_c_gw
1102  INTEGER :: count_ev_sc_gw, cut_memory, group_size_p, gw_corr_lev_tot, handle, handle3, i, &
1103  ikp_local, ispin, iter_evgw, iter_sc_gw0, j, jquad, min_bsize, mm_style, nkp, &
1104  nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, pspin, qspin, &
1105  size_p
1106  INTEGER(int_8) :: dbcsr_nflop
1107  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c
1108  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c
1109  INTEGER, DIMENSION(:), POINTER :: col_blk_size, prim_blk_sizes, &
1110  ri_blk_sizes
1111  LOGICAL :: do_apply_ic_corr_to_gw, do_gw_im_time, do_ic_model, do_kpoints_cubic_rpa, &
1112  do_periodic, do_print, do_ri_sigma_x, exit_ev_gw, first_cycle, &
1113  first_cycle_periodic_correction, my_open_shell, print_ic_values
1114  LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: has_mat_p_blocks
1115  REAL(kind=dp) :: a_scaling, alpha, dbcsr_time, e_axk, e_axk_corr, eps_filter, &
1116  eps_filter_im_time, ext_scaling, fermi_level_offset, fermi_level_offset_input, &
1117  my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old
1118  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, e_fermi, tau_tj, tau_wj, tj, &
1119  trace_qomega, vec_omega_fit_gw, wj, &
1120  wkp_w
1121  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_w_gw, weights_cos_tf_t_to_w, &
1122  weights_cos_tf_w_to_t, &
1123  weights_sin_tf_t_to_w
1124  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: eigenval_last, eigenval_scf, &
1125  vec_sigma_x_gw
1126  TYPE(cp_cfm_type) :: cfm_mat_q
1127  TYPE(cp_fm_type) :: fm_mat_q_static_bse, fm_mat_q_static_bse_gemm, fm_mat_ri_global_work, &
1128  fm_mat_s_ia_bse, fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1129  fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau
1130  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_s_gw_work, fm_mat_w, &
1131  fm_mo_coeff_occ, fm_mo_coeff_virt
1132  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_l_kpoints, fm_mat_minv_l_kpoints
1133  TYPE(dbcsr_p_type) :: mat_dm, mat_l, mat_m_p_munu_occ, &
1134  mat_m_p_munu_virt, mat_minvvminv
1135  TYPE(dbcsr_p_type), ALLOCATABLE, &
1136  DIMENSION(:, :, :) :: mat_p_omega
1137  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_im_mo_mo, &
1138  matrix_berry_re_mo_mo
1139  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_p_omega_kp
1140  TYPE(dbcsr_type), POINTER :: mat_w, mat_work
1141  TYPE(dbt_type) :: t_3c_overl_int_ao_mo
1142  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_3c_overl_int_gw_ao, &
1143  t_3c_overl_int_gw_ri, &
1144  t_3c_overl_nnp_ic, &
1145  t_3c_overl_nnp_ic_reflected
1146  TYPE(dgemm_counter_type) :: dgemm_counter
1147  TYPE(hfx_compression_type), ALLOCATABLE, &
1148  DIMENSION(:) :: t_3c_o_mo_compressed
1149  TYPE(im_time_force_type) :: force_data
1150  TYPE(rpa_grad_type) :: rpa_grad
1151  TYPE(two_dim_int_array), ALLOCATABLE, DIMENSION(:) :: t_3c_o_mo_ind
1152 
1153  CALL timeset(routinen, handle)
1154 
1155  nspins = SIZE(homo)
1156  nmo = homo(1) + virtual(1)
1157  my_open_shell = (nspins == 2)
1158 
1159  do_gw_im_time = my_do_gw .AND. do_im_time
1160  do_ri_sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x
1161  do_ic_model = mp2_env%ri_g0w0%do_ic_model
1162  print_ic_values = mp2_env%ri_g0w0%print_ic_values
1163  do_periodic = mp2_env%ri_g0w0%do_periodic
1164  do_kpoints_cubic_rpa = mp2_env%ri_rpa_im_time%do_im_time_kpoints
1165 
1166  ! For SOS-MP2 only gemm is implemented
1167  mm_style = wfc_mm_style_gemm
1168  IF (.NOT. do_ri_sos_laplace_mp2) mm_style = mp2_env%ri_rpa%mm_style
1169 
1170  IF (my_do_gw) THEN
1171  ext_scaling = 0.2_dp
1172  omega_max_fit = mp2_env%ri_g0w0%omega_max_fit
1173  fermi_level_offset_input = mp2_env%ri_g0w0%fermi_level_offset
1174  iter_evgw = mp2_env%ri_g0w0%iter_evGW
1175  iter_sc_gw0 = mp2_env%ri_g0w0%iter_sc_GW0
1176  IF ((.NOT. do_im_time)) THEN
1177  IF (iter_sc_gw0 .NE. 1 .AND. iter_evgw .NE. 1) cpabort("Mixed scGW0/evGW not implemented.")
1178  ! in case of scGW0 with the N^4 algorithm, we use the evGW code but use the DFT eigenvalues for W
1179  IF (iter_sc_gw0 .NE. 1) iter_evgw = iter_sc_gw0
1180  END IF
1181  ELSE
1182  ext_scaling = 0.0_dp
1183  iter_evgw = 1
1184  iter_sc_gw0 = 1
1185  END IF
1186 
1187  IF (do_kpoints_cubic_rpa .AND. do_ri_sos_laplace_mp2) THEN
1188  cpabort("RI-SOS-Laplace-MP2 with k-point-sampling is not implemented.")
1189  END IF
1190 
1191  do_apply_ic_corr_to_gw = .false.
1192  IF (mp2_env%ri_g0w0%ic_corr_list(1)%array(1) > 0.0_dp) do_apply_ic_corr_to_gw = .true.
1193 
1194  IF (do_im_time) THEN
1195  cpassert(do_minimax_quad .OR. do_ri_sos_laplace_mp2)
1196 
1197  group_size_p = mp2_env%ri_rpa_im_time%group_size_P
1198  cut_memory = mp2_env%ri_rpa_im_time%cut_memory
1199  eps_filter = mp2_env%ri_rpa_im_time%eps_filter
1200  eps_filter_im_time = mp2_env%ri_rpa_im_time%eps_filter* &
1201  mp2_env%ri_rpa_im_time%eps_filter_factor
1202 
1203  min_bsize = mp2_env%ri_rpa_im_time%min_bsize
1204 
1205  CALL alloc_im_time(qs_env, para_env, dimen_ri, dimen_ri_red, &
1206  num_integ_points, nspins, fm_mat_q(1), fm_mo_coeff_occ, fm_mo_coeff_virt, &
1207  fm_matrix_minv_l_kpoints, fm_matrix_l_kpoints, mat_p_global, &
1208  t_3c_o, matrix_s, kpoints, eps_filter_im_time, &
1209  cut_memory, nkp, num_cells_dm, num_3c_repl, &
1210  size_p, ikp_local, &
1211  index_to_cell_3c, &
1212  cell_to_index_3c, &
1213  col_blk_size, &
1214  do_ic_model, do_kpoints_cubic_rpa, &
1215  do_kpoints_from_gamma, do_ri_sigma_x, my_open_shell, &
1216  has_mat_p_blocks, wkp_w, &
1217  cfm_mat_q, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, &
1218  fm_mat_ri_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, &
1219  fm_mo_coeff_virt_scaled, mat_dm, mat_l, mat_m_p_munu_occ, mat_m_p_munu_virt, &
1220  mat_minvvminv, mat_p_omega, mat_p_omega_kp, mat_work, mo_coeff, &
1221  fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, homo, nmo)
1222 
1223  IF (calc_forces) CALL init_im_time_forces(force_data, fm_matrix_pq, t_3c_m, unit_nr, mp2_env, qs_env)
1224 
1225  IF (my_do_gw) THEN
1226 
1227  CALL dbcsr_get_info(mat_p_global%matrix, &
1228  row_blk_size=ri_blk_sizes)
1229 
1230  CALL dbcsr_get_info(matrix_s(1)%matrix, &
1231  row_blk_size=prim_blk_sizes)
1232 
1233  gw_corr_lev_tot = gw_corr_lev_occ(1) + gw_corr_lev_virt(1)
1234 
1235  IF (.NOT. do_kpoints_cubic_rpa) THEN
1236  CALL allocate_matrices_gw_im_time(gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, &
1237  num_integ_points, unit_nr, &
1238  ri_blk_sizes, do_ic_model, &
1239  para_env, fm_mat_w, fm_mat_q(1), &
1240  mo_coeff, &
1241  t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, &
1242  t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, &
1243  starts_array_mc, ends_array_mc, &
1244  t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, &
1245  matrix_s, mat_w, t_3c_o, &
1246  t_3c_o_compressed, t_3c_o_ind, &
1247  qs_env)
1248 
1249  END IF
1250  END IF
1251 
1252  END IF
1253 
1254  IF (do_ic_model) THEN
1255  ! image charge model only implemented for cubic scaling GW
1256  cpassert(do_gw_im_time)
1257  cpassert(.NOT. do_periodic)
1258  IF (cut_memory .NE. 1) cpabort("For IC, use MEMORY_CUT 1 in the LOW_SCALING section.")
1259  END IF
1260 
1261  ALLOCATE (e_fermi(nspins))
1262  IF (do_minimax_quad .OR. do_ri_sos_laplace_mp2) THEN
1263  do_print = .NOT. do_ic_model
1264  CALL get_minimax_grid(para_env, unit_nr, homo, eigenval, num_integ_points, do_im_time, &
1265  do_ri_sos_laplace_mp2, do_print, &
1266  tau_tj, tau_wj, qs_env, do_gw_im_time, &
1267  do_kpoints_cubic_rpa, e_fermi(1), tj, wj, &
1268  weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
1269  qs_env%mp2_env%ri_g0w0%regularization_minimax)
1270 
1271  !For sos_laplace_mp2 and low-scaling RPA, potentially need to store/retrieve the initial weights
1272  IF (qs_env%mp2_env%ri_rpa_im_time%keep_quad) THEN
1273  CALL keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, &
1274  weights_cos_tf_w_to_t, do_ri_sos_laplace_mp2, do_im_time, &
1275  num_integ_points, unit_nr, qs_env)
1276  END IF
1277  ELSE
1278  IF (calc_forces) cpabort("Forces with Clenshaw-Curtis grid not implemented.")
1279  CALL get_clenshaw_grid(para_env, para_env_rpa, unit_nr, homo, virtual, eigenval, num_integ_points, &
1280  num_integ_group, color_rpa_group, fm_mat_s, my_do_gw, &
1281  ext_scaling, a_scaling, tj, wj)
1282  END IF
1283 
1284  ! This array is needed for RPA
1285  IF (.NOT. do_ri_sos_laplace_mp2) THEN
1286  ALLOCATE (trace_qomega(dimen_ri_red))
1287  END IF
1288 
1289  IF (do_ri_sos_laplace_mp2 .AND. .NOT. do_im_time) THEN
1290  alpha = 1.0_dp
1291  ELSE IF (my_open_shell .OR. do_ri_sos_laplace_mp2) THEN
1292  alpha = 2.0_dp
1293  ELSE
1294  alpha = 4.0_dp
1295  END IF
1296 
1297  IF (my_do_gw) THEN
1298  CALL allocate_matrices_gw(vec_sigma_c_gw, color_rpa_group, dimen_nm_gw, &
1299  gw_corr_lev_occ, gw_corr_lev_virt, homo, &
1300  nmo, num_integ_group, num_integ_points, unit_nr, &
1301  gw_corr_lev_tot, num_fit_points, omega_max_fit, &
1302  do_minimax_quad, do_periodic, do_ri_sigma_x,.NOT. do_im_time, &
1303  first_cycle_periodic_correction, &
1304  a_scaling, eigenval, tj, vec_omega_fit_gw, vec_sigma_x_gw, &
1305  delta_corr, eigenval_last, eigenval_scf, vec_w_gw, &
1306  fm_mat_s_gw, fm_mat_s_gw_work, &
1307  para_env, mp2_env, kpoints, nkp, nkp_self_energy, &
1308  do_kpoints_cubic_rpa, do_kpoints_from_gamma)
1309 
1310  IF (do_bse) THEN
1311 
1312  CALL cp_fm_create(fm_mat_q_static_bse_gemm, fm_mat_q_gemm(1)%matrix_struct)
1313  CALL cp_fm_to_fm(fm_mat_q_gemm(1), fm_mat_q_static_bse_gemm)
1314  CALL cp_fm_set_all(fm_mat_q_static_bse_gemm, 0.0_dp)
1315 
1316  CALL cp_fm_create(fm_mat_q_static_bse, fm_mat_q(1)%matrix_struct)
1317  CALL cp_fm_to_fm(fm_mat_q(1), fm_mat_q_static_bse)
1318  CALL cp_fm_set_all(fm_mat_q_static_bse, 0.0_dp)
1319 
1320  END IF
1321 
1322  END IF
1323 
1324  IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_create(rpa_grad, fm_mat_q(1), &
1325  fm_mat_s, homo, virtual, mp2_env, eigenval(:, 1, :), &
1326  unit_nr, do_ri_sos_laplace_mp2)
1327 
1328  erpa = 0.0_dp
1329  IF (mp2_env%ri_rpa%do_ri_axk) e_axk = 0.0_dp
1330  first_cycle = .true.
1331  omega_old = 0.0_dp
1332  CALL dgemm_counter_init(dgemm_counter, unit_nr, mp2_env%ri_rpa%print_dgemm_info)
1333 
1334  DO count_ev_sc_gw = 1, iter_evgw
1335 
1336  dbcsr_time = 0.0_dp
1337  dbcsr_nflop = 0
1338 
1339  IF (do_ic_model) cycle
1340 
1341  ! reset some values, important when doing eigenvalue self-consistent GW
1342  IF (my_do_gw) THEN
1343  erpa = 0.0_dp
1344  vec_sigma_c_gw = z_zero
1345  first_cycle = .true.
1346  END IF
1347 
1348  ! calculate Q_PQ(it)
1349  IF (do_im_time) THEN
1350 
1351  IF (.NOT. do_kpoints_cubic_rpa) THEN
1352  DO ispin = 1, nspins
1353  e_fermi(ispin) = (eigenval(homo(ispin), 1, ispin) + eigenval(homo(ispin) + 1, 1, ispin))*0.5_dp
1354  END DO
1355  END IF
1356 
1357  tau = 0.0_dp
1358  tau_old = 0.0_dp
1359 
1360  IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(/T3,A,T66,i15)") &
1361  "MEMORY_INFO| Memory cut:", cut_memory
1362  IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,ES15.2)") &
1363  "SPARSITY_INFO| Eps filter for M virt/occ tensors:", eps_filter
1364  IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,ES15.2)") &
1365  "SPARSITY_INFO| Eps filter for P matrix:", eps_filter_im_time
1366  IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,i15)") &
1367  "SPARSITY_INFO| Minimum tensor block size:", min_bsize
1368 
1369  ! for evGW, we have to ensure that mat_P_omega is zero
1370  CALL zero_mat_p_omega(mat_p_omega(:, :, 1))
1371 
1372  ! compute the matrix Q(it) and Fourier transform it directly to mat_P_omega(iw)
1373  CALL compute_mat_p_omega(mat_p_omega(:, :, 1), fm_scaled_dm_occ_tau, &
1374  fm_scaled_dm_virt_tau, fm_mo_coeff_occ(1), fm_mo_coeff_virt(1), &
1375  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1376  mat_p_global, matrix_s, 1, &
1377  t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, &
1378  starts_array_mc, ends_array_mc, &
1379  starts_array_mc_block, ends_array_mc_block, &
1380  weights_cos_tf_t_to_w, tj, tau_tj, e_fermi(1), eps_filter, alpha, &
1381  eps_filter_im_time, eigenval(:, 1, 1), nmo, &
1382  num_integ_points, cut_memory, &
1383  unit_nr, mp2_env, para_env, &
1384  qs_env, do_kpoints_from_gamma, &
1385  index_to_cell_3c, cell_to_index_3c, &
1386  has_mat_p_blocks, do_ri_sos_laplace_mp2, &
1387  dbcsr_time, dbcsr_nflop)
1388 
1389  ! the same for open shell, use fm_mo_coeff_occ_beta and fm_mo_coeff_virt_beta
1390  IF (my_open_shell) THEN
1391  CALL zero_mat_p_omega(mat_p_omega(:, :, 2))
1392  CALL compute_mat_p_omega(mat_p_omega(:, :, 2), fm_scaled_dm_occ_tau, &
1393  fm_scaled_dm_virt_tau, fm_mo_coeff_occ(2), &
1394  fm_mo_coeff_virt(2), &
1395  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1396  mat_p_global, matrix_s, 2, &
1397  t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, &
1398  starts_array_mc, ends_array_mc, &
1399  starts_array_mc_block, ends_array_mc_block, &
1400  weights_cos_tf_t_to_w, tj, tau_tj, e_fermi(2), eps_filter, alpha, &
1401  eps_filter_im_time, eigenval(:, 1, 2), nmo, &
1402  num_integ_points, cut_memory, &
1403  unit_nr, mp2_env, para_env, &
1404  qs_env, do_kpoints_from_gamma, &
1405  index_to_cell_3c, cell_to_index_3c, &
1406  has_mat_p_blocks, do_ri_sos_laplace_mp2, &
1407  dbcsr_time, dbcsr_nflop)
1408 
1409  !For RPA, we sum up the P matrices. If no force needed, can clean-up the beta spin one
1410  IF (.NOT. do_ri_sos_laplace_mp2) THEN
1411  DO j = 1, SIZE(mat_p_omega, 2)
1412  DO i = 1, SIZE(mat_p_omega, 1)
1413  CALL dbcsr_add(mat_p_omega(i, j, 1)%matrix, mat_p_omega(i, j, 2)%matrix, 1.0_dp, 1.0_dp)
1414  IF (.NOT. calc_forces) CALL dbcsr_clear(mat_p_omega(i, j, 2)%matrix)
1415  END DO
1416  END DO
1417  END IF
1418  END IF ! my_open_shell
1419 
1420  END IF ! do im time
1421 
1422  DO jquad = 1, num_integ_points
1423 
1424  IF (modulo(jquad, num_integ_group) /= color_rpa_group) cycle
1425 
1426  CALL timeset(routinen//"_RPA_matrix_operations", handle3)
1427 
1428  IF (do_ri_sos_laplace_mp2) THEN
1429  omega = tau_tj(jquad)
1430  ELSE
1431  IF (do_minimax_quad) THEN
1432  omega = tj(jquad)
1433  ELSE
1434  omega = a_scaling/tan(tj(jquad))
1435  END IF
1436  END IF ! do_ri_sos_laplace_mp2
1437 
1438  IF (do_im_time) THEN
1439  ! in case we do imag time, we already calculated fm_mat_Q by a Fourier transform from im. time
1440 
1441  IF (.NOT. (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma)) THEN
1442 
1443  DO ispin = 1, SIZE(mat_p_omega, 3)
1444  CALL contract_p_omega_with_mat_l(mat_p_omega(jquad, 1, ispin)%matrix, mat_l%matrix, mat_work, &
1445  eps_filter_im_time, fm_mat_work, dimen_ri, dimen_ri_red, &
1446  fm_mat_minv_l_kpoints(1, 1), fm_mat_q(ispin))
1447  END DO
1448  END IF
1449 
1450  ELSE
1451 
1452  IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3, A, 1X, I3, 1X, A, 1X, I3)") &
1453  "INTEG_INFO| Started with Integration point", jquad, "of", num_integ_points
1454 
1455  IF (first_cycle .AND. count_ev_sc_gw > 1) THEN
1456  IF (iter_sc_gw0 == 1) THEN
1457  DO ispin = 1, nspins
1458  CALL remove_scaling_factor_rpa(fm_mat_s(ispin), virtual(ispin), &
1459  eigenval_last(:, 1, ispin), homo(ispin), omega_old)
1460  END DO
1461  ELSE
1462  DO ispin = 1, nspins
1463  CALL remove_scaling_factor_rpa(fm_mat_s(ispin), virtual(ispin), &
1464  eigenval_scf(:, 1, ispin), homo(ispin), omega_old)
1465  END DO
1466  END IF
1467  END IF
1468 
1469  CALL calc_mat_q(fm_mat_s(1), do_ri_sos_laplace_mp2, first_cycle, iter_sc_gw0, virtual(1), &
1470  eigenval(:, 1, 1), eigenval_scf(:, 1, 1), &
1471  homo(1), omega, omega_old, jquad, mm_style, dimen_ri_red, dimen_ia(1), alpha, fm_mat_q(1), &
1472  fm_mat_q_gemm(1), do_bse, fm_mat_q_static_bse_gemm, dgemm_counter, &
1473  num_integ_points)
1474 
1475  IF (my_open_shell) THEN
1476  CALL calc_mat_q(fm_mat_s(2), do_ri_sos_laplace_mp2, first_cycle, iter_sc_gw0, virtual(2), &
1477  eigenval(:, 1, 2), eigenval_scf(:, 1, 2), &
1478  homo(2), omega, omega_old, jquad, mm_style, &
1479  dimen_ri_red, dimen_ia(2), alpha, fm_mat_q(2), fm_mat_q_gemm(2), do_bse, &
1480  fm_mat_q_static_bse_gemm, dgemm_counter, &
1481  num_integ_points)
1482 
1483  ! For SOS-MP2 we need both matrices separately
1484  IF (.NOT. do_ri_sos_laplace_mp2) THEN
1485  CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_q(1), beta=1.0_dp, matrix_b=fm_mat_q(2))
1486  END IF
1487 
1488  END IF ! open shell
1489 
1490  END IF ! im time
1491 
1492  ! Calculate AXK energy correction
1493  IF (mp2_env%ri_rpa%do_ri_axk) THEN
1494  CALL compute_axk_ener(qs_env, fm_mat_q(1), fm_mat_q_gemm(1), dimen_ri_red, dimen_ia(1), para_env_sub, &
1495  para_env_rpa, eigenval(:, 1, 1), fm_mat_s(1), homo(1), virtual(1), omega, &
1496  mp2_env, mat_munu, unit_nr, e_axk_corr)
1497 
1498  ! Evaluate the final AXK energy correction
1499  e_axk = e_axk + e_axk_corr*wj(jquad)
1500  END IF ! do_ri_axk
1501 
1502  IF (do_ri_sos_laplace_mp2) THEN
1503  CALL sos_mp2_postprocessing(fm_mat_q, erpa, tau_wj(jquad))
1504 
1505  IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
1506  fm_mat_q, fm_mat_q_gemm, dgemm_counter, fm_mat_s, omega, homo, virtual, &
1507  eigenval(:, 1, :), tau_wj(jquad), unit_nr)
1508  ELSE
1509  IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_copy_q(fm_mat_q(1), rpa_grad)
1510 
1511  CALL q_trace_and_add_unit_matrix(dimen_ri_red, trace_qomega, fm_mat_q(1), para_env_rpa)
1512 
1513  IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
1514  CALL invert_eps_compute_w_and_erpa_kp(dimen_ri, num_integ_points, jquad, nkp, count_ev_sc_gw, para_env, &
1515  erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, &
1516  wkp_w, do_gw_im_time, do_ri_sigma_x, do_kpoints_from_gamma, &
1517  cfm_mat_q, ikp_local, &
1518  mat_p_omega(:, :, 1), mat_p_omega_kp, qs_env, eps_filter_im_time, unit_nr, &
1519  kpoints, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, &
1520  fm_mat_w, fm_mat_ri_global_work, mat_minvvminv, &
1521  fm_matrix_minv, fm_matrix_minv_vtrunc_minv)
1522  ELSE
1523  CALL compute_erpa_by_freq_int(dimen_ri_red, trace_qomega, fm_mat_q(1), para_env_rpa, erpa, wj(jquad))
1524  END IF
1525 
1526  IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
1527  fm_mat_q, fm_mat_q_gemm, dgemm_counter, fm_mat_s, omega, homo, virtual, &
1528  eigenval(:, 1, :), wj(jquad), unit_nr)
1529  END IF ! do_ri_sos_laplace_mp2
1530 
1531  ! save omega and reset the first_cycle flag
1532  first_cycle = .false.
1533  omega_old = omega
1534 
1535  CALL timestop(handle3)
1536 
1537  IF (my_do_gw) THEN
1538 
1539  CALL get_fermi_level_offset(fermi_level_offset, fermi_level_offset_input, eigenval(:, 1, :), homo)
1540 
1541  ! do_im_time = TRUE means low-scaling calculation
1542  IF (do_im_time) THEN
1543  ! only for molecules
1544  IF (.NOT. (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma)) THEN
1545  CALL compute_w_cubic_gw(fm_mat_w, fm_mat_q(1), fm_mat_work, dimen_ri, fm_mat_minv_l_kpoints, num_integ_points, &
1546  tj, tau_tj, weights_cos_tf_w_to_t, jquad, omega)
1547  END IF
1548  ELSE
1549  CALL compute_gw_self_energy(vec_sigma_c_gw, dimen_nm_gw, dimen_ri_red, gw_corr_lev_occ, &
1550  gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, num_integ_points, &
1551  do_bse, do_im_time, do_periodic, first_cycle_periodic_correction, &
1552  fermi_level_offset, &
1553  omega, eigenval(:, 1, :), delta_corr, vec_omega_fit_gw, vec_w_gw, wj, &
1554  fm_mat_q(1), fm_mat_q_static_bse, fm_mat_r_gw, fm_mat_s_gw, &
1555  fm_mat_s_gw_work, mo_coeff(1), para_env, &
1556  para_env_rpa, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, &
1557  kpoints, qs_env, mp2_env)
1558  END IF
1559  END IF
1560 
1561  IF (unit_nr > 0) CALL m_flush(unit_nr)
1562  CALL para_env%sync() ! sync to see output
1563 
1564  END DO ! jquad
1565 
1566  CALL para_env%sum(erpa)
1567 
1568  IF (.NOT. do_ri_sos_laplace_mp2) THEN
1569  erpa = erpa/(pi*2.0_dp)
1570  IF (do_minimax_quad) erpa = erpa/2.0_dp
1571  END IF
1572 
1573  IF (mp2_env%ri_rpa%do_ri_axk) THEN
1574  CALL para_env%sum(e_axk)
1575  e_axk = e_axk/(pi*2.0_dp)
1576  IF (do_minimax_quad) e_axk = e_axk/2.0_dp
1577  mp2_env%ri_rpa%ener_axk = e_axk
1578  END IF
1579 
1580  IF (calc_forces .AND. do_ri_sos_laplace_mp2 .AND. do_im_time) THEN
1581  IF (my_open_shell) THEN
1582  pspin = 1
1583  qspin = 2
1584  CALL calc_laplace_loop_forces(force_data, mat_p_omega(:, 1, :), t_3c_m, t_3c_o(1, 1), &
1585  t_3c_o_compressed(1, 1, :), t_3c_o_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1586  fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1587  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1588  starts_array_mc, ends_array_mc, starts_array_mc_block, &
1589  ends_array_mc_block, num_integ_points, nmo, eigenval(:, 1, :), &
1590  tau_tj, tau_wj, cut_memory, pspin, qspin, my_open_shell, &
1591  unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1592  pspin = 2
1593  qspin = 1
1594  CALL calc_laplace_loop_forces(force_data, mat_p_omega(:, 1, :), t_3c_m, t_3c_o(1, 1), &
1595  t_3c_o_compressed(1, 1, :), t_3c_o_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1596  fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1597  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1598  starts_array_mc, ends_array_mc, starts_array_mc_block, &
1599  ends_array_mc_block, num_integ_points, nmo, eigenval(:, 1, :), &
1600  tau_tj, tau_wj, cut_memory, pspin, qspin, my_open_shell, &
1601  unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1602 
1603  ELSE
1604  pspin = 1
1605  qspin = 1
1606  CALL calc_laplace_loop_forces(force_data, mat_p_omega(:, 1, :), t_3c_m, t_3c_o(1, 1), &
1607  t_3c_o_compressed(1, 1, :), t_3c_o_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1608  fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1609  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1610  starts_array_mc, ends_array_mc, starts_array_mc_block, &
1611  ends_array_mc_block, num_integ_points, nmo, eigenval(:, 1, :), &
1612  tau_tj, tau_wj, cut_memory, pspin, qspin, my_open_shell, &
1613  unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1614  END IF
1615  CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1616  END IF !laplace SOS-MP2
1617 
1618  IF (calc_forces .AND. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) THEN
1619  DO ispin = 1, nspins
1620  CALL calc_rpa_loop_forces(force_data, mat_p_omega(:, 1, :), t_3c_m, t_3c_o(1, 1), &
1621  t_3c_o_compressed(1, 1, :), t_3c_o_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1622  fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1623  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1624  starts_array_mc, ends_array_mc, starts_array_mc_block, &
1625  ends_array_mc_block, num_integ_points, nmo, eigenval(:, 1, :), &
1626  e_fermi(ispin), weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, tj, &
1627  wj, tau_tj, cut_memory, ispin, my_open_shell, unit_nr, dbcsr_time, &
1628  dbcsr_nflop, mp2_env, qs_env)
1629  END DO
1630  CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1631  END IF
1632 
1633  IF (do_im_time) THEN
1634 
1635  my_flop_rate = real(dbcsr_nflop, dp)/(1.0e09_dp*dbcsr_time)
1636  IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(/T3,A,T73,ES8.2)") &
1637  "PERFORMANCE| DBCSR total number of flops:", real(dbcsr_nflop*para_env%num_pe, dp)
1638  IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.2)") &
1639  "PERFORMANCE| DBCSR total execution time:", dbcsr_time
1640  IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.2)") &
1641  "PERFORMANCE| DBCSR flop rate (Gflops / MPI rank):", my_flop_rate
1642 
1643  ELSE
1644 
1645  CALL dgemm_counter_write(dgemm_counter, para_env)
1646 
1647  END IF
1648 
1649  ! GW: for low-scaling calculation: Compute self-energy Sigma(i*tau), Sigma(i*omega)
1650  ! for low-scaling and ordinary-scaling: analytic continuation from Sigma(iw) -> Sigma(w)
1651  ! and correction of quasiparticle energies e_n^GW
1652  IF (my_do_gw) THEN
1653 
1654  CALL compute_qp_energies(vec_sigma_c_gw, count_ev_sc_gw, gw_corr_lev_occ, &
1655  gw_corr_lev_tot, gw_corr_lev_virt, homo, &
1656  nmo, num_fit_points, num_integ_points, &
1657  unit_nr, do_apply_ic_corr_to_gw, do_im_time, &
1658  do_periodic, do_ri_sigma_x, first_cycle_periodic_correction, &
1659  e_fermi, eps_filter, fermi_level_offset, &
1660  delta_corr, eigenval, &
1661  eigenval_last, eigenval_scf, iter_sc_gw0, exit_ev_gw, tau_tj, tj, &
1662  vec_omega_fit_gw, vec_sigma_x_gw, mp2_env%ri_g0w0%ic_corr_list, &
1663  weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, &
1664  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_mo_coeff_occ, &
1665  fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
1666  mo_coeff(1), fm_mat_w, para_env, para_env_rpa, mat_dm, mat_minvvminv, &
1667  t_3c_o, t_3c_m, t_3c_overl_int_ao_mo, t_3c_o_compressed, t_3c_o_mo_compressed, &
1668  t_3c_o_ind, t_3c_o_mo_ind, &
1669  t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, &
1670  matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, mat_w, matrix_s, &
1671  kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_rpa, &
1672  starts_array_mc, ends_array_mc)
1673 
1674  ! if HOMO-LUMO gap differs by less than mp2_env%ri_g0w0%eps_ev_sc_iter, exit ev sc GW loop
1675  IF (exit_ev_gw) EXIT
1676 
1677  END IF ! my_do_gw if
1678 
1679  END DO ! evGW loop
1680 
1681  IF (do_ic_model) THEN
1682 
1683  IF (my_open_shell) THEN
1684 
1685  CALL calculate_ic_correction(eigenval(:, 1, 1), mat_minvvminv%matrix, &
1686  t_3c_overl_nnp_ic(1), t_3c_overl_nnp_ic_reflected(1), &
1687  gw_corr_lev_tot, &
1688  gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), unit_nr, &
1689  print_ic_values, para_env, do_alpha=.true.)
1690 
1691  CALL calculate_ic_correction(eigenval(:, 1, 2), mat_minvvminv%matrix, &
1692  t_3c_overl_nnp_ic(2), t_3c_overl_nnp_ic_reflected(2), &
1693  gw_corr_lev_tot, &
1694  gw_corr_lev_occ(2), gw_corr_lev_virt(2), homo(2), unit_nr, &
1695  print_ic_values, para_env, do_beta=.true.)
1696 
1697  ELSE
1698 
1699  CALL calculate_ic_correction(eigenval(:, 1, 1), mat_minvvminv%matrix, &
1700  t_3c_overl_nnp_ic(1), t_3c_overl_nnp_ic_reflected(1), &
1701  gw_corr_lev_tot, &
1702  gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), unit_nr, &
1703  print_ic_values, para_env)
1704 
1705  END IF
1706 
1707  END IF
1708 
1709  ! postprocessing after GW for Bethe-Salpeter
1710  IF (do_bse) THEN
1711  ! Create a copy of fm_mat_S for usage in BSE
1712  CALL cp_fm_create(fm_mat_s_ia_bse, fm_mat_s(1)%matrix_struct)
1713  CALL cp_fm_to_fm(fm_mat_s(1), fm_mat_s_ia_bse)
1714  ! Remove energy/frequency factor from 3c-Integral for BSE
1715  IF (iter_sc_gw0 == 1) THEN
1716  CALL remove_scaling_factor_rpa(fm_mat_s_ia_bse, virtual(1), &
1717  eigenval_last(:, 1, 1), homo(1), omega)
1718  ELSE
1719  CALL remove_scaling_factor_rpa(fm_mat_s_ia_bse, virtual(1), &
1720  eigenval_scf(:, 1, 1), homo(1), omega)
1721  END IF
1722  ! Main routine for all BSE postprocessing
1723  CALL start_bse_calculation(fm_mat_s_ia_bse, fm_mat_s_ij_bse, fm_mat_s_ab_bse, &
1724  fm_mat_q_static_bse, fm_mat_q_static_bse_gemm, &
1725  eigenval, homo, virtual, dimen_ri, dimen_ri_red, &
1726  gd_array, color_sub, mp2_env, unit_nr)
1727  ! Release BSE-copy of fm_mat_S
1728  CALL cp_fm_release(fm_mat_s_ia_bse)
1729  END IF
1730 
1731  IF (my_do_gw) THEN
1732  CALL deallocate_matrices_gw(fm_mat_s_gw_work, vec_w_gw, vec_sigma_c_gw, vec_omega_fit_gw, &
1733  mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw, &
1734  eigenval_last, eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
1735  kpoints, vec_sigma_x_gw,.NOT. do_im_time)
1736  END IF
1737 
1738  IF (do_im_time) THEN
1739 
1740  CALL dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, &
1741  fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, index_to_cell_3c, &
1742  cell_to_index_3c, do_ic_model, &
1743  do_kpoints_cubic_rpa, do_kpoints_from_gamma, do_ri_sigma_x, &
1744  has_mat_p_blocks, &
1745  wkp_w, cfm_mat_q, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, &
1746  fm_matrix_minv, fm_matrix_minv_vtrunc_minv, fm_mat_ri_global_work, fm_mat_work, &
1747  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_l, &
1748  mat_minvvminv, mat_p_omega, mat_p_omega_kp, &
1749  t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, mat_work, qs_env)
1750 
1751  IF (my_do_gw) THEN
1752  CALL deallocate_matrices_gw_im_time(weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, do_ic_model, &
1753  do_kpoints_cubic_rpa, fm_mat_w, &
1754  t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, &
1755  t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, &
1756  t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, &
1757  mat_w, qs_env)
1758  END IF
1759 
1760  END IF
1761 
1762  IF (.NOT. do_ri_sos_laplace_mp2) THEN
1763  DEALLOCATE (tj)
1764  DEALLOCATE (wj)
1765  DEALLOCATE (trace_qomega)
1766  END IF
1767 
1768  IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
1769  DEALLOCATE (tau_tj)
1770  DEALLOCATE (tau_wj)
1771  END IF
1772 
1773  IF (do_im_time .AND. calc_forces) THEN
1774  CALL im_time_force_release(force_data)
1775  END IF
1776 
1777  IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_finalize(rpa_grad, mp2_env, para_env_sub, para_env, &
1778  qs_env, gd_array, color_sub, do_ri_sos_laplace_mp2, &
1779  homo, virtual)
1780 
1781  CALL timestop(handle)
1782 
1783  END SUBROUTINE rpa_num_int
1784 
1785 END MODULE rpa_main
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
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::wilhelm2016b
integer, save, public wilhelm2016b
Definition: bibliography.F:43
bibliography::delben2015
integer, save, public delben2015
Definition: bibliography.F:43
bibliography::wilhelm2016a
integer, save, public wilhelm2016a
Definition: bibliography.F:43
bibliography::wilhelm2018
integer, save, public wilhelm2018
Definition: bibliography.F:43
bibliography::delben2013
integer, save, public delben2013
Definition: bibliography.F:43
bibliography::ren2013
integer, save, public ren2013
Definition: bibliography.F:43
bibliography::wilhelm2017
integer, save, public wilhelm2017
Definition: bibliography.F:43
bibliography::ren2011
integer, save, public ren2011
Definition: bibliography.F:43
bibliography::bates2013
integer, save, public bates2013
Definition: bibliography.F:43
bse_main
Main routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition: bse_main.F:14
bse_main::start_bse_calculation
subroutine, public start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm, Eigenval, homo, virtual, dimen_RI, dimen_RI_red, gd_array, color_sub, mp2_env, unit_nr)
Main subroutine managing BSE calculations.
Definition: bse_main.F:72
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_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:101
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_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_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_init
elemental subroutine, public dgemm_counter_init(dgemm_counter, unit_nr, print_info)
Initialize a dgemm_counter.
Definition: dgemm_counter_types.F:49
dgemm_counter_types::dgemm_counter_write
subroutine, public dgemm_counter_write(dgemm_counter, para_env)
calculate and print flop rates
Definition: dgemm_counter_types.F:103
group_dist_types
Types to describe group distributions.
Definition: group_dist_types.F:14
group_dist_types::maxsize
elemental integer function, public maxsize(this)
...
Definition: group_dist_types.F:270
hfx_types
Types and set/get functions for HFX.
Definition: hfx_types.F:15
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_gemm
integer, parameter, public wfc_mm_style_gemm
Definition: input_constants.F:1046
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition: kinds.F:54
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
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition: machine.F:332
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
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
minimax_exp
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition: minimax_exp.F:29
minimax_exp::check_exp_minimax_range
subroutine, public check_exp_minimax_range(k, Rc, ierr)
Check that a minimax approximation is available for given input k, Rc. ierr == 0: everything ok ierr ...
Definition: minimax_exp.F:101
mp2_grids
Routines to calculate frequency and time grids (integration points and weights) for correlation metho...
Definition: mp2_grids.F:14
mp2_grids::get_clenshaw_grid
subroutine, public get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, num_integ_group, color_rpa_group, fm_mat_S, my_do_gw, ext_scaling, a_scaling, tj, wj)
...
Definition: mp2_grids.F:314
mp2_grids::get_minimax_grid
subroutine, public get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, do_ri_sos_laplace_mp2, do_print, tau_tj, tau_wj, qs_env, do_gw_im_time, do_kpoints_cubic_RPA, e_fermi, tj, wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, regularization)
...
Definition: mp2_grids.F:75
mp2_laplace
Routines to calculate MP2 energy with laplace approach.
Definition: mp2_laplace.F:13
mp2_laplace::sos_mp2_postprocessing
subroutine, public sos_mp2_postprocessing(fm_mat_Q, Erpa, tau_wjquad)
...
Definition: mp2_laplace.F:81
mp2_ri_grad_util
Routines for calculating RI-MP2 gradients.
Definition: mp2_ri_grad_util.F:13
mp2_ri_grad_util::array2fm
subroutine, public array2fm(mat2D, fm_struct, my_start_row, my_end_row, my_start_col, my_end_col, gd_row, gd_col, group_grid_2_mepos, ngroup_row, ngroup_col, fm_mat, integ_group_size, color_group, do_release_mat)
redistribute local part of array to fm
Definition: mp2_ri_grad_util.F:495
mp2_types
Types needed for MP2 calculations.
Definition: mp2_types.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_axk
Auxiliary routines needed for RPA-AXK given blacs_env to another.
Definition: rpa_axk.F:17
rpa_axk::compute_axk_ener
subroutine, public compute_axk_ener(qs_env, fm_mat_Q, fm_mat_Q_gemm, dimen_RI, dimen_ia, para_env_sub, para_env_RPA, eig, fm_mat_S, homo, virtual, omega, mp2_env, mat_munu, unit_nr, e_axk_corr)
Main driver for RPA-AXK energies.
Definition: rpa_axk.F:98
rpa_grad
Routines to calculate RI-RPA and SOS-MP2 gradients.
Definition: rpa_grad.F:13
rpa_grad::rpa_grad_matrix_operations
subroutine, public rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, fm_mat_Q, fm_mat_Q_gemm, dgemm_counter, fm_mat_S, omega, homo, virtual, Eigenval, weight, unit_nr)
...
Definition: rpa_grad.F:746
rpa_grad::rpa_grad_needed_mem
pure subroutine, public rpa_grad_needed_mem(homo, virtual, dimen_RI, mem_per_rank, mem_per_repl, do_ri_sos_laplace_mp2)
Calculates the necessary minimum memory for the Gradient code ion MiB.
Definition: rpa_grad.F:124
rpa_grad::rpa_grad_create
subroutine, public rpa_grad_create(rpa_grad, fm_mat_Q, fm_mat_S, homo, virtual, mp2_env, Eigenval, unit_nr, do_ri_sos_laplace_mp2)
Creates the arrays of a rpa_grad_type.
Definition: rpa_grad.F:170
rpa_grad::rpa_grad_finalize
subroutine, public rpa_grad_finalize(rpa_grad, mp2_env, para_env_sub, para_env, qs_env, gd_array, color_sub, do_ri_sos_laplace_mp2, homo, virtual)
...
Definition: rpa_grad.F:2042
rpa_grad::rpa_grad_copy_q
subroutine, public rpa_grad_copy_q(fm_mat_Q, rpa_grad)
...
Definition: rpa_grad.F:721
rpa_gw_ic
Routines to calculate image charge corrections.
Definition: rpa_gw_ic.F:13
rpa_gw_ic::calculate_ic_correction
subroutine, public calculate_ic_correction(Eigenval, mat_SinvVSinv, t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, unit_nr, print_ic_values, para_env, do_alpha, do_beta)
...
Definition: rpa_gw_ic.F:58
rpa_gw_kpoints_util
Routines treating GW and RPA calculations with kpoints.
Definition: rpa_gw_kpoints_util.F:13
rpa_gw_kpoints_util::invert_eps_compute_w_and_erpa_kp
subroutine, public invert_eps_compute_w_and_erpa_kp(dimen_RI, num_integ_points, jquad, nkp, count_ev_sc_GW, para_env, Erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, wkp_W, do_gw_im_time, do_ri_Sigma_x, do_kpoints_from_Gamma, cfm_mat_Q, ikp_local, mat_P_omega, mat_P_omega_kp, qs_env, eps_filter_im_time, unit_nr, kpoints, fm_mat_Minv_L_kpoints, fm_matrix_L_kpoints, fm_mat_W, fm_mat_RI_global_work, mat_MinvVMinv, fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv)
...
Definition: rpa_gw_kpoints_util.F:136
rpa_gw_kpoints_util::get_bandstruc_and_k_dependent_mos
subroutine, public get_bandstruc_and_k_dependent_mos(qs_env, Eigenval_kp)
...
Definition: rpa_gw_kpoints_util.F:1480
rpa_gw
Routines for GW, continuous development [Jan Wilhelm].
Definition: rpa_gw.F:14
rpa_gw::get_fermi_level_offset
subroutine, public get_fermi_level_offset(fermi_level_offset, fermi_level_offset_input, Eigenval, homo)
...
Definition: rpa_gw.F:920
rpa_gw::allocate_matrices_gw
subroutine, public allocate_matrices_gw(vec_Sigma_c_gw, color_rpa_group, dimen_nm_gw, gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, num_integ_group, num_integ_points, unit_nr, gw_corr_lev_tot, num_fit_points, omega_max_fit, do_minimax_quad, do_periodic, do_ri_Sigma_x, my_do_gw, first_cycle_periodic_correction, a_scaling, Eigenval, tj, vec_omega_fit_gw, vec_Sigma_x_gw, delta_corr, Eigenval_last, Eigenval_scf, vec_W_gw, fm_mat_S_gw, fm_mat_S_gw_work, para_env, mp2_env, kpoints, nkp, nkp_self_energy, do_kpoints_cubic_RPA, do_kpoints_from_Gamma)
...
Definition: rpa_gw.F:349
rpa_gw::allocate_matrices_gw_im_time
subroutine, public allocate_matrices_gw_im_time(gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, num_integ_points, unit_nr, RI_blk_sizes, do_ic_model, para_env, fm_mat_W, fm_mat_Q, mo_coeff, t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, starts_array_mc, ends_array_mc, t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, matrix_s, mat_W, t_3c_overl_int, t_3c_O_compressed, t_3c_O_ind, qs_env)
...
Definition: rpa_gw.F:203
rpa_gw::compute_qp_energies
subroutine, public compute_qp_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, homo, nmo, num_fit_points, num_integ_points, unit_nr, do_apply_ic_corr_to_gw, do_im_time, do_periodic, do_ri_Sigma_x, first_cycle_periodic_correction, e_fermi, eps_filter, fermi_level_offset, delta_corr, Eigenval, Eigenval_last, Eigenval_scf, iter_sc_GW0, exit_ev_gw, tau_tj, tj, vec_omega_fit_gw, vec_Sigma_x_gw, ic_corr_list, weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, mo_coeff, fm_mat_W, para_env, para_env_RPA, mat_dm, mat_MinvVMinv, t_3c_O, t_3c_M, t_3c_overl_int_ao_mo, t_3c_O_compressed, t_3c_O_mo_compressed, t_3c_O_ind, t_3c_O_mo_ind, t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, mat_W, matrix_s, kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_RPA, starts_array_mc, ends_array_mc)
...
Definition: rpa_gw.F:1236
rpa_gw::compute_w_cubic_gw
subroutine, public compute_w_cubic_gw(fm_mat_W, fm_mat_Q, fm_mat_work, dimen_RI, fm_mat_L, num_integ_points, tj, tau_tj, weights_cos_tf_w_to_t, jquad, omega)
...
Definition: rpa_gw.F:962
rpa_gw::compute_gw_self_energy
subroutine, public compute_gw_self_energy(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw_corr_lev_occ, gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, num_integ_points, do_bse, do_im_time, do_periodic, first_cycle_periodic_correction, fermi_level_offset, omega, Eigenval, delta_corr, vec_omega_fit_gw, vec_W_gw, wj, fm_mat_Q, fm_mat_Q_static_bse, fm_mat_R_gw, fm_mat_S_gw, fm_mat_S_gw_work, mo_coeff, para_env, para_env_RPA, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, kpoints, qs_env, mp2_env)
...
Definition: rpa_gw.F:811
rpa_gw::deallocate_matrices_gw
subroutine, public deallocate_matrices_gw(fm_mat_S_gw_work, vec_W_gw, vec_Sigma_c_gw, vec_omega_fit_gw, vec_Sigma_x_minus_vxc_gw, Eigenval_last, Eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, kpoints, vec_Sigma_x_gw, my_do_gw)
...
Definition: rpa_gw.F:614
rpa_gw::deallocate_matrices_gw_im_time
subroutine, public deallocate_matrices_gw_im_time(weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, do_ic_model, do_kpoints_cubic_RPA, fm_mat_W, t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, mat_W, qs_env)
...
Definition: rpa_gw.F:689
rpa_im_time_force_methods
Routines needed for cubic-scaling RPA and SOS-Laplace-MP2 forces.
Definition: rpa_im_time_force_methods.F:12
rpa_im_time_force_methods::calc_rpa_loop_forces
subroutine, public calc_rpa_loop_forces(force_data, mat_P_omega, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, num_integ_points, nmo, Eigenval, e_fermi, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, tj, wj, tau_tj, cut_memory, ispin, open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
Updates the cubic-scaling RPA contribution to the forces at each quadrature point....
Definition: rpa_im_time_force_methods.F:1218
rpa_im_time_force_methods::calc_laplace_loop_forces
subroutine, public calc_laplace_loop_forces(force_data, mat_P_omega, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, num_integ_points, nmo, Eigenval, tau_tj, tau_wj, cut_memory, Pspin, Qspin, open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
Updates the cubic-scaling SOS-Laplace-MP2 contribution to the forces at each quadrature point.
Definition: rpa_im_time_force_methods.F:667
rpa_im_time_force_methods::calc_post_loop_forces
subroutine, public calc_post_loop_forces(force_data, unit_nr, qs_env)
All the forces that can be calculated after the loop on the Laplace quaradture, using terms collected...
Definition: rpa_im_time_force_methods.F:2211
rpa_im_time_force_methods::init_im_time_forces
subroutine, public init_im_time_forces(force_data, fm_matrix_PQ, t_3c_M, unit_nr, mp2_env, qs_env)
Initializes and pre-calculates all needed tensors for the forces.
Definition: rpa_im_time_force_methods.F:188
rpa_im_time_force_methods::keep_initial_quad
subroutine, public keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, do_laplace, do_im_time, num_integ_points, unit_nr, qs_env)
Overwrites the "optimal" Laplace quadrature with that of the first step.
Definition: rpa_im_time_force_methods.F:3929
rpa_im_time_force_types
Types needed for cubic-scaling RPA and SOS-Laplace-MP2 forces.
Definition: rpa_im_time_force_types.F:12
rpa_im_time_force_types::im_time_force_release
subroutine, public im_time_force_release(force_data)
Cleans everything up.
Definition: rpa_im_time_force_types.F:91
rpa_im_time
Routines for low-scaling RPA/GW with imaginary time.
Definition: rpa_im_time.F:13
rpa_im_time::compute_mat_p_omega
subroutine, public compute_mat_p_omega(mat_P_omega, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_P_global, matrix_s, ispin, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, weights_cos_tf_t_to_w, tj, tau_tj, e_fermi, eps_filter, alpha, eps_filter_im_time, Eigenval, nmo, num_integ_points, cut_memory, unit_nr, mp2_env, para_env, qs_env, do_kpoints_from_Gamma, index_to_cell_3c, cell_to_index_3c, has_mat_P_blocks, do_ri_sos_laplace_mp2, dbcsr_time, dbcsr_nflop)
...
Definition: rpa_im_time.F:139
rpa_im_time::zero_mat_p_omega
subroutine, public zero_mat_p_omega(mat_P_omega)
...
Definition: rpa_im_time.F:616
rpa_main
Routines to calculate RI-RPA energy.
Definition: rpa_main.F:17
rpa_main::rpa_ri_compute_en
subroutine, public rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, para_env, para_env_sub, color_sub, gd_array, gd_B_virtual, gd_B_all, gd_B_occ_bse, gd_B_virt_bse, mo_coeff, fm_matrix_PQ, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, kpoints, Eigenval, nmo, homo, dimen_RI, dimen_RI_red, gw_corr_lev_occ, gw_corr_lev_virt, unit_nr, do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_bse, matrix_s, mat_munu, mat_P_global, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, calc_forces)
...
Definition: rpa_main.F:181
rpa_util
Utility functions for RPA calculations.
Definition: rpa_util.F:13
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
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333