(git:6344619)
Loading...
Searching...
No Matches
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-2026 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! **************************************************************************************************
14MODULE 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_diag, ONLY: cp_cfm_power
28 USE cp_cfm_types, ONLY: cp_cfm_copy_general,&
29 cp_cfm_create,&
30 cp_cfm_release,&
31 cp_cfm_to_cfm,&
32 cp_cfm_to_fm,&
33 cp_cfm_type,&
34 cp_fm_to_cfm
35 USE cp_control_types, ONLY: dft_control_type
36 USE cp_dbcsr_api, ONLY: &
37 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
38 dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_p_type, dbcsr_release, &
39 dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
40 dbcsr_type_symmetric
41 USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
42 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
43 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
44 cp_dbcsr_dist2d_to_dist,&
45 dbcsr_allocate_matrix_set,&
46 dbcsr_deallocate_matrix_set
47 USE cp_eri_mme_interface, ONLY: cp_eri_mme_param
48 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
49 cp_fm_triangular_invert
50 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
51 USE cp_fm_diag, ONLY: choose_eigv_solver,&
52 cp_fm_power,&
53 cp_fm_syevx
54 USE cp_fm_redistribution, ONLY: cp_array_columns_to_fm
55 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
56 cp_fm_struct_release,&
57 cp_fm_struct_type
58 USE cp_fm_types, ONLY: cp_fm_copy_general,&
59 cp_fm_create,&
60 cp_fm_get_info,&
61 cp_fm_release,&
62 cp_fm_set_all,&
63 cp_fm_to_fm,&
64 cp_fm_type
65 USE distribution_2d_types, ONLY: distribution_2d_type
66 USE group_dist_types, ONLY: create_group_dist,&
67 get_group_dist,&
68 group_dist_d1_type,&
69 release_group_dist
70 USE input_constants, ONLY: do_eri_gpw,&
71 do_eri_mme,&
72 do_eri_os,&
73 do_potential_coulomb,&
74 do_potential_id,&
75 do_potential_long,&
76 do_potential_truncated
77 USE kinds, ONLY: dp
78 USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
79 USE kpoint_methods, ONLY: kpoint_init_cell_index,&
80 rskp_transform
81 USE kpoint_types, ONLY: get_kpoint_info,&
82 kpoint_type
83 USE libint_2c_3c, ONLY: compare_potential_types,&
84 libint_potential_type
85 USE machine, ONLY: m_flush
86 USE mathconstants, ONLY: gaussi,&
87 z_one,&
88 z_zero
89 USE message_passing, ONLY: mp_comm_type,&
90 mp_para_env_release,&
91 mp_para_env_type
92 USE mp2_eri, ONLY: mp2_eri_2c_integrate
93 USE mp2_eri_gpw, ONLY: mp2_eri_2c_integrate_gpw
94 USE parallel_gemm_api, ONLY: parallel_gemm
95 USE particle_methods, ONLY: get_particle_set
96 USE particle_types, ONLY: particle_type
97 USE qs_environment_types, ONLY: get_qs_env,&
98 qs_environment_type
99 USE qs_integral_utils, ONLY: basis_set_list_setup
100 USE qs_interactions, ONLY: init_interaction_radii
101 USE qs_kind_types, ONLY: get_qs_kind,&
102 qs_kind_type
103 USE qs_ks_types, ONLY: get_ks_env,&
104 set_ks_env
105 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
106 release_neighbor_list_sets
107 USE qs_tensors, ONLY: build_2c_integrals,&
108 build_2c_neighbor_lists
109
110!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
111#include "./base/base_uses.f90"
112
113 IMPLICIT NONE
114
115 PRIVATE
116
117 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_2c'
118
119 PUBLIC :: get_2c_integrals, trunc_coulomb_for_exchange, ri_2c_integral_mat, ri_2c_integral_mat_kp, &
120 inversion_of_m_and_mult_with_chol_dec_of_v
121
122CONTAINS
123
124! **************************************************************************************************
125!> \brief ...
126!> \param qs_env ...
127!> \param eri_method ...
128!> \param eri_param ...
129!> \param para_env ...
130!> \param para_env_sub ...
131!> \param mp2_memory ...
132!> \param my_Lrows ...
133!> \param my_Vrows ...
134!> \param fm_matrix_PQ ...
135!> \param ngroup ...
136!> \param color_sub ...
137!> \param dimen_RI ...
138!> \param dimen_RI_red reduced RI dimension, needed if we perform SVD instead of Cholesky
139!> \param kpoints ...
140!> \param my_group_L_size ...
141!> \param my_group_L_start ...
142!> \param my_group_L_end ...
143!> \param gd_array ...
144!> \param calc_PQ_cond_num ...
145!> \param do_svd ...
146!> \param eps_svd ...
147!> \param potential ...
148!> \param ri_metric ...
149!> \param cfm_matrix_L_kpoints ...
150!> \param cfm_matrix_Minv_L_kpoints ...
151!> \param fm_matrix_Minv_L_gamma ...
152!> \param fm_matrix_Minv ...
153!> \param fm_matrix_Minv_Vtrunc_Minv ...
154!> \param do_im_time ...
155!> \param do_kpoints ...
156!> \param mp2_eps_pgf_orb_S ...
157!> \param qs_kind_set ...
158!> \param sab_orb_sub ...
159!> \param calc_forces ...
160!> \param unit_nr ...
161! **************************************************************************************************
162 SUBROUTINE get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, mp2_memory, &
163 my_Lrows, my_Vrows, fm_matrix_PQ, ngroup, color_sub, dimen_RI, dimen_RI_red, &
164 kpoints, my_group_L_size, my_group_L_start, my_group_L_end, &
165 gd_array, calc_PQ_cond_num, do_svd, eps_svd, potential, ri_metric, &
166 cfm_matrix_L_kpoints, cfm_matrix_Minv_L_kpoints, fm_matrix_Minv_L_gamma, &
167 fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
168 do_im_time, do_kpoints, mp2_eps_pgf_orb_S, qs_kind_set, sab_orb_sub, calc_forces, unit_nr)
169
170 TYPE(qs_environment_type), POINTER :: qs_env
171 INTEGER, INTENT(IN) :: eri_method
172 TYPE(cp_eri_mme_param), POINTER :: eri_param
173 TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
174 REAL(kind=dp), INTENT(IN) :: mp2_memory
175 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
176 INTENT(OUT) :: my_lrows, my_vrows
177 TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_pq
178 INTEGER, INTENT(IN) :: ngroup, color_sub
179 INTEGER, INTENT(OUT) :: dimen_ri, dimen_ri_red
180 TYPE(kpoint_type), POINTER :: kpoints
181 INTEGER, INTENT(OUT) :: my_group_l_size, my_group_l_start, &
182 my_group_l_end
183 TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array
184 LOGICAL, INTENT(IN) :: calc_pq_cond_num, do_svd
185 REAL(kind=dp), INTENT(IN) :: eps_svd
186 TYPE(libint_potential_type) :: potential, ri_metric
187 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_matrix_l_kpoints, &
188 cfm_matrix_minv_l_kpoints
189 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :), &
190 INTENT(OUT) :: fm_matrix_minv_l_gamma
191 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_minv, &
192 fm_matrix_minv_vtrunc_minv
193 LOGICAL, INTENT(IN) :: do_im_time, do_kpoints
194 REAL(kind=dp), INTENT(IN) :: mp2_eps_pgf_orb_s
195 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
196 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
197 POINTER :: sab_orb_sub
198 LOGICAL, INTENT(IN) :: calc_forces
199 INTEGER, INTENT(IN) :: unit_nr
200
201 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_2c_integrals'
202
203 INTEGER :: handle, i_size, num_small_eigen
204 REAL(kind=dp) :: cond_num, eps_pgf_orb_old
205 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_matrix_m_kpoints
206 TYPE(cp_fm_type) :: fm_matrix_l_work, fm_matrix_m_inv_work, &
207 fm_matrix_v
208 TYPE(dft_control_type), POINTER :: dft_control
209 TYPE(libint_potential_type) :: trunc_coulomb
210 TYPE(mp_para_env_type), POINTER :: para_env_l
211
212 CALL timeset(routinen, handle)
213
214 ! calculate V and store it in fm_matrix_L_work
215 CALL compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_l, mp2_memory, &
216 fm_matrix_l_work, ngroup, color_sub, dimen_ri, &
217 my_group_l_size, my_group_l_start, my_group_l_end, &
218 gd_array, calc_pq_cond_num, cond_num, &
219 num_small_eigen, potential, sab_orb_sub, do_im_time=do_im_time)
220
221 IF (do_im_time .AND. calc_forces) THEN
222 !need a copy of the (P|Q) integrals
223 CALL cp_fm_create(fm_matrix_pq, fm_matrix_l_work%matrix_struct)
224 CALL cp_fm_to_fm(fm_matrix_l_work, fm_matrix_pq)
225 END IF
226
227 dimen_ri_red = dimen_ri
228
229 IF (compare_potential_types(ri_metric, potential) .AND. .NOT. do_im_time) THEN
230 CALL decomp_mat_l(fm_matrix_l_work, do_svd, eps_svd, num_small_eigen, cond_num, .true., gd_array, ngroup, &
231 dimen_ri, dimen_ri_red, para_env)
232 ELSE
233
234 ! RI-metric matrix (depending on RI metric: Coulomb, overlap or truncated Coulomb)
235 IF (do_im_time) THEN
236 CALL get_qs_env(qs_env, dft_control=dft_control)
237
238 ! re-init the radii to be able to generate pair lists with appropriate screening for overlap matrix
239 eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
240 dft_control%qs_control%eps_pgf_orb = mp2_eps_pgf_orb_s
241 CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
242
243 IF (do_kpoints) THEN
244 CALL ri_2c_integral_mat_kp(qs_env, cfm_matrix_m_kpoints, fm_matrix_l_work%matrix_struct, dimen_ri, &
245 ri_metric, kpoints, put_mat_ks_env=.true., &
246 regularization_ri=qs_env%mp2_env%ri_rpa_im_time%regularization_RI)
247 ELSE
248 CALL ri_2c_integral_mat(qs_env, fm_matrix_minv_l_gamma, fm_matrix_l_work%matrix_struct, dimen_ri, &
249 ri_metric, put_mat_ks_env=.true., &
250 regularization_ri=qs_env%mp2_env%ri_rpa_im_time%regularization_RI)
251 END IF
252
253 ! re-init the radii to previous values
254 dft_control%qs_control%eps_pgf_orb = eps_pgf_orb_old
255 CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
256
257 ! GPW overlap matrix
258 ELSE
259
260 CALL mp_para_env_release(para_env_l)
261 CALL release_group_dist(gd_array)
262
263 ALLOCATE (fm_matrix_minv_l_gamma(1, 1))
264
265 ! Calculate matrix of RI operator (for overlap metric: S), store it in fm_matrix_Minv_L_gamma
266 CALL compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_l, mp2_memory, &
267 fm_matrix_minv_l_gamma(1, 1), ngroup, color_sub, dimen_ri, &
268 my_group_l_size, my_group_l_start, my_group_l_end, &
269 gd_array, calc_pq_cond_num, cond_num, &
270 num_small_eigen, ri_metric, sab_orb_sub, &
271 fm_matrix_l_struct=fm_matrix_l_work%matrix_struct)
272
273 END IF
274
275 IF (do_kpoints) THEN
276
277 ! k-dependent 1/r Coulomb matrix V_PQ(k)
278 CALL compute_v_by_lattice_sum(qs_env, cfm_matrix_l_kpoints, fm_matrix_l_work%matrix_struct, kpoints)
279
280 CALL inversion_of_m_and_mult_with_chol_dec_of_v(cfm_matrix_m_kpoints, cfm_matrix_l_kpoints, &
281 cfm_matrix_minv_l_kpoints, dimen_ri, kpoints, &
282 qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S)
283 DO i_size = 1, SIZE(cfm_matrix_m_kpoints)
284 CALL cp_cfm_release(cfm_matrix_m_kpoints(i_size))
285 END DO
286 DEALLOCATE (cfm_matrix_m_kpoints)
287
288 CALL trunc_coulomb_for_exchange(qs_env, trunc_coulomb)
289
290 ! Gamma-only truncated Coulomb matrix V^tr with cutoff radius = half the unit cell size; for exchange self-energy
291 CALL ri_2c_integral_mat(qs_env, fm_matrix_minv_vtrunc_minv, fm_matrix_l_work%matrix_struct, dimen_ri, trunc_coulomb, &
292 put_mat_ks_env=.false., regularization_ri=0.0_dp)
293
294 ! Gamma-only RI-metric matrix; for computing Gamma-only/MIC self-energy
295 CALL ri_2c_integral_mat(qs_env, fm_matrix_minv, fm_matrix_l_work%matrix_struct, dimen_ri, ri_metric, &
296 put_mat_ks_env=.false., regularization_ri=0.0_dp)
297
298 CALL gamma_only_inversion_of_m_and_mult_with_chol_dec_of_vtrunc(fm_matrix_minv_vtrunc_minv, &
299 fm_matrix_minv, qs_env)
300 ELSE
301 IF (calc_forces .AND. (.NOT. do_im_time)) THEN
302 ! For the gradients, we make a copy of the inverse root of L
303 CALL cp_fm_create(fm_matrix_v, fm_matrix_l_work%matrix_struct)
304 CALL cp_fm_to_fm(fm_matrix_l_work, fm_matrix_v)
305
306 CALL decomp_mat_l(fm_matrix_v, do_svd, eps_svd, num_small_eigen, cond_num, .true., gd_array, ngroup, &
307 dimen_ri, dimen_ri_red, para_env)
308 END IF
309
310 CALL decomp_mat_l(fm_matrix_l_work, do_svd, eps_svd, num_small_eigen, cond_num, .false., gd_array, ngroup, &
311 dimen_ri, dimen_ri_red, para_env)
312
313 CALL decomp_mat_l(fm_matrix_minv_l_gamma(1, 1), .false., 0.0_dp, num_small_eigen, cond_num, .true., &
314 gd_array, ngroup, dimen_ri, dimen_ri_red, para_env)
315
316 CALL cp_fm_create(fm_matrix_m_inv_work, fm_matrix_minv_l_gamma(1, 1)%matrix_struct)
317 CALL cp_fm_set_all(fm_matrix_m_inv_work, 0.0_dp)
318
319 CALL parallel_gemm('N', 'T', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_matrix_minv_l_gamma(1, 1), &
320 fm_matrix_minv_l_gamma(1, 1), 0.0_dp, fm_matrix_m_inv_work)
321
322 IF (do_svd) THEN
323 ! We have to reset the size of fm_matrix_Minv_L_gamma
324 CALL reset_size_matrix(fm_matrix_minv_l_gamma(1, 1), dimen_ri_red, fm_matrix_l_work%matrix_struct)
325
326 ! L (=fm_matrix_Minv_L_gamma) = S^(-1)*chol_dec(V)
327 CALL parallel_gemm('T', 'N', dimen_ri, dimen_ri_red, dimen_ri, 1.0_dp, fm_matrix_m_inv_work, &
328 fm_matrix_l_work, 0.0_dp, fm_matrix_minv_l_gamma(1, 1))
329 ELSE
330 ! L (=fm_matrix_Minv_L_gamma) = S^(-1)*chol_dec(V)
331 CALL parallel_gemm('T', 'T', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_matrix_m_inv_work, &
332 fm_matrix_l_work, 0.0_dp, fm_matrix_minv_l_gamma(1, 1))
333 END IF
334
335 ! clean the S_inv matrix
336 CALL cp_fm_release(fm_matrix_m_inv_work)
337
338 END IF
339
340 IF (.NOT. do_im_time) THEN
341 IF (do_kpoints) THEN
342 CALL cp_cfm_to_fm(cfm_matrix_minv_l_kpoints(1), mtargetr=fm_matrix_l_work)
343 ELSE
344 CALL cp_fm_to_fm(fm_matrix_minv_l_gamma(1, 1), fm_matrix_l_work)
345 CALL cp_fm_release(fm_matrix_minv_l_gamma)
346 END IF
347
348 END IF
349
350 END IF
351
352 ! If the group distribution changed because of an SVD, we get the new values here
353 CALL get_group_dist(gd_array, color_sub, my_group_l_start, my_group_l_end, my_group_l_size)
354
355 IF (.NOT. do_im_time) THEN
356 IF (unit_nr > 0) THEN
357 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") "RI_INFO| Cholesky decomposition group size:", para_env_l%num_pe
358 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") "RI_INFO| Number of groups for auxiliary basis functions", ngroup
359 IF (calc_pq_cond_num .OR. do_svd) THEN
360 WRITE (unit=unit_nr, fmt="(T3,A,T67,ES14.5)") &
361 "RI_INFO| Condition number of the (P|Q):", cond_num
362 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
363 "RI_INFO| Number of non-positive Eigenvalues of (P|Q):", num_small_eigen
364 END IF
365 CALL m_flush(unit_nr)
366 END IF
367
368 ! replicate the necessary row of the L^{-1} matrix on each proc
369 CALL grep_lcols(fm_matrix_l_work, my_group_l_start, my_group_l_end, my_group_l_size, my_lrows)
370 IF (calc_forces .AND. .NOT. compare_potential_types(qs_env%mp2_env%ri_metric, &
371 qs_env%mp2_env%potential_parameter)) THEN
372 CALL grep_lcols(fm_matrix_v, my_group_l_start, my_group_l_end, my_group_l_size, my_vrows)
373 END IF
374 END IF
375 CALL cp_fm_release(fm_matrix_l_work)
376 IF (calc_forces .AND. .NOT. (do_im_time .OR. compare_potential_types(qs_env%mp2_env%ri_metric, &
377 qs_env%mp2_env%potential_parameter))) CALL cp_fm_release(fm_matrix_v)
378 CALL mp_para_env_release(para_env_l)
379
380 CALL timestop(handle)
381
382 END SUBROUTINE get_2c_integrals
383
384! **************************************************************************************************
385!> \brief ...
386!> \param fm_matrix_L ...
387!> \param do_svd ...
388!> \param eps_svd ...
389!> \param num_small_eigen ...
390!> \param cond_num ...
391!> \param do_inversion ...
392!> \param gd_array ...
393!> \param ngroup ...
394!> \param dimen_RI ...
395!> \param dimen_RI_red ...
396!> \param para_env ...
397! **************************************************************************************************
398 SUBROUTINE decomp_mat_l(fm_matrix_L, do_svd, eps_svd, num_small_eigen, cond_num, do_inversion, gd_array, ngroup, &
399 dimen_RI, dimen_RI_red, para_env)
400
401 TYPE(cp_fm_type), INTENT(INOUT) :: fm_matrix_l
402 LOGICAL, INTENT(IN) :: do_svd
403 REAL(kind=dp), INTENT(IN) :: eps_svd
404 INTEGER, INTENT(INOUT) :: num_small_eigen
405 REAL(kind=dp), INTENT(INOUT) :: cond_num
406 LOGICAL, INTENT(IN) :: do_inversion
407 TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_array
408 INTEGER, INTENT(IN) :: ngroup, dimen_ri
409 INTEGER, INTENT(INOUT) :: dimen_ri_red
410 TYPE(mp_para_env_type), INTENT(IN) :: para_env
411
412 IF (do_svd) THEN
413 CALL matrix_root_with_svd(fm_matrix_l, num_small_eigen, cond_num, eps_svd, do_inversion, para_env)
414
415 dimen_ri_red = dimen_ri - num_small_eigen
416
417 ! We changed the size of fm_matrix_L in matrix_root_with_svd.
418 ! So, we have to get new group sizes
419 CALL release_group_dist(gd_array)
420
421 CALL create_group_dist(gd_array, ngroup, dimen_ri_red)
422 ELSE
423
424 ! calculate Cholesky decomposition L (V=LL^T)
425 CALL cholesky_decomp(fm_matrix_l, dimen_ri, do_inversion=do_inversion)
426
427 IF (do_inversion) CALL invert_mat(fm_matrix_l)
428 END IF
429
430 END SUBROUTINE decomp_mat_l
431
432! **************************************************************************************************
433!> \brief ...
434!> \param qs_env ...
435!> \param cfm_matrix_L_kpoints ...
436!> \param fm_matrix_L_struct ...
437!> \param kpoints ...
438! **************************************************************************************************
439 SUBROUTINE compute_v_by_lattice_sum(qs_env, cfm_matrix_L_kpoints, fm_matrix_L_struct, kpoints)
440 TYPE(qs_environment_type), POINTER :: qs_env
441 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_matrix_l_kpoints
442 TYPE(cp_fm_struct_type), POINTER :: fm_matrix_l_struct
443 TYPE(kpoint_type), POINTER :: kpoints
444
445 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_V_by_lattice_sum'
446
447 INTEGER :: handle, i_dim, ikp, nkp, nkp_extra, &
448 nkp_orig
449 INTEGER, DIMENSION(3) :: nkp_grid_orig, periodic
450 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
451 TYPE(cell_type), POINTER :: cell
452 TYPE(cp_fm_type) :: fm_matrix_l_work
453 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_ri_aux_transl, matrix_v_ri_kp
454 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
455 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
456
457 CALL timeset(routinen, handle)
458
459 NULLIFY (matrix_s_ri_aux_transl, particle_set, cell, qs_kind_set)
460
461 CALL get_qs_env(qs_env=qs_env, &
462 matrix_s_ri_aux_kp=matrix_s_ri_aux_transl, &
463 particle_set=particle_set, &
464 cell=cell, &
465 qs_kind_set=qs_kind_set, &
466 atomic_kind_set=atomic_kind_set)
467
468 ! check that we have a 2n x 2m x 2k mesh to guarantee good convergence
469 CALL get_cell(cell=cell, periodic=periodic)
470 DO i_dim = 1, 3
471 IF (periodic(i_dim) == 1) THEN
472 cpassert(modulo(kpoints%nkp_grid(i_dim), 2) == 0)
473 END IF
474 END DO
475
476 nkp = kpoints%nkp
477
478 ALLOCATE (cfm_matrix_l_kpoints(nkp))
479 DO ikp = 1, nkp
480 CALL cp_cfm_create(cfm_matrix_l_kpoints(ikp), fm_matrix_l_struct)
481 END DO
482 CALL cp_fm_create(fm_matrix_l_work, fm_matrix_l_struct)
483
484 CALL allocate_matrix_v_ri_kp(matrix_v_ri_kp, matrix_s_ri_aux_transl, nkp)
485
486 IF (qs_env%mp2_env%ri_rpa_im_time%do_extrapolate_kpoints) THEN
487
488 nkp_orig = qs_env%mp2_env%ri_rpa_im_time%nkp_orig
489 nkp_extra = qs_env%mp2_env%ri_rpa_im_time%nkp_extra
490
491 CALL build_2c_coulomb_matrix_kp(matrix_v_ri_kp, kpoints, basis_type="RI_AUX", &
492 cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
493 atomic_kind_set=atomic_kind_set, &
494 size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
495 operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp_orig)
496
497 nkp_grid_orig = kpoints%nkp_grid
498 kpoints%nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid_extra(1:3)
499
500 CALL build_2c_coulomb_matrix_kp(matrix_v_ri_kp, kpoints, basis_type="RI_AUX", &
501 cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
502 atomic_kind_set=atomic_kind_set, &
503 size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
504 operator_type=operator_coulomb, ikp_start=nkp_orig + 1, ikp_end=nkp)
505
506 kpoints%nkp_grid = nkp_grid_orig
507
508 ELSE
509
510 CALL build_2c_coulomb_matrix_kp(matrix_v_ri_kp, kpoints, basis_type="RI_AUX", &
511 cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
512 atomic_kind_set=atomic_kind_set, &
513 size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
514 operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp)
515
516 END IF
517
518 DO ikp = 1, nkp
519 CALL copy_dbcsr_pair_to_cfm(matrix_v_ri_kp(ikp, 1)%matrix, matrix_v_ri_kp(ikp, 2)%matrix, &
520 fm_matrix_l_work, cfm_matrix_l_kpoints(ikp))
521 END DO
522 CALL cp_fm_release(fm_matrix_l_work)
523
524 CALL dbcsr_deallocate_matrix_set(matrix_v_ri_kp)
525
526 CALL timestop(handle)
527
528 END SUBROUTINE compute_v_by_lattice_sum
529
530! **************************************************************************************************
531!> \brief ...
532!> \param matrix_v_RI_kp ...
533!> \param matrix_s_RI_aux_transl ...
534!> \param nkp ...
535! **************************************************************************************************
536 SUBROUTINE allocate_matrix_v_ri_kp(matrix_v_RI_kp, matrix_s_RI_aux_transl, nkp)
537
538 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_ri_kp, matrix_s_ri_aux_transl
539 INTEGER :: nkp
540
541 INTEGER :: ikp
542
543 NULLIFY (matrix_v_ri_kp)
544 CALL dbcsr_allocate_matrix_set(matrix_v_ri_kp, nkp, 2)
545
546 DO ikp = 1, nkp
547
548 ALLOCATE (matrix_v_ri_kp(ikp, 1)%matrix)
549 CALL dbcsr_create(matrix_v_ri_kp(ikp, 1)%matrix, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
550 matrix_type=dbcsr_type_no_symmetry)
551 CALL dbcsr_reserve_all_blocks(matrix_v_ri_kp(ikp, 1)%matrix)
552 CALL dbcsr_set(matrix_v_ri_kp(ikp, 1)%matrix, 0.0_dp)
553
554 ALLOCATE (matrix_v_ri_kp(ikp, 2)%matrix)
555 CALL dbcsr_create(matrix_v_ri_kp(ikp, 2)%matrix, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
556 matrix_type=dbcsr_type_no_symmetry)
557 CALL dbcsr_reserve_all_blocks(matrix_v_ri_kp(ikp, 2)%matrix)
558 CALL dbcsr_set(matrix_v_ri_kp(ikp, 2)%matrix, 0.0_dp)
559
560 END DO
561
562 END SUBROUTINE allocate_matrix_v_ri_kp
563
564! **************************************************************************************************
565!> \brief ...
566!> \param qs_env ...
567!> \param fm_matrix_Minv_L_gamma real Gamma-point result ...
568!> \param fm_matrix_L_struct matrix distribution used for the result ...
569!> \param dimen_RI ...
570!> \param ri_metric ...
571!> \param put_mat_KS_env ...
572!> \param regularization_RI ...
573! **************************************************************************************************
574 SUBROUTINE ri_2c_integral_mat(qs_env, fm_matrix_Minv_L_gamma, fm_matrix_L_struct, dimen_RI, ri_metric, &
575 put_mat_KS_env, regularization_RI)
576
577 TYPE(qs_environment_type), POINTER :: qs_env
578 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_minv_l_gamma
579 TYPE(cp_fm_struct_type), INTENT(IN), POINTER :: fm_matrix_l_struct
580 INTEGER, INTENT(IN) :: dimen_ri
581 TYPE(libint_potential_type) :: ri_metric
582 LOGICAL, OPTIONAL :: put_mat_ks_env
583 REAL(kind=dp), OPTIONAL :: regularization_ri
584
585 CALL ri_2c_integral_mat_core(qs_env, fm_matrix_minv_l_gamma=fm_matrix_minv_l_gamma, &
586 fm_matrix_l_struct=fm_matrix_l_struct, dimen_ri=dimen_ri, ri_metric=ri_metric, &
587 do_kpoints=.false., put_mat_ks_env=put_mat_ks_env, &
588 regularization_ri=regularization_ri)
589
590 END SUBROUTINE ri_2c_integral_mat
591
592! **************************************************************************************************
593!> \brief Computes the complex k-point RI metric as native CFM matrices.
594!> \param qs_env ...
595!> \param cfm_matrix_M_kpoints ...
596!> \param fm_matrix_struct ...
597!> \param dimen_RI ...
598!> \param ri_metric ...
599!> \param kpoints ...
600!> \param put_mat_KS_env ...
601!> \param regularization_RI ...
602!> \param ikp_ext ...
603!> \param do_build_cell_index ...
604! **************************************************************************************************
605 SUBROUTINE ri_2c_integral_mat_kp(qs_env, cfm_matrix_M_kpoints, fm_matrix_struct, dimen_RI, ri_metric, &
606 kpoints, put_mat_KS_env, regularization_RI, ikp_ext, do_build_cell_index)
607
608 TYPE(qs_environment_type), POINTER :: qs_env
609 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_matrix_m_kpoints
610 TYPE(cp_fm_struct_type), INTENT(IN), POINTER :: fm_matrix_struct
611 INTEGER, INTENT(IN) :: dimen_ri
612 TYPE(libint_potential_type) :: ri_metric
613 TYPE(kpoint_type), POINTER :: kpoints
614 LOGICAL, OPTIONAL :: put_mat_ks_env
615 REAL(kind=dp), OPTIONAL :: regularization_ri
616 INTEGER, OPTIONAL :: ikp_ext
617 LOGICAL, OPTIONAL :: do_build_cell_index
618
619 CALL ri_2c_integral_mat_core(qs_env, cfm_matrix_m_kpoints=cfm_matrix_m_kpoints, &
620 fm_matrix_l_struct=fm_matrix_struct, dimen_ri=dimen_ri, ri_metric=ri_metric, &
621 do_kpoints=.true., kpoints=kpoints, put_mat_ks_env=put_mat_ks_env, &
622 regularization_ri=regularization_ri, ikp_ext=ikp_ext, &
623 do_build_cell_index=do_build_cell_index)
624
625 END SUBROUTINE ri_2c_integral_mat_kp
626
627! **************************************************************************************************
628!> \brief Shared DBCSR evaluator with representation-specific matrix output adapters.
629!> \param qs_env ...
630!> \param fm_matrix_Minv_L_gamma ...
631!> \param cfm_matrix_M_kpoints ...
632!> \param fm_matrix_L_struct ...
633!> \param dimen_RI ...
634!> \param ri_metric ...
635!> \param do_kpoints ...
636!> \param kpoints ...
637!> \param put_mat_KS_env ...
638!> \param regularization_RI ...
639!> \param ikp_ext ...
640!> \param do_build_cell_index ...
641! **************************************************************************************************
642 SUBROUTINE ri_2c_integral_mat_core(qs_env, fm_matrix_Minv_L_gamma, cfm_matrix_M_kpoints, &
643 fm_matrix_L_struct, dimen_RI, ri_metric, do_kpoints, kpoints, &
644 put_mat_KS_env, regularization_RI, ikp_ext, do_build_cell_index)
645
646 TYPE(qs_environment_type), POINTER :: qs_env
647 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :), &
648 INTENT(OUT), OPTIONAL :: fm_matrix_minv_l_gamma
649 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:), &
650 INTENT(OUT), OPTIONAL :: cfm_matrix_m_kpoints
651 TYPE(cp_fm_struct_type), INTENT(IN), POINTER :: fm_matrix_l_struct
652 INTEGER, INTENT(IN) :: dimen_ri
653 TYPE(libint_potential_type) :: ri_metric
654 LOGICAL, INTENT(IN) :: do_kpoints
655 TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
656 LOGICAL, OPTIONAL :: put_mat_ks_env
657 REAL(kind=dp), OPTIONAL :: regularization_ri
658 INTEGER, OPTIONAL :: ikp_ext
659 LOGICAL, OPTIONAL :: do_build_cell_index
660
661 CHARACTER(LEN=*), PARAMETER :: routinen = 'RI_2c_integral_mat_core'
662
663 INTEGER :: handle, ikp, ikp_for_xkp, img, natom, &
664 nimg, nkind, nkp
665 INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_ri
666 INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
667 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
668 LOGICAL :: my_do_build_cell_index, my_put_mat_ks_env
669 REAL(kind=dp) :: my_regularization_ri
670 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
671 TYPE(cp_blacs_env_type), POINTER :: blacs_env
672 TYPE(cp_cfm_type) :: cfm_matrix_s_global
673 TYPE(cp_fm_struct_type), POINTER :: fm_struct
674 TYPE(cp_fm_type) :: fm_matrix_s_global
675 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_m_kpoints
676 TYPE(dbcsr_distribution_type) :: dbcsr_dist
677 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_ri_aux_transl
678 TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: mat_2c
679 TYPE(dbcsr_type), POINTER :: cmatrix, matrix_s_ri_aux_desymm, &
680 rmatrix, tmpmat
681 TYPE(dft_control_type), POINTER :: dft_control
682 TYPE(distribution_2d_type), POINTER :: dist_2d
683 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ri
684 TYPE(mp_para_env_type), POINTER :: para_env
685 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
686 POINTER :: sab_ri
687 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
688 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
689
690 CALL timeset(routinen, handle)
691
692 IF (do_kpoints) THEN
693 cpassert(PRESENT(cfm_matrix_m_kpoints))
694 ELSE
695 cpassert(PRESENT(fm_matrix_minv_l_gamma))
696 END IF
697
698 NULLIFY (sab_ri, matrix_s_ri_aux_transl, dist_2d)
699
700 IF (PRESENT(regularization_ri)) THEN
701 my_regularization_ri = regularization_ri
702 ELSE
703 my_regularization_ri = 0.0_dp
704 END IF
705
706 IF (PRESENT(put_mat_ks_env)) THEN
707 my_put_mat_ks_env = put_mat_ks_env
708 ELSE
709 my_put_mat_ks_env = .false.
710 END IF
711
712 IF (PRESENT(do_build_cell_index)) THEN
713 my_do_build_cell_index = do_build_cell_index
714 ELSE
715 my_do_build_cell_index = .false.
716 END IF
717
718 CALL get_qs_env(qs_env=qs_env, &
719 para_env=para_env, &
720 blacs_env=blacs_env, &
721 nkind=nkind, &
722 distribution_2d=dist_2d, &
723 qs_kind_set=qs_kind_set, &
724 particle_set=particle_set, &
725 dft_control=dft_control, &
726 natom=natom)
727
728 ALLOCATE (sizes_ri(natom))
729 ALLOCATE (basis_set_ri(nkind))
730
731 CALL basis_set_list_setup(basis_set_ri, 'RI_AUX', qs_kind_set)
732 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri)
733
734 CALL build_2c_neighbor_lists(sab_ri, basis_set_ri, basis_set_ri, ri_metric, "RPA_2c_nl_RI", qs_env, &
735 sym_ij=.true., dist_2d=dist_2d)
736
737 CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
738 ALLOCATE (row_bsize(SIZE(sizes_ri)))
739 ALLOCATE (col_bsize(SIZE(sizes_ri)))
740 row_bsize(:) = sizes_ri
741 col_bsize(:) = sizes_ri
742
743 IF (do_kpoints) THEN
744 cpassert(PRESENT(kpoints))
745 IF (my_do_build_cell_index) THEN
746 CALL kpoint_init_cell_index(kpoints, sab_ri, para_env, dft_control%nimages)
747 END IF
748 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, &
749 cell_to_index=cell_to_index)
750 nimg = dft_control%nimages
751 ELSE
752 nkp = 1
753 nimg = 1
754 END IF
755
756 ALLOCATE (mat_2c(nimg))
757 CALL dbcsr_create(mat_2c(1), "(RI|RI)", dbcsr_dist, dbcsr_type_symmetric, &
758 row_bsize, col_bsize)
759 DEALLOCATE (row_bsize, col_bsize)
760
761 DO img = 2, nimg
762 CALL dbcsr_create(mat_2c(img), template=mat_2c(1))
763 END DO
764
765 CALL build_2c_integrals(mat_2c, 0.0_dp, qs_env, sab_ri, basis_set_ri, basis_set_ri, &
766 ri_metric, do_kpoints=do_kpoints, ext_kpoints=kpoints, &
767 regularization_ri=my_regularization_ri)
768
769 CALL dbcsr_distribution_release(dbcsr_dist)
770 DEALLOCATE (basis_set_ri)
771
772 IF (my_put_mat_ks_env) THEN
773 CALL get_ks_env(qs_env%ks_env, matrix_s_ri_aux_kp=matrix_s_ri_aux_transl)
774 IF (ASSOCIATED(matrix_s_ri_aux_transl)) CALL dbcsr_deallocate_matrix_set(matrix_s_ri_aux_transl)
775 END IF
776 CALL dbcsr_allocate_matrix_set(matrix_s_ri_aux_transl, 1, nimg)
777 DO img = 1, nimg
778 ALLOCATE (matrix_s_ri_aux_transl(1, img)%matrix)
779 CALL dbcsr_copy(matrix_s_ri_aux_transl(1, img)%matrix, mat_2c(img))
780 CALL dbcsr_release(mat_2c(img))
781 END DO
782
783 IF (my_put_mat_ks_env) THEN
784 CALL set_ks_env(qs_env%ks_env, matrix_s_ri_aux_kp=matrix_s_ri_aux_transl)
785 END IF
786
787 IF (PRESENT(ikp_ext)) nkp = 1
788
789 IF (do_kpoints) THEN
790 ALLOCATE (cfm_matrix_m_kpoints(nkp))
791 DO ikp = 1, nkp
792 CALL cp_cfm_create(cfm_matrix_m_kpoints(ikp), fm_matrix_l_struct)
793 END DO
794 ELSE
795 ALLOCATE (fm_matrix_m_kpoints(1, 1))
796 CALL cp_fm_create(fm_matrix_m_kpoints(1, 1), fm_matrix_l_struct)
797 CALL cp_fm_set_all(fm_matrix_m_kpoints(1, 1), 0.0_dp)
798 END IF
799
800 NULLIFY (fm_struct)
801 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_ri, &
802 ncol_global=dimen_ri, para_env=para_env)
803
804 CALL cp_fm_create(fm_matrix_s_global, fm_struct)
805 CALL cp_fm_set_all(fm_matrix_s_global, 0.0_dp)
806 IF (do_kpoints) CALL cp_cfm_create(cfm_matrix_s_global, fm_struct, set_zero=.true.)
807
808 IF (do_kpoints) THEN
809
810 ALLOCATE (rmatrix, cmatrix, tmpmat)
811 CALL dbcsr_create(rmatrix, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
812 matrix_type=dbcsr_type_symmetric)
813 CALL dbcsr_create(cmatrix, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
814 matrix_type=dbcsr_type_antisymmetric)
815 CALL dbcsr_create(tmpmat, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
816 matrix_type=dbcsr_type_no_symmetry)
817 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_ri)
818 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_ri)
819
820 DO ikp = 1, nkp
821
822 ! s matrix is not spin dependent, double the work
823 CALL dbcsr_set(rmatrix, 0.0_dp)
824 CALL dbcsr_set(cmatrix, 0.0_dp)
825
826 IF (PRESENT(ikp_ext)) THEN
827 ikp_for_xkp = ikp_ext
828 ELSE
829 ikp_for_xkp = ikp
830 END IF
831
832 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s_ri_aux_transl, ispin=1, &
833 xkp=xkp(1:3, ikp_for_xkp), cell_to_index=cell_to_index, sab_nl=sab_ri)
834
835 CALL dbcsr_set(tmpmat, 0.0_dp)
836 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
837
838 CALL cp_fm_set_all(fm_matrix_s_global, 0.0_dp)
839 CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_s_global)
840 CALL cp_fm_to_cfm(msourcer=fm_matrix_s_global, mtarget=cfm_matrix_s_global)
841
842 CALL dbcsr_set(tmpmat, 0.0_dp)
843 CALL dbcsr_desymmetrize(cmatrix, tmpmat)
844
845 CALL cp_fm_set_all(fm_matrix_s_global, 0.0_dp)
846 CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_s_global)
847 CALL cp_cfm_scale_and_add_fm(z_one, cfm_matrix_s_global, gaussi, fm_matrix_s_global)
848 CALL cp_cfm_copy_general(cfm_matrix_s_global, cfm_matrix_m_kpoints(ikp), para_env)
849
850 END DO
851
852 CALL dbcsr_deallocate_matrix(rmatrix)
853 CALL dbcsr_deallocate_matrix(cmatrix)
854 CALL dbcsr_deallocate_matrix(tmpmat)
855
856 ELSE
857
858 NULLIFY (matrix_s_ri_aux_desymm)
859 ALLOCATE (matrix_s_ri_aux_desymm)
860 CALL dbcsr_create(matrix_s_ri_aux_desymm, template=matrix_s_ri_aux_transl(1, 1)%matrix, &
861 name='S_RI non_symm', matrix_type=dbcsr_type_no_symmetry)
862
863 CALL dbcsr_desymmetrize(matrix_s_ri_aux_transl(1, 1)%matrix, matrix_s_ri_aux_desymm)
864
865 CALL copy_dbcsr_to_fm(matrix_s_ri_aux_desymm, fm_matrix_s_global)
866
867 CALL cp_fm_copy_general(fm_matrix_s_global, fm_matrix_m_kpoints(1, 1), para_env)
868
869 CALL dbcsr_deallocate_matrix(matrix_s_ri_aux_desymm)
870
871 END IF
872
873 CALL release_neighbor_list_sets(sab_ri)
874
875 CALL cp_fm_struct_release(fm_struct)
876
877 CALL cp_fm_release(fm_matrix_s_global)
878
879 IF (do_kpoints) THEN
880 CALL cp_cfm_release(cfm_matrix_s_global)
881 ELSE
882 cpassert(PRESENT(fm_matrix_minv_l_gamma))
883 CALL move_alloc(fm_matrix_m_kpoints, fm_matrix_minv_l_gamma)
884 END IF
885
886 IF (.NOT. my_put_mat_ks_env) THEN
887 CALL dbcsr_deallocate_matrix_set(matrix_s_ri_aux_transl)
888 END IF
889
890 CALL timestop(handle)
891
892 END SUBROUTINE ri_2c_integral_mat_core
893
894! **************************************************************************************************
895!> \brief Copies a real/imaginary DBCSR pair into its native complex full-matrix representation.
896!> \param matrix_re ...
897!> \param matrix_im ...
898!> \param fm_work ...
899!> \param cfm_matrix ...
900! **************************************************************************************************
901 SUBROUTINE copy_dbcsr_pair_to_cfm(matrix_re, matrix_im, fm_work, cfm_matrix)
902 TYPE(dbcsr_type), POINTER :: matrix_re, matrix_im
903 TYPE(cp_fm_type), INTENT(INOUT) :: fm_work
904 TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_matrix
905
906 CALL cp_fm_set_all(fm_work, 0.0_dp)
907 CALL copy_dbcsr_to_fm(matrix_re, fm_work)
908 CALL cp_fm_to_cfm(msourcer=fm_work, mtarget=cfm_matrix)
909
910 CALL cp_fm_set_all(fm_work, 0.0_dp)
911 CALL copy_dbcsr_to_fm(matrix_im, fm_work)
912 CALL cp_cfm_scale_and_add_fm(z_one, cfm_matrix, gaussi, fm_work)
913 END SUBROUTINE copy_dbcsr_pair_to_cfm
914
915! **************************************************************************************************
916!> \brief ...
917!> \param qs_env ...
918!> \param eri_method ...
919!> \param eri_param ...
920!> \param para_env ...
921!> \param para_env_sub ...
922!> \param para_env_L ...
923!> \param mp2_memory ...
924!> \param fm_matrix_L ...
925!> \param ngroup ...
926!> \param color_sub ...
927!> \param dimen_RI ...
928!> \param my_group_L_size ...
929!> \param my_group_L_start ...
930!> \param my_group_L_end ...
931!> \param gd_array ...
932!> \param calc_PQ_cond_num ...
933!> \param cond_num ...
934!> \param num_small_eigen ...
935!> \param potential ...
936!> \param sab_orb_sub ...
937!> \param do_im_time ...
938!> \param fm_matrix_L_struct ...
939! **************************************************************************************************
940 SUBROUTINE compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
941 fm_matrix_L, ngroup, color_sub, dimen_RI, &
942 my_group_L_size, my_group_L_start, my_group_L_end, &
943 gd_array, calc_PQ_cond_num, cond_num, num_small_eigen, potential, &
944 sab_orb_sub, do_im_time, fm_matrix_L_struct)
945
946 TYPE(qs_environment_type), POINTER :: qs_env
947 INTEGER, INTENT(IN) :: eri_method
948 TYPE(cp_eri_mme_param), POINTER :: eri_param
949 TYPE(mp_para_env_type), INTENT(IN) :: para_env
950 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
951 TYPE(mp_para_env_type), INTENT(OUT), POINTER :: para_env_l
952 REAL(kind=dp), INTENT(IN) :: mp2_memory
953 TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_l
954 INTEGER, INTENT(IN) :: ngroup, color_sub
955 INTEGER, INTENT(OUT) :: dimen_ri, my_group_l_size, &
956 my_group_l_start, my_group_l_end
957 TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array
958 LOGICAL, INTENT(IN) :: calc_pq_cond_num
959 REAL(kind=dp), INTENT(OUT) :: cond_num
960 INTEGER, INTENT(OUT) :: num_small_eigen
961 TYPE(libint_potential_type), INTENT(IN) :: potential
962 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
963 POINTER :: sab_orb_sub
964 LOGICAL, INTENT(IN), OPTIONAL :: do_im_time
965 TYPE(cp_fm_struct_type), OPTIONAL, POINTER :: fm_matrix_l_struct
966
967 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_2c_integrals'
968
969 INTEGER :: best_group_size, color_l, group_size, handle, handle2, i_global, iatom, iib, &
970 ikind, iproc, j_global, jjb, natom, ncol_local, nkind, nrow_local, nsgf, potential_type, &
971 proc_receive, proc_receive_static, proc_send_static, proc_shift, rec_l_end, rec_l_size, &
972 rec_l_start, strat_group_size
973 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
974 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
975 LOGICAL :: my_do_im_time
976 REAL(kind=dp) :: min_mem_for_qk
977 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: egen_l
978 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: l_external_col, l_local_col
979 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
980 TYPE(cp_blacs_env_type), POINTER :: blacs_env_l
981 TYPE(cp_fm_type) :: fm_matrix_l_diag
982 TYPE(group_dist_d1_type) :: gd_sub_array
983 TYPE(gto_basis_set_type), POINTER :: basis_set_a
984 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
985 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
986
987 CALL timeset(routinen, handle)
988
989 my_do_im_time = .false.
990 IF (PRESENT(do_im_time)) THEN
991 my_do_im_time = do_im_time
992 END IF
993
994 ! get stuff
995 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
996 particle_set=particle_set)
997
998 nkind = SIZE(qs_kind_set)
999 natom = SIZE(particle_set)
1000
1001 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
1002
1003 DO ikind = 1, nkind
1004 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
1005 cpassert(ASSOCIATED(basis_set_a))
1006 END DO
1007
1008 dimen_ri = 0
1009 DO iatom = 1, natom
1010 ikind = kind_of(iatom)
1011 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
1012 dimen_ri = dimen_ri + nsgf
1013 END DO
1014
1015 ! check that very small systems are not running on too many processes
1016 IF (dimen_ri < ngroup) THEN
1017 CALL cp_abort(__location__, "Product of block size and number "// &
1018 "of RI functions should not exceed total number of processes")
1019 END IF
1020
1021 CALL create_group_dist(gd_array, ngroup, dimen_ri)
1022
1023 CALL get_group_dist(gd_array, color_sub, my_group_l_start, my_group_l_end, my_group_l_size)
1024
1025 CALL timeset(routinen//"_loop_lm", handle2)
1026
1027 ALLOCATE (l_local_col(dimen_ri, my_group_l_size))
1028 l_local_col = 0.0_dp
1029
1030 potential_type = potential%potential_type
1031
1032 IF ((eri_method == do_eri_mme .OR. eri_method == do_eri_os) &
1033 .AND. potential_type == do_potential_coulomb .OR. &
1034 (eri_method == do_eri_mme .AND. potential_type == do_potential_long)) THEN
1035
1036 CALL mp2_eri_2c_integrate(eri_param, potential, para_env_sub, qs_env, &
1037 basis_type_a="RI_AUX", basis_type_b="RI_AUX", &
1038 hab=l_local_col, first_b=my_group_l_start, last_b=my_group_l_end, &
1039 eri_method=eri_method)
1040
1041 ELSE IF (eri_method == do_eri_gpw .OR. &
1042 (potential_type == do_potential_long .AND. qs_env%mp2_env%eri_method == do_eri_os) &
1043 .OR. (potential_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme)) THEN
1044
1045 CALL mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_l_start, my_group_l_end, &
1046 natom, potential, sab_orb_sub, l_local_col, kind_of)
1047
1048 ELSE
1049 cpabort("unknown ERI method")
1050 END IF
1051
1052 CALL timestop(handle2)
1053
1054 ! split the total number of proc in a subgroup of the size of ~1/10 of the
1055 ! total num of proc
1056 best_group_size = para_env%num_pe
1057
1058 strat_group_size = max(1, para_env%num_pe/10)
1059
1060 min_mem_for_qk = real(dimen_ri, kind=dp)*dimen_ri*3.0_dp*8.0_dp/1024_dp/1024_dp
1061
1062 group_size = strat_group_size - 1
1063 DO iproc = strat_group_size, para_env%num_pe
1064 group_size = group_size + 1
1065 ! check that group_size is a multiple of sub_group_size and a divisor of
1066 ! the total num of proc
1067 IF (mod(para_env%num_pe, group_size) /= 0 .OR. mod(group_size, para_env_sub%num_pe) /= 0) cycle
1068
1069 ! check for memory
1070 IF (real(group_size, kind=dp)*mp2_memory < min_mem_for_qk) cycle
1071
1072 best_group_size = group_size
1073 EXIT
1074 END DO
1075
1076 IF (my_do_im_time) THEN
1077 ! para_env_L is the para_env for im. time to avoid bug
1078 best_group_size = para_env%num_pe
1079 END IF
1080
1081 ! create the L group
1082 color_l = para_env%mepos/best_group_size
1083 ALLOCATE (para_env_l)
1084 CALL para_env_l%from_split(para_env, color_l)
1085
1086 ! create the blacs_L
1087 NULLIFY (blacs_env_l)
1088 CALL cp_blacs_env_create(blacs_env=blacs_env_l, para_env=para_env_l)
1089
1090 ! create the full matrix L defined in the L group
1091 CALL create_matrix_l(fm_matrix_l, blacs_env_l, dimen_ri, para_env_l, "fm_matrix_L", fm_matrix_l_struct)
1092
1093 IF (my_do_im_time .AND. para_env%num_pe > 1) THEN
1094
1095 CALL cp_array_columns_to_fm(fm_matrix_l, l_local_col, para_env, &
1096 my_group_l_start, my_group_l_end, &
1097 dimen_ri)
1098
1099 ELSE
1100 block
1101 TYPE(mp_comm_type) :: comm_exchange
1102
1103 CALL comm_exchange%from_split(para_env_l, para_env_sub%mepos)
1104
1105 CALL create_group_dist(gd_sub_array, my_group_l_start, &
1106 my_group_l_end, my_group_l_size, comm_exchange)
1107
1108 CALL cp_fm_get_info(matrix=fm_matrix_l, &
1109 nrow_local=nrow_local, &
1110 ncol_local=ncol_local, &
1111 row_indices=row_indices, &
1112 col_indices=col_indices)
1113
1114 DO jjb = 1, ncol_local
1115 j_global = col_indices(jjb)
1116 IF (j_global >= my_group_l_start .AND. j_global <= my_group_l_end) THEN
1117 DO iib = 1, nrow_local
1118 i_global = row_indices(iib)
1119 fm_matrix_l%local_data(iib, jjb) = l_local_col(i_global, j_global - my_group_l_start + 1)
1120 END DO
1121 END IF
1122 END DO
1123
1124 proc_send_static = modulo(comm_exchange%mepos + 1, comm_exchange%num_pe)
1125 proc_receive_static = modulo(comm_exchange%mepos - 1, comm_exchange%num_pe)
1126
1127 DO proc_shift = 1, comm_exchange%num_pe - 1
1128 proc_receive = modulo(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
1129
1130 CALL get_group_dist(gd_sub_array, proc_receive, rec_l_start, rec_l_end, rec_l_size)
1131
1132 ALLOCATE (l_external_col(dimen_ri, rec_l_size))
1133 l_external_col = 0.0_dp
1134
1135 CALL comm_exchange%sendrecv(l_local_col, proc_send_static, l_external_col, proc_receive_static)
1136
1137 DO jjb = 1, ncol_local
1138 j_global = col_indices(jjb)
1139 IF (j_global >= rec_l_start .AND. j_global <= rec_l_end) THEN
1140 DO iib = 1, nrow_local
1141 i_global = row_indices(iib)
1142 fm_matrix_l%local_data(iib, jjb) = l_external_col(i_global, j_global - rec_l_start + 1)
1143 END DO
1144 END IF
1145 END DO
1146
1147 CALL move_alloc(l_external_col, l_local_col)
1148 END DO
1149
1150 CALL release_group_dist(gd_sub_array)
1151 CALL comm_exchange%free()
1152 END block
1153 END IF
1154
1155 DEALLOCATE (l_local_col)
1156
1157 ! create the new group for the sum of the local data over all processes
1158 block
1159 TYPE(mp_comm_type) :: comm_exchange
1160 comm_exchange = fm_matrix_l%matrix_struct%context%interconnect(para_env)
1161 CALL comm_exchange%sum(fm_matrix_l%local_data)
1162 CALL comm_exchange%free()
1163 END block
1164
1165 cond_num = 1.0_dp
1166 num_small_eigen = 0
1167 IF (calc_pq_cond_num) THEN
1168 ! calculate the condition number of the (P|Q) matrix
1169 ! create a copy of the matrix
1170
1171 CALL create_matrix_l(fm_matrix_l_diag, blacs_env_l, dimen_ri, para_env_l, "fm_matrix_L_diag", fm_matrix_l_struct)
1172
1173 CALL cp_fm_to_fm(source=fm_matrix_l, destination=fm_matrix_l_diag)
1174
1175 ALLOCATE (egen_l(dimen_ri))
1176
1177 egen_l = 0.0_dp
1178 CALL cp_fm_syevx(matrix=fm_matrix_l_diag, eigenvalues=egen_l)
1179
1180 num_small_eigen = 0
1181 DO iib = 1, dimen_ri
1182 IF (abs(egen_l(iib)) < 0.001_dp) num_small_eigen = num_small_eigen + 1
1183 END DO
1184
1185 cond_num = maxval(abs(egen_l))/minval(abs(egen_l))
1186
1187 CALL cp_fm_release(fm_matrix_l_diag)
1188
1189 DEALLOCATE (egen_l)
1190 END IF
1191
1192 ! release blacs_env
1193 CALL cp_blacs_env_release(blacs_env_l)
1194
1195 CALL timestop(handle)
1196
1197 END SUBROUTINE compute_2c_integrals
1198
1199! **************************************************************************************************
1200!> \brief ...
1201!> \param matrix ...
1202!> \param num_small_evals ...
1203!> \param cond_num ...
1204!> \param eps_svd ...
1205!> \param do_inversion ...
1206!> \param para_env ...
1207! **************************************************************************************************
1208 SUBROUTINE matrix_root_with_svd(matrix, num_small_evals, cond_num, eps_svd, do_inversion, para_env)
1209 TYPE(cp_fm_type), INTENT(INOUT) :: matrix
1210 INTEGER, INTENT(OUT) :: num_small_evals
1211 REAL(kind=dp), INTENT(OUT) :: cond_num
1212 REAL(kind=dp), INTENT(IN) :: eps_svd
1213 LOGICAL, INTENT(IN) :: do_inversion
1214 TYPE(mp_para_env_type), INTENT(IN) :: para_env
1215
1216 CHARACTER(LEN=*), PARAMETER :: routinen = 'matrix_root_with_svd'
1217
1218 INTEGER :: group_size_l, handle, ii, needed_evals, &
1219 nrow, pos_max(1)
1220 INTEGER, ALLOCATABLE, DIMENSION(:) :: num_eval
1221 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
1222 TYPE(cp_fm_type) :: evecs
1223 TYPE(mp_comm_type) :: comm_exchange
1224
1225 CALL timeset(routinen, handle)
1226
1227 ! Create arrays carrying eigenvalues and eigenvectors later
1228 CALL cp_fm_get_info(matrix=matrix, nrow_global=nrow)
1229
1230 ALLOCATE (evals(nrow))
1231 evals = 0
1232
1233 CALL cp_fm_create(evecs, matrix%matrix_struct)
1234
1235 ! Perform the EVD (=SVD of a positive semidefinite matrix)
1236 CALL choose_eigv_solver(matrix, evecs, evals)
1237
1238 ! Determine the number of neglectable eigenvalues assuming that the eigenvalues are in ascending order
1239 num_small_evals = 0
1240 DO ii = 1, nrow
1241 IF (evals(ii) > eps_svd) THEN
1242 num_small_evals = ii - 1
1243 EXIT
1244 END IF
1245 END DO
1246 needed_evals = nrow - num_small_evals
1247
1248 ! Get the condition number w.r.t. considered eigenvalues
1249 cond_num = evals(nrow)/evals(num_small_evals + 1)
1250
1251 ! Determine the eigenvalues of the request matrix root or its inverse
1252 evals(1:num_small_evals) = 0.0_dp
1253 IF (do_inversion) THEN
1254 evals(num_small_evals + 1:nrow) = 1.0_dp/sqrt(evals(num_small_evals + 1:nrow))
1255 ELSE
1256 evals(num_small_evals + 1:nrow) = sqrt(evals(num_small_evals + 1:nrow))
1257 END IF
1258
1259 CALL cp_fm_column_scale(evecs, evals)
1260
1261 ! As it turns out, the results in the subgroups differ.
1262 ! Thus, we have to equilibrate the results.
1263 ! Step 1: Get a communicator connecting processes with same id within subgroups
1264 ! use that group_size_L is divisible by the total number of processes (see above)
1265 group_size_l = para_env%num_pe/matrix%matrix_struct%para_env%num_pe
1266 comm_exchange = matrix%matrix_struct%context%interconnect(para_env)
1267
1268 ALLOCATE (num_eval(0:group_size_l - 1))
1269 CALL comm_exchange%allgather(num_small_evals, num_eval)
1270
1271 num_small_evals = minval(num_eval)
1272
1273 IF (num_small_evals /= maxval(num_eval)) THEN
1274 ! Step 2: Get position of maximum value
1275 pos_max = maxloc(num_eval)
1276 num_small_evals = num_eval(pos_max(1))
1277 needed_evals = nrow - num_small_evals
1278
1279 ! Step 3: Broadcast your local data to all other processes
1280 CALL comm_exchange%bcast(evecs%local_data, pos_max(1))
1281 CALL comm_exchange%bcast(cond_num, pos_max(1))
1282 END IF
1283
1284 DEALLOCATE (num_eval)
1285
1286 CALL comm_exchange%free()
1287
1288 CALL reset_size_matrix(matrix, needed_evals, matrix%matrix_struct)
1289
1290 ! Copy the needed eigenvectors
1291 CALL cp_fm_to_fm(evecs, matrix, needed_evals, num_small_evals + 1)
1292
1293 CALL cp_fm_release(evecs)
1294
1295 CALL timestop(handle)
1296
1297 END SUBROUTINE matrix_root_with_svd
1298
1299! **************************************************************************************************
1300!> \brief ...
1301!> \param matrix ...
1302!> \param new_size ...
1303!> \param fm_struct_template ...
1304! **************************************************************************************************
1305 SUBROUTINE reset_size_matrix(matrix, new_size, fm_struct_template)
1306 TYPE(cp_fm_type), INTENT(INOUT) :: matrix
1307 INTEGER, INTENT(IN) :: new_size
1308 TYPE(cp_fm_struct_type), POINTER :: fm_struct_template
1309
1310 CHARACTER(LEN=*), PARAMETER :: routinen = 'reset_size_matrix'
1311
1312 INTEGER :: handle
1313 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1314
1315 CALL timeset(routinen, handle)
1316
1317 ! Choose the block sizes as large as in the template for the later steps
1318 NULLIFY (fm_struct)
1319 CALL cp_fm_struct_create(fm_struct, ncol_global=new_size, template_fmstruct=fm_struct_template, force_block=.true.)
1320
1321 CALL cp_fm_release(matrix)
1322
1323 CALL cp_fm_create(matrix, fm_struct)
1324 CALL cp_fm_set_all(matrix, 0.0_dp)
1325
1326 CALL cp_fm_struct_release(fm_struct)
1327
1328 CALL timestop(handle)
1329
1330 END SUBROUTINE reset_size_matrix
1331
1332! **************************************************************************************************
1333!> \brief ...
1334!> \param fm_matrix_L ...
1335!> \param dimen_RI ...
1336!> \param do_inversion ...
1337! **************************************************************************************************
1338 SUBROUTINE cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion)
1339
1340 TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_l
1341 INTEGER, INTENT(IN) :: dimen_ri
1342 LOGICAL, INTENT(IN) :: do_inversion
1343
1344 CHARACTER(LEN=*), PARAMETER :: routinen = 'cholesky_decomp'
1345
1346 INTEGER :: handle, i_global, iib, info_chol, &
1347 j_global, jjb, ncol_local, nrow_local
1348 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1349
1350 CALL timeset(routinen, handle)
1351
1352 ! do cholesky decomposition
1353 CALL cp_fm_cholesky_decompose(matrix=fm_matrix_l, n=dimen_ri, info_out=info_chol)
1354 cpassert(info_chol == 0)
1355
1356 IF (.NOT. do_inversion) THEN
1357 ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
1358 CALL cp_fm_get_info(matrix=fm_matrix_l, &
1359 nrow_local=nrow_local, &
1360 ncol_local=ncol_local, &
1361 row_indices=row_indices, &
1362 col_indices=col_indices)
1363 DO jjb = 1, ncol_local
1364 j_global = col_indices(jjb)
1365 DO iib = 1, nrow_local
1366 i_global = row_indices(iib)
1367 IF (j_global < i_global) fm_matrix_l%local_data(iib, jjb) = 0.0_dp
1368 END DO
1369 END DO
1370
1371 END IF
1372
1373 CALL timestop(handle)
1374
1375 END SUBROUTINE cholesky_decomp
1376
1377! **************************************************************************************************
1378!> \brief ...
1379!> \param cfm_matrix_M_kpoints ...
1380!> \param cfm_matrix_L_kpoints ...
1381!> \param cfm_matrix_Minv_L_kpoints ...
1382!> \param dimen_RI ...
1383!> \param kpoints ...
1384!> \param eps_eigval_S ...
1385! **************************************************************************************************
1386 SUBROUTINE inversion_of_m_and_mult_with_chol_dec_of_v(cfm_matrix_M_kpoints, cfm_matrix_L_kpoints, &
1387 cfm_matrix_Minv_L_kpoints, dimen_RI, kpoints, eps_eigval_S)
1388
1389 TYPE(cp_cfm_type), DIMENSION(:), INTENT(IN) :: cfm_matrix_m_kpoints
1390 TYPE(cp_cfm_type), DIMENSION(:), INTENT(INOUT) :: cfm_matrix_l_kpoints
1391 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:), &
1392 INTENT(INOUT) :: cfm_matrix_minv_l_kpoints
1393 INTEGER, INTENT(IN) :: dimen_ri
1394 TYPE(kpoint_type), POINTER :: kpoints
1395 REAL(kind=dp), INTENT(IN) :: eps_eigval_s
1396
1397 CHARACTER(LEN=*), PARAMETER :: routinen = 'inversion_of_M_and_mult_with_chol_dec_of_V'
1398
1399 INTEGER :: handle, ikp, nkp
1400 TYPE(cp_cfm_type) :: cfm_matrix_m_tmp
1401 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1402
1403 CALL timeset(routinen, handle)
1404
1405 matrix_struct => cfm_matrix_m_kpoints(1)%matrix_struct
1406
1407 CALL cp_cfm_create(cfm_matrix_m_tmp, matrix_struct)
1408
1409 CALL get_kpoint_info(kpoints, nkp=nkp)
1410 ALLOCATE (cfm_matrix_minv_l_kpoints(nkp))
1411
1412 DO ikp = 1, nkp
1413
1414 CALL cp_cfm_to_cfm(cfm_matrix_m_kpoints(ikp), cfm_matrix_m_tmp)
1415
1416 CALL cp_cfm_power(cfm_matrix_m_tmp, threshold=eps_eigval_s, exponent=-1.0_dp)
1417
1418 ! save L(k) = SQRT(V(k))
1419 CALL cp_cfm_power(cfm_matrix_l_kpoints(ikp), threshold=0.0_dp, exponent=0.5_dp)
1420
1421 ! get K = M^(-1)*L, where L is the Cholesky decomposition of V = L^T*L
1422 CALL cp_cfm_create(cfm_matrix_minv_l_kpoints(ikp), matrix_struct)
1423 CALL parallel_gemm("N", "C", dimen_ri, dimen_ri, dimen_ri, z_one, cfm_matrix_m_tmp, &
1424 cfm_matrix_l_kpoints(ikp), &
1425 z_zero, cfm_matrix_minv_l_kpoints(ikp))
1426
1427 END DO
1428
1429 CALL cp_cfm_release(cfm_matrix_m_tmp)
1430
1431 CALL timestop(handle)
1432
1433 END SUBROUTINE inversion_of_m_and_mult_with_chol_dec_of_v
1434
1435! **************************************************************************************************
1436!> \brief ...
1437!> \param fm_matrix_Minv_Vtrunc_Minv ...
1438!> \param fm_matrix_Minv ...
1439!> \param qs_env ...
1440! **************************************************************************************************
1441 SUBROUTINE gamma_only_inversion_of_m_and_mult_with_chol_dec_of_vtrunc(fm_matrix_Minv_Vtrunc_Minv, &
1442 fm_matrix_Minv, qs_env)
1443
1444 TYPE(cp_fm_type), DIMENSION(:, :) :: fm_matrix_minv_vtrunc_minv, &
1445 fm_matrix_minv
1446 TYPE(qs_environment_type), POINTER :: qs_env
1447
1448 CHARACTER(LEN=*), PARAMETER :: &
1449 routinen = 'Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc'
1450 REAL(kind=dp), PARAMETER :: eps_eigval_s_gamma_thre = 1.0e-18_dp
1451
1452 INTEGER :: dimen_ri, handle, ndep
1453 REAL(kind=dp) :: eps_eigval_s_gamma
1454 TYPE(cp_fm_type) :: fm_matrix_ri_metric_inv_work, fm_work
1455
1456 CALL timeset(routinen, handle)
1457
1458 CALL cp_fm_create(fm_work, fm_matrix_minv(1, 1)%matrix_struct)
1459 CALL cp_fm_set_all(fm_work, 0.0_dp)
1460
1461 CALL cp_fm_create(fm_matrix_ri_metric_inv_work, fm_matrix_minv(1, 1)%matrix_struct)
1462 CALL cp_fm_set_all(fm_matrix_ri_metric_inv_work, 0.0_dp)
1463
1464 CALL cp_fm_get_info(matrix=fm_matrix_minv(1, 1), nrow_global=dimen_ri)
1465
1466 eps_eigval_s_gamma = qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S_Gamma
1467
1468 IF (eps_eigval_s_gamma > eps_eigval_s_gamma_thre) THEN
1469
1470 ! remove small eigenvalues
1471 CALL cp_fm_power(fm_matrix_minv(1, 1), fm_matrix_ri_metric_inv_work, -0.5_dp, &
1472 eps_eigval_s_gamma, ndep)
1473
1474 ELSE
1475
1476 CALL cholesky_decomp(fm_matrix_minv(1, 1), dimen_ri, do_inversion=.true.)
1477
1478 CALL invert_mat(fm_matrix_minv(1, 1))
1479
1480 END IF
1481
1482 CALL parallel_gemm('N', 'T', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_matrix_minv(1, 1), &
1483 fm_matrix_minv(1, 1), 0.0_dp, fm_matrix_ri_metric_inv_work)
1484
1485 CALL cp_fm_to_fm(fm_matrix_ri_metric_inv_work, fm_matrix_minv(1, 1))
1486
1487 CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_matrix_ri_metric_inv_work, &
1488 fm_matrix_minv_vtrunc_minv(1, 1), 0.0_dp, fm_work)
1489
1490 CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_work, &
1491 fm_matrix_ri_metric_inv_work, 0.0_dp, fm_matrix_minv_vtrunc_minv(1, 1))
1492
1493 CALL cp_fm_release(fm_work)
1494 CALL cp_fm_release(fm_matrix_ri_metric_inv_work)
1495
1496 CALL timestop(handle)
1497
1498 END SUBROUTINE gamma_only_inversion_of_m_and_mult_with_chol_dec_of_vtrunc
1499
1500! **************************************************************************************************
1501!> \brief ...
1502!> \param qs_env ...
1503!> \param trunc_coulomb ...
1504!> \param rel_cutoff_trunc_coulomb_ri_x ...
1505!> \param cell_grid ...
1506!> \param do_BvK_cell ...
1507! **************************************************************************************************
1508 SUBROUTINE trunc_coulomb_for_exchange(qs_env, trunc_coulomb, rel_cutoff_trunc_coulomb_ri_x, &
1509 cell_grid, do_BvK_cell)
1510 TYPE(qs_environment_type), POINTER :: qs_env
1511 TYPE(libint_potential_type), OPTIONAL :: trunc_coulomb
1512 REAL(kind=dp), OPTIONAL :: rel_cutoff_trunc_coulomb_ri_x
1513 INTEGER, DIMENSION(3), OPTIONAL :: cell_grid
1514 LOGICAL, OPTIONAL :: do_bvk_cell
1515
1516 CHARACTER(LEN=*), PARAMETER :: routinen = 'trunc_coulomb_for_exchange'
1517
1518 INTEGER :: handle, i_dim
1519 INTEGER, DIMENSION(3) :: periodic
1520 LOGICAL :: my_do_bvk_cell
1521 REAL(kind=dp) :: kp_fac, kp_fac_idim, my_rel_cutoff_trunc_coulomb_ri_x, &
1522 shortest_dist_cell_planes
1523 TYPE(cell_type), POINTER :: cell
1524 TYPE(kpoint_type), POINTER :: kpoints_scf
1525
1526 CALL timeset(routinen, handle)
1527
1528 NULLIFY (cell)
1529 CALL get_qs_env(qs_env, cell=cell, kpoints=kpoints_scf)
1530 CALL get_cell(cell=cell, periodic=periodic)
1531
1532 my_do_bvk_cell = .false.
1533 IF (PRESENT(do_bvk_cell)) my_do_bvk_cell = do_bvk_cell
1534 IF (my_do_bvk_cell) THEN
1535 kp_fac = 1.0e10_dp
1536 DO i_dim = 1, 3
1537 ! look for smallest k-point mesh in periodic direction
1538 IF (periodic(i_dim) == 1) THEN
1539 IF (PRESENT(cell_grid)) THEN
1540 kp_fac_idim = real(cell_grid(i_dim), kind=dp)
1541 ELSE
1542 kp_fac_idim = real(kpoints_scf%nkp_grid(i_dim), kind=dp)
1543 END IF
1544 IF (kp_fac > kp_fac_idim) kp_fac = kp_fac_idim
1545 END IF
1546 END DO
1547 ELSE
1548 kp_fac = 1.0_dp
1549 END IF
1550
1551 shortest_dist_cell_planes = 1.0e4_dp
1552 IF (periodic(1) == 1) THEN
1553 IF (shortest_dist_cell_planes > plane_distance(1, 0, 0, cell)) THEN
1554 shortest_dist_cell_planes = plane_distance(1, 0, 0, cell)
1555 END IF
1556 END IF
1557 IF (periodic(2) == 1) THEN
1558 IF (shortest_dist_cell_planes > plane_distance(0, 1, 0, cell)) THEN
1559 shortest_dist_cell_planes = plane_distance(0, 1, 0, cell)
1560 END IF
1561 END IF
1562 IF (periodic(3) == 1) THEN
1563 IF (shortest_dist_cell_planes > plane_distance(0, 0, 1, cell)) THEN
1564 shortest_dist_cell_planes = plane_distance(0, 0, 1, cell)
1565 END IF
1566 END IF
1567
1568 IF (PRESENT(rel_cutoff_trunc_coulomb_ri_x)) THEN
1569 my_rel_cutoff_trunc_coulomb_ri_x = rel_cutoff_trunc_coulomb_ri_x
1570 ELSE
1571 my_rel_cutoff_trunc_coulomb_ri_x = qs_env%mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x
1572 END IF
1573
1574 IF (PRESENT(trunc_coulomb)) THEN
1575 trunc_coulomb%potential_type = do_potential_truncated
1576 trunc_coulomb%cutoff_radius = shortest_dist_cell_planes* &
1577 my_rel_cutoff_trunc_coulomb_ri_x* &
1578 kp_fac
1579 trunc_coulomb%filename = "t_c_g.dat"
1580 ! dummy
1581 trunc_coulomb%omega = 0.0_dp
1582 END IF
1583
1584 CALL timestop(handle)
1585
1586 END SUBROUTINE trunc_coulomb_for_exchange
1587
1588! **************************************************************************************************
1589!> \brief ...
1590!> \param fm_matrix_L ...
1591! **************************************************************************************************
1592 SUBROUTINE invert_mat(fm_matrix_L)
1593
1594 TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_l
1595
1596 CHARACTER(LEN=*), PARAMETER :: routinen = 'invert_mat'
1597
1598 INTEGER :: handle, i_global, iib, j_global, jjb, &
1599 ncol_local, nrow_local
1600 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1601
1602 CALL timeset(routinen, handle)
1603
1604 CALL cp_fm_triangular_invert(matrix_a=fm_matrix_l, uplo_tr='U')
1605
1606 ! clear the lower part of the L^{-1} matrix (just to have no surprises afterwards)
1607 CALL cp_fm_get_info(matrix=fm_matrix_l, &
1608 nrow_local=nrow_local, &
1609 ncol_local=ncol_local, &
1610 row_indices=row_indices, &
1611 col_indices=col_indices)
1612 ! assume upper triangular format
1613 DO jjb = 1, ncol_local
1614 j_global = col_indices(jjb)
1615 DO iib = 1, nrow_local
1616 i_global = row_indices(iib)
1617 IF (j_global < i_global) fm_matrix_l%local_data(iib, jjb) = 0.0_dp
1618 END DO
1619 END DO
1620
1621 CALL timestop(handle)
1622
1623 END SUBROUTINE invert_mat
1624
1625! **************************************************************************************************
1626!> \brief ...
1627!> \param fm_matrix_L ...
1628!> \param blacs_env_L ...
1629!> \param dimen_RI ...
1630!> \param para_env_L ...
1631!> \param name ...
1632!> \param fm_matrix_L_struct ...
1633! **************************************************************************************************
1634 SUBROUTINE create_matrix_l(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, name, fm_matrix_L_struct)
1635 TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_l
1636 TYPE(cp_blacs_env_type), POINTER :: blacs_env_l
1637 INTEGER, INTENT(IN) :: dimen_ri
1638 TYPE(mp_para_env_type), POINTER :: para_env_l
1639 CHARACTER(LEN=*), INTENT(IN) :: name
1640 TYPE(cp_fm_struct_type), OPTIONAL, POINTER :: fm_matrix_l_struct
1641
1642 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_matrix_L'
1643
1644 INTEGER :: handle
1645 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1646
1647 CALL timeset(routinen, handle)
1648
1649 ! create the full matrix L defined in the L group
1650 IF (PRESENT(fm_matrix_l_struct)) THEN
1651 CALL cp_fm_create(fm_matrix_l, fm_matrix_l_struct, name=name)
1652 ELSE
1653 NULLIFY (fm_struct)
1654 CALL cp_fm_struct_create(fm_struct, context=blacs_env_l, nrow_global=dimen_ri, &
1655 ncol_global=dimen_ri, para_env=para_env_l)
1656
1657 CALL cp_fm_create(fm_matrix_l, fm_struct, name=name)
1658
1659 CALL cp_fm_struct_release(fm_struct)
1660 END IF
1661
1662 CALL cp_fm_set_all(matrix=fm_matrix_l, alpha=0.0_dp)
1663
1664 CALL timestop(handle)
1665
1666 END SUBROUTINE create_matrix_l
1667
1668! **************************************************************************************************
1669!> \brief ...
1670!> \param fm_matrix_L ...
1671!> \param my_group_L_start ...
1672!> \param my_group_L_end ...
1673!> \param my_group_L_size ...
1674!> \param my_Lrows ...
1675! **************************************************************************************************
1676 SUBROUTINE grep_lcols(fm_matrix_L, &
1677 my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
1678 TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_l
1679 INTEGER, INTENT(IN) :: my_group_l_start, my_group_l_end, &
1680 my_group_l_size
1681 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
1682 INTENT(OUT) :: my_lrows
1683
1684 CHARACTER(LEN=*), PARAMETER :: routinen = 'grep_Lcols'
1685
1686 INTEGER :: dimen_ri, handle, handle2, i_global, iib, j_global, jjb, max_row_col_local, &
1687 ncol_local, ncol_rec, nrow_local, nrow_rec, proc_receive_static, proc_send_static, &
1688 proc_shift
1689 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, rec_col_row_info
1690 INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_rec, &
1691 row_indices, row_indices_rec
1692 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: local_l, rec_l
1693 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1694 POINTER :: local_l_internal
1695 TYPE(mp_para_env_type), POINTER :: para_env
1696
1697 CALL timeset(routinen, handle)
1698
1699 CALL cp_fm_get_info(matrix=fm_matrix_l, &
1700 nrow_local=nrow_local, &
1701 ncol_local=ncol_local, &
1702 row_indices=row_indices, &
1703 col_indices=col_indices, &
1704 nrow_global=dimen_ri, &
1705 local_data=local_l_internal, &
1706 para_env=para_env)
1707
1708 ALLOCATE (my_lrows(dimen_ri, my_group_l_size))
1709 my_lrows = 0.0_dp
1710
1711 ALLOCATE (local_l(nrow_local, ncol_local))
1712 local_l(:, :) = local_l_internal(1:nrow_local, 1:ncol_local)
1713
1714 max_row_col_local = max(nrow_local, ncol_local)
1715 CALL para_env%max(max_row_col_local)
1716
1717 ALLOCATE (local_col_row_info(0:max_row_col_local, 2))
1718 local_col_row_info = 0
1719 ! 0,1 nrows
1720 local_col_row_info(0, 1) = nrow_local
1721 local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local)
1722 ! 0,2 ncols
1723 local_col_row_info(0, 2) = ncol_local
1724 local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local)
1725
1726 ALLOCATE (rec_col_row_info(0:max_row_col_local, 2))
1727
1728 ! accumulate data on my_Lrows starting from myself
1729 DO jjb = 1, ncol_local
1730 j_global = col_indices(jjb)
1731 IF (j_global >= my_group_l_start .AND. j_global <= my_group_l_end) THEN
1732 DO iib = 1, nrow_local
1733 i_global = row_indices(iib)
1734 my_lrows(i_global, j_global - my_group_l_start + 1) = local_l(iib, jjb)
1735 END DO
1736 END IF
1737 END DO
1738
1739 proc_send_static = modulo(para_env%mepos + 1, para_env%num_pe)
1740 proc_receive_static = modulo(para_env%mepos - 1, para_env%num_pe)
1741
1742 CALL timeset(routinen//"_comm", handle2)
1743
1744 DO proc_shift = 1, para_env%num_pe - 1
1745 ! first exchange information on the local data
1746 rec_col_row_info = 0
1747 CALL para_env%sendrecv(local_col_row_info, proc_send_static, rec_col_row_info, proc_receive_static)
1748 nrow_rec = rec_col_row_info(0, 1)
1749 ncol_rec = rec_col_row_info(0, 2)
1750
1751 ALLOCATE (row_indices_rec(nrow_rec))
1752 row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
1753
1754 ALLOCATE (col_indices_rec(ncol_rec))
1755 col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
1756
1757 ALLOCATE (rec_l(nrow_rec, ncol_rec))
1758 rec_l = 0.0_dp
1759
1760 ! then send and receive the real data
1761 CALL para_env%sendrecv(local_l, proc_send_static, rec_l, proc_receive_static)
1762
1763 ! accumulate the received data on my_Lrows
1764 DO jjb = 1, ncol_rec
1765 j_global = col_indices_rec(jjb)
1766 IF (j_global >= my_group_l_start .AND. j_global <= my_group_l_end) THEN
1767 DO iib = 1, nrow_rec
1768 i_global = row_indices_rec(iib)
1769 my_lrows(i_global, j_global - my_group_l_start + 1) = rec_l(iib, jjb)
1770 END DO
1771 END IF
1772 END DO
1773
1774 local_col_row_info(:, :) = rec_col_row_info
1775 CALL move_alloc(rec_l, local_l)
1776
1777 DEALLOCATE (col_indices_rec)
1778 DEALLOCATE (row_indices_rec)
1779 END DO
1780 CALL timestop(handle2)
1781
1782 DEALLOCATE (local_col_row_info)
1783 DEALLOCATE (rec_col_row_info)
1784 DEALLOCATE (local_l)
1785
1786 CALL timestop(handle)
1787
1788 END SUBROUTINE grep_lcols
1789
1790END 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:126
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:95
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:103
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
group_dist_types::create_group_dist
Definition group_dist_types.F:39
group_dist_types::get_group_dist
Definition group_dist_types.F:47
group_dist_types::release_group_dist
Definition group_dist_types.F:53
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
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:234
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:325
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:274
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:353
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_power
subroutine, public cp_cfm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
Raise the positive spectral part of a Hermitian matrix to a real power. Eigenmodes with eigenvalues b...
Definition cp_cfm_diag.F:239
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_copy_general
subroutine, public cp_cfm_copy_general(source, destination, para_env)
General copy of a cfm matrix to another cfm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition cp_cfm_types.F:905
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:845
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
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:793
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:643
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_reserve_all_blocks
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
Definition cp_dbcsr_contrib.F:300
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::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:233
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:696
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:1691
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:1969
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:62
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:1144
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:262
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:814
cp_fm_redistribution
Redistribution of full matrices and distributed arrays.
Definition cp_fm_redistribution.F:11
cp_fm_redistribution::cp_array_columns_to_fm
subroutine, public cp_array_columns_to_fm(matrix, array, para_env, first_col, last_col, nrows)
Redistribute a column-distributed array into an existing full matrix.
Definition cp_fm_redistribution.F:354
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:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
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:1569
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:1045
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, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:166
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:18
input_constants::do_eri_os
integer, parameter, public do_eri_os
Definition input_constants.F:1338
input_constants::do_eri_mme
integer, parameter, public do_eri_mme
Definition input_constants.F:1338
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:856
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:856
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:856
input_constants::do_potential_long
integer, parameter, public do_potential_long
Definition input_constants.F:856
input_constants::do_eri_gpw
integer, parameter, public do_eri_gpw
Definition input_constants.F:1338
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:76
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:1138
kpoint_methods::kpoint_init_cell_index
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
Definition kpoint_methods.F:1004
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:585
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 Coulomb operator types.
Definition libint_2c_3c.F:886
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:127
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:143
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
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:2163
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:167
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:119
mp2_ri_2c
Framework for 2c-integrals for RI.
Definition mp2_ri_2c.F:14
mp2_ri_2c::ri_2c_integral_mat
subroutine, public ri_2c_integral_mat(qs_env, fm_matrix_minv_l_gamma, fm_matrix_l_struct, dimen_ri, ri_metric, put_mat_ks_env, regularization_ri)
...
Definition mp2_ri_2c.F:576
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, eps_svd, potential, ri_metric, cfm_matrix_l_kpoints, cfm_matrix_minv_l_kpoints, fm_matrix_minv_l_gamma, 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:169
mp2_ri_2c::trunc_coulomb_for_exchange
subroutine, public trunc_coulomb_for_exchange(qs_env, trunc_coulomb, rel_cutoff_trunc_coulomb_ri_x, cell_grid, do_bvk_cell)
...
Definition mp2_ri_2c.F:1510
mp2_ri_2c::ri_2c_integral_mat_kp
subroutine, public ri_2c_integral_mat_kp(qs_env, cfm_matrix_m_kpoints, fm_matrix_struct, dimen_ri, ri_metric, kpoints, put_mat_ks_env, regularization_ri, ikp_ext, do_build_cell_index)
Computes the complex k-point RI metric as native CFM matrices.
Definition mp2_ri_2c.F:607
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(cfm_matrix_m_kpoints, cfm_matrix_l_kpoints, cfm_matrix_minv_l_kpoints, dimen_ri, kpoints, eps_eigval_s)
...
Definition mp2_ri_2c.F:1388
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, ncgf)
Get the components of a particle set.
Definition particle_methods.F:123
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, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
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:77
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, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, 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_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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, 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_vhxc, 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, xcint_weights, 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_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:591
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, 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_vhxc, 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, xcint_weights, 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_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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:343
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:1096
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors.F:11
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:3156
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:144
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_eri_mme_interface::cp_eri_mme_param
Definition cp_eri_mme_interface.F:79
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_types.F:62
group_dist_types::group_dist_d1_type
Definition group_dist_types.F:35
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
message_passing::mp_comm_type
Definition message_passing.F:179
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67