1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines to calculate and distribute 2c- and 3c- integrals for RI 8!> \par History 9!> 06.2012 created [Mauro Del Ben] 10!> 03.2019 separated from mp2_ri_gpw [Frederick Stein] 11! ************************************************************************************************** 12MODULE mp2_integrals 13 USE ai_contraction_sphi, ONLY: ab_contract,& 14 abc_contract 15 USE ai_overlap3, ONLY: overlap3 16 USE atomic_kind_types, ONLY: atomic_kind_type,& 17 get_atomic_kind_set 18 USE basis_set_types, ONLY: gto_basis_set_p_type,& 19 gto_basis_set_type 20 USE bibliography, ONLY: DelBen2013,& 21 cite_reference 22 USE cell_types, ONLY: cell_type,& 23 get_cell 24 USE cp_blacs_env, ONLY: cp_blacs_env_type 25 USE cp_control_types, ONLY: dft_control_type 26 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 27 cp_dbcsr_m_by_n_from_template,& 28 dbcsr_allocate_matrix_set,& 29 dbcsr_deallocate_matrix_set 30 USE cp_eri_mme_interface, ONLY: cp_eri_mme_param,& 31 cp_eri_mme_set_params 32 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 33 cp_fm_struct_release,& 34 cp_fm_struct_type 35 USE cp_fm_types, ONLY: cp_fm_create,& 36 cp_fm_get_info,& 37 cp_fm_p_type,& 38 cp_fm_release,& 39 cp_fm_type 40 USE cp_para_env, ONLY: cp_para_env_release 41 USE cp_para_types, ONLY: cp_para_env_type 42 USE dbcsr_api, ONLY: & 43 dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, & 44 dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, & 45 dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, & 46 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, & 47 dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_release_p, dbcsr_reserve_blocks, & 48 dbcsr_scalar, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_real_8 49 USE dbcsr_tensor_api, ONLY: & 50 dbcsr_t_contract, dbcsr_t_copy, dbcsr_t_create, dbcsr_t_destroy, & 51 dbcsr_t_distribution_destroy, dbcsr_t_distribution_new, dbcsr_t_distribution_type, & 52 dbcsr_t_filter, dbcsr_t_get_block, dbcsr_t_get_stored_coordinates, & 53 dbcsr_t_mp_environ_pgrid, dbcsr_t_pgrid_create, dbcsr_t_pgrid_destroy, dbcsr_t_pgrid_type, & 54 dbcsr_t_put_block, dbcsr_t_reserve_blocks, dbcsr_t_split_blocks, dbcsr_t_type 55 USE group_dist_types, ONLY: create_group_dist,& 56 get_group_dist,& 57 group_dist_d1_type,& 58 release_group_dist 59 USE input_constants, ONLY: do_eri_gpw,& 60 do_eri_mme,& 61 do_eri_os,& 62 do_potential_coulomb,& 63 do_potential_id,& 64 do_potential_long,& 65 do_potential_short,& 66 do_potential_truncated 67 USE kinds, ONLY: dp 68 USE kpoint_methods, ONLY: kpoint_init_cell_index 69 USE kpoint_types, ONLY: get_kpoint_info,& 70 kpoint_type 71 USE libint_2c_3c, ONLY: libint_potential_type 72 USE machine, ONLY: m_flush 73 USE message_passing, ONLY: mp_alltoall,& 74 mp_cart_create,& 75 mp_dims_create,& 76 mp_environ,& 77 mp_max,& 78 mp_sendrecv,& 79 mp_sync 80 USE molecule_kind_types, ONLY: molecule_kind_type 81 USE mp2_eri, ONLY: mp2_eri_3c_integrate 82 USE mp2_eri_gpw, ONLY: cleanup_gpw,& 83 mp2_eri_3c_integrate_gpw,& 84 prepare_gpw 85 USE mp2_ri_2c, ONLY: get_2c_integrals 86 USE mp2_types, ONLY: integ_mat_buffer_type,& 87 one_dim_int_array,& 88 two_dim_real_array 89 USE orbital_pointers, ONLY: ncoset 90 USE particle_methods, ONLY: get_particle_set 91 USE particle_types, ONLY: particle_type 92 USE pw_env_types, ONLY: pw_env_type 93 USE pw_poisson_types, ONLY: pw_poisson_type 94 USE pw_pool_types, ONLY: pw_pool_type 95 USE pw_types, ONLY: pw_p_type 96 USE qs_environment_types, ONLY: get_qs_env,& 97 qs_environment_type,& 98 set_qs_env 99 USE qs_integral_utils, ONLY: basis_set_list_setup 100 USE qs_kind_types, ONLY: get_qs_kind,& 101 qs_kind_type 102 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 103 neighbor_list_iterate,& 104 neighbor_list_iterator_create,& 105 neighbor_list_iterator_p_type,& 106 neighbor_list_iterator_release,& 107 neighbor_list_set_p_type 108 USE qs_tensors, ONLY: build_3c_integrals,& 109 build_3c_neighbor_lists,& 110 contiguous_tensor_dist,& 111 neighbor_list_3c_destroy 112 USE qs_tensors_types, ONLY: cyclic_tensor_dist,& 113 distribution_3d_create,& 114 distribution_3d_type,& 115 neighbor_list_3c_type,& 116 split_block_sizes 117 USE rpa_communication, ONLY: communicate_buffer 118 USE task_list_types, ONLY: task_list_type 119 USE util, ONLY: get_limit 120 121!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num 122#include "./base/base_uses.f90" 123 124 IMPLICIT NONE 125 126 PRIVATE 127 128 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_integrals' 129 130 PUBLIC :: mp2_ri_gpw_compute_in 131 132CONTAINS 133 134! ************************************************************************************************** 135!> \brief with ri mp2 gpw 136!> \param BIb_C ... 137!> \param BIb_C_gw ... 138!> \param BIb_C_bse_ij ... 139!> \param BIb_C_bse_ab ... 140!> \param gd_array ... 141!> \param gd_B_virtual ... 142!> \param dimen_RI ... 143!> \param dimen_RI_red ... 144!> \param qs_env ... 145!> \param para_env ... 146!> \param para_env_sub ... 147!> \param color_sub ... 148!> \param cell ... 149!> \param particle_set ... 150!> \param atomic_kind_set ... 151!> \param qs_kind_set ... 152!> \param mo_coeff ... 153!> \param fm_matrix_L_RI_metric ... 154!> \param nmo ... 155!> \param homo ... 156!> \param mat_munu ... 157!> \param mat_munu_mao_occ_virt ... 158!> \param mat_munu_mao_virt_occ ... 159!> \param sab_orb_sub ... 160!> \param sab_orb_all ... 161!> \param mo_coeff_o ... 162!> \param mo_coeff_v ... 163!> \param mo_coeff_all ... 164!> \param mo_coeff_gw ... 165!> \param eps_filter ... 166!> \param unit_nr ... 167!> \param mp2_memory ... 168!> \param calc_PQ_cond_num ... 169!> \param calc_forces ... 170!> \param blacs_env_sub ... 171!> \param my_do_gw ... 172!> \param do_bse ... 173!> \param gd_B_all ... 174!> \param starts_array_mc_t ... 175!> \param ends_array_mc_t ... 176!> \param gw_corr_lev_occ ... 177!> \param gw_corr_lev_virt ... 178!> \param do_im_time ... 179!> \param do_mao ... 180!> \param do_kpoints_cubic_RPA ... 181!> \param kpoints ... 182!> \param mat_3c_overl_int ... 183!> \param do_dbcsr_t ... 184!> \param t_3c_overl_int ... 185!> \param t_3c_M ... 186!> \param t_3c_O ... 187!> \param mat_3c_overl_int_mao_for_occ ... 188!> \param mat_3c_overl_int_mao_for_virt ... 189!> \param mao_coeff_occ ... 190!> \param mao_coeff_virt ... 191!> \param ri_metric ... 192!> \param gd_B_occ_bse ... 193!> \param gd_B_virt_bse ... 194!> \param BIb_C_beta ... 195!> \param BIb_C_gw_beta ... 196!> \param gd_B_virtual_beta ... 197!> \param homo_beta ... 198!> \param mo_coeff_o_beta ... 199!> \param mo_coeff_v_beta ... 200!> \param mo_coeff_all_beta ... 201!> \param mo_coeff_gw_beta ... 202!> \author Mauro Del Ben 203! ************************************************************************************************** 204 SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd_array, gd_B_virtual, & 205 dimen_RI, dimen_RI_red, qs_env, para_env, para_env_sub, color_sub, & 206 cell, particle_set, atomic_kind_set, qs_kind_set, mo_coeff, & 207 fm_matrix_L_RI_metric, nmo, homo, & 208 mat_munu, mat_munu_mao_occ_virt, mat_munu_mao_virt_occ, & 209 sab_orb_sub, sab_orb_all, mo_coeff_o, mo_coeff_v, mo_coeff_all, & 210 mo_coeff_gw, eps_filter, unit_nr, & 211 mp2_memory, calc_PQ_cond_num, calc_forces, blacs_env_sub, my_do_gw, do_bse, & 212 gd_B_all, starts_array_mc_t, ends_array_mc_t, & 213 gw_corr_lev_occ, gw_corr_lev_virt, & 214 do_im_time, do_mao, do_kpoints_cubic_RPA, kpoints, & 215 mat_3c_overl_int, do_dbcsr_t, t_3c_overl_int, t_3c_M, t_3c_O, mat_3c_overl_int_mao_for_occ, & 216 mat_3c_overl_int_mao_for_virt, mao_coeff_occ, mao_coeff_virt, & 217 ri_metric, gd_B_occ_bse, gd_B_virt_bse, BIb_C_beta, BIb_C_gw_beta, & 218 gd_B_virtual_beta, homo_beta, mo_coeff_o_beta, mo_coeff_v_beta, & 219 mo_coeff_all_beta, mo_coeff_gw_beta) 220 221 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & 222 INTENT(OUT) :: BIb_C, BIb_C_gw, BIb_C_bse_ij, & 223 BIb_C_bse_ab 224 TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array, gd_B_virtual 225 INTEGER, INTENT(OUT) :: dimen_RI, dimen_RI_red 226 TYPE(qs_environment_type), POINTER :: qs_env 227 TYPE(cp_para_env_type), POINTER :: para_env, para_env_sub 228 INTEGER, INTENT(IN) :: color_sub 229 TYPE(cell_type), POINTER :: cell 230 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 231 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 232 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 233 TYPE(cp_fm_type), POINTER :: mo_coeff 234 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_matrix_L_RI_metric 235 INTEGER, INTENT(IN) :: nmo, homo 236 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu, mat_munu_mao_occ_virt, & 237 mat_munu_mao_virt_occ 238 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 239 INTENT(IN), POINTER :: sab_orb_sub, sab_orb_all 240 TYPE(dbcsr_type), POINTER :: mo_coeff_o, mo_coeff_v, mo_coeff_all, & 241 mo_coeff_gw 242 REAL(KIND=dp), INTENT(IN) :: eps_filter 243 INTEGER, INTENT(IN) :: unit_nr 244 REAL(KIND=dp), INTENT(IN) :: mp2_memory 245 LOGICAL, INTENT(IN) :: calc_PQ_cond_num, calc_forces 246 TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub 247 LOGICAL, INTENT(IN) :: my_do_gw, do_bse 248 TYPE(group_dist_d1_type), INTENT(OUT) :: gd_B_all 249 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: starts_array_mc_t, ends_array_mc_t 250 INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt 251 LOGICAL, INTENT(IN) :: do_im_time, do_mao, do_kpoints_cubic_RPA 252 TYPE(kpoint_type), POINTER :: kpoints 253 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int 254 LOGICAL, INTENT(IN) :: do_dbcsr_t 255 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_overl_int 256 TYPE(dbcsr_t_type), INTENT(OUT) :: t_3c_M 257 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :), & 258 INTENT(OUT) :: t_3c_O 259 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_mao_for_occ, & 260 mat_3c_overl_int_mao_for_virt 261 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coeff_occ, mao_coeff_virt 262 TYPE(libint_potential_type), INTENT(IN) :: ri_metric 263 TYPE(group_dist_d1_type), INTENT(OUT) :: gd_B_occ_bse, gd_B_virt_bse 264 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & 265 INTENT(OUT), OPTIONAL :: BIb_C_beta, BIb_C_gw_beta 266 TYPE(group_dist_d1_type), INTENT(OUT), OPTIONAL :: gd_B_virtual_beta 267 INTEGER, INTENT(IN), OPTIONAL :: homo_beta 268 TYPE(dbcsr_type), OPTIONAL, POINTER :: mo_coeff_o_beta, mo_coeff_v_beta, & 269 mo_coeff_all_beta, mo_coeff_gw_beta 270 271 CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_ri_gpw_compute_in', & 272 routineP = moduleN//':'//routineN 273 274 INTEGER :: color_sub_3c, cut_memory, cut_RI, eri_method, gw_corr_lev_total, handle, handle2, & 275 handle4, i, i_counter, itmp(2), j, LLL, max_row_col_local, max_row_col_local_beta, & 276 max_row_col_local_gw, max_row_col_local_occ_bse, max_row_col_local_virt_bse, & 277 my_B_all_end, my_B_all_size, my_B_all_start, my_B_occ_bse_end, my_B_occ_bse_size, & 278 my_B_occ_bse_start, my_B_size, my_B_size_beta, my_B_virt_bse_end, my_B_virt_bse_size, & 279 my_B_virt_bse_start, my_B_virtual_end, my_B_virtual_end_beta, my_B_virtual_start, & 280 my_B_virtual_start_beta, my_group_L_end, my_group_L_size, my_group_L_start, natom, ngroup 281 INTEGER :: ngroup_im_time, ngroup_im_time_P, nimg, nkind, num_diff_blk, num_small_eigen, & 282 potential_type, ri_metric_type, sqrt_max_bsize, virtual, virtual_beta 283 INTEGER, ALLOCATABLE, DIMENSION(:) :: local_atoms_for_mao_basis, my_group_L_ends_im_time, & 284 my_group_L_sizes_im_time, my_group_L_starts_im_time, sizes_AO, sizes_AO_split, sizes_RI, & 285 sizes_RI_split, sub_proc_map 286 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, local_col_row_info_beta, & 287 local_col_row_info_gw, local_col_row_info_occ_bse, local_col_row_info_virt_bse 288 INTEGER, DIMENSION(3) :: pdims_t3c, periodic 289 LOGICAL :: do_alpha_beta, do_gpw, & 290 do_kpoints_from_Gamma, do_svd, & 291 impose_split, memory_info 292 REAL(KIND=dp) :: cond_num, cutoff_old, eps_svd, & 293 mem_for_iaK, omega_pot, & 294 relative_cutoff_old 295 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old 296 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: my_Lrows 297 TYPE(cp_eri_mme_param), POINTER :: eri_param 298 TYPE(cp_fm_type), POINTER :: fm_BIb_bse_ab, fm_BIb_bse_ij, fm_BIb_gw, & 299 fm_BIb_gw_beta, fm_BIb_jb, & 300 fm_BIb_jb_beta, fm_matrix_L 301 TYPE(cp_para_env_type), POINTER :: para_env_L 302 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_munu_local_L 303 TYPE(dbcsr_t_pgrid_type) :: pgrid_t3c_M, pgrid_t3c_overl 304 TYPE(dbcsr_t_type) :: t_3c_overl_int_template 305 TYPE(dbcsr_type) :: matrix_bse_ab, matrix_bse_anu, matrix_bse_ij, matrix_bse_inu, & 306 matrix_ia_jb, matrix_ia_jb_beta, matrix_ia_jnu, matrix_ia_jnu_beta, matrix_in_jm, & 307 matrix_in_jm_beta, matrix_in_jnu, matrix_in_jnu_beta 308 TYPE(dft_control_type), POINTER :: dft_control 309 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ao, basis_set_ri_aux 310 TYPE(neighbor_list_3c_type) :: nl_3c 311 TYPE(pw_env_type), POINTER :: pw_env_sub 312 TYPE(pw_p_type) :: pot_g, psi_L, rho_g, rho_r 313 TYPE(pw_poisson_type), POINTER :: poisson_env 314 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 315 TYPE(task_list_type), POINTER :: task_list_sub 316 TYPE(two_dim_real_array), ALLOCATABLE, & 317 DIMENSION(:) :: mao_coeff_repl_occ, mao_coeff_repl_virt 318 319 CALL timeset(routineN, handle) 320 321 CALL cite_reference(DelBen2013) 322 323 do_alpha_beta = .FALSE. 324 IF (PRESENT(BIb_C_beta) .AND. & 325 PRESENT(gd_B_virtual_beta) .AND. & 326 PRESENT(homo_beta) .AND. & 327 PRESENT(mo_coeff_o_beta) .AND. & 328 PRESENT(mo_coeff_v_beta)) do_alpha_beta = .TRUE. 329 330 IF ((ri_metric%potential_type == do_potential_short .OR. ri_metric%potential_type == do_potential_truncated) & 331 .AND. .NOT. do_im_time) THEN 332 CPABORT("TRUNCATED and SHORTRANGE RI metrics not yet implemented") 333 ENDIF 334 335 virtual = nmo - homo 336 gw_corr_lev_total = gw_corr_lev_virt + gw_corr_lev_occ 337 338 eri_method = qs_env%mp2_env%eri_method 339 eri_param => qs_env%mp2_env%eri_mme_param 340 do_svd = qs_env%mp2_env%do_svd 341 eps_svd = qs_env%mp2_env%eps_svd 342 potential_type = qs_env%mp2_env%potential_parameter%potential_type 343 ri_metric_type = ri_metric%potential_type 344 omega_pot = qs_env%mp2_env%potential_parameter%omega 345 346 ! whether we need gpw integrals (plus pw stuff) 347 do_gpw = (eri_method == do_eri_gpw) .OR. & 348 ((potential_type == do_potential_long .OR. ri_metric_type == do_potential_long) & 349 .AND. qs_env%mp2_env%eri_method == do_eri_os) & 350 .OR. (ri_metric_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme) 351 352 IF (do_svd .AND. calc_forces) THEN 353 CPABORT("SVD not implemented for forces.!") 354 END IF 355 356 do_kpoints_from_Gamma = SUM(qs_env%mp2_env%ri_rpa_im_time%kp_grid) > 0 357 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 358 CALL get_qs_env(qs_env=qs_env, & 359 kpoints=kpoints) 360 END IF 361 IF (do_kpoints_from_Gamma) THEN 362 CALL compute_kpoints(qs_env, kpoints, unit_nr) 363 END IF 364 365 ! in case of imaginary time, we do not need the intermediate matrices 366 IF (.NOT. do_im_time) THEN 367 368 CALL create_intermediate_matrices(matrix_ia_jnu, matrix_ia_jb, mo_coeff_o, virtual, homo, & 369 fm_BIb_jb, "fm_BIb_jb", max_row_col_local, & 370 blacs_env_sub, para_env_sub, local_col_row_info) 371 372 CALL create_group_dist(gd_B_virtual, para_env_sub%num_pe, virtual) 373 CALL get_group_dist(gd_B_virtual, para_env_sub%mepos, my_B_virtual_start, my_B_virtual_end, my_B_size) 374 375 IF (do_alpha_beta) THEN 376 377 virtual_beta = nmo - homo_beta 378 379 CALL create_intermediate_matrices(matrix_ia_jnu_beta, matrix_ia_jb_beta, mo_coeff_o_beta, & 380 virtual_beta, homo_beta, & 381 fm_BIb_jb_beta, "fm_BIb_jb_beta", & 382 max_row_col_local_beta, & 383 blacs_env_sub, para_env_sub, local_col_row_info_beta) 384 385 CALL create_group_dist(gd_B_virtual_beta, para_env_sub%num_pe, virtual_beta) 386 CALL get_group_dist(gd_B_virtual_beta, para_env_sub%mepos, my_B_virtual_start_beta, my_B_virtual_end_beta, & 387 my_B_size_beta) 388 389 END IF 390 391 ! in the case of G0W0, we need (K|nm), n,m may be occ or virt (m restricted to corrected levels) 392 IF (my_do_gw) THEN 393 394 CALL create_intermediate_matrices(matrix_in_jnu, matrix_in_jm, mo_coeff_gw, & 395 nmo, gw_corr_lev_total, & 396 fm_BIb_gw, "fm_BIb_gw", & 397 max_row_col_local_gw, & 398 blacs_env_sub, para_env_sub, local_col_row_info_gw) 399 400 CALL create_group_dist(gd_B_all, para_env_sub%num_pe, nmo) 401 CALL get_group_dist(gd_B_all, para_env_sub%mepos, my_B_all_start, my_B_all_end, my_B_all_size) 402 403 IF (do_alpha_beta) THEN 404 ! deallocate local_col_row_info_gw, otherwise it gets twice allocated in create_intermediate_m 405 DEALLOCATE (local_col_row_info_gw) 406 CALL create_intermediate_matrices(matrix_in_jnu_beta, matrix_in_jm_beta, mo_coeff_gw_beta, & 407 nmo, gw_corr_lev_total, & 408 fm_BIb_gw_beta, "fm_BIb_gw_beta", & 409 max_row_col_local_gw, & 410 blacs_env_sub, para_env_sub, local_col_row_info_gw) 411 412 ! we don"t need parallelization arrays for beta since the matrix sizes of B_nm^P is the same 413 ! for the beta case and therefore the parallelization of beta is the same than for alpha 414 415 END IF 416 END IF 417 END IF ! not for imaginary time 418 419 IF (do_bse) THEN 420 421 CPASSERT(my_do_gw) 422 CPASSERT(.NOT. do_im_time) 423 ! GPW integrals have to be implemented later 424 CPASSERT(.NOT. (eri_method == do_eri_gpw)) 425 426 ! virt x virt matrices 427 CALL create_intermediate_matrices(matrix_bse_anu, matrix_bse_ab, mo_coeff_v, virtual, virtual, & 428 fm_BIb_bse_ab, "fm_BIb_bse_ab", max_row_col_local_virt_bse, & 429 blacs_env_sub, para_env_sub, local_col_row_info_virt_bse) 430 431 CALL create_group_dist(gd_B_virt_bse, para_env_sub%num_pe, virtual) 432 CALL get_group_dist(gd_B_virt_bse, para_env_sub%mepos, my_B_virt_bse_start, my_B_virt_bse_end, my_B_virt_bse_size) 433 434 ! occ x occ matrices 435 CALL create_intermediate_matrices(matrix_bse_inu, matrix_bse_ij, mo_coeff_o, homo, homo, & 436 fm_BIb_bse_ij, "fm_BIb_bse_ij", max_row_col_local_occ_bse, & 437 blacs_env_sub, para_env_sub, local_col_row_info_occ_bse) 438 439 CALL create_group_dist(gd_B_occ_bse, para_env_sub%num_pe, homo) 440 CALL get_group_dist(gd_B_occ_bse, para_env_sub%mepos, my_B_occ_bse_start, my_B_occ_bse_end, my_B_occ_bse_size) 441 442 END IF 443 444 ngroup = para_env%num_pe/para_env_sub%num_pe 445 446 ! Preparations for MME method to compute ERIs 447 IF (qs_env%mp2_env%eri_method .EQ. do_eri_mme) THEN 448 ! cell might have changed, so we need to reset parameters 449 CALL cp_eri_mme_set_params(eri_param, cell, qs_kind_set, basis_type_1="ORB", basis_type_2="RI_AUX", para_env=para_env) 450 ENDIF 451 452 CALL get_cell(cell=cell, periodic=periodic) 453 ! for minimax Ewald summation, full periodicity is required 454 IF (eri_method == do_eri_mme) THEN 455 CPASSERT(periodic(1) == 1 .AND. periodic(2) == 1 .AND. periodic(3) == 1) 456 END IF 457 458 IF (do_svd .AND. (do_kpoints_from_Gamma .OR. do_kpoints_cubic_RPA)) THEN 459 CPABORT("SVD with kpoints not implemented yet!") 460 END IF 461 462 CALL get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, & 463 fm_matrix_L, ngroup, color_sub, dimen_RI, dimen_RI_red, kpoints, mo_coeff, & 464 my_group_L_size, my_group_L_start, my_group_L_end, & 465 gd_array, calc_PQ_cond_num .AND. .NOT. do_svd, cond_num, do_svd, eps_svd, & 466 num_small_eigen, qs_env%mp2_env%potential_parameter, ri_metric, & 467 fm_matrix_L_RI_metric, & 468 do_im_time, do_kpoints_from_Gamma .OR. do_kpoints_cubic_RPA, qs_env%mp2_env%mp2_gpw%eps_pgf_orb_S, & 469 atomic_kind_set, qs_kind_set, sab_orb_sub) 470 471 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 472 "RI_INFO| Cholesky decomposition group size:", para_env_L%num_pe 473 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 474 "RI_INFO| Number of groups for auxiliary basis functions", ngroup 475 IF (calc_PQ_cond_num .OR. do_svd) THEN 476 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T67,ES14.5)") & 477 "RI_INFO| Condition number of the (P|Q):", cond_num 478 IF (do_svd) THEN 479 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,ES8.1,A,T75,i6)") & 480 "RI_INFO| Number of neglected Eigenvalues of (P|Q) smaller than ", eps_svd, ":", num_small_eigen 481 ELSE 482 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,ES14.1,A,T75,i6)") & 483 "RI_INFO| Number of Eigenvalue of (P|Q) smaller ", eps_svd, ":", num_small_eigen 484 END IF 485 END IF 486 487 IF (.NOT. do_im_time) THEN 488 ! replicate the necessary row of the L^{-1} matrix on each proc 489 CALL grep_Lcols(para_env_L, dimen_RI, fm_matrix_L, & 490 my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows) 491 END IF 492 ! clean the L^{-1} matrix 493 CALL cp_fm_release(fm_matrix_L) 494 495 ! in case of imag. time we need the para_env_L later 496 IF (.NOT. do_im_time) THEN 497 CALL cp_para_env_release(para_env_L) 498 END IF 499 500 IF (calc_forces) THEN 501 ! we need (P|Q)^(-1/2) for future use, just save it 502 ! in a fully (home made) distributed way 503 itmp = get_limit(dimen_RI, para_env_sub%num_pe, para_env_sub%mepos) 504 lll = itmp(2) - itmp(1) + 1 505 ALLOCATE (qs_env%mp2_env%ri_grad%PQ_half(lll, my_group_L_size)) 506 qs_env%mp2_env%ri_grad%PQ_half(:, :) = my_Lrows(itmp(1):itmp(2), 1:my_group_L_size) 507 END IF 508 509 IF (unit_nr > 0) THEN 510 WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 511 "RI_INFO| Occupied basis set size:", homo, & 512 "RI_INFO| Virtual basis set size:", virtual, & 513 "RI_INFO| Auxiliary basis set size:", dimen_RI 514 IF (do_svd) THEN 515 WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 516 "RI_INFO| Reduced auxiliary basis set size:", dimen_RI_red 517 END IF 518 519 mem_for_iaK = dimen_RI*REAL(homo, KIND=dp)*virtual*8.0_dp/(1024_dp**2) 520 IF (do_alpha_beta) mem_for_iaK = mem_for_iaK + & 521 dimen_RI*REAL(homo_beta, KIND=dp)*(nmo - homo_beta)*8.0_dp/(1024_dp**2) 522 523 WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for (ia|K) integrals:', & 524 mem_for_iaK, ' MiB' 525 IF (my_do_gw .AND. .NOT. do_im_time) THEN 526 mem_for_iaK = dimen_RI*REAL(nmo, KIND=dp)*gw_corr_lev_total*8.0_dp/(1024_dp**2) 527 528 WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for G0W0-(nm|K) integrals:', & 529 mem_for_iaK, ' MiB' 530 END IF 531 CALL m_flush(unit_nr) 532 ENDIF 533 534 CALL mp_sync(para_env%group) ! sync to see memory output 535 536 ! in case we do imaginary time, we need the overlap matrix (alpha beta P) 537 IF (.NOT. do_im_time) THEN 538 539 ALLOCATE (sub_proc_map(-para_env_sub%num_pe:2*para_env_sub%num_pe - 1)) 540 DO i = 0, para_env_sub%num_pe - 1 541 sub_proc_map(i) = i 542 sub_proc_map(-i - 1) = para_env_sub%num_pe - i - 1 543 sub_proc_map(para_env_sub%num_pe + i) = i 544 END DO 545 546 ! array that will store the (ia|K) integrals 547 ALLOCATE (BIb_C(my_group_L_size, my_B_size, homo)) 548 BIb_C = 0.0_dp 549 550 IF (do_alpha_beta) THEN 551 ALLOCATE (BIb_C_beta(my_group_L_size, my_B_size_beta, homo_beta)) 552 BIb_C_beta = 0.0_dp 553 END IF 554 555 ! in the case of GW, we also need (nm|K) 556 IF (my_do_gw) THEN 557 558 ALLOCATE (BIb_C_gw(my_group_L_size, my_B_all_size, gw_corr_lev_total)) 559 BIb_C_gw = 0.0_dp 560 IF (do_alpha_beta) THEN 561 ALLOCATE (BIb_C_gw_beta(my_group_L_size, my_B_all_size, gw_corr_lev_total)) 562 BIb_C_gw_beta = 0.0_dp 563 END IF 564 565 END IF 566 567 IF (do_bse) THEN 568 569 ALLOCATE (BIb_C_bse_ij(my_group_L_size, my_B_occ_bse_size, homo)) 570 BIb_C_bse_ij = 0.0_dp 571 572 ALLOCATE (BIb_C_bse_ab(my_group_L_size, my_B_virt_bse_size, virtual)) 573 BIb_C_bse_ab = 0.0_dp 574 575 END IF 576 577 CALL timeset(routineN//"_loop", handle2) 578 579 IF (eri_method == do_eri_mme .AND. & 580 (ri_metric%potential_type == do_potential_coulomb .OR. ri_metric%potential_type == do_potential_long) .OR. & 581 eri_method == do_eri_os .AND. ri_metric%potential_type == do_potential_coulomb) THEN 582 583 NULLIFY (mat_munu_local_L) 584 ALLOCATE (mat_munu_local_L(my_group_L_size)) 585 DO LLL = 1, my_group_L_size 586 NULLIFY (mat_munu_local_L(LLL)%matrix) 587 ALLOCATE (mat_munu_local_L(LLL)%matrix) 588 CALL dbcsr_copy(mat_munu_local_L(LLL)%matrix, mat_munu%matrix) 589 CALL dbcsr_set(mat_munu_local_L(LLL)%matrix, 0.0_dp) 590 ENDDO 591 CALL mp2_eri_3c_integrate(eri_param, ri_metric%potential_type, ri_metric%omega, para_env_sub, qs_env, & 592 first_c=my_group_L_start, last_c=my_group_L_end, & 593 mat_ab=mat_munu_local_L, & 594 basis_type_a="ORB", basis_type_b="ORB", & 595 basis_type_c="RI_AUX", & 596 sab_nl=sab_orb_sub, eri_method=eri_method) 597 598 DO LLL = 1, my_group_L_size 599 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_ia_jnu, matrix_ia_jb, & 600 fm_BIb_jb, BIb_C(LLL, 1:my_B_size, 1:homo), & 601 mo_coeff_o, mo_coeff_v, eps_filter, max_row_col_local, sub_proc_map, & 602 local_col_row_info, my_B_virtual_end, my_B_virtual_start, "alpha") 603 ENDDO 604 CALL contract_B_L(BIb_C, my_Lrows, gd_B_virtual%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 605 ngroup, color_sub, para_env%group, para_env_sub) 606 607 IF (do_alpha_beta) THEN 608 609 DO LLL = 1, my_group_L_size 610 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_ia_jnu_beta, matrix_ia_jb_beta, & 611 fm_BIb_jb_beta, BIb_C_beta(LLL, 1:my_B_size_beta, 1:homo_beta), & 612 mo_coeff_o_beta, mo_coeff_v_beta, eps_filter, max_row_col_local_beta, sub_proc_map, & 613 local_col_row_info_beta, my_B_virtual_end_beta, my_B_virtual_start_beta, "beta") 614 ENDDO 615 CALL contract_B_L(BIb_C_beta, my_Lrows, gd_B_virtual_beta%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 616 ngroup, color_sub, para_env%group, para_env_sub) 617 618 ENDIF 619 620 IF (my_do_gw) THEN 621 622 DO LLL = 1, my_group_L_size 623 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_in_jnu, matrix_in_jm, & 624 fm_BIb_gw, BIb_C_gw(LLL, 1:my_B_all_size, 1:gw_corr_lev_total), & 625 mo_coeff_gw, mo_coeff_all, eps_filter, max_row_col_local_gw, sub_proc_map, & 626 local_col_row_info_gw, my_B_all_end, my_B_all_start, "gw_alpha") 627 ENDDO 628 CALL contract_B_L(BIb_C_gw, my_Lrows, gd_B_all%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 629 ngroup, color_sub, para_env%group, para_env_sub) 630 631 IF (do_alpha_beta) THEN 632 633 DO LLL = 1, my_group_L_size 634 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_in_jnu_beta, matrix_in_jm_beta, & 635 fm_BIb_gw_beta, BIb_C_gw_beta(LLL, 1:my_B_all_size, 1:gw_corr_lev_total), & 636 mo_coeff_gw_beta, mo_coeff_all_beta, eps_filter, max_row_col_local_gw, & 637 sub_proc_map, local_col_row_info_gw, my_B_all_end, my_B_all_start, "gw_beta") 638 ENDDO 639 CALL contract_B_L(BIb_C_gw_beta, my_Lrows, gd_B_all%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 640 ngroup, color_sub, para_env%group, para_env_sub) 641 642 ENDIF 643 ENDIF 644 645 IF (do_bse) THEN 646 647 ! B^ab_P matrix elements for BSE 648 DO LLL = 1, my_group_L_size 649 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_bse_anu, matrix_bse_ab, & 650 fm_BIb_bse_ab, BIb_C_bse_ab(LLL, 1:my_B_virt_bse_size, 1:virtual), & 651 mo_coeff_v, mo_coeff_v, eps_filter, max_row_col_local_virt_bse, & 652 sub_proc_map, local_col_row_info_virt_bse, my_B_all_end, my_B_all_start, "bse_ab") 653 ENDDO 654 CALL contract_B_L(BIb_C_bse_ab, my_Lrows, gd_B_virt_bse%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 655 ngroup, color_sub, para_env%group, para_env_sub) 656 657 ! B^ij_P matrix elements for BSE 658 DO LLL = 1, my_group_L_size 659 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_bse_inu, matrix_bse_ij, & 660 fm_BIb_bse_ij, BIb_C_bse_ij(LLL, 1:my_B_occ_bse_size, 1:homo), & 661 mo_coeff_o, mo_coeff_o, eps_filter, max_row_col_local_occ_bse, & 662 sub_proc_map, local_col_row_info_occ_bse, & 663 my_B_occ_bse_end, my_B_occ_bse_start, "bse_ij") 664 ENDDO 665 CALL contract_B_L(BIb_C_bse_ij, my_Lrows, gd_B_occ_bse%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 666 ngroup, color_sub, para_env%group, para_env_sub) 667 668 END IF 669 670 DO LLL = 1, my_group_L_size 671 CALL dbcsr_release_p(mat_munu_local_L(LLL)%matrix) 672 ENDDO 673 DEALLOCATE (mat_munu_local_L) 674 675 ELSE IF (do_gpw) THEN 676 677 CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, & 678 auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub) 679 680 DO i_counter = 1, my_group_L_size 681 682 CALL mp2_eri_3c_integrate_gpw(mo_coeff, psi_L, rho_g, atomic_kind_set, qs_kind_set, cell, dft_control, & 683 particle_set, pw_env_sub, my_Lrows(:, i_counter), poisson_env, rho_r, pot_g, & 684 ri_metric%potential_type, ri_metric%omega, mat_munu, qs_env, task_list_sub) 685 686 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, matrix_ia_jnu, matrix_ia_jb, & 687 fm_BIb_jb, BIb_C(i_counter, 1:my_B_size, 1:homo), & 688 mo_coeff_o, mo_coeff_v, eps_filter, max_row_col_local, sub_proc_map, & 689 local_col_row_info, my_B_virtual_end, my_B_virtual_start, "alpha") 690 691 IF (do_alpha_beta) THEN 692 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, matrix_ia_jnu_beta, matrix_ia_jb_beta, & 693 fm_BIb_jb_beta, BIb_C_beta(i_counter, 1:my_B_size_beta, 1:homo_beta), & 694 mo_coeff_o_beta, mo_coeff_v_beta, eps_filter, max_row_col_local_beta, sub_proc_map, & 695 local_col_row_info_beta, my_B_virtual_end_beta, my_B_virtual_start_beta, "beta") 696 697 ENDIF 698 699 IF (my_do_gw) THEN 700 ! transform (K|mu nu) to (K|nm), n corresponds to corrected GW levels, m is in nmo 701 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, matrix_in_jnu, matrix_in_jm, & 702 fm_BIb_gw, BIb_C_gw(i_counter, 1:my_B_all_size, 1:gw_corr_lev_total), & 703 mo_coeff_gw, mo_coeff_all, eps_filter, max_row_col_local_gw, sub_proc_map, & 704 local_col_row_info_gw, my_B_all_end, my_B_all_start, "gw_alpha") 705 706 ! the same for beta 707 IF (do_alpha_beta) THEN 708 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, matrix_in_jnu_beta, matrix_in_jm_beta, & 709 fm_BIb_gw_beta, BIb_C_gw_beta(i_counter, 1:my_B_all_size, 1:gw_corr_lev_total), & 710 mo_coeff_gw_beta, mo_coeff_all_beta, eps_filter, max_row_col_local_gw, & 711 sub_proc_map, local_col_row_info_gw, my_B_all_end, my_B_all_start, "gw_beta") 712 713 END IF 714 END IF 715 716 END DO 717 718 CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, pw_env_sub, & 719 task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L) 720 ELSE 721 CPABORT("Integration method not implemented!") 722 END IF 723 724 CALL timestop(handle2) 725 726 DEALLOCATE (my_Lrows) 727 728 CALL cp_fm_release(fm_BIb_jb) 729 DEALLOCATE (local_col_row_info) 730 731 CALL dbcsr_release(matrix_ia_jnu) 732 CALL dbcsr_release(matrix_ia_jb) 733 IF (do_alpha_beta) THEN 734 CALL dbcsr_release(matrix_ia_jnu_beta) 735 CALL dbcsr_release(matrix_ia_jb_beta) 736 CALL cp_fm_release(fm_BIb_jb_beta) 737 DEALLOCATE (local_col_row_info_beta) 738 END IF 739 740 IF (my_do_gw) THEN 741 CALL dbcsr_release(matrix_in_jnu) 742 CALL dbcsr_release(matrix_in_jm) 743 CALL cp_fm_release(fm_BIb_gw) 744 DEALLOCATE (local_col_row_info_gw) 745 IF (do_alpha_beta) THEN 746 CALL dbcsr_release(matrix_in_jnu_beta) 747 CALL dbcsr_release(matrix_in_jm_beta) 748 CALL cp_fm_release(fm_BIb_gw_beta) 749 END IF 750 END IF 751 752 IF (do_bse) THEN 753 CALL dbcsr_release(matrix_bse_anu) 754 CALL dbcsr_release(matrix_bse_ab) 755 CALL cp_fm_release(fm_BIb_bse_ab) 756 CALL dbcsr_release(matrix_bse_inu) 757 CALL dbcsr_release(matrix_bse_ij) 758 CALL cp_fm_release(fm_BIb_bse_ij) 759 END IF 760 761 DEALLOCATE (sub_proc_map) 762 763 ELSE 764 765 cut_memory = qs_env%mp2_env%ri_rpa_im_time%cut_memory 766 767 ngroup_im_time = para_env%num_pe/qs_env%mp2_env%ri_rpa_im_time%group_size_3c 768 ngroup_im_time_P = para_env%num_pe/qs_env%mp2_env%ri_rpa_im_time%group_size_p 769 impose_split = .NOT. qs_env%mp2_env%ri_rpa_im_time%group_size_internal 770 771 color_sub_3c = para_env%mepos/qs_env%mp2_env%ri_rpa_im_time%group_size_3c 772 773 CALL setup_group_L_im_time(my_group_L_starts_im_time, my_group_L_ends_im_time, my_group_L_sizes_im_time, & 774 dimen_RI, ngroup_im_time, cut_memory, & 775 cut_RI, & 776 color_sub_3c, qs_env) 777 778 memory_info = qs_env%mp2_env%ri_rpa_im_time%memory_info 779 780 IF (do_dbcsr_t) THEN 781 ! we need 3 tensors: 782 ! 1) t_3c_overl_int: 3c overlap integrals, optimized for easy access to integral blocks 783 ! (atomic blocks) 784 ! 2) t_3c_O: 3c overlap integrals, optimized for contraction (split blocks) 785 ! 3) t_3c_M: tensor M, optimized for contraction 786 787 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, dft_control=dft_control) 788 789 IF (.NOT. impose_split) THEN 790 pdims_t3c = 0 791 CALL dbcsr_t_pgrid_create(para_env%group, pdims_t3c, pgrid_t3c_overl, map1_2d=[1, 2], map2_2d=[3]) 792 ELSE 793 pgrid_t3c_overl = get_pgrid_from_ngroup(para_env, ngroup_im_time, map1_2d=[1, 2], map2_2d=[3]) 794 ENDIF 795 796 ! set up basis 797 ALLOCATE (sizes_RI(natom), sizes_AO(natom)) 798 ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind)) 799 CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set) 800 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_ri_aux) 801 CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set) 802 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_ao) 803 804 CALL create_tensor_O_3c(t_3c_overl_int_template, pgrid_t3c_overl, & 805 sizes_AO, sizes_AO, sizes_RI, & 806 ri_metric=ri_metric, nl_3c=nl_3c, qs_env=qs_env, & 807 basis_ao_1=basis_set_ao, basis_ao_2=basis_set_ao, basis_RI=basis_set_ri_aux, & 808 do_kpoints=do_kpoints_cubic_RPA) 809 810 ! init k points 811 IF (do_kpoints_cubic_RPA) THEN 812 ! set up new kpoint type with periodic images according to eps_grid from MP2 section 813 ! instead of eps_pgf_orb from QS section 814 CALL kpoint_init_cell_index(kpoints, nl_3c%jk_list, para_env, dft_control) 815 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 816 "3C_OVERLAP_INTEGRALS_INFO| Number of periodic images considered:", dft_control%nimages 817 818 nimg = dft_control%nimages 819 ELSE 820 nimg = 1 821 ENDIF 822 823 ALLOCATE (t_3c_overl_int(nimg, nimg)) 824 825 DO i = 1, SIZE(t_3c_overl_int, 1) 826 DO j = 1, SIZE(t_3c_overl_int, 2) 827 CALL dbcsr_t_create(t_3c_overl_int_template, t_3c_overl_int(i, j)) 828 ENDDO 829 ENDDO 830 831 CALL dbcsr_t_destroy(t_3c_overl_int_template) 832 833 ! split blocks to improve load balancing for tensor contraction 834 sqrt_max_bsize = qs_env%mp2_env%ri_rpa_im_time%max_bsize_sqrt 835 836 CALL split_block_sizes(sizes_AO, sizes_AO_split, sqrt_max_bsize) 837 CALL split_block_sizes(sizes_RI, sizes_RI_split, sqrt_max_bsize) 838 839 IF (.NOT. impose_split) THEN 840 pdims_t3c = 0 841 CALL dbcsr_t_pgrid_create(para_env%group, pdims_t3c, pgrid_t3c_M, map1_2d=[1], map2_2d=[2, 3]) 842 ELSE 843 pgrid_t3c_M = get_pgrid_from_ngroup(para_env, ngroup_im_time_P, map1_2d=[1], map2_2d=[2, 3]) 844 ENDIF 845 846 CALL create_tensor_M_3c(t_3c_M, pgrid_t3c_M, & 847 cut_memory, & 848 sizes_AO_split, sizes_RI, starts_array_mc_t, ends_array_mc_t) 849 850 CALL dbcsr_t_pgrid_destroy(pgrid_t3c_M) 851 852 ALLOCATE (t_3c_O(SIZE(t_3c_overl_int, 1), SIZE(t_3c_overl_int, 2))) 853 CALL create_tensor_O_3c(t_3c_O(1, 1), pgrid_t3c_overl, sizes_AO_split, & 854 sizes_AO, sizes_RI_split) 855 856 CALL dbcsr_t_pgrid_destroy(pgrid_t3c_overl) 857 DO i = 1, SIZE(t_3c_O, 1) 858 DO j = 1, SIZE(t_3c_O, 2) 859 IF (i > 1 .OR. j > 1) CALL dbcsr_t_create(t_3c_O(1, 1), t_3c_O(i, j)) 860 ENDDO 861 ENDDO 862 863 CALL build_3c_integrals(t_3c_overl_int, & 864 qs_env%mp2_env%mp2_gpw%eps_filter, & 865 qs_env, & 866 nl_3c, & 867 basis_i=basis_set_ri_aux, & 868 basis_j=basis_set_ao, basis_k=basis_set_ao, & 869 potential_parameter=ri_metric, & 870 do_kpoints=do_kpoints_cubic_RPA) 871 872 DEALLOCATE (basis_set_ri_aux, basis_set_ao) 873 ! copy integral tensor t_3c_overl_int to t_3c_O tensor optimized for contraction 874 CALL timeset(routineN//"_copy_3c", handle4) 875 DO i = 1, SIZE(t_3c_overl_int, 1) 876 DO j = 1, SIZE(t_3c_overl_int, 2) 877 878 CALL dbcsr_t_copy(t_3c_overl_int(i, j), t_3c_O(i, j), order=[1, 3, 2], move_data=.FALSE.) 879 CALL dbcsr_t_filter(t_3c_O(i, j), qs_env%mp2_env%mp2_gpw%eps_filter) 880 ENDDO 881 ENDDO 882 CALL timestop(handle4) 883 884 CALL neighbor_list_3c_destroy(nl_3c) 885 ELSE 886 CALL reserve_blocks_3c(mat_3c_overl_int, & 887 mat_munu, qs_env, & 888 dimen_RI, cut_RI, unit_nr, & 889 my_group_L_starts_im_time, my_group_L_ends_im_time, & 890 my_group_L_sizes_im_time, & 891 sab_orb_sub, sab_orb_all, para_env, num_diff_blk) 892 893 ! if we do modified atomic orbitals for the primary basis (one MAO basis for D^occ one for D^virt, then 894 ! we also need different 3-center overlap matrix elements! 895 IF (do_mao) THEN 896 897 natom = SIZE(particle_set) 898 ALLOCATE (local_atoms_for_mao_basis(natom)) 899 local_atoms_for_mao_basis = 0 900 901 CALL reserve_blocks_3c(mat_3c_overl_int_mao_for_occ, mat_munu_mao_occ_virt, qs_env, & 902 dimen_RI, cut_RI, unit_nr, & 903 my_group_L_starts_im_time, my_group_L_ends_im_time, & 904 my_group_L_sizes_im_time, sab_orb_sub, sab_orb_all, para_env, num_diff_blk, & 905 local_atoms_for_mao_basis=local_atoms_for_mao_basis) 906 907 CALL reserve_blocks_3c(mat_3c_overl_int_mao_for_virt, mat_munu_mao_virt_occ, qs_env, & 908 dimen_RI, cut_RI, unit_nr, & 909 my_group_L_starts_im_time, my_group_L_ends_im_time, & 910 my_group_L_sizes_im_time, sab_orb_sub, sab_orb_all, para_env, num_diff_blk, & 911 local_atoms_for_mao_basis=local_atoms_for_mao_basis) 912 913 CALL replicate_mao_coeff(mao_coeff_repl_occ, mao_coeff_occ, local_atoms_for_mao_basis, natom, para_env) 914 915 CALL replicate_mao_coeff(mao_coeff_repl_virt, mao_coeff_virt, local_atoms_for_mao_basis, natom, para_env) 916 917 END IF 918 919 CALL compute_3c_overlap_int(mat_3c_overl_int, & 920 mat_3c_overl_int_mao_for_occ, & 921 mat_3c_overl_int_mao_for_virt, & 922 mat_munu, mat_munu_mao_occ_virt, & 923 qs_env, my_group_L_starts_im_time, & 924 my_group_L_ends_im_time, my_group_L_sizes_im_time, & 925 dimen_RI, cut_RI, sab_orb_sub, sab_orb_all, do_mao, & 926 mao_coeff_repl_occ, mao_coeff_repl_virt, & 927 num_diff_blk) 928 ENDIF 929 930 CALL cp_para_env_release(para_env_L) 931 932 IF (do_mao) THEN 933 CALL clean_up(mao_coeff_repl_occ, mao_coeff_repl_virt, local_atoms_for_mao_basis, natom) 934 END IF 935 936 END IF 937 938 CALL timestop(handle) 939 940 END SUBROUTINE mp2_ri_gpw_compute_in 941 942! ************************************************************************************************** 943!> \brief heuristic to get a global 3d process grid that is consistent with a number of process subgroups 944!> such that balanced (ideally square) 2d grids on subgroups are obtained 945!> \param para_env ... 946!> \param ngroup number of process groups 947!> \param map1_2d mapping of 3d grid to 2 dimensions (see dbcsr_t_pgrid_create) 948!> \param map2_2d mapping of 3d grid to 2 dimensions (see dbcsr_t_pgrid_create) 949!> \return process grid object compatible with DBCSR tensors 950! ************************************************************************************************** 951 FUNCTION get_pgrid_from_ngroup(para_env, ngroup, map1_2d, map2_2d) RESULT(pgrid) 952 TYPE(cp_para_env_type), INTENT(IN) :: para_env 953 INTEGER, INTENT(IN) :: ngroup 954 INTEGER, DIMENSION(:), INTENT(IN) :: map1_2d, map2_2d 955 TYPE(dbcsr_t_pgrid_type) :: pgrid 956 957 INTEGER :: dimsplit, nproc_sub 958 INTEGER, DIMENSION(2) :: pdims_2_2d, pdims_sub_2d 959 INTEGER, DIMENSION(3) :: pdims_t3c 960 961 nproc_sub = para_env%num_pe/ngroup 962 pdims_sub_2d = 0 963 CALL mp_dims_create(nproc_sub, pdims_sub_2d) 964 pdims_2_2d = 0 965 CALL mp_dims_create(pdims_sub_2d(2)*ngroup, pdims_2_2d) 966 IF (SIZE(map1_2d) == 1) THEN 967 dimsplit = 2 968 pdims_t3c(map1_2d(1)) = pdims_sub_2d(1) 969 pdims_t3c(map2_2d) = pdims_2_2d 970 ELSEIF (SIZE(map2_2d) == 1) THEN 971 dimsplit = 1 972 pdims_t3c(map2_2d(1)) = pdims_sub_2d(1) 973 pdims_t3c(map1_2d) = pdims_2_2d 974 ELSE 975 CPABORT("map1_2d and map2_2d not compatible with a 3d grid") 976 ENDIF 977 CALL dbcsr_t_pgrid_create(para_env%group, pdims_t3c, pgrid, map1_2d=map1_2d, map2_2d=map2_2d, & 978 nsplit=ngroup, dimsplit=dimsplit) 979 END FUNCTION 980 981! ************************************************************************************************** 982!> \brief ... 983!> \param t_3c_M ... 984!> \param pgrid_t3c_M ... 985!> \param mem_cut ... 986!> \param sizes_AO ... 987!> \param sizes_RI ... 988!> \param starts_array_mc ... 989!> \param ends_array_mc ... 990! ************************************************************************************************** 991 SUBROUTINE create_tensor_M_3c(t_3c_M, pgrid_t3c_M, mem_cut, sizes_AO, sizes_RI, & 992 starts_array_mc, ends_array_mc) 993 TYPE(dbcsr_t_type), INTENT(OUT) :: t_3c_M 994 TYPE(dbcsr_t_pgrid_type), INTENT(IN) :: pgrid_t3c_M 995 INTEGER, INTENT(IN) :: mem_cut 996 INTEGER, DIMENSION(:), INTENT(IN) :: sizes_AO, sizes_RI 997 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: starts_array_mc, ends_array_mc 998 999 INTEGER :: bsum, imem, size_AO, size_AO_cut, size_RI 1000 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_ao_1, dist_ao_2, dist_RI, & 1001 ends_array_mc_block, & 1002 starts_array_mc_block 1003 INTEGER, DIMENSION(3) :: pcoord, pdims 1004 TYPE(dbcsr_t_distribution_type) :: dist 1005 1006 CALL dbcsr_t_mp_environ_pgrid(pgrid_t3c_M, pdims, pcoord) 1007 1008 size_RI = SIZE(sizes_RI) 1009 size_AO = SIZE(sizes_AO) 1010 1011 ALLOCATE (dist_RI(size_RI)) 1012 CALL cyclic_tensor_dist(size_RI, pdims(1), sizes_RI, dist_RI) 1013 1014 ALLOCATE (starts_array_mc(mem_cut)) 1015 ALLOCATE (ends_array_mc(mem_cut)) 1016 ALLOCATE (starts_array_mc_block(mem_cut)) 1017 ALLOCATE (ends_array_mc_block(mem_cut)) 1018 1019 IF (size_AO < mem_cut) THEN 1020 CPABORT("Use smaller MEMORY_CUT") 1021 ENDIF 1022 1023 CALL contiguous_tensor_dist(size_AO, mem_cut, sizes_AO, limits_start=starts_array_mc_block, limits_end=ends_array_mc_block) 1024 1025 bsum = 0 1026 DO imem = 1, mem_cut 1027 starts_array_mc(imem) = bsum + 1 1028 bsum = bsum + SUM(sizes_AO(starts_array_mc_block(imem):ends_array_mc_block(imem))) 1029 ends_array_mc(imem) = bsum 1030 ENDDO 1031 1032 ALLOCATE (dist_ao_1(size_AO)) 1033 ALLOCATE (dist_ao_2(size_AO)) 1034 DO imem = 1, mem_cut 1035 1036 size_AO_cut = ends_array_mc_block(imem) - starts_array_mc_block(imem) + 1 1037 1038 IF (size_AO_cut < MIN(pdims(2), pdims(3))) THEN 1039 CPABORT("use smaller MEMORY_CUT or use less MPI ranks") 1040 ENDIF 1041 CALL cyclic_tensor_dist(size_AO_cut, pdims(2), sizes_AO(starts_array_mc_block(imem):ends_array_mc_block(imem)), & 1042 dist=dist_ao_1(starts_array_mc_block(imem):ends_array_mc_block(imem))) 1043 CALL cyclic_tensor_dist(size_AO_cut, pdims(3), sizes_AO(starts_array_mc_block(imem):ends_array_mc_block(imem)), & 1044 dist=dist_ao_2(starts_array_mc_block(imem):ends_array_mc_block(imem))) 1045 ENDDO 1046 1047 CALL dbcsr_t_distribution_new(dist, pgrid_t3c_M, [1], [2, 3], dist_RI, dist_ao_1, dist_ao_2) 1048 CALL dbcsr_t_create(t_3c_M, "M (RI | AO AO)", dist, [1], [2, 3], dbcsr_type_real_8, sizes_RI, & 1049 sizes_AO, sizes_AO) 1050 CALL dbcsr_t_distribution_destroy(dist) 1051 1052 END SUBROUTINE create_tensor_M_3c 1053 1054! ************************************************************************************************** 1055!> \brief ... 1056!> \param t_3c_O Create tensor with overlap integrals, optionally create 3c neighbor list matching tensor 1057!> distribution 1058!> \param mp_comm_t3c ... 1059!> \param ao_sizes_1 block sizes of first AO index 1060!> \param ao_sizes_2 block sizes of second AO index 1061!> \param sizes_RI block sizes of RI index 1062!> \param ri_metric ... 1063!> \param nl_3c 3c-neighborlist - if present, the sizes arguments must refer to atomic blocks 1064!> \param qs_env needed for nl_3c 1065!> \param basis_ao_1 needed for nl_3c 1066!> \param basis_ao_2 needed for nl_3c 1067!> \param basis_RI needed for nl_3c 1068!> \param do_kpoints ... 1069! ************************************************************************************************** 1070 SUBROUTINE create_tensor_O_3c(t_3c_O, mp_comm_t3c, ao_sizes_1, ao_sizes_2, sizes_RI, & 1071 ri_metric, nl_3c, qs_env, basis_ao_1, basis_ao_2, basis_RI, & 1072 do_kpoints) 1073 TYPE(dbcsr_t_type), INTENT(OUT) :: t_3c_O 1074 TYPE(dbcsr_t_pgrid_type), INTENT(IN) :: mp_comm_t3c 1075 INTEGER, DIMENSION(:), INTENT(IN) :: ao_sizes_1, ao_sizes_2, sizes_RI 1076 TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: ri_metric 1077 TYPE(neighbor_list_3c_type), INTENT(OUT), OPTIONAL :: nl_3c 1078 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env 1079 TYPE(gto_basis_set_p_type), DIMENSION(:), & 1080 OPTIONAL, POINTER :: basis_ao_1, basis_ao_2, basis_RI 1081 LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints 1082 1083 INTEGER :: mp_comm_t3c_2, nkind, size_AO_1, & 1084 size_AO_2, size_RI 1085 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_ao_1, dist_ao_2, dist_RI 1086 INTEGER, DIMENSION(3) :: pcoor, pcoord, pdims 1087 LOGICAL :: kp 1088 TYPE(dbcsr_t_distribution_type) :: dist 1089 TYPE(distribution_3d_type) :: dist_3d 1090 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1091 1092 IF (PRESENT(do_kpoints)) THEN 1093 kp = do_kpoints 1094 ELSE 1095 kp = .FALSE. 1096 ENDIF 1097 1098 CALL dbcsr_t_mp_environ_pgrid(mp_comm_t3c, pdims, pcoord) 1099 1100 size_AO_1 = SIZE(ao_sizes_1) 1101 size_AO_2 = SIZE(ao_sizes_2) 1102 size_RI = SIZE(sizes_RI) 1103 1104 ALLOCATE (dist_RI(size_RI)) 1105 ALLOCATE (dist_ao_1(size_AO_1)) 1106 ALLOCATE (dist_ao_2(size_AO_2)) 1107 1108 CALL cyclic_tensor_dist(size_RI, pdims(1), sizes_RI, dist_RI) 1109 CALL cyclic_tensor_dist(size_AO_1, pdims(2), ao_sizes_1, dist_ao_1) 1110 CALL cyclic_tensor_dist(size_AO_2, pdims(3), ao_sizes_2, dist_ao_2) 1111 1112 IF (PRESENT(nl_3c)) THEN 1113 CPASSERT(PRESENT(qs_env)) 1114 CPASSERT(PRESENT(ri_metric)) 1115 CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set) 1116 CPASSERT(PRESENT(basis_RI)) 1117 CPASSERT(PRESENT(basis_ao_1)) 1118 CPASSERT(PRESENT(basis_ao_2)) 1119 CALL dbcsr_t_mp_environ_pgrid(mp_comm_t3c, pdims, pcoor) 1120 CALL mp_cart_create(mp_comm_t3c%mp_comm_2d, 3, pdims, pcoord, mp_comm_t3c_2) 1121 CALL distribution_3d_create(dist_3d, dist_RI, dist_ao_1, dist_ao_2, & 1122 nkind, particle_set, mp_comm_t3c_2, own_comm=.TRUE.) 1123 1124 CALL build_3c_neighbor_lists(nl_3c, basis_RI, basis_ao_1, basis_ao_2, & 1125 dist_3d, ri_metric, "RPA_3c_nl", qs_env, & 1126 sym_jk=.NOT. kp, own_dist=.TRUE.) 1127 ENDIF 1128 1129 CALL dbcsr_t_distribution_new(dist, mp_comm_t3c, [1, 2], [3], dist_RI, dist_ao_1, dist_ao_2) 1130 CALL dbcsr_t_create(t_3c_O, "O (RI AO | AO)", dist, [1, 2], [3], dbcsr_type_real_8, sizes_RI, & 1131 ao_sizes_1, ao_sizes_2) 1132 CALL dbcsr_t_distribution_destroy(dist) 1133 1134 END SUBROUTINE create_tensor_O_3c 1135 1136! ************************************************************************************************** 1137!> \brief Contract (P|ai) = (R|P) x (R|ai) 1138!> \param BIb_C (R|ai) 1139!> \param my_Lrows (R|P) 1140!> \param sizes_B number of a (virtual) indices per subgroup process 1141!> \param sizes_L number of P / R (auxiliary) indices per subgroup 1142!> \param blk_size ... 1143!> \param ngroup how many subgroups (NG) 1144!> \param igroup subgroup color 1145!> \param mp_comm communicator 1146!> \param para_env_sub ... 1147! ************************************************************************************************** 1148 SUBROUTINE contract_B_L(BIb_C, my_Lrows, sizes_B, sizes_L, blk_size, ngroup, igroup, mp_comm, para_env_sub) 1149 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: BIb_C 1150 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: my_Lrows 1151 INTEGER, DIMENSION(:), INTENT(IN) :: sizes_B, sizes_L 1152 INTEGER, DIMENSION(2), INTENT(IN) :: blk_size 1153 INTEGER, INTENT(IN) :: ngroup, igroup, mp_comm 1154 TYPE(cp_para_env_type), POINTER :: para_env_sub 1155 1156 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_B_L', routineP = moduleN//':'//routineN 1157 LOGICAL, PARAMETER :: debug = .FALSE. 1158 1159 INTEGER :: check_proc, handle, i, iend, ii, ioff, & 1160 iproc, iproc_glob, istart, loc_a, & 1161 loc_P, nproc, nproc_glob 1162 INTEGER, ALLOCATABLE, DIMENSION(:) :: block_ind_L_P, block_ind_L_R 1163 INTEGER, DIMENSION(1) :: dist_B_i, map_B_1, map_L_1, map_L_2, & 1164 sizes_i 1165 INTEGER, DIMENSION(2) :: map_B_2, pdims_L 1166 INTEGER, DIMENSION(3) :: pdims_B 1167 LOGICAL :: found 1168 INTEGER, DIMENSION(ngroup) :: dist_L_P, dist_L_R 1169 INTEGER, DIMENSION(para_env_sub%num_pe) :: dist_B_a 1170 TYPE(dbcsr_t_distribution_type) :: dist_B, dist_L 1171 TYPE(dbcsr_t_pgrid_type) :: mp_comm_B, mp_comm_L 1172 TYPE(dbcsr_t_type) :: tB_in, tB_in_split, tB_out, & 1173 tB_out_split, tL, tL_split 1174 1175 CALL timeset(routineN, handle) 1176 1177 sizes_i(1) = SIZE(BIb_C, 3) 1178 1179 nproc = para_env_sub%num_pe ! number of processes per subgroup (Nw) 1180 iproc = para_env_sub%mepos ! subgroup-local process ID 1181 1182 ! Total number of processes and global process ID 1183 CALL mp_environ(nproc_glob, iproc_glob, mp_comm) 1184 1185 ! local block index for R/P and a 1186 loc_P = igroup + 1; loc_a = iproc + 1 1187 1188 CPASSERT(SIZE(sizes_L) .EQ. ngroup) 1189 CPASSERT(SIZE(sizes_B) .EQ. nproc) 1190 CPASSERT(sizes_L(loc_P) .EQ. SIZE(BIb_C, 1)) 1191 CPASSERT(sizes_L(loc_P) .EQ. SIZE(my_Lrows, 2)) 1192 CPASSERT(sizes_B(loc_a) .EQ. SIZE(BIb_C, 2)) 1193 1194 ! Tensor distributions as follows: 1195 ! Process grid NG x Nw 1196 ! Each process has coordinates (np, nw) 1197 ! tB_in: (R|ai): R distributed (np), a distributed (nw) 1198 ! tB_out: (P|ai): P distributed (np), a distributed (nw) 1199 ! tL: (R|P): R distributed (nw), P distributed (np) 1200 1201 ! define mappings between tensor index and matrix index: 1202 ! (R|ai) and (P|ai): 1203 map_B_1 = [1] ! index 1 (R or P) maps to 1st matrix index (np distributed) 1204 map_B_2 = [2, 3] ! indices 2, 3 (a, i) map to 2nd matrix index (nw distributed) 1205 ! (R|P): 1206 map_L_1 = [2] ! index 2 (P) maps to 1st matrix index (np distributed) 1207 map_L_2 = [1] ! index 1 (R) maps to 2nd matrix index (nw distributed) 1208 1209 ! derive nd process grid that is compatible with distributions and 2d process grid 1210 ! (R|ai) / (P|ai) on process grid NG x Nw x 1 1211 ! (R|P) on process grid NG x Nw 1212 pdims_B = [ngroup, nproc, 1] 1213 pdims_L = [nproc, ngroup] 1214 1215 CALL dbcsr_t_pgrid_create(mp_comm, pdims_B, mp_comm_B) 1216 CALL dbcsr_t_pgrid_create(mp_comm, pdims_L, mp_comm_L) 1217 1218 ! setup distribution vectors such that distribution matches parallel data layout of BIb_C and my_Lrows 1219 dist_B_i = [0] 1220 dist_B_a = (/(i, i=0, nproc - 1)/) 1221 dist_L_R = (/(MODULO(i, nproc), i=0, ngroup - 1)/) ! R index is replicated in my_Lrows, we impose a cyclic distribution 1222 dist_L_P = (/(i, i=0, ngroup - 1)/) 1223 1224 ! create distributions and tensors 1225 CALL dbcsr_t_distribution_new(dist_B, mp_comm_B, map_B_1, map_B_2, dist_L_P, dist_B_a, dist_B_i) 1226 CALL dbcsr_t_distribution_new(dist_L, mp_comm_L, map_L_1, map_L_2, dist_L_R, dist_L_P) 1227 1228 CALL dbcsr_t_create(tB_in, "(R|ai)", dist_B, map_B_1, map_B_2, dbcsr_type_real_8, sizes_L, sizes_B, sizes_i) 1229 CALL dbcsr_t_create(tB_out, "(P|ai)", dist_B, map_B_1, map_B_2, dbcsr_type_real_8, sizes_L, sizes_B, sizes_i) 1230 CALL dbcsr_t_create(tL, "(R|P)", dist_L, map_L_1, map_L_2, dbcsr_type_real_8, sizes_L, sizes_L) 1231 1232 IF (debug) THEN 1233 ! check that tensor distribution is correct 1234 CALL dbcsr_t_get_stored_coordinates(tB_in, [loc_P, loc_a, 1], check_proc) 1235 CPASSERT(check_proc .EQ. iproc_glob) 1236 ENDIF 1237 1238 ! reserve (R|ai) block 1239 CALL dbcsr_t_reserve_blocks(tB_in, [loc_P], [loc_a], [1]) 1240 1241 ! reserve (R|P) blocks 1242 ! in my_Lrows, R index is replicated. For (R|P), we distribute quadratic blocks cyclically over 1243 ! the processes in a subgroup. 1244 ! There are NG blocks, so each process holds at most NG/Nw+1 blocks. 1245 ALLOCATE (block_ind_L_R(ngroup/nproc + 1)) 1246 ALLOCATE (block_ind_L_P(ngroup/nproc + 1)) 1247 block_ind_L_R(:) = 0; block_ind_L_P(:) = 0 1248 ii = 0 1249 DO i = 1, ngroup 1250 CALL dbcsr_t_get_stored_coordinates(tL, [i, loc_P], check_proc) 1251 IF (check_proc == iproc_glob) THEN 1252 ii = ii + 1 1253 block_ind_L_R(ii) = i 1254 block_ind_L_P(ii) = loc_P 1255 ENDIF 1256 ENDDO 1257 CALL dbcsr_t_reserve_blocks(tL, block_ind_L_R(1:ii), block_ind_L_P(1:ii)) 1258 1259 ! insert (R|ai) block 1260 CALL dbcsr_t_put_block(tB_in, [loc_P, loc_a, 1], SHAPE(BIb_C), BIb_C) 1261 1262 ! insert (R|P) blocks 1263 ioff = 0 1264 DO i = 1, ngroup 1265 istart = ioff + 1; iend = ioff + sizes_L(i) 1266 ioff = ioff + sizes_L(i) 1267 CALL dbcsr_t_get_stored_coordinates(tL, [i, loc_P], check_proc) 1268 IF (check_proc == iproc_glob) THEN 1269 CALL dbcsr_t_put_block(tL, [i, loc_P], [sizes_L(i), sizes_L(loc_P)], my_Lrows(istart:iend, :)) 1270 ENDIF 1271 ENDDO 1272 1273 CALL dbcsr_t_split_blocks(tB_in, tB_in_split, [blk_size(2), blk_size(1), blk_size(1)]) 1274 CALL dbcsr_t_split_blocks(tL, tL_split, [blk_size(2), blk_size(2)]) 1275 CALL dbcsr_t_split_blocks(tB_out, tB_out_split, [blk_size(2), blk_size(1), blk_size(1)]) 1276 1277 ! contract 1278 CALL dbcsr_t_contract(alpha=dbcsr_scalar(1.0_dp), tensor_1=tB_in_split, tensor_2=tL_split, & 1279 beta=dbcsr_scalar(0.0_dp), tensor_3=tB_out_split, & 1280 contract_1=[1], notcontract_1=[2, 3], & 1281 contract_2=[1], notcontract_2=[2], & 1282 map_1=[2, 3], map_2=[1], optimize_dist=.TRUE.) 1283 1284 ! retrieve local block of contraction result (P|ai) 1285 CALL dbcsr_t_copy(tB_out_split, tB_out) 1286 1287 CALL dbcsr_t_get_block(tB_out, [loc_P, loc_a, 1], SHAPE(BIb_C), BIb_C, found) 1288 CPASSERT(found) 1289 1290 ! cleanup 1291 CALL dbcsr_t_destroy(tB_in) 1292 CALL dbcsr_t_destroy(tB_in_split) 1293 CALL dbcsr_t_destroy(tB_out) 1294 CALL dbcsr_t_destroy(tB_out_split) 1295 CALL dbcsr_t_destroy(tL) 1296 CALL dbcsr_t_destroy(tL_split) 1297 1298 CALL dbcsr_t_distribution_destroy(dist_B) 1299 CALL dbcsr_t_distribution_destroy(dist_L) 1300 1301 CALL dbcsr_t_pgrid_destroy(mp_comm_B) 1302 CALL dbcsr_t_pgrid_destroy(mp_comm_L) 1303 1304 CALL timestop(handle) 1305 1306 END SUBROUTINE contract_B_L 1307 1308! ************************************************************************************************** 1309!> \brief Encapsulate building of intermediate matrices matrix_ia_jnu(_beta 1310!> matrix_ia_jb(_beta),fm_BIb_jb(_beta),matrix_in_jnu(for G0W0) and 1311!> fm_BIb_all(for G0W0) 1312!> \param matrix_ia_jnu ... 1313!> \param matrix_ia_jb ... 1314!> \param mo_coeff_templ ... 1315!> \param size_1 ... 1316!> \param size_2 ... 1317!> \param fm_BIb_jb ... 1318!> \param matrix_name_2 ... 1319!> \param max_row_col_local ... 1320!> \param blacs_env_sub ... 1321!> \param para_env_sub ... 1322!> \param local_col_row_info ... 1323!> \author Jan Wilhelm 1324! ************************************************************************************************** 1325 SUBROUTINE create_intermediate_matrices(matrix_ia_jnu, matrix_ia_jb, mo_coeff_templ, size_1, size_2, & 1326 fm_BIb_jb, matrix_name_2, max_row_col_local, & 1327 blacs_env_sub, para_env_sub, local_col_row_info) 1328 1329 TYPE(dbcsr_type), INTENT(OUT) :: matrix_ia_jnu, matrix_ia_jb 1330 TYPE(dbcsr_type), INTENT(INOUT) :: mo_coeff_templ 1331 INTEGER, INTENT(IN) :: size_1, size_2 1332 TYPE(cp_fm_type), POINTER :: fm_BIb_jb 1333 CHARACTER(LEN=*), INTENT(IN) :: matrix_name_2 1334 INTEGER, INTENT(OUT) :: max_row_col_local 1335 TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub 1336 TYPE(cp_para_env_type), POINTER :: para_env_sub 1337 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: local_col_row_info 1338 1339 CHARACTER(LEN=*), PARAMETER :: routineN = 'create_intermediate_matrices', & 1340 routineP = moduleN//':'//routineN 1341 1342 INTEGER :: handle, ncol_local, nfullcols_total, & 1343 nfullrows_total, nrow_local 1344 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 1345 TYPE(cp_fm_struct_type), POINTER :: fm_struct 1346 1347 CALL timeset(routineN, handle) 1348 1349 ! initialize and create the matrix (K|jnu) 1350 CALL dbcsr_create(matrix_ia_jnu, template=mo_coeff_templ) 1351 1352 ! Allocate Sparse matrices: (K|jb) 1353 CALL cp_dbcsr_m_by_n_from_template(matrix_ia_jb, template=mo_coeff_templ, m=size_2, n=size_1, & 1354 sym=dbcsr_type_no_symmetry) 1355 1356 ! set all to zero in such a way that the memory is actually allocated 1357 CALL dbcsr_set(matrix_ia_jnu, 0.0_dp) 1358 CALL dbcsr_set(matrix_ia_jb, 0.0_dp) 1359 1360 ! create the analogous of matrix_ia_jb in fm type 1361 NULLIFY (fm_BIb_jb) 1362 NULLIFY (fm_struct) 1363 CALL dbcsr_get_info(matrix_ia_jb, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total) 1364 CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, nrow_global=nfullrows_total, & 1365 ncol_global=nfullcols_total, para_env=para_env_sub) 1366 CALL cp_fm_create(fm_BIb_jb, fm_struct, name=matrix_name_2) 1367 1368 CALL copy_dbcsr_to_fm(matrix_ia_jb, fm_BIb_jb) 1369 CALL cp_fm_struct_release(fm_struct) 1370 1371 CALL cp_fm_get_info(matrix=fm_BIb_jb, & 1372 nrow_local=nrow_local, & 1373 ncol_local=ncol_local, & 1374 row_indices=row_indices, & 1375 col_indices=col_indices) 1376 1377 max_row_col_local = MAX(nrow_local, ncol_local) 1378 CALL mp_max(max_row_col_local, para_env_sub%group) 1379 1380 ALLOCATE (local_col_row_info(0:max_row_col_local, 2)) 1381 local_col_row_info = 0 1382 ! 0,1 nrows 1383 local_col_row_info(0, 1) = nrow_local 1384 local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local) 1385 ! 0,2 ncols 1386 local_col_row_info(0, 2) = ncol_local 1387 local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local) 1388 1389 CALL timestop(handle) 1390 1391 END SUBROUTINE create_intermediate_matrices 1392 1393! ************************************************************************************************** 1394!> \brief Encapsulate ERI postprocessing: AO to MO transformation and store in B matrix. 1395!> \param para_env ... 1396!> \param mat_munu ... 1397!> \param matrix_ia_jnu ... 1398!> \param matrix_ia_jb ... 1399!> \param fm_BIb_jb ... 1400!> \param BIb_jb ... 1401!> \param mo_coeff_o ... 1402!> \param mo_coeff_v ... 1403!> \param eps_filter ... 1404!> \param max_row_col_local ... 1405!> \param sub_proc_map ... 1406!> \param local_col_row_info ... 1407!> \param my_B_end ... 1408!> \param my_B_start ... 1409!> \param descr ... 1410! ************************************************************************************************** 1411 SUBROUTINE ao_to_mo_and_store_B(para_env, mat_munu, matrix_ia_jnu, matrix_ia_jb, fm_BIb_jb, BIb_jb, & 1412 mo_coeff_o, mo_coeff_v, eps_filter, & 1413 max_row_col_local, sub_proc_map, local_col_row_info, & 1414 my_B_end, my_B_start, descr) 1415 TYPE(cp_para_env_type), POINTER :: para_env 1416 TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu 1417 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ia_jnu, matrix_ia_jb 1418 TYPE(cp_fm_type), POINTER :: fm_BIb_jb 1419 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: BIb_jb 1420 TYPE(dbcsr_type), POINTER :: mo_coeff_o, mo_coeff_v 1421 REAL(KIND=dp), INTENT(IN) :: eps_filter 1422 INTEGER, INTENT(IN) :: max_row_col_local 1423 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: sub_proc_map 1424 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: local_col_row_info 1425 INTEGER, INTENT(IN) :: my_B_end, my_B_start 1426 CHARACTER(LEN=*), INTENT(IN) :: descr 1427 1428 CHARACTER(LEN=*), PARAMETER :: routineN = 'ao_to_mo_and_store_B' 1429 1430 INTEGER :: handle1, handle2 1431 1432 CALL timeset(routineN//"_mult_"//descr, handle1) 1433 1434 CALL dbcsr_multiply("N", "N", 1.0_dp, mat_munu%matrix, mo_coeff_o, & 1435 0.0_dp, matrix_ia_jnu, filter_eps=eps_filter) 1436 CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_ia_jnu, mo_coeff_v, & 1437 0.0_dp, matrix_ia_jb, filter_eps=eps_filter) 1438 CALL timestop(handle1) 1439 1440 CALL timeset(routineN//"_E_Ex_"//descr, handle2) 1441 CALL copy_dbcsr_to_fm(matrix_ia_jb, fm_BIb_jb) 1442 1443 IF (.NOT. (descr .EQ. "bse_ab")) THEN 1444 1445 CALL grep_my_integrals(para_env, fm_BIb_jb, BIb_jb, max_row_col_local, & 1446 sub_proc_map, local_col_row_info, & 1447 my_B_end, my_B_start) 1448 1449 END IF 1450 1451 CALL timestop(handle2) 1452 END SUBROUTINE ao_to_mo_and_store_B 1453 1454! ************************************************************************************************** 1455!> \brief ... 1456!> \param qs_env ... 1457!> \param kpoints ... 1458!> \param unit_nr ... 1459! ************************************************************************************************** 1460 SUBROUTINE compute_kpoints(qs_env, kpoints, unit_nr) 1461 1462 TYPE(qs_environment_type), POINTER :: qs_env 1463 TYPE(kpoint_type), POINTER :: kpoints 1464 INTEGER :: unit_nr 1465 1466 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_kpoints', & 1467 routineP = moduleN//':'//routineN 1468 1469 INTEGER :: handle, i, i_dim, ix, iy, iz, nkp, & 1470 num_dim 1471 INTEGER, DIMENSION(3) :: nkp_grid, periodic 1472 TYPE(cell_type), POINTER :: cell 1473 TYPE(cp_para_env_type), POINTER :: para_env 1474 TYPE(dft_control_type), POINTER :: dft_control 1475 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1476 POINTER :: sab_orb 1477 1478 CALL timeset(routineN, handle) 1479 1480 NULLIFY (cell, dft_control, para_env) 1481 CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env, dft_control=dft_control, sab_orb=sab_orb) 1482 CALL get_cell(cell=cell, periodic=periodic) 1483 1484 ! general because we augment a Monkhorst-Pack mesh by additional points in the BZ 1485 kpoints%kp_scheme = "GENERAL" 1486 kpoints%symmetry = .FALSE. 1487 kpoints%verbose = .FALSE. 1488 kpoints%full_grid = .FALSE. 1489 kpoints%use_real_wfn = .FALSE. 1490 kpoints%eps_geo = 1.e-6_dp 1491 kpoints%full_grid = .TRUE. 1492 nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid(1:3) 1493 kpoints%nkp_grid(1:3) = nkp_grid(1:3) 1494 1495 num_dim = periodic(1) + periodic(2) + periodic(3) 1496 1497 DO i_dim = 1, 3 1498 IF (periodic(i_dim) == 1) THEN 1499 CPASSERT(MODULO(kpoints%nkp_grid(i_dim), 2) == 0) 1500 END IF 1501 END DO 1502 1503 IF (num_dim == 3) THEN 1504 nkp = nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2 1505 ELSE IF (num_dim == 2) THEN 1506 nkp = MAX(nkp_grid(1)*nkp_grid(2)/2, nkp_grid(1)*nkp_grid(3)/2, nkp_grid(2)*nkp_grid(3)/2) 1507 ELSE IF (num_dim == 1) THEN 1508 nkp = MAX(nkp_grid(1)/2, nkp_grid(2)/2, nkp_grid(3)/2) 1509 END IF 1510 1511 kpoints%nkp = nkp 1512 1513 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 1514 "KPOINT_INFO| Number of kpoints for V and W:", nkp 1515 1516 ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp)) 1517 kpoints%wkp(:) = 1.0_dp/REAL(nkp, KIND=dp) 1518 1519 i = 0 1520 DO ix = 1, nkp_grid(1) 1521 DO iy = 1, nkp_grid(2) 1522 DO iz = 1, nkp_grid(3) 1523 1524 i = i + 1 1525 IF (i > nkp) CYCLE 1526 1527 kpoints%xkp(1, i) = REAL(2*ix - nkp_grid(1) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(1), KIND=dp)) 1528 kpoints%xkp(2, i) = REAL(2*iy - nkp_grid(2) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(2), KIND=dp)) 1529 kpoints%xkp(3, i) = REAL(2*iz - nkp_grid(3) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(3), KIND=dp)) 1530 1531 END DO 1532 END DO 1533 END DO 1534 1535 CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, dft_control) 1536 1537 CALL set_qs_env(qs_env, kpoints=kpoints) 1538 1539 CALL timestop(handle) 1540 1541 END SUBROUTINE compute_kpoints 1542 1543! ************************************************************************************************** 1544!> \brief ... 1545!> \param para_env ... 1546!> \param dimen_RI ... 1547!> \param fm_matrix_L ... 1548!> \param my_group_L_start ... 1549!> \param my_group_L_end ... 1550!> \param my_group_L_size ... 1551!> \param my_Lrows ... 1552! ************************************************************************************************** 1553 SUBROUTINE grep_Lcols(para_env, dimen_RI, fm_matrix_L, & 1554 my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows) 1555 TYPE(cp_para_env_type), POINTER :: para_env 1556 INTEGER, INTENT(IN) :: dimen_RI 1557 TYPE(cp_fm_type), POINTER :: fm_matrix_L 1558 INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, & 1559 my_group_L_size 1560 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1561 INTENT(OUT) :: my_Lrows 1562 1563 CHARACTER(LEN=*), PARAMETER :: routineN = 'grep_Lcols', routineP = moduleN//':'//routineN 1564 1565 INTEGER :: handle, handle2, i_global, iiB, j_global, jjB, max_row_col_local, ncol_local, & 1566 ncol_rec, nrow_local, nrow_rec, proc_receive, proc_receive_static, proc_send, & 1567 proc_send_static, proc_shift 1568 INTEGER, ALLOCATABLE, DIMENSION(:) :: proc_map 1569 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, rec_col_row_info 1570 INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_rec, & 1571 row_indices, row_indices_rec 1572 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_L, local_L_internal, rec_L 1573 1574 CALL timeset(routineN, handle) 1575 1576 ALLOCATE (my_Lrows(dimen_RI, my_group_L_size)) 1577 my_Lrows = 0.0_dp 1578 1579 ! proc_map, vector that replicates the processor numbers also 1580 ! for negative and positive number > num_pe 1581 ! needed to know which is the processor, to respect to another one, 1582 ! for a given shift 1583 ALLOCATE (proc_map(-para_env%num_pe:2*para_env%num_pe - 1)) 1584 DO iiB = 0, para_env%num_pe - 1 1585 proc_map(iiB) = iiB 1586 proc_map(-iiB - 1) = para_env%num_pe - iiB - 1 1587 proc_map(para_env%num_pe + iiB) = iiB 1588 END DO 1589 1590 CALL cp_fm_get_info(matrix=fm_matrix_L, & 1591 nrow_local=nrow_local, & 1592 ncol_local=ncol_local, & 1593 row_indices=row_indices, & 1594 col_indices=col_indices, & 1595 local_data=local_L_internal) 1596 1597 ALLOCATE (local_L(nrow_local, ncol_local)) 1598 local_L = local_L_internal(1:nrow_local, 1:ncol_local) 1599 1600 max_row_col_local = MAX(nrow_local, ncol_local) 1601 CALL mp_max(max_row_col_local, para_env%group) 1602 1603 ALLOCATE (local_col_row_info(0:max_row_col_local, 2)) 1604 local_col_row_info = 0 1605 ! 0,1 nrows 1606 local_col_row_info(0, 1) = nrow_local 1607 local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local) 1608 ! 0,2 ncols 1609 local_col_row_info(0, 2) = ncol_local 1610 local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local) 1611 1612 ALLOCATE (rec_col_row_info(0:max_row_col_local, 2)) 1613 1614 ! accumulate data on my_Lrows starting from myself 1615 DO jjB = 1, ncol_local 1616 j_global = col_indices(jjB) 1617 IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN 1618 DO iiB = 1, nrow_local 1619 i_global = row_indices(iiB) 1620 my_Lrows(i_global, j_global - my_group_L_start + 1) = local_L(iiB, jjB) 1621 END DO 1622 END IF 1623 END DO 1624 1625 proc_send_static = proc_map(para_env%mepos + 1) 1626 proc_receive_static = proc_map(para_env%mepos - 1) 1627 1628 CALL timeset(routineN//"_comm", handle2) 1629 1630 DO proc_shift = 1, para_env%num_pe - 1 1631 proc_send = proc_map(para_env%mepos + proc_shift) 1632 proc_receive = proc_map(para_env%mepos - proc_shift) 1633 1634 ! first exchange information on the local data 1635 rec_col_row_info = 0 1636 CALL mp_sendrecv(local_col_row_info, proc_send_static, rec_col_row_info, proc_receive_static, para_env%group) 1637 nrow_rec = rec_col_row_info(0, 1) 1638 ncol_rec = rec_col_row_info(0, 2) 1639 1640 ALLOCATE (row_indices_rec(nrow_rec)) 1641 row_indices_rec = rec_col_row_info(1:nrow_rec, 1) 1642 1643 ALLOCATE (col_indices_rec(ncol_rec)) 1644 col_indices_rec = rec_col_row_info(1:ncol_rec, 2) 1645 1646 ALLOCATE (rec_L(nrow_rec, ncol_rec)) 1647 rec_L = 0.0_dp 1648 1649 ! then send and receive the real data 1650 CALL mp_sendrecv(local_L, proc_send_static, rec_L, proc_receive_static, para_env%group) 1651 1652 ! accumulate the received data on my_Lrows 1653 DO jjB = 1, ncol_rec 1654 j_global = col_indices_rec(jjB) 1655 IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN 1656 DO iiB = 1, nrow_rec 1657 i_global = row_indices_rec(iiB) 1658 my_Lrows(i_global, j_global - my_group_L_start + 1) = rec_L(iiB, jjB) 1659 END DO 1660 END IF 1661 END DO 1662 1663 local_col_row_info(:, :) = rec_col_row_info 1664 DEALLOCATE (local_L) 1665 ALLOCATE (local_L(nrow_rec, ncol_rec)) 1666 local_L = rec_L 1667 1668 DEALLOCATE (col_indices_rec) 1669 DEALLOCATE (row_indices_rec) 1670 DEALLOCATE (rec_L) 1671 END DO 1672 CALL timestop(handle2) 1673 1674 DEALLOCATE (local_col_row_info) 1675 DEALLOCATE (rec_col_row_info) 1676 DEALLOCATE (proc_map) 1677 DEALLOCATE (local_L) 1678 1679 CALL timestop(handle) 1680 1681 END SUBROUTINE grep_Lcols 1682 1683! ************************************************************************************************** 1684!> \brief ... 1685!> \param para_env_sub ... 1686!> \param fm_BIb_jb ... 1687!> \param BIb_jb ... 1688!> \param max_row_col_local ... 1689!> \param proc_map ... 1690!> \param local_col_row_info ... 1691!> \param my_B_virtual_end ... 1692!> \param my_B_virtual_start ... 1693! ************************************************************************************************** 1694 SUBROUTINE grep_my_integrals(para_env_sub, fm_BIb_jb, BIb_jb, max_row_col_local, & 1695 proc_map, local_col_row_info, & 1696 my_B_virtual_end, my_B_virtual_start) 1697 TYPE(cp_para_env_type), POINTER :: para_env_sub 1698 TYPE(cp_fm_type), POINTER :: fm_BIb_jb 1699 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: BIb_jb 1700 INTEGER, INTENT(IN) :: max_row_col_local 1701 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: proc_map 1702 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: local_col_row_info 1703 INTEGER, INTENT(IN) :: my_B_virtual_end, my_B_virtual_start 1704 1705 CHARACTER(LEN=*), PARAMETER :: routineN = 'grep_my_integrals', & 1706 routineP = moduleN//':'//routineN 1707 1708 INTEGER :: i_global, iiB, j_global, jjB, ncol_rec, & 1709 nrow_rec, proc_receive, proc_send, & 1710 proc_shift 1711 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: rec_col_row_info 1712 INTEGER, DIMENSION(:), POINTER :: col_indices_rec, row_indices_rec 1713 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_BI, rec_BI 1714 1715 ALLOCATE (rec_col_row_info(0:max_row_col_local, 2)) 1716 1717 rec_col_row_info(:, :) = local_col_row_info 1718 1719 nrow_rec = rec_col_row_info(0, 1) 1720 ncol_rec = rec_col_row_info(0, 2) 1721 1722 ALLOCATE (row_indices_rec(nrow_rec)) 1723 row_indices_rec = rec_col_row_info(1:nrow_rec, 1) 1724 1725 ALLOCATE (col_indices_rec(ncol_rec)) 1726 col_indices_rec = rec_col_row_info(1:ncol_rec, 2) 1727 1728 ! accumulate data on BIb_jb buffer starting from myself 1729 DO jjB = 1, ncol_rec 1730 j_global = col_indices_rec(jjB) 1731 IF (j_global >= my_B_virtual_start .AND. j_global <= my_B_virtual_end) THEN 1732 DO iiB = 1, nrow_rec 1733 i_global = row_indices_rec(iiB) 1734 BIb_jb(j_global - my_B_virtual_start + 1, i_global) = fm_BIb_jb%local_data(iiB, jjB) 1735 END DO 1736 END IF 1737 END DO 1738 1739 DEALLOCATE (row_indices_rec) 1740 DEALLOCATE (col_indices_rec) 1741 1742 IF (para_env_sub%num_pe > 1) THEN 1743 ALLOCATE (local_BI(nrow_rec, ncol_rec)) 1744 local_BI(1:nrow_rec, 1:ncol_rec) = fm_BIb_jb%local_data(1:nrow_rec, 1:ncol_rec) 1745 1746 DO proc_shift = 1, para_env_sub%num_pe - 1 1747 proc_send = proc_map(para_env_sub%mepos + proc_shift) 1748 proc_receive = proc_map(para_env_sub%mepos - proc_shift) 1749 1750 ! first exchange information on the local data 1751 rec_col_row_info = 0 1752 CALL mp_sendrecv(local_col_row_info, proc_send, rec_col_row_info, proc_receive, para_env_sub%group) 1753 nrow_rec = rec_col_row_info(0, 1) 1754 ncol_rec = rec_col_row_info(0, 2) 1755 1756 ALLOCATE (row_indices_rec(nrow_rec)) 1757 row_indices_rec = rec_col_row_info(1:nrow_rec, 1) 1758 1759 ALLOCATE (col_indices_rec(ncol_rec)) 1760 col_indices_rec = rec_col_row_info(1:ncol_rec, 2) 1761 1762 ALLOCATE (rec_BI(nrow_rec, ncol_rec)) 1763 rec_BI = 0.0_dp 1764 1765 ! then send and receive the real data 1766 CALL mp_sendrecv(local_BI, proc_send, rec_BI, proc_receive, para_env_sub%group) 1767 1768 ! accumulate the received data on BIb_jb buffer 1769 DO jjB = 1, ncol_rec 1770 j_global = col_indices_rec(jjB) 1771 IF (j_global >= my_B_virtual_start .AND. j_global <= my_B_virtual_end) THEN 1772 DO iiB = 1, nrow_rec 1773 i_global = row_indices_rec(iiB) 1774 BIb_jb(j_global - my_B_virtual_start + 1, i_global) = rec_BI(iiB, jjB) 1775 END DO 1776 END IF 1777 END DO 1778 1779 DEALLOCATE (col_indices_rec) 1780 DEALLOCATE (row_indices_rec) 1781 DEALLOCATE (rec_BI) 1782 END DO 1783 1784 DEALLOCATE (local_BI) 1785 END IF 1786 1787 DEALLOCATE (rec_col_row_info) 1788 1789 END SUBROUTINE grep_my_integrals 1790 1791! ************************************************************************************************** 1792!> \brief compute three center overlap integrals and write them to mat_3_overl_int 1793!> \param mat_3c_overl_int ... 1794!> \param mat_3c_overl_int_mao_for_occ ... 1795!> \param mat_3c_overl_int_mao_for_virt ... 1796!> \param mat_munu ... 1797!> \param mat_munu_mao_occ_virt ... 1798!> \param qs_env ... 1799!> \param my_group_L_starts_im_time ... 1800!> \param my_group_L_ends_im_time ... 1801!> \param my_group_L_sizes_im_time ... 1802!> \param dimen_RI ... 1803!> \param cut_RI ... 1804!> \param sab_orb_sub ... 1805!> \param sab_orb_all ... 1806!> \param do_mao ... 1807!> \param mao_coeff_repl_occ ... 1808!> \param mao_coeff_repl_virt ... 1809!> \param num_diff_blk ... 1810! ************************************************************************************************** 1811 SUBROUTINE compute_3c_overlap_int(mat_3c_overl_int, mat_3c_overl_int_mao_for_occ, & 1812 mat_3c_overl_int_mao_for_virt, mat_munu, mat_munu_mao_occ_virt, & 1813 qs_env, my_group_L_starts_im_time, & 1814 my_group_L_ends_im_time, my_group_L_sizes_im_time, dimen_RI, & 1815 cut_RI, sab_orb_sub, sab_orb_all, do_mao, mao_coeff_repl_occ, & 1816 mao_coeff_repl_virt, num_diff_blk) 1817 1818 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int, & 1819 mat_3c_overl_int_mao_for_occ, & 1820 mat_3c_overl_int_mao_for_virt 1821 TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu, mat_munu_mao_occ_virt 1822 TYPE(qs_environment_type), POINTER :: qs_env 1823 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: my_group_L_starts_im_time, & 1824 my_group_L_ends_im_time, & 1825 my_group_L_sizes_im_time 1826 INTEGER, INTENT(IN) :: dimen_RI, cut_RI 1827 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1828 POINTER :: sab_orb_sub, sab_orb_all 1829 LOGICAL, INTENT(IN) :: do_mao 1830 TYPE(two_dim_real_array), ALLOCATABLE, & 1831 DIMENSION(:), INTENT(IN) :: mao_coeff_repl_occ, mao_coeff_repl_virt 1832 INTEGER, INTENT(IN) :: num_diff_blk 1833 1834 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_overlap_int', & 1835 routineP = moduleN//':'//routineN 1836 1837 INTEGER :: acell, atom_RI, atom_RI_end, atom_RI_start, bcell, block_end_col, block_end_row, & 1838 block_start_col, block_start_row, col, col_basis_size_mao_occ, col_basis_size_mao_virt, & 1839 end_col_from_LLL_mao_occ, end_col_from_LLL_mao_virt, end_row_from_LLL_mao_occ, & 1840 end_row_from_LLL_mao_virt, handle, handle2, i_cut_RI, i_img, iatom, & 1841 iatom_basis_size_mao_occ, iatom_basis_size_mao_virt, iatom_outer, iblk, iblk_send, ikind, & 1842 img_col, img_row, iset, isgf_RI, isgfa, j_img, jatom, jatom_basis_size_mao_occ, & 1843 jatom_basis_size_mao_virt, jatom_outer, jkind, jset, jsgfb, kind_RI, LLL, LLL_set_end 1844 INTEGER :: LLL_set_start, maxval_1, maxval_2, maxval_3, mepos, minval_1, minval_2, minval_3, & 1845 my_group_L_end, my_group_L_size, my_group_L_start, natom, nco_RI, ncoa, ncob, nimg, & 1846 nkind, nset_RI, nseta, nsetb, nthread, offset_col_from_LLL, offset_row_from_LLL, row, & 1847 row_basis_size_mao_occ, row_basis_size_mao_virt, set_RI, sgf_RI, sgfa, sgfb, & 1848 start_col_from_LLL_mao_occ, start_col_from_LLL_mao_virt, start_row_from_LLL_mao_occ, & 1849 start_row_from_LLL_mao_virt 1850 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_RI_index, kind_of 1851 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: indices_for_uncommon_blocks 1852 INTEGER, DIMENSION(3) :: bcellvec, cell_vec, outer_cell_vec, & 1853 size_cell_to_index 1854 INTEGER, DIMENSION(:), POINTER :: blk_sizes_mao_occ, blk_sizes_mao_virt, col_blk_sizes, & 1855 l_max_RI, l_min_RI, la_max, la_min, lb_max, lb_min, npgf_RI, npgfa, npgfb, nsgf_RI, & 1856 nsgfa, nsgfb, row_blk_end, row_blk_sizes, row_blk_start 1857 INTEGER, DIMENSION(:, :), POINTER :: first_sgf_RI, first_sgfa, first_sgfb 1858 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 1859 LOGICAL :: do_kpoints_cubic_RPA, found 1860 REAL(KIND=dp) :: dab, daRI, dbRI 1861 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: block, block_mao_for_occ, & 1862 block_mao_for_occ_transposed, block_mao_for_virt, block_mao_for_virt_transposed, & 1863 block_transposed, temp_mat_mao_occ_virt, temp_mat_mao_virt_occ 1864 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: s_abRI, s_abRI_contr, & 1865 s_abRI_contr_transposed 1866 REAL(KIND=dp), DIMENSION(3) :: rab, rab_outer, raRI, rbRI 1867 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_RI 1868 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dummy_block, rpgf_RI, rpgfa, rpgfb, & 1869 sphi_a, sphi_b, sphi_RI, zet_RI, zeta, & 1870 zetb 1871 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1872 TYPE(cell_type), POINTER :: cell 1873 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_tmp 1874 TYPE(dft_control_type), POINTER :: dft_control 1875 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_RI_tmp 1876 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, basis_set_RI 1877 TYPE(kpoint_type), POINTER :: kpoints 1878 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1879 TYPE(neighbor_list_iterator_p_type), & 1880 DIMENSION(:), POINTER :: nl_iterator, nl_iterator_outer 1881 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1882 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1883 TYPE(two_dim_real_array), ALLOCATABLE, & 1884 DIMENSION(:) :: blocks_to_send_around 1885 1886 CALL timeset(routineN, handle) 1887 1888 NULLIFY (qs_kind_set, atomic_kind_set, cell, particle_set, molecule_kind_set, basis_set_RI, basis_set_a, basis_set_b, & 1889 nl_iterator) 1890 1891 ! get stuff 1892 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, & 1893 cell=cell, particle_set=particle_set, molecule_kind_set=molecule_kind_set, & 1894 nkind=nkind, natom=natom, kpoints=kpoints, dft_control=dft_control) 1895 1896 do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints 1897 1898 IF (do_kpoints_cubic_RPA) THEN 1899 ! No MAOs kpoints 1900 CPASSERT(.NOT. do_mao) 1901 nimg = dft_control%nimages 1902 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index) 1903 size_cell_to_index(1) = SIZE(cell_to_index, 1) 1904 size_cell_to_index(2) = SIZE(cell_to_index, 2) 1905 size_cell_to_index(3) = SIZE(cell_to_index, 3) 1906 ! 1: row, 2: col, 3: i_img, 4: j_img, 5: cut_RI, 6: proc_to_send 1907 ALLOCATE (indices_for_uncommon_blocks(num_diff_blk, 6)) 1908 indices_for_uncommon_blocks(:, :) = 0 1909 ALLOCATE (blocks_to_send_around(num_diff_blk)) 1910 iblk_send = 0 1911 minval_1 = MINVAL(kpoints%index_to_cell(1, :)) 1912 maxval_1 = MAXVAL(kpoints%index_to_cell(1, :)) 1913 minval_2 = MINVAL(kpoints%index_to_cell(2, :)) 1914 maxval_2 = MAXVAL(kpoints%index_to_cell(2, :)) 1915 minval_3 = MINVAL(kpoints%index_to_cell(3, :)) 1916 maxval_3 = MAXVAL(kpoints%index_to_cell(3, :)) 1917 ELSE 1918 nimg = 1 1919 END IF 1920 1921 ALLOCATE (row_blk_start(natom)) 1922 ALLOCATE (row_blk_end(natom)) 1923 1924 ALLOCATE (basis_set_RI_tmp(nkind)) 1925 CALL basis_set_list_setup(basis_set_RI_tmp, "RI_AUX", qs_kind_set) 1926 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, & 1927 basis=basis_set_RI_tmp) 1928 DEALLOCATE (basis_set_RI_tmp) 1929 1930 ALLOCATE (atom_from_RI_index(dimen_RI)) 1931 1932 DO LLL = 1, dimen_RI 1933 1934 DO iatom = 1, natom 1935 1936 IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN 1937 1938 atom_from_RI_index(LLL) = iatom 1939 1940 END IF 1941 1942 END DO 1943 1944 END DO 1945 1946 CALL dbcsr_get_info(mat_munu%matrix, & 1947 row_blk_size=row_blk_sizes, & 1948 col_blk_size=col_blk_sizes) 1949 1950 IF (do_mao) THEN 1951 1952 CALL dbcsr_get_info(mat_munu_mao_occ_virt%matrix, & 1953 row_blk_size=blk_sizes_mao_occ, & 1954 col_blk_size=blk_sizes_mao_virt) 1955 1956 END IF 1957 1958 ALLOCATE (kind_of(natom)) 1959 1960 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of) 1961 1962 DO i_cut_RI = 1, cut_RI 1963 1964 my_group_L_start = my_group_L_starts_im_time(i_cut_RI) 1965 my_group_L_end = my_group_L_ends_im_time(i_cut_RI) 1966 my_group_L_size = my_group_L_sizes_im_time(i_cut_RI) 1967 1968 atom_RI_start = atom_from_RI_index(my_group_L_start) 1969 atom_RI_end = atom_from_RI_index(my_group_L_end) 1970 1971 DO atom_RI = atom_RI_start, atom_RI_end 1972 1973 kind_RI = kind_of(atom_RI) 1974 CALL get_qs_kind(qs_kind=qs_kind_set(kind_RI), basis_set=basis_set_RI, basis_type="RI_AUX") 1975 1976 first_sgf_RI => basis_set_RI%first_sgf 1977 l_max_RI => basis_set_RI%lmax 1978 l_min_RI => basis_set_RI%lmin 1979 npgf_RI => basis_set_RI%npgf 1980 nset_RI = basis_set_RI%nset 1981 nsgf_RI => basis_set_RI%nsgf_set 1982 rpgf_RI => basis_set_RI%pgf_radius 1983 set_radius_RI => basis_set_RI%set_radius 1984 sphi_RI => basis_set_RI%sphi 1985 zet_RI => basis_set_RI%zet 1986 1987 CALL neighbor_list_iterator_create(nl_iterator_outer, sab_orb_all) 1988 DO WHILE (neighbor_list_iterate(nl_iterator_outer) == 0) 1989 1990 CALL get_iterator_info(nl_iterator_outer, & 1991 iatom=iatom_outer, jatom=jatom_outer, r=rab_outer, & 1992 cell=outer_cell_vec) 1993 1994 IF (iatom_outer .NE. atom_RI .AND. jatom_outer .NE. atom_RI) CYCLE 1995 1996 nthread = 1 1997!$ nthread = omp_get_max_threads() 1998 1999 CALL neighbor_list_iterator_create(nl_iterator, sab_orb_sub, nthread=nthread) 2000 2001!$OMP PARALLEL DEFAULT(NONE) & 2002!$OMP SHARED (nthread,atom_RI,first_sgf_RI,l_max_RI,l_min_RI,npgf_RI,nset_RI,nsgf_RI,rpgf_RI,set_radius_RI,sphi_RI,zet_RI,& 2003!$OMP iatom_outer,jatom_outer,rab_outer,kind_of,my_group_L_start,row_blk_start,& 2004!$OMP mat_3c_overl_int,nl_iterator,qs_kind_set,ncoset,my_group_L_size,row_blk_sizes,& 2005!$OMP col_blk_sizes, i_cut_RI, do_mao, blk_sizes_mao_occ, blk_sizes_mao_virt, outer_cell_vec, & 2006!$OMP mao_coeff_repl_occ, mao_coeff_repl_virt, do_kpoints_cubic_RPA, & 2007!$OMP mat_3c_overl_int_mao_for_occ, mat_3c_overl_int_mao_for_virt, kpoints, cell_to_index, nimg, & 2008!$OMP minval_1, minval_2, minval_3, maxval_1, maxval_2, maxval_3) & 2009!$OMP PRIVATE (mepos,iatom,jatom,rab,set_RI,raRI,rbRI,dab,daRI,dbRI,ikind,first_sgfa,la_max,la_min,npgfa,nseta,nsgfa,rpgfa,& 2010!$OMP set_radius_a,sphi_a,zeta,jkind,first_sgfb,lb_max,lb_min,npgfb,nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,& 2011!$OMP iset,jset,LLL_set_start,nco_RI,ncoa,ncob,sgf_RI,sgfa,sgfb,isgf_RI,isgfa,jsgfb,LLL,row,col,& 2012!$OMP block,block_transposed,LLL_set_end, start_row_from_LLL_mao_occ, start_col_from_LLL_mao_occ, & 2013!$OMP start_row_from_LLL_mao_virt, start_col_from_LLL_mao_virt, end_row_from_LLL_mao_occ, & 2014!$OMP end_col_from_LLL_mao_occ, end_row_from_LLL_mao_virt, end_col_from_LLL_mao_virt, & 2015!$OMP row_basis_size_mao_occ, col_basis_size_mao_occ, row_basis_size_mao_virt, col_basis_size_mao_virt, & 2016!$OMP offset_row_from_LLL,offset_col_from_LLL,block_start_row,block_end_row,block_start_col,block_end_col,& 2017!$OMP handle2,s_abRI,s_abRI_contr,s_abRI_contr_transposed,basis_set_a,basis_set_b, jatom_basis_size_mao_occ, & 2018!$OMP iatom_basis_size_mao_virt, jatom_basis_size_mao_virt, temp_mat_mao_occ_virt, block_mao_for_occ, & 2019!$OMP block_mao_for_occ_transposed, block_mao_for_virt, block_mao_for_virt_transposed, temp_mat_mao_virt_occ, & 2020!$OMP iatom_basis_size_mao_occ, found, dummy_block, cell_vec, acell, bcell, bcellvec, & 2021!$OMP img_row, img_col, i_img, j_img) 2022 2023 mepos = 0 2024!$ mepos = omp_get_thread_num() 2025 2026 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 2027 2028 CALL get_iterator_info(nl_iterator, mepos=mepos, & 2029 iatom=iatom, jatom=jatom, r=rab, & 2030 cell=cell_vec) 2031 2032 DO i_img = 1, nimg 2033 2034 DO j_img = 1, nimg 2035 2036 IF (atom_RI .EQ. iatom_outer .AND. & 2037 iatom .NE. jatom_outer .AND. & 2038 jatom .NE. jatom_outer) & 2039 CYCLE 2040 2041 IF (atom_RI .EQ. jatom_outer .AND. & 2042 iatom .NE. iatom_outer .AND. & 2043 jatom .NE. iatom_outer) & 2044 CYCLE 2045 2046 IF (iatom .EQ. iatom_outer .AND. & 2047 atom_RI .EQ. jatom_outer) THEN 2048 2049 IF (do_kpoints_cubic_RPA) THEN 2050 2051 IF (outer_cell_vec(1) < minval_1) CYCLE 2052 IF (outer_cell_vec(1) > maxval_1) CYCLE 2053 IF (outer_cell_vec(2) < minval_2) CYCLE 2054 IF (outer_cell_vec(2) > maxval_2) CYCLE 2055 IF (outer_cell_vec(3) < minval_3) CYCLE 2056 IF (outer_cell_vec(3) > maxval_3) CYCLE 2057 2058 acell = cell_to_index(-outer_cell_vec(1), -outer_cell_vec(2), -outer_cell_vec(3)) 2059 2060 IF (.NOT. (acell == i_img)) CYCLE 2061 2062 bcellvec = -outer_cell_vec + cell_vec 2063 2064 IF (bcellvec(1) < minval_1) CYCLE 2065 IF (bcellvec(1) > maxval_1) CYCLE 2066 IF (bcellvec(2) < minval_2) CYCLE 2067 IF (bcellvec(2) > maxval_2) CYCLE 2068 IF (bcellvec(3) < minval_3) CYCLE 2069 IF (bcellvec(3) > maxval_3) CYCLE 2070 2071 bcell = cell_to_index(bcellvec(1), bcellvec(2), bcellvec(3)) 2072 2073 IF (.NOT. (bcell == j_img)) CYCLE 2074 2075 END IF 2076 2077 raRI = rab_outer 2078 rbRI = raRI - rab 2079 2080 ELSE IF (atom_RI .EQ. iatom_outer .AND. & 2081 iatom .EQ. jatom_outer) THEN 2082 2083 IF (do_kpoints_cubic_RPA) THEN 2084 ! for kpoints we have the non-symmetric neighbor list 2085 CYCLE 2086 END IF 2087 2088 raRI = -rab_outer 2089 rbRI = raRI - rab 2090 2091 ELSE IF (jatom .EQ. iatom_outer .AND. & 2092 atom_RI .EQ. jatom_outer) THEN 2093 2094 IF (do_kpoints_cubic_RPA) THEN 2095 ! for kpoints we have the non-symmetric neighbor list 2096 CYCLE 2097 END IF 2098 2099 rbRI = rab_outer 2100 raRI = rbRI + rab 2101 2102 ELSE IF (atom_RI .EQ. iatom_outer .AND. & 2103 jatom .EQ. jatom_outer) THEN 2104 2105 IF (do_kpoints_cubic_RPA) THEN 2106 ! for kpoints we have the non-symmetric neighbor list 2107 CYCLE 2108 END IF 2109 2110 rbRI = -rab_outer 2111 raRI = rbRI + rab 2112 2113 END IF 2114 2115 DO set_RI = 1, nset_RI 2116 2117 dab = SQRT(SUM(rab**2)) 2118 daRI = SQRT(SUM(raRI**2)) 2119 dbRI = SQRT(SUM(rbRI**2)) 2120 2121 ikind = kind_of(iatom) 2122 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a) 2123 first_sgfa => basis_set_a%first_sgf 2124 la_max => basis_set_a%lmax 2125 la_min => basis_set_a%lmin 2126 npgfa => basis_set_a%npgf 2127 nseta = basis_set_a%nset 2128 nsgfa => basis_set_a%nsgf_set 2129 rpgfa => basis_set_a%pgf_radius 2130 set_radius_a => basis_set_a%set_radius 2131 sphi_a => basis_set_a%sphi 2132 zeta => basis_set_a%zet 2133 2134 jkind = kind_of(jatom) 2135 CALL get_qs_kind(qs_kind=qs_kind_set(jkind), basis_set=basis_set_b) 2136 first_sgfb => basis_set_b%first_sgf 2137 lb_max => basis_set_b%lmax 2138 lb_min => basis_set_b%lmin 2139 npgfb => basis_set_b%npgf 2140 nsetb = basis_set_b%nset 2141 nsgfb => basis_set_b%nsgf_set 2142 rpgfb => basis_set_b%pgf_radius 2143 set_radius_b => basis_set_b%set_radius 2144 sphi_b => basis_set_b%sphi 2145 zetb => basis_set_b%zet 2146 2147 DO iset = 1, nseta 2148 2149 IF (set_radius_a(iset) + set_radius_RI(set_RI) < daRI) CYCLE 2150 2151 DO jset = 1, nsetb 2152 2153 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 2154 IF (set_radius_b(jset) + set_radius_RI(set_RI) < dbRI) CYCLE 2155 2156 nco_RI = npgf_RI(set_RI)*ncoset(l_max_RI(set_RI)) 2157 ncoa = npgfa(iset)*ncoset(la_max(iset)) 2158 ncob = npgfb(jset)*ncoset(lb_max(jset)) 2159 2160 sgf_RI = first_sgf_RI(1, set_RI) 2161 sgfa = first_sgfa(1, iset) 2162 sgfb = first_sgfb(1, jset) 2163 2164 LLL_set_start = 1 + sgf_RI - 1 - my_group_L_start + 1 + row_blk_start(atom_RI) - 1 2165 2166 IF (LLL_set_start > my_group_L_size) CYCLE 2167 2168 LLL_set_end = nsgf_RI(set_RI) + sgf_RI - 1 - my_group_L_start + 1 + row_blk_start(atom_RI) - 1 2169 2170 IF (LLL_set_end < 1) CYCLE 2171 2172 IF (ncoa*ncob*nco_RI > 0) THEN 2173 2174 ALLOCATE (s_abRI(ncoa, ncob, nco_RI)) 2175 s_abRI(:, :, :) = 0.0_dp 2176 2177 CALL overlap3(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 2178 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), & 2179 l_max_RI(set_RI), npgf_RI(set_RI), zet_RI(:, set_RI), rpgf_RI(:, set_RI), & 2180 l_min_RI(set_RI), rab, dab, raRI, daRI, rbRI, dbRI, s_abRI) 2181 2182 ALLOCATE (s_abRI_contr(nsgfa(iset), nsgfb(jset), nsgf_RI(set_RI))) 2183 2184 CALL abc_contract(s_abRI_contr, s_abRI, & 2185 sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_RI(:, sgf_RI:), & 2186 ncoa, ncob, nco_RI, nsgfa(iset), nsgfb(jset), nsgf_RI(set_RI)) 2187 2188 ALLOCATE (s_abRI_contr_transposed(nsgfb(jset), nsgfa(iset), nsgf_RI(set_RI))) 2189 2190 DO isgf_RI = 1, nsgf_RI(set_RI) 2191 DO isgfa = 1, nsgfa(iset) 2192 DO jsgfb = 1, nsgfb(jset) 2193 s_abRI_contr_transposed(jsgfb, isgfa, isgf_RI) = & 2194 s_abRI_contr(isgfa, jsgfb, isgf_RI) 2195 END DO 2196 END DO 2197 END DO 2198 2199 DO isgf_RI = 1, nsgf_RI(set_RI) 2200 2201 LLL = isgf_RI + sgf_RI - 1 - my_group_L_start + 1 + row_blk_start(atom_RI) - 1 2202 2203 IF (LLL < 1) CYCLE 2204 IF (LLL > my_group_L_size) CYCLE 2205 2206 IF (iatom < jatom) THEN 2207 2208 row = iatom 2209 col = jatom 2210 img_row = i_img 2211 img_col = j_img 2212 2213 ELSE IF (iatom >= jatom) THEN 2214 2215 row = jatom 2216 col = iatom 2217 img_row = j_img 2218 img_col = i_img 2219 2220 END IF 2221 2222 ALLOCATE (block(row_blk_sizes(row), col_blk_sizes(col)*my_group_L_size)) 2223 block = 0.0_dp 2224 2225 ALLOCATE (block_transposed(col_blk_sizes(col), row_blk_sizes(row)*my_group_L_size)) 2226 block_transposed = 0.0_dp 2227 2228 offset_row_from_LLL = (LLL - 1)*row_blk_sizes(row) 2229 offset_col_from_LLL = (LLL - 1)*col_blk_sizes(col) 2230 2231 IF (iatom < jatom) THEN 2232 2233 block_start_row = sgfa 2234 block_end_row = sgfa + nsgfa(iset) - 1 2235 block_start_col = sgfb 2236 block_end_col = sgfb + nsgfb(jset) - 1 2237 2238 ! factor 0.5 is necessary due to double filling due to double iterate loop 2239 block(block_start_row:block_end_row, & 2240 block_start_col + offset_col_from_LLL:block_end_col + offset_col_from_LLL) = & 2241 0.5_dp*s_abRI_contr(:, :, isgf_RI) 2242 2243 block_transposed(block_start_col:block_end_col, & 2244 block_start_row + offset_row_from_LLL: & 2245 block_end_row + offset_row_from_LLL) = & 2246 0.5_dp*s_abRI_contr_transposed(:, :, isgf_RI) 2247 2248 ELSE IF (iatom > jatom) THEN 2249 2250 block_start_row = sgfb 2251 block_end_row = sgfb + nsgfb(jset) - 1 2252 block_start_col = sgfa 2253 block_end_col = sgfa + nsgfa(iset) - 1 2254 2255 ! factor 0.5 is necessary due to double filling due to double iterate loop 2256 block(block_start_row:block_end_row, & 2257 block_start_col + offset_col_from_LLL:block_end_col + offset_col_from_LLL) = & 2258 0.5_dp*s_abRI_contr_transposed(:, :, isgf_RI) 2259 2260 block_transposed(block_start_col:block_end_col, & 2261 block_start_row + offset_row_from_LLL: & 2262 block_end_row + offset_row_from_LLL) = & 2263 0.5_dp*s_abRI_contr(:, :, isgf_RI) 2264 2265 ELSE IF (iatom .EQ. jatom) THEN 2266 2267 block_start_row = sgfa 2268 block_end_row = sgfa + nsgfa(iset) - 1 2269 block_start_col = sgfb + offset_col_from_LLL 2270 block_end_col = sgfb + nsgfb(jset) - 1 + offset_col_from_LLL 2271 2272 block(block_start_row:block_end_row, & 2273 block_start_col:block_end_col) = s_abRI_contr(:, :, isgf_RI) 2274 2275 END IF 2276 2277 CALL timeset(routineN//"_put_dbcsr", handle2) 2278 2279!$OMP CRITICAL(putblock) 2280 NULLIFY (dummy_block) 2281 CALL dbcsr_get_block_p(matrix=mat_3c_overl_int(i_cut_RI, img_row, img_col)%matrix, & 2282 row=row, col=col, block=dummy_block, found=found) 2283 IF (found) THEN 2284 CALL dbcsr_put_block(mat_3c_overl_int(i_cut_RI, img_row, img_col)%matrix, & 2285 row=row, col=col, block=block, summation=.TRUE.) 2286 END IF 2287!$OMP END CRITICAL(putblock) 2288 2289 IF (row .NE. col) THEN 2290!$OMP CRITICAL(putblock) 2291 NULLIFY (dummy_block) 2292 CALL dbcsr_get_block_p(matrix=mat_3c_overl_int(i_cut_RI, img_col, img_row)%matrix, & 2293 row=col, col=row, block=dummy_block, found=found) 2294 IF (found) THEN 2295 CALL dbcsr_put_block(mat_3c_overl_int(i_cut_RI, img_col, img_row)%matrix, & 2296 row=col, col=row, block=block_transposed, summation=.TRUE.) 2297 END IF 2298!$OMP END CRITICAL(putblock) 2299 END IF 2300 2301 DEALLOCATE (block, block_transposed) 2302 2303 IF (do_mao) THEN 2304 2305 iatom_basis_size_mao_occ = blk_sizes_mao_occ(iatom) 2306 jatom_basis_size_mao_occ = blk_sizes_mao_occ(jatom) 2307 iatom_basis_size_mao_virt = blk_sizes_mao_virt(iatom) 2308 jatom_basis_size_mao_virt = blk_sizes_mao_virt(jatom) 2309 2310 row_basis_size_mao_occ = blk_sizes_mao_occ(row) 2311 col_basis_size_mao_occ = blk_sizes_mao_occ(col) 2312 row_basis_size_mao_virt = blk_sizes_mao_virt(row) 2313 col_basis_size_mao_virt = blk_sizes_mao_virt(col) 2314 2315 ALLOCATE (temp_mat_mao_occ_virt(iatom_basis_size_mao_occ, jatom_basis_size_mao_virt)) 2316 ALLOCATE (temp_mat_mao_virt_occ(iatom_basis_size_mao_virt, jatom_basis_size_mao_occ)) 2317 2318 ! transform from the primary Gaussian basis into the MAO basis 2319 CALL ab_contract(temp_mat_mao_occ_virt(:, :), s_abRI_contr(:, :, isgf_RI), & 2320 mao_coeff_repl_occ(iatom)%array(sgfa:sgfa + nsgfa(iset) - 1, :), & 2321 mao_coeff_repl_virt(jatom)%array(sgfb:sgfb + nsgfb(jset) - 1, :), & 2322 nsgfa(iset), nsgfb(jset), & 2323 iatom_basis_size_mao_occ, jatom_basis_size_mao_virt) 2324 2325 CALL ab_contract(temp_mat_mao_virt_occ(:, :), s_abRI_contr(:, :, isgf_RI), & 2326 mao_coeff_repl_virt(iatom)%array(sgfa:sgfa + nsgfa(iset) - 1, :), & 2327 mao_coeff_repl_occ(jatom)%array(sgfb:sgfb + nsgfb(jset) - 1, :), & 2328 nsgfa(iset), nsgfb(jset), & 2329 iatom_basis_size_mao_virt, jatom_basis_size_mao_occ) 2330 2331 start_row_from_LLL_mao_occ = (LLL - 1)*row_basis_size_mao_occ + 1 2332 start_col_from_LLL_mao_occ = (LLL - 1)*col_basis_size_mao_occ + 1 2333 2334 start_row_from_LLL_mao_virt = (LLL - 1)*row_basis_size_mao_virt + 1 2335 start_col_from_LLL_mao_virt = (LLL - 1)*col_basis_size_mao_virt + 1 2336 2337 end_row_from_LLL_mao_occ = LLL*row_basis_size_mao_occ 2338 end_col_from_LLL_mao_occ = LLL*col_basis_size_mao_occ 2339 2340 end_row_from_LLL_mao_virt = LLL*row_basis_size_mao_virt 2341 end_col_from_LLL_mao_virt = LLL*col_basis_size_mao_virt 2342 2343 ALLOCATE (block_mao_for_occ(row_basis_size_mao_occ, & 2344 col_basis_size_mao_virt*my_group_L_size)) 2345 block_mao_for_occ = 0.0_dp 2346 2347 ALLOCATE (block_mao_for_occ_transposed(col_basis_size_mao_occ, & 2348 row_basis_size_mao_virt*my_group_L_size)) 2349 block_mao_for_occ_transposed = 0.0_dp 2350 2351 ALLOCATE (block_mao_for_virt(row_basis_size_mao_virt, & 2352 col_basis_size_mao_occ*my_group_L_size)) 2353 block_mao_for_virt = 0.0_dp 2354 2355 ALLOCATE (block_mao_for_virt_transposed(col_basis_size_mao_virt, & 2356 row_basis_size_mao_occ*my_group_L_size)) 2357 block_mao_for_virt_transposed = 0.0_dp 2358 2359 IF (iatom < jatom) THEN 2360 2361 block_mao_for_occ(:, start_col_from_LLL_mao_virt:end_col_from_LLL_mao_virt) = & 2362 0.5_dp*temp_mat_mao_occ_virt(:, :) 2363 2364 block_mao_for_occ_transposed(:, start_row_from_LLL_mao_virt: & 2365 end_row_from_LLL_mao_virt) = & 2366 0.5_dp*TRANSPOSE(temp_mat_mao_virt_occ(:, :)) 2367 2368 block_mao_for_virt(:, start_col_from_LLL_mao_occ:end_col_from_LLL_mao_occ) = & 2369 0.5_dp*temp_mat_mao_virt_occ(:, :) 2370 2371 block_mao_for_virt_transposed(:, start_row_from_LLL_mao_occ: & 2372 end_row_from_LLL_mao_occ) = & 2373 0.5_dp*TRANSPOSE(temp_mat_mao_occ_virt(:, :)) 2374 2375 ELSE IF (iatom > jatom) THEN 2376 2377 block_mao_for_occ(:, start_col_from_LLL_mao_virt:end_col_from_LLL_mao_virt) = & 2378 0.5_dp*TRANSPOSE(temp_mat_mao_virt_occ(:, :)) 2379 2380 block_mao_for_occ_transposed(:, start_row_from_LLL_mao_virt: & 2381 end_row_from_LLL_mao_virt) = & 2382 0.5_dp*temp_mat_mao_occ_virt(:, :) 2383 2384 block_mao_for_virt(:, start_col_from_LLL_mao_occ:end_col_from_LLL_mao_occ) = & 2385 0.5_dp*TRANSPOSE(temp_mat_mao_occ_virt(:, :)) 2386 2387 block_mao_for_virt_transposed(:, start_row_from_LLL_mao_occ: & 2388 end_row_from_LLL_mao_occ) = & 2389 0.5_dp*temp_mat_mao_virt_occ(:, :) 2390 2391 ELSE IF (iatom .EQ. jatom) THEN 2392 2393 block_mao_for_occ(:, start_col_from_LLL_mao_virt:end_col_from_LLL_mao_virt) = & 2394 temp_mat_mao_occ_virt(:, :) 2395 2396 block_mao_for_virt(:, start_col_from_LLL_mao_occ:end_col_from_LLL_mao_occ) = & 2397 temp_mat_mao_virt_occ(:, :) 2398 2399 END IF 2400 2401!$OMP CRITICAL(putblock) 2402 CALL dbcsr_put_block(mat_3c_overl_int_mao_for_occ(i_cut_RI, i_img, j_img)%matrix, & 2403 row=row, col=col, block=block_mao_for_occ, summation=.TRUE.) 2404!$OMP END CRITICAL(putblock) 2405 2406!$OMP CRITICAL(putblock) 2407 CALL dbcsr_put_block(mat_3c_overl_int_mao_for_virt(i_cut_RI, i_img, j_img)%matrix, & 2408 row=row, col=col, block=block_mao_for_virt, summation=.TRUE.) 2409!$OMP END CRITICAL(putblock) 2410 2411 IF (row .NE. col) THEN 2412!$OMP CRITICAL(putblock) 2413 CALL dbcsr_put_block(mat_3c_overl_int_mao_for_occ(i_cut_RI, j_img, i_img)%matrix, & 2414 row=col, col=row, block=block_mao_for_occ_transposed, & 2415 summation=.TRUE.) 2416!$OMP END CRITICAL(putblock) 2417 2418!$OMP CRITICAL(putblock) 2419 CALL dbcsr_put_block(mat_3c_overl_int_mao_for_virt(i_cut_RI, j_img, i_img)%matrix, & 2420 row=col, col=row, block=block_mao_for_virt_transposed, & 2421 summation=.TRUE.) 2422!$OMP END CRITICAL(putblock) 2423 2424 END IF 2425 2426 DEALLOCATE (temp_mat_mao_occ_virt, temp_mat_mao_virt_occ, block_mao_for_occ, & 2427 block_mao_for_occ_transposed, block_mao_for_virt, & 2428 block_mao_for_virt_transposed) 2429 2430 END IF 2431 2432 CALL timestop(handle2) 2433 2434 END DO ! single RI basis functions 2435 2436 DEALLOCATE (s_abRI, s_abRI_contr, s_abRI_contr_transposed) 2437 2438 END IF ! number of triples > 0 2439 2440 END DO ! jset 2441 2442 END DO ! iset 2443 2444 END DO ! set RI 2445 2446 END DO ! j_img 2447 2448 END DO ! i_img 2449 2450 END DO ! inner neighbor list 2451 2452!$OMP END PARALLEL 2453 2454 CALL neighbor_list_iterator_release(nl_iterator) 2455 2456 END DO ! outer neighbor list 2457 2458 CALL neighbor_list_iterator_release(nl_iterator_outer) 2459 2460 END DO ! atom_RI 2461 2462 END DO ! RI index ranges 2463 2464 ! when doing kpoints, there might be some zero blocks 2465 IF (do_kpoints_cubic_RPA) THEN 2466 2467 NULLIFY (mat_3c_overl_int_tmp) 2468 CALL dbcsr_allocate_matrix_set(mat_3c_overl_int_tmp, cut_RI, nimg, nimg) 2469 2470 DO i_img = 1, nimg 2471 DO j_img = 1, nimg 2472 DO i_cut_RI = 1, cut_RI 2473 2474 CALL dbcsr_filter(mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix, 1.0E-30_dp) 2475 2476 ALLOCATE (mat_3c_overl_int_tmp(i_cut_RI, i_img, j_img)%matrix) 2477 CALL dbcsr_create(mat_3c_overl_int_tmp(i_cut_RI, i_img, j_img)%matrix, & 2478 template=mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix) 2479 CALL dbcsr_copy(mat_3c_overl_int_tmp(i_cut_RI, i_img, j_img)%matrix, & 2480 mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix) 2481 CALL dbcsr_set(mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix, 0.0_dp) 2482 CALL dbcsr_filter(mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix, 1.0E-30_dp) 2483 CALL dbcsr_complete_redistribute(mat_3c_overl_int_tmp(i_cut_RI, i_img, j_img)%matrix, & 2484 mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix, summation=.TRUE.) 2485 2486 END DO 2487 END DO 2488 END DO 2489 2490 DO iblk = 1, num_diff_blk 2491 DEALLOCATE (blocks_to_send_around(iblk)%array) 2492 END DO 2493 2494 DEALLOCATE (indices_for_uncommon_blocks, blocks_to_send_around) 2495 CALL dbcsr_deallocate_matrix_set(mat_3c_overl_int_tmp) 2496 2497 END IF 2498 2499 DEALLOCATE (atom_from_RI_index) 2500 2501 DEALLOCATE (row_blk_start, row_blk_end) 2502 2503 CALL timestop(handle) 2504 2505 END SUBROUTINE compute_3c_overlap_int 2506 2507! ************************************************************************************************** 2508!> \brief compute three center overlap integrals and write them to mat_3_overl_int 2509!> \param mat_3c_overl_int ... 2510!> \param mat_munu ... 2511!> \param qs_env ... 2512!> \param dimen_RI ... 2513!> \param cut_RI ... 2514!> \param unit_nr ... 2515!> \param my_group_L_starts_im_time ... 2516!> \param my_group_L_ends_im_time ... 2517!> \param my_group_L_sizes_im_time ... 2518!> \param sab_orb_sub ... 2519!> \param sab_orb_all ... 2520!> \param para_env ... 2521!> \param num_diff_blk ... 2522!> \param local_atoms_for_mao_basis ... 2523! ************************************************************************************************** 2524 SUBROUTINE reserve_blocks_3c(mat_3c_overl_int, mat_munu, qs_env, dimen_RI, & 2525 cut_RI, unit_nr, my_group_L_starts_im_time, my_group_L_ends_im_time, & 2526 my_group_L_sizes_im_time, sab_orb_sub, sab_orb_all, para_env, & 2527 num_diff_blk, local_atoms_for_mao_basis) 2528 2529 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int 2530 TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu 2531 TYPE(qs_environment_type), POINTER :: qs_env 2532 INTEGER, INTENT(IN) :: dimen_RI, cut_RI, unit_nr 2533 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: my_group_L_starts_im_time, & 2534 my_group_L_ends_im_time, & 2535 my_group_L_sizes_im_time 2536 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 2537 POINTER :: sab_orb_sub, sab_orb_all 2538 TYPE(cp_para_env_type), POINTER :: para_env 2539 INTEGER, INTENT(OUT) :: num_diff_blk 2540 INTEGER, ALLOCATABLE, DIMENSION(:), & 2541 INTENT(INOUT), OPTIONAL :: local_atoms_for_mao_basis 2542 2543 CHARACTER(LEN=*), PARAMETER :: routineN = 'reserve_blocks_3c', & 2544 routineP = moduleN//':'//routineN 2545 2546 INTEGER :: acell, atom_RI, atom_RI_end, atom_RI_start, bcell, blk, blk_cnt, col, handle, & 2547 handle2, i_cut_RI, i_img, i_img_outer, i_loop, iatom, iatom_outer, ikind, iset, isgf_RI, & 2548 j_img, j_img_outer, jatom, jatom_outer, jkind, jset, kind_RI, LLL, LLL_set_end, & 2549 LLL_set_start, maxval_1, maxval_2, maxval_3, minval_1, minval_2, minval_3, & 2550 my_group_L_end, my_group_L_size, my_group_L_start, natom, nblkrows_total, nco_RI, ncoa, & 2551 ncob, nimg, nkind, nset_RI, nseta, nsetb, num_loops, row, set_RI, sgf_RI, sgfa, sgfb 2552 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_RI_index, kind_of, tmp 2553 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: old_col, old_row 2554 INTEGER, DIMENSION(3) :: cell_vec, outer_cell_vec 2555 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, col_blk_sizes_scaled, l_max_RI, l_min_RI, & 2556 la_max, la_min, lb_max, lb_min, npgf_RI, npgfa, npgfb, nsgf_RI, nsgfa, nsgfb, & 2557 row_blk_end, row_blk_offset, row_blk_sizes, row_blk_sizes_scaled, row_blk_start 2558 INTEGER, DIMENSION(:, :), POINTER :: first_sgf_RI, first_sgfa, first_sgfb 2559 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 2560 LOGICAL :: do_kpoints_cubic_RPA, new_block 2561 REAL(KIND=dp) :: dab, daRI, dbRI 2562 REAL(KIND=dp), DIMENSION(3) :: rab, rab_outer, raRI, rbRI 2563 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_RI 2564 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_RI, rpgfa, rpgfb, sphi_a, sphi_b, & 2565 sphi_RI, zet_RI, zeta, zetb 2566 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2567 TYPE(cell_type), POINTER :: cell 2568 TYPE(dbcsr_distribution_type) :: dist_sub 2569 TYPE(dft_control_type), POINTER :: dft_control 2570 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_RI_tmp 2571 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, basis_set_RI 2572 TYPE(kpoint_type), POINTER :: kpoints 2573 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 2574 TYPE(neighbor_list_iterator_p_type), & 2575 DIMENSION(:), POINTER :: nl_iterator, nl_iterator_outer 2576 TYPE(one_dim_int_array), ALLOCATABLE, & 2577 DIMENSION(:, :) :: cols_to_alloc, rows_to_alloc 2578 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 2579 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 2580 2581 CALL timeset(routineN, handle) 2582 2583 NULLIFY (qs_kind_set, atomic_kind_set, cell, particle_set, molecule_kind_set, basis_set_RI, basis_set_a, basis_set_b, & 2584 nl_iterator) 2585 2586 do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints 2587 2588 ! get stuff 2589 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, dft_control=dft_control, & 2590 cell=cell, particle_set=particle_set, molecule_kind_set=molecule_kind_set, & 2591 nkind=nkind, natom=natom, kpoints=kpoints) 2592 2593 IF (do_kpoints_cubic_RPA) THEN 2594 ! set up new kpoint type with periodic images according to eps_grid from MP2 section 2595 ! instead of eps_pgf_orb from QS section 2596 CALL kpoint_init_cell_index(kpoints, sab_orb_sub, para_env, dft_control) 2597 2598 nimg = dft_control%nimages 2599 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index) 2600 2601 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 2602 "3C_OVERLAP_INTEGRALS_INFO| Number of periodic images considered:", nimg 2603 2604 num_diff_blk = 0 2605 ELSE 2606 nimg = 1 2607 END IF 2608 2609 ! allocate matrix for storing 3c overlap integrals 2610 NULLIFY (mat_3c_overl_int) 2611 CALL dbcsr_allocate_matrix_set(mat_3c_overl_int, cut_RI, nimg, nimg) 2612 2613 ALLOCATE (row_blk_start(natom)) 2614 ALLOCATE (row_blk_end(natom)) 2615 2616 ALLOCATE (basis_set_RI_tmp(nkind)) 2617 CALL basis_set_list_setup(basis_set_RI_tmp, "RI_AUX", qs_kind_set) 2618 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, & 2619 basis=basis_set_RI_tmp) 2620 DEALLOCATE (basis_set_RI_tmp) 2621 2622 ALLOCATE (atom_from_RI_index(dimen_RI)) 2623 2624 DO LLL = 1, dimen_RI 2625 2626 DO iatom = 1, natom 2627 2628 IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN 2629 2630 atom_from_RI_index(LLL) = iatom 2631 2632 END IF 2633 2634 END DO 2635 2636 END DO 2637 2638 CALL dbcsr_get_info(mat_munu%matrix, & 2639 row_blk_offset=row_blk_offset, & 2640 row_blk_size=row_blk_sizes, & 2641 nblkrows_total=nblkrows_total) 2642 2643 IF (do_kpoints_cubic_RPA) THEN 2644 2645 minval_1 = MINVAL(kpoints%index_to_cell(1, :)) 2646 maxval_1 = MAXVAL(kpoints%index_to_cell(1, :)) 2647 minval_2 = MINVAL(kpoints%index_to_cell(2, :)) 2648 maxval_2 = MAXVAL(kpoints%index_to_cell(2, :)) 2649 minval_3 = MINVAL(kpoints%index_to_cell(3, :)) 2650 maxval_3 = MAXVAL(kpoints%index_to_cell(3, :)) 2651 2652 END IF 2653 2654 ALLOCATE (kind_of(natom)) 2655 ALLOCATE (rows_to_alloc(nimg, nimg)) 2656 ALLOCATE (cols_to_alloc(nimg, nimg)) 2657 2658 ALLOCATE (old_row(nimg, nimg)) 2659 old_row = -1 2660 ALLOCATE (old_col(nimg, nimg)) 2661 old_col = -1 2662 2663 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of) 2664 2665 DO i_cut_RI = 1, cut_RI 2666 2667 DO i_img = 1, nimg 2668 DO j_img = 1, nimg 2669 ALLOCATE (rows_to_alloc(i_img, j_img)%array(1)) 2670 rows_to_alloc(i_img, j_img)%array(1) = 0 2671 ALLOCATE (cols_to_alloc(i_img, j_img)%array(1)) 2672 cols_to_alloc(i_img, j_img)%array(1) = 0 2673 END DO 2674 END DO 2675 2676 my_group_L_start = my_group_L_starts_im_time(i_cut_RI) 2677 my_group_L_end = my_group_L_ends_im_time(i_cut_RI) 2678 my_group_L_size = my_group_L_sizes_im_time(i_cut_RI) 2679 2680 atom_RI_start = atom_from_RI_index(my_group_L_start) 2681 atom_RI_end = atom_from_RI_index(my_group_L_end) 2682 2683 DO atom_RI = atom_RI_start, atom_RI_end 2684 2685 kind_RI = kind_of(atom_RI) 2686 CALL get_qs_kind(qs_kind=qs_kind_set(kind_RI), basis_set=basis_set_RI, basis_type="RI_AUX") 2687 2688 first_sgf_RI => basis_set_RI%first_sgf 2689 2690 l_max_RI => basis_set_RI%lmax 2691 l_min_RI => basis_set_RI%lmin 2692 npgf_RI => basis_set_RI%npgf 2693 nset_RI = basis_set_RI%nset 2694 nsgf_RI => basis_set_RI%nsgf_set 2695 rpgf_RI => basis_set_RI%pgf_radius 2696 set_radius_RI => basis_set_RI%set_radius 2697 sphi_RI => basis_set_RI%sphi 2698 zet_RI => basis_set_RI%zet 2699 2700 CALL neighbor_list_iterator_create(nl_iterator_outer, sab_orb_all) 2701 DO WHILE (neighbor_list_iterate(nl_iterator_outer) == 0) 2702 2703 CALL get_iterator_info(nl_iterator_outer, & 2704 iatom=iatom_outer, jatom=jatom_outer, r=rab_outer, cell=outer_cell_vec) 2705 2706 IF (iatom_outer .NE. atom_RI .AND. jatom_outer .NE. atom_RI) CYCLE 2707 2708 CALL neighbor_list_iterator_create(nl_iterator, sab_orb_sub) 2709 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 2710 2711 CALL get_iterator_info(nl_iterator, & 2712 iatom=iatom, jatom=jatom, r=rab, cell=cell_vec) 2713 2714 DO i_img_outer = 1, nimg 2715 2716 DO j_img_outer = 1, nimg 2717 2718 IF (i_img_outer > j_img_outer) CYCLE 2719 IF (i_img_outer < j_img_outer) THEN 2720 num_loops = 2 2721 ELSE IF (i_img_outer == j_img_outer) THEN 2722 num_loops = 1 2723 END IF 2724 2725 DO i_loop = 1, num_loops 2726 2727 IF (i_loop == 1) THEN 2728 i_img = i_img_outer 2729 j_img = j_img_outer 2730 ELSE IF (i_loop == 2) THEN 2731 i_img = j_img_outer 2732 j_img = i_img_outer 2733 END IF 2734 2735 IF (atom_RI .EQ. iatom_outer .AND. & 2736 iatom .NE. jatom_outer .AND. & 2737 jatom .NE. jatom_outer) & 2738 CYCLE 2739 2740 IF (atom_RI .EQ. jatom_outer .AND. & 2741 iatom .NE. iatom_outer .AND. & 2742 jatom .NE. iatom_outer) & 2743 CYCLE 2744 2745 IF (iatom .EQ. iatom_outer .AND. & 2746 atom_RI .EQ. jatom_outer) THEN 2747 2748 IF (do_kpoints_cubic_RPA) THEN 2749 2750 IF (outer_cell_vec(1) < minval_1) CYCLE 2751 IF (outer_cell_vec(1) > maxval_1) CYCLE 2752 IF (outer_cell_vec(2) < minval_2) CYCLE 2753 IF (outer_cell_vec(2) > maxval_2) CYCLE 2754 IF (outer_cell_vec(3) < minval_3) CYCLE 2755 IF (outer_cell_vec(3) > maxval_3) CYCLE 2756 2757 IF (-outer_cell_vec(1) + cell_vec(1) < minval_1) CYCLE 2758 IF (-outer_cell_vec(1) + cell_vec(1) > maxval_1) CYCLE 2759 IF (-outer_cell_vec(2) + cell_vec(2) < minval_2) CYCLE 2760 IF (-outer_cell_vec(2) + cell_vec(2) > maxval_2) CYCLE 2761 IF (-outer_cell_vec(3) + cell_vec(3) < minval_3) CYCLE 2762 IF (-outer_cell_vec(3) + cell_vec(3) > maxval_3) CYCLE 2763 2764 acell = cell_to_index(-outer_cell_vec(1), -outer_cell_vec(2), -outer_cell_vec(3)) 2765 IF (.NOT. (acell == i_img)) CYCLE 2766 bcell = cell_to_index(-outer_cell_vec(1) + cell_vec(1), & 2767 -outer_cell_vec(2) + cell_vec(2), & 2768 -outer_cell_vec(3) + cell_vec(3)) 2769 IF (.NOT. (bcell == j_img)) CYCLE 2770 2771 END IF 2772 2773 raRI = rab_outer 2774 rbRI = raRI - rab 2775 2776 ELSE IF (atom_RI .EQ. iatom_outer .AND. & 2777 iatom .EQ. jatom_outer) THEN 2778 2779 IF (do_kpoints_cubic_RPA) THEN 2780 ! for kpoints we have the non-symmetric neighbor list 2781 CYCLE 2782 END IF 2783 2784 raRI = -rab_outer 2785 rbRI = raRI - rab 2786 2787 ELSE IF (jatom .EQ. iatom_outer .AND. & 2788 atom_RI .EQ. jatom_outer) THEN 2789 2790 IF (do_kpoints_cubic_RPA) THEN 2791 ! for kpoints we have the non-symmetric neighbor list 2792 CYCLE 2793 END IF 2794 2795 rbRI = rab_outer 2796 raRI = rbRI + rab 2797 2798 ELSE IF (atom_RI .EQ. iatom_outer .AND. & 2799 jatom .EQ. jatom_outer) THEN 2800 2801 IF (do_kpoints_cubic_RPA) THEN 2802 ! for kpoints we have the non-symmetric neighbor list 2803 CYCLE 2804 END IF 2805 2806 rbRI = -rab_outer 2807 raRI = rbRI + rab 2808 2809 END IF 2810 2811 IF (PRESENT(local_atoms_for_mao_basis)) THEN 2812 local_atoms_for_mao_basis(iatom) = iatom 2813 local_atoms_for_mao_basis(jatom) = jatom 2814 END IF 2815 2816 dab = SQRT(SUM(rab**2)) 2817 daRI = SQRT(SUM(raRI**2)) 2818 dbRI = SQRT(SUM(rbRI**2)) 2819 2820 ikind = kind_of(iatom) 2821 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a) 2822 first_sgfa => basis_set_a%first_sgf 2823 la_max => basis_set_a%lmax 2824 la_min => basis_set_a%lmin 2825 npgfa => basis_set_a%npgf 2826 nseta = basis_set_a%nset 2827 nsgfa => basis_set_a%nsgf_set 2828 rpgfa => basis_set_a%pgf_radius 2829 set_radius_a => basis_set_a%set_radius 2830 sphi_a => basis_set_a%sphi 2831 zeta => basis_set_a%zet 2832 2833 jkind = kind_of(jatom) 2834 CALL get_qs_kind(qs_kind=qs_kind_set(jkind), basis_set=basis_set_b) 2835 first_sgfb => basis_set_b%first_sgf 2836 lb_max => basis_set_b%lmax 2837 lb_min => basis_set_b%lmin 2838 npgfb => basis_set_b%npgf 2839 nsetb = basis_set_b%nset 2840 nsgfb => basis_set_b%nsgf_set 2841 rpgfb => basis_set_b%pgf_radius 2842 set_radius_b => basis_set_b%set_radius 2843 sphi_b => basis_set_b%sphi 2844 zetb => basis_set_b%zet 2845 2846 DO set_RI = 1, nset_RI 2847 2848 DO iset = 1, nseta 2849 2850 IF (set_radius_a(iset) + set_radius_RI(set_RI) < daRI) CYCLE 2851 2852 DO jset = 1, nsetb 2853 2854 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 2855 IF (set_radius_b(jset) + set_radius_RI(set_RI) < dbRI) CYCLE 2856 2857 nco_RI = npgf_RI(set_RI)*ncoset(l_max_RI(set_RI)) 2858 ncoa = npgfa(iset)*ncoset(la_max(iset)) 2859 ncob = npgfb(jset)*ncoset(lb_max(jset)) 2860 2861 sgf_RI = first_sgf_RI(1, set_RI) 2862 sgfa = first_sgfa(1, iset) 2863 sgfb = first_sgfb(1, jset) 2864 2865 LLL_set_start = 1 + sgf_RI - 1 - my_group_L_start + 1 + row_blk_start(atom_RI) - 1 2866 2867 IF (LLL_set_start > my_group_L_size) CYCLE 2868 2869 LLL_set_end = nsgf_RI(set_RI) + sgf_RI - 1 - my_group_L_start + 1 + row_blk_start(atom_RI) - 1 2870 2871 IF (LLL_set_end < 1) CYCLE 2872 2873 IF (ncoa*ncob*nco_RI > 0) THEN 2874 2875 DO isgf_RI = 1, nsgf_RI(set_RI) 2876 2877 LLL = isgf_RI + sgf_RI - 1 - my_group_L_start + 1 + row_blk_start(atom_RI) - 1 2878 2879 IF (LLL < 1) CYCLE 2880 IF (LLL > my_group_L_size) CYCLE 2881 2882 IF (iatom < jatom) THEN 2883 2884 row = iatom 2885 col = jatom 2886 2887 ELSE IF (iatom >= jatom) THEN 2888 2889 row = jatom 2890 col = iatom 2891 2892 END IF 2893 2894 IF (SIZE(rows_to_alloc(i_img_outer, j_img_outer)%array) == 1 .AND. & 2895 rows_to_alloc(i_img_outer, j_img_outer)%array(1) == 0) THEN 2896 2897 IF (row == col) THEN 2898 2899 rows_to_alloc(i_img_outer, j_img_outer)%array = row 2900 cols_to_alloc(i_img_outer, j_img_outer)%array = col 2901 2902 ELSE 2903 2904 DEALLOCATE (rows_to_alloc(i_img_outer, j_img_outer)%array) 2905 DEALLOCATE (cols_to_alloc(i_img_outer, j_img_outer)%array) 2906 2907 ALLOCATE (rows_to_alloc(i_img_outer, j_img_outer)%array(2)) 2908 ALLOCATE (cols_to_alloc(i_img_outer, j_img_outer)%array(2)) 2909 2910 rows_to_alloc(i_img_outer, j_img_outer)%array(1) = row 2911 rows_to_alloc(i_img_outer, j_img_outer)%array(2) = col 2912 2913 cols_to_alloc(i_img_outer, j_img_outer)%array(1) = col 2914 cols_to_alloc(i_img_outer, j_img_outer)%array(2) = row 2915 2916 END IF 2917 2918 old_row(i_img_outer, j_img_outer) = row 2919 old_col(i_img_outer, j_img_outer) = col 2920 2921 ELSE 2922 2923 new_block = .FALSE. 2924 2925 IF ((old_row(i_img_outer, j_img_outer) .NE. row) .OR. & 2926 (old_col(i_img_outer, j_img_outer) .NE. col)) THEN 2927 2928 blk_cnt = SIZE(rows_to_alloc(i_img_outer, j_img_outer)%array) 2929 2930 new_block = .TRUE. 2931 2932 DO blk = 1, blk_cnt 2933 2934 IF (rows_to_alloc(i_img_outer, j_img_outer)%array(blk) == row .AND. & 2935 cols_to_alloc(i_img_outer, j_img_outer)%array(blk) == col) THEN 2936 2937 new_block = .FALSE. 2938 EXIT 2939 2940 END IF 2941 2942 END DO 2943 2944 IF (new_block .AND. row == col) THEN 2945 2946 ALLOCATE (tmp(blk_cnt)) 2947 tmp(1:blk_cnt) = rows_to_alloc(i_img_outer, j_img_outer)%array(1:blk_cnt) 2948 DEALLOCATE (rows_to_alloc(i_img_outer, j_img_outer)%array) 2949 ALLOCATE (rows_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 1)) 2950 rows_to_alloc(i_img_outer, j_img_outer)%array(1:blk_cnt) = tmp(1:blk_cnt) 2951 rows_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 1) = row 2952 tmp(1:blk_cnt) = cols_to_alloc(i_img_outer, j_img_outer)%array(1:blk_cnt) 2953 DEALLOCATE (cols_to_alloc(i_img_outer, j_img_outer)%array) 2954 ALLOCATE (cols_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 1)) 2955 cols_to_alloc(i_img_outer, j_img_outer)%array(1:blk_cnt) = tmp(1:blk_cnt) 2956 cols_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 1) = col 2957 DEALLOCATE (tmp) 2958 2959 old_row(i_img_outer, j_img_outer) = row 2960 old_col(i_img_outer, j_img_outer) = col 2961 2962 ELSE IF (new_block .AND. row .NE. col) THEN 2963 2964 ALLOCATE (tmp(blk_cnt)) 2965 tmp(1:blk_cnt) = rows_to_alloc(i_img_outer, j_img_outer)%array(1:blk_cnt) 2966 DEALLOCATE (rows_to_alloc(i_img_outer, j_img_outer)%array) 2967 ALLOCATE (rows_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 2)) 2968 rows_to_alloc(i_img_outer, j_img_outer)%array(1:blk_cnt) = tmp(1:blk_cnt) 2969 rows_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 1) = row 2970 rows_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 2) = col 2971 tmp(1:blk_cnt) = cols_to_alloc(i_img_outer, j_img_outer)%array(1:blk_cnt) 2972 DEALLOCATE (cols_to_alloc(i_img_outer, j_img_outer)%array) 2973 ALLOCATE (cols_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 2)) 2974 cols_to_alloc(i_img_outer, j_img_outer)%array(1:blk_cnt) = tmp(1:blk_cnt) 2975 cols_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 1) = col 2976 cols_to_alloc(i_img_outer, j_img_outer)%array(blk_cnt + 2) = row 2977 DEALLOCATE (tmp) 2978 2979 old_row(i_img_outer, j_img_outer) = row 2980 old_col(i_img_outer, j_img_outer) = col 2981 2982 END IF 2983 2984 END IF 2985 2986 END IF 2987 2988 END DO ! single RI basis functions 2989 2990 END IF ! number of triples > 0 2991 2992 END DO ! jset 2993 2994 END DO ! iset 2995 2996 END DO ! i_loop for i_img and j_img for kpoints 2997 2998 END DO ! j_img 2999 3000 END DO ! i_img 3001 3002 END DO ! set RI 3003 3004 END DO ! inner neighbor list 3005 3006 CALL neighbor_list_iterator_release(nl_iterator) 3007 3008 END DO ! outer neighbor list 3009 3010 CALL neighbor_list_iterator_release(nl_iterator_outer) 3011 3012 END DO ! atom_RI 3013 3014 DO i_img = 1, nimg 3015 3016 DO j_img = 1, nimg 3017 3018 IF (i_img > j_img) CYCLE 3019 3020 CALL timeset(routineN//"_reserve_dbcsr", handle2) 3021 3022 blk_cnt = SIZE(rows_to_alloc(i_img, j_img)%array) 3023 3024 CALL dbcsr_get_info(mat_munu%matrix, & 3025 row_blk_size=row_blk_sizes, & 3026 col_blk_size=col_blk_sizes, & 3027 distribution=dist_sub) 3028 3029 ALLOCATE (col_blk_sizes_scaled(SIZE(col_blk_sizes))) 3030 col_blk_sizes_scaled(:) = col_blk_sizes(:)*my_group_L_size 3031 3032 ALLOCATE (mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix) 3033 3034 CALL dbcsr_create(matrix=mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix, & 3035 name="(munuP)", & 3036 dist=dist_sub, & 3037 matrix_type=dbcsr_type_no_symmetry, & 3038 row_blk_size=row_blk_sizes, & 3039 col_blk_size=col_blk_sizes_scaled, & 3040 nze=0) 3041 3042 IF (rows_to_alloc(i_img, j_img)%array(1) .NE. 0) THEN 3043 3044 ! that should be okay for the kpoints sinceif block (row, col) is allocated, also (col, row) 3045 ! afterwards, the matrix will be filtered anyway 3046 CALL dbcsr_reserve_blocks(mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix, & 3047 rows=rows_to_alloc(i_img, j_img)%array(1:blk_cnt), & 3048 cols=cols_to_alloc(i_img, j_img)%array(1:blk_cnt)) 3049 3050 END IF 3051 3052 CALL dbcsr_finalize(mat_3c_overl_int(i_cut_RI, i_img, j_img)%matrix) 3053 3054 IF (i_img .NE. j_img) THEN 3055 3056 ALLOCATE (row_blk_sizes_scaled(SIZE(row_blk_sizes))) 3057 row_blk_sizes_scaled(:) = row_blk_sizes(:)*my_group_L_size 3058 3059 ALLOCATE (mat_3c_overl_int(i_cut_RI, j_img, i_img)%matrix) 3060 3061 CALL dbcsr_create(matrix=mat_3c_overl_int(i_cut_RI, j_img, i_img)%matrix, & 3062 name="(munuP)", & 3063 dist=dist_sub, & 3064 matrix_type=dbcsr_type_no_symmetry, & 3065 row_blk_size=col_blk_sizes, & 3066 col_blk_size=row_blk_sizes_scaled, & 3067 nze=0) 3068 3069 DEALLOCATE (row_blk_sizes_scaled) 3070 3071 IF (rows_to_alloc(i_img, j_img)%array(1) .NE. 0) THEN 3072 CALL dbcsr_reserve_blocks(mat_3c_overl_int(i_cut_RI, j_img, i_img)%matrix, & 3073 rows=cols_to_alloc(i_img, j_img)%array(1:blk_cnt), & 3074 cols=rows_to_alloc(i_img, j_img)%array(1:blk_cnt)) 3075 END IF 3076 CALL dbcsr_finalize(mat_3c_overl_int(i_cut_RI, j_img, i_img)%matrix) 3077 END IF 3078 3079 DEALLOCATE (col_blk_sizes_scaled) 3080 3081 CALL timestop(handle2) 3082 3083 END DO ! j_img 3084 3085 END DO ! i_img 3086 3087 DO i_img = 1, nimg 3088 DO j_img = 1, nimg 3089 DEALLOCATE (rows_to_alloc(i_img, j_img)%array) 3090 DEALLOCATE (cols_to_alloc(i_img, j_img)%array) 3091 END DO 3092 END DO 3093 3094 END DO ! cut_RI 3095 3096 DEALLOCATE (rows_to_alloc, cols_to_alloc) 3097 DEALLOCATE (atom_from_RI_index) 3098 DEALLOCATE (row_blk_start, row_blk_end) 3099 3100 CALL timestop(handle) 3101 3102 END SUBROUTINE reserve_blocks_3c 3103 3104! ************************************************************************************************** 3105!> \brief ... 3106!> \param mao_coeff_repl ... 3107!> \param mao_coeff ... 3108!> \param local_atoms_for_mao_basis ... 3109!> \param natom ... 3110!> \param para_env ... 3111! ************************************************************************************************** 3112 SUBROUTINE replicate_mao_coeff(mao_coeff_repl, mao_coeff, local_atoms_for_mao_basis, natom, para_env) 3113 TYPE(two_dim_real_array), ALLOCATABLE, & 3114 DIMENSION(:) :: mao_coeff_repl 3115 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coeff 3116 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: local_atoms_for_mao_basis 3117 INTEGER, INTENT(IN) :: natom 3118 TYPE(cp_para_env_type), POINTER :: para_env 3119 3120 CHARACTER(LEN=*), PARAMETER :: routineN = 'replicate_mao_coeff', & 3121 routineP = moduleN//':'//routineN 3122 3123 INTEGER :: atom_rec, block_size, col, col_size, end_msg, handle, iatom, imepos, mao_size, & 3124 num_atom_from_imepos, prim_size, row, row_size, start_msg 3125 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_counter, entry_counter, & 3126 num_entries_atoms_rec, & 3127 num_entries_atoms_send, sizes_rec, & 3128 sizes_send 3129 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes_mao, row_blk_sizes_prim 3130 INTEGER, DIMENSION(:, :), POINTER :: req_array 3131 REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block 3132 TYPE(dbcsr_iterator_type) :: iter 3133 TYPE(integ_mat_buffer_type), ALLOCATABLE, & 3134 DIMENSION(:) :: buffer_rec, buffer_send 3135 3136 CALL timeset(routineN, handle) 3137 3138 CALL dbcsr_get_info(mao_coeff(1)%matrix, & 3139 row_blk_size=row_blk_sizes_prim, & 3140 col_blk_size=col_blk_sizes_mao) 3141 3142 ALLOCATE (mao_coeff_repl(1:natom)) 3143 3144 ALLOCATE (num_entries_atoms_send(0:2*para_env%num_pe - 1)) 3145 num_entries_atoms_send(:) = 0 3146 3147 ALLOCATE (num_entries_atoms_rec(0:2*para_env%num_pe - 1)) 3148 num_entries_atoms_rec(:) = 0 3149 3150 DO iatom = 1, natom 3151 3152 ! check whether local atom is there 3153 IF (local_atoms_for_mao_basis(iatom) == iatom) THEN 3154 ! row index belongs to primary basis, col index belongs to mao basis 3155 prim_size = row_blk_sizes_prim(iatom) 3156 mao_size = col_blk_sizes_mao(iatom) 3157 ALLOCATE (mao_coeff_repl(iatom)%array(prim_size, mao_size)) 3158 3159 CALL dbcsr_get_stored_coordinates(mao_coeff(1)%matrix, iatom, iatom, imepos) 3160 3161 num_entries_atoms_rec(2*imepos) = num_entries_atoms_rec(2*imepos) + prim_size*mao_size 3162 num_entries_atoms_rec(2*imepos + 1) = num_entries_atoms_rec(2*imepos + 1) + 1 3163 3164 END IF 3165 3166 END DO 3167 3168 IF (para_env%num_pe > 1) THEN 3169 CALL mp_alltoall(num_entries_atoms_rec, num_entries_atoms_send, 2, para_env%group) 3170 ELSE 3171 num_entries_atoms_send(0:1) = num_entries_atoms_rec(0:1) 3172 END IF 3173 3174 ALLOCATE (buffer_rec(0:para_env%num_pe - 1)) 3175 ALLOCATE (buffer_send(0:para_env%num_pe - 1)) 3176 3177 ! allocate data message and corresponding indices 3178 DO imepos = 0, para_env%num_pe - 1 3179 3180 ALLOCATE (buffer_rec(imepos)%msg(num_entries_atoms_rec(2*imepos))) 3181 buffer_rec(imepos)%msg = 0.0_dp 3182 3183 ALLOCATE (buffer_send(imepos)%msg(num_entries_atoms_send(2*imepos))) 3184 buffer_send(imepos)%msg = 0.0_dp 3185 3186 ALLOCATE (buffer_rec(imepos)%indx(num_entries_atoms_rec(2*imepos + 1), 6)) 3187 buffer_rec(imepos)%indx = 0 3188 3189 ALLOCATE (buffer_send(imepos)%indx(num_entries_atoms_send(2*imepos + 1), 6)) 3190 buffer_send(imepos)%indx = 0 3191 3192 END DO 3193 3194 ALLOCATE (entry_counter(0:para_env%num_pe - 1)) 3195 entry_counter(:) = 0 3196 3197 ALLOCATE (atom_counter(0:para_env%num_pe - 1)) 3198 atom_counter = 0 3199 3200 DO iatom = 1, natom 3201 3202 ! check whether local atom is there 3203 IF (local_atoms_for_mao_basis(iatom) == iatom) THEN 3204 3205 CALL dbcsr_get_stored_coordinates(mao_coeff(1)%matrix, iatom, iatom, imepos) 3206 3207 atom_counter(imepos) = atom_counter(imepos) + 1 3208 3209 buffer_rec(imepos)%indx(atom_counter(imepos), 1) = iatom 3210 3211 END IF 3212 3213 END DO 3214 3215 ALLOCATE (req_array(1:para_env%num_pe, 4)) 3216 3217 ALLOCATE (sizes_rec(0:para_env%num_pe - 1)) 3218 ALLOCATE (sizes_send(0:para_env%num_pe - 1)) 3219 3220 DO imepos = 0, para_env%num_pe - 1 3221 3222 sizes_send(imepos) = num_entries_atoms_send(2*imepos) 3223 sizes_rec(imepos) = num_entries_atoms_rec(2*imepos) 3224 3225 END DO 3226 3227 ! change rec and send! 3228 CALL communicate_buffer(para_env, sizes_send, sizes_rec, buffer_send, buffer_rec, req_array, do_msg=.FALSE.) 3229 3230 atom_counter(:) = 0 3231 3232 CALL dbcsr_iterator_start(iter, mao_coeff(1)%matrix) 3233 3234 DO WHILE (dbcsr_iterator_blocks_left(iter)) 3235 3236 CALL dbcsr_iterator_next_block(iter, row, col, data_block, & 3237 row_size=row_size, col_size=col_size) 3238 3239 ! for the mao trafo matrix, only the diagonal blocks should be allocated 3240 CPASSERT(row == col) 3241 3242 DO imepos = 0, para_env%num_pe - 1 3243 3244 IF (ANY(buffer_send(imepos)%indx(:, 1) == row)) THEN 3245 3246 block_size = row_size*col_size 3247 start_msg = entry_counter(imepos) + 1 3248 end_msg = entry_counter(imepos) + block_size 3249 3250 buffer_send(imepos)%msg(start_msg:end_msg) = RESHAPE(data_block(1:row_size, 1:col_size), (/block_size/)) 3251 3252 entry_counter(imepos) = entry_counter(imepos) + block_size 3253 atom_counter(imepos) = atom_counter(imepos) + 1 3254 3255 buffer_send(imepos)%indx(atom_counter(imepos), 2) = row 3256 buffer_send(imepos)%indx(atom_counter(imepos), 3) = row_size 3257 buffer_send(imepos)%indx(atom_counter(imepos), 4) = col_size 3258 buffer_send(imepos)%indx(atom_counter(imepos), 5) = start_msg 3259 buffer_send(imepos)%indx(atom_counter(imepos), 6) = end_msg 3260 3261 END IF 3262 3263 END DO 3264 3265 END DO ! blocks 3266 3267 CALL dbcsr_iterator_stop(iter) 3268 3269 ! communicate the mao coefficients 3270 CALL communicate_buffer(para_env, sizes_rec, sizes_send, buffer_rec, buffer_send, req_array) 3271 3272 DEALLOCATE (req_array, sizes_rec, sizes_send) 3273 3274 ! fill mao_coeff_repl 3275 DO imepos = 0, para_env%num_pe - 1 3276 3277 num_atom_from_imepos = num_entries_atoms_rec(2*imepos + 1) 3278 3279 DO atom_rec = 1, num_atom_from_imepos 3280 3281 iatom = buffer_rec(imepos)%indx(atom_rec, 2) 3282 row_size = buffer_rec(imepos)%indx(atom_rec, 3) 3283 col_size = buffer_rec(imepos)%indx(atom_rec, 4) 3284 start_msg = buffer_rec(imepos)%indx(atom_rec, 5) 3285 end_msg = buffer_rec(imepos)%indx(atom_rec, 6) 3286 3287 mao_coeff_repl(iatom)%array(1:row_size, 1:col_size) = & 3288 RESHAPE(buffer_rec(imepos)%msg(start_msg:end_msg), (/row_size, col_size/)) 3289 3290 END DO 3291 3292 END DO 3293 3294 DEALLOCATE (num_entries_atoms_send, num_entries_atoms_rec, entry_counter, atom_counter) 3295 3296 DO imepos = 0, para_env%num_pe - 1 3297 DEALLOCATE (buffer_rec(imepos)%msg) 3298 DEALLOCATE (buffer_rec(imepos)%indx) 3299 DEALLOCATE (buffer_send(imepos)%msg) 3300 DEALLOCATE (buffer_send(imepos)%indx) 3301 END DO 3302 3303 DEALLOCATE (buffer_rec, buffer_send) 3304 3305 CALL timestop(handle) 3306 3307 END SUBROUTINE replicate_mao_coeff 3308 3309!*************************************************************************************************** 3310! ************************************************************************************************** 3311!> \brief ... 3312!> \param mao_coeff_repl_occ ... 3313!> \param mao_coeff_repl_virt ... 3314!> \param local_atoms_for_mao_basis ... 3315!> \param natom ... 3316! ************************************************************************************************** 3317 SUBROUTINE clean_up(mao_coeff_repl_occ, mao_coeff_repl_virt, local_atoms_for_mao_basis, natom) 3318 TYPE(two_dim_real_array), ALLOCATABLE, & 3319 DIMENSION(:), INTENT(INOUT) :: mao_coeff_repl_occ, mao_coeff_repl_virt 3320 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: local_atoms_for_mao_basis 3321 INTEGER, INTENT(IN) :: natom 3322 3323 INTEGER :: iatom 3324 3325 DO iatom = 1, natom 3326 ! check whether local atom is there 3327 IF (local_atoms_for_mao_basis(iatom) == iatom) THEN 3328 DEALLOCATE (mao_coeff_repl_occ(iatom)%array) 3329 END IF 3330 END DO 3331 3332 DO iatom = 1, natom 3333 ! check whether local atom is there 3334 IF (local_atoms_for_mao_basis(iatom) == iatom) THEN 3335 DEALLOCATE (mao_coeff_repl_virt(iatom)%array) 3336 END IF 3337 END DO 3338 3339 DEALLOCATE (mao_coeff_repl_occ, mao_coeff_repl_virt) 3340 3341 DEALLOCATE (local_atoms_for_mao_basis) 3342 3343 END SUBROUTINE clean_up 3344 3345! ************************************************************************************************** 3346!> \brief ... 3347!> \param my_group_L_starts_im_time ... 3348!> \param my_group_L_ends_im_time ... 3349!> \param my_group_L_sizes_im_time ... 3350!> \param dimen_RI ... 3351!> \param ngroup ... 3352!> \param cut_memory ... 3353!> \param cut_RI ... 3354!> \param color_sub ... 3355!> \param qs_env ... 3356! ************************************************************************************************** 3357 SUBROUTINE setup_group_L_im_time(my_group_L_starts_im_time, my_group_L_ends_im_time, my_group_L_sizes_im_time, & 3358 dimen_RI, ngroup, cut_memory, cut_RI, color_sub, qs_env) 3359 3360 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: my_group_L_starts_im_time, & 3361 my_group_L_ends_im_time, & 3362 my_group_L_sizes_im_time 3363 INTEGER, INTENT(IN) :: dimen_RI, ngroup, cut_memory 3364 INTEGER, INTENT(OUT) :: cut_RI 3365 INTEGER, INTENT(IN) :: color_sub 3366 TYPE(qs_environment_type), POINTER :: qs_env 3367 3368 CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_group_L_im_time', & 3369 routineP = moduleN//':'//routineN 3370 3371 INTEGER :: cut_RI_old, handle, iblk_RI, icut, & 3372 itmp(2), natom, nblks_RI, nkind, & 3373 size_RI, start_RI 3374 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_group_L_ends_im_time_tmp, & 3375 my_group_L_starts_im_time_tmp, & 3376 row_blk_end_RI, row_blk_offset_RI 3377 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes_RI 3378 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 3379 TYPE(group_dist_d1_type) :: gd_array_mem_cut 3380 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ri_aux 3381 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 3382 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 3383 3384 CALL timeset(routineN, handle) 3385 3386 CALL create_group_dist(gd_array_mem_cut, cut_memory, dimen_RI) 3387 3388 ALLOCATE (my_group_L_starts_im_time(1:cut_memory)) 3389 my_group_L_starts_im_time = 0 3390 ALLOCATE (my_group_L_ends_im_time(1:cut_memory)) 3391 my_group_L_ends_im_time = 0 3392 3393 DO icut = 0, cut_memory - 1 3394 3395 CALL get_group_dist(gd_array_mem_cut, icut, sizes=size_RI, starts=start_RI) 3396 3397 itmp = get_limit(size_RI, ngroup, color_sub) 3398 3399 my_group_L_starts_im_time(icut + 1) = itmp(1) + start_RI - 1 3400 my_group_L_ends_im_time(icut + 1) = itmp(2) + start_RI - 1 3401 3402 END DO 3403 3404 CALL get_qs_env(qs_env, & 3405 qs_kind_set=qs_kind_set, & 3406 atomic_kind_set=atomic_kind_set, & 3407 particle_set=particle_set) 3408 3409 natom = SIZE(particle_set) 3410 nkind = SIZE(atomic_kind_set) 3411 3412 ALLOCATE (row_blk_sizes_RI(natom)) 3413 3414 ALLOCATE (basis_set_ri_aux(nkind)) 3415 CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set) 3416 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes_RI, basis=basis_set_ri_aux) 3417 DEALLOCATE (basis_set_ri_aux) 3418 3419 ! check whether my_group_L_im_time covers different atoms. In this case, cut my_group_L_im_time 3420 ! especially improves memory a bit 3421 nblks_RI = SIZE(row_blk_sizes_RI) 3422 3423 ALLOCATE (row_blk_offset_RI(nblks_RI)) 3424 ALLOCATE (row_blk_end_RI(nblks_RI)) 3425 3426 row_blk_offset_RI(1) = 1 3427 DO iblk_RI = 2, nblks_RI 3428 row_blk_offset_RI(iblk_RI) = row_blk_offset_RI(iblk_RI - 1) + row_blk_sizes_RI(iblk_RI - 1) 3429 END DO 3430 3431 row_blk_end_RI(nblks_RI) = dimen_RI 3432 DO iblk_RI = 1, nblks_RI - 1 3433 row_blk_end_RI(iblk_RI) = row_blk_offset_RI(iblk_RI + 1) - 1 3434 END DO 3435 3436 cut_RI = cut_memory 3437 cut_RI_old = 0 3438 3439 DO WHILE (cut_RI .NE. cut_RI_old) 3440 3441 cut_RI_old = cut_RI 3442 3443 DO iblk_RI = 1, nblks_RI 3444 3445 DO icut = 1, cut_RI_old 3446 3447 IF (row_blk_offset_RI(iblk_RI) >= my_group_L_starts_im_time(icut) .AND. & 3448 row_blk_offset_RI(iblk_RI) <= my_group_L_ends_im_time(icut) .AND. & 3449 row_blk_end_RI(iblk_RI) < my_group_L_ends_im_time(icut)) THEN 3450 3451 cut_RI = cut_RI + 1 3452 3453 ALLOCATE (my_group_L_starts_im_time_tmp(cut_RI - 1)) 3454 ALLOCATE (my_group_L_ends_im_time_tmp(cut_RI - 1)) 3455 3456 my_group_L_starts_im_time_tmp(:) = my_group_L_starts_im_time(:) 3457 my_group_L_ends_im_time_tmp(:) = my_group_L_ends_im_time(:) 3458 3459 DEALLOCATE (my_group_L_starts_im_time) 3460 DEALLOCATE (my_group_L_ends_im_time) 3461 3462 ALLOCATE (my_group_L_starts_im_time(cut_RI)) 3463 ALLOCATE (my_group_L_ends_im_time(cut_RI)) 3464 3465 my_group_L_starts_im_time(1:icut) = my_group_L_starts_im_time_tmp(1:icut) 3466 my_group_L_starts_im_time(icut + 1) = row_blk_offset_RI(iblk_RI + 1) 3467 IF (cut_RI >= icut + 2) THEN 3468 my_group_L_starts_im_time(icut + 2:cut_RI) = my_group_L_starts_im_time_tmp(icut + 1:cut_RI - 1) 3469 END IF 3470 3471 IF (icut - 1 >= 1) THEN 3472 my_group_L_ends_im_time(1:icut - 1) = my_group_L_ends_im_time_tmp(1:icut - 1) 3473 END IF 3474 my_group_L_ends_im_time(icut) = row_blk_offset_RI(iblk_RI + 1) - 1 3475 my_group_L_ends_im_time(icut + 1:cut_RI) = my_group_L_ends_im_time_tmp(icut:cut_RI - 1) 3476 3477 DEALLOCATE (my_group_L_starts_im_time_tmp) 3478 DEALLOCATE (my_group_L_ends_im_time_tmp) 3479 3480 END IF 3481 3482 END DO 3483 3484 END DO 3485 3486 END DO 3487 3488 ALLOCATE (my_group_L_sizes_im_time(cut_RI)) 3489 my_group_L_sizes_im_time(:) = my_group_L_ends_im_time(:) - my_group_L_starts_im_time(:) + 1 3490 3491 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time_util(cut_memory)) 3492 3493 qs_env%mp2_env%ri_rpa_im_time_util(1)%cut_RI = cut_RI 3494 3495 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time_util(1)%my_group_L_starts_im_time(1:cut_RI)) 3496 qs_env%mp2_env%ri_rpa_im_time_util(1)%my_group_L_starts_im_time(:) = my_group_L_starts_im_time(:) 3497 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time_util(1)%my_group_L_ends_im_time(1:cut_RI)) 3498 qs_env%mp2_env%ri_rpa_im_time_util(1)%my_group_L_ends_im_time(:) = my_group_L_ends_im_time(:) 3499 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time_util(1)%my_group_L_sizes_im_time(1:cut_RI)) 3500 qs_env%mp2_env%ri_rpa_im_time_util(1)%my_group_L_sizes_im_time(:) = my_group_L_sizes_im_time(:) 3501 3502 CALL release_group_dist(gd_array_mem_cut) 3503 DEALLOCATE (row_blk_sizes_RI) 3504 3505 CALL timestop(handle) 3506 3507 END SUBROUTINE setup_group_L_im_time 3508 3509END MODULE mp2_integrals 3510