1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utility functions for RPA calculations 8!> \par History 9!> 06.2019 Moved from rpa_ri_gpw.F [Frederick Stein] 10! ************************************************************************************************** 11MODULE rpa_util 12 USE cell_types, ONLY: cell_type,& 13 get_cell 14 USE cp_blacs_env, ONLY: cp_blacs_env_create,& 15 cp_blacs_env_release,& 16 cp_blacs_env_type 17 USE cp_cfm_types, ONLY: cp_cfm_create,& 18 cp_cfm_release,& 19 cp_cfm_set_all,& 20 cp_cfm_type 21 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 22 copy_fm_to_dbcsr,& 23 dbcsr_allocate_matrix_set,& 24 dbcsr_deallocate_matrix_set 25 USE cp_fm_basic_linalg, ONLY: cp_fm_syrk,& 26 cp_fm_transpose 27 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose 28 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 29 cp_fm_struct_release,& 30 cp_fm_struct_type 31 USE cp_fm_types, ONLY: cp_fm_copy_general,& 32 cp_fm_create,& 33 cp_fm_get_info,& 34 cp_fm_p_type,& 35 cp_fm_release,& 36 cp_fm_set_all,& 37 cp_fm_to_fm,& 38 cp_fm_type 39 USE cp_gemm_interface, ONLY: cp_gemm 40 USE cp_para_env, ONLY: cp_para_env_create,& 41 cp_para_env_release 42 USE cp_para_types, ONLY: cp_para_env_type 43 USE dbcsr_api, ONLY: & 44 dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_filter, dbcsr_get_info, & 45 dbcsr_get_occupation, dbcsr_init_p, 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_release, dbcsr_release_p, dbcsr_set, dbcsr_type 48 USE dbcsr_tensor_api, ONLY: dbcsr_t_destroy,& 49 dbcsr_t_type 50 USE input_constants, ONLY: wfc_mm_style_gemm,& 51 wfc_mm_style_syrk 52 USE kinds, ONLY: dp 53 USE kpoint_types, ONLY: get_kpoint_info,& 54 kpoint_type 55 USE machine, ONLY: m_walltime 56 USE mathconstants, ONLY: z_zero 57 USE message_passing, ONLY: mp_comm_split_direct,& 58 mp_sum 59 USE mp2_laplace, ONLY: calc_fm_mat_S_laplace 60 USE mp2_types, ONLY: integ_mat_buffer_type,& 61 mp2_type,& 62 two_dim_int_array 63 USE qs_environment_types, ONLY: get_qs_env,& 64 qs_environment_type 65 USE rpa_communication, ONLY: fm_redistribute 66 USE rpa_gw_kpoints, ONLY: compute_wkp_W 67 USE rpa_im_time, ONLY: setup_mat_for_mem_cut_3c 68#include "./base/base_uses.f90" 69 70 IMPLICIT NONE 71 72 PRIVATE 73 74 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_util' 75 76 PUBLIC :: RPA_postprocessing_nokp, alloc_im_time, calc_mat_Q, RPA_postprocessing_start, get_mat_3c_overl_int_cut, & 77 dealloc_im_time, contract_P_omega_with_mat_L 78 79CONTAINS 80 81! ************************************************************************************************** 82!> \brief ... 83!> \param qs_env ... 84!> \param mp2_env ... 85!> \param para_env ... 86!> \param para_env_sub ... 87!> \param dimen_RI ... 88!> \param dimen_RI_red ... 89!> \param num_integ_points ... 90!> \param fm_mat_Q ... 91!> \param do_mao ... 92!> \param fm_mo_coeff_occ ... 93!> \param fm_mo_coeff_virt ... 94!> \param fm_matrix_L_RI_metric ... 95!> \param mat_munu ... 96!> \param mat_dm_occ_local ... 97!> \param mat_dm_virt_local ... 98!> \param mat_P_global ... 99!> \param mat_M ... 100!> \param mat_3c_overl_int ... 101!> \param mat_3c_overl_int_mao_for_occ ... 102!> \param mat_3c_overl_int_mao_for_virt ... 103!> \param do_dbcsr_t ... 104!> \param t_3c_O ... 105!> \param matrix_s ... 106!> \param kpoints ... 107!> \param eps_filter ... 108!> \param eps_filter_im_time ... 109!> \param do_ri_sos_laplace_mp2 ... 110!> \param cut_RI ... 111!> \param cut_memory ... 112!> \param nkp ... 113!> \param num_cells_dm ... 114!> \param num_3c_repl ... 115!> \param size_P ... 116!> \param n_group_col ... 117!> \param n_group_row ... 118!> \param ikp_local ... 119!> \param mepos_P_from_RI_row ... 120!> \param my_group_L_sizes_im_time ... 121!> \param my_group_L_starts_im_time ... 122!> \param row_from_LLL ... 123!> \param ends_array_prim_col ... 124!> \param ends_array_prim_fullcol ... 125!> \param ends_array_prim_fullrow ... 126!> \param ends_array_prim_row ... 127!> \param index_to_cell_3c ... 128!> \param sizes_array_prim_col ... 129!> \param sizes_array_prim_row ... 130!> \param starts_array_prim_col ... 131!> \param cell_to_index_3c ... 132!> \param non_zero_blocks_3c ... 133!> \param starts_array_prim_fullrow ... 134!> \param starts_array_prim_row ... 135!> \param starts_array_prim_fullcol ... 136!> \param col_blk_size ... 137!> \param ends_array_cm ... 138!> \param ends_array_cm_mao_occ ... 139!> \param ends_array_cm_mao_virt ... 140!> \param row_blk_offset ... 141!> \param row_blk_size ... 142!> \param starts_array_cm ... 143!> \param starts_array_cm_mao_occ ... 144!> \param starts_array_cm_mao_virt ... 145!> \param do_ic_model ... 146!> \param do_kpoints_cubic_RPA ... 147!> \param do_kpoints_from_Gamma ... 148!> \param do_ri_Sigma_x ... 149!> \param my_open_shell ... 150!> \param cycle_due_to_sparse_dm ... 151!> \param multiply_needed_occ ... 152!> \param multiply_needed_virt ... 153!> \param has_mat_P_blocks ... 154!> \param buffer_mat_M ... 155!> \param wkp_W ... 156!> \param cfm_mat_Q ... 157!> \param fm_mat_L ... 158!> \param fm_mat_RI_global_work ... 159!> \param fm_mat_work ... 160!> \param fm_mo_coeff_occ_scaled ... 161!> \param fm_mo_coeff_virt_scaled ... 162!> \param mat_dm ... 163!> \param mat_L ... 164!> \param mat_M_P_munu_occ ... 165!> \param mat_M_P_munu_virt ... 166!> \param mat_P_global_copy ... 167!> \param mat_SinvVSinv ... 168!> \param mat_M_mu_Pnu_occ ... 169!> \param mat_M_mu_Pnu_virt ... 170!> \param mat_P_omega ... 171!> \param mat_P_omega_kp ... 172!> \param mat_dm_loc_occ_cut ... 173!> \param mat_dm_loc_virt_cut ... 174!> \param mat_dm_loc_occ ... 175!> \param mat_dm_loc_virt ... 176!> \param mat_work ... 177!> \param offset_combi_block ... 178!> \param mat_P_omega_beta ... 179! ************************************************************************************************** 180 SUBROUTINE alloc_im_time(qs_env, mp2_env, para_env, para_env_sub, dimen_RI, dimen_RI_red, num_integ_points, & 181 fm_mat_Q, do_mao, fm_mo_coeff_occ, fm_mo_coeff_virt, & 182 fm_matrix_L_RI_metric, mat_munu, mat_dm_occ_local, mat_dm_virt_local, mat_P_global, mat_M, & 183 mat_3c_overl_int, mat_3c_overl_int_mao_for_occ, mat_3c_overl_int_mao_for_virt, do_dbcsr_t, & 184 t_3c_O, matrix_s, kpoints, eps_filter, eps_filter_im_time, & 185 do_ri_sos_laplace_mp2, cut_RI, cut_memory, nkp, num_cells_dm, num_3c_repl, & 186 size_P, n_group_col, n_group_row, ikp_local, mepos_P_from_RI_row, & 187 my_group_L_sizes_im_time, my_group_L_starts_im_time, row_from_LLL, ends_array_prim_col, & 188 ends_array_prim_fullcol, ends_array_prim_fullrow, ends_array_prim_row, index_to_cell_3c, & 189 sizes_array_prim_col, sizes_array_prim_row, starts_array_prim_col, cell_to_index_3c, & 190 non_zero_blocks_3c, starts_array_prim_fullrow, starts_array_prim_row, & 191 starts_array_prim_fullcol, col_blk_size, ends_array_cm, ends_array_cm_mao_occ, ends_array_cm_mao_virt, & 192 row_blk_offset, row_blk_size, starts_array_cm, starts_array_cm_mao_occ, & 193 starts_array_cm_mao_virt, do_ic_model, do_kpoints_cubic_RPA, & 194 do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, cycle_due_to_sparse_dm, multiply_needed_occ, & 195 multiply_needed_virt, has_mat_P_blocks, buffer_mat_M, wkp_W, & 196 cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, & 197 fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, mat_P_global_copy, & 198 mat_SinvVSinv, mat_M_mu_Pnu_occ, mat_M_mu_Pnu_virt, mat_P_omega, mat_P_omega_kp, mat_dm_loc_occ_cut, & 199 mat_dm_loc_virt_cut, mat_dm_loc_occ, mat_dm_loc_virt, mat_work, offset_combi_block, & 200 mat_P_omega_beta) 201 202 TYPE(qs_environment_type), POINTER :: qs_env 203 TYPE(mp2_type), POINTER :: mp2_env 204 TYPE(cp_para_env_type), POINTER :: para_env, para_env_sub 205 INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, num_integ_points 206 TYPE(cp_fm_type), POINTER :: fm_mat_Q 207 LOGICAL, INTENT(IN) :: do_mao 208 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_occ, fm_mo_coeff_virt 209 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_matrix_L_RI_metric 210 TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu, mat_dm_occ_local, & 211 mat_dm_virt_local, mat_P_global, mat_M 212 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int, & 213 mat_3c_overl_int_mao_for_occ, & 214 mat_3c_overl_int_mao_for_virt 215 LOGICAL, INTENT(IN) :: do_dbcsr_t 216 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :), & 217 INTENT(INOUT) :: t_3c_O 218 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 219 TYPE(kpoint_type), POINTER :: kpoints 220 REAL(KIND=dp), INTENT(IN) :: eps_filter, eps_filter_im_time 221 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2 222 INTEGER, INTENT(IN) :: cut_RI, cut_memory 223 INTEGER, INTENT(OUT) :: nkp, num_cells_dm, num_3c_repl, size_P, & 224 n_group_col, n_group_row 225 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ikp_local, mepos_P_from_RI_row, & 226 my_group_L_sizes_im_time, & 227 my_group_L_starts_im_time, row_from_LLL 228 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: ends_array_prim_col, & 229 ends_array_prim_fullcol, ends_array_prim_fullrow, ends_array_prim_row, index_to_cell_3c, & 230 sizes_array_prim_col, sizes_array_prim_row, starts_array_prim_col 231 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & 232 INTENT(OUT) :: cell_to_index_3c, non_zero_blocks_3c 233 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: starts_array_prim_fullrow, & 234 starts_array_prim_row, & 235 starts_array_prim_fullcol 236 INTEGER, DIMENSION(:), POINTER :: col_blk_size, ends_array_cm, ends_array_cm_mao_occ, & 237 ends_array_cm_mao_virt, row_blk_offset, row_blk_size, starts_array_cm, & 238 starts_array_cm_mao_occ, starts_array_cm_mao_virt 239 LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA, & 240 do_kpoints_from_Gamma, do_ri_Sigma_x, & 241 my_open_shell 242 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :), & 243 INTENT(OUT) :: cycle_due_to_sparse_dm, & 244 multiply_needed_occ, & 245 multiply_needed_virt 246 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), & 247 INTENT(OUT) :: has_mat_P_blocks 248 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 249 INTENT(OUT) :: buffer_mat_M 250 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 251 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 252 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L 253 TYPE(cp_fm_type), POINTER :: fm_mat_RI_global_work, fm_mat_work, & 254 fm_mo_coeff_occ_scaled, & 255 fm_mo_coeff_virt_scaled 256 TYPE(dbcsr_p_type), INTENT(OUT) :: mat_dm, mat_L, mat_M_P_munu_occ, & 257 mat_M_P_munu_virt, mat_P_global_copy, & 258 mat_SinvVSinv 259 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_M_mu_Pnu_occ, mat_M_mu_Pnu_virt 260 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega, mat_P_omega_kp 261 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_dm_loc_occ_cut, mat_dm_loc_virt_cut 262 TYPE(dbcsr_type), POINTER :: mat_dm_loc_occ, mat_dm_loc_virt, mat_work 263 TYPE(two_dim_int_array), ALLOCATABLE, & 264 DIMENSION(:, :), INTENT(OUT) :: offset_combi_block 265 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_beta 266 267 CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_im_time', routineP = moduleN//':'//routineN 268 269 INTEGER :: cell_grid_dm(3), col_start_local, color_sub_col, color_sub_row, first_ikp_local, & 270 handle, i_cell, i_cut_RI, i_dim, i_mem, j_cell, j_mem, LLL, n_local_col, n_local_row, & 271 nblkrows_total, periodic(3), row, row_start_local 272 TYPE(cell_type), POINTER :: cell 273 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_sub_kp 274 275 CALL timeset(routineN, handle) 276 277 ALLOCATE (my_group_L_starts_im_time(cut_RI)) 278 my_group_L_starts_im_time(:) = mp2_env%ri_rpa_im_time_util(1)%my_group_L_starts_im_time 279 ALLOCATE (my_group_L_sizes_im_time(cut_RI)) 280 my_group_L_sizes_im_time(:) = mp2_env%ri_rpa_im_time_util(1)%my_group_L_sizes_im_time 281 282 IF (.NOT. do_dbcsr_t) THEN 283 num_3c_repl = SIZE(mat_3c_overl_int, 3) 284 ELSE 285 num_3c_repl = SIZE(t_3c_O, 2) 286 ENDIF 287 288 IF (.NOT. do_dbcsr_t) THEN 289 290 DO i_cut_RI = 1, cut_RI 291 292 DO i_cell = 1, num_3c_repl 293 DO j_cell = 1, num_3c_repl 294 295 CALL dbcsr_filter(mat_3c_overl_int(i_cut_RI, i_cell, j_cell)%matrix, & 296 eps_filter) 297 298 END DO 299 END DO 300 301 END DO 302 303 CALL get_non_zero_blocks_3c(mat_3c_overl_int, para_env_sub, cut_RI, non_zero_blocks_3c) 304 305 ENDIF 306 307 IF (do_kpoints_cubic_RPA) THEN 308 ! we always use an odd number of image cells 309 ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical 310 DO i_dim = 1, 3 311 cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1 312 END DO 313 num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3) 314 ALLOCATE (index_to_cell_3c(3, SIZE(kpoints%index_to_cell, 2))) 315 CPASSERT(SIZE(kpoints%index_to_cell, 1) == 3) 316 index_to_cell_3c(:, :) = kpoints%index_to_cell(:, :) 317 ALLOCATE (cell_to_index_3c(LBOUND(kpoints%cell_to_index, 1):UBOUND(kpoints%cell_to_index, 1), & 318 LBOUND(kpoints%cell_to_index, 2):UBOUND(kpoints%cell_to_index, 2), & 319 LBOUND(kpoints%cell_to_index, 3):UBOUND(kpoints%cell_to_index, 3))) 320 cell_to_index_3c(:, :, :) = kpoints%cell_to_index(:, :, :) 321 322 ELSE 323 ALLOCATE (index_to_cell_3c(3, 1)) 324 index_to_cell_3c(:, 1) = 0 325 ALLOCATE (cell_to_index_3c(0:0, 0:0, 0:0)) 326 cell_to_index_3c(0, 0, 0) = 1 327 num_cells_dm = 1 328 END IF 329 330 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 331 332 CALL get_sub_para_kp(fm_struct_sub_kp, para_env, kpoints%nkp, & 333 dimen_RI, ikp_local, first_ikp_local, do_kpoints_cubic_RPA) 334 335 NULLIFY (cfm_mat_Q) 336 CALL cp_cfm_create(cfm_mat_Q, fm_struct_sub_kp) 337 CALL cp_cfm_set_all(cfm_mat_Q, z_zero) 338 ELSE 339 first_ikp_local = 1 340 END IF 341 342 IF (.NOT. do_dbcsr_t) THEN 343 CALL allocate_dm_loc(mat_dm_loc_occ, mat_dm_occ_local, mat_dm_loc_occ_cut, cut_RI, cut_memory, num_cells_dm) 344 345 CALL allocate_dm_loc(mat_dm_loc_virt, mat_dm_virt_local, mat_dm_loc_virt_cut, cut_RI, cut_memory, num_cells_dm) 346 347 CALL dbcsr_set(mat_munu%matrix, 0.0_dp) 348 CALL dbcsr_filter(mat_munu%matrix, 1.0_dp) 349 350 IF (.NOT. do_mao) THEN 351 mat_3c_overl_int_mao_for_occ => mat_3c_overl_int 352 mat_3c_overl_int_mao_for_virt => mat_3c_overl_int 353 END IF 354 355 CALL allocate_mat_M_P_munu(mat_M_P_munu_occ, mat_M, mat_M_mu_Pnu_occ, cut_RI, mat_3c_overl_int_mao_for_occ) 356 357 CALL allocate_mat_M_P_munu(mat_M_P_munu_virt, mat_M, mat_M_mu_Pnu_virt, cut_RI, mat_3c_overl_int_mao_for_virt) 358 359 ENDIF 360 361 ! if we do kpoints, mat_P has a kpoint and mat_P_omega has the inted 362 ! mat_P(tau, kpoint) 363 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 364 365 NULLIFY (cell) 366 CALL get_qs_env(qs_env, cell=cell) 367 CALL get_cell(cell=cell, periodic=periodic) 368 369 CALL get_kpoint_info(kpoints, nkp=nkp) 370 ! compute k-point weights such that functions 1/k^2, 1/k and const function are 371 ! integrated correctly 372 CALL compute_wkp_W(wkp_W, kpoints, cell%hmat, cell%h_inv, & 373 qs_env%mp2_env%ri_rpa_im_time%exp_kpoints, periodic) 374 ELSE 375 nkp = 1 376 END IF 377 378 IF (do_kpoints_cubic_RPA) THEN 379 size_P = MAX(num_cells_dm/2 + 1, nkp) 380 ELSE IF (do_kpoints_from_Gamma) THEN 381 size_P = MAX(3**(periodic(1) + periodic(2) + periodic(3)), nkp) 382 ELSE 383 size_P = 1 384 END IF 385 386 CALL alloc_mat_P_omega(mat_P_omega, num_integ_points, size_P, mat_P_global%matrix) 387 388 IF (my_open_shell .AND. do_ri_sos_laplace_mp2) THEN 389 CALL alloc_mat_P_omega(mat_P_omega_beta, num_integ_points, size_P, mat_P_global%matrix) 390 END IF 391 392 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 393 CALL alloc_mat_P_omega(mat_P_omega_kp, 2, size_P, mat_P_global%matrix) 394 END IF 395 396 IF (.NOT. do_dbcsr_t) THEN 397 NULLIFY (mat_P_global_copy%matrix) 398 ALLOCATE (mat_P_global_copy%matrix) 399 CALL dbcsr_create(mat_P_global_copy%matrix, template=mat_P_global%matrix) 400 CALL dbcsr_copy(mat_P_global_copy%matrix, mat_P_global%matrix) 401 402 n_group_row = mp2_env%ri_rpa_im_time_util(1)%n_group_row 403 ALLOCATE (sizes_array_prim_row(0:n_group_row - 1, cut_memory)) 404 DO i_mem = 1, cut_memory 405 sizes_array_prim_row(:, i_mem) = mp2_env%ri_rpa_im_time_util(i_mem)%sizes_array_prim_row(:) 406 END DO 407 ALLOCATE (starts_array_prim_row(0:n_group_row - 1, cut_memory)) 408 DO i_mem = 1, cut_memory 409 starts_array_prim_row(:, i_mem) = mp2_env%ri_rpa_im_time_util(i_mem)%starts_array_prim_row(:) 410 END DO 411 ALLOCATE (ends_array_prim_row(0:n_group_row - 1, cut_memory)) 412 DO i_mem = 1, cut_memory 413 ends_array_prim_row(:, i_mem) = mp2_env%ri_rpa_im_time_util(i_mem)%ends_array_prim_row(:) 414 END DO 415 ALLOCATE (starts_array_prim_fullrow(0:n_group_row - 1, cut_memory)) 416 DO i_mem = 1, cut_memory 417 starts_array_prim_fullrow(:, i_mem) = mp2_env%ri_rpa_im_time_util(i_mem)%starts_array_prim_fullrow(:) 418 END DO 419 ALLOCATE (ends_array_prim_fullrow(0:n_group_row - 1, cut_memory)) 420 DO i_mem = 1, cut_memory 421 ends_array_prim_fullrow(:, i_mem) = mp2_env%ri_rpa_im_time_util(i_mem)%ends_array_prim_fullrow(:) 422 END DO 423 424 n_group_col = mp2_env%ri_rpa_im_time_util(1)%n_group_col 425 ALLOCATE (sizes_array_prim_col(0:n_group_col - 1, cut_memory)) 426 DO j_mem = 1, cut_memory 427 sizes_array_prim_col(:, j_mem) = mp2_env%ri_rpa_im_time_util(j_mem)%sizes_array_prim_col(:) 428 END DO 429 ALLOCATE (starts_array_prim_col(0:n_group_col - 1, cut_memory)) 430 DO j_mem = 1, cut_memory 431 starts_array_prim_col(:, j_mem) = mp2_env%ri_rpa_im_time_util(j_mem)%starts_array_prim_col(:) 432 END DO 433 ALLOCATE (ends_array_prim_col(0:n_group_col - 1, cut_memory)) 434 DO j_mem = 1, cut_memory 435 ends_array_prim_col(:, j_mem) = mp2_env%ri_rpa_im_time_util(j_mem)%ends_array_prim_col(:) 436 END DO 437 438 ALLOCATE (starts_array_prim_fullcol(0:n_group_col - 1, cut_memory)) 439 DO j_mem = 1, cut_memory 440 starts_array_prim_fullcol(:, j_mem) = mp2_env%ri_rpa_im_time_util(j_mem)%starts_array_prim_fullcol(:) 441 END DO 442 ALLOCATE (ends_array_prim_fullcol(0:n_group_col - 1, cut_memory)) 443 DO j_mem = 1, cut_memory 444 ends_array_prim_fullcol(:, j_mem) = mp2_env%ri_rpa_im_time_util(j_mem)%ends_array_prim_fullcol(:) 445 END DO 446 447 ALLOCATE (offset_combi_block(cut_memory, cut_memory)) 448 449 color_sub_row = mp2_env%ri_rpa_im_time_util(1)%color_sub_row 450 color_sub_col = mp2_env%ri_rpa_im_time_util(1)%color_sub_col 451 452 DO i_mem = 1, cut_memory 453 DO j_mem = 1, cut_memory 454 455 n_local_row = sizes_array_prim_row(color_sub_row, i_mem) 456 row_start_local = starts_array_prim_row(color_sub_row, i_mem) 457 458 n_local_col = sizes_array_prim_col(color_sub_col, j_mem) 459 col_start_local = starts_array_prim_col(color_sub_col, j_mem) 460 461 ALLOCATE (offset_combi_block(i_mem, j_mem)%array(row_start_local:row_start_local + n_local_row - 1, & 462 col_start_local:col_start_local + n_local_col - 1)) 463 offset_combi_block(i_mem, j_mem)%array(:, :) = & 464 mp2_env%ri_rpa_im_time_2d_util(i_mem, j_mem)%offset_combi_block(:, :) 465 466 END DO 467 END DO 468 469 NULLIFY (starts_array_cm, ends_array_cm) 470 starts_array_cm => mp2_env%ri_rpa_im_time%starts_array_cm 471 ends_array_cm => mp2_env%ri_rpa_im_time%ends_array_cm 472 473 NULLIFY (starts_array_cm_mao_occ, starts_array_cm_mao_virt, ends_array_cm_mao_occ, ends_array_cm_mao_virt) 474 IF (do_mao) THEN 475 starts_array_cm_mao_occ => mp2_env%ri_rpa_im_time%starts_array_cm_mao_occ 476 starts_array_cm_mao_virt => mp2_env%ri_rpa_im_time%starts_array_cm_mao_virt 477 ends_array_cm_mao_occ => mp2_env%ri_rpa_im_time%ends_array_cm_mao_occ 478 ends_array_cm_mao_virt => mp2_env%ri_rpa_im_time%ends_array_cm_mao_virt 479 ELSE 480 starts_array_cm_mao_occ => starts_array_cm 481 starts_array_cm_mao_virt => starts_array_cm 482 ends_array_cm_mao_occ => ends_array_cm 483 ends_array_cm_mao_virt => ends_array_cm 484 END IF 485 486 ! allocatable arrays for fast filling of dbcsr matrices 487 CALL dbcsr_get_info(mat_M_P_munu_occ%matrix, row_blk_size=row_blk_size, & 488 col_blk_size=col_blk_size, nblkrows_total=nblkrows_total) 489 ALLOCATE (buffer_mat_M(MAXVAL(row_blk_size), MAXVAL(col_blk_size))) 490 491 ALLOCATE (mepos_P_from_RI_row(nblkrows_total)) 492 mepos_P_from_RI_row(:) = mp2_env%ri_rpa_im_time_util(1)%mepos_P_from_RI_row(:) 493 494 ! allocatable arrays for fast filling of dbcsr matrices 495 CALL dbcsr_get_info(mat_P_global%matrix, row_blk_size=row_blk_size, col_blk_size=col_blk_size) 496 497 ALLOCATE (row_from_LLL(dimen_RI)) 498 row_from_LLL = 0 499 500 CALL dbcsr_get_info(mat_M_P_munu_occ%matrix, & 501 nblkrows_total=nblkrows_total, & 502 row_blk_offset=row_blk_offset, & 503 row_blk_size=row_blk_size) 504 505 DO LLL = 1, dimen_RI 506 DO row = 1, nblkrows_total 507 IF (row_blk_offset(row) <= LLL .AND. LLL < row_blk_offset(row) + row_blk_size(row)) THEN 508 row_from_LLL(LLL) = row 509 END IF 510 END DO 511 END DO 512 ENDIF 513 514 NULLIFY (fm_mo_coeff_occ_scaled) 515 CALL cp_fm_create(fm_mo_coeff_occ_scaled, fm_mo_coeff_occ%matrix_struct) 516 CALL cp_fm_to_fm(fm_mo_coeff_occ, fm_mo_coeff_occ_scaled) 517 CALL cp_fm_set_all(matrix=fm_mo_coeff_occ_scaled, alpha=0.0_dp) 518 519 NULLIFY (fm_mo_coeff_virt_scaled) 520 CALL cp_fm_create(fm_mo_coeff_virt_scaled, fm_mo_coeff_virt%matrix_struct) 521 CALL cp_fm_to_fm(fm_mo_coeff_virt, fm_mo_coeff_virt_scaled) 522 CALL cp_fm_set_all(matrix=fm_mo_coeff_virt_scaled, alpha=0.0_dp) 523 524 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 525 NULLIFY (fm_mat_RI_global_work) 526 CALL cp_fm_create(fm_mat_RI_global_work, fm_matrix_L_RI_metric(1, 1)%matrix%matrix_struct) 527 CALL cp_fm_to_fm(fm_matrix_L_RI_metric(1, 1)%matrix, fm_mat_RI_global_work) 528 CALL cp_fm_set_all(fm_mat_RI_global_work, 0.0_dp) 529 END IF 530 531 IF (.NOT. do_dbcsr_t) THEN 532 ALLOCATE (cycle_due_to_sparse_dm(cut_memory, cut_memory, num_integ_points)) 533 cycle_due_to_sparse_dm = .FALSE. 534 535 ALLOCATE (multiply_needed_occ(cut_memory, cut_RI, num_cells_dm)) 536 multiply_needed_occ = .TRUE. 537 538 ALLOCATE (multiply_needed_virt(cut_memory, cut_RI, num_cells_dm)) 539 multiply_needed_virt = .TRUE. 540 ENDIF 541 542 ALLOCATE (has_mat_P_blocks(num_cells_dm/2 + 1, cut_memory, cut_memory, num_3c_repl, num_3c_repl)) 543 has_mat_P_blocks = .TRUE. 544 545 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 546 CALL reorder_mat_L(fm_mat_L, fm_matrix_L_RI_metric, fm_mat_Q%matrix_struct, para_env, mat_L, & 547 mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp) 548 549 CALL cp_fm_struct_release(fm_struct_sub_kp) 550 ELSE 551 CALL reorder_mat_L(fm_mat_L, fm_matrix_L_RI_metric, fm_mat_Q%matrix_struct, para_env, mat_L, & 552 mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local) 553 END IF 554 555 ! Create Scalapack working matrix for the contraction with the metric 556 IF (dimen_RI == dimen_RI_red) THEN 557 NULLIFY (fm_mat_work) 558 CALL cp_fm_create(fm_mat_work, fm_mat_Q%matrix_struct) 559 CALL cp_fm_to_fm(fm_mat_Q, fm_mat_work) 560 CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp) 561 562 ELSE 563 NULLIFY (fm_struct) 564 CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_mat_Q%matrix_struct, & 565 ncol_global=dimen_RI_red, nrow_global=dimen_RI) 566 567 NULLIFY (fm_mat_work) 568 CALL cp_fm_create(fm_mat_work, fm_struct) 569 CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp) 570 571 CALL cp_fm_struct_release(fm_struct) 572 573 END IF 574 575 ! Then its DBCSR counter part 576 IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN 577 CALL dbcsr_get_info(mat_L%matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size) 578 579 ! Create mat_work having the shape of the transposed of mat_L (compare with contract_P_omega_with_mat_L) 580 NULLIFY (mat_work) 581 ALLOCATE (mat_work) 582 CALL dbcsr_create(mat_work, template=mat_L%matrix, row_blk_size=col_blk_size, col_blk_size=row_blk_size) 583 END IF 584 585 IF (do_ri_Sigma_x .OR. do_ic_model) THEN 586 587 NULLIFY (mat_SinvVSinv%matrix) 588 ALLOCATE (mat_SinvVSinv%matrix) 589 CALL dbcsr_create(mat_SinvVSinv%matrix, template=mat_P_global%matrix) 590 CALL dbcsr_set(mat_SinvVSinv%matrix, 0.0_dp) 591 592 ! for kpoints we compute SinvVSinv later with kpoints 593 IF (.NOT. do_kpoints_from_Gamma) THEN 594 595 ! get the Coulomb matrix for Sigma_x = G*V 596 CALL dbcsr_multiply("T", "N", 1.0_dp, mat_L%matrix, mat_L%matrix, & 597 0.0_dp, mat_SinvVSinv%matrix, filter_eps=eps_filter_im_time) 598 599 END IF 600 601 END IF 602 603 IF (do_ri_Sigma_x) THEN 604 605 NULLIFY (mat_dm%matrix) 606 ALLOCATE (mat_dm%matrix) 607 CALL dbcsr_create(mat_dm%matrix, template=matrix_s(1)%matrix) 608 609 END IF 610 611 CALL timestop(handle) 612 613 END SUBROUTINE alloc_im_time 614 615! ************************************************************************************************** 616!> \brief ... 617!> \param mat_dm_loc ... 618!> \param mat_dm_local ... 619!> \param mat_dm_loc_cut ... 620!> \param cut_RI ... 621!> \param cut_memory ... 622!> \param num_cells_dm ... 623! ************************************************************************************************** 624 SUBROUTINE allocate_dm_loc(mat_dm_loc, mat_dm_local, mat_dm_loc_cut, cut_RI, cut_memory, num_cells_dm) 625 TYPE(dbcsr_type), POINTER :: mat_dm_loc 626 TYPE(dbcsr_p_type), INTENT(IN) :: mat_dm_local 627 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_dm_loc_cut 628 INTEGER, INTENT(IN) :: cut_RI, cut_memory, num_cells_dm 629 630 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_dm_loc', & 631 routineP = moduleN//':'//routineN 632 633 INTEGER :: handle, i_cell, i_cut_RI, i_mem 634 635 CALL timeset(routineN, handle) 636 637 NULLIFY (mat_dm_loc) 638 CALL dbcsr_init_p(mat_dm_loc) 639 CALL dbcsr_desymmetrize(mat_dm_local%matrix, mat_dm_loc) 640 641 NULLIFY (mat_dm_loc_cut) 642 CALL dbcsr_allocate_matrix_set(mat_dm_loc_cut, cut_RI, cut_memory, num_cells_dm) 643 644 DO i_mem = 1, cut_memory 645 DO i_cut_RI = 1, cut_RI 646 DO i_cell = 1, num_cells_dm 647 648 ALLOCATE (mat_dm_loc_cut(i_cut_RI, i_mem, i_cell)%matrix) 649 CALL dbcsr_create(matrix=mat_dm_loc_cut(i_cut_RI, i_mem, i_cell)%matrix, & 650 template=mat_dm_loc) 651 652 END DO 653 END DO 654 END DO 655 656 CALL timestop(handle) 657 658 END SUBROUTINE allocate_dm_loc 659 660! ************************************************************************************************** 661!> \brief ... 662!> \param mat_M_P_munu ... 663!> \param mat_M ... 664!> \param mat_M_mu_Pnu ... 665!> \param cut_RI ... 666!> \param mat_3c_overl_int_mao ... 667! ************************************************************************************************** 668 SUBROUTINE allocate_mat_M_P_munu(mat_M_P_munu, mat_M, mat_M_mu_Pnu, cut_RI, mat_3c_overl_int_mao) 669 TYPE(dbcsr_p_type), INTENT(OUT) :: mat_M_P_munu 670 TYPE(dbcsr_p_type), INTENT(IN) :: mat_M 671 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_M_mu_Pnu 672 INTEGER, INTENT(IN) :: cut_RI 673 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_mao 674 675 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_mat_M_P_munu', & 676 routineP = moduleN//':'//routineN 677 678 INTEGER :: handle, i_cut_RI 679 680 CALL timeset(routineN, handle) 681 682 NULLIFY (mat_M_P_munu%matrix) 683 ALLOCATE (mat_M_P_munu%matrix) 684 CALL dbcsr_create(mat_M_P_munu%matrix, template=mat_M%matrix) 685 686 NULLIFY (mat_M_mu_Pnu) 687 CALL dbcsr_allocate_matrix_set(mat_M_mu_Pnu, cut_RI) 688 DO i_cut_RI = 1, cut_RI 689 ALLOCATE (mat_M_mu_Pnu(i_cut_RI)%matrix) 690 CALL dbcsr_create(matrix=mat_M_mu_Pnu(i_cut_RI)%matrix, & 691 template=mat_3c_overl_int_mao(i_cut_RI, 1, 1)%matrix) 692 END DO 693 694 CALL timestop(handle) 695 696 END SUBROUTINE allocate_mat_M_P_munu 697 698! ************************************************************************************************** 699!> \brief ... 700!> \param mat_P_omega ... 701!> \param num_integ_points ... 702!> \param size_P ... 703!> \param template ... 704! ************************************************************************************************** 705 SUBROUTINE alloc_mat_P_omega(mat_P_omega, num_integ_points, size_P, template) 706 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega 707 INTEGER, INTENT(IN) :: num_integ_points, size_P 708 TYPE(dbcsr_type), POINTER :: template 709 710 CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_mat_P_omega', & 711 routineP = moduleN//':'//routineN 712 713 INTEGER :: handle, i_kp, jquad 714 715 CALL timeset(routineN, handle) 716 717 NULLIFY (mat_P_omega) 718 CALL dbcsr_allocate_matrix_set(mat_P_omega, num_integ_points, size_P) 719 DO jquad = 1, num_integ_points 720 DO i_kp = 1, size_P 721 ALLOCATE (mat_P_omega(jquad, i_kp)%matrix) 722 CALL dbcsr_create(matrix=mat_P_omega(jquad, i_kp)%matrix, & 723 template=template) 724 CALL dbcsr_set(mat_P_omega(jquad, i_kp)%matrix, 0.0_dp) 725 END DO 726 END DO 727 728 CALL timestop(handle) 729 730 END SUBROUTINE alloc_mat_P_omega 731 732! ************************************************************************************************** 733!> \brief ... 734!> \param fm_mat_L ... 735!> \param fm_matrix_L_RI_metric ... 736!> \param fm_struct_template ... 737!> \param para_env ... 738!> \param mat_L ... 739!> \param mat_template ... 740!> \param dimen_RI ... 741!> \param dimen_RI_red ... 742!> \param first_ikp_local ... 743!> \param ikp_local ... 744!> \param fm_struct_sub_kp ... 745! ************************************************************************************************** 746 SUBROUTINE reorder_mat_L(fm_mat_L, fm_matrix_L_RI_metric, fm_struct_template, para_env, mat_L, mat_template, & 747 dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp) 748 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L, fm_matrix_L_RI_metric 749 TYPE(cp_fm_struct_type), POINTER :: fm_struct_template 750 TYPE(cp_para_env_type), POINTER :: para_env 751 TYPE(dbcsr_p_type), INTENT(OUT) :: mat_L 752 TYPE(dbcsr_type), INTENT(IN) :: mat_template 753 INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, first_ikp_local 754 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: ikp_local 755 TYPE(cp_fm_struct_type), OPTIONAL, POINTER :: fm_struct_sub_kp 756 757 CHARACTER(LEN=*), PARAMETER :: routineN = 'reorder_mat_L', routineP = moduleN//':'//routineN 758 759 INTEGER :: handle, i_size, j_size, nblk 760 INTEGER, DIMENSION(:), POINTER :: row_blk_size 761 LOGICAL :: do_kpoints 762 TYPE(cp_blacs_env_type), POINTER :: blacs_env 763 TYPE(cp_fm_struct_type), POINTER :: fm_struct 764 TYPE(cp_fm_type), POINTER :: fm_mat_L_transposed, fmdummy 765 766 CALL timeset(routineN, handle) 767 768 do_kpoints = .FALSE. 769 IF (PRESENT(ikp_local) .AND. PRESENT(fm_struct_sub_kp)) THEN 770 do_kpoints = .TRUE. 771 END IF 772 773 ! Get the fm_struct for fm_mat_L 774 NULLIFY (fm_struct) 775 IF (dimen_RI == dimen_RI_red) THEN 776 fm_struct => fm_struct_template 777 ELSE 778 ! The template is assumed to be square such that we need a new fm_struct if dimensions are not equal 779 CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI_red, ncol_global=dimen_RI, template_fmstruct=fm_struct_template) 780 END IF 781 782 ! Start to allocate the new full matrix 783 ALLOCATE (fm_mat_L(SIZE(fm_matrix_L_RI_metric, 1), SIZE(fm_matrix_L_RI_metric, 2))) 784 DO i_size = 1, SIZE(fm_matrix_L_RI_metric, 1) 785 DO j_size = 1, SIZE(fm_matrix_L_RI_metric, 2) 786 IF (do_kpoints) THEN 787 IF (ANY(ikp_local(:) == i_size)) THEN 788 CALL cp_fm_create(fm_mat_L(i_size, j_size)%matrix, fm_struct_sub_kp) 789 CALL cp_fm_set_all(fm_mat_L(i_size, j_size)%matrix, 0.0_dp) 790 END IF 791 ELSE 792 CALL cp_fm_create(fm_mat_L(i_size, j_size)%matrix, fm_struct) 793 CALL cp_fm_set_all(fm_mat_L(i_size, j_size)%matrix, 0.0_dp) 794 END IF 795 END DO 796 END DO 797 798 ! For the transposed matric we need a different fm_struct 799 IF (dimen_RI == dimen_RI_red) THEN 800 fm_struct => fm_mat_L(first_ikp_local, 1)%matrix%matrix_struct 801 ELSE 802 CALL cp_fm_struct_release(fm_struct) 803 804 ! Create a fm_struct with transposed sizes 805 NULLIFY (fm_struct) 806 CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI, ncol_global=dimen_RI_red, & 807 template_fmstruct=fm_mat_L(first_ikp_local, 1)%matrix%matrix_struct) !, force_block=.TRUE.) 808 END IF 809 810 ! Allocate buffer matrix 811 NULLIFY (fm_mat_L_transposed) 812 CALL cp_fm_create(fm_mat_L_transposed, fm_struct) 813 CALL cp_fm_set_all(matrix=fm_mat_L_transposed, alpha=0.0_dp) 814 815 IF (dimen_RI /= dimen_RI_red) CALL cp_fm_struct_release(fm_struct) 816 817 CALL cp_fm_get_info(fm_mat_L_transposed, context=blacs_env) 818 819 ! For k-points copy matrices of your group 820 ! Without kpoints, transpose matrix 821 ! without kpoints, the size of fm_mat_L is 1x1. with kpoints, the size is N_kpoints x 2 (2 for real/complex) 822 DO i_size = 1, SIZE(fm_matrix_L_RI_metric, 1) 823 DO j_size = 1, SIZE(fm_matrix_L_RI_metric, 2) 824 IF (do_kpoints) THEN 825 IF (ANY(ikp_local(:) == i_size)) THEN 826 CALL cp_fm_copy_general(fm_matrix_L_RI_metric(i_size, j_size)%matrix, fm_mat_L_transposed, para_env) 827 CALL cp_fm_to_fm(fm_mat_L_transposed, fm_mat_L(i_size, j_size)%matrix) 828 ELSE 829 NULLIFY (fmdummy) 830 CALL cp_fm_copy_general(fm_matrix_L_RI_metric(i_size, j_size)%matrix, fmdummy, para_env) 831 END IF 832 ELSE 833 CALL cp_fm_copy_general(fm_matrix_L_RI_metric(i_size, j_size)%matrix, fm_mat_L_transposed, blacs_env%para_env) 834 CALL cp_fm_transpose(fm_mat_L_transposed, fm_mat_L(i_size, j_size)%matrix) 835 END IF 836 END DO 837 END DO 838 839 ! Release old matrix 840 DO i_size = 1, SIZE(fm_matrix_L_RI_metric, 1) 841 DO j_size = 1, SIZE(fm_matrix_L_RI_metric, 2) 842 CALL cp_fm_release(fm_matrix_L_RI_metric(i_size, j_size)%matrix) 843 END DO 844 END DO 845 DEALLOCATE (fm_matrix_L_RI_metric) 846 ! Release buffer 847 CALL cp_fm_release(fm_mat_L_transposed) 848 849 ! Create sparse variant of L 850 NULLIFY (mat_L%matrix) 851 ALLOCATE (mat_L%matrix) 852 IF (dimen_RI == dimen_RI_red) THEN 853 CALL dbcsr_create(mat_L%matrix, template=mat_template) 854 ELSE 855 CALL dbcsr_get_info(mat_template, nblkrows_total=nblk) 856 857 CALL calculate_equal_blk_size(row_blk_size, dimen_RI_red, nblk) 858 859 CALL dbcsr_create(mat_L%matrix, template=mat_template, row_blk_size=row_blk_size) 860 861 DEALLOCATE (row_blk_size) 862 END IF 863 IF (.NOT. (do_kpoints)) THEN 864 CALL copy_fm_to_dbcsr(fm_mat_L(1, 1)%matrix, mat_L%matrix) 865 END IF 866 867 CALL timestop(handle) 868 869 END SUBROUTINE reorder_mat_L 870 871! ************************************************************************************************** 872!> \brief ... 873!> \param blk_size_new ... 874!> \param dimen_RI_red ... 875!> \param nblk ... 876! ************************************************************************************************** 877 SUBROUTINE calculate_equal_blk_size(blk_size_new, dimen_RI_red, nblk) 878 INTEGER, DIMENSION(:), POINTER :: blk_size_new 879 INTEGER, INTENT(IN) :: dimen_RI_red, nblk 880 881 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_equal_blk_size', & 882 routineP = moduleN//':'//routineN 883 884 INTEGER :: col_per_blk, remainder 885 886 NULLIFY (blk_size_new) 887 ALLOCATE (blk_size_new(nblk)) 888 889 remainder = MOD(dimen_RI_red, nblk) 890 col_per_blk = dimen_RI_red/nblk 891 892 ! Determine a new distribution for the columns (corresponding to the number of columns) 893 IF (remainder > 0) blk_size_new(1:remainder) = col_per_blk + 1 894 blk_size_new(remainder + 1:nblk) = col_per_blk 895 896 END SUBROUTINE calculate_equal_blk_size 897 898! ************************************************************************************************** 899!> \brief ... 900!> \param para_env_sub ... 901!> \param do_mao ... 902!> \param do_gw_im_time ... 903!> \param mat_3c_overl_int ... 904!> \param mat_3c_overl_int_mao_for_occ ... 905!> \param mat_3c_overl_int_mao_for_virt ... 906!> \param do_dbcsr_t ... 907!> \param eps_filter ... 908!> \param cut_RI ... 909!> \param cut_memory ... 910!> \param num_3c_repl ... 911!> \param my_group_L_sizes_im_time ... 912!> \param non_zero_blocks_3c_cut_col ... 913!> \param ends_array_cm ... 914!> \param ends_array_cm_mao_occ ... 915!> \param ends_array_cm_mao_virt ... 916!> \param starts_array_cm ... 917!> \param starts_array_cm_mao_occ ... 918!> \param starts_array_cm_mao_virt ... 919!> \param do_kpoints_cubic_RPA ... 920!> \param needed_cutRI_mem_R1vec_R2vec_for_kp ... 921!> \param mat_3c_overl_int_cut ... 922!> \param mat_3c_overl_int_mao_for_occ_cut ... 923!> \param mat_3c_overl_int_mao_for_virt_cut ... 924! ************************************************************************************************** 925 SUBROUTINE get_mat_3c_overl_int_cut(para_env_sub, do_mao, do_gw_im_time, mat_3c_overl_int, mat_3c_overl_int_mao_for_occ, & 926 mat_3c_overl_int_mao_for_virt, do_dbcsr_t, eps_filter, cut_RI, & 927 cut_memory, num_3c_repl, my_group_L_sizes_im_time, & 928 non_zero_blocks_3c_cut_col, ends_array_cm, ends_array_cm_mao_occ, ends_array_cm_mao_virt, & 929 starts_array_cm, starts_array_cm_mao_occ, starts_array_cm_mao_virt, & 930 do_kpoints_cubic_RPA, needed_cutRI_mem_R1vec_R2vec_for_kp, & 931 mat_3c_overl_int_cut, mat_3c_overl_int_mao_for_occ_cut, mat_3c_overl_int_mao_for_virt_cut) 932 TYPE(cp_para_env_type), POINTER :: para_env_sub 933 LOGICAL, INTENT(IN) :: do_mao, do_gw_im_time 934 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int, & 935 mat_3c_overl_int_mao_for_occ, & 936 mat_3c_overl_int_mao_for_virt 937 LOGICAL, INTENT(IN) :: do_dbcsr_t 938 REAL(KIND=dp), INTENT(IN) :: eps_filter 939 INTEGER, INTENT(IN) :: cut_RI, cut_memory 940 INTEGER, INTENT(OUT) :: num_3c_repl 941 INTEGER, DIMENSION(:), INTENT(IN) :: my_group_L_sizes_im_time 942 INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :), & 943 INTENT(OUT) :: non_zero_blocks_3c_cut_col 944 INTEGER, DIMENSION(:), POINTER :: ends_array_cm, ends_array_cm_mao_occ, & 945 ends_array_cm_mao_virt, starts_array_cm, starts_array_cm_mao_occ, starts_array_cm_mao_virt 946 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA 947 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :), INTENT(OUT) :: & 948 needed_cutRI_mem_R1vec_R2vec_for_kp 949 TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: mat_3c_overl_int_cut, & 950 mat_3c_overl_int_mao_for_occ_cut, mat_3c_overl_int_mao_for_virt_cut 951 952 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_3c_overl_int_cut', & 953 routineP = moduleN//':'//routineN 954 955 INTEGER :: handle 956 957 CALL timeset(routineN, handle) 958 959 ! This splitting allows to keep allocation for gw im.time outside from remaining allocation routines retaining slim interfaces 960 IF (.NOT. do_dbcsr_t) THEN 961 CALL setup_mat_for_mem_cut_3c(mat_3c_overl_int_cut, mat_3c_overl_int, cut_memory, cut_RI, & 962 starts_array_cm, ends_array_cm, my_group_L_sizes_im_time, eps_filter, & 963 do_kpoints_cubic_RPA, do_gw_im_time, num_3c_repl) 964 ENDIF 965 966 IF (do_mao) THEN 967 968 ! in setup_mat_for_mem_cut_3c, one deallocates mat_3c_overl_int, therefore for the beginning, deallocate 969 ! the mao 3c overlap integrals here 970 CALL setup_mat_for_mem_cut_3c(mat_3c_overl_int_mao_for_occ_cut, mat_3c_overl_int_mao_for_occ, & 971 cut_memory, cut_RI, starts_array_cm_mao_virt, ends_array_cm_mao_virt, & 972 my_group_L_sizes_im_time, eps_filter, do_kpoints_cubic_RPA, do_gw_im_time, & 973 num_3c_repl) 974 CALL setup_mat_for_mem_cut_3c(mat_3c_overl_int_mao_for_virt_cut, mat_3c_overl_int_mao_for_virt, & 975 cut_memory, cut_RI, starts_array_cm_mao_occ, ends_array_cm_mao_occ, & 976 my_group_L_sizes_im_time, eps_filter, do_kpoints_cubic_RPA, do_gw_im_time, & 977 num_3c_repl) 978 979 ELSE 980 IF (.NOT. do_dbcsr_t) THEN 981 982 mat_3c_overl_int_mao_for_occ_cut => mat_3c_overl_int_cut 983 mat_3c_overl_int_mao_for_virt_cut => mat_3c_overl_int_cut 984 ENDIF 985 END IF 986 987 IF (.NOT. do_dbcsr_t) THEN 988 CALL get_non_zero_blocks_3c_cut_col(mat_3c_overl_int_cut, para_env_sub, cut_RI, & 989 cut_memory, non_zero_blocks_3c_cut_col) 990 991 CALL check_sparsity_arrays_for_kp(needed_cutRI_mem_R1vec_R2vec_for_kp, & 992 mat_3c_overl_int_cut, cut_RI, cut_memory, para_env_sub, & 993 do_kpoints_cubic_RPA) 994 ENDIF 995 996 CALL timestop(handle) 997 998 END SUBROUTINE get_mat_3c_overl_int_cut 999 1000! ************************************************************************************************** 1001!> \brief ... 1002!> \param fm_mat_S ... 1003!> \param do_ri_sos_laplace_mp2 ... 1004!> \param first_cycle ... 1005!> \param count_ev_sc_GW ... 1006!> \param virtual ... 1007!> \param Eigenval ... 1008!> \param Eigenval_last ... 1009!> \param homo ... 1010!> \param omega ... 1011!> \param omega_old ... 1012!> \param jquad ... 1013!> \param mm_style ... 1014!> \param dimen_RI ... 1015!> \param dimen_ia ... 1016!> \param alpha ... 1017!> \param fm_mat_Q ... 1018!> \param fm_mat_Q_gemm ... 1019!> \param para_env_RPA ... 1020!> \param do_bse ... 1021!> \param fm_mat_Q_static_bse_gemm ... 1022!> \param RPA_proc_map ... 1023!> \param buffer_rec ... 1024!> \param buffer_send ... 1025!> \param number_of_send ... 1026!> \param map_send_size ... 1027!> \param map_rec_size ... 1028!> \param local_size_source ... 1029!> \param my_num_dgemm_call ... 1030!> \param my_flop_rate ... 1031! ************************************************************************************************** 1032 SUBROUTINE calc_mat_Q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, count_ev_sc_GW, virtual, Eigenval, Eigenval_last, & 1033 homo, omega, omega_old, jquad, mm_style, dimen_RI, dimen_ia, alpha, fm_mat_Q, fm_mat_Q_gemm, & 1034 para_env_RPA, do_bse, fm_mat_Q_static_bse_gemm, RPA_proc_map, buffer_rec, buffer_send, number_of_send, & 1035 map_send_size, map_rec_size, local_size_source, my_num_dgemm_call, my_flop_rate) 1036 TYPE(cp_fm_type), POINTER :: fm_mat_S 1037 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, first_cycle 1038 INTEGER, INTENT(IN) :: count_ev_sc_GW, virtual 1039 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval, Eigenval_last 1040 INTEGER, INTENT(IN) :: homo 1041 REAL(KIND=dp), INTENT(IN) :: omega, omega_old 1042 INTEGER, INTENT(IN) :: jquad, mm_style, dimen_RI, dimen_ia 1043 REAL(KIND=dp), INTENT(IN) :: alpha 1044 TYPE(cp_fm_type), POINTER :: fm_mat_Q, fm_mat_Q_gemm 1045 TYPE(cp_para_env_type), POINTER :: para_env_RPA 1046 LOGICAL, INTENT(IN) :: do_bse 1047 TYPE(cp_fm_type), POINTER :: fm_mat_Q_static_bse_gemm 1048 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: RPA_proc_map 1049 TYPE(integ_mat_buffer_type), ALLOCATABLE, & 1050 DIMENSION(:), INTENT(INOUT) :: buffer_rec, buffer_send 1051 INTEGER, INTENT(IN) :: number_of_send 1052 INTEGER, DIMENSION(0:para_env_RPA%num_pe-1), & 1053 INTENT(IN) :: map_send_size, map_rec_size 1054 INTEGER, DIMENSION(2, 0:para_env_RPA%num_pe-1), & 1055 INTENT(IN) :: local_size_source 1056 INTEGER, INTENT(INOUT) :: my_num_dgemm_call 1057 REAL(KIND=dp), INTENT(INOUT) :: my_flop_rate 1058 1059 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mat_Q', routineP = moduleN//':'//routineN 1060 1061 INTEGER :: handle 1062 1063 CALL timeset(routineN, handle) 1064 1065 IF (do_ri_sos_laplace_mp2) THEN 1066 ! the first index of tau_tj starts with 0 (see mp2_weights) 1067 CALL calc_fm_mat_S_laplace(fm_mat_S, first_cycle, homo, virtual, Eigenval, omega, omega_old) 1068 1069 ELSE 1070 CALL calc_fm_mat_S_rpa(fm_mat_S, first_cycle, count_ev_sc_GW, virtual, Eigenval, Eigenval_last, & 1071 homo, omega, omega_old) 1072 END IF 1073 1074 CALL contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, para_env_RPA, & 1075 my_num_dgemm_call, fm_mat_Q, RPA_proc_map, buffer_rec, buffer_send, number_of_send, & 1076 map_send_size, map_rec_size, local_size_source, my_flop_rate) 1077 1078 IF (do_bse .AND. jquad == 1) THEN 1079 CALL cp_fm_to_fm(fm_mat_Q_gemm, fm_mat_Q_static_bse_gemm) 1080 END IF 1081 CALL timestop(handle) 1082 1083 END SUBROUTINE calc_mat_Q 1084 1085! ************************************************************************************************** 1086!> \brief ... 1087!> \param fm_mat_S ... 1088!> \param first_cycle ... 1089!> \param count_ev_sc_GW ... 1090!> \param virtual ... 1091!> \param Eigenval ... 1092!> \param Eigenval_last ... 1093!> \param homo ... 1094!> \param omega ... 1095!> \param omega_old ... 1096! ************************************************************************************************** 1097 SUBROUTINE calc_fm_mat_S_rpa(fm_mat_S, first_cycle, count_ev_sc_GW, virtual, Eigenval, Eigenval_last, homo, & 1098 omega, omega_old) 1099 TYPE(cp_fm_type), POINTER :: fm_mat_S 1100 LOGICAL, INTENT(IN) :: first_cycle 1101 INTEGER, INTENT(IN) :: count_ev_sc_GW, virtual 1102 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval, Eigenval_last 1103 INTEGER, INTENT(IN) :: homo 1104 REAL(KIND=dp), INTENT(IN) :: omega, omega_old 1105 1106 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_fm_mat_S_rpa', & 1107 routineP = moduleN//':'//routineN 1108 1109 INTEGER :: avirt, handle, i_global, iiB, iocc, & 1110 j_global, jjB, ncol_local, nrow_local 1111 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 1112 REAL(KIND=dp) :: eigen_diff 1113 1114 CALL timeset(routineN, handle) 1115 1116 ! get info of fm_mat_S 1117 CALL cp_fm_get_info(matrix=fm_mat_S, & 1118 nrow_local=nrow_local, & 1119 ncol_local=ncol_local, & 1120 row_indices=row_indices, & 1121 col_indices=col_indices) 1122 1123 ! remove eigenvalue part of S matrix from the last eigenvalue self-c. cycle 1124 IF (first_cycle .AND. count_ev_sc_GW > 1) THEN 1125!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,iocc,avirt,eigen_diff,i_global,j_global) & 1126!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Eigenval_last,fm_mat_S,virtual,homo,omega_old) 1127 DO jjB = 1, ncol_local 1128 j_global = col_indices(jjB) 1129 DO iiB = 1, nrow_local 1130 i_global = row_indices(iiB) 1131 1132 iocc = MAX(1, i_global - 1)/virtual + 1 1133 avirt = i_global - (iocc - 1)*virtual 1134 eigen_diff = Eigenval_last(avirt + homo) - Eigenval_last(iocc) 1135 1136 fm_mat_S%local_data(iiB, jjB) = fm_mat_S%local_data(iiB, jjB)/ & 1137 SQRT(eigen_diff/(eigen_diff**2 + omega_old**2)) 1138 1139 END DO 1140 END DO 1141 1142 END IF 1143 1144 ! update G matrix with the new value of omega 1145 IF (first_cycle) THEN 1146 ! In this case just update the matrix (symmetric form) with 1147 ! SQRT((epsi_a-epsi_i)/((epsi_a-epsi_i)**2+omega**2)) 1148!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,iocc,avirt,eigen_diff,i_global,j_global) & 1149!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Eigenval,fm_mat_S,virtual,homo,omega) 1150 DO jjB = 1, ncol_local 1151 j_global = col_indices(jjB) 1152 DO iiB = 1, nrow_local 1153 i_global = row_indices(iiB) 1154 1155 iocc = MAX(1, i_global - 1)/virtual + 1 1156 avirt = i_global - (iocc - 1)*virtual 1157 eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc) 1158 1159 fm_mat_S%local_data(iiB, jjB) = fm_mat_S%local_data(iiB, jjB)* & 1160 SQRT(eigen_diff/(eigen_diff**2 + omega**2)) 1161 1162 END DO 1163 END DO 1164 ELSE 1165 ! In this case the update has to remove the old omega component thus 1166 ! SQRT(((epsi_a-epsi_i)**2+omega_old**2)/((epsi_a-epsi_i)**2+omega**2)) 1167!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,iocc,avirt,eigen_diff,i_global,j_global) & 1168!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Eigenval,& 1169!$OMP fm_mat_S,virtual,homo,omega,omega_old) 1170 DO jjB = 1, ncol_local 1171 j_global = col_indices(jjB) 1172 DO iiB = 1, nrow_local 1173 i_global = row_indices(iiB) 1174 1175 iocc = MAX(1, i_global - 1)/virtual + 1 1176 avirt = i_global - (iocc - 1)*virtual 1177 eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc) 1178 1179 fm_mat_S%local_data(iiB, jjB) = fm_mat_S%local_data(iiB, jjB)* & 1180 SQRT((eigen_diff**2 + omega_old**2)/(eigen_diff**2 + omega**2)) 1181 1182 END DO 1183 END DO 1184 END IF 1185 1186 CALL timestop(handle) 1187 1188 END SUBROUTINE calc_fm_mat_S_rpa 1189 1190! ************************************************************************************************** 1191!> \brief ... 1192!> \param mm_style ... 1193!> \param dimen_RI ... 1194!> \param dimen_ia ... 1195!> \param alpha ... 1196!> \param fm_mat_S ... 1197!> \param fm_mat_Q_gemm ... 1198!> \param para_env_RPA ... 1199!> \param my_num_dgemm_call ... 1200!> \param fm_mat_Q ... 1201!> \param RPA_proc_map ... 1202!> \param buffer_rec ... 1203!> \param buffer_send ... 1204!> \param number_of_send ... 1205!> \param map_send_size ... 1206!> \param map_rec_size ... 1207!> \param local_size_source ... 1208!> \param my_flop_rate ... 1209! ************************************************************************************************** 1210 SUBROUTINE contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, para_env_RPA, my_num_dgemm_call, & 1211 fm_mat_Q, RPA_proc_map, buffer_rec, buffer_send, number_of_send, & 1212 map_send_size, map_rec_size, local_size_source, my_flop_rate) 1213 1214 INTEGER, INTENT(IN) :: mm_style, dimen_RI, dimen_ia 1215 REAL(KIND=dp), INTENT(IN) :: alpha 1216 TYPE(cp_fm_type), POINTER :: fm_mat_S, fm_mat_Q_gemm 1217 TYPE(cp_para_env_type), POINTER :: para_env_RPA 1218 INTEGER, INTENT(INOUT) :: my_num_dgemm_call 1219 TYPE(cp_fm_type), POINTER :: fm_mat_Q 1220 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: RPA_proc_map 1221 TYPE(integ_mat_buffer_type), DIMENSION(:), & 1222 INTENT(INOUT) :: buffer_rec, buffer_send 1223 INTEGER, INTENT(IN) :: number_of_send 1224 INTEGER, DIMENSION(0:para_env_RPA%num_pe-1), & 1225 INTENT(IN) :: map_send_size, map_rec_size 1226 INTEGER, DIMENSION(2, 0:para_env_RPA%num_pe-1), & 1227 INTENT(IN) :: local_size_source 1228 REAL(KIND=dp), INTENT(INOUT) :: my_flop_rate 1229 1230 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_S_to_Q', & 1231 routineP = moduleN//':'//routineN 1232 1233 INTEGER :: handle 1234 REAL(KIND=dp) :: actual_flop_rate, t_end, t_start 1235 1236 CALL timeset(routineN, handle) 1237 1238 t_start = m_walltime() 1239 SELECT CASE (mm_style) 1240 CASE (wfc_mm_style_gemm) 1241 ! waste-fully computes the full symmetrix matrix, but maybe faster than cp_fm_syrk for optimized cp_fm_gemm 1242 CALL cp_gemm(transa="T", transb="N", m=dimen_RI, n=dimen_RI, k=dimen_ia, alpha=alpha, & 1243 matrix_a=fm_mat_S, matrix_b=fm_mat_S, beta=0.0_dp, & 1244 matrix_c=fm_mat_Q_gemm) 1245 CASE (wfc_mm_style_syrk) 1246 ! will only compute the upper half of the matrix, which is fine, since we only use it for cholesky later 1247 CALL cp_fm_syrk(uplo='U', trans='T', k=dimen_ia, alpha=alpha, matrix_a=fm_mat_S, & 1248 ia=1, ja=1, beta=0.0_dp, matrix_c=fm_mat_Q_gemm) 1249 CASE DEFAULT 1250 CPABORT("") 1251 END SELECT 1252 t_end = m_walltime() 1253 actual_flop_rate = 2.0_dp*REAL(dimen_ia, KIND=dp)*dimen_RI*REAL(dimen_RI, KIND=dp)/(MAX(0.01_dp, t_end - t_start)) 1254 IF (para_env_RPA%mepos == 0) my_flop_rate = my_flop_rate + actual_flop_rate 1255 my_num_dgemm_call = my_num_dgemm_call + 1 1256 1257 ! copy/redistribute fm_mat_Q_gemm to fm_mat_Q 1258 CALL cp_fm_set_all(matrix=fm_mat_Q, alpha=0.0_dp) 1259 CALL fm_redistribute(fm_mat_Q_gemm, fm_mat_Q, RPA_proc_map, buffer_rec, buffer_send, & 1260 number_of_send, map_send_size, map_rec_size, local_size_source, para_env_RPA) 1261 1262 CALL timestop(handle) 1263 1264 END SUBROUTINE contract_S_to_Q 1265 1266! ************************************************************************************************** 1267!> \brief ... 1268!> \param dimen_RI ... 1269!> \param trace_Qomega ... 1270!> \param fm_mat_Q ... 1271!> \param para_env_RPA ... 1272! ************************************************************************************************** 1273 SUBROUTINE RPA_postprocessing_start(dimen_RI, trace_Qomega, fm_mat_Q, para_env_RPA) 1274 1275 INTEGER, INTENT(IN) :: dimen_RI 1276 REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(OUT) :: trace_Qomega 1277 TYPE(cp_fm_type), POINTER :: fm_mat_Q 1278 TYPE(cp_para_env_type), POINTER :: para_env_RPA 1279 1280 CHARACTER(LEN=*), PARAMETER :: routineN = 'RPA_postprocessing_start', & 1281 routineP = moduleN//':'//routineN 1282 1283 INTEGER :: handle, i_global, iiB, j_global, jjB, & 1284 ncol_local, nrow_local 1285 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 1286 1287 CALL timeset(routineN, handle) 1288 1289 CALL cp_fm_get_info(matrix=fm_mat_Q, & 1290 nrow_local=nrow_local, & 1291 ncol_local=ncol_local, & 1292 row_indices=row_indices, & 1293 col_indices=col_indices) 1294 1295 ! calculate the trace of Q and add 1 on the diagonal 1296 trace_Qomega = 0.0_dp 1297!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 1298!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,trace_Qomega,fm_mat_Q,dimen_RI) 1299 DO jjB = 1, ncol_local 1300 j_global = col_indices(jjB) 1301 DO iiB = 1, nrow_local 1302 i_global = row_indices(iiB) 1303 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 1304 trace_Qomega(i_global) = fm_mat_Q%local_data(iiB, jjB) 1305 fm_mat_Q%local_data(iiB, jjB) = fm_mat_Q%local_data(iiB, jjB) + 1.0_dp 1306 END IF 1307 END DO 1308 END DO 1309 CALL mp_sum(trace_Qomega, para_env_RPA%group) 1310 1311 CALL timestop(handle) 1312 1313 END SUBROUTINE RPA_postprocessing_start 1314 1315! ************************************************************************************************** 1316!> \brief ... 1317!> \param dimen_RI ... 1318!> \param trace_Qomega ... 1319!> \param fm_mat_Q ... 1320!> \param para_env_RPA ... 1321!> \param Erpa ... 1322!> \param wjquad ... 1323! ************************************************************************************************** 1324 SUBROUTINE RPA_postprocessing_nokp(dimen_RI, trace_Qomega, fm_mat_Q, para_env_RPA, Erpa, wjquad) 1325 1326 INTEGER, INTENT(IN) :: dimen_RI 1327 REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(IN) :: trace_Qomega 1328 TYPE(cp_fm_type), POINTER :: fm_mat_Q 1329 TYPE(cp_para_env_type), POINTER :: para_env_RPA 1330 REAL(KIND=dp), INTENT(INOUT) :: Erpa 1331 REAL(KIND=dp), INTENT(IN) :: wjquad 1332 1333 CHARACTER(LEN=*), PARAMETER :: routineN = 'RPA_postprocessing_nokp', & 1334 routineP = moduleN//':'//routineN 1335 1336 INTEGER :: handle, i_global, iiB, info_chol, & 1337 j_global, jjB, ncol_local, nrow_local 1338 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 1339 REAL(KIND=dp) :: FComega 1340 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Q_log 1341 1342 CALL timeset(routineN, handle) 1343 1344 CALL cp_fm_get_info(matrix=fm_mat_Q, & 1345 nrow_local=nrow_local, & 1346 ncol_local=ncol_local, & 1347 row_indices=row_indices, & 1348 col_indices=col_indices) 1349 1350 ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition 1351 CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q, n=dimen_RI, info_out=info_chol) 1352 CPASSERT(info_chol == 0) 1353 1354 ALLOCATE (Q_log(dimen_RI)) 1355 Q_log = 0.0_dp 1356!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 1357!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Q_log,fm_mat_Q,dimen_RI) 1358 DO jjB = 1, ncol_local 1359 j_global = col_indices(jjB) 1360 DO iiB = 1, nrow_local 1361 i_global = row_indices(iiB) 1362 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 1363 Q_log(i_global) = 2.0_dp*LOG(fm_mat_Q%local_data(iiB, jjB)) 1364 END IF 1365 END DO 1366 END DO 1367 CALL mp_sum(Q_log, para_env_RPA%group) 1368 1369 FComega = 0.0_dp 1370 DO iiB = 1, dimen_RI 1371 IF (MODULO(iiB, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE 1372 FComega = FComega + (Q_log(iiB) - trace_Qomega(iiB))/2.0_dp 1373 END DO 1374 Erpa = Erpa + FComega*wjquad 1375 1376 DEALLOCATE (Q_log) 1377 1378 CALL timestop(handle) 1379 1380 END SUBROUTINE RPA_postprocessing_nokp 1381 1382! ************************************************************************************************** 1383!> \brief ... 1384!> \param mat_3c_overl_int ... 1385!> \param para_env_sub ... 1386!> \param cut_RI ... 1387!> \param non_zero_blocks_3c ... 1388! ************************************************************************************************** 1389 SUBROUTINE get_non_zero_blocks_3c(mat_3c_overl_int, para_env_sub, cut_RI, non_zero_blocks_3c) 1390 1391 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int 1392 TYPE(cp_para_env_type), POINTER :: para_env_sub 1393 INTEGER, INTENT(IN) :: cut_RI 1394 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & 1395 INTENT(OUT) :: non_zero_blocks_3c 1396 1397 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_non_zero_blocks_3c', & 1398 routineP = moduleN//':'//routineN 1399 1400 INTEGER :: blk, block_counter, col, handle, i_cell, i_cut_RI, iblk, imepos, j_cell, & 1401 maxlength, maxlength_tmp, nblkrows_total, num_cells_3c, row 1402 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: non_zero_blocks_3c_tmp 1403 REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block 1404 TYPE(dbcsr_iterator_type) :: iter 1405 1406 CALL timeset(routineN, handle) 1407 1408 num_cells_3c = SIZE(mat_3c_overl_int, 2) 1409 CPASSERT(num_cells_3c == SIZE(mat_3c_overl_int, 3)) 1410 1411 CALL dbcsr_get_info(mat_3c_overl_int(1, 1, 1)%matrix, nblkrows_total=nblkrows_total) 1412 1413 ALLOCATE (non_zero_blocks_3c_tmp(1:cut_RI, 1:nblkrows_total, 0:(para_env_sub%num_pe - 1))) 1414 non_zero_blocks_3c_tmp = 0 1415 1416 DO i_cut_RI = 1, cut_RI 1417 1418 DO i_cell = 1, num_cells_3c 1419 1420 DO j_cell = 1, num_cells_3c 1421 1422 CALL dbcsr_iterator_start(iter, mat_3c_overl_int(i_cut_RI, i_cell, j_cell)%matrix) 1423 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1424 CALL dbcsr_iterator_next_block(iter, row, col, data_block, blk) 1425 1426 non_zero_blocks_3c_tmp(i_cut_RI, row, para_env_sub%mepos) = 1 1427 1428 ENDDO 1429 1430 CALL dbcsr_iterator_stop(iter) 1431 1432 END DO 1433 1434 END DO 1435 1436 END DO 1437 1438 CALL mp_sum(non_zero_blocks_3c_tmp, para_env_sub%group) 1439 1440 maxlength = 0 1441 1442 DO imepos = 0, para_env_sub%num_pe - 1 1443 DO i_cut_RI = 1, cut_RI 1444 maxlength_tmp = 0 1445 DO iblk = 1, nblkrows_total 1446 IF (non_zero_blocks_3c_tmp(i_cut_RI, iblk, imepos) .NE. 0) THEN 1447 maxlength_tmp = maxlength_tmp + 1 1448 END IF 1449 END DO 1450 IF (maxlength_tmp > maxlength) THEN 1451 maxlength = maxlength_tmp 1452 END IF 1453 END DO 1454 END DO 1455 1456 ! save memory with smaller non_zero_blocks_3c 1457 ALLOCATE (non_zero_blocks_3c(1:cut_RI, 1:maxlength, 0:(para_env_sub%num_pe - 1))) 1458 non_zero_blocks_3c = 0 1459 1460 DO imepos = 0, para_env_sub%num_pe - 1 1461 DO i_cut_RI = 1, cut_RI 1462 block_counter = 0 1463 DO iblk = 1, nblkrows_total 1464 IF (non_zero_blocks_3c_tmp(i_cut_RI, iblk, imepos) .NE. 0) THEN 1465 block_counter = block_counter + 1 1466 non_zero_blocks_3c(i_cut_RI, block_counter, imepos) = iblk 1467 END IF 1468 END DO 1469 END DO 1470 END DO 1471 1472 DEALLOCATE (non_zero_blocks_3c_tmp) 1473 1474 CALL timestop(handle) 1475 1476 END SUBROUTINE get_non_zero_blocks_3c 1477 1478! ************************************************************************************************** 1479!> \brief ... 1480!> \param mat_3c_overl_int_cut_col ... 1481!> \param para_env_sub ... 1482!> \param cut_RI ... 1483!> \param cut_memory ... 1484!> \param non_zero_blocks_3c_cut ... 1485! ************************************************************************************************** 1486 SUBROUTINE get_non_zero_blocks_3c_cut_col(mat_3c_overl_int_cut_col, para_env_sub, cut_RI, cut_memory, & 1487 non_zero_blocks_3c_cut) 1488 1489 TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: mat_3c_overl_int_cut_col 1490 TYPE(cp_para_env_type), POINTER :: para_env_sub 1491 INTEGER, INTENT(IN) :: cut_RI, cut_memory 1492 INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :), & 1493 INTENT(OUT) :: non_zero_blocks_3c_cut 1494 1495 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_non_zero_blocks_3c_cut_col', & 1496 routineP = moduleN//':'//routineN 1497 1498 INTEGER :: blk, block_counter, col, handle, i_cell, i_cut_RI, i_mem, iblk, imepos, j_cell, & 1499 maxlength, maxlength_tmp, nblkrows_total, nblkrows_total_max, num_3c_repl, row 1500 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: non_zero_blocks_3c_tmp 1501 REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block 1502 TYPE(dbcsr_iterator_type) :: iter 1503 1504 CALL timeset(routineN, handle) 1505 1506 nblkrows_total_max = 0 1507 1508 num_3c_repl = SIZE(mat_3c_overl_int_cut_col, 3) 1509 ! quickly check that replica of fourth component matches the one of third component 1510 CPASSERT(num_3c_repl == SIZE(mat_3c_overl_int_cut_col, 4)) 1511 1512 DO i_cut_RI = 1, cut_RI 1513 DO i_mem = 1, cut_memory 1514 CALL dbcsr_get_info(mat_3c_overl_int_cut_col(i_cut_RI, i_mem, 1, 1)%matrix, & 1515 nblkrows_total=nblkrows_total) 1516 IF (nblkrows_total > nblkrows_total_max) THEN 1517 nblkrows_total_max = nblkrows_total 1518 END IF 1519 END DO 1520 END DO 1521 1522 ALLOCATE (non_zero_blocks_3c_tmp(1:cut_RI, 1:nblkrows_total_max, 0:(para_env_sub%num_pe - 1))) 1523 non_zero_blocks_3c_tmp = 0 1524 1525 maxlength = 0 1526 1527 ! first, determine maxlength 1528 DO i_mem = 1, cut_memory 1529 1530 DO i_cut_RI = 1, cut_RI 1531 1532 DO i_cell = 1, num_3c_repl 1533 1534 DO j_cell = 1, num_3c_repl 1535 1536 CALL dbcsr_iterator_start(iter, mat_3c_overl_int_cut_col(i_cut_RI, i_mem, i_cell, j_cell)%matrix) 1537 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1538 CALL dbcsr_iterator_next_block(iter, row, col, data_block, blk) 1539 1540 non_zero_blocks_3c_tmp(i_cut_RI, col, para_env_sub%mepos) = 1 1541 1542 ENDDO 1543 1544 CALL dbcsr_iterator_stop(iter) 1545 1546 END DO 1547 1548 END DO 1549 1550 END DO ! cut_RI 1551 1552 CALL mp_sum(non_zero_blocks_3c_tmp, para_env_sub%group) 1553 1554 maxlength_tmp = 0 1555 1556 DO imepos = 0, para_env_sub%num_pe - 1 1557 DO i_cut_RI = 1, cut_RI 1558 maxlength_tmp = 0 1559 DO iblk = 1, nblkrows_total 1560 IF (non_zero_blocks_3c_tmp(i_cut_RI, iblk, imepos) .NE. 0) THEN 1561 maxlength_tmp = maxlength_tmp + 1 1562 END IF 1563 END DO 1564 IF (maxlength_tmp > maxlength) THEN 1565 maxlength = maxlength_tmp 1566 END IF 1567 END DO 1568 END DO 1569 1570 non_zero_blocks_3c_tmp = 0 1571 1572 END DO ! i_mem 1573 ! end determine maxlength 1574 1575 ! save memory with a smaller non_zero_blocks_3c_cut 1576 ALLOCATE (non_zero_blocks_3c_cut(1:cut_RI, 1:maxlength, 0:(para_env_sub%num_pe - 1), 1:cut_memory)) 1577 non_zero_blocks_3c_cut = 0 1578 1579 ! now, fill non_zero_blocks_3c_cut 1580 DO i_mem = 1, cut_memory 1581 1582 DO i_cut_RI = 1, cut_RI 1583 1584 DO i_cell = 1, num_3c_repl 1585 1586 DO j_cell = 1, num_3c_repl 1587 1588 CALL dbcsr_iterator_start(iter, mat_3c_overl_int_cut_col(i_cut_RI, i_mem, i_cell, j_cell)%matrix) 1589 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1590 CALL dbcsr_iterator_next_block(iter, row, col, data_block, blk) 1591 1592 non_zero_blocks_3c_tmp(i_cut_RI, col, para_env_sub%mepos) = 1 1593 1594 ENDDO 1595 1596 CALL dbcsr_iterator_stop(iter) 1597 1598 END DO 1599 END DO 1600 1601 END DO ! cut_RI 1602 1603 CALL mp_sum(non_zero_blocks_3c_tmp, para_env_sub%group) 1604 1605 DO imepos = 0, para_env_sub%num_pe - 1 1606 DO i_cut_RI = 1, cut_RI 1607 block_counter = 0 1608 DO iblk = 1, nblkrows_total 1609 IF (non_zero_blocks_3c_tmp(i_cut_RI, iblk, imepos) .NE. 0) THEN 1610 block_counter = block_counter + 1 1611 non_zero_blocks_3c_cut(i_cut_RI, block_counter, imepos, i_mem) = iblk 1612 END IF 1613 END DO 1614 END DO 1615 END DO 1616 1617 non_zero_blocks_3c_tmp = 0 1618 1619 END DO ! i_mem 1620 ! end fill non_zero_blocks_3c_cut 1621 1622 DEALLOCATE (non_zero_blocks_3c_tmp) 1623 1624 CALL timestop(handle) 1625 1626 END SUBROUTINE get_non_zero_blocks_3c_cut_col 1627 1628! ************************************************************************************************** 1629!> \brief ... 1630!> \param needed_cutRI_mem_R1vec_R2vec_for_kp ... 1631!> \param mat_3c_overl_int_cut ... 1632!> \param cut_RI ... 1633!> \param cut_memory ... 1634!> \param para_env_sub ... 1635!> \param do_kpoints_cubic_RPA ... 1636! ************************************************************************************************** 1637 SUBROUTINE check_sparsity_arrays_for_kp(needed_cutRI_mem_R1vec_R2vec_for_kp, & 1638 mat_3c_overl_int_cut, cut_RI, cut_memory, para_env_sub, & 1639 do_kpoints_cubic_RPA) 1640 1641 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :), INTENT(OUT) :: & 1642 needed_cutRI_mem_R1vec_R2vec_for_kp 1643 TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: mat_3c_overl_int_cut 1644 INTEGER, INTENT(IN) :: cut_RI, cut_memory 1645 TYPE(cp_para_env_type), POINTER :: para_env_sub 1646 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA 1647 1648 CHARACTER(LEN=*), PARAMETER :: routineN = 'check_sparsity_arrays_for_kp', & 1649 routineP = moduleN//':'//routineN 1650 1651 INTEGER :: handle, i_cell, i_cut_RI, i_mem, j_cell, & 1652 num_cells_3c 1653 REAL(KIND=dp) :: occ_local, occ_local_sum_j_cell 1654 1655 CALL timeset(routineN, handle) 1656 1657 num_cells_3c = SIZE(mat_3c_overl_int_cut, 3) 1658 1659 ALLOCATE (needed_cutRI_mem_R1vec_R2vec_for_kp(cut_RI, cut_memory, 0:num_cells_3c, 0:num_cells_3c)) 1660 needed_cutRI_mem_R1vec_R2vec_for_kp(:, :, :, :) = .TRUE. 1661 1662 IF (do_kpoints_cubic_RPA) THEN 1663 1664 DO i_mem = 1, cut_memory 1665 1666 DO i_cell = 1, num_cells_3c 1667 1668 occ_local_sum_j_cell = 0.0_dp 1669 1670 DO i_cut_RI = 1, cut_RI 1671 1672 DO j_cell = 1, num_cells_3c 1673 1674 occ_local = dbcsr_get_occupation(mat_3c_overl_int_cut(i_cut_RI, i_mem, & 1675 j_cell, i_cell)%matrix) 1676 1677 CALL mp_sum(occ_local, para_env_sub%group) 1678 1679 IF (occ_local < 1E-100_dp) THEN 1680 1681 needed_cutRI_mem_R1vec_R2vec_for_kp(i_cut_RI, i_mem, j_cell, i_cell) = .FALSE. 1682 1683 END IF 1684 1685 occ_local_sum_j_cell = occ_local_sum_j_cell + occ_local 1686 1687 END DO ! j_cell 1688 1689 END DO ! i_cut_RI 1690 1691 END DO ! i_cell 1692 END DO ! i_mem 1693 1694 END IF 1695 1696 CALL timestop(handle) 1697 1698 END SUBROUTINE check_sparsity_arrays_for_kp 1699 1700! ************************************************************************************************** 1701!> \brief ... 1702!> \param fm_struct_sub_kp ... 1703!> \param para_env ... 1704!> \param nkp ... 1705!> \param dimen_RI ... 1706!> \param ikp_local ... 1707!> \param first_ikp_local ... 1708!> \param do_kpoints_cubic_RPA ... 1709! ************************************************************************************************** 1710 SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, nkp, dimen_RI, & 1711 ikp_local, first_ikp_local, do_kpoints_cubic_RPA) 1712 TYPE(cp_fm_struct_type), POINTER :: fm_struct_sub_kp 1713 TYPE(cp_para_env_type), POINTER :: para_env 1714 INTEGER, INTENT(IN) :: nkp, dimen_RI 1715 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ikp_local 1716 INTEGER, INTENT(OUT) :: first_ikp_local 1717 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA 1718 1719 CHARACTER(len=*), PARAMETER :: routineN = 'get_sub_para_kp', & 1720 routineP = moduleN//':'//routineN 1721 1722 INTEGER :: color_sub_kp, comm_sub_kp, handle, ikp, & 1723 num_proc_per_kp 1724 TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub_kp 1725 TYPE(cp_para_env_type), POINTER :: para_env_sub_kp 1726 1727 CALL timeset(routineN, handle) 1728 1729 IF (nkp > para_env%num_pe .OR. do_kpoints_cubic_RPA) THEN 1730 ! we have all kpoints on every processpr 1731 num_proc_per_kp = para_env%num_pe 1732 ELSE 1733 ! we have only one kpoint per group 1734 num_proc_per_kp = para_env%num_pe/nkp 1735 END IF 1736 1737 color_sub_kp = MOD(para_env%mepos/num_proc_per_kp, nkp) 1738 CALL mp_comm_split_direct(para_env%group, comm_sub_kp, color_sub_kp) 1739 NULLIFY (para_env_sub_kp) 1740 CALL cp_para_env_create(para_env_sub_kp, comm_sub_kp) 1741 NULLIFY (blacs_env_sub_kp) 1742 CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp) 1743 1744 NULLIFY (fm_struct_sub_kp) 1745 CALL cp_fm_struct_create(fm_struct_sub_kp, context=blacs_env_sub_kp, nrow_global=dimen_RI, & 1746 ncol_global=dimen_RI, para_env=para_env_sub_kp) 1747 1748 CALL cp_para_env_release(para_env_sub_kp) 1749 1750 CALL cp_blacs_env_release(blacs_env_sub_kp) 1751 1752 ALLOCATE (ikp_local(nkp)) 1753 ikp_local = 0 1754 first_ikp_local = 1 1755 DO ikp = 1, nkp 1756 IF (nkp > para_env%num_pe .OR. do_kpoints_cubic_RPA .OR. ikp == color_sub_kp + 1) THEN 1757 ikp_local(ikp) = ikp 1758 first_ikp_local = ikp 1759 END IF 1760 END DO 1761 1762 CALL timestop(handle) 1763 1764 END SUBROUTINE get_sub_para_kp 1765 1766! ************************************************************************************************** 1767!> \brief ... 1768!> \param do_mao ... 1769!> \param do_dbcsr_t ... 1770!> \param do_ri_sos_laplace_mp2 ... 1771!> \param fm_mo_coeff_occ ... 1772!> \param fm_mo_coeff_virt ... 1773!> \param fm_scaled_dm_occ_tau ... 1774!> \param fm_scaled_dm_virt_tau ... 1775!> \param ikp_local ... 1776!> \param row_from_LLL ... 1777!> \param index_to_cell_3c ... 1778!> \param cell_to_index_3c ... 1779!> \param non_zero_blocks_3c ... 1780!> \param non_zero_blocks_3c_cut_col ... 1781!> \param do_ic_model ... 1782!> \param do_kpoints_cubic_RPA ... 1783!> \param do_kpoints_from_Gamma ... 1784!> \param do_ri_Sigma_x ... 1785!> \param cycle_due_to_sparse_dm ... 1786!> \param multiply_needed_occ ... 1787!> \param multiply_needed_virt ... 1788!> \param needed_cutRI_mem_R1vec_R2vec_for_kp ... 1789!> \param has_mat_P_blocks ... 1790!> \param buffer_mat_M ... 1791!> \param wkp_W ... 1792!> \param cfm_mat_Q ... 1793!> \param fm_mat_L ... 1794!> \param fm_mat_RI_global_work ... 1795!> \param fm_mat_work ... 1796!> \param fm_mo_coeff_occ_scaled ... 1797!> \param fm_mo_coeff_virt_scaled ... 1798!> \param mat_dm ... 1799!> \param mat_L ... 1800!> \param mat_M_P_munu_occ ... 1801!> \param mat_M_P_munu_virt ... 1802!> \param mat_P_global_copy ... 1803!> \param mat_SinvVSinv ... 1804!> \param mat_M_mu_Pnu_occ ... 1805!> \param mat_M_mu_Pnu_virt ... 1806!> \param mat_P_omega ... 1807!> \param mat_P_omega_kp ... 1808!> \param mat_dm_loc_occ_cut ... 1809!> \param mat_dm_loc_virt_cut ... 1810!> \param mat_3c_overl_int_cut ... 1811!> \param mat_3c_overl_int_mao_for_occ_cut ... 1812!> \param mat_3c_overl_int_mao_for_virt_cut ... 1813!> \param t_3c_overl_int ... 1814!> \param t_3c_M ... 1815!> \param t_3c_O ... 1816!> \param mat_dm_loc_occ ... 1817!> \param mat_dm_loc_virt ... 1818!> \param mat_work ... 1819!> \param fm_mo_coeff_occ_beta ... 1820!> \param fm_mo_coeff_virt_beta ... 1821!> \param mat_P_omega_beta ... 1822! ************************************************************************************************** 1823 SUBROUTINE dealloc_im_time(do_mao, do_dbcsr_t, do_ri_sos_laplace_mp2, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_scaled_dm_occ_tau, & 1824 fm_scaled_dm_virt_tau, ikp_local, row_from_LLL, index_to_cell_3c, & 1825 cell_to_index_3c, non_zero_blocks_3c, non_zero_blocks_3c_cut_col, do_ic_model, & 1826 do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, cycle_due_to_sparse_dm, & 1827 multiply_needed_occ, multiply_needed_virt, needed_cutRI_mem_R1vec_R2vec_for_kp, has_mat_P_blocks, & 1828 buffer_mat_M, wkp_W, cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, & 1829 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, & 1830 mat_P_global_copy, mat_SinvVSinv, mat_M_mu_Pnu_occ, mat_M_mu_Pnu_virt, mat_P_omega, mat_P_omega_kp, & 1831 mat_dm_loc_occ_cut, mat_dm_loc_virt_cut, mat_3c_overl_int_cut, mat_3c_overl_int_mao_for_occ_cut, & 1832 mat_3c_overl_int_mao_for_virt_cut, t_3c_overl_int, t_3c_M, t_3c_O, & 1833 mat_dm_loc_occ, mat_dm_loc_virt, mat_work, & 1834 fm_mo_coeff_occ_beta, fm_mo_coeff_virt_beta, mat_P_omega_beta) 1835 1836 LOGICAL, INTENT(IN) :: do_mao, do_dbcsr_t, do_ri_sos_laplace_mp2 1837 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_occ, fm_mo_coeff_virt, & 1838 fm_scaled_dm_occ_tau, & 1839 fm_scaled_dm_virt_tau 1840 INTEGER, DIMENSION(:), INTENT(IN) :: ikp_local 1841 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: row_from_LLL 1842 INTEGER, ALLOCATABLE, DIMENSION(:, :), & 1843 INTENT(INOUT) :: index_to_cell_3c 1844 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & 1845 INTENT(INOUT) :: cell_to_index_3c, non_zero_blocks_3c 1846 INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :), & 1847 INTENT(INOUT) :: non_zero_blocks_3c_cut_col 1848 LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA, & 1849 do_kpoints_from_Gamma, do_ri_Sigma_x 1850 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :), & 1851 INTENT(INOUT) :: cycle_due_to_sparse_dm, & 1852 multiply_needed_occ, & 1853 multiply_needed_virt 1854 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :), INTENT(INOUT) :: & 1855 needed_cutRI_mem_R1vec_R2vec_for_kp 1856 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), & 1857 INTENT(INOUT) :: has_mat_P_blocks 1858 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1859 INTENT(INOUT) :: buffer_mat_M 1860 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 1861 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 1862 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L 1863 TYPE(cp_fm_type), POINTER :: fm_mat_RI_global_work, fm_mat_work, & 1864 fm_mo_coeff_occ_scaled, & 1865 fm_mo_coeff_virt_scaled 1866 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_dm, mat_L, mat_M_P_munu_occ, & 1867 mat_M_P_munu_virt, mat_P_global_copy, & 1868 mat_SinvVSinv 1869 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_M_mu_Pnu_occ, mat_M_mu_Pnu_virt 1870 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega, mat_P_omega_kp 1871 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_dm_loc_occ_cut, mat_dm_loc_virt_cut 1872 TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: mat_3c_overl_int_cut, & 1873 mat_3c_overl_int_mao_for_occ_cut, mat_3c_overl_int_mao_for_virt_cut 1874 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_overl_int 1875 TYPE(dbcsr_t_type) :: t_3c_M 1876 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O 1877 TYPE(dbcsr_type), POINTER :: mat_dm_loc_occ, mat_dm_loc_virt, mat_work 1878 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_mo_coeff_occ_beta, & 1879 fm_mo_coeff_virt_beta 1880 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & 1881 POINTER :: mat_P_omega_beta 1882 1883 CHARACTER(LEN=*), PARAMETER :: routineN = 'dealloc_im_time', & 1884 routineP = moduleN//':'//routineN 1885 1886 INTEGER :: handle, i_size, j_size 1887 LOGICAL :: my_open_shell 1888 1889 CALL timeset(routineN, handle) 1890 1891 my_open_shell = .FALSE. 1892 IF (PRESENT(fm_mo_coeff_occ_beta) .AND. PRESENT(fm_mo_coeff_virt_beta) .AND. PRESENT(mat_P_omega_beta)) THEN 1893 my_open_shell = .TRUE. 1894 END IF 1895 1896 CALL cp_fm_release(fm_scaled_dm_occ_tau) 1897 CALL cp_fm_release(fm_scaled_dm_virt_tau) 1898 CALL cp_fm_release(fm_mo_coeff_occ) 1899 CALL cp_fm_release(fm_mo_coeff_virt) 1900 CALL cp_fm_release(fm_mo_coeff_occ_scaled) 1901 CALL cp_fm_release(fm_mo_coeff_virt_scaled) 1902 1903 IF (my_open_shell) THEN 1904 CALL cp_fm_release(fm_mo_coeff_occ_beta) 1905 CALL cp_fm_release(fm_mo_coeff_virt_beta) 1906 END IF 1907 1908 DO i_size = 1, SIZE(fm_mat_L, 1) 1909 DO j_size = 1, SIZE(fm_mat_L, 2) 1910 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 1911 IF (ANY(ikp_local(:) == i_size)) THEN 1912 CALL cp_fm_release(fm_mat_L(i_size, j_size)%matrix) 1913 END IF 1914 ELSE 1915 CALL cp_fm_release(fm_mat_L(i_size, j_size)%matrix) 1916 END IF 1917 END DO 1918 END DO 1919 DEALLOCATE (fm_mat_L) 1920 CALL cp_fm_release(fm_mat_work) 1921 1922 IF (.NOT. do_dbcsr_t) THEN 1923 CALL dbcsr_release_p(mat_dm_loc_occ) 1924 CALL dbcsr_release_p(mat_dm_loc_virt) 1925 1926 CALL dbcsr_deallocate_matrix_set(mat_dm_loc_occ_cut) 1927 CALL dbcsr_deallocate_matrix_set(mat_dm_loc_virt_cut) 1928 1929 CALL dbcsr_release(mat_M_P_munu_occ%matrix) 1930 DEALLOCATE (mat_M_P_munu_occ%matrix) 1931 1932 CALL dbcsr_release(mat_M_P_munu_virt%matrix) 1933 DEALLOCATE (mat_M_P_munu_virt%matrix) 1934 1935 CALL dbcsr_release(mat_P_global_copy%matrix) 1936 DEALLOCATE (mat_P_global_copy%matrix) 1937 ENDIF 1938 1939 IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN 1940 CALL dbcsr_release(mat_work) 1941 DEALLOCATE (mat_work) 1942 END IF 1943 1944 CALL dbcsr_release(mat_L%matrix) 1945 DEALLOCATE (mat_L%matrix) 1946 IF (do_ri_Sigma_x .OR. do_ic_model) THEN 1947 CALL dbcsr_release(mat_SinvVSinv%matrix) 1948 DEALLOCATE (mat_SinvVSinv%matrix) 1949 END IF 1950 IF (do_ri_Sigma_x) THEN 1951 CALL dbcsr_release(mat_dm%matrix) 1952 DEALLOCATE (mat_dm%matrix) 1953 END IF 1954 1955 DEALLOCATE (index_to_cell_3c, cell_to_index_3c) 1956 1957 IF (.NOT. do_dbcsr_t) THEN 1958 CALL dbcsr_deallocate_matrix_set(mat_M_mu_Pnu_occ) 1959 CALL dbcsr_deallocate_matrix_set(mat_M_mu_Pnu_virt) 1960 ENDIF 1961 1962 CALL dbcsr_deallocate_matrix_set(mat_P_omega) 1963 IF (my_open_shell .AND. do_ri_sos_laplace_mp2) CALL dbcsr_deallocate_matrix_set(mat_P_omega_beta) 1964 1965 IF (.NOT. do_dbcsr_t) THEN 1966 CALL dbcsr_deallocate_matrix_set(mat_3c_overl_int_cut) 1967 1968 IF (do_mao) THEN 1969 CALL dbcsr_deallocate_matrix_set(mat_3c_overl_int_mao_for_occ_cut) 1970 CALL dbcsr_deallocate_matrix_set(mat_3c_overl_int_mao_for_virt_cut) 1971 END IF 1972 ELSE 1973 DO i_size = 1, SIZE(t_3c_O, 1) 1974 DO j_size = 1, SIZE(t_3c_O, 2) 1975 CALL dbcsr_t_destroy(t_3c_O(i_size, j_size)) 1976 CALL dbcsr_t_destroy(t_3c_overl_int(i_size, j_size)) 1977 ENDDO 1978 ENDDO 1979 1980 DEALLOCATE (t_3c_O, t_3c_overl_int) 1981 CALL dbcsr_t_destroy(t_3c_M) 1982 ENDIF 1983 1984 IF (.NOT. do_dbcsr_t) THEN 1985 DEALLOCATE (buffer_mat_M) 1986 DEALLOCATE (non_zero_blocks_3c) 1987 DEALLOCATE (non_zero_blocks_3c_cut_col) 1988 DEALLOCATE (cycle_due_to_sparse_dm, multiply_needed_occ, multiply_needed_virt) 1989 DEALLOCATE (row_from_LLL) 1990 DEALLOCATE (needed_cutRI_mem_R1vec_R2vec_for_kp) 1991 ENDIF 1992 DEALLOCATE (has_mat_P_blocks) 1993 1994 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 1995 CALL cp_cfm_release(cfm_mat_Q) 1996 CALL cp_fm_release(fm_mat_RI_global_work) 1997 CALL dbcsr_deallocate_matrix_set(mat_P_omega_kp) 1998 DEALLOCATE (wkp_W) 1999 END IF 2000 2001 CALL timestop(handle) 2002 2003 END SUBROUTINE dealloc_im_time 2004 2005! ************************************************************************************************** 2006!> \brief ... 2007!> \param mat_P_omega ... 2008!> \param mat_L ... 2009!> \param mat_work ... 2010!> \param eps_filter_im_time ... 2011!> \param fm_mat_work ... 2012!> \param dimen_RI ... 2013!> \param dimen_RI_red ... 2014!> \param fm_mat_L ... 2015!> \param fm_mat_Q ... 2016! ************************************************************************************************** 2017 SUBROUTINE contract_P_omega_with_mat_L(mat_P_omega, mat_L, mat_work, eps_filter_im_time, fm_mat_work, dimen_RI, & 2018 dimen_RI_red, fm_mat_L, fm_mat_Q) 2019 2020 TYPE(dbcsr_type), POINTER :: mat_P_omega, mat_L, mat_work 2021 REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time 2022 TYPE(cp_fm_type), POINTER :: fm_mat_work 2023 INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red 2024 TYPE(cp_fm_type), POINTER :: fm_mat_L, fm_mat_Q 2025 2026 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_P_omega_with_mat_L', & 2027 routineP = moduleN//':'//routineN 2028 2029 INTEGER :: handle 2030 2031 CALL timeset(routineN, handle) 2032 2033 ! multiplication with RI metric/Coulomb operator 2034 CALL dbcsr_multiply("N", "T", 1.0_dp, mat_P_omega, mat_L, & 2035 0.0_dp, mat_work, filter_eps=eps_filter_im_time) 2036 2037 CALL copy_dbcsr_to_fm(mat_work, fm_mat_work) 2038 2039 CALL cp_gemm('N', 'N', dimen_RI_red, dimen_RI_red, dimen_RI, 1.0_dp, fm_mat_L, fm_mat_work, & 2040 0.0_dp, fm_mat_Q) 2041 2042 ! Reset mat_work to save memory 2043 CALL dbcsr_set(mat_work, 0.0_dp) 2044 CALL dbcsr_filter(mat_work, 1.0_dp) 2045 2046 CALL timestop(handle) 2047 2048 END SUBROUTINE contract_P_omega_with_mat_L 2049 2050END MODULE rpa_util 2051