1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 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: dbcsr_create,& 44 dbcsr_filter,& 45 dbcsr_get_info,& 46 dbcsr_multiply,& 47 dbcsr_p_type,& 48 dbcsr_release,& 49 dbcsr_set,& 50 dbcsr_type 51 USE dbcsr_tensor_api, ONLY: dbcsr_t_destroy,& 52 dbcsr_t_type 53 USE input_constants, ONLY: wfc_mm_style_gemm,& 54 wfc_mm_style_syrk 55 USE kinds, ONLY: dp 56 USE kpoint_types, ONLY: get_kpoint_info,& 57 kpoint_type 58 USE machine, ONLY: m_walltime 59 USE mathconstants, ONLY: z_zero 60 USE message_passing, ONLY: mp_comm_split_direct,& 61 mp_sum 62 USE mp2_laplace, ONLY: calc_fm_mat_S_laplace 63 USE mp2_types, ONLY: integ_mat_buffer_type 64 USE qs_environment_types, ONLY: get_qs_env,& 65 qs_environment_type 66 USE rpa_communication, ONLY: fm_redistribute 67 USE rpa_gw_kpoints, ONLY: compute_wkp_W 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, & 77 dealloc_im_time, contract_P_omega_with_mat_L 78 79CONTAINS 80 81! ************************************************************************************************** 82!> \brief ... 83!> \param qs_env ... 84!> \param para_env ... 85!> \param dimen_RI ... 86!> \param dimen_RI_red ... 87!> \param num_integ_points ... 88!> \param fm_mat_Q ... 89!> \param fm_mo_coeff_occ ... 90!> \param fm_mo_coeff_virt ... 91!> \param fm_matrix_L_RI_metric ... 92!> \param mat_P_global ... 93!> \param t_3c_O ... 94!> \param matrix_s ... 95!> \param kpoints ... 96!> \param eps_filter_im_time ... 97!> \param do_ri_sos_laplace_mp2 ... 98!> \param cut_memory ... 99!> \param nkp ... 100!> \param num_cells_dm ... 101!> \param num_3c_repl ... 102!> \param size_P ... 103!> \param ikp_local ... 104!> \param index_to_cell_3c ... 105!> \param cell_to_index_3c ... 106!> \param col_blk_size ... 107!> \param do_ic_model ... 108!> \param do_kpoints_cubic_RPA ... 109!> \param do_kpoints_from_Gamma ... 110!> \param do_ri_Sigma_x ... 111!> \param my_open_shell ... 112!> \param has_mat_P_blocks ... 113!> \param wkp_W ... 114!> \param cfm_mat_Q ... 115!> \param fm_mat_L ... 116!> \param fm_mat_RI_global_work ... 117!> \param fm_mat_work ... 118!> \param fm_mo_coeff_occ_scaled ... 119!> \param fm_mo_coeff_virt_scaled ... 120!> \param mat_dm ... 121!> \param mat_L ... 122!> \param mat_M_P_munu_occ ... 123!> \param mat_M_P_munu_virt ... 124!> \param mat_P_global_copy ... 125!> \param mat_SinvVSinv ... 126!> \param mat_P_omega ... 127!> \param mat_P_omega_kp ... 128!> \param mat_work ... 129!> \param mat_P_omega_beta ... 130! ************************************************************************************************** 131 SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_points, & 132 fm_mat_Q, fm_mo_coeff_occ, fm_mo_coeff_virt, & 133 fm_matrix_L_RI_metric, mat_P_global, & 134 t_3c_O, matrix_s, kpoints, eps_filter_im_time, & 135 do_ri_sos_laplace_mp2, cut_memory, nkp, num_cells_dm, num_3c_repl, & 136 size_P, ikp_local, & 137 index_to_cell_3c, & 138 cell_to_index_3c, & 139 col_blk_size, & 140 do_ic_model, do_kpoints_cubic_RPA, & 141 do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, & 142 has_mat_P_blocks, wkp_W, & 143 cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, & 144 fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, mat_P_global_copy, & 145 mat_SinvVSinv, mat_P_omega, mat_P_omega_kp, & 146 mat_work, & 147 mat_P_omega_beta) 148 149 TYPE(qs_environment_type), POINTER :: qs_env 150 TYPE(cp_para_env_type), POINTER :: para_env 151 INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, num_integ_points 152 TYPE(cp_fm_type), POINTER :: fm_mat_Q, fm_mo_coeff_occ, & 153 fm_mo_coeff_virt 154 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_matrix_L_RI_metric 155 TYPE(dbcsr_p_type), INTENT(IN) :: mat_P_global 156 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :), & 157 INTENT(INOUT) :: t_3c_O 158 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 159 TYPE(kpoint_type), POINTER :: kpoints 160 REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time 161 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2 162 INTEGER, INTENT(IN) :: cut_memory 163 INTEGER, INTENT(OUT) :: nkp, num_cells_dm, num_3c_repl, size_P 164 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ikp_local 165 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell_3c 166 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & 167 INTENT(OUT) :: cell_to_index_3c 168 INTEGER, DIMENSION(:), POINTER :: col_blk_size 169 LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA, & 170 do_kpoints_from_Gamma, do_ri_Sigma_x, & 171 my_open_shell 172 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), & 173 INTENT(OUT) :: has_mat_P_blocks 174 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 175 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 176 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L 177 TYPE(cp_fm_type), POINTER :: fm_mat_RI_global_work, fm_mat_work, & 178 fm_mo_coeff_occ_scaled, & 179 fm_mo_coeff_virt_scaled 180 TYPE(dbcsr_p_type), INTENT(OUT) :: mat_dm, mat_L, mat_M_P_munu_occ, & 181 mat_M_P_munu_virt, mat_P_global_copy, & 182 mat_SinvVSinv 183 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega, mat_P_omega_kp 184 TYPE(dbcsr_type), POINTER :: mat_work 185 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_beta 186 187 CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_im_time', routineP = moduleN//':'//routineN 188 189 INTEGER :: cell_grid_dm(3), first_ikp_local, & 190 handle, i_dim, periodic(3) 191 INTEGER, DIMENSION(:), POINTER :: row_blk_size 192 TYPE(cell_type), POINTER :: cell 193 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_sub_kp 194 195 CALL timeset(routineN, handle) 196 197 num_3c_repl = SIZE(t_3c_O, 2) 198 199 IF (do_kpoints_cubic_RPA) THEN 200 ! we always use an odd number of image cells 201 ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical 202 DO i_dim = 1, 3 203 cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1 204 END DO 205 num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3) 206 ALLOCATE (index_to_cell_3c(3, SIZE(kpoints%index_to_cell, 2))) 207 CPASSERT(SIZE(kpoints%index_to_cell, 1) == 3) 208 index_to_cell_3c(:, :) = kpoints%index_to_cell(:, :) 209 ALLOCATE (cell_to_index_3c(LBOUND(kpoints%cell_to_index, 1):UBOUND(kpoints%cell_to_index, 1), & 210 LBOUND(kpoints%cell_to_index, 2):UBOUND(kpoints%cell_to_index, 2), & 211 LBOUND(kpoints%cell_to_index, 3):UBOUND(kpoints%cell_to_index, 3))) 212 cell_to_index_3c(:, :, :) = kpoints%cell_to_index(:, :, :) 213 214 ELSE 215 ALLOCATE (index_to_cell_3c(3, 1)) 216 index_to_cell_3c(:, 1) = 0 217 ALLOCATE (cell_to_index_3c(0:0, 0:0, 0:0)) 218 cell_to_index_3c(0, 0, 0) = 1 219 num_cells_dm = 1 220 END IF 221 222 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 223 224 CALL get_sub_para_kp(fm_struct_sub_kp, para_env, kpoints%nkp, & 225 dimen_RI, ikp_local, first_ikp_local, do_kpoints_cubic_RPA) 226 227 NULLIFY (cfm_mat_Q) 228 CALL cp_cfm_create(cfm_mat_Q, fm_struct_sub_kp) 229 CALL cp_cfm_set_all(cfm_mat_Q, z_zero) 230 ELSE 231 first_ikp_local = 1 232 END IF 233 234 ! if we do kpoints, mat_P has a kpoint and mat_P_omega has the inted 235 ! mat_P(tau, kpoint) 236 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 237 238 NULLIFY (cell) 239 CALL get_qs_env(qs_env, cell=cell) 240 CALL get_cell(cell=cell, periodic=periodic) 241 242 CALL get_kpoint_info(kpoints, nkp=nkp) 243 ! compute k-point weights such that functions 1/k^2, 1/k and const function are 244 ! integrated correctly 245 CALL compute_wkp_W(wkp_W, kpoints, cell%hmat, cell%h_inv, & 246 qs_env%mp2_env%ri_rpa_im_time%exp_kpoints, periodic) 247 ELSE 248 nkp = 1 249 END IF 250 251 IF (do_kpoints_cubic_RPA) THEN 252 size_P = MAX(num_cells_dm/2 + 1, nkp) 253 ELSE IF (do_kpoints_from_Gamma) THEN 254 size_P = MAX(3**(periodic(1) + periodic(2) + periodic(3)), nkp) 255 ELSE 256 size_P = 1 257 END IF 258 259 CALL alloc_mat_P_omega(mat_P_omega, num_integ_points, size_P, mat_P_global%matrix) 260 261 IF (my_open_shell .AND. do_ri_sos_laplace_mp2) THEN 262 CALL alloc_mat_P_omega(mat_P_omega_beta, num_integ_points, size_P, mat_P_global%matrix) 263 END IF 264 265 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 266 CALL alloc_mat_P_omega(mat_P_omega_kp, 2, size_P, mat_P_global%matrix) 267 END IF 268 269 NULLIFY (fm_mo_coeff_occ_scaled) 270 CALL cp_fm_create(fm_mo_coeff_occ_scaled, fm_mo_coeff_occ%matrix_struct) 271 CALL cp_fm_to_fm(fm_mo_coeff_occ, fm_mo_coeff_occ_scaled) 272 CALL cp_fm_set_all(matrix=fm_mo_coeff_occ_scaled, alpha=0.0_dp) 273 274 NULLIFY (fm_mo_coeff_virt_scaled) 275 CALL cp_fm_create(fm_mo_coeff_virt_scaled, fm_mo_coeff_virt%matrix_struct) 276 CALL cp_fm_to_fm(fm_mo_coeff_virt, fm_mo_coeff_virt_scaled) 277 CALL cp_fm_set_all(matrix=fm_mo_coeff_virt_scaled, alpha=0.0_dp) 278 279 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 280 NULLIFY (fm_mat_RI_global_work) 281 CALL cp_fm_create(fm_mat_RI_global_work, fm_matrix_L_RI_metric(1, 1)%matrix%matrix_struct) 282 CALL cp_fm_to_fm(fm_matrix_L_RI_metric(1, 1)%matrix, fm_mat_RI_global_work) 283 CALL cp_fm_set_all(fm_mat_RI_global_work, 0.0_dp) 284 END IF 285 286 ALLOCATE (has_mat_P_blocks(num_cells_dm/2 + 1, cut_memory, cut_memory, num_3c_repl, num_3c_repl)) 287 has_mat_P_blocks = .TRUE. 288 289 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 290 CALL reorder_mat_L(fm_mat_L, fm_matrix_L_RI_metric, fm_mat_Q%matrix_struct, para_env, mat_L, & 291 mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp) 292 293 CALL cp_fm_struct_release(fm_struct_sub_kp) 294 ELSE 295 CALL reorder_mat_L(fm_mat_L, fm_matrix_L_RI_metric, fm_mat_Q%matrix_struct, para_env, mat_L, & 296 mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local) 297 END IF 298 299 ! Create Scalapack working matrix for the contraction with the metric 300 IF (dimen_RI == dimen_RI_red) THEN 301 NULLIFY (fm_mat_work) 302 CALL cp_fm_create(fm_mat_work, fm_mat_Q%matrix_struct) 303 CALL cp_fm_to_fm(fm_mat_Q, fm_mat_work) 304 CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp) 305 306 ELSE 307 NULLIFY (fm_struct) 308 CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_mat_Q%matrix_struct, & 309 ncol_global=dimen_RI_red, nrow_global=dimen_RI) 310 311 NULLIFY (fm_mat_work) 312 CALL cp_fm_create(fm_mat_work, fm_struct) 313 CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp) 314 315 CALL cp_fm_struct_release(fm_struct) 316 317 END IF 318 319 ! Then its DBCSR counter part 320 IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN 321 CALL dbcsr_get_info(mat_L%matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size) 322 323 ! Create mat_work having the shape of the transposed of mat_L (compare with contract_P_omega_with_mat_L) 324 NULLIFY (mat_work) 325 ALLOCATE (mat_work) 326 CALL dbcsr_create(mat_work, template=mat_L%matrix, row_blk_size=col_blk_size, col_blk_size=row_blk_size) 327 END IF 328 329 IF (do_ri_Sigma_x .OR. do_ic_model) THEN 330 331 NULLIFY (mat_SinvVSinv%matrix) 332 ALLOCATE (mat_SinvVSinv%matrix) 333 CALL dbcsr_create(mat_SinvVSinv%matrix, template=mat_P_global%matrix) 334 CALL dbcsr_set(mat_SinvVSinv%matrix, 0.0_dp) 335 336 ! for kpoints we compute SinvVSinv later with kpoints 337 IF (.NOT. do_kpoints_from_Gamma) THEN 338 339 ! get the Coulomb matrix for Sigma_x = G*V 340 CALL dbcsr_multiply("T", "N", 1.0_dp, mat_L%matrix, mat_L%matrix, & 341 0.0_dp, mat_SinvVSinv%matrix, filter_eps=eps_filter_im_time) 342 343 END IF 344 345 END IF 346 347 IF (do_ri_Sigma_x) THEN 348 349 NULLIFY (mat_dm%matrix) 350 ALLOCATE (mat_dm%matrix) 351 CALL dbcsr_create(mat_dm%matrix, template=matrix_s(1)%matrix) 352 353 END IF 354 355 CALL timestop(handle) 356 357 END SUBROUTINE alloc_im_time 358 359! ************************************************************************************************** 360!> \brief ... 361!> \param mat_P_omega ... 362!> \param num_integ_points ... 363!> \param size_P ... 364!> \param template ... 365! ************************************************************************************************** 366 SUBROUTINE alloc_mat_P_omega(mat_P_omega, num_integ_points, size_P, template) 367 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega 368 INTEGER, INTENT(IN) :: num_integ_points, size_P 369 TYPE(dbcsr_type), POINTER :: template 370 371 CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_mat_P_omega', & 372 routineP = moduleN//':'//routineN 373 374 INTEGER :: handle, i_kp, jquad 375 376 CALL timeset(routineN, handle) 377 378 NULLIFY (mat_P_omega) 379 CALL dbcsr_allocate_matrix_set(mat_P_omega, num_integ_points, size_P) 380 DO jquad = 1, num_integ_points 381 DO i_kp = 1, size_P 382 ALLOCATE (mat_P_omega(jquad, i_kp)%matrix) 383 CALL dbcsr_create(matrix=mat_P_omega(jquad, i_kp)%matrix, & 384 template=template) 385 CALL dbcsr_set(mat_P_omega(jquad, i_kp)%matrix, 0.0_dp) 386 END DO 387 END DO 388 389 CALL timestop(handle) 390 391 END SUBROUTINE alloc_mat_P_omega 392 393! ************************************************************************************************** 394!> \brief ... 395!> \param fm_mat_L ... 396!> \param fm_matrix_L_RI_metric ... 397!> \param fm_struct_template ... 398!> \param para_env ... 399!> \param mat_L ... 400!> \param mat_template ... 401!> \param dimen_RI ... 402!> \param dimen_RI_red ... 403!> \param first_ikp_local ... 404!> \param ikp_local ... 405!> \param fm_struct_sub_kp ... 406! ************************************************************************************************** 407 SUBROUTINE reorder_mat_L(fm_mat_L, fm_matrix_L_RI_metric, fm_struct_template, para_env, mat_L, mat_template, & 408 dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp) 409 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L, fm_matrix_L_RI_metric 410 TYPE(cp_fm_struct_type), POINTER :: fm_struct_template 411 TYPE(cp_para_env_type), POINTER :: para_env 412 TYPE(dbcsr_p_type), INTENT(OUT) :: mat_L 413 TYPE(dbcsr_type), INTENT(IN) :: mat_template 414 INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, first_ikp_local 415 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: ikp_local 416 TYPE(cp_fm_struct_type), OPTIONAL, POINTER :: fm_struct_sub_kp 417 418 CHARACTER(LEN=*), PARAMETER :: routineN = 'reorder_mat_L', routineP = moduleN//':'//routineN 419 420 INTEGER :: handle, i_size, j_size, nblk 421 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size 422 LOGICAL :: do_kpoints 423 TYPE(cp_blacs_env_type), POINTER :: blacs_env 424 TYPE(cp_fm_struct_type), POINTER :: fm_struct 425 TYPE(cp_fm_type), POINTER :: fm_mat_L_transposed, fmdummy 426 427 CALL timeset(routineN, handle) 428 429 do_kpoints = .FALSE. 430 IF (PRESENT(ikp_local) .AND. PRESENT(fm_struct_sub_kp)) THEN 431 do_kpoints = .TRUE. 432 END IF 433 434 ! Get the fm_struct for fm_mat_L 435 NULLIFY (fm_struct) 436 IF (dimen_RI == dimen_RI_red) THEN 437 fm_struct => fm_struct_template 438 ELSE 439 ! The template is assumed to be square such that we need a new fm_struct if dimensions are not equal 440 CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI_red, ncol_global=dimen_RI, template_fmstruct=fm_struct_template) 441 END IF 442 443 ! Start to allocate the new full matrix 444 ALLOCATE (fm_mat_L(SIZE(fm_matrix_L_RI_metric, 1), SIZE(fm_matrix_L_RI_metric, 2))) 445 DO i_size = 1, SIZE(fm_matrix_L_RI_metric, 1) 446 DO j_size = 1, SIZE(fm_matrix_L_RI_metric, 2) 447 IF (do_kpoints) THEN 448 IF (ANY(ikp_local(:) == i_size)) THEN 449 CALL cp_fm_create(fm_mat_L(i_size, j_size)%matrix, fm_struct_sub_kp) 450 CALL cp_fm_set_all(fm_mat_L(i_size, j_size)%matrix, 0.0_dp) 451 END IF 452 ELSE 453 CALL cp_fm_create(fm_mat_L(i_size, j_size)%matrix, fm_struct) 454 CALL cp_fm_set_all(fm_mat_L(i_size, j_size)%matrix, 0.0_dp) 455 END IF 456 END DO 457 END DO 458 459 ! For the transposed matric we need a different fm_struct 460 IF (dimen_RI == dimen_RI_red) THEN 461 fm_struct => fm_mat_L(first_ikp_local, 1)%matrix%matrix_struct 462 ELSE 463 CALL cp_fm_struct_release(fm_struct) 464 465 ! Create a fm_struct with transposed sizes 466 NULLIFY (fm_struct) 467 CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI, ncol_global=dimen_RI_red, & 468 template_fmstruct=fm_mat_L(first_ikp_local, 1)%matrix%matrix_struct) !, force_block=.TRUE.) 469 END IF 470 471 ! Allocate buffer matrix 472 NULLIFY (fm_mat_L_transposed) 473 CALL cp_fm_create(fm_mat_L_transposed, fm_struct) 474 CALL cp_fm_set_all(matrix=fm_mat_L_transposed, alpha=0.0_dp) 475 476 IF (dimen_RI /= dimen_RI_red) CALL cp_fm_struct_release(fm_struct) 477 478 CALL cp_fm_get_info(fm_mat_L_transposed, context=blacs_env) 479 480 ! For k-points copy matrices of your group 481 ! Without kpoints, transpose matrix 482 ! without kpoints, the size of fm_mat_L is 1x1. with kpoints, the size is N_kpoints x 2 (2 for real/complex) 483 DO i_size = 1, SIZE(fm_matrix_L_RI_metric, 1) 484 DO j_size = 1, SIZE(fm_matrix_L_RI_metric, 2) 485 IF (do_kpoints) THEN 486 IF (ANY(ikp_local(:) == i_size)) THEN 487 CALL cp_fm_copy_general(fm_matrix_L_RI_metric(i_size, j_size)%matrix, fm_mat_L_transposed, para_env) 488 CALL cp_fm_to_fm(fm_mat_L_transposed, fm_mat_L(i_size, j_size)%matrix) 489 ELSE 490 NULLIFY (fmdummy) 491 CALL cp_fm_copy_general(fm_matrix_L_RI_metric(i_size, j_size)%matrix, fmdummy, para_env) 492 END IF 493 ELSE 494 CALL cp_fm_copy_general(fm_matrix_L_RI_metric(i_size, j_size)%matrix, fm_mat_L_transposed, blacs_env%para_env) 495 CALL cp_fm_transpose(fm_mat_L_transposed, fm_mat_L(i_size, j_size)%matrix) 496 END IF 497 END DO 498 END DO 499 500 ! Release old matrix 501 DO i_size = 1, SIZE(fm_matrix_L_RI_metric, 1) 502 DO j_size = 1, SIZE(fm_matrix_L_RI_metric, 2) 503 CALL cp_fm_release(fm_matrix_L_RI_metric(i_size, j_size)%matrix) 504 END DO 505 END DO 506 DEALLOCATE (fm_matrix_L_RI_metric) 507 ! Release buffer 508 CALL cp_fm_release(fm_mat_L_transposed) 509 510 ! Create sparse variant of L 511 NULLIFY (mat_L%matrix) 512 ALLOCATE (mat_L%matrix) 513 IF (dimen_RI == dimen_RI_red) THEN 514 CALL dbcsr_create(mat_L%matrix, template=mat_template) 515 ELSE 516 CALL dbcsr_get_info(mat_template, nblkrows_total=nblk, col_blk_size=col_blk_size) 517 518 CALL calculate_equal_blk_size(row_blk_size, dimen_RI_red, nblk) 519 520 CALL dbcsr_create(mat_L%matrix, template=mat_template, row_blk_size=row_blk_size, col_blk_size=col_blk_size) 521 522 DEALLOCATE (row_blk_size) 523 END IF 524 IF (.NOT. (do_kpoints)) THEN 525 CALL copy_fm_to_dbcsr(fm_mat_L(1, 1)%matrix, mat_L%matrix) 526 END IF 527 528 CALL timestop(handle) 529 530 END SUBROUTINE reorder_mat_L 531 532! ************************************************************************************************** 533!> \brief ... 534!> \param blk_size_new ... 535!> \param dimen_RI_red ... 536!> \param nblk ... 537! ************************************************************************************************** 538 SUBROUTINE calculate_equal_blk_size(blk_size_new, dimen_RI_red, nblk) 539 INTEGER, DIMENSION(:), POINTER :: blk_size_new 540 INTEGER, INTENT(IN) :: dimen_RI_red, nblk 541 542 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_equal_blk_size', & 543 routineP = moduleN//':'//routineN 544 545 INTEGER :: col_per_blk, remainder 546 547 NULLIFY (blk_size_new) 548 ALLOCATE (blk_size_new(nblk)) 549 550 remainder = MOD(dimen_RI_red, nblk) 551 col_per_blk = dimen_RI_red/nblk 552 553 ! Determine a new distribution for the columns (corresponding to the number of columns) 554 IF (remainder > 0) blk_size_new(1:remainder) = col_per_blk + 1 555 blk_size_new(remainder + 1:nblk) = col_per_blk 556 557 END SUBROUTINE calculate_equal_blk_size 558 559! ************************************************************************************************** 560!> \brief ... 561!> \param fm_mat_S ... 562!> \param do_ri_sos_laplace_mp2 ... 563!> \param first_cycle ... 564!> \param count_ev_sc_GW ... 565!> \param iter_sc_GW0 ... 566!> \param virtual ... 567!> \param Eigenval ... 568!> \param Eigenval_last ... 569!> \param Eigenval_scf ... 570!> \param homo ... 571!> \param omega ... 572!> \param omega_old ... 573!> \param jquad ... 574!> \param mm_style ... 575!> \param dimen_RI ... 576!> \param dimen_ia ... 577!> \param alpha ... 578!> \param fm_mat_Q ... 579!> \param fm_mat_Q_gemm ... 580!> \param para_env_RPA ... 581!> \param do_bse ... 582!> \param fm_mat_Q_static_bse_gemm ... 583!> \param RPA_proc_map ... 584!> \param buffer_rec ... 585!> \param buffer_send ... 586!> \param number_of_send ... 587!> \param map_send_size ... 588!> \param map_rec_size ... 589!> \param local_size_source ... 590!> \param my_num_dgemm_call ... 591!> \param my_flop_rate ... 592! ************************************************************************************************** 593 SUBROUTINE calc_mat_Q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, count_ev_sc_GW, iter_sc_GW0, virtual, & 594 Eigenval, Eigenval_last, Eigenval_scf, & 595 homo, omega, omega_old, jquad, mm_style, dimen_RI, dimen_ia, alpha, fm_mat_Q, fm_mat_Q_gemm, & 596 para_env_RPA, do_bse, fm_mat_Q_static_bse_gemm, RPA_proc_map, buffer_rec, buffer_send, number_of_send, & 597 map_send_size, map_rec_size, local_size_source, my_num_dgemm_call, my_flop_rate) 598 TYPE(cp_fm_type), POINTER :: fm_mat_S 599 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, first_cycle 600 INTEGER, INTENT(IN) :: count_ev_sc_GW, iter_sc_GW0, virtual 601 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval, Eigenval_last, Eigenval_scf 602 INTEGER, INTENT(IN) :: homo 603 REAL(KIND=dp), INTENT(IN) :: omega, omega_old 604 INTEGER, INTENT(IN) :: jquad, mm_style, dimen_RI, dimen_ia 605 REAL(KIND=dp), INTENT(IN) :: alpha 606 TYPE(cp_fm_type), POINTER :: fm_mat_Q, fm_mat_Q_gemm 607 TYPE(cp_para_env_type), POINTER :: para_env_RPA 608 LOGICAL, INTENT(IN) :: do_bse 609 TYPE(cp_fm_type), POINTER :: fm_mat_Q_static_bse_gemm 610 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: RPA_proc_map 611 TYPE(integ_mat_buffer_type), ALLOCATABLE, & 612 DIMENSION(:), INTENT(INOUT) :: buffer_rec, buffer_send 613 INTEGER, INTENT(IN) :: number_of_send 614 INTEGER, DIMENSION(0:para_env_RPA%num_pe-1), & 615 INTENT(IN) :: map_send_size, map_rec_size 616 INTEGER, DIMENSION(2, 0:para_env_RPA%num_pe-1), & 617 INTENT(IN) :: local_size_source 618 INTEGER, INTENT(INOUT) :: my_num_dgemm_call 619 REAL(KIND=dp), INTENT(INOUT) :: my_flop_rate 620 621 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mat_Q', routineP = moduleN//':'//routineN 622 623 INTEGER :: handle 624 625 CALL timeset(routineN, handle) 626 627 IF (do_ri_sos_laplace_mp2) THEN 628 ! the first index of tau_tj starts with 0 (see mp2_weights) 629 CALL calc_fm_mat_S_laplace(fm_mat_S, first_cycle, homo, virtual, Eigenval, omega, omega_old) 630 ELSE 631 IF (iter_sc_GW0 == 1) THEN 632 CALL calc_fm_mat_S_rpa(fm_mat_S, first_cycle, count_ev_sc_GW, virtual, Eigenval, Eigenval_last, & 633 homo, omega, omega_old) 634 ELSE 635 CALL calc_fm_mat_S_rpa(fm_mat_S, first_cycle, count_ev_sc_GW, virtual, Eigenval_scf, Eigenval_scf, & 636 homo, omega, omega_old) 637 END IF 638 END IF 639 640 CALL contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, para_env_RPA, & 641 my_num_dgemm_call, fm_mat_Q, RPA_proc_map, buffer_rec, buffer_send, number_of_send, & 642 map_send_size, map_rec_size, local_size_source, my_flop_rate) 643 644 IF (do_bse .AND. jquad == 1) THEN 645 CALL cp_fm_to_fm(fm_mat_Q_gemm, fm_mat_Q_static_bse_gemm) 646 END IF 647 CALL timestop(handle) 648 649 END SUBROUTINE calc_mat_Q 650 651! ************************************************************************************************** 652!> \brief ... 653!> \param fm_mat_S ... 654!> \param first_cycle ... 655!> \param count_ev_sc_GW ... 656!> \param virtual ... 657!> \param Eigenval ... 658!> \param Eigenval_last ... 659!> \param homo ... 660!> \param omega ... 661!> \param omega_old ... 662! ************************************************************************************************** 663 SUBROUTINE calc_fm_mat_S_rpa(fm_mat_S, first_cycle, count_ev_sc_GW, virtual, Eigenval, Eigenval_last, homo, & 664 omega, omega_old) 665 TYPE(cp_fm_type), POINTER :: fm_mat_S 666 LOGICAL, INTENT(IN) :: first_cycle 667 INTEGER, INTENT(IN) :: count_ev_sc_GW, virtual 668 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval, Eigenval_last 669 INTEGER, INTENT(IN) :: homo 670 REAL(KIND=dp), INTENT(IN) :: omega, omega_old 671 672 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_fm_mat_S_rpa', & 673 routineP = moduleN//':'//routineN 674 675 INTEGER :: avirt, handle, i_global, iiB, iocc, & 676 j_global, jjB, ncol_local, nrow_local 677 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 678 REAL(KIND=dp) :: eigen_diff 679 680 CALL timeset(routineN, handle) 681 682 ! get info of fm_mat_S 683 CALL cp_fm_get_info(matrix=fm_mat_S, & 684 nrow_local=nrow_local, & 685 ncol_local=ncol_local, & 686 row_indices=row_indices, & 687 col_indices=col_indices) 688 689 ! remove eigenvalue part of S matrix from the last eigenvalue self-c. cycle 690 IF (first_cycle .AND. count_ev_sc_GW > 1) THEN 691!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,iocc,avirt,eigen_diff,i_global,j_global) & 692!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Eigenval_last,fm_mat_S,virtual,homo,omega_old) 693 DO jjB = 1, ncol_local 694 j_global = col_indices(jjB) 695 DO iiB = 1, nrow_local 696 i_global = row_indices(iiB) 697 698 iocc = MAX(1, i_global - 1)/virtual + 1 699 avirt = i_global - (iocc - 1)*virtual 700 eigen_diff = Eigenval_last(avirt + homo) - Eigenval_last(iocc) 701 702 fm_mat_S%local_data(iiB, jjB) = fm_mat_S%local_data(iiB, jjB)/ & 703 SQRT(eigen_diff/(eigen_diff**2 + omega_old**2)) 704 705 END DO 706 END DO 707 708 END IF 709 710 ! update G matrix with the new value of omega 711 IF (first_cycle) THEN 712 ! In this case just update the matrix (symmetric form) with 713 ! SQRT((epsi_a-epsi_i)/((epsi_a-epsi_i)**2+omega**2)) 714!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,iocc,avirt,eigen_diff,i_global,j_global) & 715!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Eigenval,fm_mat_S,virtual,homo,omega) 716 DO jjB = 1, ncol_local 717 j_global = col_indices(jjB) 718 DO iiB = 1, nrow_local 719 i_global = row_indices(iiB) 720 721 iocc = MAX(1, i_global - 1)/virtual + 1 722 avirt = i_global - (iocc - 1)*virtual 723 eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc) 724 725 fm_mat_S%local_data(iiB, jjB) = fm_mat_S%local_data(iiB, jjB)* & 726 SQRT(eigen_diff/(eigen_diff**2 + omega**2)) 727 728 END DO 729 END DO 730 ELSE 731 ! In this case the update has to remove the old omega component thus 732 ! SQRT(((epsi_a-epsi_i)**2+omega_old**2)/((epsi_a-epsi_i)**2+omega**2)) 733!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,iocc,avirt,eigen_diff,i_global,j_global) & 734!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Eigenval,& 735!$OMP fm_mat_S,virtual,homo,omega,omega_old) 736 DO jjB = 1, ncol_local 737 j_global = col_indices(jjB) 738 DO iiB = 1, nrow_local 739 i_global = row_indices(iiB) 740 741 iocc = MAX(1, i_global - 1)/virtual + 1 742 avirt = i_global - (iocc - 1)*virtual 743 eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc) 744 745 fm_mat_S%local_data(iiB, jjB) = fm_mat_S%local_data(iiB, jjB)* & 746 SQRT((eigen_diff**2 + omega_old**2)/(eigen_diff**2 + omega**2)) 747 748 END DO 749 END DO 750 END IF 751 752 CALL timestop(handle) 753 754 END SUBROUTINE calc_fm_mat_S_rpa 755 756! ************************************************************************************************** 757!> \brief ... 758!> \param mm_style ... 759!> \param dimen_RI ... 760!> \param dimen_ia ... 761!> \param alpha ... 762!> \param fm_mat_S ... 763!> \param fm_mat_Q_gemm ... 764!> \param para_env_RPA ... 765!> \param my_num_dgemm_call ... 766!> \param fm_mat_Q ... 767!> \param RPA_proc_map ... 768!> \param buffer_rec ... 769!> \param buffer_send ... 770!> \param number_of_send ... 771!> \param map_send_size ... 772!> \param map_rec_size ... 773!> \param local_size_source ... 774!> \param my_flop_rate ... 775! ************************************************************************************************** 776 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, & 777 fm_mat_Q, RPA_proc_map, buffer_rec, buffer_send, number_of_send, & 778 map_send_size, map_rec_size, local_size_source, my_flop_rate) 779 780 INTEGER, INTENT(IN) :: mm_style, dimen_RI, dimen_ia 781 REAL(KIND=dp), INTENT(IN) :: alpha 782 TYPE(cp_fm_type), POINTER :: fm_mat_S, fm_mat_Q_gemm 783 TYPE(cp_para_env_type), POINTER :: para_env_RPA 784 INTEGER, INTENT(INOUT) :: my_num_dgemm_call 785 TYPE(cp_fm_type), POINTER :: fm_mat_Q 786 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: RPA_proc_map 787 TYPE(integ_mat_buffer_type), DIMENSION(:), & 788 INTENT(INOUT) :: buffer_rec, buffer_send 789 INTEGER, INTENT(IN) :: number_of_send 790 INTEGER, DIMENSION(0:para_env_RPA%num_pe-1), & 791 INTENT(IN) :: map_send_size, map_rec_size 792 INTEGER, DIMENSION(2, 0:para_env_RPA%num_pe-1), & 793 INTENT(IN) :: local_size_source 794 REAL(KIND=dp), INTENT(INOUT) :: my_flop_rate 795 796 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_S_to_Q', & 797 routineP = moduleN//':'//routineN 798 799 INTEGER :: handle 800 REAL(KIND=dp) :: actual_flop_rate, t_end, t_start 801 802 CALL timeset(routineN, handle) 803 804 t_start = m_walltime() 805 SELECT CASE (mm_style) 806 CASE (wfc_mm_style_gemm) 807 ! waste-fully computes the full symmetrix matrix, but maybe faster than cp_fm_syrk for optimized cp_fm_gemm 808 CALL cp_gemm(transa="T", transb="N", m=dimen_RI, n=dimen_RI, k=dimen_ia, alpha=alpha, & 809 matrix_a=fm_mat_S, matrix_b=fm_mat_S, beta=0.0_dp, & 810 matrix_c=fm_mat_Q_gemm) 811 CASE (wfc_mm_style_syrk) 812 ! will only compute the upper half of the matrix, which is fine, since we only use it for cholesky later 813 CALL cp_fm_syrk(uplo='U', trans='T', k=dimen_ia, alpha=alpha, matrix_a=fm_mat_S, & 814 ia=1, ja=1, beta=0.0_dp, matrix_c=fm_mat_Q_gemm) 815 CASE DEFAULT 816 CPABORT("") 817 END SELECT 818 t_end = m_walltime() 819 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)) 820 IF (para_env_RPA%mepos == 0) my_flop_rate = my_flop_rate + actual_flop_rate 821 my_num_dgemm_call = my_num_dgemm_call + 1 822 823 ! copy/redistribute fm_mat_Q_gemm to fm_mat_Q 824 CALL cp_fm_set_all(matrix=fm_mat_Q, alpha=0.0_dp) 825 CALL fm_redistribute(fm_mat_Q_gemm, fm_mat_Q, RPA_proc_map, buffer_rec, buffer_send, & 826 number_of_send, map_send_size, map_rec_size, local_size_source, para_env_RPA) 827 828 CALL timestop(handle) 829 830 END SUBROUTINE contract_S_to_Q 831 832! ************************************************************************************************** 833!> \brief ... 834!> \param dimen_RI ... 835!> \param trace_Qomega ... 836!> \param fm_mat_Q ... 837!> \param para_env_RPA ... 838! ************************************************************************************************** 839 SUBROUTINE RPA_postprocessing_start(dimen_RI, trace_Qomega, fm_mat_Q, para_env_RPA) 840 841 INTEGER, INTENT(IN) :: dimen_RI 842 REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(OUT) :: trace_Qomega 843 TYPE(cp_fm_type), POINTER :: fm_mat_Q 844 TYPE(cp_para_env_type), POINTER :: para_env_RPA 845 846 CHARACTER(LEN=*), PARAMETER :: routineN = 'RPA_postprocessing_start', & 847 routineP = moduleN//':'//routineN 848 849 INTEGER :: handle, i_global, iiB, j_global, jjB, & 850 ncol_local, nrow_local 851 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 852 853 CALL timeset(routineN, handle) 854 855 CALL cp_fm_get_info(matrix=fm_mat_Q, & 856 nrow_local=nrow_local, & 857 ncol_local=ncol_local, & 858 row_indices=row_indices, & 859 col_indices=col_indices) 860 861 ! calculate the trace of Q and add 1 on the diagonal 862 trace_Qomega = 0.0_dp 863!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 864!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,trace_Qomega,fm_mat_Q,dimen_RI) 865 DO jjB = 1, ncol_local 866 j_global = col_indices(jjB) 867 DO iiB = 1, nrow_local 868 i_global = row_indices(iiB) 869 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 870 trace_Qomega(i_global) = fm_mat_Q%local_data(iiB, jjB) 871 fm_mat_Q%local_data(iiB, jjB) = fm_mat_Q%local_data(iiB, jjB) + 1.0_dp 872 END IF 873 END DO 874 END DO 875 CALL mp_sum(trace_Qomega, para_env_RPA%group) 876 877 CALL timestop(handle) 878 879 END SUBROUTINE RPA_postprocessing_start 880 881! ************************************************************************************************** 882!> \brief ... 883!> \param dimen_RI ... 884!> \param trace_Qomega ... 885!> \param fm_mat_Q ... 886!> \param para_env_RPA ... 887!> \param Erpa ... 888!> \param wjquad ... 889! ************************************************************************************************** 890 SUBROUTINE RPA_postprocessing_nokp(dimen_RI, trace_Qomega, fm_mat_Q, para_env_RPA, Erpa, wjquad) 891 892 INTEGER, INTENT(IN) :: dimen_RI 893 REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(IN) :: trace_Qomega 894 TYPE(cp_fm_type), POINTER :: fm_mat_Q 895 TYPE(cp_para_env_type), POINTER :: para_env_RPA 896 REAL(KIND=dp), INTENT(INOUT) :: Erpa 897 REAL(KIND=dp), INTENT(IN) :: wjquad 898 899 CHARACTER(LEN=*), PARAMETER :: routineN = 'RPA_postprocessing_nokp', & 900 routineP = moduleN//':'//routineN 901 902 INTEGER :: handle, i_global, iiB, info_chol, & 903 j_global, jjB, ncol_local, nrow_local 904 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 905 REAL(KIND=dp) :: FComega 906 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Q_log 907 908 CALL timeset(routineN, handle) 909 910 CALL cp_fm_get_info(matrix=fm_mat_Q, & 911 nrow_local=nrow_local, & 912 ncol_local=ncol_local, & 913 row_indices=row_indices, & 914 col_indices=col_indices) 915 916 ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition 917 CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q, n=dimen_RI, info_out=info_chol) 918 IF (info_chol .NE. 0) THEN 919 CALL cp_warn(__LOCATION__, & 920 "The Cholesky decomposition before inverting the RPA matrix / dielectric "// & 921 "function failed. "// & 922 "In case of low-scaling RPA/GW, decreasing EPS_FILTER in the &LOW_SCALING "// & 923 "section might "// & 924 "increase the overall accuracy making the matrix positive definite. "// & 925 "Code will abort.") 926 END IF 927 928 CPASSERT(info_chol == 0) 929 930 ALLOCATE (Q_log(dimen_RI)) 931 Q_log = 0.0_dp 932!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 933!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Q_log,fm_mat_Q,dimen_RI) 934 DO jjB = 1, ncol_local 935 j_global = col_indices(jjB) 936 DO iiB = 1, nrow_local 937 i_global = row_indices(iiB) 938 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 939 Q_log(i_global) = 2.0_dp*LOG(fm_mat_Q%local_data(iiB, jjB)) 940 END IF 941 END DO 942 END DO 943 CALL mp_sum(Q_log, para_env_RPA%group) 944 945 FComega = 0.0_dp 946 DO iiB = 1, dimen_RI 947 IF (MODULO(iiB, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE 948 FComega = FComega + (Q_log(iiB) - trace_Qomega(iiB))/2.0_dp 949 END DO 950 Erpa = Erpa + FComega*wjquad 951 952 DEALLOCATE (Q_log) 953 954 CALL timestop(handle) 955 956 END SUBROUTINE RPA_postprocessing_nokp 957 958! ************************************************************************************************** 959!> \brief ... 960!> \param fm_struct_sub_kp ... 961!> \param para_env ... 962!> \param nkp ... 963!> \param dimen_RI ... 964!> \param ikp_local ... 965!> \param first_ikp_local ... 966!> \param do_kpoints_cubic_RPA ... 967! ************************************************************************************************** 968 SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, nkp, dimen_RI, & 969 ikp_local, first_ikp_local, do_kpoints_cubic_RPA) 970 TYPE(cp_fm_struct_type), POINTER :: fm_struct_sub_kp 971 TYPE(cp_para_env_type), POINTER :: para_env 972 INTEGER, INTENT(IN) :: nkp, dimen_RI 973 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ikp_local 974 INTEGER, INTENT(OUT) :: first_ikp_local 975 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA 976 977 CHARACTER(len=*), PARAMETER :: routineN = 'get_sub_para_kp', & 978 routineP = moduleN//':'//routineN 979 980 INTEGER :: color_sub_kp, comm_sub_kp, handle, ikp, & 981 num_proc_per_kp 982 TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub_kp 983 TYPE(cp_para_env_type), POINTER :: para_env_sub_kp 984 985 CALL timeset(routineN, handle) 986 987 IF (nkp > para_env%num_pe .OR. do_kpoints_cubic_RPA) THEN 988 ! we have all kpoints on every processpr 989 num_proc_per_kp = para_env%num_pe 990 ELSE 991 ! we have only one kpoint per group 992 num_proc_per_kp = para_env%num_pe/nkp 993 END IF 994 995 color_sub_kp = MOD(para_env%mepos/num_proc_per_kp, nkp) 996 CALL mp_comm_split_direct(para_env%group, comm_sub_kp, color_sub_kp) 997 NULLIFY (para_env_sub_kp) 998 CALL cp_para_env_create(para_env_sub_kp, comm_sub_kp) 999 NULLIFY (blacs_env_sub_kp) 1000 CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp) 1001 1002 NULLIFY (fm_struct_sub_kp) 1003 CALL cp_fm_struct_create(fm_struct_sub_kp, context=blacs_env_sub_kp, nrow_global=dimen_RI, & 1004 ncol_global=dimen_RI, para_env=para_env_sub_kp) 1005 1006 CALL cp_para_env_release(para_env_sub_kp) 1007 1008 CALL cp_blacs_env_release(blacs_env_sub_kp) 1009 1010 ALLOCATE (ikp_local(nkp)) 1011 ikp_local = 0 1012 first_ikp_local = 1 1013 DO ikp = 1, nkp 1014 IF (nkp > para_env%num_pe .OR. do_kpoints_cubic_RPA .OR. ikp == color_sub_kp + 1) THEN 1015 ikp_local(ikp) = ikp 1016 first_ikp_local = ikp 1017 END IF 1018 END DO 1019 1020 CALL timestop(handle) 1021 1022 END SUBROUTINE get_sub_para_kp 1023 1024! ************************************************************************************************** 1025!> \brief ... 1026!> \param do_ri_sos_laplace_mp2 ... 1027!> \param fm_mo_coeff_occ ... 1028!> \param fm_mo_coeff_virt ... 1029!> \param fm_scaled_dm_occ_tau ... 1030!> \param fm_scaled_dm_virt_tau ... 1031!> \param ikp_local ... 1032!> \param index_to_cell_3c ... 1033!> \param cell_to_index_3c ... 1034!> \param do_ic_model ... 1035!> \param do_kpoints_cubic_RPA ... 1036!> \param do_kpoints_from_Gamma ... 1037!> \param do_ri_Sigma_x ... 1038!> \param has_mat_P_blocks ... 1039!> \param wkp_W ... 1040!> \param cfm_mat_Q ... 1041!> \param fm_mat_L ... 1042!> \param fm_mat_RI_global_work ... 1043!> \param fm_mat_work ... 1044!> \param fm_mo_coeff_occ_scaled ... 1045!> \param fm_mo_coeff_virt_scaled ... 1046!> \param mat_dm ... 1047!> \param mat_L ... 1048!> \param mat_SinvVSinv ... 1049!> \param mat_P_omega ... 1050!> \param mat_P_omega_kp ... 1051!> \param t_3c_M ... 1052!> \param t_3c_O ... 1053!> \param mat_work ... 1054!> \param fm_mo_coeff_occ_beta ... 1055!> \param fm_mo_coeff_virt_beta ... 1056!> \param mat_P_omega_beta ... 1057! ************************************************************************************************** 1058 SUBROUTINE dealloc_im_time(do_ri_sos_laplace_mp2, fm_mo_coeff_occ, fm_mo_coeff_virt, & 1059 fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, ikp_local, index_to_cell_3c, & 1060 cell_to_index_3c, do_ic_model, & 1061 do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, & 1062 has_mat_P_blocks, & 1063 wkp_W, cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, & 1064 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_L, & 1065 mat_SinvVSinv, mat_P_omega, mat_P_omega_kp, & 1066 t_3c_M, t_3c_O, & 1067 mat_work, & 1068 fm_mo_coeff_occ_beta, fm_mo_coeff_virt_beta, mat_P_omega_beta) 1069 1070 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2 1071 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_occ, fm_mo_coeff_virt, & 1072 fm_scaled_dm_occ_tau, & 1073 fm_scaled_dm_virt_tau 1074 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: ikp_local 1075 INTEGER, ALLOCATABLE, DIMENSION(:, :), & 1076 INTENT(INOUT) :: index_to_cell_3c 1077 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & 1078 INTENT(INOUT) :: cell_to_index_3c 1079 LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA, & 1080 do_kpoints_from_Gamma, do_ri_Sigma_x 1081 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), & 1082 INTENT(INOUT) :: has_mat_P_blocks 1083 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 1084 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 1085 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L 1086 TYPE(cp_fm_type), POINTER :: fm_mat_RI_global_work, fm_mat_work, & 1087 fm_mo_coeff_occ_scaled, & 1088 fm_mo_coeff_virt_scaled 1089 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_dm, mat_L, mat_SinvVSinv 1090 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega, mat_P_omega_kp 1091 TYPE(dbcsr_t_type) :: t_3c_M 1092 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O 1093 TYPE(dbcsr_type), POINTER :: mat_work 1094 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_mo_coeff_occ_beta, & 1095 fm_mo_coeff_virt_beta 1096 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & 1097 POINTER :: mat_P_omega_beta 1098 1099 CHARACTER(LEN=*), PARAMETER :: routineN = 'dealloc_im_time', & 1100 routineP = moduleN//':'//routineN 1101 1102 INTEGER :: handle, i_size, j_size 1103 LOGICAL :: my_open_shell 1104 1105 CALL timeset(routineN, handle) 1106 1107 my_open_shell = .FALSE. 1108 IF (PRESENT(fm_mo_coeff_occ_beta) .AND. PRESENT(fm_mo_coeff_virt_beta) .AND. PRESENT(mat_P_omega_beta)) THEN 1109 my_open_shell = .TRUE. 1110 END IF 1111 1112 CALL cp_fm_release(fm_scaled_dm_occ_tau) 1113 CALL cp_fm_release(fm_scaled_dm_virt_tau) 1114 CALL cp_fm_release(fm_mo_coeff_occ) 1115 CALL cp_fm_release(fm_mo_coeff_virt) 1116 CALL cp_fm_release(fm_mo_coeff_occ_scaled) 1117 CALL cp_fm_release(fm_mo_coeff_virt_scaled) 1118 1119 IF (my_open_shell) THEN 1120 CALL cp_fm_release(fm_mo_coeff_occ_beta) 1121 CALL cp_fm_release(fm_mo_coeff_virt_beta) 1122 END IF 1123 1124 DO i_size = 1, SIZE(fm_mat_L, 1) 1125 DO j_size = 1, SIZE(fm_mat_L, 2) 1126 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 1127 IF (ANY(ikp_local(:) == i_size)) THEN 1128 CALL cp_fm_release(fm_mat_L(i_size, j_size)%matrix) 1129 END IF 1130 ELSE 1131 CALL cp_fm_release(fm_mat_L(i_size, j_size)%matrix) 1132 END IF 1133 END DO 1134 END DO 1135 DEALLOCATE (fm_mat_L) 1136 CALL cp_fm_release(fm_mat_work) 1137 1138 IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN 1139 CALL dbcsr_release(mat_work) 1140 DEALLOCATE (mat_work) 1141 END IF 1142 1143 CALL dbcsr_release(mat_L%matrix) 1144 DEALLOCATE (mat_L%matrix) 1145 IF (do_ri_Sigma_x .OR. do_ic_model) THEN 1146 CALL dbcsr_release(mat_SinvVSinv%matrix) 1147 DEALLOCATE (mat_SinvVSinv%matrix) 1148 END IF 1149 IF (do_ri_Sigma_x) THEN 1150 CALL dbcsr_release(mat_dm%matrix) 1151 DEALLOCATE (mat_dm%matrix) 1152 END IF 1153 1154 DEALLOCATE (index_to_cell_3c, cell_to_index_3c) 1155 1156 CALL dbcsr_deallocate_matrix_set(mat_P_omega) 1157 IF (my_open_shell .AND. do_ri_sos_laplace_mp2) CALL dbcsr_deallocate_matrix_set(mat_P_omega_beta) 1158 1159 DO i_size = 1, SIZE(t_3c_O, 1) 1160 DO j_size = 1, SIZE(t_3c_O, 2) 1161 CALL dbcsr_t_destroy(t_3c_O(i_size, j_size)) 1162 ENDDO 1163 ENDDO 1164 1165 DEALLOCATE (t_3c_O) 1166 CALL dbcsr_t_destroy(t_3c_M) 1167 1168 DEALLOCATE (has_mat_P_blocks) 1169 1170 IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN 1171 CALL cp_cfm_release(cfm_mat_Q) 1172 CALL cp_fm_release(fm_mat_RI_global_work) 1173 CALL dbcsr_deallocate_matrix_set(mat_P_omega_kp) 1174 DEALLOCATE (wkp_W) 1175 END IF 1176 1177 CALL timestop(handle) 1178 1179 END SUBROUTINE dealloc_im_time 1180 1181! ************************************************************************************************** 1182!> \brief ... 1183!> \param mat_P_omega ... 1184!> \param mat_L ... 1185!> \param mat_work ... 1186!> \param eps_filter_im_time ... 1187!> \param fm_mat_work ... 1188!> \param dimen_RI ... 1189!> \param dimen_RI_red ... 1190!> \param fm_mat_L ... 1191!> \param fm_mat_Q ... 1192! ************************************************************************************************** 1193 SUBROUTINE contract_P_omega_with_mat_L(mat_P_omega, mat_L, mat_work, eps_filter_im_time, fm_mat_work, dimen_RI, & 1194 dimen_RI_red, fm_mat_L, fm_mat_Q) 1195 1196 TYPE(dbcsr_type), POINTER :: mat_P_omega, mat_L, mat_work 1197 REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time 1198 TYPE(cp_fm_type), POINTER :: fm_mat_work 1199 INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red 1200 TYPE(cp_fm_type), POINTER :: fm_mat_L, fm_mat_Q 1201 1202 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_P_omega_with_mat_L', & 1203 routineP = moduleN//':'//routineN 1204 1205 INTEGER :: handle 1206 1207 CALL timeset(routineN, handle) 1208 1209 ! multiplication with RI metric/Coulomb operator 1210 CALL dbcsr_multiply("N", "T", 1.0_dp, mat_P_omega, mat_L, & 1211 0.0_dp, mat_work, filter_eps=eps_filter_im_time) 1212 1213 CALL copy_dbcsr_to_fm(mat_work, fm_mat_work) 1214 1215 CALL cp_gemm('N', 'N', dimen_RI_red, dimen_RI_red, dimen_RI, 1.0_dp, fm_mat_L, fm_mat_work, & 1216 0.0_dp, fm_mat_Q) 1217 1218 ! Reset mat_work to save memory 1219 CALL dbcsr_set(mat_work, 0.0_dp) 1220 CALL dbcsr_filter(mat_work, 1.0_dp) 1221 1222 CALL timestop(handle) 1223 1224 END SUBROUTINE contract_P_omega_with_mat_L 1225 1226END MODULE rpa_util 1227