1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 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 atomic_kind_types, ONLY: atomic_kind_type 14 USE basis_set_types, ONLY: gto_basis_set_p_type 15 USE bibliography, ONLY: DelBen2013,& 16 cite_reference 17 USE cell_types, ONLY: cell_type,& 18 get_cell 19 USE cp_blacs_env, ONLY: cp_blacs_env_type 20 USE cp_control_types, ONLY: dft_control_type 21 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 22 cp_dbcsr_m_by_n_from_template 23 USE cp_eri_mme_interface, ONLY: cp_eri_mme_param,& 24 cp_eri_mme_set_params 25 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 26 cp_fm_struct_release,& 27 cp_fm_struct_type 28 USE cp_fm_types, ONLY: cp_fm_create,& 29 cp_fm_get_info,& 30 cp_fm_p_type,& 31 cp_fm_release,& 32 cp_fm_type 33 USE cp_para_env, ONLY: cp_para_env_release 34 USE cp_para_types, ONLY: cp_para_env_type 35 USE cp_units, ONLY: cp_unit_from_cp2k 36 USE dbcsr_api, ONLY: & 37 dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, dbcsr_p_type, dbcsr_release, & 38 dbcsr_release_p, dbcsr_scalar, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, & 39 dbcsr_type_real_8 40 USE dbcsr_tensor_api, ONLY: & 41 dbcsr_t_clear, dbcsr_t_contract, dbcsr_t_copy, dbcsr_t_create, dbcsr_t_destroy, & 42 dbcsr_t_distribution_destroy, dbcsr_t_distribution_new, dbcsr_t_distribution_type, & 43 dbcsr_t_filter, dbcsr_t_get_block, dbcsr_t_get_stored_coordinates, & 44 dbcsr_t_mp_environ_pgrid, dbcsr_t_pgrid_create, dbcsr_t_pgrid_create_expert, & 45 dbcsr_t_pgrid_destroy, dbcsr_t_pgrid_type, dbcsr_t_put_block, dbcsr_t_reserve_blocks, & 46 dbcsr_t_split_blocks, dbcsr_t_type 47 USE group_dist_types, ONLY: create_group_dist,& 48 get_group_dist,& 49 group_dist_d1_type 50 USE input_constants, ONLY: do_eri_gpw,& 51 do_eri_mme,& 52 do_eri_os,& 53 do_potential_coulomb,& 54 do_potential_id,& 55 do_potential_long,& 56 do_potential_short,& 57 do_potential_truncated 58 USE kinds, ONLY: dp 59 USE kpoint_methods, ONLY: kpoint_init_cell_index 60 USE kpoint_types, ONLY: kpoint_type 61 USE libint_2c_3c, ONLY: libint_potential_type 62 USE machine, ONLY: m_flush 63 USE message_passing, ONLY: mp_cart_create,& 64 mp_dims_create,& 65 mp_environ,& 66 mp_max,& 67 mp_sendrecv,& 68 mp_sync 69 USE mp2_eri, ONLY: mp2_eri_3c_integrate 70 USE mp2_eri_gpw, ONLY: cleanup_gpw,& 71 mp2_eri_3c_integrate_gpw,& 72 prepare_gpw 73 USE mp2_ri_2c, ONLY: get_2c_integrals 74 USE particle_methods, ONLY: get_particle_set 75 USE particle_types, ONLY: particle_type 76 USE pw_env_types, ONLY: pw_env_type 77 USE pw_poisson_types, ONLY: pw_poisson_type 78 USE pw_pool_types, ONLY: pw_pool_type 79 USE pw_types, ONLY: pw_p_type 80 USE qs_environment_types, ONLY: get_qs_env,& 81 qs_environment_type,& 82 set_qs_env 83 USE qs_integral_utils, ONLY: basis_set_list_setup 84 USE qs_kind_types, ONLY: qs_kind_type 85 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 86 USE qs_tensors, ONLY: build_3c_integrals,& 87 build_3c_neighbor_lists,& 88 neighbor_list_3c_destroy 89 USE qs_tensors_types, ONLY: create_3c_tensor,& 90 create_tensor_batches,& 91 distribution_3d_create,& 92 distribution_3d_type,& 93 neighbor_list_3c_type,& 94 pgf_block_sizes 95 USE task_list_types, ONLY: task_list_type 96 USE util, ONLY: get_limit 97 98!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num 99#include "./base/base_uses.f90" 100 101 IMPLICIT NONE 102 103 PRIVATE 104 105 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_integrals' 106 107 PUBLIC :: mp2_ri_gpw_compute_in 108 109CONTAINS 110 111! ************************************************************************************************** 112!> \brief with ri mp2 gpw 113!> \param BIb_C ... 114!> \param BIb_C_gw ... 115!> \param BIb_C_bse_ij ... 116!> \param BIb_C_bse_ab ... 117!> \param gd_array ... 118!> \param gd_B_virtual ... 119!> \param dimen_RI ... 120!> \param dimen_RI_red ... 121!> \param qs_env ... 122!> \param para_env ... 123!> \param para_env_sub ... 124!> \param color_sub ... 125!> \param cell ... 126!> \param particle_set ... 127!> \param atomic_kind_set ... 128!> \param qs_kind_set ... 129!> \param mo_coeff ... 130!> \param fm_matrix_L_RI_metric ... 131!> \param nmo ... 132!> \param homo ... 133!> \param mat_munu ... 134!> \param sab_orb_sub ... 135!> \param mo_coeff_o ... 136!> \param mo_coeff_v ... 137!> \param mo_coeff_all ... 138!> \param mo_coeff_gw ... 139!> \param eps_filter ... 140!> \param unit_nr ... 141!> \param mp2_memory ... 142!> \param calc_PQ_cond_num ... 143!> \param calc_forces ... 144!> \param blacs_env_sub ... 145!> \param my_do_gw ... 146!> \param do_bse ... 147!> \param gd_B_all ... 148!> \param starts_array_mc ... 149!> \param ends_array_mc ... 150!> \param starts_array_mc_block ... 151!> \param ends_array_mc_block ... 152!> \param gw_corr_lev_occ ... 153!> \param gw_corr_lev_virt ... 154!> \param do_im_time ... 155!> \param do_kpoints_cubic_RPA ... 156!> \param kpoints ... 157!> \param t_3c_M ... 158!> \param t_3c_O ... 159!> \param ri_metric ... 160!> \param gd_B_occ_bse ... 161!> \param gd_B_virt_bse ... 162!> \param BIb_C_beta ... 163!> \param BIb_C_gw_beta ... 164!> \param gd_B_virtual_beta ... 165!> \param homo_beta ... 166!> \param mo_coeff_o_beta ... 167!> \param mo_coeff_v_beta ... 168!> \param mo_coeff_all_beta ... 169!> \param mo_coeff_gw_beta ... 170!> \author Mauro Del Ben 171! ************************************************************************************************** 172 SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd_array, gd_B_virtual, & 173 dimen_RI, dimen_RI_red, qs_env, para_env, para_env_sub, color_sub, & 174 cell, particle_set, atomic_kind_set, qs_kind_set, mo_coeff, & 175 fm_matrix_L_RI_metric, nmo, homo, & 176 mat_munu, & 177 sab_orb_sub, mo_coeff_o, mo_coeff_v, mo_coeff_all, & 178 mo_coeff_gw, eps_filter, unit_nr, & 179 mp2_memory, calc_PQ_cond_num, calc_forces, blacs_env_sub, my_do_gw, do_bse, & 180 gd_B_all, starts_array_mc, ends_array_mc, & 181 starts_array_mc_block, ends_array_mc_block, & 182 gw_corr_lev_occ, gw_corr_lev_virt, & 183 do_im_time, do_kpoints_cubic_RPA, kpoints, & 184 t_3c_M, t_3c_O, & 185 ri_metric, gd_B_occ_bse, gd_B_virt_bse, BIb_C_beta, BIb_C_gw_beta, & 186 gd_B_virtual_beta, homo_beta, mo_coeff_o_beta, mo_coeff_v_beta, & 187 mo_coeff_all_beta, mo_coeff_gw_beta) 188 189 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & 190 INTENT(OUT) :: BIb_C, BIb_C_gw, BIb_C_bse_ij, & 191 BIb_C_bse_ab 192 TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array, gd_B_virtual 193 INTEGER, INTENT(OUT) :: dimen_RI, dimen_RI_red 194 TYPE(qs_environment_type), POINTER :: qs_env 195 TYPE(cp_para_env_type), POINTER :: para_env, para_env_sub 196 INTEGER, INTENT(IN) :: color_sub 197 TYPE(cell_type), POINTER :: cell 198 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 199 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 200 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 201 TYPE(cp_fm_type), POINTER :: mo_coeff 202 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_matrix_L_RI_metric 203 INTEGER, INTENT(IN) :: nmo, homo 204 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu 205 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 206 INTENT(IN), POINTER :: sab_orb_sub 207 TYPE(dbcsr_type), POINTER :: mo_coeff_o, mo_coeff_v, mo_coeff_all, & 208 mo_coeff_gw 209 REAL(KIND=dp), INTENT(IN) :: eps_filter 210 INTEGER, INTENT(IN) :: unit_nr 211 REAL(KIND=dp), INTENT(IN) :: mp2_memory 212 LOGICAL, INTENT(IN) :: calc_PQ_cond_num, calc_forces 213 TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub 214 LOGICAL, INTENT(IN) :: my_do_gw, do_bse 215 TYPE(group_dist_d1_type), INTENT(OUT) :: gd_B_all 216 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: starts_array_mc, ends_array_mc, & 217 starts_array_mc_block, & 218 ends_array_mc_block 219 INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt 220 LOGICAL, INTENT(IN) :: do_im_time, do_kpoints_cubic_RPA 221 TYPE(kpoint_type), POINTER :: kpoints 222 TYPE(dbcsr_t_type), INTENT(OUT) :: t_3c_M 223 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :), & 224 INTENT(OUT) :: t_3c_O 225 TYPE(libint_potential_type), INTENT(IN) :: ri_metric 226 TYPE(group_dist_d1_type), INTENT(OUT) :: gd_B_occ_bse, gd_B_virt_bse 227 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & 228 INTENT(OUT), OPTIONAL :: BIb_C_beta, BIb_C_gw_beta 229 TYPE(group_dist_d1_type), INTENT(OUT), OPTIONAL :: gd_B_virtual_beta 230 INTEGER, INTENT(IN), OPTIONAL :: homo_beta 231 TYPE(dbcsr_type), OPTIONAL, POINTER :: mo_coeff_o_beta, mo_coeff_v_beta, & 232 mo_coeff_all_beta, mo_coeff_gw_beta 233 234 CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_ri_gpw_compute_in', & 235 routineP = moduleN//':'//routineN 236 237 INTEGER :: cm, cut_memory, cut_memory_int, eri_method, gw_corr_lev_total, handle, handle2, & 238 handle4, i, i_counter, itmp(2), j, jcell, kcell, LLL, max_row_col_local, & 239 max_row_col_local_beta, max_row_col_local_gw, max_row_col_local_occ_bse, & 240 max_row_col_local_virt_bse, min_bsize, mp_comm_t3c_2, my_B_all_end, my_B_all_size, & 241 my_B_all_start, my_B_occ_bse_end, my_B_occ_bse_size, my_B_occ_bse_start, my_B_size, & 242 my_B_size_beta, my_B_virt_bse_end, my_B_virt_bse_size, my_B_virt_bse_start, & 243 my_B_virtual_end, my_B_virtual_end_beta, my_B_virtual_start, my_B_virtual_start_beta, & 244 my_group_L_end 245 INTEGER :: my_group_L_size, my_group_L_start, natom, ngroup, ngroup_im_time, & 246 ngroup_im_time_P, nimg, nkind, num_small_eigen, potential_type, ri_metric_type, virtual, & 247 virtual_beta 248 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_AO_1, dist_AO_2, dist_RI, & 249 ends_array_mc_block_int, ends_array_mc_int, sizes_AO, sizes_AO_split, sizes_RI, & 250 sizes_RI_split, starts_array_mc_block_int, starts_array_mc_int, sub_proc_map 251 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, local_col_row_info_beta, & 252 local_col_row_info_gw, local_col_row_info_occ_bse, local_col_row_info_virt_bse 253 INTEGER, DIMENSION(3) :: pcoord, pdims, pdims_t3c, periodic 254 LOGICAL :: do_alpha_beta, do_gpw, & 255 do_kpoints_from_Gamma, do_svd, & 256 impose_split, memory_info 257 REAL(KIND=dp) :: cond_num, cutoff_old, eps_svd, & 258 mem_for_iaK, omega_pot, rc_ang, & 259 relative_cutoff_old 260 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old 261 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: my_Lrows 262 TYPE(cp_eri_mme_param), POINTER :: eri_param 263 TYPE(cp_fm_type), POINTER :: fm_BIb_bse_ab, fm_BIb_bse_ij, fm_BIb_gw, & 264 fm_BIb_gw_beta, fm_BIb_jb, & 265 fm_BIb_jb_beta, fm_matrix_L 266 TYPE(cp_para_env_type), POINTER :: para_env_L 267 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_munu_local_L 268 TYPE(dbcsr_t_pgrid_type) :: pgrid_t3c_M, pgrid_t3c_overl 269 TYPE(dbcsr_t_type) :: t_3c_overl_int_template, t_3c_tmp 270 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_overl_int 271 TYPE(dbcsr_type) :: matrix_bse_ab, matrix_bse_anu, matrix_bse_ij, matrix_bse_inu, & 272 matrix_ia_jb, matrix_ia_jb_beta, matrix_ia_jnu, matrix_ia_jnu_beta, matrix_in_jm, & 273 matrix_in_jm_beta, matrix_in_jnu, matrix_in_jnu_beta 274 TYPE(dft_control_type), POINTER :: dft_control 275 TYPE(distribution_3d_type) :: dist_3d 276 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ao, basis_set_ri_aux 277 TYPE(neighbor_list_3c_type) :: nl_3c 278 TYPE(pw_env_type), POINTER :: pw_env_sub 279 TYPE(pw_p_type) :: pot_g, psi_L, rho_g, rho_r 280 TYPE(pw_poisson_type), POINTER :: poisson_env 281 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 282 TYPE(task_list_type), POINTER :: task_list_sub 283 284 CALL timeset(routineN, handle) 285 286 CALL cite_reference(DelBen2013) 287 288 do_alpha_beta = .FALSE. 289 IF (PRESENT(BIb_C_beta) .AND. & 290 PRESENT(gd_B_virtual_beta) .AND. & 291 PRESENT(homo_beta) .AND. & 292 PRESENT(mo_coeff_o_beta) .AND. & 293 PRESENT(mo_coeff_v_beta)) do_alpha_beta = .TRUE. 294 295 IF ((ri_metric%potential_type == do_potential_short .OR. ri_metric%potential_type == do_potential_truncated) & 296 .AND. .NOT. do_im_time) THEN 297 CPABORT("TRUNCATED and SHORTRANGE RI metrics not yet implemented") 298 ENDIF 299 300 virtual = nmo - homo 301 gw_corr_lev_total = gw_corr_lev_virt + gw_corr_lev_occ 302 303 eri_method = qs_env%mp2_env%eri_method 304 eri_param => qs_env%mp2_env%eri_mme_param 305 do_svd = qs_env%mp2_env%do_svd 306 eps_svd = qs_env%mp2_env%eps_svd 307 potential_type = qs_env%mp2_env%potential_parameter%potential_type 308 ri_metric_type = ri_metric%potential_type 309 omega_pot = qs_env%mp2_env%potential_parameter%omega 310 311 ! whether we need gpw integrals (plus pw stuff) 312 do_gpw = (eri_method == do_eri_gpw) .OR. & 313 ((potential_type == do_potential_long .OR. ri_metric_type == do_potential_long) & 314 .AND. qs_env%mp2_env%eri_method == do_eri_os) & 315 .OR. (ri_metric_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme) 316 317 IF (do_svd .AND. calc_forces) THEN 318 CPABORT("SVD not implemented for forces.!") 319 END IF 320 321 do_kpoints_from_Gamma = SUM(qs_env%mp2_env%ri_rpa_im_time%kp_grid) > 0 322 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 323 CALL get_qs_env(qs_env=qs_env, & 324 kpoints=kpoints) 325 END IF 326 IF (do_kpoints_from_Gamma) THEN 327 CALL compute_kpoints(qs_env, kpoints, unit_nr) 328 END IF 329 330 ! in case of imaginary time, we do not need the intermediate matrices 331 IF (.NOT. do_im_time) THEN 332 333 CALL create_intermediate_matrices(matrix_ia_jnu, matrix_ia_jb, mo_coeff_o, virtual, homo, & 334 fm_BIb_jb, "fm_BIb_jb", max_row_col_local, & 335 blacs_env_sub, para_env_sub, local_col_row_info) 336 337 CALL create_group_dist(gd_B_virtual, para_env_sub%num_pe, virtual) 338 CALL get_group_dist(gd_B_virtual, para_env_sub%mepos, my_B_virtual_start, my_B_virtual_end, my_B_size) 339 340 IF (do_alpha_beta) THEN 341 342 virtual_beta = nmo - homo_beta 343 344 CALL create_intermediate_matrices(matrix_ia_jnu_beta, matrix_ia_jb_beta, mo_coeff_o_beta, & 345 virtual_beta, homo_beta, & 346 fm_BIb_jb_beta, "fm_BIb_jb_beta", & 347 max_row_col_local_beta, & 348 blacs_env_sub, para_env_sub, local_col_row_info_beta) 349 350 CALL create_group_dist(gd_B_virtual_beta, para_env_sub%num_pe, virtual_beta) 351 CALL get_group_dist(gd_B_virtual_beta, para_env_sub%mepos, my_B_virtual_start_beta, my_B_virtual_end_beta, & 352 my_B_size_beta) 353 354 END IF 355 356 ! in the case of G0W0, we need (K|nm), n,m may be occ or virt (m restricted to corrected levels) 357 IF (my_do_gw) THEN 358 359 CALL create_intermediate_matrices(matrix_in_jnu, matrix_in_jm, mo_coeff_gw, & 360 nmo, gw_corr_lev_total, & 361 fm_BIb_gw, "fm_BIb_gw", & 362 max_row_col_local_gw, & 363 blacs_env_sub, para_env_sub, local_col_row_info_gw) 364 365 CALL create_group_dist(gd_B_all, para_env_sub%num_pe, nmo) 366 CALL get_group_dist(gd_B_all, para_env_sub%mepos, my_B_all_start, my_B_all_end, my_B_all_size) 367 368 IF (do_alpha_beta) THEN 369 ! deallocate local_col_row_info_gw, otherwise it gets twice allocated in create_intermediate_m 370 DEALLOCATE (local_col_row_info_gw) 371 CALL create_intermediate_matrices(matrix_in_jnu_beta, matrix_in_jm_beta, mo_coeff_gw_beta, & 372 nmo, gw_corr_lev_total, & 373 fm_BIb_gw_beta, "fm_BIb_gw_beta", & 374 max_row_col_local_gw, & 375 blacs_env_sub, para_env_sub, local_col_row_info_gw) 376 377 ! we don"t need parallelization arrays for beta since the matrix sizes of B_nm^P is the same 378 ! for the beta case and therefore the parallelization of beta is the same than for alpha 379 380 END IF 381 END IF 382 END IF ! not for imaginary time 383 384 IF (do_bse) THEN 385 386 CPASSERT(my_do_gw) 387 CPASSERT(.NOT. do_im_time) 388 ! GPW integrals have to be implemented later 389 CPASSERT(.NOT. (eri_method == do_eri_gpw)) 390 391 ! virt x virt matrices 392 CALL create_intermediate_matrices(matrix_bse_anu, matrix_bse_ab, mo_coeff_v, virtual, virtual, & 393 fm_BIb_bse_ab, "fm_BIb_bse_ab", max_row_col_local_virt_bse, & 394 blacs_env_sub, para_env_sub, local_col_row_info_virt_bse) 395 396 CALL create_group_dist(gd_B_virt_bse, para_env_sub%num_pe, virtual) 397 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) 398 399 ! occ x occ matrices 400 CALL create_intermediate_matrices(matrix_bse_inu, matrix_bse_ij, mo_coeff_o, homo, homo, & 401 fm_BIb_bse_ij, "fm_BIb_bse_ij", max_row_col_local_occ_bse, & 402 blacs_env_sub, para_env_sub, local_col_row_info_occ_bse) 403 404 CALL create_group_dist(gd_B_occ_bse, para_env_sub%num_pe, homo) 405 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) 406 407 END IF 408 409 ngroup = para_env%num_pe/para_env_sub%num_pe 410 411 ! Preparations for MME method to compute ERIs 412 IF (qs_env%mp2_env%eri_method .EQ. do_eri_mme) THEN 413 ! cell might have changed, so we need to reset parameters 414 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) 415 ENDIF 416 417 CALL get_cell(cell=cell, periodic=periodic) 418 ! for minimax Ewald summation, full periodicity is required 419 IF (eri_method == do_eri_mme) THEN 420 CPASSERT(periodic(1) == 1 .AND. periodic(2) == 1 .AND. periodic(3) == 1) 421 END IF 422 423 IF (do_svd .AND. (do_kpoints_from_Gamma .OR. do_kpoints_cubic_RPA)) THEN 424 CPABORT("SVD with kpoints not implemented yet!") 425 END IF 426 427 CALL get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, & 428 fm_matrix_L, ngroup, color_sub, dimen_RI, dimen_RI_red, kpoints, mo_coeff, & 429 my_group_L_size, my_group_L_start, my_group_L_end, & 430 gd_array, calc_PQ_cond_num .AND. .NOT. do_svd, cond_num, do_svd, eps_svd, & 431 num_small_eigen, qs_env%mp2_env%potential_parameter, ri_metric, & 432 fm_matrix_L_RI_metric, & 433 do_im_time, do_kpoints_from_Gamma .OR. do_kpoints_cubic_RPA, qs_env%mp2_env%mp2_gpw%eps_pgf_orb_S, & 434 qs_kind_set, sab_orb_sub) 435 436 IF (unit_nr > 0) THEN 437 ASSOCIATE (ri_metric=>qs_env%mp2_env%ri_metric) 438 SELECT CASE (ri_metric%potential_type) 439 CASE (do_potential_coulomb) 440 WRITE (unit_nr, FMT="(/T3,A,T74,A)") & 441 "RI_INFO| RI metric: ", "COULOMB" 442 CASE (do_potential_short) 443 WRITE (unit_nr, FMT="(T3,A,T71,A)") & 444 "RI_INFO| RI metric: ", "SHORTRANGE" 445 WRITE (unit_nr, '(T3,A,T61,F20.10)') & 446 "RI_INFO| Omega: ", ri_metric%omega 447 rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom") 448 WRITE (unit_nr, '(T3,A,T61,F20.10)') & 449 "RI_INFO| Cutoff Radius [angstrom]: ", rc_ang 450 CASE (do_potential_long) 451 WRITE (unit_nr, FMT="(T3,A,T72,A)") & 452 "RI_INFO| RI metric: ", "LONGRANGE" 453 WRITE (unit_nr, '(T3,A,T61,F20.10)') & 454 "RI_INFO| Omega: ", ri_metric%omega 455 CASE (do_potential_id) 456 WRITE (unit_nr, FMT="(T3,A,T73,A)") & 457 "RI_INFO| RI metric: ", "IDENTITY" 458 CASE (do_potential_truncated) 459 WRITE (unit_nr, FMT="(T3,A,T72,A)") & 460 "RI_INFO| RI metric: ", "TRUNCATED" 461 rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom") 462 WRITE (unit_nr, '(T3,A,T61,F20.10)') & 463 "RI_INFO| Cutoff Radius [angstrom]: ", rc_ang 464 END SELECT 465 END ASSOCIATE 466 ENDIF 467 468 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 469 "RI_INFO| Cholesky decomposition group size:", para_env_L%num_pe 470 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 471 "RI_INFO| Number of groups for auxiliary basis functions", ngroup 472 IF (calc_PQ_cond_num .OR. do_svd) THEN 473 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T67,ES14.5)") & 474 "RI_INFO| Condition number of the (P|Q):", cond_num 475 IF (do_svd) THEN 476 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,ES8.1,A,T75,i6)") & 477 "RI_INFO| Number of neglected Eigenvalues of (P|Q) smaller than ", eps_svd, ":", num_small_eigen 478 ELSE 479 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,ES14.1,A,T75,i6)") & 480 "RI_INFO| Number of Eigenvalue of (P|Q) smaller ", eps_svd, ":", num_small_eigen 481 END IF 482 END IF 483 484 IF (.NOT. do_im_time) THEN 485 ! replicate the necessary row of the L^{-1} matrix on each proc 486 CALL grep_Lcols(para_env_L, dimen_RI, fm_matrix_L, & 487 my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows) 488 END IF 489 ! clean the L^{-1} matrix 490 CALL cp_fm_release(fm_matrix_L) 491 492 ! in case of imag. time we need the para_env_L later 493 IF (.NOT. do_im_time) THEN 494 CALL cp_para_env_release(para_env_L) 495 END IF 496 497 IF (calc_forces) THEN 498 ! we need (P|Q)^(-1/2) for future use, just save it 499 ! in a fully (home made) distributed way 500 itmp = get_limit(dimen_RI, para_env_sub%num_pe, para_env_sub%mepos) 501 lll = itmp(2) - itmp(1) + 1 502 ALLOCATE (qs_env%mp2_env%ri_grad%PQ_half(lll, my_group_L_size)) 503 qs_env%mp2_env%ri_grad%PQ_half(:, :) = my_Lrows(itmp(1):itmp(2), 1:my_group_L_size) 504 END IF 505 506 IF (unit_nr > 0) THEN 507 WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 508 "RI_INFO| Occupied basis set size:", homo, & 509 "RI_INFO| Virtual basis set size:", virtual, & 510 "RI_INFO| Auxiliary basis set size:", dimen_RI 511 IF (do_svd) THEN 512 WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 513 "RI_INFO| Reduced auxiliary basis set size:", dimen_RI_red 514 END IF 515 516 mem_for_iaK = dimen_RI*REAL(homo, KIND=dp)*virtual*8.0_dp/(1024_dp**2) 517 IF (do_alpha_beta) mem_for_iaK = mem_for_iaK + & 518 dimen_RI*REAL(homo_beta, KIND=dp)*(nmo - homo_beta)*8.0_dp/(1024_dp**2) 519 520 WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for (ia|K) integrals:', & 521 mem_for_iaK, ' MiB' 522 IF (my_do_gw .AND. .NOT. do_im_time) THEN 523 mem_for_iaK = dimen_RI*REAL(nmo, KIND=dp)*gw_corr_lev_total*8.0_dp/(1024_dp**2) 524 525 WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for G0W0-(nm|K) integrals:', & 526 mem_for_iaK, ' MiB' 527 END IF 528 CALL m_flush(unit_nr) 529 ENDIF 530 531 CALL mp_sync(para_env%group) ! sync to see memory output 532 533 ! in case we do imaginary time, we need the overlap tensor (alpha beta P) or trunc. Coulomb tensor 534 IF (.NOT. do_im_time) THEN 535 536 ALLOCATE (sub_proc_map(-para_env_sub%num_pe:2*para_env_sub%num_pe - 1)) 537 DO i = 0, para_env_sub%num_pe - 1 538 sub_proc_map(i) = i 539 sub_proc_map(-i - 1) = para_env_sub%num_pe - i - 1 540 sub_proc_map(para_env_sub%num_pe + i) = i 541 END DO 542 543 ! array that will store the (ia|K) integrals 544 ALLOCATE (BIb_C(my_group_L_size, my_B_size, homo)) 545 BIb_C = 0.0_dp 546 547 IF (do_alpha_beta) THEN 548 ALLOCATE (BIb_C_beta(my_group_L_size, my_B_size_beta, homo_beta)) 549 BIb_C_beta = 0.0_dp 550 END IF 551 552 ! in the case of GW, we also need (nm|K) 553 IF (my_do_gw) THEN 554 555 ALLOCATE (BIb_C_gw(my_group_L_size, my_B_all_size, gw_corr_lev_total)) 556 BIb_C_gw = 0.0_dp 557 IF (do_alpha_beta) THEN 558 ALLOCATE (BIb_C_gw_beta(my_group_L_size, my_B_all_size, gw_corr_lev_total)) 559 BIb_C_gw_beta = 0.0_dp 560 END IF 561 562 END IF 563 564 IF (do_bse) THEN 565 566 ALLOCATE (BIb_C_bse_ij(my_group_L_size, my_B_occ_bse_size, homo)) 567 BIb_C_bse_ij = 0.0_dp 568 569 ALLOCATE (BIb_C_bse_ab(my_group_L_size, my_B_virt_bse_size, virtual)) 570 BIb_C_bse_ab = 0.0_dp 571 572 END IF 573 574 CALL timeset(routineN//"_loop", handle2) 575 576 IF (eri_method == do_eri_mme .AND. & 577 (ri_metric%potential_type == do_potential_coulomb .OR. ri_metric%potential_type == do_potential_long) .OR. & 578 eri_method == do_eri_os .AND. ri_metric%potential_type == do_potential_coulomb) THEN 579 580 NULLIFY (mat_munu_local_L) 581 ALLOCATE (mat_munu_local_L(my_group_L_size)) 582 DO LLL = 1, my_group_L_size 583 NULLIFY (mat_munu_local_L(LLL)%matrix) 584 ALLOCATE (mat_munu_local_L(LLL)%matrix) 585 CALL dbcsr_copy(mat_munu_local_L(LLL)%matrix, mat_munu%matrix) 586 CALL dbcsr_set(mat_munu_local_L(LLL)%matrix, 0.0_dp) 587 ENDDO 588 CALL mp2_eri_3c_integrate(eri_param, ri_metric%potential_type, ri_metric%omega, para_env_sub, qs_env, & 589 first_c=my_group_L_start, last_c=my_group_L_end, & 590 mat_ab=mat_munu_local_L, & 591 basis_type_a="ORB", basis_type_b="ORB", & 592 basis_type_c="RI_AUX", & 593 sab_nl=sab_orb_sub, eri_method=eri_method) 594 595 DO LLL = 1, my_group_L_size 596 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_ia_jnu, matrix_ia_jb, & 597 fm_BIb_jb, BIb_C(LLL, 1:my_B_size, 1:homo), & 598 mo_coeff_o, mo_coeff_v, eps_filter, max_row_col_local, sub_proc_map, & 599 local_col_row_info, my_B_virtual_end, my_B_virtual_start, "alpha") 600 ENDDO 601 CALL contract_B_L(BIb_C, my_Lrows, gd_B_virtual%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 602 ngroup, color_sub, para_env%group, para_env_sub) 603 604 IF (do_alpha_beta) THEN 605 606 DO LLL = 1, my_group_L_size 607 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_ia_jnu_beta, matrix_ia_jb_beta, & 608 fm_BIb_jb_beta, BIb_C_beta(LLL, 1:my_B_size_beta, 1:homo_beta), & 609 mo_coeff_o_beta, mo_coeff_v_beta, eps_filter, max_row_col_local_beta, sub_proc_map, & 610 local_col_row_info_beta, my_B_virtual_end_beta, my_B_virtual_start_beta, "beta") 611 ENDDO 612 CALL contract_B_L(BIb_C_beta, my_Lrows, gd_B_virtual_beta%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 613 ngroup, color_sub, para_env%group, para_env_sub) 614 615 ENDIF 616 617 IF (my_do_gw) THEN 618 619 DO LLL = 1, my_group_L_size 620 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_in_jnu, matrix_in_jm, & 621 fm_BIb_gw, BIb_C_gw(LLL, 1:my_B_all_size, 1:gw_corr_lev_total), & 622 mo_coeff_gw, mo_coeff_all, eps_filter, max_row_col_local_gw, sub_proc_map, & 623 local_col_row_info_gw, my_B_all_end, my_B_all_start, "gw_alpha") 624 ENDDO 625 CALL contract_B_L(BIb_C_gw, my_Lrows, gd_B_all%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 626 ngroup, color_sub, para_env%group, para_env_sub) 627 628 IF (do_alpha_beta) THEN 629 630 DO LLL = 1, my_group_L_size 631 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_in_jnu_beta, matrix_in_jm_beta, & 632 fm_BIb_gw_beta, BIb_C_gw_beta(LLL, 1:my_B_all_size, 1:gw_corr_lev_total), & 633 mo_coeff_gw_beta, mo_coeff_all_beta, eps_filter, max_row_col_local_gw, & 634 sub_proc_map, local_col_row_info_gw, my_B_all_end, my_B_all_start, "gw_beta") 635 ENDDO 636 CALL contract_B_L(BIb_C_gw_beta, my_Lrows, gd_B_all%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 637 ngroup, color_sub, para_env%group, para_env_sub) 638 639 ENDIF 640 ENDIF 641 642 IF (do_bse) THEN 643 644 ! B^ab_P matrix elements for BSE 645 DO LLL = 1, my_group_L_size 646 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_bse_anu, matrix_bse_ab, & 647 fm_BIb_bse_ab, BIb_C_bse_ab(LLL, 1:my_B_virt_bse_size, 1:virtual), & 648 mo_coeff_v, mo_coeff_v, eps_filter, max_row_col_local_virt_bse, & 649 sub_proc_map, local_col_row_info_virt_bse, my_B_all_end, my_B_all_start, "bse_ab") 650 ENDDO 651 CALL contract_B_L(BIb_C_bse_ab, my_Lrows, gd_B_virt_bse%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 652 ngroup, color_sub, para_env%group, para_env_sub) 653 654 ! B^ij_P matrix elements for BSE 655 DO LLL = 1, my_group_L_size 656 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), matrix_bse_inu, matrix_bse_ij, & 657 fm_BIb_bse_ij, BIb_C_bse_ij(LLL, 1:my_B_occ_bse_size, 1:homo), & 658 mo_coeff_o, mo_coeff_o, eps_filter, max_row_col_local_occ_bse, & 659 sub_proc_map, local_col_row_info_occ_bse, & 660 my_B_occ_bse_end, my_B_occ_bse_start, "bse_ij") 661 ENDDO 662 CALL contract_B_L(BIb_C_bse_ij, my_Lrows, gd_B_occ_bse%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, & 663 ngroup, color_sub, para_env%group, para_env_sub) 664 665 END IF 666 667 DO LLL = 1, my_group_L_size 668 CALL dbcsr_release_p(mat_munu_local_L(LLL)%matrix) 669 ENDDO 670 DEALLOCATE (mat_munu_local_L) 671 672 ELSE IF (do_gpw) THEN 673 674 CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, & 675 auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub) 676 677 DO i_counter = 1, my_group_L_size 678 679 CALL mp2_eri_3c_integrate_gpw(mo_coeff, psi_L, rho_g, atomic_kind_set, qs_kind_set, cell, dft_control, & 680 particle_set, pw_env_sub, my_Lrows(:, i_counter), poisson_env, rho_r, pot_g, & 681 ri_metric%potential_type, ri_metric%omega, mat_munu, qs_env, task_list_sub) 682 683 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, matrix_ia_jnu, matrix_ia_jb, & 684 fm_BIb_jb, BIb_C(i_counter, 1:my_B_size, 1:homo), & 685 mo_coeff_o, mo_coeff_v, eps_filter, max_row_col_local, sub_proc_map, & 686 local_col_row_info, my_B_virtual_end, my_B_virtual_start, "alpha") 687 688 IF (do_alpha_beta) THEN 689 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, matrix_ia_jnu_beta, matrix_ia_jb_beta, & 690 fm_BIb_jb_beta, BIb_C_beta(i_counter, 1:my_B_size_beta, 1:homo_beta), & 691 mo_coeff_o_beta, mo_coeff_v_beta, eps_filter, max_row_col_local_beta, sub_proc_map, & 692 local_col_row_info_beta, my_B_virtual_end_beta, my_B_virtual_start_beta, "beta") 693 694 ENDIF 695 696 IF (my_do_gw) THEN 697 ! transform (K|mu nu) to (K|nm), n corresponds to corrected GW levels, m is in nmo 698 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, matrix_in_jnu, matrix_in_jm, & 699 fm_BIb_gw, BIb_C_gw(i_counter, 1:my_B_all_size, 1:gw_corr_lev_total), & 700 mo_coeff_gw, mo_coeff_all, eps_filter, max_row_col_local_gw, sub_proc_map, & 701 local_col_row_info_gw, my_B_all_end, my_B_all_start, "gw_alpha") 702 703 ! the same for beta 704 IF (do_alpha_beta) THEN 705 CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, matrix_in_jnu_beta, matrix_in_jm_beta, & 706 fm_BIb_gw_beta, BIb_C_gw_beta(i_counter, 1:my_B_all_size, 1:gw_corr_lev_total), & 707 mo_coeff_gw_beta, mo_coeff_all_beta, eps_filter, max_row_col_local_gw, & 708 sub_proc_map, local_col_row_info_gw, my_B_all_end, my_B_all_start, "gw_beta") 709 710 END IF 711 END IF 712 713 END DO 714 715 CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, pw_env_sub, & 716 task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L) 717 ELSE 718 CPABORT("Integration method not implemented!") 719 END IF 720 721 CALL timestop(handle2) 722 723 DEALLOCATE (my_Lrows) 724 725 CALL cp_fm_release(fm_BIb_jb) 726 DEALLOCATE (local_col_row_info) 727 728 CALL dbcsr_release(matrix_ia_jnu) 729 CALL dbcsr_release(matrix_ia_jb) 730 IF (do_alpha_beta) THEN 731 CALL dbcsr_release(matrix_ia_jnu_beta) 732 CALL dbcsr_release(matrix_ia_jb_beta) 733 CALL cp_fm_release(fm_BIb_jb_beta) 734 DEALLOCATE (local_col_row_info_beta) 735 END IF 736 737 IF (my_do_gw) THEN 738 CALL dbcsr_release(matrix_in_jnu) 739 CALL dbcsr_release(matrix_in_jm) 740 CALL cp_fm_release(fm_BIb_gw) 741 DEALLOCATE (local_col_row_info_gw) 742 IF (do_alpha_beta) THEN 743 CALL dbcsr_release(matrix_in_jnu_beta) 744 CALL dbcsr_release(matrix_in_jm_beta) 745 CALL cp_fm_release(fm_BIb_gw_beta) 746 END IF 747 END IF 748 749 IF (do_bse) THEN 750 CALL dbcsr_release(matrix_bse_anu) 751 CALL dbcsr_release(matrix_bse_ab) 752 CALL cp_fm_release(fm_BIb_bse_ab) 753 CALL dbcsr_release(matrix_bse_inu) 754 CALL dbcsr_release(matrix_bse_ij) 755 CALL cp_fm_release(fm_BIb_bse_ij) 756 END IF 757 758 DEALLOCATE (sub_proc_map) 759 760 ! imag. time = low-scaling SOS-MP2, RPA, GW 761 ELSE 762 763 impose_split = .NOT. qs_env%mp2_env%ri_rpa_im_time%group_size_internal 764 765 IF (impose_split) THEN 766 ngroup_im_time = para_env%num_pe/qs_env%mp2_env%ri_rpa_im_time%group_size_3c 767 ngroup_im_time_P = para_env%num_pe/qs_env%mp2_env%ri_rpa_im_time%group_size_p 768 ENDIF 769 770 memory_info = qs_env%mp2_env%ri_rpa_im_time%memory_info 771 772 ! we need 3 tensors: 773 ! 1) t_3c_overl_int: 3c overlap integrals, optimized for easy access to integral blocks 774 ! (atomic blocks) 775 ! 2) t_3c_O: 3c overlap integrals, optimized for contraction (split blocks) 776 ! 3) t_3c_M: tensor M, optimized for contraction 777 778 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, dft_control=dft_control) 779 780 IF (.NOT. impose_split) THEN 781 pdims_t3c = 0 782 CALL dbcsr_t_pgrid_create(para_env%group, pdims_t3c, pgrid_t3c_overl) 783 ELSE 784 pgrid_t3c_overl = get_pgrid_from_ngroup(para_env, ngroup_im_time, map1_2d=[1, 2], map2_2d=[3]) 785 ENDIF 786 787 ! set up basis 788 ALLOCATE (sizes_RI(natom), sizes_AO(natom)) 789 ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind)) 790 CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set) 791 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_ri_aux) 792 CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set) 793 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_ao) 794 795 cut_memory_int = qs_env%mp2_env%ri_rpa_im_time%cut_memory 796 CALL create_tensor_batches(sizes_AO, cut_memory_int, starts_array_mc_int, ends_array_mc_int, & 797 starts_array_mc_block_int, ends_array_mc_block_int) 798 799 DEALLOCATE (starts_array_mc_int, ends_array_mc_int) 800 801 CALL create_3c_tensor(t_3c_overl_int_template, dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_overl, & 802 sizes_RI, sizes_AO, sizes_AO, map1=[1, 2], map2=[3], & 803 starts_array_block_2=starts_array_mc_block_int, & 804 ends_array_block_2=ends_array_mc_block_int, & 805 name="O (RI AO | AO)") 806 807 CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set) 808 CALL dbcsr_t_mp_environ_pgrid(pgrid_t3c_overl, pdims, pcoord) 809 CALL mp_cart_create(pgrid_t3c_overl%mp_comm_2d, 3, pdims, pcoord, mp_comm_t3c_2) 810 CALL distribution_3d_create(dist_3d, dist_RI, dist_AO_1, dist_AO_2, & 811 nkind, particle_set, mp_comm_t3c_2, own_comm=.TRUE.) 812 DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2) 813 814 CALL build_3c_neighbor_lists(nl_3c, basis_set_ri_aux, basis_set_ao, basis_set_ao, & 815 dist_3d, ri_metric, "RPA_3c_nl", qs_env, & 816 sym_jk=.NOT. do_kpoints_cubic_RPA, own_dist=.TRUE.) 817 818 ! init k points 819 IF (do_kpoints_cubic_RPA) THEN 820 ! set up new kpoint type with periodic images according to eps_grid from MP2 section 821 ! instead of eps_pgf_orb from QS section 822 CALL kpoint_init_cell_index(kpoints, nl_3c%jk_list, para_env, dft_control) 823 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 824 "3C_OVERLAP_INTEGRALS_INFO| Number of periodic images considered:", dft_control%nimages 825 826 nimg = dft_control%nimages 827 ELSE 828 nimg = 1 829 ENDIF 830 831 ALLOCATE (t_3c_overl_int(nimg, nimg)) 832 833 DO i = 1, SIZE(t_3c_overl_int, 1) 834 DO j = 1, SIZE(t_3c_overl_int, 2) 835 CALL dbcsr_t_create(t_3c_overl_int_template, t_3c_overl_int(i, j)) 836 ENDDO 837 ENDDO 838 839 CALL dbcsr_t_destroy(t_3c_overl_int_template) 840 841 ! split blocks to improve load balancing for tensor contraction 842 min_bsize = qs_env%mp2_env%ri_rpa_im_time%min_bsize 843 844 CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, min_bsize, sizes_AO_split) 845 CALL pgf_block_sizes(atomic_kind_set, basis_set_ri_aux, min_bsize, sizes_RI_split) 846 847 IF (.NOT. impose_split) THEN 848 pdims_t3c = 0 849 CALL dbcsr_t_pgrid_create(para_env%group, pdims_t3c, pgrid_t3c_M) 850 ELSE 851 pgrid_t3c_M = get_pgrid_from_ngroup(para_env, ngroup_im_time_P, map1_2d=[1], map2_2d=[2, 3]) 852 ENDIF 853 854 ASSOCIATE (cut_memory=>qs_env%mp2_env%ri_rpa_im_time%cut_memory) 855 CALL create_tensor_batches(sizes_AO_split, cut_memory, starts_array_mc, ends_array_mc, & 856 starts_array_mc_block, ends_array_mc_block) 857 END ASSOCIATE 858 cut_memory = qs_env%mp2_env%ri_rpa_im_time%cut_memory 859 860 CALL create_3c_tensor(t_3c_M, dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_M, & 861 sizes_RI_split, sizes_AO_split, sizes_AO_split, & 862 map1=[1], map2=[2, 3], & 863 starts_array_block_2=starts_array_mc_block, & 864 ends_array_block_2=ends_array_mc_block, & 865 starts_array_block_3=starts_array_mc_block, & 866 ends_array_block_3=ends_array_mc_block, & 867 name="M (RI | AO AO)") 868 869 DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2) 870 871 CALL dbcsr_t_pgrid_destroy(pgrid_t3c_M) 872 873 ALLOCATE (t_3c_O(SIZE(t_3c_overl_int, 1), SIZE(t_3c_overl_int, 2))) 874 CALL create_3c_tensor(t_3c_O(1, 1), dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_overl, & 875 sizes_RI_split, sizes_AO_split, sizes_AO_split, & 876 map1=[1, 2], map2=[3], & 877 starts_array_block_2=starts_array_mc_block, ends_array_block_2=ends_array_mc_block, & 878 name="O (RI AO | AO)") 879 DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2) 880 881 CALL dbcsr_t_pgrid_destroy(pgrid_t3c_overl) 882 DO i = 1, SIZE(t_3c_O, 1) 883 DO j = 1, SIZE(t_3c_O, 2) 884 IF (i > 1 .OR. j > 1) CALL dbcsr_t_create(t_3c_O(1, 1), t_3c_O(i, j)) 885 ENDDO 886 ENDDO 887 888 ! build integrals in batches and copy to optimized format 889 ! note: integrals are stored in terms of atomic blocks. To avoid a memory bottleneck, 890 ! integrals are calculated in batches and copied to optimized format with subatomic blocks 891 892 DO cm = 1, cut_memory_int 893 CALL build_3c_integrals(t_3c_overl_int, & 894 qs_env%mp2_env%ri_rpa_im_time%eps_filter/2, & 895 qs_env, & 896 nl_3c, & 897 basis_i=basis_set_ri_aux, & 898 basis_j=basis_set_ao, basis_k=basis_set_ao, & 899 potential_parameter=ri_metric, & 900 do_kpoints=do_kpoints_cubic_RPA, & 901 bounds_j=[starts_array_mc_block_int(cm), ends_array_mc_block_int(cm)], desymmetrize=.FALSE.) 902 CALL timeset(routineN//"_copy_3c", handle4) 903 ! copy integral tensor t_3c_overl_int to t_3c_O tensor optimized for contraction 904 DO i = 1, SIZE(t_3c_overl_int, 1) 905 DO j = 1, SIZE(t_3c_overl_int, 2) 906 907 CALL dbcsr_t_copy(t_3c_overl_int(i, j), t_3c_O(i, j), order=[1, 3, 2], & 908 summation=.TRUE., move_data=.TRUE.) 909 CALL dbcsr_t_clear(t_3c_overl_int(i, j)) 910 CALL dbcsr_t_filter(t_3c_O(i, j), qs_env%mp2_env%ri_rpa_im_time%eps_filter/2) 911 ENDDO 912 ENDDO 913 914 CALL timestop(handle4) 915 ENDDO 916 917 DO i = 1, SIZE(t_3c_overl_int, 1) 918 DO j = 1, SIZE(t_3c_overl_int, 2) 919 CALL dbcsr_t_destroy(t_3c_overl_int(i, j)) 920 ENDDO 921 ENDDO 922 DEALLOCATE (t_3c_overl_int) 923 924 CALL timeset(routineN//"_copy_3c", handle4) 925 ! desymmetrize 926 CALL dbcsr_t_create(t_3c_O(1, 1), t_3c_tmp) 927 DO jcell = 1, nimg 928 DO kcell = 1, jcell 929 CALL dbcsr_t_copy(t_3c_O(jcell, kcell), t_3c_tmp) 930 CALL dbcsr_t_copy(t_3c_tmp, t_3c_O(kcell, jcell), order=[1, 3, 2], summation=.TRUE., move_data=.TRUE.) 931 CALL dbcsr_t_filter(t_3c_O(kcell, jcell), qs_env%mp2_env%ri_rpa_im_time%eps_filter) 932 ENDDO 933 ENDDO 934 DO jcell = 1, nimg 935 DO kcell = jcell + 1, nimg 936 CALL dbcsr_t_copy(t_3c_O(jcell, kcell), t_3c_tmp) 937 CALL dbcsr_t_copy(t_3c_tmp, t_3c_O(kcell, jcell), order=[1, 3, 2], summation=.FALSE., move_data=.TRUE.) 938 CALL dbcsr_t_filter(t_3c_O(kcell, jcell), qs_env%mp2_env%ri_rpa_im_time%eps_filter) 939 ENDDO 940 ENDDO 941 CALL dbcsr_t_destroy(t_3c_tmp) 942 943 CALL timestop(handle4) 944 945 DEALLOCATE (basis_set_ri_aux, basis_set_ao) 946 947 CALL neighbor_list_3c_destroy(nl_3c) 948 949 CALL cp_para_env_release(para_env_L) 950 951 END IF 952 953 CALL timestop(handle) 954 955 END SUBROUTINE mp2_ri_gpw_compute_in 956 957! ************************************************************************************************** 958!> \brief heuristic to get a global 3d process grid that is consistent with a number of process subgroups 959!> such that balanced (ideally square) 2d grids on subgroups are obtained 960!> \param para_env ... 961!> \param ngroup number of process groups 962!> \param map1_2d mapping of 3d grid to 2 dimensions (see dbcsr_t_pgrid_create) 963!> \param map2_2d mapping of 3d grid to 2 dimensions (see dbcsr_t_pgrid_create) 964!> \return process grid object compatible with DBCSR tensors 965! ************************************************************************************************** 966 FUNCTION get_pgrid_from_ngroup(para_env, ngroup, map1_2d, map2_2d) RESULT(pgrid) 967 TYPE(cp_para_env_type), INTENT(IN) :: para_env 968 INTEGER, INTENT(IN) :: ngroup 969 INTEGER, DIMENSION(:), INTENT(IN) :: map1_2d, map2_2d 970 TYPE(dbcsr_t_pgrid_type) :: pgrid 971 972 INTEGER :: dimsplit, nproc_sub 973 INTEGER, DIMENSION(2) :: pdims_2_2d, pdims_sub_2d 974 INTEGER, DIMENSION(3) :: pdims_t3c 975 976 nproc_sub = para_env%num_pe/ngroup 977 pdims_sub_2d = 0 978 CALL mp_dims_create(nproc_sub, pdims_sub_2d) 979 pdims_2_2d = 0 980 CALL mp_dims_create(pdims_sub_2d(2)*ngroup, pdims_2_2d) 981 IF (SIZE(map1_2d) == 1) THEN 982 dimsplit = 2 983 pdims_t3c(map1_2d(1)) = pdims_sub_2d(1) 984 pdims_t3c(map2_2d) = pdims_2_2d 985 ELSEIF (SIZE(map2_2d) == 1) THEN 986 dimsplit = 1 987 pdims_t3c(map2_2d(1)) = pdims_sub_2d(1) 988 pdims_t3c(map1_2d) = pdims_2_2d 989 ELSE 990 CPABORT("map1_2d and map2_2d not compatible with a 3d grid") 991 ENDIF 992 CALL dbcsr_t_pgrid_create_expert(para_env%group, pdims_t3c, pgrid, & 993 map1_2d=map1_2d, map2_2d=map2_2d, & 994 nsplit=ngroup, dimsplit=dimsplit) 995 END FUNCTION 996 997! ************************************************************************************************** 998!> \brief Contract (P|ai) = (R|P) x (R|ai) 999!> \param BIb_C (R|ai) 1000!> \param my_Lrows (R|P) 1001!> \param sizes_B number of a (virtual) indices per subgroup process 1002!> \param sizes_L number of P / R (auxiliary) indices per subgroup 1003!> \param blk_size ... 1004!> \param ngroup how many subgroups (NG) 1005!> \param igroup subgroup color 1006!> \param mp_comm communicator 1007!> \param para_env_sub ... 1008! ************************************************************************************************** 1009 SUBROUTINE contract_B_L(BIb_C, my_Lrows, sizes_B, sizes_L, blk_size, ngroup, igroup, mp_comm, para_env_sub) 1010 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: BIb_C 1011 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: my_Lrows 1012 INTEGER, DIMENSION(:), INTENT(IN) :: sizes_B, sizes_L 1013 INTEGER, DIMENSION(2), INTENT(IN) :: blk_size 1014 INTEGER, INTENT(IN) :: ngroup, igroup, mp_comm 1015 TYPE(cp_para_env_type), POINTER :: para_env_sub 1016 1017 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_B_L', routineP = moduleN//':'//routineN 1018 LOGICAL, PARAMETER :: debug = .FALSE. 1019 1020 INTEGER :: check_proc, handle, i, iend, ii, ioff, & 1021 iproc, iproc_glob, istart, loc_a, & 1022 loc_P, nproc, nproc_glob 1023 INTEGER, ALLOCATABLE, DIMENSION(:) :: block_ind_L_P, block_ind_L_R 1024 INTEGER, DIMENSION(1) :: dist_B_i, map_B_1, map_L_1, map_L_2, & 1025 sizes_i 1026 INTEGER, DIMENSION(2) :: map_B_2, pdims_L 1027 INTEGER, DIMENSION(3) :: pdims_B 1028 LOGICAL :: found 1029 INTEGER, DIMENSION(ngroup) :: dist_L_P, dist_L_R 1030 INTEGER, DIMENSION(para_env_sub%num_pe) :: dist_B_a 1031 TYPE(dbcsr_t_distribution_type) :: dist_B, dist_L 1032 TYPE(dbcsr_t_pgrid_type) :: mp_comm_B, mp_comm_L 1033 TYPE(dbcsr_t_type) :: tB_in, tB_in_split, tB_out, & 1034 tB_out_split, tL, tL_split 1035 1036 CALL timeset(routineN, handle) 1037 1038 sizes_i(1) = SIZE(BIb_C, 3) 1039 1040 nproc = para_env_sub%num_pe ! number of processes per subgroup (Nw) 1041 iproc = para_env_sub%mepos ! subgroup-local process ID 1042 1043 ! Total number of processes and global process ID 1044 CALL mp_environ(nproc_glob, iproc_glob, mp_comm) 1045 1046 ! local block index for R/P and a 1047 loc_P = igroup + 1; loc_a = iproc + 1 1048 1049 CPASSERT(SIZE(sizes_L) .EQ. ngroup) 1050 CPASSERT(SIZE(sizes_B) .EQ. nproc) 1051 CPASSERT(sizes_L(loc_P) .EQ. SIZE(BIb_C, 1)) 1052 CPASSERT(sizes_L(loc_P) .EQ. SIZE(my_Lrows, 2)) 1053 CPASSERT(sizes_B(loc_a) .EQ. SIZE(BIb_C, 2)) 1054 1055 ! Tensor distributions as follows: 1056 ! Process grid NG x Nw 1057 ! Each process has coordinates (np, nw) 1058 ! tB_in: (R|ai): R distributed (np), a distributed (nw) 1059 ! tB_out: (P|ai): P distributed (np), a distributed (nw) 1060 ! tL: (R|P): R distributed (nw), P distributed (np) 1061 1062 ! define mappings between tensor index and matrix index: 1063 ! (R|ai) and (P|ai): 1064 map_B_1 = [1] ! index 1 (R or P) maps to 1st matrix index (np distributed) 1065 map_B_2 = [2, 3] ! indices 2, 3 (a, i) map to 2nd matrix index (nw distributed) 1066 ! (R|P): 1067 map_L_1 = [2] ! index 2 (P) maps to 1st matrix index (np distributed) 1068 map_L_2 = [1] ! index 1 (R) maps to 2nd matrix index (nw distributed) 1069 1070 ! derive nd process grid that is compatible with distributions and 2d process grid 1071 ! (R|ai) / (P|ai) on process grid NG x Nw x 1 1072 ! (R|P) on process grid NG x Nw 1073 pdims_B = [ngroup, nproc, 1] 1074 pdims_L = [nproc, ngroup] 1075 1076 CALL dbcsr_t_pgrid_create(mp_comm, pdims_B, mp_comm_B) 1077 CALL dbcsr_t_pgrid_create(mp_comm, pdims_L, mp_comm_L) 1078 1079 ! setup distribution vectors such that distribution matches parallel data layout of BIb_C and my_Lrows 1080 dist_B_i = [0] 1081 dist_B_a = (/(i, i=0, nproc - 1)/) 1082 dist_L_R = (/(MODULO(i, nproc), i=0, ngroup - 1)/) ! R index is replicated in my_Lrows, we impose a cyclic distribution 1083 dist_L_P = (/(i, i=0, ngroup - 1)/) 1084 1085 ! create distributions and tensors 1086 CALL dbcsr_t_distribution_new(dist_B, mp_comm_B, dist_L_P, dist_B_a, dist_B_i) 1087 CALL dbcsr_t_distribution_new(dist_L, mp_comm_L, dist_L_R, dist_L_P) 1088 1089 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) 1090 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) 1091 CALL dbcsr_t_create(tL, "(R|P)", dist_L, map_L_1, map_L_2, dbcsr_type_real_8, sizes_L, sizes_L) 1092 1093 IF (debug) THEN 1094 ! check that tensor distribution is correct 1095 CALL dbcsr_t_get_stored_coordinates(tB_in, [loc_P, loc_a, 1], check_proc) 1096 CPASSERT(check_proc .EQ. iproc_glob) 1097 ENDIF 1098 1099 ! reserve (R|ai) block 1100 CALL dbcsr_t_reserve_blocks(tB_in, [loc_P], [loc_a], [1]) 1101 1102 ! reserve (R|P) blocks 1103 ! in my_Lrows, R index is replicated. For (R|P), we distribute quadratic blocks cyclically over 1104 ! the processes in a subgroup. 1105 ! There are NG blocks, so each process holds at most NG/Nw+1 blocks. 1106 ALLOCATE (block_ind_L_R(ngroup/nproc + 1)) 1107 ALLOCATE (block_ind_L_P(ngroup/nproc + 1)) 1108 block_ind_L_R(:) = 0; block_ind_L_P(:) = 0 1109 ii = 0 1110 DO i = 1, ngroup 1111 CALL dbcsr_t_get_stored_coordinates(tL, [i, loc_P], check_proc) 1112 IF (check_proc == iproc_glob) THEN 1113 ii = ii + 1 1114 block_ind_L_R(ii) = i 1115 block_ind_L_P(ii) = loc_P 1116 ENDIF 1117 ENDDO 1118 CALL dbcsr_t_reserve_blocks(tL, block_ind_L_R(1:ii), block_ind_L_P(1:ii)) 1119 1120 ! insert (R|ai) block 1121 CALL dbcsr_t_put_block(tB_in, [loc_P, loc_a, 1], SHAPE(BIb_C), BIb_C) 1122 1123 ! insert (R|P) blocks 1124 ioff = 0 1125 DO i = 1, ngroup 1126 istart = ioff + 1; iend = ioff + sizes_L(i) 1127 ioff = ioff + sizes_L(i) 1128 CALL dbcsr_t_get_stored_coordinates(tL, [i, loc_P], check_proc) 1129 IF (check_proc == iproc_glob) THEN 1130 CALL dbcsr_t_put_block(tL, [i, loc_P], [sizes_L(i), sizes_L(loc_P)], my_Lrows(istart:iend, :)) 1131 ENDIF 1132 ENDDO 1133 1134 CALL dbcsr_t_split_blocks(tB_in, tB_in_split, [blk_size(2), blk_size(1), blk_size(1)]) 1135 CALL dbcsr_t_split_blocks(tL, tL_split, [blk_size(2), blk_size(2)]) 1136 CALL dbcsr_t_split_blocks(tB_out, tB_out_split, [blk_size(2), blk_size(1), blk_size(1)]) 1137 1138 ! contract 1139 CALL dbcsr_t_contract(alpha=dbcsr_scalar(1.0_dp), tensor_1=tB_in_split, tensor_2=tL_split, & 1140 beta=dbcsr_scalar(0.0_dp), tensor_3=tB_out_split, & 1141 contract_1=[1], notcontract_1=[2, 3], & 1142 contract_2=[1], notcontract_2=[2], & 1143 map_1=[2, 3], map_2=[1], optimize_dist=.TRUE.) 1144 1145 ! retrieve local block of contraction result (P|ai) 1146 CALL dbcsr_t_copy(tB_out_split, tB_out) 1147 1148 CALL dbcsr_t_get_block(tB_out, [loc_P, loc_a, 1], SHAPE(BIb_C), BIb_C, found) 1149 CPASSERT(found) 1150 1151 ! cleanup 1152 CALL dbcsr_t_destroy(tB_in) 1153 CALL dbcsr_t_destroy(tB_in_split) 1154 CALL dbcsr_t_destroy(tB_out) 1155 CALL dbcsr_t_destroy(tB_out_split) 1156 CALL dbcsr_t_destroy(tL) 1157 CALL dbcsr_t_destroy(tL_split) 1158 1159 CALL dbcsr_t_distribution_destroy(dist_B) 1160 CALL dbcsr_t_distribution_destroy(dist_L) 1161 1162 CALL dbcsr_t_pgrid_destroy(mp_comm_B) 1163 CALL dbcsr_t_pgrid_destroy(mp_comm_L) 1164 1165 CALL timestop(handle) 1166 1167 END SUBROUTINE contract_B_L 1168 1169! ************************************************************************************************** 1170!> \brief Encapsulate building of intermediate matrices matrix_ia_jnu(_beta 1171!> matrix_ia_jb(_beta),fm_BIb_jb(_beta),matrix_in_jnu(for G0W0) and 1172!> fm_BIb_all(for G0W0) 1173!> \param matrix_ia_jnu ... 1174!> \param matrix_ia_jb ... 1175!> \param mo_coeff_templ ... 1176!> \param size_1 ... 1177!> \param size_2 ... 1178!> \param fm_BIb_jb ... 1179!> \param matrix_name_2 ... 1180!> \param max_row_col_local ... 1181!> \param blacs_env_sub ... 1182!> \param para_env_sub ... 1183!> \param local_col_row_info ... 1184!> \author Jan Wilhelm 1185! ************************************************************************************************** 1186 SUBROUTINE create_intermediate_matrices(matrix_ia_jnu, matrix_ia_jb, mo_coeff_templ, size_1, size_2, & 1187 fm_BIb_jb, matrix_name_2, max_row_col_local, & 1188 blacs_env_sub, para_env_sub, local_col_row_info) 1189 1190 TYPE(dbcsr_type), INTENT(OUT) :: matrix_ia_jnu, matrix_ia_jb 1191 TYPE(dbcsr_type), INTENT(INOUT) :: mo_coeff_templ 1192 INTEGER, INTENT(IN) :: size_1, size_2 1193 TYPE(cp_fm_type), POINTER :: fm_BIb_jb 1194 CHARACTER(LEN=*), INTENT(IN) :: matrix_name_2 1195 INTEGER, INTENT(OUT) :: max_row_col_local 1196 TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub 1197 TYPE(cp_para_env_type), POINTER :: para_env_sub 1198 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: local_col_row_info 1199 1200 CHARACTER(LEN=*), PARAMETER :: routineN = 'create_intermediate_matrices', & 1201 routineP = moduleN//':'//routineN 1202 1203 INTEGER :: handle, ncol_local, nfullcols_total, & 1204 nfullrows_total, nrow_local 1205 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 1206 TYPE(cp_fm_struct_type), POINTER :: fm_struct 1207 1208 CALL timeset(routineN, handle) 1209 1210 ! initialize and create the matrix (K|jnu) 1211 CALL dbcsr_create(matrix_ia_jnu, template=mo_coeff_templ) 1212 1213 ! Allocate Sparse matrices: (K|jb) 1214 CALL cp_dbcsr_m_by_n_from_template(matrix_ia_jb, template=mo_coeff_templ, m=size_2, n=size_1, & 1215 sym=dbcsr_type_no_symmetry) 1216 1217 ! set all to zero in such a way that the memory is actually allocated 1218 CALL dbcsr_set(matrix_ia_jnu, 0.0_dp) 1219 CALL dbcsr_set(matrix_ia_jb, 0.0_dp) 1220 1221 ! create the analogous of matrix_ia_jb in fm type 1222 NULLIFY (fm_BIb_jb) 1223 NULLIFY (fm_struct) 1224 CALL dbcsr_get_info(matrix_ia_jb, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total) 1225 CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, nrow_global=nfullrows_total, & 1226 ncol_global=nfullcols_total, para_env=para_env_sub) 1227 CALL cp_fm_create(fm_BIb_jb, fm_struct, name=matrix_name_2) 1228 1229 CALL copy_dbcsr_to_fm(matrix_ia_jb, fm_BIb_jb) 1230 CALL cp_fm_struct_release(fm_struct) 1231 1232 CALL cp_fm_get_info(matrix=fm_BIb_jb, & 1233 nrow_local=nrow_local, & 1234 ncol_local=ncol_local, & 1235 row_indices=row_indices, & 1236 col_indices=col_indices) 1237 1238 max_row_col_local = MAX(nrow_local, ncol_local) 1239 CALL mp_max(max_row_col_local, para_env_sub%group) 1240 1241 ALLOCATE (local_col_row_info(0:max_row_col_local, 2)) 1242 local_col_row_info = 0 1243 ! 0,1 nrows 1244 local_col_row_info(0, 1) = nrow_local 1245 local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local) 1246 ! 0,2 ncols 1247 local_col_row_info(0, 2) = ncol_local 1248 local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local) 1249 1250 CALL timestop(handle) 1251 1252 END SUBROUTINE create_intermediate_matrices 1253 1254! ************************************************************************************************** 1255!> \brief Encapsulate ERI postprocessing: AO to MO transformation and store in B matrix. 1256!> \param para_env ... 1257!> \param mat_munu ... 1258!> \param matrix_ia_jnu ... 1259!> \param matrix_ia_jb ... 1260!> \param fm_BIb_jb ... 1261!> \param BIb_jb ... 1262!> \param mo_coeff_o ... 1263!> \param mo_coeff_v ... 1264!> \param eps_filter ... 1265!> \param max_row_col_local ... 1266!> \param sub_proc_map ... 1267!> \param local_col_row_info ... 1268!> \param my_B_end ... 1269!> \param my_B_start ... 1270!> \param descr ... 1271! ************************************************************************************************** 1272 SUBROUTINE ao_to_mo_and_store_B(para_env, mat_munu, matrix_ia_jnu, matrix_ia_jb, fm_BIb_jb, BIb_jb, & 1273 mo_coeff_o, mo_coeff_v, eps_filter, & 1274 max_row_col_local, sub_proc_map, local_col_row_info, & 1275 my_B_end, my_B_start, descr) 1276 TYPE(cp_para_env_type), POINTER :: para_env 1277 TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu 1278 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ia_jnu, matrix_ia_jb 1279 TYPE(cp_fm_type), POINTER :: fm_BIb_jb 1280 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: BIb_jb 1281 TYPE(dbcsr_type), POINTER :: mo_coeff_o, mo_coeff_v 1282 REAL(KIND=dp), INTENT(IN) :: eps_filter 1283 INTEGER, INTENT(IN) :: max_row_col_local 1284 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: sub_proc_map 1285 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: local_col_row_info 1286 INTEGER, INTENT(IN) :: my_B_end, my_B_start 1287 CHARACTER(LEN=*), INTENT(IN) :: descr 1288 1289 CHARACTER(LEN=*), PARAMETER :: routineN = 'ao_to_mo_and_store_B' 1290 1291 INTEGER :: handle1, handle2 1292 1293 CALL timeset(routineN//"_mult_"//descr, handle1) 1294 1295 CALL dbcsr_multiply("N", "N", 1.0_dp, mat_munu%matrix, mo_coeff_o, & 1296 0.0_dp, matrix_ia_jnu, filter_eps=eps_filter) 1297 CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_ia_jnu, mo_coeff_v, & 1298 0.0_dp, matrix_ia_jb, filter_eps=eps_filter) 1299 CALL timestop(handle1) 1300 1301 CALL timeset(routineN//"_E_Ex_"//descr, handle2) 1302 CALL copy_dbcsr_to_fm(matrix_ia_jb, fm_BIb_jb) 1303 1304 IF (.NOT. (descr .EQ. "bse_ab")) THEN 1305 1306 CALL grep_my_integrals(para_env, fm_BIb_jb, BIb_jb, max_row_col_local, & 1307 sub_proc_map, local_col_row_info, & 1308 my_B_end, my_B_start) 1309 1310 END IF 1311 1312 CALL timestop(handle2) 1313 END SUBROUTINE ao_to_mo_and_store_B 1314 1315! ************************************************************************************************** 1316!> \brief ... 1317!> \param qs_env ... 1318!> \param kpoints ... 1319!> \param unit_nr ... 1320! ************************************************************************************************** 1321 SUBROUTINE compute_kpoints(qs_env, kpoints, unit_nr) 1322 1323 TYPE(qs_environment_type), POINTER :: qs_env 1324 TYPE(kpoint_type), POINTER :: kpoints 1325 INTEGER :: unit_nr 1326 1327 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_kpoints', & 1328 routineP = moduleN//':'//routineN 1329 1330 INTEGER :: handle, i, i_dim, ix, iy, iz, nkp, & 1331 num_dim 1332 INTEGER, DIMENSION(3) :: nkp_grid, periodic 1333 TYPE(cell_type), POINTER :: cell 1334 TYPE(cp_para_env_type), POINTER :: para_env 1335 TYPE(dft_control_type), POINTER :: dft_control 1336 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1337 POINTER :: sab_orb 1338 1339 CALL timeset(routineN, handle) 1340 1341 NULLIFY (cell, dft_control, para_env) 1342 CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env, dft_control=dft_control, sab_orb=sab_orb) 1343 CALL get_cell(cell=cell, periodic=periodic) 1344 1345 ! general because we augment a Monkhorst-Pack mesh by additional points in the BZ 1346 kpoints%kp_scheme = "GENERAL" 1347 kpoints%symmetry = .FALSE. 1348 kpoints%verbose = .FALSE. 1349 kpoints%full_grid = .FALSE. 1350 kpoints%use_real_wfn = .FALSE. 1351 kpoints%eps_geo = 1.e-6_dp 1352 kpoints%full_grid = .TRUE. 1353 nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid(1:3) 1354 kpoints%nkp_grid(1:3) = nkp_grid(1:3) 1355 1356 num_dim = periodic(1) + periodic(2) + periodic(3) 1357 1358 DO i_dim = 1, 3 1359 IF (periodic(i_dim) == 1) THEN 1360 CPASSERT(MODULO(kpoints%nkp_grid(i_dim), 2) == 0) 1361 END IF 1362 END DO 1363 1364 IF (num_dim == 3) THEN 1365 nkp = nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2 1366 ELSE IF (num_dim == 2) THEN 1367 nkp = MAX(nkp_grid(1)*nkp_grid(2)/2, nkp_grid(1)*nkp_grid(3)/2, nkp_grid(2)*nkp_grid(3)/2) 1368 ELSE IF (num_dim == 1) THEN 1369 nkp = MAX(nkp_grid(1)/2, nkp_grid(2)/2, nkp_grid(3)/2) 1370 END IF 1371 1372 kpoints%nkp = nkp 1373 1374 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 1375 "KPOINT_INFO| Number of kpoints for V and W:", nkp 1376 1377 ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp)) 1378 kpoints%wkp(:) = 1.0_dp/REAL(nkp, KIND=dp) 1379 1380 i = 0 1381 DO ix = 1, nkp_grid(1) 1382 DO iy = 1, nkp_grid(2) 1383 DO iz = 1, nkp_grid(3) 1384 1385 i = i + 1 1386 IF (i > nkp) CYCLE 1387 1388 kpoints%xkp(1, i) = REAL(2*ix - nkp_grid(1) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(1), KIND=dp)) 1389 kpoints%xkp(2, i) = REAL(2*iy - nkp_grid(2) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(2), KIND=dp)) 1390 kpoints%xkp(3, i) = REAL(2*iz - nkp_grid(3) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(3), KIND=dp)) 1391 1392 END DO 1393 END DO 1394 END DO 1395 1396 CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, dft_control) 1397 1398 CALL set_qs_env(qs_env, kpoints=kpoints) 1399 1400 CALL timestop(handle) 1401 1402 END SUBROUTINE compute_kpoints 1403 1404! ************************************************************************************************** 1405!> \brief ... 1406!> \param para_env ... 1407!> \param dimen_RI ... 1408!> \param fm_matrix_L ... 1409!> \param my_group_L_start ... 1410!> \param my_group_L_end ... 1411!> \param my_group_L_size ... 1412!> \param my_Lrows ... 1413! ************************************************************************************************** 1414 SUBROUTINE grep_Lcols(para_env, dimen_RI, fm_matrix_L, & 1415 my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows) 1416 TYPE(cp_para_env_type), POINTER :: para_env 1417 INTEGER, INTENT(IN) :: dimen_RI 1418 TYPE(cp_fm_type), POINTER :: fm_matrix_L 1419 INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, & 1420 my_group_L_size 1421 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1422 INTENT(OUT) :: my_Lrows 1423 1424 CHARACTER(LEN=*), PARAMETER :: routineN = 'grep_Lcols', routineP = moduleN//':'//routineN 1425 1426 INTEGER :: handle, handle2, i_global, iiB, j_global, jjB, max_row_col_local, ncol_local, & 1427 ncol_rec, nrow_local, nrow_rec, proc_receive, proc_receive_static, proc_send, & 1428 proc_send_static, proc_shift 1429 INTEGER, ALLOCATABLE, DIMENSION(:) :: proc_map 1430 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, rec_col_row_info 1431 INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_rec, & 1432 row_indices, row_indices_rec 1433 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_L, local_L_internal, rec_L 1434 1435 CALL timeset(routineN, handle) 1436 1437 ALLOCATE (my_Lrows(dimen_RI, my_group_L_size)) 1438 my_Lrows = 0.0_dp 1439 1440 ! proc_map, vector that replicates the processor numbers also 1441 ! for negative and positive number > num_pe 1442 ! needed to know which is the processor, to respect to another one, 1443 ! for a given shift 1444 ALLOCATE (proc_map(-para_env%num_pe:2*para_env%num_pe - 1)) 1445 DO iiB = 0, para_env%num_pe - 1 1446 proc_map(iiB) = iiB 1447 proc_map(-iiB - 1) = para_env%num_pe - iiB - 1 1448 proc_map(para_env%num_pe + iiB) = iiB 1449 END DO 1450 1451 CALL cp_fm_get_info(matrix=fm_matrix_L, & 1452 nrow_local=nrow_local, & 1453 ncol_local=ncol_local, & 1454 row_indices=row_indices, & 1455 col_indices=col_indices, & 1456 local_data=local_L_internal) 1457 1458 ALLOCATE (local_L(nrow_local, ncol_local)) 1459 local_L = local_L_internal(1:nrow_local, 1:ncol_local) 1460 1461 max_row_col_local = MAX(nrow_local, ncol_local) 1462 CALL mp_max(max_row_col_local, para_env%group) 1463 1464 ALLOCATE (local_col_row_info(0:max_row_col_local, 2)) 1465 local_col_row_info = 0 1466 ! 0,1 nrows 1467 local_col_row_info(0, 1) = nrow_local 1468 local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local) 1469 ! 0,2 ncols 1470 local_col_row_info(0, 2) = ncol_local 1471 local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local) 1472 1473 ALLOCATE (rec_col_row_info(0:max_row_col_local, 2)) 1474 1475 ! accumulate data on my_Lrows starting from myself 1476 DO jjB = 1, ncol_local 1477 j_global = col_indices(jjB) 1478 IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN 1479 DO iiB = 1, nrow_local 1480 i_global = row_indices(iiB) 1481 my_Lrows(i_global, j_global - my_group_L_start + 1) = local_L(iiB, jjB) 1482 END DO 1483 END IF 1484 END DO 1485 1486 proc_send_static = proc_map(para_env%mepos + 1) 1487 proc_receive_static = proc_map(para_env%mepos - 1) 1488 1489 CALL timeset(routineN//"_comm", handle2) 1490 1491 DO proc_shift = 1, para_env%num_pe - 1 1492 proc_send = proc_map(para_env%mepos + proc_shift) 1493 proc_receive = proc_map(para_env%mepos - proc_shift) 1494 1495 ! first exchange information on the local data 1496 rec_col_row_info = 0 1497 CALL mp_sendrecv(local_col_row_info, proc_send_static, rec_col_row_info, proc_receive_static, para_env%group) 1498 nrow_rec = rec_col_row_info(0, 1) 1499 ncol_rec = rec_col_row_info(0, 2) 1500 1501 ALLOCATE (row_indices_rec(nrow_rec)) 1502 row_indices_rec = rec_col_row_info(1:nrow_rec, 1) 1503 1504 ALLOCATE (col_indices_rec(ncol_rec)) 1505 col_indices_rec = rec_col_row_info(1:ncol_rec, 2) 1506 1507 ALLOCATE (rec_L(nrow_rec, ncol_rec)) 1508 rec_L = 0.0_dp 1509 1510 ! then send and receive the real data 1511 CALL mp_sendrecv(local_L, proc_send_static, rec_L, proc_receive_static, para_env%group) 1512 1513 ! accumulate the received data on my_Lrows 1514 DO jjB = 1, ncol_rec 1515 j_global = col_indices_rec(jjB) 1516 IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN 1517 DO iiB = 1, nrow_rec 1518 i_global = row_indices_rec(iiB) 1519 my_Lrows(i_global, j_global - my_group_L_start + 1) = rec_L(iiB, jjB) 1520 END DO 1521 END IF 1522 END DO 1523 1524 local_col_row_info(:, :) = rec_col_row_info 1525 DEALLOCATE (local_L) 1526 ALLOCATE (local_L(nrow_rec, ncol_rec)) 1527 local_L = rec_L 1528 1529 DEALLOCATE (col_indices_rec) 1530 DEALLOCATE (row_indices_rec) 1531 DEALLOCATE (rec_L) 1532 END DO 1533 CALL timestop(handle2) 1534 1535 DEALLOCATE (local_col_row_info) 1536 DEALLOCATE (rec_col_row_info) 1537 DEALLOCATE (proc_map) 1538 DEALLOCATE (local_L) 1539 1540 CALL timestop(handle) 1541 1542 END SUBROUTINE grep_Lcols 1543 1544! ************************************************************************************************** 1545!> \brief ... 1546!> \param para_env_sub ... 1547!> \param fm_BIb_jb ... 1548!> \param BIb_jb ... 1549!> \param max_row_col_local ... 1550!> \param proc_map ... 1551!> \param local_col_row_info ... 1552!> \param my_B_virtual_end ... 1553!> \param my_B_virtual_start ... 1554! ************************************************************************************************** 1555 SUBROUTINE grep_my_integrals(para_env_sub, fm_BIb_jb, BIb_jb, max_row_col_local, & 1556 proc_map, local_col_row_info, & 1557 my_B_virtual_end, my_B_virtual_start) 1558 TYPE(cp_para_env_type), POINTER :: para_env_sub 1559 TYPE(cp_fm_type), POINTER :: fm_BIb_jb 1560 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: BIb_jb 1561 INTEGER, INTENT(IN) :: max_row_col_local 1562 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: proc_map 1563 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: local_col_row_info 1564 INTEGER, INTENT(IN) :: my_B_virtual_end, my_B_virtual_start 1565 1566 CHARACTER(LEN=*), PARAMETER :: routineN = 'grep_my_integrals', & 1567 routineP = moduleN//':'//routineN 1568 1569 INTEGER :: i_global, iiB, j_global, jjB, ncol_rec, & 1570 nrow_rec, proc_receive, proc_send, & 1571 proc_shift 1572 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: rec_col_row_info 1573 INTEGER, DIMENSION(:), POINTER :: col_indices_rec, row_indices_rec 1574 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_BI, rec_BI 1575 1576 ALLOCATE (rec_col_row_info(0:max_row_col_local, 2)) 1577 1578 rec_col_row_info(:, :) = local_col_row_info 1579 1580 nrow_rec = rec_col_row_info(0, 1) 1581 ncol_rec = rec_col_row_info(0, 2) 1582 1583 ALLOCATE (row_indices_rec(nrow_rec)) 1584 row_indices_rec = rec_col_row_info(1:nrow_rec, 1) 1585 1586 ALLOCATE (col_indices_rec(ncol_rec)) 1587 col_indices_rec = rec_col_row_info(1:ncol_rec, 2) 1588 1589 ! accumulate data on BIb_jb buffer starting from myself 1590 DO jjB = 1, ncol_rec 1591 j_global = col_indices_rec(jjB) 1592 IF (j_global >= my_B_virtual_start .AND. j_global <= my_B_virtual_end) THEN 1593 DO iiB = 1, nrow_rec 1594 i_global = row_indices_rec(iiB) 1595 BIb_jb(j_global - my_B_virtual_start + 1, i_global) = fm_BIb_jb%local_data(iiB, jjB) 1596 END DO 1597 END IF 1598 END DO 1599 1600 DEALLOCATE (row_indices_rec) 1601 DEALLOCATE (col_indices_rec) 1602 1603 IF (para_env_sub%num_pe > 1) THEN 1604 ALLOCATE (local_BI(nrow_rec, ncol_rec)) 1605 local_BI(1:nrow_rec, 1:ncol_rec) = fm_BIb_jb%local_data(1:nrow_rec, 1:ncol_rec) 1606 1607 DO proc_shift = 1, para_env_sub%num_pe - 1 1608 proc_send = proc_map(para_env_sub%mepos + proc_shift) 1609 proc_receive = proc_map(para_env_sub%mepos - proc_shift) 1610 1611 ! first exchange information on the local data 1612 rec_col_row_info = 0 1613 CALL mp_sendrecv(local_col_row_info, proc_send, rec_col_row_info, proc_receive, para_env_sub%group) 1614 nrow_rec = rec_col_row_info(0, 1) 1615 ncol_rec = rec_col_row_info(0, 2) 1616 1617 ALLOCATE (row_indices_rec(nrow_rec)) 1618 row_indices_rec = rec_col_row_info(1:nrow_rec, 1) 1619 1620 ALLOCATE (col_indices_rec(ncol_rec)) 1621 col_indices_rec = rec_col_row_info(1:ncol_rec, 2) 1622 1623 ALLOCATE (rec_BI(nrow_rec, ncol_rec)) 1624 rec_BI = 0.0_dp 1625 1626 ! then send and receive the real data 1627 CALL mp_sendrecv(local_BI, proc_send, rec_BI, proc_receive, para_env_sub%group) 1628 1629 ! accumulate the received data on BIb_jb buffer 1630 DO jjB = 1, ncol_rec 1631 j_global = col_indices_rec(jjB) 1632 IF (j_global >= my_B_virtual_start .AND. j_global <= my_B_virtual_end) THEN 1633 DO iiB = 1, nrow_rec 1634 i_global = row_indices_rec(iiB) 1635 BIb_jb(j_global - my_B_virtual_start + 1, i_global) = rec_BI(iiB, jjB) 1636 END DO 1637 END IF 1638 END DO 1639 1640 DEALLOCATE (col_indices_rec) 1641 DEALLOCATE (row_indices_rec) 1642 DEALLOCATE (rec_BI) 1643 END DO 1644 1645 DEALLOCATE (local_BI) 1646 END IF 1647 1648 DEALLOCATE (rec_col_row_info) 1649 1650 END SUBROUTINE grep_my_integrals 1651 1652END MODULE mp2_integrals 1653