(git:34ef472)
mp2_ri_2c.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 Framework for 2c-integrals for RI
10 !> \par History
11 !> 06.2012 created [Mauro Del Ben]
12 !> 03.2019 separated from mp2_ri_gpw [Frederick Stein]
13 ! **************************************************************************************************
14 MODULE mp2_ri_2c
15  USE atomic_kind_types, ONLY: atomic_kind_type,&
16  get_atomic_kind_set
17  USE basis_set_types, ONLY: gto_basis_set_p_type,&
18  gto_basis_set_type
19  USE cell_types, ONLY: cell_type,&
20  get_cell,&
21  plane_distance
22  USE constants_operator, ONLY: operator_coulomb
23  USE cp_blacs_env, ONLY: cp_blacs_env_create,&
24  cp_blacs_env_release,&
25  cp_blacs_env_type
26  USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm
27  USE cp_cfm_types, ONLY: cp_cfm_create,&
28  cp_cfm_release,&
29  cp_cfm_to_fm,&
30  cp_cfm_type
31  USE cp_control_types, ONLY: dft_control_type
32  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
33  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
34  cp_dbcsr_dist2d_to_dist,&
35  dbcsr_allocate_matrix_set,&
36  dbcsr_deallocate_matrix_set
37  USE cp_eri_mme_interface, ONLY: cp_eri_mme_param
38  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
39  cp_fm_triangular_invert
40  USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
41  USE cp_fm_diag, ONLY: choose_eigv_solver,&
42  cp_fm_power,&
43  cp_fm_syevx
44  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
45  cp_fm_struct_release,&
46  cp_fm_struct_type
47  USE cp_fm_types, ONLY: cp_fm_copy_general,&
48  cp_fm_create,&
49  cp_fm_get_info,&
50  cp_fm_indxg2p,&
51  cp_fm_release,&
52  cp_fm_set_all,&
53  cp_fm_to_fm,&
54  cp_fm_type
55  USE dbcsr_api, ONLY: &
56  dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
57  dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_p_type, dbcsr_release, &
58  dbcsr_reserve_all_blocks, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
59  dbcsr_type_no_symmetry, dbcsr_type_symmetric
60  USE distribution_2d_types, ONLY: distribution_2d_type
61  USE group_dist_types, ONLY: create_group_dist,&
62  get_group_dist,&
63  group_dist_d1_type,&
64  release_group_dist
65  USE input_constants, ONLY: do_eri_gpw,&
66  do_eri_mme,&
67  do_eri_os,&
68  do_potential_coulomb,&
69  do_potential_id,&
70  do_potential_long,&
71  do_potential_truncated
72  USE kinds, ONLY: dp
73  USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
74  USE kpoint_methods, ONLY: rskp_transform
75  USE kpoint_types, ONLY: get_kpoint_info,&
76  kpoint_type
77  USE libint_2c_3c, ONLY: compare_potential_types,&
78  libint_potential_type
79  USE machine, ONLY: m_flush
80  USE message_passing, ONLY: mp_comm_type,&
81  mp_para_env_release,&
82  mp_para_env_type,&
83  mp_request_type
84  USE mp2_eri, ONLY: mp2_eri_2c_integrate
85  USE mp2_eri_gpw, ONLY: mp2_eri_2c_integrate_gpw
86  USE mp2_types, ONLY: integ_mat_buffer_type
87  USE parallel_gemm_api, ONLY: parallel_gemm
88  USE particle_methods, ONLY: get_particle_set
89  USE particle_types, ONLY: particle_type
90  USE qs_environment_types, ONLY: get_qs_env,&
91  qs_environment_type
92  USE qs_integral_utils, ONLY: basis_set_list_setup
93  USE qs_interactions, ONLY: init_interaction_radii
94  USE qs_kind_types, ONLY: get_qs_kind,&
95  qs_kind_type
96  USE qs_ks_types, ONLY: get_ks_env,&
97  set_ks_env
98  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
99  release_neighbor_list_sets
100  USE qs_tensors, ONLY: build_2c_integrals,&
101  build_2c_neighbor_lists
102  USE rpa_communication, ONLY: communicate_buffer
103  USE rpa_gw_kpoints_util, ONLY: cp_cfm_power
104 
105 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
106 #include "./base/base_uses.f90"
107 
108  IMPLICIT NONE
109 
110  PRIVATE
111 
112  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_2c'
113 
114  PUBLIC :: get_2c_integrals, setup_trunc_coulomb_pot_for_exchange_self_energy, ri_2c_integral_mat, &
115  inversion_of_m_and_mult_with_chol_dec_of_v
116 
117 CONTAINS
118 
119 ! **************************************************************************************************
120 !> \brief ...
121 !> \param qs_env ...
122 !> \param eri_method ...
123 !> \param eri_param ...
124 !> \param para_env ...
125 !> \param para_env_sub ...
126 !> \param mp2_memory ...
127 !> \param my_Lrows ...
128 !> \param my_Vrows ...
129 !> \param fm_matrix_PQ ...
130 !> \param ngroup ...
131 !> \param color_sub ...
132 !> \param dimen_RI ...
133 !> \param dimen_RI_red reduced RI dimension, needed if we perform SVD instead of Cholesky
134 !> \param kpoints ...
135 !> \param my_group_L_size ...
136 !> \param my_group_L_start ...
137 !> \param my_group_L_end ...
138 !> \param gd_array ...
139 !> \param calc_PQ_cond_num ...
140 !> \param do_svd ...
141 !> \param potential ...
142 !> \param ri_metric ...
143 !> \param fm_matrix_L_kpoints ...
144 !> \param fm_matrix_Minv_L_kpoints ...
145 !> \param fm_matrix_Minv ...
146 !> \param fm_matrix_Minv_Vtrunc_Minv ...
147 !> \param do_im_time ...
148 !> \param do_kpoints ...
149 !> \param mp2_eps_pgf_orb_S ...
150 !> \param qs_kind_set ...
151 !> \param sab_orb_sub ...
152 !> \param calc_forces ...
153 !> \param unit_nr ...
154 ! **************************************************************************************************
155  SUBROUTINE get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, mp2_memory, &
156  my_Lrows, my_Vrows, fm_matrix_PQ, ngroup, color_sub, dimen_RI, dimen_RI_red, &
157  kpoints, my_group_L_size, my_group_L_start, my_group_L_end, &
158  gd_array, calc_PQ_cond_num, do_svd, potential, ri_metric, &
159  fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
160  fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
161  do_im_time, do_kpoints, mp2_eps_pgf_orb_S, qs_kind_set, sab_orb_sub, calc_forces, unit_nr)
162 
163  TYPE(qs_environment_type), POINTER :: qs_env
164  INTEGER, INTENT(IN) :: eri_method
165  TYPE(cp_eri_mme_param), POINTER :: eri_param
166  TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
167  REAL(kind=dp), INTENT(IN) :: mp2_memory
168  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
169  INTENT(OUT) :: my_lrows, my_vrows
170  TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_pq
171  INTEGER, INTENT(IN) :: ngroup, color_sub
172  INTEGER, INTENT(OUT) :: dimen_ri, dimen_ri_red
173  TYPE(kpoint_type), POINTER :: kpoints
174  INTEGER, INTENT(OUT) :: my_group_l_size, my_group_l_start, &
175  my_group_l_end
176  TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array
177  LOGICAL, INTENT(IN) :: calc_pq_cond_num, do_svd
178  TYPE(libint_potential_type) :: potential, ri_metric
179  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_l_kpoints, &
180  fm_matrix_minv_l_kpoints, &
181  fm_matrix_minv, &
182  fm_matrix_minv_vtrunc_minv
183  LOGICAL, INTENT(IN) :: do_im_time, do_kpoints
184  REAL(kind=dp), INTENT(IN) :: mp2_eps_pgf_orb_s
185  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
186  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
187  POINTER :: sab_orb_sub
188  LOGICAL, INTENT(IN) :: calc_forces
189  INTEGER, INTENT(IN) :: unit_nr
190 
191  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_2c_integrals'
192 
193  INTEGER :: handle, num_small_eigen
194  REAL(kind=dp) :: cond_num, eps_pgf_orb_old
195  TYPE(cp_fm_type) :: fm_matrix_l_work, fm_matrix_m_inv_work, &
196  fm_matrix_v
197  TYPE(dft_control_type), POINTER :: dft_control
198  TYPE(libint_potential_type) :: trunc_coulomb
199  TYPE(mp_para_env_type), POINTER :: para_env_l
200 
201  CALL timeset(routinen, handle)
202 
203  ! calculate V and store it in fm_matrix_L_work
204  CALL compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_l, mp2_memory, &
205  fm_matrix_l_work, ngroup, color_sub, dimen_ri, &
206  my_group_l_size, my_group_l_start, my_group_l_end, &
207  gd_array, calc_pq_cond_num, cond_num, &
208  num_small_eigen, potential, sab_orb_sub, do_im_time=do_im_time)
209 
210  IF (do_im_time .AND. calc_forces) THEN
211  !need a copy of the (P|Q) integrals
212  CALL cp_fm_create(fm_matrix_pq, fm_matrix_l_work%matrix_struct)
213  CALL cp_fm_to_fm(fm_matrix_l_work, fm_matrix_pq)
214  END IF
215 
216  dimen_ri_red = dimen_ri
217 
218  IF (compare_potential_types(ri_metric, potential) .AND. .NOT. do_im_time) THEN
219  CALL decomp_mat_l(fm_matrix_l_work, do_svd, num_small_eigen, cond_num, .true., gd_array, ngroup, &
220  dimen_ri, dimen_ri_red, para_env)
221  ELSE
222 
223  ! RI-metric matrix (depending on RI metric: Coulomb, overlap or truncated Coulomb)
224  IF (do_im_time) THEN
225  CALL get_qs_env(qs_env, dft_control=dft_control)
226 
227  ! re-init the radii to be able to generate pair lists with appropriate screening for overlap matrix
228  eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
229  dft_control%qs_control%eps_pgf_orb = mp2_eps_pgf_orb_s
230  CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
231 
232  CALL ri_2c_integral_mat(qs_env, fm_matrix_minv_l_kpoints, fm_matrix_l_work, dimen_ri, ri_metric, &
233  do_kpoints, kpoints, put_mat_ks_env=.true., &
234  regularization_ri=qs_env%mp2_env%ri_rpa_im_time%regularization_RI)
235 
236  ! re-init the radii to previous values
237  dft_control%qs_control%eps_pgf_orb = eps_pgf_orb_old
238  CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
239 
240  ! GPW overlap matrix
241  ELSE
242 
243  CALL mp_para_env_release(para_env_l)
244  CALL release_group_dist(gd_array)
245 
246  ALLOCATE (fm_matrix_minv_l_kpoints(1, 1))
247 
248  ! Calculate matrix of RI operator (for overlap metric: S), store it in fm_matrix_Minv_L_kpoints
249  CALL compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_l, mp2_memory, &
250  fm_matrix_minv_l_kpoints(1, 1), ngroup, color_sub, dimen_ri, &
251  my_group_l_size, my_group_l_start, my_group_l_end, &
252  gd_array, calc_pq_cond_num, cond_num, &
253  num_small_eigen, ri_metric, sab_orb_sub, &
254  fm_matrix_l_extern=fm_matrix_l_work)
255 
256  END IF
257 
258  IF (do_kpoints) THEN
259 
260  ! k-dependent 1/r Coulomb matrix V_PQ(k)
261  CALL compute_v_by_lattice_sum(qs_env, fm_matrix_l_kpoints, fm_matrix_minv_l_kpoints, kpoints)
262 
263  CALL inversion_of_m_and_mult_with_chol_dec_of_v(fm_matrix_minv_l_kpoints, fm_matrix_l_kpoints, dimen_ri, &
264  kpoints, qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S)
265 
266  CALL setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env, trunc_coulomb)
267 
268  ! Gamma-only truncated Coulomb matrix V^tr with cutoff radius = half the unit cell size; for exchange self-energy
269  CALL ri_2c_integral_mat(qs_env, fm_matrix_minv_vtrunc_minv, fm_matrix_l_work, dimen_ri, trunc_coulomb, &
270  do_kpoints=.false., kpoints=kpoints, put_mat_ks_env=.false., regularization_ri=0.0_dp)
271 
272  ! Gamma-only RI-metric matrix; for computing Gamma-only/MIC self-energy
273  CALL ri_2c_integral_mat(qs_env, fm_matrix_minv, fm_matrix_l_work, dimen_ri, ri_metric, &
274  do_kpoints=.false., kpoints=kpoints, put_mat_ks_env=.false., regularization_ri=0.0_dp)
275 
276  CALL gamma_only_inversion_of_m_and_mult_with_chol_dec_of_vtrunc(fm_matrix_minv_vtrunc_minv, &
277  fm_matrix_minv, qs_env)
278  ELSE
279  IF (calc_forces .AND. (.NOT. do_im_time)) THEN
280  ! For the gradients, we make a copy of the inverse root of L
281  CALL cp_fm_create(fm_matrix_v, fm_matrix_l_work%matrix_struct)
282  CALL cp_fm_to_fm(fm_matrix_l_work, fm_matrix_v)
283 
284  CALL decomp_mat_l(fm_matrix_v, do_svd, num_small_eigen, cond_num, .true., gd_array, ngroup, &
285  dimen_ri, dimen_ri_red, para_env)
286  END IF
287 
288  CALL decomp_mat_l(fm_matrix_l_work, do_svd, num_small_eigen, cond_num, .false., gd_array, ngroup, &
289  dimen_ri, dimen_ri_red, para_env)
290 
291  CALL decomp_mat_l(fm_matrix_minv_l_kpoints(1, 1), .false., num_small_eigen, cond_num, .true., &
292  gd_array, ngroup, dimen_ri, dimen_ri_red, para_env)
293 
294  CALL cp_fm_create(fm_matrix_m_inv_work, fm_matrix_minv_l_kpoints(1, 1)%matrix_struct)
295  CALL cp_fm_set_all(fm_matrix_m_inv_work, 0.0_dp)
296 
297  CALL parallel_gemm('N', 'T', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_matrix_minv_l_kpoints(1, 1), &
298  fm_matrix_minv_l_kpoints(1, 1), 0.0_dp, fm_matrix_m_inv_work)
299 
300  IF (do_svd) THEN
301  ! We have to reset the size of fm_matrix_Minv_L_kpoints
302  CALL reset_size_matrix(fm_matrix_minv_l_kpoints(1, 1), dimen_ri_red, fm_matrix_l_work%matrix_struct)
303 
304  ! L (=fm_matrix_Minv_L_kpoints) = S^(-1)*chol_dec(V)
305  CALL parallel_gemm('T', 'N', dimen_ri, dimen_ri_red, dimen_ri, 1.0_dp, fm_matrix_m_inv_work, &
306  fm_matrix_l_work, 0.0_dp, fm_matrix_minv_l_kpoints(1, 1))
307  ELSE
308  ! L (=fm_matrix_Minv_L_kpoints) = S^(-1)*chol_dec(V)
309  CALL parallel_gemm('T', 'T', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_matrix_m_inv_work, &
310  fm_matrix_l_work, 0.0_dp, fm_matrix_minv_l_kpoints(1, 1))
311  END IF
312 
313  ! clean the S_inv matrix
314  CALL cp_fm_release(fm_matrix_m_inv_work)
315  END IF
316 
317  IF (.NOT. do_im_time) THEN
318 
319  CALL cp_fm_to_fm(fm_matrix_minv_l_kpoints(1, 1), fm_matrix_l_work)
320  CALL cp_fm_release(fm_matrix_minv_l_kpoints)
321 
322  END IF
323 
324  END IF
325 
326  ! If the group distribution changed because of an SVD, we get the new values here
327  CALL get_group_dist(gd_array, color_sub, my_group_l_start, my_group_l_end, my_group_l_size)
328 
329  IF (.NOT. do_im_time) THEN
330  IF (unit_nr > 0) THEN
331  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") "RI_INFO| Cholesky decomposition group size:", para_env_l%num_pe
332  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") "RI_INFO| Number of groups for auxiliary basis functions", ngroup
333  IF (calc_pq_cond_num .OR. do_svd) THEN
334  WRITE (unit=unit_nr, fmt="(T3,A,T67,ES14.5)") &
335  "RI_INFO| Condition number of the (P|Q):", cond_num
336  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
337  "RI_INFO| Number of non-positive Eigenvalues of (P|Q):", num_small_eigen
338  END IF
339  CALL m_flush(unit_nr)
340  END IF
341 
342  ! replicate the necessary row of the L^{-1} matrix on each proc
343  CALL grep_lcols(fm_matrix_l_work, my_group_l_start, my_group_l_end, my_group_l_size, my_lrows)
344  IF (calc_forces .AND. .NOT. compare_potential_types(qs_env%mp2_env%ri_metric, &
345  qs_env%mp2_env%potential_parameter)) THEN
346  CALL grep_lcols(fm_matrix_v, my_group_l_start, my_group_l_end, my_group_l_size, my_vrows)
347  END IF
348  END IF
349  CALL cp_fm_release(fm_matrix_l_work)
350  IF (calc_forces .AND. .NOT. (do_im_time .OR. compare_potential_types(qs_env%mp2_env%ri_metric, &
351  qs_env%mp2_env%potential_parameter))) CALL cp_fm_release(fm_matrix_v)
352  CALL mp_para_env_release(para_env_l)
353 
354  CALL timestop(handle)
355 
356  END SUBROUTINE get_2c_integrals
357 
358 ! **************************************************************************************************
359 !> \brief ...
360 !> \param fm_matrix_L ...
361 !> \param do_svd ...
362 !> \param num_small_eigen ...
363 !> \param cond_num ...
364 !> \param do_inversion ...
365 !> \param gd_array ...
366 !> \param ngroup ...
367 !> \param dimen_RI ...
368 !> \param dimen_RI_red ...
369 !> \param para_env ...
370 ! **************************************************************************************************
371  SUBROUTINE decomp_mat_l(fm_matrix_L, do_svd, num_small_eigen, cond_num, do_inversion, gd_array, ngroup, &
372  dimen_RI, dimen_RI_red, para_env)
373 
374  TYPE(cp_fm_type), INTENT(INOUT) :: fm_matrix_l
375  LOGICAL, INTENT(IN) :: do_svd
376  INTEGER, INTENT(INOUT) :: num_small_eigen
377  REAL(kind=dp), INTENT(INOUT) :: cond_num
378  LOGICAL, INTENT(IN) :: do_inversion
379  TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_array
380  INTEGER, INTENT(IN) :: ngroup, dimen_ri
381  INTEGER, INTENT(INOUT) :: dimen_ri_red
382  TYPE(mp_para_env_type), INTENT(IN) :: para_env
383 
384  IF (do_svd) THEN
385  CALL matrix_root_with_svd(fm_matrix_l, num_small_eigen, cond_num, do_inversion, para_env)
386 
387  dimen_ri_red = dimen_ri - num_small_eigen
388 
389  ! We changed the size of fm_matrix_L in matrix_root_with_svd.
390  ! So, we have to get new group sizes
391  CALL release_group_dist(gd_array)
392 
393  CALL create_group_dist(gd_array, ngroup, dimen_ri_red)
394  ELSE
395 
396  ! calculate Cholesky decomposition L (V=LL^T)
397  CALL cholesky_decomp(fm_matrix_l, dimen_ri, do_inversion=do_inversion)
398 
399  IF (do_inversion) CALL invert_mat(fm_matrix_l)
400  END IF
401 
402  END SUBROUTINE decomp_mat_l
403 
404 ! **************************************************************************************************
405 !> \brief ...
406 !> \param qs_env ...
407 !> \param fm_matrix_L_kpoints ...
408 !> \param fm_matrix_Minv_L_kpoints ...
409 !> \param kpoints ...
410 ! **************************************************************************************************
411  SUBROUTINE compute_v_by_lattice_sum(qs_env, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, kpoints)
412  TYPE(qs_environment_type), POINTER :: qs_env
413  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_l_kpoints, &
414  fm_matrix_minv_l_kpoints
415  TYPE(kpoint_type), POINTER :: kpoints
416 
417  CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_V_by_lattice_sum'
418 
419  INTEGER :: handle, i_dim, i_real_imag, ikp, nkp, &
420  nkp_extra, nkp_orig
421  INTEGER, DIMENSION(3) :: nkp_grid_orig, periodic
422  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
423  TYPE(cell_type), POINTER :: cell
424  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_ri_aux_transl, matrix_v_ri_kp
425  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
426  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
427 
428  CALL timeset(routinen, handle)
429 
430  NULLIFY (matrix_s_ri_aux_transl, particle_set, cell, qs_kind_set)
431 
432  CALL get_qs_env(qs_env=qs_env, &
433  matrix_s_ri_aux_kp=matrix_s_ri_aux_transl, &
434  particle_set=particle_set, &
435  cell=cell, &
436  qs_kind_set=qs_kind_set, &
437  atomic_kind_set=atomic_kind_set)
438 
439  ! check that we have a 2n x 2m x 2k mesh to guarantee good convergence
440  CALL get_cell(cell=cell, periodic=periodic)
441  DO i_dim = 1, 3
442  IF (periodic(i_dim) == 1) THEN
443  cpassert(modulo(kpoints%nkp_grid(i_dim), 2) == 0)
444  END IF
445  END DO
446 
447  nkp = kpoints%nkp
448 
449  ALLOCATE (fm_matrix_l_kpoints(nkp, 2))
450  DO ikp = 1, nkp
451  DO i_real_imag = 1, 2
452  CALL cp_fm_create(fm_matrix_l_kpoints(ikp, i_real_imag), &
453  fm_matrix_minv_l_kpoints(1, i_real_imag)%matrix_struct)
454  CALL cp_fm_set_all(fm_matrix_l_kpoints(ikp, i_real_imag), 0.0_dp)
455  END DO
456  END DO
457 
458  CALL allocate_matrix_v_ri_kp(matrix_v_ri_kp, matrix_s_ri_aux_transl, nkp)
459 
460  IF (qs_env%mp2_env%ri_rpa_im_time%do_extrapolate_kpoints) THEN
461 
462  nkp_orig = qs_env%mp2_env%ri_rpa_im_time%nkp_orig
463  nkp_extra = qs_env%mp2_env%ri_rpa_im_time%nkp_extra
464 
465  CALL build_2c_coulomb_matrix_kp(matrix_v_ri_kp, kpoints, basis_type="RI_AUX", &
466  cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
467  atomic_kind_set=atomic_kind_set, &
468  size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
469  operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp_orig)
470 
471  nkp_grid_orig = kpoints%nkp_grid
472  kpoints%nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid_extra(1:3)
473 
474  CALL build_2c_coulomb_matrix_kp(matrix_v_ri_kp, kpoints, basis_type="RI_AUX", &
475  cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
476  atomic_kind_set=atomic_kind_set, &
477  size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
478  operator_type=operator_coulomb, ikp_start=nkp_orig + 1, ikp_end=nkp)
479 
480  kpoints%nkp_grid = nkp_grid_orig
481 
482  ELSE
483 
484  CALL build_2c_coulomb_matrix_kp(matrix_v_ri_kp, kpoints, basis_type="RI_AUX", &
485  cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
486  atomic_kind_set=atomic_kind_set, &
487  size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
488  operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp)
489 
490  END IF
491 
492  DO ikp = 1, nkp
493 
494  CALL copy_dbcsr_to_fm(matrix_v_ri_kp(ikp, 1)%matrix, fm_matrix_l_kpoints(ikp, 1))
495  CALL copy_dbcsr_to_fm(matrix_v_ri_kp(ikp, 2)%matrix, fm_matrix_l_kpoints(ikp, 2))
496 
497  END DO
498 
499  CALL dbcsr_deallocate_matrix_set(matrix_v_ri_kp)
500 
501  CALL timestop(handle)
502 
503  END SUBROUTINE compute_v_by_lattice_sum
504 
505 ! **************************************************************************************************
506 !> \brief ...
507 !> \param matrix_v_RI_kp ...
508 !> \param matrix_s_RI_aux_transl ...
509 !> \param nkp ...
510 ! **************************************************************************************************
511  SUBROUTINE allocate_matrix_v_ri_kp(matrix_v_RI_kp, matrix_s_RI_aux_transl, nkp)
512 
513  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_ri_kp, matrix_s_ri_aux_transl
514  INTEGER :: nkp
515 
516  INTEGER :: ikp
517 
518  NULLIFY (matrix_v_ri_kp)
519  CALL dbcsr_allocate_matrix_set(matrix_v_ri_kp, nkp, 2)
520 
521  DO ikp = 1, nkp
522 
523  ALLOCATE (matrix_v_ri_kp(ikp, 1)%matrix)
524  CALL dbcsr_create(matrix_v_ri_kp(ikp, 1)%matrix, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
525  matrix_type=dbcsr_type_no_symmetry)
526  CALL dbcsr_reserve_all_blocks(matrix_v_ri_kp(ikp, 1)%matrix)
527  CALL dbcsr_set(matrix_v_ri_kp(ikp, 1)%matrix, 0.0_dp)
528 
529  ALLOCATE (matrix_v_ri_kp(ikp, 2)%matrix)
530  CALL dbcsr_create(matrix_v_ri_kp(ikp, 2)%matrix, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
531  matrix_type=dbcsr_type_no_symmetry)
532  CALL dbcsr_reserve_all_blocks(matrix_v_ri_kp(ikp, 2)%matrix)
533  CALL dbcsr_set(matrix_v_ri_kp(ikp, 2)%matrix, 0.0_dp)
534 
535  END DO
536 
537  END SUBROUTINE allocate_matrix_v_ri_kp
538 
539 ! **************************************************************************************************
540 !> \brief ...
541 !> \param qs_env ...
542 !> \param fm_matrix_Minv_L_kpoints ...
543 !> \param fm_matrix_L ...
544 !> \param dimen_RI ...
545 !> \param ri_metric ...
546 !> \param do_kpoints ...
547 !> \param kpoints ...
548 !> \param put_mat_KS_env ...
549 !> \param regularization_RI ...
550 !> \param ikp_ext ...
551 ! **************************************************************************************************
552  SUBROUTINE ri_2c_integral_mat(qs_env, fm_matrix_Minv_L_kpoints, fm_matrix_L, dimen_RI, ri_metric, &
553  do_kpoints, kpoints, put_mat_KS_env, regularization_RI, ikp_ext)
554 
555  TYPE(qs_environment_type), POINTER :: qs_env
556  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_minv_l_kpoints
557  TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_l
558  INTEGER, INTENT(IN) :: dimen_ri
559  TYPE(libint_potential_type) :: ri_metric
560  LOGICAL, INTENT(IN) :: do_kpoints
561  TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
562  LOGICAL, OPTIONAL :: put_mat_ks_env
563  REAL(kind=dp), OPTIONAL :: regularization_ri
564  INTEGER, OPTIONAL :: ikp_ext
565 
566  CHARACTER(LEN=*), PARAMETER :: routinen = 'RI_2c_integral_mat'
567 
568  INTEGER :: handle, i_real_imag, i_size, ikp, &
569  ikp_for_xkp, img, n_real_imag, natom, &
570  nimg, nkind, nkp
571  INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_ri
572  INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
573  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
574  LOGICAL :: my_put_mat_ks_env
575  REAL(kind=dp) :: my_regularization_ri
576  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
577  TYPE(cp_blacs_env_type), POINTER :: blacs_env
578  TYPE(cp_fm_struct_type), POINTER :: fm_struct
579  TYPE(cp_fm_type) :: fm_matrix_s_global
580  TYPE(dbcsr_distribution_type) :: dbcsr_dist
581  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_ri_aux_transl
582  TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: mat_2c
583  TYPE(dbcsr_type), POINTER :: cmatrix, matrix_s_ri_aux_desymm, &
584  rmatrix, tmpmat
585  TYPE(dft_control_type), POINTER :: dft_control
586  TYPE(distribution_2d_type), POINTER :: dist_2d
587  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ri
588  TYPE(mp_para_env_type), POINTER :: para_env
589  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
590  POINTER :: sab_ri
591  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
592  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
593 
594  CALL timeset(routinen, handle)
595 
596  NULLIFY (sab_ri, matrix_s_ri_aux_transl, dist_2d)
597 
598  IF (PRESENT(regularization_ri)) THEN
599  my_regularization_ri = regularization_ri
600  ELSE
601  my_regularization_ri = 0.0_dp
602  END IF
603 
604  IF (PRESENT(put_mat_ks_env)) THEN
605  my_put_mat_ks_env = put_mat_ks_env
606  ELSE
607  my_put_mat_ks_env = .false.
608  END IF
609 
610  CALL get_qs_env(qs_env=qs_env, &
611  para_env=para_env, &
612  blacs_env=blacs_env, &
613  nkind=nkind, &
614  distribution_2d=dist_2d, &
615  qs_kind_set=qs_kind_set, &
616  particle_set=particle_set, &
617  dft_control=dft_control, &
618  natom=natom)
619 
620  ALLOCATE (sizes_ri(natom))
621  ALLOCATE (basis_set_ri(nkind))
622 
623  CALL basis_set_list_setup(basis_set_ri, 'RI_AUX', qs_kind_set)
624  CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri)
625 
626  CALL build_2c_neighbor_lists(sab_ri, basis_set_ri, basis_set_ri, ri_metric, "RPA_2c_nl_RI", qs_env, &
627  sym_ij=.true., dist_2d=dist_2d)
628 
629  CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
630  ALLOCATE (row_bsize(SIZE(sizes_ri)))
631  ALLOCATE (col_bsize(SIZE(sizes_ri)))
632  row_bsize(:) = sizes_ri
633  col_bsize(:) = sizes_ri
634 
635  IF (do_kpoints) THEN
636  nimg = dft_control%nimages
637  ELSE
638  nimg = 1
639  END IF
640 
641  ALLOCATE (mat_2c(nimg))
642  CALL dbcsr_create(mat_2c(1), "(RI|RI)", dbcsr_dist, dbcsr_type_symmetric, &
643  row_bsize, col_bsize)
644  DEALLOCATE (row_bsize, col_bsize)
645 
646  DO img = 2, nimg
647  CALL dbcsr_create(mat_2c(img), template=mat_2c(1))
648  END DO
649 
650  CALL build_2c_integrals(mat_2c, 0.0_dp, qs_env, sab_ri, basis_set_ri, basis_set_ri, &
651  ri_metric, do_kpoints=do_kpoints, ext_kpoints=kpoints, &
652  regularization_ri=my_regularization_ri)
653 
654  CALL dbcsr_distribution_release(dbcsr_dist)
655  DEALLOCATE (basis_set_ri)
656 
657  IF (my_put_mat_ks_env) THEN
658  CALL get_ks_env(qs_env%ks_env, matrix_s_ri_aux_kp=matrix_s_ri_aux_transl)
659  IF (ASSOCIATED(matrix_s_ri_aux_transl)) CALL dbcsr_deallocate_matrix_set(matrix_s_ri_aux_transl)
660  END IF
661  CALL dbcsr_allocate_matrix_set(matrix_s_ri_aux_transl, 1, nimg)
662  DO img = 1, nimg
663  ALLOCATE (matrix_s_ri_aux_transl(1, img)%matrix)
664  CALL dbcsr_copy(matrix_s_ri_aux_transl(1, img)%matrix, mat_2c(img))
665  CALL dbcsr_release(mat_2c(img))
666  END DO
667 
668  IF (my_put_mat_ks_env) THEN
669  CALL set_ks_env(qs_env%ks_env, matrix_s_ri_aux_kp=matrix_s_ri_aux_transl)
670  END IF
671 
672  IF (do_kpoints) THEN
673  cpassert(PRESENT(kpoints))
674  CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, &
675  cell_to_index=cell_to_index)
676  n_real_imag = 2
677  ELSE
678  nkp = 1
679  n_real_imag = 1
680  END IF
681 
682  IF (PRESENT(ikp_ext)) nkp = 1
683 
684  ALLOCATE (fm_matrix_minv_l_kpoints(nkp, n_real_imag))
685  DO i_size = 1, nkp
686  DO i_real_imag = 1, n_real_imag
687  CALL cp_fm_create(fm_matrix_minv_l_kpoints(i_size, i_real_imag), fm_matrix_l%matrix_struct)
688  CALL cp_fm_set_all(fm_matrix_minv_l_kpoints(i_size, i_real_imag), 0.0_dp)
689  END DO
690  END DO
691 
692  NULLIFY (fm_struct)
693  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_ri, &
694  ncol_global=dimen_ri, para_env=para_env)
695 
696  CALL cp_fm_create(fm_matrix_s_global, fm_struct)
697  CALL cp_fm_set_all(fm_matrix_s_global, 0.0_dp)
698 
699  IF (do_kpoints) THEN
700 
701  ALLOCATE (rmatrix, cmatrix, tmpmat)
702  CALL dbcsr_create(rmatrix, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
703  matrix_type=dbcsr_type_symmetric)
704  CALL dbcsr_create(cmatrix, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
705  matrix_type=dbcsr_type_antisymmetric)
706  CALL dbcsr_create(tmpmat, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
707  matrix_type=dbcsr_type_no_symmetry)
708  CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_ri)
709  CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_ri)
710 
711  DO ikp = 1, nkp
712 
713  ! s matrix is not spin dependent, double the work
714  CALL dbcsr_set(rmatrix, 0.0_dp)
715  CALL dbcsr_set(cmatrix, 0.0_dp)
716 
717  IF (PRESENT(ikp_ext)) THEN
718  ikp_for_xkp = ikp_ext
719  ELSE
720  ikp_for_xkp = ikp
721  END IF
722 
723  CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s_ri_aux_transl, ispin=1, &
724  xkp=xkp(1:3, ikp_for_xkp), cell_to_index=cell_to_index, sab_nl=sab_ri)
725 
726  CALL dbcsr_set(tmpmat, 0.0_dp)
727  CALL dbcsr_desymmetrize(rmatrix, tmpmat)
728 
729  CALL cp_fm_set_all(fm_matrix_s_global, 0.0_dp)
730  CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_s_global)
731  CALL cp_fm_copy_general(fm_matrix_s_global, fm_matrix_minv_l_kpoints(ikp, 1), para_env)
732 
733  CALL dbcsr_set(tmpmat, 0.0_dp)
734  CALL dbcsr_desymmetrize(cmatrix, tmpmat)
735 
736  CALL cp_fm_set_all(fm_matrix_s_global, 0.0_dp)
737  CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_s_global)
738  CALL cp_fm_copy_general(fm_matrix_s_global, fm_matrix_minv_l_kpoints(ikp, 2), para_env)
739 
740  END DO
741 
742  CALL dbcsr_deallocate_matrix(rmatrix)
743  CALL dbcsr_deallocate_matrix(cmatrix)
744  CALL dbcsr_deallocate_matrix(tmpmat)
745 
746  ELSE
747 
748  NULLIFY (matrix_s_ri_aux_desymm)
749  ALLOCATE (matrix_s_ri_aux_desymm)
750  CALL dbcsr_create(matrix_s_ri_aux_desymm, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
751  name='S_RI non_symm', matrix_type=dbcsr_type_no_symmetry)
752 
753  CALL dbcsr_desymmetrize(matrix_s_ri_aux_transl(1, 1)%matrix, matrix_s_ri_aux_desymm)
754 
755  CALL copy_dbcsr_to_fm(matrix_s_ri_aux_desymm, fm_matrix_s_global)
756 
757  CALL cp_fm_copy_general(fm_matrix_s_global, fm_matrix_minv_l_kpoints(1, 1), para_env)
758 
759  CALL dbcsr_deallocate_matrix(matrix_s_ri_aux_desymm)
760 
761  END IF
762 
763  CALL release_neighbor_list_sets(sab_ri)
764 
765  CALL cp_fm_struct_release(fm_struct)
766 
767  CALL cp_fm_release(fm_matrix_s_global)
768 
769  IF (.NOT. my_put_mat_ks_env) THEN
770  CALL dbcsr_deallocate_matrix_set(matrix_s_ri_aux_transl)
771  END IF
772 
773  CALL timestop(handle)
774 
775  END SUBROUTINE ri_2c_integral_mat
776 
777 ! **************************************************************************************************
778 !> \brief ...
779 !> \param qs_env ...
780 !> \param eri_method ...
781 !> \param eri_param ...
782 !> \param para_env ...
783 !> \param para_env_sub ...
784 !> \param para_env_L ...
785 !> \param mp2_memory ...
786 !> \param fm_matrix_L ...
787 !> \param ngroup ...
788 !> \param color_sub ...
789 !> \param dimen_RI ...
790 !> \param my_group_L_size ...
791 !> \param my_group_L_start ...
792 !> \param my_group_L_end ...
793 !> \param gd_array ...
794 !> \param calc_PQ_cond_num ...
795 !> \param cond_num ...
796 !> \param num_small_eigen ...
797 !> \param potential ...
798 !> \param sab_orb_sub ...
799 !> \param do_im_time ...
800 !> \param fm_matrix_L_extern ...
801 ! **************************************************************************************************
802  SUBROUTINE compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
803  fm_matrix_L, ngroup, color_sub, dimen_RI, &
804  my_group_L_size, my_group_L_start, my_group_L_end, &
805  gd_array, calc_PQ_cond_num, cond_num, num_small_eigen, potential, &
806  sab_orb_sub, do_im_time, fm_matrix_L_extern)
807 
808  TYPE(qs_environment_type), POINTER :: qs_env
809  INTEGER, INTENT(IN) :: eri_method
810  TYPE(cp_eri_mme_param), POINTER :: eri_param
811  TYPE(mp_para_env_type), INTENT(IN) :: para_env
812  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
813  TYPE(mp_para_env_type), INTENT(OUT), POINTER :: para_env_l
814  REAL(kind=dp), INTENT(IN) :: mp2_memory
815  TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_l
816  INTEGER, INTENT(IN) :: ngroup, color_sub
817  INTEGER, INTENT(OUT) :: dimen_ri, my_group_l_size, &
818  my_group_l_start, my_group_l_end
819  TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array
820  LOGICAL, INTENT(IN) :: calc_pq_cond_num
821  REAL(kind=dp), INTENT(OUT) :: cond_num
822  INTEGER, INTENT(OUT) :: num_small_eigen
823  TYPE(libint_potential_type), INTENT(IN) :: potential
824  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
825  POINTER :: sab_orb_sub
826  LOGICAL, INTENT(IN), OPTIONAL :: do_im_time
827  TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_matrix_l_extern
828 
829  CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_2c_integrals'
830 
831  INTEGER :: best_group_size, color_l, group_size, handle, handle2, i_global, iatom, iib, &
832  ikind, iproc, j_global, jjb, natom, ncol_local, nkind, nrow_local, nsgf, potential_type, &
833  proc_receive, proc_receive_static, proc_send_static, proc_shift, rec_l_end, rec_l_size, &
834  rec_l_start, strat_group_size
835  INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
836  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
837  LOGICAL :: my_do_im_time
838  REAL(kind=dp) :: min_mem_for_qk
839  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: egen_l
840  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: l_external_col, l_local_col
841  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
842  TYPE(cp_blacs_env_type), POINTER :: blacs_env_l
843  TYPE(cp_fm_type) :: fm_matrix_l_diag
844  TYPE(group_dist_d1_type) :: gd_sub_array
845  TYPE(gto_basis_set_type), POINTER :: basis_set_a
846  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
847  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
848 
849  CALL timeset(routinen, handle)
850 
851  my_do_im_time = .false.
852  IF (PRESENT(do_im_time)) THEN
853  my_do_im_time = do_im_time
854  END IF
855 
856  ! get stuff
857  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
858  particle_set=particle_set)
859 
860  nkind = SIZE(qs_kind_set)
861  natom = SIZE(particle_set)
862 
863  CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
864 
865  DO ikind = 1, nkind
866  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
867  cpassert(ASSOCIATED(basis_set_a))
868  END DO
869 
870  dimen_ri = 0
871  DO iatom = 1, natom
872  ikind = kind_of(iatom)
873  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
874  dimen_ri = dimen_ri + nsgf
875  END DO
876 
877  ! check that very small systems are not running on too many processes
878  IF (dimen_ri < ngroup) THEN
879  CALL cp_abort(__location__, "Product of block size and number "// &
880  "of RI functions should not exceed total number of processes")
881  END IF
882 
883  CALL create_group_dist(gd_array, ngroup, dimen_ri)
884 
885  CALL get_group_dist(gd_array, color_sub, my_group_l_start, my_group_l_end, my_group_l_size)
886 
887  CALL timeset(routinen//"_loop_lm", handle2)
888 
889  ALLOCATE (l_local_col(dimen_ri, my_group_l_size))
890  l_local_col = 0.0_dp
891 
892  potential_type = potential%potential_type
893 
894  IF ((eri_method .EQ. do_eri_mme .OR. eri_method .EQ. do_eri_os) &
895  .AND. potential_type .EQ. do_potential_coulomb .OR. &
896  (eri_method .EQ. do_eri_mme .AND. potential_type .EQ. do_potential_long)) THEN
897 
898  CALL mp2_eri_2c_integrate(eri_param, potential, para_env_sub, qs_env, &
899  basis_type_a="RI_AUX", basis_type_b="RI_AUX", &
900  hab=l_local_col, first_b=my_group_l_start, last_b=my_group_l_end, &
901  eri_method=eri_method)
902 
903  ELSEIF (eri_method .EQ. do_eri_gpw .OR. &
904  (potential_type == do_potential_long .AND. qs_env%mp2_env%eri_method == do_eri_os) &
905  .OR. (potential_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme)) THEN
906 
907  CALL mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_l_start, my_group_l_end, &
908  natom, potential, sab_orb_sub, l_local_col, kind_of)
909 
910  ELSE
911  cpabort("unknown ERI method")
912  END IF
913 
914  CALL timestop(handle2)
915 
916  ! split the total number of proc in a subgroup of the size of ~1/10 of the
917  ! total num of proc
918  best_group_size = para_env%num_pe
919 
920  strat_group_size = max(1, para_env%num_pe/10)
921 
922  min_mem_for_qk = real(dimen_ri, kind=dp)*dimen_ri*3.0_dp*8.0_dp/1024_dp/1024_dp
923 
924  group_size = strat_group_size - 1
925  DO iproc = strat_group_size, para_env%num_pe
926  group_size = group_size + 1
927  ! check that group_size is a multiple of sub_group_size and a divisor of
928  ! the total num of proc
929  IF (mod(para_env%num_pe, group_size) /= 0 .OR. mod(group_size, para_env_sub%num_pe) /= 0) cycle
930 
931  ! check for memory
932  IF (real(group_size, kind=dp)*mp2_memory < min_mem_for_qk) cycle
933 
934  best_group_size = group_size
935  EXIT
936  END DO
937 
938  IF (my_do_im_time) THEN
939  ! para_env_L is the para_env for im. time to avoid bug
940  best_group_size = para_env%num_pe
941  END IF
942 
943  ! create the L group
944  color_l = para_env%mepos/best_group_size
945  ALLOCATE (para_env_l)
946  CALL para_env_l%from_split(para_env, color_l)
947 
948  ! create the blacs_L
949  NULLIFY (blacs_env_l)
950  CALL cp_blacs_env_create(blacs_env=blacs_env_l, para_env=para_env_l)
951 
952  ! create the full matrix L defined in the L group
953  CALL create_matrix_l(fm_matrix_l, blacs_env_l, dimen_ri, para_env_l, "fm_matrix_L", fm_matrix_l_extern)
954 
955  IF (my_do_im_time .AND. para_env%num_pe > 1) THEN
956 
957  CALL fill_fm_l_from_l_loc_non_blocking(fm_matrix_l, l_local_col, para_env, &
958  my_group_l_start, my_group_l_end, &
959  dimen_ri)
960 
961  ELSE
962  block
963  TYPE(mp_comm_type) :: comm_exchange
964 
965  CALL comm_exchange%from_split(para_env_l, para_env_sub%mepos)
966 
967  CALL create_group_dist(gd_sub_array, my_group_l_start, &
968  my_group_l_end, my_group_l_size, comm_exchange)
969 
970  CALL cp_fm_get_info(matrix=fm_matrix_l, &
971  nrow_local=nrow_local, &
972  ncol_local=ncol_local, &
973  row_indices=row_indices, &
974  col_indices=col_indices)
975 
976  DO jjb = 1, ncol_local
977  j_global = col_indices(jjb)
978  IF (j_global >= my_group_l_start .AND. j_global <= my_group_l_end) THEN
979  DO iib = 1, nrow_local
980  i_global = row_indices(iib)
981  fm_matrix_l%local_data(iib, jjb) = l_local_col(i_global, j_global - my_group_l_start + 1)
982  END DO
983  END IF
984  END DO
985 
986  proc_send_static = modulo(comm_exchange%mepos + 1, comm_exchange%num_pe)
987  proc_receive_static = modulo(comm_exchange%mepos - 1, comm_exchange%num_pe)
988 
989  DO proc_shift = 1, comm_exchange%num_pe - 1
990  proc_receive = modulo(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
991 
992  CALL get_group_dist(gd_sub_array, proc_receive, rec_l_start, rec_l_end, rec_l_size)
993 
994  ALLOCATE (l_external_col(dimen_ri, rec_l_size))
995  l_external_col = 0.0_dp
996 
997  CALL comm_exchange%sendrecv(l_local_col, proc_send_static, l_external_col, proc_receive_static)
998 
999  DO jjb = 1, ncol_local
1000  j_global = col_indices(jjb)
1001  IF (j_global >= rec_l_start .AND. j_global <= rec_l_end) THEN
1002  DO iib = 1, nrow_local
1003  i_global = row_indices(iib)
1004  fm_matrix_l%local_data(iib, jjb) = l_external_col(i_global, j_global - rec_l_start + 1)
1005  END DO
1006  END IF
1007  END DO
1008 
1009  CALL move_alloc(l_external_col, l_local_col)
1010  END DO
1011 
1012  CALL release_group_dist(gd_sub_array)
1013  CALL comm_exchange%free()
1014  END block
1015  END IF
1016 
1017  DEALLOCATE (l_local_col)
1018 
1019  ! create the new group for the sum of the local data over all processes
1020  block
1021  TYPE(mp_comm_type) :: comm_exchange
1022  comm_exchange = fm_matrix_l%matrix_struct%context%interconnect(para_env)
1023  CALL comm_exchange%sum(fm_matrix_l%local_data)
1024  CALL comm_exchange%free()
1025  END block
1026 
1027  cond_num = 1.0_dp
1028  num_small_eigen = 0
1029  IF (calc_pq_cond_num) THEN
1030  ! calculate the condition number of the (P|Q) matrix
1031  ! create a copy of the matrix
1032 
1033  CALL create_matrix_l(fm_matrix_l_diag, blacs_env_l, dimen_ri, para_env_l, "fm_matrix_L_diag", fm_matrix_l_extern)
1034 
1035  CALL cp_fm_to_fm(source=fm_matrix_l, destination=fm_matrix_l_diag)
1036 
1037  ALLOCATE (egen_l(dimen_ri))
1038 
1039  egen_l = 0.0_dp
1040  CALL cp_fm_syevx(matrix=fm_matrix_l_diag, eigenvalues=egen_l)
1041 
1042  num_small_eigen = 0
1043  DO iib = 1, dimen_ri
1044  IF (abs(egen_l(iib)) < 0.001_dp) num_small_eigen = num_small_eigen + 1
1045  END DO
1046 
1047  cond_num = maxval(abs(egen_l))/minval(abs(egen_l))
1048 
1049  CALL cp_fm_release(fm_matrix_l_diag)
1050 
1051  DEALLOCATE (egen_l)
1052  END IF
1053 
1054  ! release blacs_env
1055  CALL cp_blacs_env_release(blacs_env_l)
1056 
1057  CALL timestop(handle)
1058 
1059  END SUBROUTINE compute_2c_integrals
1060 
1061 ! **************************************************************************************************
1062 !> \brief ...
1063 !> \param matrix ...
1064 !> \param num_small_evals ...
1065 !> \param cond_num ...
1066 !> \param do_inversion ...
1067 !> \param para_env ...
1068 ! **************************************************************************************************
1069  SUBROUTINE matrix_root_with_svd(matrix, num_small_evals, cond_num, do_inversion, para_env)
1070  TYPE(cp_fm_type), INTENT(INOUT) :: matrix
1071  INTEGER, INTENT(OUT) :: num_small_evals
1072  REAL(kind=dp), INTENT(OUT) :: cond_num
1073  LOGICAL, INTENT(IN) :: do_inversion
1074  TYPE(mp_para_env_type), INTENT(IN) :: para_env
1075 
1076  CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_root_with_svd'
1077 
1078  INTEGER :: group_size_l, handle, ii, needed_evals, &
1079  nrow, pos_max(1)
1080  INTEGER, ALLOCATABLE, DIMENSION(:) :: num_eval
1081  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
1082  TYPE(cp_fm_type) :: evecs
1083  TYPE(mp_comm_type) :: comm_exchange
1084 
1085  CALL timeset(routinen, handle)
1086 
1087  ! Create arrays carrying eigenvalues and eigenvectors later
1088  CALL cp_fm_get_info(matrix=matrix, nrow_global=nrow)
1089 
1090  ALLOCATE (evals(nrow))
1091  evals = 0
1092 
1093  CALL cp_fm_create(evecs, matrix%matrix_struct)
1094 
1095  ! Perform the EVD (=SVD of a positive semidefinite matrix)
1096  CALL choose_eigv_solver(matrix, evecs, evals)
1097 
1098  ! Determine the number of neglectable eigenvalues assuming that the eigenvalues are in ascending order
1099  num_small_evals = 0
1100  DO ii = 1, nrow
1101  IF (evals(ii) > 0.0_dp) THEN
1102  num_small_evals = ii - 1
1103  EXIT
1104  END IF
1105  END DO
1106  needed_evals = nrow - num_small_evals
1107 
1108  ! Get the condition number w.r.t. considered eigenvalues
1109  cond_num = evals(nrow)/evals(num_small_evals + 1)
1110 
1111  ! Determine the eigenvalues of the request matrix root or its inverse
1112  evals(1:num_small_evals) = 0.0_dp
1113  IF (do_inversion) THEN
1114  evals(num_small_evals + 1:nrow) = 1.0_dp/sqrt(evals(num_small_evals + 1:nrow))
1115  ELSE
1116  evals(num_small_evals + 1:nrow) = sqrt(evals(num_small_evals + 1:nrow))
1117  END IF
1118 
1119  CALL cp_fm_column_scale(evecs, evals)
1120 
1121  ! As it turns out, the results in the subgroups differ.
1122  ! Thus, we have to equilibrate the results.
1123  ! Step 1: Get a communicator connecting processes with same id within subgroups
1124  ! use that group_size_L is divisible by the total number of processes (see above)
1125  group_size_l = para_env%num_pe/matrix%matrix_struct%para_env%num_pe
1126  comm_exchange = matrix%matrix_struct%context%interconnect(para_env)
1127 
1128  ALLOCATE (num_eval(0:group_size_l - 1))
1129  CALL comm_exchange%allgather(num_small_evals, num_eval)
1130 
1131  num_small_evals = minval(num_eval)
1132 
1133  IF (num_small_evals /= maxval(num_eval)) THEN
1134  ! Step 2: Get position of maximum value
1135  pos_max = maxloc(num_eval)
1136  num_small_evals = num_eval(pos_max(1))
1137  needed_evals = nrow - num_small_evals
1138 
1139  ! Step 3: Broadcast your local data to all other processes
1140  CALL comm_exchange%bcast(evecs%local_data, pos_max(1))
1141  CALL comm_exchange%bcast(cond_num, pos_max(1))
1142  END IF
1143 
1144  DEALLOCATE (num_eval)
1145 
1146  CALL comm_exchange%free()
1147 
1148  CALL reset_size_matrix(matrix, needed_evals, matrix%matrix_struct)
1149 
1150  ! Copy the needed eigenvectors
1151  CALL cp_fm_to_fm(evecs, matrix, needed_evals, num_small_evals + 1)
1152 
1153  CALL cp_fm_release(evecs)
1154 
1155  CALL timestop(handle)
1156 
1157  END SUBROUTINE matrix_root_with_svd
1158 
1159 ! **************************************************************************************************
1160 !> \brief ...
1161 !> \param matrix ...
1162 !> \param new_size ...
1163 !> \param fm_struct_template ...
1164 ! **************************************************************************************************
1165  SUBROUTINE reset_size_matrix(matrix, new_size, fm_struct_template)
1166  TYPE(cp_fm_type), INTENT(INOUT) :: matrix
1167  INTEGER, INTENT(IN) :: new_size
1168  TYPE(cp_fm_struct_type), POINTER :: fm_struct_template
1169 
1170  CHARACTER(LEN=*), PARAMETER :: routinen = 'reset_size_matrix'
1171 
1172  INTEGER :: handle
1173  TYPE(cp_fm_struct_type), POINTER :: fm_struct
1174 
1175  CALL timeset(routinen, handle)
1176 
1177  ! Choose the block sizes as large as in the template for the later steps
1178  NULLIFY (fm_struct)
1179  CALL cp_fm_struct_create(fm_struct, ncol_global=new_size, template_fmstruct=fm_struct_template, force_block=.true.)
1180 
1181  CALL cp_fm_release(matrix)
1182 
1183  CALL cp_fm_create(matrix, fm_struct)
1184  CALL cp_fm_set_all(matrix, 0.0_dp)
1185 
1186  CALL cp_fm_struct_release(fm_struct)
1187 
1188  CALL timestop(handle)
1189 
1190  END SUBROUTINE reset_size_matrix
1191 
1192 ! **************************************************************************************************
1193 !> \brief ...
1194 !> \param fm_matrix_L ...
1195 !> \param L_local_col ...
1196 !> \param para_env ...
1197 !> \param my_group_L_start ...
1198 !> \param my_group_L_end ...
1199 !> \param dimen_RI ...
1200 ! **************************************************************************************************
1201  SUBROUTINE fill_fm_l_from_l_loc_non_blocking(fm_matrix_L, L_local_col, para_env, my_group_L_start, &
1202  my_group_L_end, dimen_RI)
1203  TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_l
1204  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
1205  INTENT(IN) :: l_local_col
1206  TYPE(mp_para_env_type), INTENT(IN) :: para_env
1207  INTEGER, INTENT(IN) :: my_group_l_start, my_group_l_end, &
1208  dimen_ri
1209 
1210  CHARACTER(LEN=*), PARAMETER :: routinen = 'fill_fm_L_from_L_loc_non_blocking'
1211 
1212  INTEGER :: dummy_proc, handle, handle2, i_entry_rec, i_row, i_row_global, iproc, j_col, &
1213  j_col_global, lll, mmm, ncol_block, ncol_local, npcol, nprow, nrow_block, nrow_local, &
1214  proc_send, send_pcol, send_prow
1215  INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, &
1216  num_entries_send
1217  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1218  TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1219  DIMENSION(:) :: buffer_rec, buffer_send
1220  TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
1221 
1222  CALL timeset(routinen, handle)
1223 
1224  CALL timeset(routinen//"_1", handle2)
1225 
1226  ! get info for the dest
1227  CALL cp_fm_get_info(matrix=fm_matrix_l, &
1228  nrow_local=nrow_local, &
1229  ncol_local=ncol_local, &
1230  row_indices=row_indices, &
1231  col_indices=col_indices, &
1232  nrow_block=nrow_block, &
1233  ncol_block=ncol_block)
1234 
1235  nprow = fm_matrix_l%matrix_struct%context%num_pe(1)
1236  npcol = fm_matrix_l%matrix_struct%context%num_pe(2)
1237 
1238  ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
1239  ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
1240 
1241  num_entries_rec(:) = 0
1242  num_entries_send(:) = 0
1243 
1244  dummy_proc = 0
1245 
1246  ! get the process, where the elements from L_local_col have to be sent
1247  DO lll = 1, dimen_ri
1248 
1249  send_prow = cp_fm_indxg2p(lll, nrow_block, dummy_proc, &
1250  fm_matrix_l%matrix_struct%first_p_pos(1), nprow)
1251 
1252  DO mmm = my_group_l_start, my_group_l_end
1253 
1254  send_pcol = cp_fm_indxg2p(mmm, ncol_block, dummy_proc, &
1255  fm_matrix_l%matrix_struct%first_p_pos(2), npcol)
1256 
1257  proc_send = fm_matrix_l%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
1258 
1259  num_entries_send(proc_send) = num_entries_send(proc_send) + 1
1260 
1261  END DO
1262 
1263  END DO
1264 
1265  CALL timestop(handle2)
1266 
1267  CALL timeset(routinen//"_2", handle2)
1268 
1269  CALL para_env%alltoall(num_entries_send, num_entries_rec, 1)
1270 
1271  CALL timestop(handle2)
1272 
1273  CALL timeset(routinen//"_3", handle2)
1274 
1275  ! allocate buffers to send the entries and the information of the entries
1276  ALLOCATE (buffer_rec(0:para_env%num_pe - 1))
1277  ALLOCATE (buffer_send(0:para_env%num_pe - 1))
1278 
1279  ! allocate data message and corresponding indices
1280  DO iproc = 0, para_env%num_pe - 1
1281 
1282  ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc)))
1283  buffer_rec(iproc)%msg = 0.0_dp
1284 
1285  END DO
1286 
1287  CALL timestop(handle2)
1288 
1289  CALL timeset(routinen//"_4", handle2)
1290 
1291  DO iproc = 0, para_env%num_pe - 1
1292 
1293  ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc)))
1294  buffer_send(iproc)%msg = 0.0_dp
1295 
1296  END DO
1297 
1298  CALL timestop(handle2)
1299 
1300  CALL timeset(routinen//"_5", handle2)
1301 
1302  DO iproc = 0, para_env%num_pe - 1
1303 
1304  ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2))
1305  buffer_rec(iproc)%indx = 0
1306 
1307  END DO
1308 
1309  CALL timestop(handle2)
1310 
1311  CALL timeset(routinen//"_6", handle2)
1312 
1313  DO iproc = 0, para_env%num_pe - 1
1314 
1315  ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2))
1316  buffer_send(iproc)%indx = 0
1317 
1318  END DO
1319 
1320  CALL timestop(handle2)
1321 
1322  CALL timeset(routinen//"_7", handle2)
1323 
1324  ALLOCATE (entry_counter(0:para_env%num_pe - 1))
1325  entry_counter(:) = 0
1326 
1327  ! get the process, where the elements from L_local_col have to be sent and
1328  ! write the data and the info to buffer_send
1329  DO lll = 1, dimen_ri
1330 
1331  send_prow = cp_fm_indxg2p(lll, nrow_block, dummy_proc, &
1332  fm_matrix_l%matrix_struct%first_p_pos(1), nprow)
1333 
1334  DO mmm = my_group_l_start, my_group_l_end
1335 
1336  send_pcol = cp_fm_indxg2p(mmm, ncol_block, dummy_proc, &
1337  fm_matrix_l%matrix_struct%first_p_pos(2), npcol)
1338 
1339  proc_send = fm_matrix_l%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
1340 
1341  entry_counter(proc_send) = entry_counter(proc_send) + 1
1342 
1343  buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
1344  l_local_col(lll, mmm - my_group_l_start + 1)
1345 
1346  buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = lll
1347  buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = mmm
1348 
1349  END DO
1350 
1351  END DO
1352 
1353  ALLOCATE (req_array(1:para_env%num_pe, 4))
1354 
1355  CALL timestop(handle2)
1356 
1357  CALL timeset(routinen//"_8", handle2)
1358 
1359  ! communicate the buffer
1360  CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, &
1361  buffer_send, req_array)
1362 
1363  fm_matrix_l%local_data = 0.0_dp
1364 
1365  CALL timestop(handle2)
1366 
1367  CALL timeset(routinen//"_9", handle2)
1368 
1369  ! fill fm_matrix_L with the entries from buffer_rec
1370  DO iproc = 0, para_env%num_pe - 1
1371 
1372  DO i_entry_rec = 1, num_entries_rec(iproc)
1373 
1374  DO i_row = 1, nrow_local
1375 
1376  i_row_global = row_indices(i_row)
1377 
1378  DO j_col = 1, ncol_local
1379 
1380  j_col_global = col_indices(j_col)
1381 
1382  IF (i_row_global == buffer_rec(iproc)%indx(i_entry_rec, 1) .AND. &
1383  j_col_global == buffer_rec(iproc)%indx(i_entry_rec, 2)) THEN
1384 
1385  fm_matrix_l%local_data(i_row, j_col) = buffer_rec(iproc)%msg(i_entry_rec)
1386 
1387  END IF
1388 
1389  END DO
1390 
1391  END DO
1392 
1393  END DO
1394 
1395  END DO
1396 
1397  CALL timestop(handle2)
1398 
1399  CALL timeset(routinen//"_10", handle2)
1400 
1401  DO iproc = 0, para_env%num_pe - 1
1402  DEALLOCATE (buffer_rec(iproc)%msg)
1403  DEALLOCATE (buffer_rec(iproc)%indx)
1404  DEALLOCATE (buffer_send(iproc)%msg)
1405  DEALLOCATE (buffer_send(iproc)%indx)
1406  END DO
1407 
1408  DEALLOCATE (buffer_rec, buffer_send)
1409  DEALLOCATE (req_array)
1410  DEALLOCATE (entry_counter)
1411  DEALLOCATE (num_entries_rec, num_entries_send)
1412 
1413  CALL timestop(handle2)
1414 
1415  CALL timestop(handle)
1416 
1417  END SUBROUTINE fill_fm_l_from_l_loc_non_blocking
1418 
1419 ! **************************************************************************************************
1420 !> \brief ...
1421 !> \param fm_matrix_L ...
1422 !> \param dimen_RI ...
1423 !> \param do_inversion ...
1424 ! **************************************************************************************************
1425  SUBROUTINE cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion)
1426 
1427  TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_l
1428  INTEGER, INTENT(IN) :: dimen_ri
1429  LOGICAL, INTENT(IN) :: do_inversion
1430 
1431  CHARACTER(LEN=*), PARAMETER :: routinen = 'cholesky_decomp'
1432 
1433  INTEGER :: handle, i_global, iib, info_chol, &
1434  j_global, jjb, ncol_local, nrow_local
1435  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1436 
1437  CALL timeset(routinen, handle)
1438 
1439  ! do cholesky decomposition
1440  CALL cp_fm_cholesky_decompose(matrix=fm_matrix_l, n=dimen_ri, info_out=info_chol)
1441  cpassert(info_chol == 0)
1442 
1443  IF (.NOT. do_inversion) THEN
1444  ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
1445  CALL cp_fm_get_info(matrix=fm_matrix_l, &
1446  nrow_local=nrow_local, &
1447  ncol_local=ncol_local, &
1448  row_indices=row_indices, &
1449  col_indices=col_indices)
1450  DO iib = 1, nrow_local
1451  i_global = row_indices(iib)
1452  DO jjb = 1, ncol_local
1453  j_global = col_indices(jjb)
1454  IF (j_global < i_global) fm_matrix_l%local_data(iib, jjb) = 0.0_dp
1455  END DO
1456  END DO
1457 
1458  END IF
1459 
1460  CALL timestop(handle)
1461 
1462  END SUBROUTINE cholesky_decomp
1463 
1464 ! **************************************************************************************************
1465 !> \brief ...
1466 !> \param fm_matrix_Minv_L_kpoints ...
1467 !> \param fm_matrix_L_kpoints ...
1468 !> \param dimen_RI ...
1469 !> \param kpoints ...
1470 !> \param eps_eigval_S ...
1471 ! **************************************************************************************************
1472  SUBROUTINE inversion_of_m_and_mult_with_chol_dec_of_v(fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, &
1473  dimen_RI, kpoints, eps_eigval_S)
1474 
1475  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: fm_matrix_minv_l_kpoints, &
1476  fm_matrix_l_kpoints
1477  INTEGER, INTENT(IN) :: dimen_ri
1478  TYPE(kpoint_type), POINTER :: kpoints
1479  REAL(kind=dp), INTENT(IN) :: eps_eigval_s
1480 
1481  CHARACTER(LEN=*), PARAMETER :: routinen = 'inversion_of_M_and_mult_with_chol_dec_of_V'
1482  COMPLEX(KIND=dp), PARAMETER :: cone = cmplx(1.0_dp, 0.0_dp, kind=dp), &
1483  czero = cmplx(0.0_dp, 0.0_dp, kind=dp), ione = cmplx(0.0_dp, 1.0_dp, kind=dp)
1484 
1485  INTEGER :: handle, ikp, nkp
1486  TYPE(cp_cfm_type) :: cfm_matrix_k_tmp, cfm_matrix_m_tmp, &
1487  cfm_matrix_v_tmp, cfm_matrix_vtrunc_tmp
1488  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1489 
1490  CALL timeset(routinen, handle)
1491 
1492  CALL cp_fm_get_info(fm_matrix_minv_l_kpoints(1, 1), matrix_struct=matrix_struct)
1493 
1494  CALL cp_cfm_create(cfm_matrix_m_tmp, matrix_struct)
1495  CALL cp_cfm_create(cfm_matrix_v_tmp, matrix_struct)
1496  CALL cp_cfm_create(cfm_matrix_k_tmp, matrix_struct)
1497  CALL cp_cfm_create(cfm_matrix_vtrunc_tmp, matrix_struct)
1498 
1499  CALL get_kpoint_info(kpoints, nkp=nkp)
1500 
1501  DO ikp = 1, nkp
1502 
1503  CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_m_tmp, cone, fm_matrix_minv_l_kpoints(ikp, 1))
1504  CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_m_tmp, ione, fm_matrix_minv_l_kpoints(ikp, 2))
1505 
1506  CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_v_tmp, cone, fm_matrix_l_kpoints(ikp, 1))
1507  CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_v_tmp, ione, fm_matrix_l_kpoints(ikp, 2))
1508 
1509  CALL cp_cfm_power(cfm_matrix_m_tmp, threshold=eps_eigval_s, exponent=-1.0_dp)
1510 
1511  CALL cp_cfm_power(cfm_matrix_v_tmp, threshold=0.0_dp, exponent=0.5_dp)
1512 
1513  ! save L(k) = SQRT(V(k))
1514  CALL cp_cfm_to_fm(cfm_matrix_v_tmp, fm_matrix_l_kpoints(ikp, 1), fm_matrix_l_kpoints(ikp, 2))
1515 
1516  ! get K = M^(-1)*L, where L is the Cholesky decomposition of V = L^T*L
1517  CALL parallel_gemm("N", "C", dimen_ri, dimen_ri, dimen_ri, cone, cfm_matrix_m_tmp, cfm_matrix_v_tmp, &
1518  czero, cfm_matrix_k_tmp)
1519 
1520  ! move
1521  CALL cp_cfm_to_fm(cfm_matrix_k_tmp, fm_matrix_minv_l_kpoints(ikp, 1), fm_matrix_minv_l_kpoints(ikp, 2))
1522 
1523  END DO
1524 
1525  CALL cp_cfm_release(cfm_matrix_m_tmp)
1526  CALL cp_cfm_release(cfm_matrix_v_tmp)
1527  CALL cp_cfm_release(cfm_matrix_k_tmp)
1528  CALL cp_cfm_release(cfm_matrix_vtrunc_tmp)
1529 
1530  CALL timestop(handle)
1531 
1532  END SUBROUTINE inversion_of_m_and_mult_with_chol_dec_of_v
1533 
1534 ! **************************************************************************************************
1535 !> \brief ...
1536 !> \param fm_matrix_Minv_Vtrunc_Minv ...
1537 !> \param fm_matrix_Minv ...
1538 !> \param qs_env ...
1539 ! **************************************************************************************************
1540  SUBROUTINE gamma_only_inversion_of_m_and_mult_with_chol_dec_of_vtrunc(fm_matrix_Minv_Vtrunc_Minv, &
1541  fm_matrix_Minv, qs_env)
1542 
1543  TYPE(cp_fm_type), DIMENSION(:, :) :: fm_matrix_minv_vtrunc_minv, &
1544  fm_matrix_minv
1545  TYPE(qs_environment_type), POINTER :: qs_env
1546 
1547  CHARACTER(LEN=*), PARAMETER :: &
1548  routinen = 'Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc'
1549 
1550  INTEGER :: dimen_ri, handle, ndep
1551  REAL(kind=dp) :: eps_eigval_s_gamma
1552  TYPE(cp_fm_type) :: fm_matrix_ri_metric_inv_work, fm_work
1553 
1554  CALL timeset(routinen, handle)
1555 
1556  CALL cp_fm_create(fm_work, fm_matrix_minv(1, 1)%matrix_struct)
1557  CALL cp_fm_set_all(fm_work, 0.0_dp)
1558 
1559  CALL cp_fm_create(fm_matrix_ri_metric_inv_work, fm_matrix_minv(1, 1)%matrix_struct)
1560  CALL cp_fm_set_all(fm_matrix_ri_metric_inv_work, 0.0_dp)
1561 
1562  CALL cp_fm_get_info(matrix=fm_matrix_minv(1, 1), nrow_global=dimen_ri)
1563 
1564  eps_eigval_s_gamma = qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S_Gamma
1565 
1566  IF (eps_eigval_s_gamma > 1.0e-18) THEN
1567 
1568  ! remove small eigenvalues
1569  CALL cp_fm_power(fm_matrix_minv(1, 1), fm_matrix_ri_metric_inv_work, -0.5_dp, &
1570  eps_eigval_s_gamma, ndep)
1571 
1572  ELSE
1573 
1574  CALL cholesky_decomp(fm_matrix_minv(1, 1), dimen_ri, do_inversion=.true.)
1575 
1576  CALL invert_mat(fm_matrix_minv(1, 1))
1577 
1578  END IF
1579 
1580  CALL parallel_gemm('N', 'T', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_matrix_minv(1, 1), &
1581  fm_matrix_minv(1, 1), 0.0_dp, fm_matrix_ri_metric_inv_work)
1582 
1583  CALL cp_fm_to_fm(fm_matrix_ri_metric_inv_work, fm_matrix_minv(1, 1))
1584 
1585  CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_matrix_ri_metric_inv_work, &
1586  fm_matrix_minv_vtrunc_minv(1, 1), 0.0_dp, fm_work)
1587 
1588  CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_work, &
1589  fm_matrix_ri_metric_inv_work, 0.0_dp, fm_matrix_minv_vtrunc_minv(1, 1))
1590 
1591  CALL cp_fm_release(fm_work)
1592  CALL cp_fm_release(fm_matrix_ri_metric_inv_work)
1593 
1594  CALL timestop(handle)
1595 
1596  END SUBROUTINE gamma_only_inversion_of_m_and_mult_with_chol_dec_of_vtrunc
1597 
1598 ! **************************************************************************************************
1599 !> \brief ...
1600 !> \param qs_env ...
1601 !> \param trunc_coulomb ...
1602 !> \param rel_cutoff_trunc_coulomb_ri_x ...
1603 ! **************************************************************************************************
1604  SUBROUTINE setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env, trunc_coulomb, &
1605  rel_cutoff_trunc_coulomb_ri_x)
1606  TYPE(qs_environment_type), POINTER :: qs_env
1607  TYPE(libint_potential_type), OPTIONAL :: trunc_coulomb
1608  REAL(kind=dp), OPTIONAL :: rel_cutoff_trunc_coulomb_ri_x
1609 
1610  CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_trunc_coulomb_pot_for_exchange_self_energy'
1611 
1612  INTEGER :: handle
1613  INTEGER, DIMENSION(3) :: periodic
1614  REAL(kind=dp) :: my_rel_cutoff_trunc_coulomb_ri_x, shortest_dist_cell_planes
1615  TYPE(cell_type), POINTER :: cell
1616 
1617  CALL timeset(routinen, handle)
1618 
1619  NULLIFY (cell)
1620  CALL get_qs_env(qs_env, cell=cell)
1621 
1622  CALL get_cell(cell=cell, periodic=periodic)
1623 
1624  shortest_dist_cell_planes = 1.0e4_dp
1625  IF (periodic(1) == 1) THEN
1626  IF (shortest_dist_cell_planes > plane_distance(1, 0, 0, cell)) THEN
1627  shortest_dist_cell_planes = plane_distance(1, 0, 0, cell)
1628  END IF
1629  END IF
1630  IF (periodic(2) == 1) THEN
1631  IF (shortest_dist_cell_planes > plane_distance(0, 1, 0, cell)) THEN
1632  shortest_dist_cell_planes = plane_distance(0, 1, 0, cell)
1633  END IF
1634  END IF
1635  IF (periodic(3) == 1) THEN
1636  IF (shortest_dist_cell_planes > plane_distance(0, 0, 1, cell)) THEN
1637  shortest_dist_cell_planes = plane_distance(0, 0, 1, cell)
1638  END IF
1639  END IF
1640 
1641  IF (PRESENT(rel_cutoff_trunc_coulomb_ri_x)) THEN
1642  my_rel_cutoff_trunc_coulomb_ri_x = rel_cutoff_trunc_coulomb_ri_x
1643  ELSE
1644  my_rel_cutoff_trunc_coulomb_ri_x = qs_env%mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x
1645  END IF
1646 
1647  IF (PRESENT(trunc_coulomb)) THEN
1648  trunc_coulomb%potential_type = do_potential_truncated
1649  trunc_coulomb%cutoff_radius = shortest_dist_cell_planes*my_rel_cutoff_trunc_coulomb_ri_x
1650  trunc_coulomb%filename = "t_c_g.dat"
1651  ! dummy
1652  trunc_coulomb%omega = 0.0_dp
1653  END IF
1654 
1655  CALL timestop(handle)
1656 
1657  END SUBROUTINE setup_trunc_coulomb_pot_for_exchange_self_energy
1658 
1659 ! **************************************************************************************************
1660 !> \brief ...
1661 !> \param fm_matrix_L ...
1662 ! **************************************************************************************************
1663  SUBROUTINE invert_mat(fm_matrix_L)
1664 
1665  TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_l
1666 
1667  CHARACTER(LEN=*), PARAMETER :: routinen = 'invert_mat'
1668 
1669  INTEGER :: handle, i_global, iib, j_global, jjb, &
1670  ncol_local, nrow_local
1671  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1672 
1673  CALL timeset(routinen, handle)
1674 
1675  CALL cp_fm_triangular_invert(matrix_a=fm_matrix_l, uplo_tr='U')
1676 
1677  ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
1678  CALL cp_fm_get_info(matrix=fm_matrix_l, &
1679  nrow_local=nrow_local, &
1680  ncol_local=ncol_local, &
1681  row_indices=row_indices, &
1682  col_indices=col_indices)
1683  DO iib = 1, nrow_local
1684  i_global = row_indices(iib)
1685  DO jjb = 1, ncol_local
1686  j_global = col_indices(jjb)
1687  IF (j_global < i_global) fm_matrix_l%local_data(iib, jjb) = 0.0_dp
1688  END DO
1689  END DO
1690 
1691  CALL timestop(handle)
1692 
1693  END SUBROUTINE invert_mat
1694 
1695 ! **************************************************************************************************
1696 !> \brief ...
1697 !> \param fm_matrix_L ...
1698 !> \param blacs_env_L ...
1699 !> \param dimen_RI ...
1700 !> \param para_env_L ...
1701 !> \param name ...
1702 !> \param fm_matrix_L_extern ...
1703 ! **************************************************************************************************
1704  SUBROUTINE create_matrix_l(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, name, fm_matrix_L_extern)
1705  TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_l
1706  TYPE(cp_blacs_env_type), POINTER :: blacs_env_l
1707  INTEGER, INTENT(IN) :: dimen_ri
1708  TYPE(mp_para_env_type), POINTER :: para_env_l
1709  CHARACTER(LEN=*), INTENT(IN) :: name
1710  TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_matrix_l_extern
1711 
1712  CHARACTER(LEN=*), PARAMETER :: routinen = 'create_matrix_L'
1713 
1714  INTEGER :: handle
1715  TYPE(cp_fm_struct_type), POINTER :: fm_struct
1716 
1717  CALL timeset(routinen, handle)
1718 
1719  ! create the full matrix L defined in the L group
1720  IF (PRESENT(fm_matrix_l_extern)) THEN
1721  CALL cp_fm_create(fm_matrix_l, fm_matrix_l_extern%matrix_struct, name=name)
1722  ELSE
1723  NULLIFY (fm_struct)
1724  CALL cp_fm_struct_create(fm_struct, context=blacs_env_l, nrow_global=dimen_ri, &
1725  ncol_global=dimen_ri, para_env=para_env_l)
1726 
1727  CALL cp_fm_create(fm_matrix_l, fm_struct, name=name)
1728 
1729  CALL cp_fm_struct_release(fm_struct)
1730  END IF
1731 
1732  CALL cp_fm_set_all(matrix=fm_matrix_l, alpha=0.0_dp)
1733 
1734  CALL timestop(handle)
1735 
1736  END SUBROUTINE create_matrix_l
1737 
1738 ! **************************************************************************************************
1739 !> \brief ...
1740 !> \param fm_matrix_L ...
1741 !> \param my_group_L_start ...
1742 !> \param my_group_L_end ...
1743 !> \param my_group_L_size ...
1744 !> \param my_Lrows ...
1745 ! **************************************************************************************************
1746  SUBROUTINE grep_lcols(fm_matrix_L, &
1747  my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
1748  TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_l
1749  INTEGER, INTENT(IN) :: my_group_l_start, my_group_l_end, &
1750  my_group_l_size
1751  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
1752  INTENT(OUT) :: my_lrows
1753 
1754  CHARACTER(LEN=*), PARAMETER :: routinen = 'grep_Lcols'
1755 
1756  INTEGER :: dimen_ri, handle, handle2, i_global, iib, j_global, jjb, max_row_col_local, &
1757  ncol_local, ncol_rec, nrow_local, nrow_rec, proc_receive_static, proc_send_static, &
1758  proc_shift
1759  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, rec_col_row_info
1760  INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_rec, &
1761  row_indices, row_indices_rec
1762  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: local_l, rec_l
1763  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1764  POINTER :: local_l_internal
1765  TYPE(mp_para_env_type), POINTER :: para_env
1766 
1767  CALL timeset(routinen, handle)
1768 
1769  CALL cp_fm_get_info(matrix=fm_matrix_l, &
1770  nrow_local=nrow_local, &
1771  ncol_local=ncol_local, &
1772  row_indices=row_indices, &
1773  col_indices=col_indices, &
1774  nrow_global=dimen_ri, &
1775  local_data=local_l_internal, &
1776  para_env=para_env)
1777 
1778  ALLOCATE (my_lrows(dimen_ri, my_group_l_size))
1779  my_lrows = 0.0_dp
1780 
1781  ALLOCATE (local_l(nrow_local, ncol_local))
1782  local_l(:, :) = local_l_internal(1:nrow_local, 1:ncol_local)
1783 
1784  max_row_col_local = max(nrow_local, ncol_local)
1785  CALL para_env%max(max_row_col_local)
1786 
1787  ALLOCATE (local_col_row_info(0:max_row_col_local, 2))
1788  local_col_row_info = 0
1789  ! 0,1 nrows
1790  local_col_row_info(0, 1) = nrow_local
1791  local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local)
1792  ! 0,2 ncols
1793  local_col_row_info(0, 2) = ncol_local
1794  local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local)
1795 
1796  ALLOCATE (rec_col_row_info(0:max_row_col_local, 2))
1797 
1798  ! accumulate data on my_Lrows starting from myself
1799  DO jjb = 1, ncol_local
1800  j_global = col_indices(jjb)
1801  IF (j_global >= my_group_l_start .AND. j_global <= my_group_l_end) THEN
1802  DO iib = 1, nrow_local
1803  i_global = row_indices(iib)
1804  my_lrows(i_global, j_global - my_group_l_start + 1) = local_l(iib, jjb)
1805  END DO
1806  END IF
1807  END DO
1808 
1809  proc_send_static = modulo(para_env%mepos + 1, para_env%num_pe)
1810  proc_receive_static = modulo(para_env%mepos - 1, para_env%num_pe)
1811 
1812  CALL timeset(routinen//"_comm", handle2)
1813 
1814  DO proc_shift = 1, para_env%num_pe - 1
1815  ! first exchange information on the local data
1816  rec_col_row_info = 0
1817  CALL para_env%sendrecv(local_col_row_info, proc_send_static, rec_col_row_info, proc_receive_static)
1818  nrow_rec = rec_col_row_info(0, 1)
1819  ncol_rec = rec_col_row_info(0, 2)
1820 
1821  ALLOCATE (row_indices_rec(nrow_rec))
1822  row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
1823 
1824  ALLOCATE (col_indices_rec(ncol_rec))
1825  col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
1826 
1827  ALLOCATE (rec_l(nrow_rec, ncol_rec))
1828  rec_l = 0.0_dp
1829 
1830  ! then send and receive the real data
1831  CALL para_env%sendrecv(local_l, proc_send_static, rec_l, proc_receive_static)
1832 
1833  ! accumulate the received data on my_Lrows
1834  DO jjb = 1, ncol_rec
1835  j_global = col_indices_rec(jjb)
1836  IF (j_global >= my_group_l_start .AND. j_global <= my_group_l_end) THEN
1837  DO iib = 1, nrow_rec
1838  i_global = row_indices_rec(iib)
1839  my_lrows(i_global, j_global - my_group_l_start + 1) = rec_l(iib, jjb)
1840  END DO
1841  END IF
1842  END DO
1843 
1844  local_col_row_info(:, :) = rec_col_row_info
1845  CALL move_alloc(rec_l, local_l)
1846 
1847  DEALLOCATE (col_indices_rec)
1848  DEALLOCATE (row_indices_rec)
1849  END DO
1850  CALL timestop(handle2)
1851 
1852  DEALLOCATE (local_col_row_info)
1853  DEALLOCATE (rec_col_row_info)
1854  DEALLOCATE (local_l)
1855 
1856  CALL timestop(handle)
1857 
1858  END SUBROUTINE grep_lcols
1859 
1860 END MODULE mp2_ri_2c
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
basis_set_types
Definition: basis_set_types.F:15
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition: cell_types.F:195
cell_types::plane_distance
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition: cell_types.F:252
constants_operator
constants for the different operators of the 2c-integrals
Definition: constants_operator.F:14
constants_operator::operator_coulomb
integer, parameter, public operator_coulomb
Definition: constants_operator.F:22
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_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition: cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add_fm
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Definition: cp_cfm_basic_linalg.F:355
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition: cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition: cp_cfm_types.F:159
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition: cp_cfm_types.F:765
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_dist2d_to_dist
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
Definition: cp_dbcsr_operations.F:611
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_eri_mme_interface
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
Definition: cp_eri_mme_interface.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_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1882
cp_fm_basic_linalg::cp_fm_triangular_invert
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
Definition: cp_fm_basic_linalg.F:2193
cp_fm_cholesky
various cholesky decomposition related routines
Definition: cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition: cp_fm_cholesky.F:50
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition: cp_fm_diag.F:896
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition: cp_fm_diag.F:208
cp_fm_diag::cp_fm_syevx
subroutine, public cp_fm_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx)
compute eigenvalues and optionally eigenvectors of a real symmetric matrix using scalapack....
Definition: cp_fm_diag.F:657
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition: cp_fm_types.F:1538
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_indxg2p
integer function, public cp_fm_indxg2p(INDXGLOB, NB, IPROC, ISRCPROC, NPROCS)
wrapper to scalapack function INDXG2P that computes the process coordinate which possesses the entry ...
Definition: cp_fm_types.F:2466
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
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition: distribution_2d_types.F:15
group_dist_types
Types to describe group distributions.
Definition: group_dist_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_eri_os
integer, parameter, public do_eri_os
Definition: input_constants.F:1226
input_constants::do_eri_mme
integer, parameter, public do_eri_mme
Definition: input_constants.F:1226
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition: input_constants.F:801
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition: input_constants.F:801
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition: input_constants.F:801
input_constants::do_potential_long
integer, parameter, public do_potential_long
Definition: input_constants.F:801
input_constants::do_eri_gpw
integer, parameter, public do_eri_gpw
Definition: input_constants.F:1226
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kpoint_coulomb_2c
Routines to compute the Coulomb integral V_(alpha beta)(k) for a k-point k using lattice summation in...
Definition: kpoint_coulomb_2c.F:15
kpoint_coulomb_2c::build_2c_coulomb_matrix_kp
subroutine, public build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, atomic_kind_set, size_lattice_sum, operator_type, ikp_start, ikp_end)
...
Definition: kpoint_coulomb_2c.F:69
kpoint_methods
Routines needed for kpoint calculation.
Definition: kpoint_methods.F:15
kpoint_methods::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition: kpoint_methods.F:780
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition: libint_2c_3c.F:14
libint_2c_3c::compare_potential_types
pure logical function, public compare_potential_types(potential1, potential2)
Helper function to compare libint_potential_types.
Definition: libint_2c_3c.F:963
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
message_passing::mp_para_env_release
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Definition: message_passing.F:2027
mp2_eri_gpw
Routines to calculate 2c- and 3c-integrals for RI with GPW.
Definition: mp2_eri_gpw.F:13
mp2_eri_gpw::mp2_eri_2c_integrate_gpw
subroutine, public mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_L_start, my_group_L_end, natom, potential_parameter, sab_orb_sub, L_local_col, kind_of)
Integrates the potential of an RI function.
Definition: mp2_eri_gpw.F:170
mp2_eri
Interface to direct methods for electron repulsion integrals for MP2.
Definition: mp2_eri.F:12
mp2_eri::mp2_eri_2c_integrate
subroutine, public mp2_eri_2c_integrate(param, potential_parameter, para_env, qs_env, basis_type_a, basis_type_b, hab, first_b, last_b, eri_method, pab, force_a, force_b, hdab, hadb, reflection_z_a, reflection_z_b, do_reflection_a, do_reflection_b)
high-level integration routine for 2c integrals over CP2K basis sets. Contiguous column-wise distribu...
Definition: mp2_eri.F:115
mp2_ri_2c
Framework for 2c-integrals for RI.
Definition: mp2_ri_2c.F:14
mp2_ri_2c::setup_trunc_coulomb_pot_for_exchange_self_energy
subroutine, public setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env, trunc_coulomb, rel_cutoff_trunc_coulomb_ri_x)
...
Definition: mp2_ri_2c.F:1606
mp2_ri_2c::inversion_of_m_and_mult_with_chol_dec_of_v
subroutine, public inversion_of_m_and_mult_with_chol_dec_of_v(fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, dimen_RI, kpoints, eps_eigval_S)
...
Definition: mp2_ri_2c.F:1474
mp2_ri_2c::get_2c_integrals
subroutine, public get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, mp2_memory, my_Lrows, my_Vrows, fm_matrix_PQ, ngroup, color_sub, dimen_RI, dimen_RI_red, kpoints, my_group_L_size, my_group_L_start, my_group_L_end, gd_array, calc_PQ_cond_num, do_svd, potential, ri_metric, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, do_im_time, do_kpoints, mp2_eps_pgf_orb_S, qs_kind_set, sab_orb_sub, calc_forces, unit_nr)
...
Definition: mp2_ri_2c.F:162
mp2_ri_2c::ri_2c_integral_mat
subroutine, public ri_2c_integral_mat(qs_env, fm_matrix_Minv_L_kpoints, fm_matrix_L, dimen_RI, ri_metric, do_kpoints, kpoints, put_mat_KS_env, regularization_RI, ikp_ext)
...
Definition: mp2_ri_2c.F:554
mp2_types
Types needed for MP2 calculations.
Definition: mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition: particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition: particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
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
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition: qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition: qs_integral_utils.F:149
qs_interactions
Calculate the interaction radii for the operator matrix calculation.
Definition: qs_interactions.F:17
qs_interactions::init_interaction_radii
subroutine, public init_interaction_radii(qs_control, qs_kind_set)
Initialize all the atomic kind radii for a given threshold value.
Definition: qs_interactions.F:75
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition: qs_ks_types.F:330
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition: qs_ks_types.F:567
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition: qs_neighbor_list_types.F:1097
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition: qs_tensors.F:11
qs_tensors::build_2c_neighbor_lists
subroutine, public build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, sym_ij, molecular, dist_2d, pot_to_rad)
Build 2-center neighborlists adapted to different operators This mainly wraps build_neighbor_lists fo...
Definition: qs_tensors.F:143
qs_tensors::build_2c_integrals
subroutine, public build_2c_integrals(t2c, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints, do_hfx_kpoints, ext_kpoints, regularization_RI)
...
Definition: qs_tensors.F:3406
rpa_communication
Auxiliary routines necessary to redistribute an fm_matrix from a given blacs_env to another.
Definition: rpa_communication.F:14
rpa_communication::communicate_buffer
subroutine, public communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array, do_indx, do_msg)
...
Definition: rpa_communication.F:457
rpa_gw_kpoints_util
Routines treating GW and RPA calculations with kpoints.
Definition: rpa_gw_kpoints_util.F:13
rpa_gw_kpoints_util::cp_cfm_power
subroutine, public cp_cfm_power(matrix, threshold, exponent, min_eigval)
...
Definition: rpa_gw_kpoints_util.F:277