1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for kpoint treatment in GW 8!> \par History 9!> 04.2019 created [Jan Wilhelm] 10! ************************************************************************************************** 11MODULE rpa_gw_kpoints 12 USE basis_set_types, ONLY: gto_basis_set_p_type 13 USE cell_types, ONLY: cell_type,& 14 get_cell,& 15 pbc 16 USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_invert,& 17 cp_cfm_gemm,& 18 cp_cfm_scale_and_add,& 19 cp_cfm_scale_and_add_fm,& 20 cp_cfm_transpose 21 USE cp_cfm_types, ONLY: cp_cfm_create,& 22 cp_cfm_get_info,& 23 cp_cfm_p_type,& 24 cp_cfm_release,& 25 cp_cfm_set_all,& 26 cp_cfm_to_fm,& 27 cp_cfm_type 28 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,& 29 dbcsr_allocate_matrix_set,& 30 dbcsr_deallocate_matrix_set 31 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add 32 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 33 cp_fm_struct_release,& 34 cp_fm_struct_type 35 USE cp_fm_types, ONLY: cp_fm_copy_general,& 36 cp_fm_create,& 37 cp_fm_p_type,& 38 cp_fm_release,& 39 cp_fm_set_all,& 40 cp_fm_set_element,& 41 cp_fm_to_fm,& 42 cp_fm_type 43 USE cp_gemm_interface, ONLY: cp_gemm 44 USE cp_para_types, ONLY: cp_para_env_type 45 USE dbcsr_api, ONLY: & 46 dbcsr_add, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_dot, dbcsr_filter, & 47 dbcsr_get_block_p, dbcsr_get_num_blocks, dbcsr_init_p, dbcsr_iterator_blocks_left, & 48 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 49 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_release_p, dbcsr_reserve_all_blocks, & 50 dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry 51 USE input_constants, ONLY: gw_read_exx 52 USE kinds, ONLY: dp,& 53 int_8 54 USE kpoint_types, ONLY: get_kpoint_info,& 55 kpoint_type 56 USE mathconstants, ONLY: gaussi,& 57 twopi,& 58 z_one,& 59 z_zero 60 USE mathlib, ONLY: invmat 61 USE message_passing, ONLY: mp_sum 62 USE particle_methods, ONLY: get_particle_set 63 USE particle_types, ONLY: particle_type 64 USE qs_environment_types, ONLY: get_qs_env,& 65 qs_environment_type 66 USE qs_integral_utils, ONLY: basis_set_list_setup 67 USE qs_kind_types, ONLY: qs_kind_type 68 USE rpa_gw_im_time_util, ONLY: fill_mat_3c_overl_int_gw,& 69 replicate_mat_to_subgroup_simple 70 USE rpa_im_time, ONLY: compute_transl_dm 71#include "./base/base_uses.f90" 72 73 IMPLICIT NONE 74 75 PRIVATE 76 77 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_kpoints' 78 79 PUBLIC :: compute_Wc_real_space_tau_GW, compute_Wc_kp_tau_GW, & 80 compute_wkp_W, compute_self_energy_im_time_gw_kp 81 82CONTAINS 83 84! ************************************************************************************************** 85!> \brief ... 86!> \param fm_mat_W_tau ... 87!> \param cfm_mat_Q ... 88!> \param fm_mat_L_re ... 89!> \param fm_mat_L_im ... 90!> \param dimen_RI ... 91!> \param num_integ_points ... 92!> \param jquad ... 93!> \param ikp ... 94!> \param tj ... 95!> \param tau_tj ... 96!> \param weights_cos_tf_w_to_t ... 97!> \param ikp_local ... 98!> \param para_env ... 99!> \param kpoints ... 100!> \param qs_env ... 101!> \param wkp_W ... 102!> \param mat_SinvVSinv ... 103!> \param do_W_and_not_V ... 104! ************************************************************************************************** 105 SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm_mat_L_im, & 106 dimen_RI, num_integ_points, jquad, & 107 ikp, tj, tau_tj, weights_cos_tf_w_to_t, ikp_local, & 108 para_env, kpoints, qs_env, wkp_W, mat_SinvVSinv, do_W_and_not_V) 109 110 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fm_mat_W_tau 111 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 112 TYPE(cp_fm_type), POINTER :: fm_mat_L_re, fm_mat_L_im 113 INTEGER :: dimen_RI, num_integ_points, jquad, ikp 114 REAL(KIND=dp), DIMENSION(:) :: tj 115 REAL(KIND=dp), DIMENSION(0:num_integ_points) :: tau_tj 116 REAL(KIND=dp), DIMENSION(:, :) :: weights_cos_tf_w_to_t 117 INTEGER, DIMENSION(:) :: ikp_local 118 TYPE(cp_para_env_type), POINTER :: para_env 119 TYPE(kpoint_type), POINTER :: kpoints 120 TYPE(qs_environment_type), POINTER :: qs_env 121 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 122 TYPE(dbcsr_p_type) :: mat_SinvVSinv 123 LOGICAL :: do_W_and_not_V 124 125 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Wc_real_space_tau_GW', & 126 routineP = moduleN//':'//routineN 127 128 INTEGER :: handle, handle2, i_global, iatom, iatom_old, icell, iiB, iquad, irow, j_global, & 129 jatom, jatom_old, jcol, jjB, jkp, LLL, natom, ncol_local, nkind, nkp, nrow_local, & 130 num_cells, xcell, ycell, zcell 131 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_RI_index 132 INTEGER, DIMENSION(:), POINTER :: col_indices, row_blk_end, row_blk_start, & 133 row_indices 134 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell 135 LOGICAL :: do_V_and_not_W 136 REAL(KIND=dp) :: abs_rab_cell, arg, contribution, coskl, cutoff_exp, d_0, omega, sinkl, & 137 sum_exp, sum_exp_k_im, sum_exp_k_re, tau, weight, weight_im, weight_re 138 REAL(KIND=dp), DIMENSION(3) :: cell_vector, rab_cell_i 139 REAL(KIND=dp), DIMENSION(3, 3) :: hmat 140 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp 141 REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp 142 TYPE(cell_type), POINTER :: cell 143 TYPE(cp_cfm_type), POINTER :: cfm_mat_L, cfm_mat_work, cfm_mat_work_2 144 TYPE(cp_fm_type), POINTER :: fm_dummy, fm_mat_work_global, & 145 fm_mat_work_local 146 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_RI_tmp 147 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 148 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 149 150 CALL timeset(routineN, handle) 151 152 CALL timeset(routineN//"_1", handle2) 153 154 NULLIFY (cfm_mat_work) 155 CALL cp_cfm_create(cfm_mat_work, cfm_mat_Q%matrix_struct) 156 CALL cp_cfm_set_all(cfm_mat_work, z_zero) 157 158 NULLIFY (cfm_mat_work_2) 159 CALL cp_cfm_create(cfm_mat_work_2, cfm_mat_Q%matrix_struct) 160 CALL cp_cfm_set_all(cfm_mat_work_2, z_zero) 161 162 NULLIFY (cfm_mat_L) 163 CALL cp_cfm_create(cfm_mat_L, cfm_mat_Q%matrix_struct) 164 CALL cp_cfm_set_all(cfm_mat_L, z_zero) 165 166 ! Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L 167 CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_L, z_one, fm_mat_L_re) 168 CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_L, gaussi, fm_mat_L_im) 169 170 NULLIFY (fm_mat_work_global) 171 CALL cp_fm_create(fm_mat_work_global, fm_mat_W_tau(1)%matrix%matrix_struct) 172 CALL cp_fm_set_all(fm_mat_work_global, 0.0_dp) 173 174 NULLIFY (fm_mat_work_local) 175 CALL cp_fm_create(fm_mat_work_local, cfm_mat_Q%matrix_struct) 176 CALL cp_fm_set_all(fm_mat_work_local, 0.0_dp) 177 178 CALL timestop(handle2) 179 180 IF (do_W_and_not_V) THEN 181 182 CALL timeset(routineN//"_2", handle2) 183 184 ! calculate [1+Q(iw')]^-1 185 CALL cp_cfm_cholesky_invert(cfm_mat_Q) 186 187 ! symmetrize the result 188 CALL own_cfm_upper_to_full(cfm_mat_Q, cfm_mat_work) 189 190 ! subtract exchange part by subtracing identity matrix from epsilon 191 CALL cp_cfm_get_info(matrix=cfm_mat_Q, & 192 nrow_local=nrow_local, & 193 ncol_local=ncol_local, & 194 row_indices=row_indices, & 195 col_indices=col_indices) 196 197 DO jjB = 1, ncol_local 198 j_global = col_indices(jjB) 199 DO iiB = 1, nrow_local 200 i_global = row_indices(iiB) 201 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 202 cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB) - z_one 203 END IF 204 END DO 205 END DO 206 207 CALL timestop(handle2) 208 209 CALL timeset(routineN//"_3.1", handle2) 210 211 ! work = epsilon(iw,k)*L^H(k) 212 CALL cp_cfm_gemm('N', 'C', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_Q, cfm_mat_L, & 213 z_zero, cfm_mat_work) 214 215 ! W(iw,k) = L(k)*work 216 CALL cp_cfm_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_work, & 217 z_zero, cfm_mat_work_2) 218 219 CALL timestop(handle2) 220 221 ELSE 222 223 ! S^-1(k)V(k)S^-1(k) = L(k)*L(k)^H 224 CALL cp_cfm_gemm('N', 'C', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_L, & 225 z_zero, cfm_mat_work_2) 226 227 END IF 228 229 CALL timeset(routineN//"_4", handle2) 230 231 CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp) 232 index_to_cell => kpoints%index_to_cell 233 num_cells = SIZE(index_to_cell, 2) 234 d_0 = qs_env%mp2_env%ri_rpa_im_time%cutoff 235 cutoff_exp = 10000.0_dp 236 CALL cp_cfm_set_all(cfm_mat_work, z_zero) 237 238 NULLIFY (qs_kind_set, cell, particle_set) 239 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell, natom=natom, nkind=nkind, & 240 particle_set=particle_set) 241 242 ALLOCATE (row_blk_start(natom)) 243 ALLOCATE (row_blk_end(natom)) 244 ALLOCATE (basis_set_RI_tmp(nkind)) 245 CALL basis_set_list_setup(basis_set_RI_tmp, "RI_AUX", qs_kind_set) 246 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, & 247 basis=basis_set_RI_tmp) 248 DEALLOCATE (basis_set_RI_tmp) 249 ALLOCATE (atom_from_RI_index(dimen_RI)) 250 DO LLL = 1, dimen_RI 251 DO iatom = 1, natom 252 IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN 253 atom_from_RI_index(LLL) = iatom 254 END IF 255 END DO 256 END DO 257 CALL get_cell(cell=cell, h=hmat) 258 iatom_old = 0 259 jatom_old = 0 260 261 CALL cp_cfm_get_info(matrix=cfm_mat_Q, & 262 nrow_local=nrow_local, & 263 ncol_local=ncol_local, & 264 row_indices=row_indices, & 265 col_indices=col_indices) 266 267 DO irow = 1, nrow_local 268 DO jcol = 1, ncol_local 269 270 iatom = atom_from_RI_index(row_indices(irow)) 271 jatom = atom_from_RI_index(col_indices(jcol)) 272 273 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN 274 275 sum_exp = 0.0_dp 276 sum_exp_k_re = 0.0_dp 277 sum_exp_k_im = 0.0_dp 278 279 DO icell = 1, num_cells 280 281 xcell = index_to_cell(1, icell) 282 ycell = index_to_cell(2, icell) 283 zcell = index_to_cell(3, icell) 284 285 arg = REAL(xcell, dp)*xkp(1, ikp) + REAL(ycell, dp)*xkp(2, ikp) + REAL(zcell, dp)*xkp(3, ikp) 286 287 coskl = wkp_W(ikp)*COS(twopi*arg) 288 sinkl = wkp_W(ikp)*SIN(twopi*arg) 289 290 cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, icell), dp)) 291 292 rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - & 293 (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3)) 294 295 abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2) 296 297 IF (abs_rab_cell/d_0 < cutoff_exp) THEN 298 sum_exp = sum_exp + EXP(-abs_rab_cell/d_0) 299 sum_exp_k_re = sum_exp_k_re + EXP(-abs_rab_cell/d_0)*coskl 300 sum_exp_k_im = sum_exp_k_im + EXP(-abs_rab_cell/d_0)*sinkl 301 END IF 302 303 END DO 304 305 weight_re = sum_exp_k_re/sum_exp 306 weight_im = sum_exp_k_im/sum_exp 307 308 iatom_old = iatom 309 jatom_old = jatom 310 311 END IF 312 313 contribution = weight_re*REAL(cfm_mat_work_2%local_data(irow, jcol)) + & 314 weight_im*AIMAG(cfm_mat_work_2%local_data(irow, jcol)) 315 316 fm_mat_work_local%local_data(irow, jcol) = fm_mat_work_local%local_data(irow, jcol) + contribution 317 318 END DO 319 END DO 320 321 CALL timestop(handle2) 322 323 CALL timeset(routineN//"_5", handle2) 324 325 IF (do_W_and_not_V) THEN 326 327 IF (SUM(ikp_local) > nkp) THEN 328 329 CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env) 330 331 DO iquad = 1, num_integ_points 332 333 omega = tj(jquad) 334 tau = tau_tj(iquad) 335 weight = weights_cos_tf_w_to_t(iquad, jquad)*COS(tau*omega) 336 337 IF (jquad == 1 .AND. ikp == 1) THEN 338 CALL cp_fm_set_all(matrix=fm_mat_W_tau(iquad)%matrix, alpha=0.0_dp) 339 END IF 340 341 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_W_tau(iquad)%matrix, beta=weight, matrix_b=fm_mat_work_global) 342 343 END DO 344 345 ELSE 346 347 DO jkp = 1, nkp 348 349 IF (ANY(ikp_local(:) == jkp)) THEN 350 CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env) 351 ELSE 352 NULLIFY (fm_dummy) 353 CALL cp_fm_copy_general(fm_dummy, fm_mat_work_global, para_env) 354 END IF 355 356 DO iquad = 1, num_integ_points 357 358 omega = tj(jquad) 359 tau = tau_tj(iquad) 360 weight = weights_cos_tf_w_to_t(iquad, jquad)*COS(tau*omega) 361 362 IF (jquad == 1 .AND. jkp == 1) THEN 363 CALL cp_fm_set_all(matrix=fm_mat_W_tau(iquad)%matrix, alpha=0.0_dp) 364 END IF 365 366 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_W_tau(iquad)%matrix, beta=weight, & 367 matrix_b=fm_mat_work_global) 368 369 END DO 370 371 END DO 372 373 END IF 374 375 END IF 376 377 do_V_and_not_W = .NOT. do_W_and_not_V 378 IF (do_V_and_not_W) THEN 379 380 IF (SUM(ikp_local) > nkp) THEN 381 CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env) 382 CALL fm_mat_work_global_to_mat_SinvVSinv(mat_SinvVSinv, fm_mat_work_global) 383 ELSE 384 DO jkp = 1, nkp 385 IF (ANY(ikp_local(:) == jkp)) THEN 386 CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env) 387 ELSE 388 NULLIFY (fm_dummy) 389 CALL cp_fm_copy_general(fm_dummy, fm_mat_work_global, para_env) 390 END IF 391 CALL fm_mat_work_global_to_mat_SinvVSinv(mat_SinvVSinv, fm_mat_work_global) 392 END DO 393 END IF 394 END IF 395 396 CALL cp_cfm_release(cfm_mat_work) 397 CALL cp_cfm_release(cfm_mat_work_2) 398 CALL cp_cfm_release(cfm_mat_L) 399 CALL cp_fm_release(fm_mat_work_global) 400 CALL cp_fm_release(fm_mat_work_local) 401 DEALLOCATE (atom_from_RI_index) 402 DEALLOCATE (row_blk_start) 403 DEALLOCATE (row_blk_end) 404 405 CALL timestop(handle2) 406 407 CALL timestop(handle) 408 409 END SUBROUTINE 410 411! ************************************************************************************************** 412!> \brief ... 413!> \param mat_SinvVSinv ... 414!> \param fm_mat_work_global ... 415! ************************************************************************************************** 416 SUBROUTINE fm_mat_work_global_to_mat_SinvVSinv(mat_SinvVSinv, fm_mat_work_global) 417 418 TYPE(dbcsr_p_type) :: mat_SinvVSinv 419 TYPE(cp_fm_type), POINTER :: fm_mat_work_global 420 421 CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_mat_work_global_to_mat_SinvVSinv', & 422 routineP = moduleN//':'//routineN 423 424 INTEGER :: handle 425 TYPE(dbcsr_p_type) :: mat_work 426 427 CALL timeset(routineN, handle) 428 429 NULLIFY (mat_work%matrix) 430 ALLOCATE (mat_work%matrix) 431 CALL dbcsr_create(mat_work%matrix, template=mat_SinvVSinv%matrix) 432 433 CALL copy_fm_to_dbcsr(fm_mat_work_global, mat_work%matrix, keep_sparsity=.FALSE.) 434 435 CALL dbcsr_add(mat_SinvVSinv%matrix, mat_work%matrix, 1.0_dp, 1.0_dp) 436 437 CALL dbcsr_release(mat_work%matrix) 438 DEALLOCATE (mat_work%matrix) 439 440 CALL timestop(handle) 441 442 END SUBROUTINE 443 444! ************************************************************************************************** 445!> \brief ... 446!> \param cfm_mat_W_kp_tau ... 447!> \param cfm_mat_Q ... 448!> \param fm_mat_L_re ... 449!> \param fm_mat_L_im ... 450!> \param dimen_RI ... 451!> \param num_integ_points ... 452!> \param jquad ... 453!> \param ikp ... 454!> \param tj ... 455!> \param tau_tj ... 456!> \param weights_cos_tf_w_to_t ... 457! ************************************************************************************************** 458 SUBROUTINE compute_Wc_kp_tau_GW(cfm_mat_W_kp_tau, cfm_mat_Q, fm_mat_L_re, fm_mat_L_im, & 459 dimen_RI, num_integ_points, jquad, & 460 ikp, tj, tau_tj, weights_cos_tf_w_to_t) 461 462 TYPE(cp_cfm_p_type), DIMENSION(:, :), POINTER :: cfm_mat_W_kp_tau 463 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 464 TYPE(cp_fm_type), POINTER :: fm_mat_L_re, fm_mat_L_im 465 INTEGER :: dimen_RI, num_integ_points, jquad, ikp 466 REAL(KIND=dp), DIMENSION(:) :: tj 467 REAL(KIND=dp), DIMENSION(0:num_integ_points) :: tau_tj 468 REAL(KIND=dp), DIMENSION(:, :) :: weights_cos_tf_w_to_t 469 470 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Wc_kp_tau_GW', & 471 routineP = moduleN//':'//routineN 472 473 INTEGER :: handle, handle2, i_global, iiB, iquad, & 474 j_global, jjB, ncol_local, nrow_local 475 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 476 REAL(KIND=dp) :: omega, tau, weight 477 TYPE(cp_cfm_type), POINTER :: cfm_mat_L, cfm_mat_work 478 479 CALL timeset(routineN, handle) 480 481 NULLIFY (cfm_mat_work) 482 CALL cp_cfm_create(cfm_mat_work, fm_mat_L_re%matrix_struct) 483 CALL cp_cfm_set_all(cfm_mat_work, z_zero) 484 485 NULLIFY (cfm_mat_L) 486 CALL cp_cfm_create(cfm_mat_L, fm_mat_L_re%matrix_struct) 487 CALL cp_cfm_set_all(cfm_mat_L, z_zero) 488 489 CALL timeset(routineN//"_cholesky_inv", handle2) 490 491 ! calculate [1+Q(iw')]^-1 492 CALL cp_cfm_cholesky_invert(cfm_mat_Q) 493 494 ! symmetrize the result 495 CALL own_cfm_upper_to_full(cfm_mat_Q, cfm_mat_work) 496 497 ! subtract exchange part by subtracing identity matrix from epsilon 498 CALL cp_cfm_get_info(matrix=cfm_mat_Q, & 499 nrow_local=nrow_local, & 500 ncol_local=ncol_local, & 501 row_indices=row_indices, & 502 col_indices=col_indices) 503 504 DO jjB = 1, ncol_local 505 j_global = col_indices(jjB) 506 DO iiB = 1, nrow_local 507 i_global = row_indices(iiB) 508 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 509 cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB) - z_one 510 END IF 511 END DO 512 END DO 513 514 CALL timestop(handle2) 515 516 ! Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L 517 CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_L, z_one, fm_mat_L_re) 518 CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_L, gaussi, fm_mat_L_im) 519 520 ! work = epsilon(iw,k)*L^H(k) 521 CALL cp_cfm_gemm('N', 'C', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_Q, cfm_mat_L, & 522 z_zero, cfm_mat_work) 523 524 ! W(iw,k) = L(k)*work 525 CALL cp_cfm_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_work, & 526 z_zero, cfm_mat_Q) 527 528 DO iquad = 1, num_integ_points 529 omega = tj(jquad) 530 tau = tau_tj(iquad) 531 weight = weights_cos_tf_w_to_t(iquad, jquad)*COS(tau*omega) 532 CALL cp_cfm_scale_and_add(alpha=z_one, matrix_a=cfm_mat_W_kp_tau(ikp, iquad)%matrix, & 533 beta=CMPLX(weight, KIND=dp), matrix_b=cfm_mat_Q) 534 END DO 535 536 CALL cp_cfm_release(cfm_mat_work) 537 CALL cp_cfm_release(cfm_mat_L) 538 539 CALL timestop(handle) 540 541 END SUBROUTINE 542 543! ************************************************************************************************** 544!> \brief ... 545!> \param wkp_W ... 546!> \param kpoints ... 547!> \param h_mat ... 548!> \param h_inv ... 549!> \param exp_kpoints ... 550!> \param periodic ... 551! ************************************************************************************************** 552 SUBROUTINE compute_wkp_W(wkp_W, kpoints, h_mat, h_inv, exp_kpoints, periodic) 553 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 554 TYPE(kpoint_type), POINTER :: kpoints 555 REAL(KIND=dp), DIMENSION(3, 3) :: h_mat, h_inv 556 REAL(KIND=dp) :: exp_kpoints 557 INTEGER, DIMENSION(3) :: periodic 558 559 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_wkp_W', routineP = moduleN//':'//routineN 560 561 INTEGER :: handle, i_dim, i_x, ikp, info, j_y, k_z, & 562 n_x, n_y, n_z, nkp, nsuperfine, & 563 num_lin_eqs 564 REAL(KIND=dp) :: a_vec_dot_k_vec, integral, k_sq, weight 565 REAL(KIND=dp), DIMENSION(3) :: a_vec, k_vec, x_vec 566 REAL(KIND=dp), DIMENSION(:), POINTER :: right_side, wkp 567 REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix_lin_eqs, xkp 568 569 CALL timeset(routineN, handle) 570 571 CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp) 572 573 ! we determine the kpoint weights of the Monkhors Pack mesh new 574 ! such that the functions 1/k^2, 1/k and const are integrated exactly 575 ! in the Brillouin zone 576 ! this is done by minimizing sum_i |w_i|^2 where w_i are the weights of 577 ! the i-th kpoint under the following constraints: 578 ! 1) 1/k^2, 1/k and const are integrated exactly 579 ! 2) the kpoint weights of kpoints with identical absolute value are 580 ! the same, of e.g. (1/8,3/8,3/8) same weight as for (3/8,1/8,3/8) 581 ! for 1d and 2d materials: we use normal Monkhorst-Pack mesh, checked 582 ! by SUM(periodic) == 3 583 584 IF (exp_kpoints < 2.0_dp .AND. SUM(periodic) == 3) THEN 585 586 ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid 587 nsuperfine = 500 588 integral = 0.0_dp 589 IF (exp_kpoints > 0.0_dp) exp_kpoints = -2.0_dp 590 591 ! actually, there is the factor *det_3x3(h_inv) missing to account for the 592 ! integration volume but for wkp det_3x3(h_inv) is needed 593 weight = 2.0_dp/(REAL(nsuperfine, dp))**3 594 DO i_x = 1, nsuperfine 595 DO j_y = 1, nsuperfine 596 DO k_z = 1, nsuperfine/2 597 598 x_vec = (/REAL(i_x - nsuperfine/2, dp) - 0.5_dp, & 599 REAL(j_y - nsuperfine/2, dp) - 0.5_dp, & 600 REAL(k_z - nsuperfine/2, dp) - 0.5_dp/)/ & 601 REAL(nsuperfine, dp) 602 k_vec = MATMUL(h_inv(1:3, 1:3), x_vec) 603 k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 604 integral = integral + weight*k_sq**(exp_kpoints*0.5_dp) 605 END DO 606 END DO 607 END DO 608 609 num_lin_eqs = nkp + 2 610 611 ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs)) 612 matrix_lin_eqs(:, :) = 0.0_dp 613 614 DO ikp = 1, nkp 615 616 k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp)) 617 k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 618 619 matrix_lin_eqs(ikp, ikp) = 2.0_dp 620 matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp 621 matrix_lin_eqs(ikp, nkp + 2) = 1.0_dp*k_sq**(exp_kpoints*0.5_dp) 622 623 matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp 624 matrix_lin_eqs(nkp + 2, ikp) = 1.0_dp*k_sq**(exp_kpoints*0.5_dp) 625 626 END DO 627 628 CALL invmat(matrix_lin_eqs, info) 629 ! check whether inversion was successfull 630 CPASSERT(info == 0) 631 632 ALLOCATE (wkp_W(num_lin_eqs)) 633 634 ALLOCATE (right_side(num_lin_eqs)) 635 right_side = 0.0_dp 636 right_side(nkp + 1) = 1.0_dp 637 right_side(nkp + 2) = integral 638 639 wkp_W(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side) 640 641 DEALLOCATE (matrix_lin_eqs, right_side) 642 643 ELSE IF (exp_kpoints < 2.0_dp .AND. SUM(periodic) == 1) THEN 644 645 ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid 646 nsuperfine = 5000 647 integral = 0.0_dp 648 649 ! actually, there is the factor *det_3x3(h_inv) missing to account for the 650 ! integration volume but for wkp det_3x3(h_inv) is needed 651 weight = 1.0_dp/REAL(nsuperfine, dp) 652 IF (periodic(1) == 1) THEN 653 n_x = nsuperfine 654 ELSE 655 n_x = 1 656 END IF 657 IF (periodic(2) == 1) THEN 658 n_y = nsuperfine 659 ELSE 660 n_y = 1 661 END IF 662 IF (periodic(3) == 1) THEN 663 n_z = nsuperfine 664 ELSE 665 n_z = 1 666 END IF 667 668 a_vec = MATMUL(h_mat(1:3, 1:3), & 669 (/REAL(periodic(1), dp), REAL(periodic(2), dp), REAL(periodic(3), dp)/)) 670 671 DO i_x = 1, n_x 672 DO j_y = 1, n_y 673 DO k_z = 1, n_z 674 675 x_vec = (/REAL(i_x - nsuperfine/2, dp) - 0.5_dp, & 676 REAL(j_y - nsuperfine/2, dp) - 0.5_dp, & 677 REAL(k_z - nsuperfine/2, dp) - 0.5_dp/)/ & 678 REAL(nsuperfine, dp) 679 680 DO i_dim = 1, 3 681 IF (periodic(i_dim) == 0) THEN 682 x_vec(i_dim) = 0.0_dp 683 END IF 684 END DO 685 686 k_vec = MATMUL(h_inv(1:3, 1:3), x_vec) 687 a_vec_dot_k_vec = a_vec(1)*k_vec(1) + a_vec(2)*k_vec(2) + a_vec(3)*k_vec(3) 688 integral = integral + weight*LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) 689 END DO 690 END DO 691 END DO 692 693 num_lin_eqs = nkp + 2 694 695 ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs)) 696 matrix_lin_eqs(:, :) = 0.0_dp 697 698 DO ikp = 1, nkp 699 700 k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp)) 701 k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 702 703 matrix_lin_eqs(ikp, ikp) = 2.0_dp 704 matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp 705 706 a_vec_dot_k_vec = a_vec(1)*k_vec(1) + a_vec(2)*k_vec(2) + a_vec(3)*k_vec(3) 707 matrix_lin_eqs(ikp, nkp + 2) = LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) 708 709 matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp 710 matrix_lin_eqs(nkp + 2, ikp) = LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) 711 712 END DO 713 714 CALL invmat(matrix_lin_eqs, info) 715 ! check whether inversion was successfull 716 CPASSERT(info == 0) 717 718 ALLOCATE (wkp_W(num_lin_eqs)) 719 720 ALLOCATE (right_side(num_lin_eqs)) 721 right_side = 0.0_dp 722 right_side(nkp + 1) = 1.0_dp 723 right_side(nkp + 2) = integral 724 725 wkp_W(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side) 726 727 DEALLOCATE (matrix_lin_eqs, right_side) 728 729 ELSE 730 731 ALLOCATE (wkp_W(nkp)) 732 wkp_W(:) = wkp(:) 733 734 END IF 735 736 CALL timestop(handle) 737 738 END SUBROUTINE 739 740! ************************************************************************************************** 741!> \brief ... 742!> \param cfm_mat_Q ... 743!> \param cfm_mat_work ... 744! ************************************************************************************************** 745 SUBROUTINE own_cfm_upper_to_full(cfm_mat_Q, cfm_mat_work) 746 747 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q, cfm_mat_work 748 749 CHARACTER(LEN=*), PARAMETER :: routineN = 'own_cfm_upper_to_full', & 750 routineP = moduleN//':'//routineN 751 752 INTEGER :: handle, i_global, iiB, j_global, jjB, & 753 ncol_local, nrow_local 754 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 755 756 CALL timeset(routineN, handle) 757 758 ! get info of fm_mat_Q 759 CALL cp_cfm_get_info(matrix=cfm_mat_Q, & 760 nrow_local=nrow_local, & 761 ncol_local=ncol_local, & 762 row_indices=row_indices, & 763 col_indices=col_indices) 764 765 DO jjB = 1, ncol_local 766 j_global = col_indices(jjB) 767 DO iiB = 1, nrow_local 768 i_global = row_indices(iiB) 769 IF (j_global < i_global) THEN 770 cfm_mat_Q%local_data(iiB, jjB) = z_zero 771 END IF 772 IF (j_global == i_global) THEN 773 cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB)/(2.0_dp, 0.0_dp) 774 END IF 775 END DO 776 END DO 777 778 CALL cp_cfm_transpose(cfm_mat_Q, 'C', cfm_mat_work) 779 780 CALL cp_cfm_scale_and_add(z_one, cfm_mat_Q, z_one, cfm_mat_work) 781 782 CALL timestop(handle) 783 784 END SUBROUTINE 785 786! ************************************************************************************************** 787!> \brief ... 788!> \param vec_Sigma_c_gw ... 789!> \param vec_Sigma_x_gw ... 790!> \param vec_Sigma_x_minus_vxc_gw ... 791!> \param mat_3c_overl_int ... 792!> \param cell_to_index_3c ... 793!> \param index_to_cell_3c ... 794!> \param num_cells_dm ... 795!> \param kpoints ... 796!> \param unit_nr ... 797!> \param gw_corr_lev_tot ... 798!> \param num_3c_repl ... 799!> \param nkp_self_energy ... 800!> \param num_fit_points ... 801!> \param RI_blk_sizes ... 802!> \param prim_blk_sizes ... 803!> \param matrix_s ... 804!> \param para_env ... 805!> \param para_env_sub ... 806!> \param gw_corr_lev_occ ... 807!> \param gw_corr_lev_virt ... 808!> \param dimen_RI ... 809!> \param homo ... 810!> \param nmo ... 811!> \param cut_RI ... 812!> \param mat_dm_virt_local ... 813!> \param row_from_LLL ... 814!> \param my_group_L_starts_im_time ... 815!> \param my_group_L_sizes_im_time ... 816!> \param cfm_mat_Q ... 817!> \param cfm_mat_W_kp_tau ... 818!> \param qs_env ... 819!> \param e_fermi ... 820!> \param eps_filter ... 821!> \param tj ... 822!> \param tau_tj ... 823!> \param weights_sin_tf_t_to_w ... 824!> \param weights_cos_tf_t_to_w ... 825!> \param num_integ_points ... 826!> \param stabilize_exp ... 827!> \param fm_mat_L ... 828!> \param wkp_W ... 829! ************************************************************************************************** 830 SUBROUTINE compute_self_energy_im_time_gw_kp(vec_Sigma_c_gw, vec_Sigma_x_gw, vec_Sigma_x_minus_vxc_gw, & 831 mat_3c_overl_int, cell_to_index_3c, index_to_cell_3c, & 832 num_cells_dm, kpoints, & 833 unit_nr, gw_corr_lev_tot, num_3c_repl, nkp_self_energy, num_fit_points, & 834 RI_blk_sizes, prim_blk_sizes, matrix_s, para_env, para_env_sub, & 835 gw_corr_lev_occ, gw_corr_lev_virt, dimen_RI, homo, nmo, cut_RI, mat_dm_virt_local, & 836 row_from_LLL, my_group_L_starts_im_time, my_group_L_sizes_im_time, & 837 cfm_mat_Q, cfm_mat_W_kp_tau, qs_env, e_fermi, eps_filter, tj, tau_tj, & 838 weights_sin_tf_t_to_w, weights_cos_tf_t_to_w, & 839 num_integ_points, stabilize_exp, fm_mat_L, wkp_W) 840 841 COMPLEX(KIND=dp), DIMENSION(:, :, :) :: vec_Sigma_c_gw 842 REAL(KIND=dp), DIMENSION(:, :) :: vec_Sigma_x_gw, vec_Sigma_x_minus_vxc_gw 843 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int 844 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c 845 INTEGER, DIMENSION(:, :) :: index_to_cell_3c 846 INTEGER :: num_cells_dm 847 TYPE(kpoint_type), POINTER :: kpoints 848 INTEGER :: unit_nr, gw_corr_lev_tot, num_3c_repl, & 849 nkp_self_energy, num_fit_points 850 INTEGER, DIMENSION(:), POINTER :: RI_blk_sizes, prim_blk_sizes 851 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 852 TYPE(cp_para_env_type), POINTER :: para_env, para_env_sub 853 INTEGER :: gw_corr_lev_occ, gw_corr_lev_virt, & 854 dimen_RI, homo, nmo, cut_RI 855 TYPE(dbcsr_p_type) :: mat_dm_virt_local 856 INTEGER, DIMENSION(:) :: row_from_LLL, my_group_L_starts_im_time, & 857 my_group_L_sizes_im_time 858 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 859 TYPE(cp_cfm_p_type), DIMENSION(:, :), POINTER :: cfm_mat_W_kp_tau 860 TYPE(qs_environment_type), POINTER :: qs_env 861 REAL(KIND=dp) :: e_fermi, eps_filter 862 REAL(KIND=dp), DIMENSION(:) :: tj 863 INTEGER :: num_integ_points 864 REAL(KIND=dp), DIMENSION(:, :) :: weights_cos_tf_t_to_w, & 865 weights_sin_tf_t_to_w 866 REAL(KIND=dp), DIMENSION(0:num_integ_points) :: tau_tj 867 REAL(KIND=dp) :: stabilize_exp 868 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L 869 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 870 871 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_self_energy_im_time_gw_kp', & 872 routineP = moduleN//':'//routineN 873 874 INTEGER :: handle, ikp, maxcell, & 875 num_cells_R1_plus_S2, num_cells_R2, & 876 start_jquad 877 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_R1, & 878 index_to_cell_R1_plus_S2, & 879 index_to_cell_R2 880 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_R2 881 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 882 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: has_blocks_mat_dm_occ_S, & 883 has_blocks_mat_dm_virt_S 884 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: cycle_R1_plus_S2_n_level, & 885 has_3c_blocks_im, has_3c_blocks_re 886 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :) :: cycle_R1_S2_n_level 887 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: are_flops_I_T_R1_plus_S2_S1_n_level 888 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_I_muP_occ_im, mat_I_muP_occ_re, & 889 mat_I_muP_virt_im, mat_I_muP_virt_re 890 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_S, mat_dm_virt_S, mat_W_R 891 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_gw_kp_im, & 892 mat_3c_overl_int_gw_kp_re 893 TYPE(dbcsr_type), POINTER :: mat_contr_gf_occ_im, mat_contr_gf_occ_re, mat_contr_gf_virt_im, & 894 mat_contr_gf_virt_re, mat_contr_W_im, mat_contr_W_re 895 896 CALL timeset(routineN, handle) 897 898 ! R2 is index on W^R2 899 CALL compute_cell_vec_for_R2(index_to_cell_R2, cell_to_index_R2, num_cells_R2, unit_nr, & 900 index_to_cell_dm, qs_env) 901 902 ! R1+S2 is index on 3c integral 903 CALL compute_cell_vec_for_R1_plus_S2_or_R1(index_to_cell_R1_plus_S2, cell_to_index_3c, & 904 num_cells_R1_plus_S2, kpoints, & 905 unit_nr, maxcell, num_cells_dm) 906 907 CALL allocate_mat_3c_gw(mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 908 gw_corr_lev_tot, num_3c_repl, num_cells_dm, num_cells_R1_plus_S2, num_integ_points, & 909 dimen_RI, nmo, RI_blk_sizes, prim_blk_sizes, matrix_s, cfm_mat_Q, & 910 mat_contr_gf_occ_re, mat_contr_gf_occ_im, mat_contr_gf_virt_re, mat_contr_gf_virt_im, & 911 mat_contr_W_re, mat_contr_W_im, mat_I_muP_occ_re, mat_I_muP_virt_re, & 912 mat_I_muP_occ_im, mat_I_muP_virt_im, has_3c_blocks_re, has_3c_blocks_im, & 913 cycle_R1_S2_n_level, & 914 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, cycle_R1_plus_S2_n_level, & 915 are_flops_I_T_R1_plus_S2_S1_n_level) 916 917 ! get W^R2 918 CALL trafo_W_from_k_to_R(index_to_cell_R2, mat_W_R, mat_3c_overl_int_gw_kp_re, cfm_mat_W_kp_tau, & 919 kpoints, RI_blk_sizes, fm_mat_L, dimen_RI, qs_env, & 920 start_jquad, wkp_W) 921 ! get G^S1 922 CALL compute_G_real_space(mat_dm_occ_S, mat_dm_virt_S, qs_env, num_integ_points, stabilize_exp, & 923 e_fermi, eps_filter, tau_tj, num_cells_dm, index_to_cell_dm, & 924 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, para_env) 925 926 DO ikp = 1, nkp_self_energy 927 928 CALL get_mat_3c_gw_kp(mat_3c_overl_int, ikp, & 929 mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 930 kpoints, para_env, para_env_sub, matrix_s, mat_dm_virt_local, & 931 gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, cut_RI, num_3c_repl, & 932 row_from_LLL, my_group_L_starts_im_time, & 933 my_group_L_sizes_im_time, has_3c_blocks_re, has_3c_blocks_im) 934 935 CALL cell_sum_self_ener(vec_Sigma_c_gw(:, :, ikp), vec_Sigma_x_gw(:, ikp), & 936 vec_Sigma_x_minus_vxc_gw(:, ikp), & 937 mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 938 mat_dm_occ_S, mat_dm_virt_S, mat_W_R, & 939 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 940 mat_contr_gf_virt_re, mat_contr_gf_virt_im, mat_contr_W_re, mat_contr_W_im, & 941 mat_I_muP_occ_re, mat_I_muP_virt_re, mat_I_muP_occ_im, mat_I_muP_virt_im, & 942 index_to_cell_R1, index_to_cell_R2, index_to_cell_3c, & 943 index_to_cell_dm, index_to_cell_R1_plus_S2, cell_to_index_3c, & 944 cell_to_index_R2, gw_corr_lev_occ, gw_corr_lev_tot, & 945 homo, num_integ_points, num_fit_points, tj, tau_tj, ikp, & 946 has_3c_blocks_re, has_3c_blocks_im, & 947 cycle_R1_S2_n_level, has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, & 948 cycle_R1_plus_S2_n_level, kpoints, & 949 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, eps_filter, para_env, & 950 are_flops_I_T_R1_plus_S2_S1_n_level, start_jquad) 951 952 END DO 953 954 CALL clean_up_self_energy_kp(mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 955 index_to_cell_R2, index_to_cell_R1_plus_S2, & 956 mat_W_R, mat_dm_occ_S, mat_dm_virt_S, & 957 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 958 mat_contr_gf_virt_re, mat_contr_gf_virt_im, mat_contr_W_re, mat_contr_W_im, & 959 mat_I_muP_occ_re, mat_I_muP_virt_re, mat_I_muP_occ_im, mat_I_muP_virt_im, & 960 has_3c_blocks_re, has_3c_blocks_im, cycle_R1_S2_n_level, & 961 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, cycle_R1_plus_S2_n_level, & 962 are_flops_I_T_R1_plus_S2_S1_n_level) 963 964 CALL timestop(handle) 965 966 END SUBROUTINE 967 968! ************************************************************************************************** 969!> \brief ... 970!> \param vec_Sigma_c_gw ... 971!> \param vec_Sigma_x_gw ... 972!> \param vec_Sigma_x_minus_vxc_gw ... 973!> \param mat_3c_overl_int_gw_kp_re ... 974!> \param mat_3c_overl_int_gw_kp_im ... 975!> \param mat_dm_occ_S ... 976!> \param mat_dm_virt_S ... 977!> \param mat_W_R ... 978!> \param mat_contr_gf_occ_re ... 979!> \param mat_contr_gf_occ_im ... 980!> \param mat_contr_gf_virt_re ... 981!> \param mat_contr_gf_virt_im ... 982!> \param mat_contr_W_re ... 983!> \param mat_contr_W_im ... 984!> \param mat_I_muP_occ_re ... 985!> \param mat_I_muP_virt_re ... 986!> \param mat_I_muP_occ_im ... 987!> \param mat_I_muP_virt_im ... 988!> \param index_to_cell_R1 ... 989!> \param index_to_cell_R2 ... 990!> \param index_to_cell_3c ... 991!> \param index_to_cell_dm ... 992!> \param index_to_cell_R1_plus_S2 ... 993!> \param cell_to_index_3c ... 994!> \param cell_to_index_R2 ... 995!> \param gw_corr_lev_occ ... 996!> \param gw_corr_lev_tot ... 997!> \param homo ... 998!> \param num_integ_points ... 999!> \param num_fit_points ... 1000!> \param tj ... 1001!> \param tau_tj ... 1002!> \param ikp ... 1003!> \param has_3c_blocks_re ... 1004!> \param has_3c_blocks_im ... 1005!> \param cycle_R1_S2_n_level ... 1006!> \param has_blocks_mat_dm_occ_S ... 1007!> \param has_blocks_mat_dm_virt_S ... 1008!> \param cycle_R1_plus_S2_n_level ... 1009!> \param kpoints ... 1010!> \param weights_cos_tf_t_to_w ... 1011!> \param weights_sin_tf_t_to_w ... 1012!> \param eps_filter ... 1013!> \param para_env ... 1014!> \param are_flops_I_T_R1_plus_S2_S1_n_level ... 1015!> \param start_jquad ... 1016! ************************************************************************************************** 1017 SUBROUTINE cell_sum_self_ener(vec_Sigma_c_gw, vec_Sigma_x_gw, vec_Sigma_x_minus_vxc_gw, & 1018 mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 1019 mat_dm_occ_S, mat_dm_virt_S, mat_W_R, & 1020 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 1021 mat_contr_gf_virt_re, mat_contr_gf_virt_im, mat_contr_W_re, mat_contr_W_im, & 1022 mat_I_muP_occ_re, mat_I_muP_virt_re, mat_I_muP_occ_im, mat_I_muP_virt_im, & 1023 index_to_cell_R1, index_to_cell_R2, index_to_cell_3c, & 1024 index_to_cell_dm, index_to_cell_R1_plus_S2, cell_to_index_3c, & 1025 cell_to_index_R2, gw_corr_lev_occ, gw_corr_lev_tot, & 1026 homo, num_integ_points, num_fit_points, tj, tau_tj, ikp, & 1027 has_3c_blocks_re, has_3c_blocks_im, & 1028 cycle_R1_S2_n_level, has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, & 1029 cycle_R1_plus_S2_n_level, kpoints, & 1030 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, eps_filter, para_env, & 1031 are_flops_I_T_R1_plus_S2_S1_n_level, start_jquad) 1032 1033 COMPLEX(KIND=dp), DIMENSION(:, :) :: vec_Sigma_c_gw 1034 REAL(KIND=dp), DIMENSION(:) :: vec_Sigma_x_gw, vec_Sigma_x_minus_vxc_gw 1035 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_gw_kp_re, & 1036 mat_3c_overl_int_gw_kp_im 1037 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_S, mat_dm_virt_S, mat_W_R 1038 TYPE(dbcsr_type), POINTER :: mat_contr_gf_occ_re, mat_contr_gf_occ_im, mat_contr_gf_virt_re, & 1039 mat_contr_gf_virt_im, mat_contr_W_re, mat_contr_W_im 1040 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_I_muP_occ_re, mat_I_muP_virt_re, & 1041 mat_I_muP_occ_im, mat_I_muP_virt_im 1042 INTEGER, DIMENSION(:, :) :: index_to_cell_R1, index_to_cell_R2, & 1043 index_to_cell_3c, index_to_cell_dm, & 1044 index_to_cell_R1_plus_S2 1045 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c, cell_to_index_R2 1046 INTEGER :: gw_corr_lev_occ, gw_corr_lev_tot, homo, & 1047 num_integ_points, num_fit_points 1048 REAL(KIND=dp), DIMENSION(:) :: tj 1049 REAL(KIND=dp), DIMENSION(0:num_integ_points) :: tau_tj 1050 INTEGER :: ikp 1051 LOGICAL, DIMENSION(:, :, :) :: has_3c_blocks_re, has_3c_blocks_im 1052 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :) :: cycle_R1_S2_n_level 1053 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: has_blocks_mat_dm_occ_S, & 1054 has_blocks_mat_dm_virt_S 1055 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: cycle_R1_plus_S2_n_level 1056 TYPE(kpoint_type), POINTER :: kpoints 1057 REAL(KIND=dp), DIMENSION(:, :) :: weights_cos_tf_t_to_w, & 1058 weights_sin_tf_t_to_w 1059 REAL(KIND=dp) :: eps_filter 1060 TYPE(cp_para_env_type), POINTER :: para_env 1061 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: are_flops_I_T_R1_plus_S2_S1_n_level 1062 INTEGER :: start_jquad 1063 1064 CHARACTER(LEN=*), PARAMETER :: routineN = 'cell_sum_self_ener', & 1065 routineP = moduleN//':'//routineN 1066 1067 COMPLEX(KIND=dp) :: im_unit 1068 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_Sigma_c_gw_cos_omega, & 1069 vec_Sigma_c_gw_cos_tau, vec_Sigma_c_gw_neg_tau, vec_Sigma_c_gw_pos_tau, & 1070 vec_Sigma_c_gw_sin_omega, vec_Sigma_c_gw_sin_tau 1071 INTEGER :: bound_1, bound_2, bound_3, bound_4, bound_5, bound_6, handle, i_cell_R1_plus_S2, & 1072 i_cell_S2, iquad, jquad, n_level_gw, num_cells_3c, num_cells_dm, num_cells_R1, & 1073 num_cells_R1_plus_S2, num_cells_R2, x_cell_R1, x_cell_R1_plus_S2, y_cell_R1, & 1074 y_cell_R1_plus_S2, z_cell_R1, z_cell_R1_plus_S2 1075 REAL(KIND=dp) :: omega, tau, weight_cos, weight_sin 1076 1077 CALL timeset(routineN, handle) 1078 1079 ALLOCATE (vec_Sigma_c_gw_pos_tau(gw_corr_lev_tot, num_integ_points)) 1080 vec_Sigma_c_gw_pos_tau = z_zero 1081 ALLOCATE (vec_Sigma_c_gw_neg_tau(gw_corr_lev_tot, num_integ_points)) 1082 vec_Sigma_c_gw_neg_tau = z_zero 1083 ALLOCATE (vec_Sigma_c_gw_cos_tau(gw_corr_lev_tot, num_integ_points)) 1084 vec_Sigma_c_gw_cos_tau = z_zero 1085 ALLOCATE (vec_Sigma_c_gw_sin_tau(gw_corr_lev_tot, num_integ_points)) 1086 vec_Sigma_c_gw_sin_tau = z_zero 1087 1088 ALLOCATE (vec_Sigma_c_gw_cos_omega(gw_corr_lev_tot, num_integ_points)) 1089 vec_Sigma_c_gw_cos_omega = z_zero 1090 ALLOCATE (vec_Sigma_c_gw_sin_omega(gw_corr_lev_tot, num_integ_points)) 1091 vec_Sigma_c_gw_sin_omega = z_zero 1092 1093 num_cells_R1 = SIZE(index_to_cell_R1, 2) 1094 num_cells_R2 = SIZE(index_to_cell_R2, 2) 1095 num_cells_dm = SIZE(index_to_cell_dm, 2) 1096 num_cells_R1_plus_S2 = SIZE(index_to_cell_R1_plus_S2, 2) 1097 num_cells_3c = SIZE(index_to_cell_3c, 2) 1098 1099 bound_1 = LBOUND(cell_to_index_3c, 1) 1100 bound_2 = UBOUND(cell_to_index_3c, 1) 1101 bound_3 = LBOUND(cell_to_index_3c, 2) 1102 bound_4 = UBOUND(cell_to_index_3c, 2) 1103 bound_5 = LBOUND(cell_to_index_3c, 3) 1104 bound_6 = UBOUND(cell_to_index_3c, 3) 1105 1106 ! jquad = 0 corresponds to exact exchange self-energy 1107 DO jquad = start_jquad, num_integ_points 1108 1109 DO i_cell_R1_plus_S2 = 1, num_cells_R1_plus_S2 1110 1111 x_cell_R1_plus_S2 = index_to_cell_R1_plus_S2(1, i_cell_R1_plus_S2) 1112 y_cell_R1_plus_S2 = index_to_cell_R1_plus_S2(2, i_cell_R1_plus_S2) 1113 z_cell_R1_plus_S2 = index_to_cell_R1_plus_S2(3, i_cell_R1_plus_S2) 1114 1115 tau = tau_tj(jquad) 1116 1117 DO n_level_gw = 1, gw_corr_lev_tot 1118 1119 IF (cycle_R1_plus_S2_n_level(i_cell_R1_plus_S2, n_level_gw, jquad)) CYCLE 1120 1121 CALL compute_I_muP_T_R1_plus_S2(i_cell_R1_plus_S2, x_cell_R1_plus_S2, y_cell_R1_plus_S2, z_cell_R1_plus_S2, & 1122 mat_I_muP_occ_re, mat_I_muP_virt_re, mat_I_muP_occ_im, mat_I_muP_virt_im, & 1123 mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 1124 mat_dm_occ_S, mat_dm_virt_S, jquad, n_level_gw, cell_to_index_3c, & 1125 index_to_cell_dm, has_3c_blocks_re, has_3c_blocks_im, & 1126 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, num_cells_dm, para_env, & 1127 are_flops_I_T_R1_plus_S2_S1_n_level, eps_filter) 1128 1129 DO i_cell_S2 = 1, num_cells_3c 1130 1131 IF (cycle_R1_S2_n_level(i_cell_R1_plus_S2, i_cell_S2, n_level_gw, jquad)) CYCLE 1132 1133 x_cell_R1 = x_cell_R1_plus_S2 - index_to_cell_3c(1, i_cell_S2) 1134 y_cell_R1 = y_cell_R1_plus_S2 - index_to_cell_3c(2, i_cell_S2) 1135 z_cell_R1 = z_cell_R1_plus_S2 - index_to_cell_3c(3, i_cell_S2) 1136 1137 CALL trafo_I_T_R1_plus_S2_to_M_R1_S2(mat_I_muP_occ_re, mat_I_muP_virt_re, & 1138 mat_I_muP_occ_im, mat_I_muP_virt_im, & 1139 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 1140 mat_contr_gf_virt_re, mat_contr_gf_virt_im, & 1141 index_to_cell_3c, kpoints, ikp, & 1142 x_cell_R1, y_cell_R1, z_cell_R1) 1143 1144 ! perform multiplication with W 1145 CALL mult_3c_with_W(mat_W_R, mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 1146 mat_contr_W_re, mat_contr_W_im, n_level_gw, jquad, num_cells_3c, & 1147 x_cell_R1, y_cell_R1, z_cell_R1, x_cell_R1_plus_S2, y_cell_R1_plus_S2, & 1148 z_cell_R1_plus_S2, index_to_cell_3c, cell_to_index_3c, & 1149 cell_to_index_R2, has_3c_blocks_re, has_3c_blocks_im) 1150 1151 ! perform contraction to self-energy and do the FT from im. time to im. frequency 1152 CALL trace_for_self_ener(mat_contr_W_re, mat_contr_W_im, n_level_gw, jquad, & 1153 gw_corr_lev_occ, homo, & 1154 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 1155 mat_contr_gf_virt_re, mat_contr_gf_virt_im, & 1156 vec_Sigma_c_gw_pos_tau, vec_Sigma_c_gw_neg_tau, & 1157 vec_Sigma_x_gw, ikp, & 1158 cycle_R1_S2_n_level, cycle_R1_plus_S2_n_level, & 1159 eps_filter, i_cell_R1_plus_S2, i_cell_S2, & 1160 start_jquad) 1161 1162 END DO ! R1 (or S2=(R1+S2)-R1 1163 END DO ! n_level_gw 1164 END DO ! jquad 1165 END DO ! R1+S2 1166 1167 vec_Sigma_x_minus_vxc_gw(:) = vec_Sigma_x_minus_vxc_gw(:) + vec_Sigma_x_gw(:) 1168 1169 im_unit = (0.0_dp, 1.0_dp) 1170 1171 vec_Sigma_c_gw_cos_tau(:, 1:num_integ_points) = 0.5_dp*(vec_Sigma_c_gw_pos_tau(:, 1:num_integ_points) + & 1172 vec_Sigma_c_gw_neg_tau(:, 1:num_integ_points)) 1173 vec_Sigma_c_gw_sin_tau(:, 1:num_integ_points) = 0.5_dp*(vec_Sigma_c_gw_pos_tau(:, 1:num_integ_points) - & 1174 vec_Sigma_c_gw_neg_tau(:, 1:num_integ_points)) 1175 1176 ! Fourier transform from time to frequency 1177 DO jquad = 1, num_integ_points 1178 1179 DO iquad = 1, num_integ_points 1180 1181 omega = tj(jquad) 1182 tau = tau_tj(iquad) 1183 weight_cos = weights_cos_tf_t_to_w(jquad, iquad)*COS(omega*tau) 1184 weight_sin = weights_sin_tf_t_to_w(jquad, iquad)*SIN(omega*tau) 1185 1186 vec_Sigma_c_gw_cos_omega(:, jquad) = vec_Sigma_c_gw_cos_omega(:, jquad) + & 1187 weight_cos*vec_Sigma_c_gw_cos_tau(:, iquad) 1188 1189 vec_Sigma_c_gw_sin_omega(:, jquad) = vec_Sigma_c_gw_sin_omega(:, jquad) + & 1190 weight_sin*vec_Sigma_c_gw_sin_tau(:, iquad) 1191 1192 END DO 1193 1194 END DO 1195 1196 ! for occupied levels, we need the correlation self-energy for negative omega. Therefore, weight_sin 1197 ! should be computed with -omega, which results in an additional minus for vec_Sigma_c_gw_sin_omega: 1198 vec_Sigma_c_gw_sin_omega(1:gw_corr_lev_occ, :) = -vec_Sigma_c_gw_sin_omega(1:gw_corr_lev_occ, :) 1199 1200 ! the third index k-point is already absorbed when calling the subroutine 1201 vec_Sigma_c_gw(:, 1:num_fit_points) = vec_Sigma_c_gw_cos_omega(:, 1:num_fit_points) + & 1202 im_unit*vec_Sigma_c_gw_sin_omega(:, 1:num_fit_points) 1203 1204 DEALLOCATE (vec_Sigma_c_gw_pos_tau) 1205 DEALLOCATE (vec_Sigma_c_gw_neg_tau) 1206 DEALLOCATE (vec_Sigma_c_gw_cos_tau) 1207 DEALLOCATE (vec_Sigma_c_gw_sin_tau) 1208 DEALLOCATE (vec_Sigma_c_gw_cos_omega) 1209 DEALLOCATE (vec_Sigma_c_gw_sin_omega) 1210 1211 CALL timestop(handle) 1212 1213 END SUBROUTINE 1214 1215! ************************************************************************************************** 1216!> \brief ... 1217!> \param mat_contr_W_re ... 1218!> \param mat_contr_W_im ... 1219!> \param n_level_gw ... 1220!> \param jquad ... 1221!> \param gw_corr_lev_occ ... 1222!> \param homo ... 1223!> \param mat_contr_gf_occ_re ... 1224!> \param mat_contr_gf_occ_im ... 1225!> \param mat_contr_gf_virt_re ... 1226!> \param mat_contr_gf_virt_im ... 1227!> \param vec_Sigma_c_gw_pos_tau ... 1228!> \param vec_Sigma_c_gw_neg_tau ... 1229!> \param vec_Sigma_x_gw ... 1230!> \param ikp ... 1231!> \param cycle_R1_S2_n_level ... 1232!> \param cycle_R1_plus_S2_n_level ... 1233!> \param eps_filter ... 1234!> \param i_cell_R1_plus_S2 ... 1235!> \param i_cell_S2 ... 1236!> \param start_jquad ... 1237! ************************************************************************************************** 1238 SUBROUTINE trace_for_self_ener(mat_contr_W_re, mat_contr_W_im, n_level_gw, jquad, & 1239 gw_corr_lev_occ, homo, & 1240 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 1241 mat_contr_gf_virt_re, mat_contr_gf_virt_im, & 1242 vec_Sigma_c_gw_pos_tau, vec_Sigma_c_gw_neg_tau, & 1243 vec_Sigma_x_gw, ikp, & 1244 cycle_R1_S2_n_level, cycle_R1_plus_S2_n_level, & 1245 eps_filter, i_cell_R1_plus_S2, i_cell_S2, & 1246 start_jquad) 1247 1248 TYPE(dbcsr_type), POINTER :: mat_contr_W_re, mat_contr_W_im 1249 INTEGER :: n_level_gw, jquad, gw_corr_lev_occ, homo 1250 TYPE(dbcsr_type), POINTER :: mat_contr_gf_occ_re, & 1251 mat_contr_gf_occ_im, & 1252 mat_contr_gf_virt_re, & 1253 mat_contr_gf_virt_im 1254 COMPLEX(KIND=dp), DIMENSION(:, :) :: vec_Sigma_c_gw_pos_tau, & 1255 vec_Sigma_c_gw_neg_tau 1256 REAL(KIND=dp), DIMENSION(:) :: vec_Sigma_x_gw 1257 INTEGER :: ikp 1258 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :) :: cycle_R1_S2_n_level 1259 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: cycle_R1_plus_S2_n_level 1260 REAL(KIND=dp) :: eps_filter 1261 INTEGER :: i_cell_R1_plus_S2, i_cell_S2, start_jquad 1262 1263 CHARACTER(LEN=*), PARAMETER :: routineN = 'trace_for_self_ener', & 1264 routineP = moduleN//':'//routineN 1265 1266 INTEGER :: handle, n_level_gw_ref 1267 REAL(KIND=dp) :: trace_neg_tau_im_1, trace_neg_tau_im_2, trace_neg_tau_re_1, & 1268 trace_neg_tau_re_2, trace_pos_tau_im_1, trace_pos_tau_im_2, trace_pos_tau_re_1, & 1269 trace_pos_tau_re_2 1270 1271 CALL timeset(routineN, handle) 1272 1273 CALL dbcsr_dot(mat_contr_gf_occ_re, & 1274 mat_contr_W_re, & 1275 trace_neg_tau_re_1) 1276 1277 CALL dbcsr_dot(mat_contr_gf_occ_im, & 1278 mat_contr_W_im, & 1279 trace_neg_tau_re_2) 1280 1281 CALL dbcsr_dot(mat_contr_gf_occ_re, & 1282 mat_contr_W_im, & 1283 trace_neg_tau_im_1) 1284 1285 CALL dbcsr_dot(mat_contr_gf_occ_im, & 1286 mat_contr_W_re, & 1287 trace_neg_tau_im_2) 1288 1289 IF (ABS(trace_neg_tau_re_1) + ABS(trace_neg_tau_re_2) + ABS(trace_neg_tau_im_1) + & 1290 ABS(trace_neg_tau_im_2) < eps_filter) THEN 1291 cycle_R1_S2_n_level(i_cell_R1_plus_S2, i_cell_S2, n_level_gw, jquad) = .TRUE. 1292 END IF 1293 1294 IF (ikp == 1 .AND. jquad == start_jquad .AND. i_cell_S2 == 1) THEN 1295 cycle_R1_plus_S2_n_level(i_cell_R1_plus_S2, n_level_gw, jquad) = .TRUE. 1296 END IF 1297 1298 IF (ABS(trace_neg_tau_re_1) + ABS(trace_neg_tau_re_2) + ABS(trace_neg_tau_im_1) + & 1299 ABS(trace_neg_tau_im_2) > eps_filter) THEN 1300 cycle_R1_plus_S2_n_level(i_cell_R1_plus_S2, n_level_gw, jquad) = .FALSE. 1301 END IF 1302 1303 IF (jquad == 0) THEN 1304 1305 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 1306 1307 vec_Sigma_x_gw(n_level_gw_ref) = vec_Sigma_x_gw(n_level_gw_ref) + trace_neg_tau_re_1 - trace_neg_tau_re_2 1308 1309 ELSE 1310 1311 vec_Sigma_c_gw_neg_tau(n_level_gw, jquad) = vec_Sigma_c_gw_neg_tau(n_level_gw, jquad) + & 1312 CMPLX(trace_neg_tau_re_1 - trace_neg_tau_re_2, & 1313 trace_neg_tau_im_1 + trace_neg_tau_im_2, dp) 1314 1315 END IF 1316 1317 CALL dbcsr_dot(mat_contr_gf_virt_re, & 1318 mat_contr_W_re, & 1319 trace_pos_tau_re_1) 1320 1321 CALL dbcsr_dot(mat_contr_gf_virt_im, & 1322 mat_contr_W_im, & 1323 trace_pos_tau_re_2) 1324 1325 CALL dbcsr_dot(mat_contr_gf_virt_re, & 1326 mat_contr_W_im, & 1327 trace_pos_tau_im_1) 1328 1329 CALL dbcsr_dot(mat_contr_gf_virt_im, & 1330 mat_contr_W_re, & 1331 trace_pos_tau_im_2) 1332 1333 IF (jquad > 0) THEN 1334 1335 vec_Sigma_c_gw_pos_tau(n_level_gw, jquad) = vec_Sigma_c_gw_pos_tau(n_level_gw, jquad) + & 1336 CMPLX(trace_pos_tau_re_1 - trace_pos_tau_re_2, & 1337 trace_pos_tau_im_1 + trace_pos_tau_im_2, dp) 1338 1339 END IF 1340 1341 CALL timestop(handle) 1342 1343 END SUBROUTINE 1344 1345! ************************************************************************************************** 1346!> \brief ... 1347!> \param mat_W_R ... 1348!> \param mat_3c_overl_int_gw_kp_re ... 1349!> \param mat_3c_overl_int_gw_kp_im ... 1350!> \param mat_contr_W_re ... 1351!> \param mat_contr_W_im ... 1352!> \param n_level_gw ... 1353!> \param jquad ... 1354!> \param num_cells_3c ... 1355!> \param x_cell_R1 ... 1356!> \param y_cell_R1 ... 1357!> \param z_cell_R1 ... 1358!> \param x_cell_R1_plus_S2 ... 1359!> \param y_cell_R1_plus_S2 ... 1360!> \param z_cell_R1_plus_S2 ... 1361!> \param index_to_cell_3c ... 1362!> \param cell_to_index_3c ... 1363!> \param cell_to_index_R2 ... 1364!> \param has_3c_blocks_re ... 1365!> \param has_3c_blocks_im ... 1366! ************************************************************************************************** 1367 SUBROUTINE mult_3c_with_W(mat_W_R, mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 1368 mat_contr_W_re, mat_contr_W_im, n_level_gw, jquad, num_cells_3c, & 1369 x_cell_R1, y_cell_R1, z_cell_R1, & 1370 x_cell_R1_plus_S2, y_cell_R1_plus_S2, z_cell_R1_plus_S2, & 1371 index_to_cell_3c, cell_to_index_3c, & 1372 cell_to_index_R2, has_3c_blocks_re, has_3c_blocks_im) 1373 1374 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_W_R 1375 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_gw_kp_re, & 1376 mat_3c_overl_int_gw_kp_im 1377 TYPE(dbcsr_type), POINTER :: mat_contr_W_re, mat_contr_W_im 1378 INTEGER :: n_level_gw, jquad, num_cells_3c, x_cell_R1, y_cell_R1, z_cell_R1, & 1379 x_cell_R1_plus_S2, y_cell_R1_plus_S2, z_cell_R1_plus_S2 1380 INTEGER, DIMENSION(:, :) :: index_to_cell_3c 1381 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c, cell_to_index_R2 1382 LOGICAL, DIMENSION(:, :, :) :: has_3c_blocks_re, has_3c_blocks_im 1383 1384 CHARACTER(LEN=*), PARAMETER :: routineN = 'mult_3c_with_W', routineP = moduleN//':'//routineN 1385 1386 INTEGER :: bound_1, bound_2, bound_3, bound_4, bound_5, bound_6, bound_R2_1, bound_R2_2, & 1387 bound_R2_3, bound_R2_4, bound_R2_5, bound_R2_6, handle, i_cell_R1_minus_R2, & 1388 i_cell_R1_plus_S2_minus_R2, i_cell_R2, x_cell_R1_plus_S2_minus_R2, x_cell_R2, & 1389 y_cell_R1_plus_S2_minus_R2, y_cell_R2, z_cell_R1_plus_S2_minus_R2, z_cell_R2 1390 1391 CALL timeset(routineN, handle) 1392 1393 bound_1 = LBOUND(cell_to_index_3c, 1) 1394 bound_2 = UBOUND(cell_to_index_3c, 1) 1395 bound_3 = LBOUND(cell_to_index_3c, 2) 1396 bound_4 = UBOUND(cell_to_index_3c, 2) 1397 bound_5 = LBOUND(cell_to_index_3c, 3) 1398 bound_6 = UBOUND(cell_to_index_3c, 3) 1399 1400 bound_R2_1 = LBOUND(cell_to_index_R2, 1) 1401 bound_R2_2 = UBOUND(cell_to_index_R2, 1) 1402 bound_R2_3 = LBOUND(cell_to_index_R2, 2) 1403 bound_R2_4 = UBOUND(cell_to_index_R2, 2) 1404 bound_R2_5 = LBOUND(cell_to_index_R2, 3) 1405 bound_R2_6 = UBOUND(cell_to_index_R2, 3) 1406 1407 CALL dbcsr_set(mat_contr_W_re, 0.0_dp) 1408 CALL dbcsr_set(mat_contr_W_im, 0.0_dp) 1409 1410 DO i_cell_R1_minus_R2 = 1, num_cells_3c 1411 1412 x_cell_R2 = x_cell_R1 - index_to_cell_3c(1, i_cell_R1_minus_R2) 1413 y_cell_R2 = y_cell_R1 - index_to_cell_3c(2, i_cell_R1_minus_R2) 1414 z_cell_R2 = z_cell_R1 - index_to_cell_3c(3, i_cell_R1_minus_R2) 1415 1416 IF (x_cell_R2 < bound_R2_1 .OR. & 1417 x_cell_R2 > bound_R2_2 .OR. & 1418 y_cell_R2 < bound_R2_3 .OR. & 1419 y_cell_R2 > bound_R2_4 .OR. & 1420 z_cell_R2 < bound_R2_5 .OR. & 1421 z_cell_R2 > bound_R2_6) THEN 1422 1423 CYCLE 1424 1425 END IF 1426 1427 i_cell_R2 = cell_to_index_R2(x_cell_R2, y_cell_R2, z_cell_R2) 1428 1429 x_cell_R1_plus_S2_minus_R2 = x_cell_R1_plus_S2 - x_cell_R2 1430 y_cell_R1_plus_S2_minus_R2 = y_cell_R1_plus_S2 - y_cell_R2 1431 z_cell_R1_plus_S2_minus_R2 = z_cell_R1_plus_S2 - z_cell_R2 1432 1433 IF (x_cell_R1_plus_S2_minus_R2 < bound_1 .OR. & 1434 x_cell_R1_plus_S2_minus_R2 > bound_2 .OR. & 1435 y_cell_R1_plus_S2_minus_R2 < bound_3 .OR. & 1436 y_cell_R1_plus_S2_minus_R2 > bound_4 .OR. & 1437 z_cell_R1_plus_S2_minus_R2 < bound_5 .OR. & 1438 z_cell_R1_plus_S2_minus_R2 > bound_6) THEN 1439 1440 CYCLE 1441 1442 END IF 1443 1444 i_cell_R1_plus_S2_minus_R2 = cell_to_index_3c(x_cell_R1_plus_S2_minus_R2, & 1445 y_cell_R1_plus_S2_minus_R2, & 1446 z_cell_R1_plus_S2_minus_R2) 1447 1448 IF (i_cell_R1_plus_S2_minus_R2 == 0) CYCLE 1449 1450 IF (has_3c_blocks_re(n_level_gw, i_cell_R1_minus_R2, i_cell_R1_plus_S2_minus_R2)) THEN 1451 1452 CALL dbcsr_multiply("N", "N", 1.0_dp, mat_W_R(i_cell_R2, jquad)%matrix, & 1453 mat_3c_overl_int_gw_kp_re(n_level_gw, i_cell_R1_minus_R2, i_cell_R1_plus_S2_minus_R2)%matrix, & 1454 1.0_dp, mat_contr_W_re) 1455 1456 END IF 1457 1458 IF (has_3c_blocks_im(n_level_gw, i_cell_R1_minus_R2, i_cell_R1_plus_S2_minus_R2)) THEN 1459 1460 CALL dbcsr_multiply("N", "N", 1.0_dp, mat_W_R(i_cell_R2, jquad)%matrix, & 1461 mat_3c_overl_int_gw_kp_im(n_level_gw, i_cell_R1_minus_R2, i_cell_R1_plus_S2_minus_R2)%matrix, & 1462 1.0_dp, mat_contr_W_im) 1463 1464 END IF 1465 1466 END DO 1467 1468 CALL timestop(handle) 1469 1470 END SUBROUTINE 1471 1472! ************************************************************************************************** 1473!> \brief ... 1474!> \param mat_I_muP_occ_re ... 1475!> \param mat_I_muP_virt_re ... 1476!> \param mat_I_muP_occ_im ... 1477!> \param mat_I_muP_virt_im ... 1478!> \param mat_contr_gf_occ_re ... 1479!> \param mat_contr_gf_occ_im ... 1480!> \param mat_contr_gf_virt_re ... 1481!> \param mat_contr_gf_virt_im ... 1482!> \param index_to_cell_3c ... 1483!> \param kpoints ... 1484!> \param ikp ... 1485!> \param x_cell_R1 ... 1486!> \param y_cell_R1 ... 1487!> \param z_cell_R1 ... 1488! ************************************************************************************************** 1489 SUBROUTINE trafo_I_T_R1_plus_S2_to_M_R1_S2(mat_I_muP_occ_re, mat_I_muP_virt_re, & 1490 mat_I_muP_occ_im, mat_I_muP_virt_im, & 1491 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 1492 mat_contr_gf_virt_re, mat_contr_gf_virt_im, & 1493 index_to_cell_3c, kpoints, ikp, & 1494 x_cell_R1, y_cell_R1, z_cell_R1) 1495 1496 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_I_muP_occ_re, mat_I_muP_virt_re, & 1497 mat_I_muP_occ_im, mat_I_muP_virt_im 1498 TYPE(dbcsr_type), POINTER :: mat_contr_gf_occ_re, & 1499 mat_contr_gf_occ_im, & 1500 mat_contr_gf_virt_re, & 1501 mat_contr_gf_virt_im 1502 INTEGER, DIMENSION(:, :) :: index_to_cell_3c 1503 TYPE(kpoint_type), POINTER :: kpoints 1504 INTEGER :: ikp, x_cell_R1, y_cell_R1, z_cell_R1 1505 1506 CHARACTER(LEN=*), PARAMETER :: routineN = 'trafo_I_T_R1_plus_S2_to_M_R1_S2', & 1507 routineP = moduleN//':'//routineN 1508 1509 INTEGER :: handle, i_cell_T, num_cells_3c, xcell, & 1510 ycell, zcell 1511 REAL(KIND=dp) :: arg, coskl, sinkl 1512 REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp 1513 1514 CALL timeset(routineN, handle) 1515 1516 num_cells_3c = SIZE(mat_I_muP_occ_re) 1517 CALL get_kpoint_info(kpoints, xkp=xkp) 1518 1519 CALL dbcsr_set(mat_contr_gf_occ_re, 0.0_dp) 1520 CALL dbcsr_set(mat_contr_gf_occ_im, 0.0_dp) 1521 CALL dbcsr_set(mat_contr_gf_virt_re, 0.0_dp) 1522 CALL dbcsr_set(mat_contr_gf_virt_im, 0.0_dp) 1523 1524 DO i_cell_T = 1, num_cells_3c 1525 1526 xcell = index_to_cell_3c(1, i_cell_T) - x_cell_R1 1527 ycell = index_to_cell_3c(2, i_cell_T) - y_cell_R1 1528 zcell = index_to_cell_3c(3, i_cell_T) - z_cell_R1 1529 arg = REAL(xcell, dp)*xkp(1, ikp) + REAL(ycell, dp)*xkp(2, ikp) + REAL(zcell, dp)*xkp(3, ikp) 1530 coskl = COS(twopi*arg) 1531 sinkl = SIN(twopi*arg) 1532 1533 CALL dbcsr_scale_and_add_local(mat_contr_gf_occ_re, mat_I_muP_occ_re(i_cell_T)%matrix, coskl) 1534 CALL dbcsr_scale_and_add_local(mat_contr_gf_occ_re, mat_I_muP_occ_im(i_cell_T)%matrix, -sinkl) 1535 CALL dbcsr_scale_and_add_local(mat_contr_gf_occ_im, mat_I_muP_occ_re(i_cell_T)%matrix, sinkl) 1536 CALL dbcsr_scale_and_add_local(mat_contr_gf_occ_im, mat_I_muP_occ_im(i_cell_T)%matrix, coskl) 1537 1538 CALL dbcsr_scale_and_add_local(mat_contr_gf_virt_re, mat_I_muP_virt_re(i_cell_T)%matrix, coskl) 1539 CALL dbcsr_scale_and_add_local(mat_contr_gf_virt_re, mat_I_muP_virt_im(i_cell_T)%matrix, -sinkl) 1540 CALL dbcsr_scale_and_add_local(mat_contr_gf_virt_im, mat_I_muP_virt_re(i_cell_T)%matrix, sinkl) 1541 CALL dbcsr_scale_and_add_local(mat_contr_gf_virt_im, mat_I_muP_virt_im(i_cell_T)%matrix, coskl) 1542 1543 END DO 1544 1545 CALL timestop(handle) 1546 1547 END SUBROUTINE 1548 1549! ************************************************************************************************** 1550!> \brief ... 1551!> \param mat_A ... 1552!> \param mat_B ... 1553!> \param beta ... 1554! ************************************************************************************************** 1555 SUBROUTINE dbcsr_scale_and_add_local(mat_A, mat_B, beta) 1556 TYPE(dbcsr_type), POINTER :: mat_A, mat_B 1557 REAL(KIND=dp) :: beta 1558 1559 INTEGER :: col, row 1560 LOGICAL :: found 1561 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_to_compute, data_block 1562 TYPE(dbcsr_iterator_type) :: iter 1563 1564 CALL dbcsr_iterator_start(iter, mat_B) 1565 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1566 1567 CALL dbcsr_iterator_next_block(iter, row, col, data_block) 1568 1569 NULLIFY (block_to_compute) 1570 CALL dbcsr_get_block_p(matrix=mat_A, & 1571 row=row, col=col, block=block_to_compute, found=found) 1572 1573 CPASSERT(found) 1574 1575 block_to_compute(:, :) = block_to_compute(:, :) + beta*data_block(:, :) 1576 1577 END DO 1578 1579 CALL dbcsr_iterator_stop(iter) 1580 1581 END SUBROUTINE 1582 1583! ************************************************************************************************** 1584!> \brief ... 1585!> \param i_cell_R1_plus_S2 ... 1586!> \param x_cell_R1_plus_S2 ... 1587!> \param y_cell_R1_plus_S2 ... 1588!> \param z_cell_R1_plus_S2 ... 1589!> \param mat_I_muP_occ_re ... 1590!> \param mat_I_muP_virt_re ... 1591!> \param mat_I_muP_occ_im ... 1592!> \param mat_I_muP_virt_im ... 1593!> \param mat_3c_overl_int_gw_kp_re ... 1594!> \param mat_3c_overl_int_gw_kp_im ... 1595!> \param mat_dm_occ_S ... 1596!> \param mat_dm_virt_S ... 1597!> \param jquad ... 1598!> \param n_level_gw ... 1599!> \param cell_to_index_3c ... 1600!> \param index_to_cell_dm ... 1601!> \param has_3c_blocks_re ... 1602!> \param has_3c_blocks_im ... 1603!> \param has_blocks_mat_dm_occ_S ... 1604!> \param has_blocks_mat_dm_virt_S ... 1605!> \param num_cells_dm ... 1606!> \param para_env ... 1607!> \param are_flops_I_T_R1_plus_S2_S1_n_level ... 1608!> \param eps_filter ... 1609! ************************************************************************************************** 1610 SUBROUTINE compute_I_muP_T_R1_plus_S2(i_cell_R1_plus_S2, x_cell_R1_plus_S2, y_cell_R1_plus_S2, z_cell_R1_plus_S2, & 1611 mat_I_muP_occ_re, mat_I_muP_virt_re, mat_I_muP_occ_im, mat_I_muP_virt_im, & 1612 mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 1613 mat_dm_occ_S, mat_dm_virt_S, jquad, n_level_gw, cell_to_index_3c, & 1614 index_to_cell_dm, has_3c_blocks_re, has_3c_blocks_im, & 1615 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, num_cells_dm, para_env, & 1616 are_flops_I_T_R1_plus_S2_S1_n_level, eps_filter) 1617 1618 INTEGER :: i_cell_R1_plus_S2, x_cell_R1_plus_S2, & 1619 y_cell_R1_plus_S2, z_cell_R1_plus_S2 1620 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_I_muP_occ_re, mat_I_muP_virt_re, & 1621 mat_I_muP_occ_im, mat_I_muP_virt_im 1622 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_gw_kp_re, & 1623 mat_3c_overl_int_gw_kp_im 1624 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_S, mat_dm_virt_S 1625 INTEGER :: jquad, n_level_gw 1626 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c 1627 INTEGER, DIMENSION(:, :) :: index_to_cell_dm 1628 LOGICAL, DIMENSION(:, :, :) :: has_3c_blocks_re, has_3c_blocks_im 1629 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: has_blocks_mat_dm_occ_S, & 1630 has_blocks_mat_dm_virt_S 1631 INTEGER :: num_cells_dm 1632 TYPE(cp_para_env_type), POINTER :: para_env 1633 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: are_flops_I_T_R1_plus_S2_S1_n_level 1634 REAL(KIND=dp) :: eps_filter 1635 1636 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_I_muP_T_R1_plus_S2', & 1637 routineP = moduleN//':'//routineN 1638 1639 INTEGER :: bound_1, bound_2, bound_3, bound_4, bound_5, bound_6, handle, i_cell_S1, & 1640 i_cell_T, index_2, num_cells_3c, x_cell_2, y_cell_2, z_cell_2 1641 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:, :) :: flops_occ_im, flops_occ_re, & 1642 flops_virt_im, flops_virt_re 1643 1644 CALL timeset(routineN, handle) 1645 1646 num_cells_3c = SIZE(mat_I_muP_occ_re) 1647 1648 bound_1 = LBOUND(cell_to_index_3c, 1) 1649 bound_2 = UBOUND(cell_to_index_3c, 1) 1650 bound_3 = LBOUND(cell_to_index_3c, 2) 1651 bound_4 = UBOUND(cell_to_index_3c, 2) 1652 bound_5 = LBOUND(cell_to_index_3c, 3) 1653 bound_6 = UBOUND(cell_to_index_3c, 3) 1654 1655 ALLOCATE (flops_occ_re(num_cells_3c, num_cells_dm)) 1656 flops_occ_re = 0 1657 ALLOCATE (flops_occ_im(num_cells_3c, num_cells_dm)) 1658 flops_occ_im = 0 1659 ALLOCATE (flops_virt_re(num_cells_3c, num_cells_dm)) 1660 flops_virt_re = 0 1661 ALLOCATE (flops_virt_im(num_cells_3c, num_cells_dm)) 1662 flops_virt_im = 0 1663 1664 DO i_cell_T = 1, num_cells_3c 1665 1666 CALL dbcsr_set(mat_I_muP_occ_re(i_cell_T)%matrix, 0.0_dp) 1667 CALL dbcsr_set(mat_I_muP_occ_im(i_cell_T)%matrix, 0.0_dp) 1668 CALL dbcsr_set(mat_I_muP_virt_re(i_cell_T)%matrix, 0.0_dp) 1669 CALL dbcsr_set(mat_I_muP_virt_im(i_cell_T)%matrix, 0.0_dp) 1670 1671 DO i_cell_S1 = 1, num_cells_dm 1672 1673 x_cell_2 = x_cell_R1_plus_S2 + index_to_cell_dm(1, i_cell_S1) 1674 y_cell_2 = y_cell_R1_plus_S2 + index_to_cell_dm(2, i_cell_S1) 1675 z_cell_2 = z_cell_R1_plus_S2 + index_to_cell_dm(3, i_cell_S1) 1676 1677 IF (x_cell_2 < bound_1 .OR. & 1678 x_cell_2 > bound_2 .OR. & 1679 y_cell_2 < bound_3 .OR. & 1680 y_cell_2 > bound_4 .OR. & 1681 z_cell_2 < bound_5 .OR. & 1682 z_cell_2 > bound_6) THEN 1683 1684 CYCLE 1685 1686 END IF 1687 1688 index_2 = cell_to_index_3c(x_cell_2, y_cell_2, z_cell_2) 1689 IF (index_2 == 0) CYCLE 1690 1691 IF (has_3c_blocks_re(n_level_gw, i_cell_T, index_2) .AND. & 1692 has_blocks_mat_dm_occ_S(jquad, i_cell_S1) .AND. & 1693 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 1)) THEN 1694 1695 ! the occ Gf has no minus, but already include the minus from Sigma = -GW 1696 CALL dbcsr_multiply("N", "N", -1.0_dp, & 1697 mat_3c_overl_int_gw_kp_re(n_level_gw, i_cell_T, index_2)%matrix, & 1698 mat_dm_occ_S(jquad, i_cell_S1)%matrix, & 1699 1.0_dp, mat_I_muP_occ_re(i_cell_T)%matrix, flop=flops_occ_re(i_cell_T, i_cell_S1), & 1700 filter_eps=eps_filter) 1701 1702 END IF 1703 1704 IF (has_3c_blocks_im(n_level_gw, i_cell_T, index_2) .AND. & 1705 has_blocks_mat_dm_occ_S(jquad, i_cell_S1) .AND. & 1706 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 2)) THEN 1707 1708 ! other sign as real part because there is a complexe conjugate on the 3c integral 1709 CALL dbcsr_multiply("N", "N", 1.0_dp, & 1710 mat_3c_overl_int_gw_kp_im(n_level_gw, i_cell_T, index_2)%matrix, & 1711 mat_dm_occ_S(jquad, i_cell_S1)%matrix, & 1712 1.0_dp, mat_I_muP_occ_im(i_cell_T)%matrix, flop=flops_occ_im(i_cell_T, i_cell_S1), & 1713 filter_eps=eps_filter) 1714 1715 END IF 1716 1717 IF (has_3c_blocks_re(n_level_gw, i_cell_T, index_2) .AND. & 1718 has_blocks_mat_dm_virt_S(jquad, i_cell_S1) .AND. & 1719 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 3) .AND. & 1720 jquad > 0) THEN 1721 1722 ! the virt Gf has a minus, but already include the minus from Sigma = -GW 1723 CALL dbcsr_multiply("N", "N", 1.0_dp, & 1724 mat_3c_overl_int_gw_kp_re(n_level_gw, i_cell_T, index_2)%matrix, & 1725 mat_dm_virt_S(jquad, i_cell_S1)%matrix, & 1726 1.0_dp, mat_I_muP_virt_re(i_cell_T)%matrix, flop=flops_virt_re(i_cell_T, i_cell_S1), & 1727 filter_eps=eps_filter) 1728 1729 END IF 1730 1731 IF (has_3c_blocks_im(n_level_gw, i_cell_T, index_2) .AND. & 1732 has_blocks_mat_dm_virt_S(jquad, i_cell_S1) .AND. & 1733 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 4) .AND. & 1734 jquad > 0) THEN 1735 1736 ! other sign as real part because there is a complexe conjugate on the 3c integral 1737 CALL dbcsr_multiply("N", "N", -1.0_dp, & 1738 mat_3c_overl_int_gw_kp_im(n_level_gw, i_cell_T, index_2)%matrix, & 1739 mat_dm_virt_S(jquad, i_cell_S1)%matrix, & 1740 1.0_dp, mat_I_muP_virt_im(i_cell_T)%matrix, flop=flops_virt_im(i_cell_T, i_cell_S1), & 1741 filter_eps=eps_filter) 1742 1743 END IF 1744 1745 END DO 1746 1747 END DO 1748 1749 CALL mp_sum(flops_occ_re, para_env%group) 1750 CALL mp_sum(flops_occ_im, para_env%group) 1751 CALL mp_sum(flops_virt_re, para_env%group) 1752 CALL mp_sum(flops_virt_im, para_env%group) 1753 1754 DO i_cell_T = 1, num_cells_3c 1755 DO i_cell_S1 = 1, num_cells_dm 1756 1757 IF (flops_occ_re(i_cell_T, i_cell_S1) > 0) THEN 1758 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 1) = .TRUE. 1759 ELSE 1760 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 1) = .FALSE. 1761 END IF 1762 1763 IF (flops_occ_im(i_cell_T, i_cell_S1) > 0) THEN 1764 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 2) = .TRUE. 1765 ELSE 1766 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 2) = .FALSE. 1767 END IF 1768 1769 IF (flops_virt_re(i_cell_T, i_cell_S1) > 0 .OR. jquad == 0) THEN 1770 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 3) = .TRUE. 1771 ELSE 1772 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 3) = .FALSE. 1773 END IF 1774 1775 IF (flops_virt_im(i_cell_T, i_cell_S1) > 0 .OR. jquad == 0) THEN 1776 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 4) = .TRUE. 1777 ELSE 1778 are_flops_I_T_R1_plus_S2_S1_n_level(i_cell_T, i_cell_R1_plus_S2, i_cell_S1, n_level_gw, 4) = .FALSE. 1779 END IF 1780 1781 END DO 1782 END DO 1783 1784 DEALLOCATE (flops_occ_re, flops_occ_im, flops_virt_re, flops_virt_im) 1785 1786 CALL timestop(handle) 1787 1788 END SUBROUTINE 1789 1790! ************************************************************************************************** 1791!> \brief ... 1792!> \param mat_dm_occ_S ... 1793!> \param mat_dm_virt_S ... 1794!> \param qs_env ... 1795!> \param num_integ_points ... 1796!> \param stabilize_exp ... 1797!> \param e_fermi ... 1798!> \param eps_filter ... 1799!> \param tau_tj ... 1800!> \param num_cells_dm ... 1801!> \param index_to_cell_dm ... 1802!> \param has_blocks_mat_dm_occ_S ... 1803!> \param has_blocks_mat_dm_virt_S ... 1804!> \param para_env ... 1805! ************************************************************************************************** 1806 SUBROUTINE compute_G_real_space(mat_dm_occ_S, mat_dm_virt_S, qs_env, num_integ_points, stabilize_exp, & 1807 e_fermi, eps_filter, tau_tj, num_cells_dm, index_to_cell_dm, & 1808 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, para_env) 1809 1810 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_S, mat_dm_virt_S 1811 TYPE(qs_environment_type), POINTER :: qs_env 1812 INTEGER :: num_integ_points 1813 REAL(KIND=dp) :: stabilize_exp, e_fermi, eps_filter 1814 REAL(KIND=dp), DIMENSION(0:num_integ_points) :: tau_tj 1815 INTEGER :: num_cells_dm 1816 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 1817 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: has_blocks_mat_dm_occ_S, & 1818 has_blocks_mat_dm_virt_S 1819 TYPE(cp_para_env_type), POINTER :: para_env 1820 1821 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_G_real_space', & 1822 routineP = moduleN//':'//routineN 1823 1824 INTEGER :: handle, i_cell, ispin, jquad, nblks 1825 REAL(KIND=dp) :: tau 1826 1827 CALL timeset(routineN, handle) 1828 1829 ispin = 1 1830 1831 ! get denity matrix for exchange self-energy 1832 tau = 0.0_dp 1833 CALL compute_transl_dm(mat_dm_occ_S, qs_env, ispin, num_integ_points, 0, e_fermi, tau, & 1834 stabilize_exp, eps_filter, num_cells_dm, index_to_cell_dm, & 1835 remove_occ=.FALSE., remove_virt=.TRUE., first_jquad=0) 1836 1837 DO i_cell = 1, num_cells_dm 1838 1839 nblks = dbcsr_get_num_blocks(mat_dm_occ_S(0, i_cell)%matrix) 1840 CALL mp_sum(nblks, para_env%group) 1841 IF (nblks == 0) has_blocks_mat_dm_occ_S(0, i_cell) = .FALSE. 1842 IF (nblks > 0) has_blocks_mat_dm_occ_S(0, i_cell) = .TRUE. 1843 1844 END DO 1845 1846 DO jquad = 1, num_integ_points 1847 1848 tau = tau_tj(jquad) 1849 1850 CALL compute_transl_dm(mat_dm_occ_S, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & 1851 stabilize_exp, eps_filter, num_cells_dm, index_to_cell_dm, & 1852 remove_occ=.FALSE., remove_virt=.TRUE., first_jquad=0) 1853 1854 CALL compute_transl_dm(mat_dm_virt_S, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & 1855 stabilize_exp, eps_filter, num_cells_dm, index_to_cell_dm, & 1856 remove_occ=.TRUE., remove_virt=.FALSE., first_jquad=1) 1857 1858 DO i_cell = 1, num_cells_dm 1859 1860 nblks = dbcsr_get_num_blocks(mat_dm_occ_S(jquad, i_cell)%matrix) 1861 CALL mp_sum(nblks, para_env%group) 1862 IF (nblks == 0) has_blocks_mat_dm_occ_S(jquad, i_cell) = .FALSE. 1863 IF (nblks > 0) has_blocks_mat_dm_occ_S(jquad, i_cell) = .TRUE. 1864 1865 nblks = dbcsr_get_num_blocks(mat_dm_virt_S(jquad, i_cell)%matrix) 1866 CALL mp_sum(nblks, para_env%group) 1867 IF (nblks == 0) has_blocks_mat_dm_virt_S(jquad, i_cell) = .FALSE. 1868 IF (nblks > 0) has_blocks_mat_dm_virt_S(jquad, i_cell) = .TRUE. 1869 1870 END DO 1871 1872 END DO 1873 1874 CALL timestop(handle) 1875 1876 END SUBROUTINE 1877 1878! ************************************************************************************************** 1879!> \brief ... 1880!> \param mat_3c_overl_int_gw_kp_re ... 1881!> \param mat_3c_overl_int_gw_kp_im ... 1882!> \param gw_corr_lev_tot ... 1883!> \param num_3c_repl ... 1884!> \param num_cells_dm ... 1885!> \param num_cells_R1_plus_S2 ... 1886!> \param num_integ_points ... 1887!> \param dimen_RI ... 1888!> \param nmo ... 1889!> \param RI_blk_sizes ... 1890!> \param prim_blk_sizes ... 1891!> \param matrix_s ... 1892!> \param cfm_mat_Q ... 1893!> \param mat_contr_gf_occ_re ... 1894!> \param mat_contr_gf_occ_im ... 1895!> \param mat_contr_gf_virt_re ... 1896!> \param mat_contr_gf_virt_im ... 1897!> \param mat_contr_W_re ... 1898!> \param mat_contr_W_im ... 1899!> \param mat_I_muP_occ_re ... 1900!> \param mat_I_muP_virt_re ... 1901!> \param mat_I_muP_occ_im ... 1902!> \param mat_I_muP_virt_im ... 1903!> \param has_3c_blocks_re ... 1904!> \param has_3c_blocks_im ... 1905!> \param cycle_R1_S2_n_level ... 1906!> \param has_blocks_mat_dm_occ_S ... 1907!> \param has_blocks_mat_dm_virt_S ... 1908!> \param cycle_R1_plus_S2_n_level ... 1909!> \param are_flops_I_T_R1_plus_S2_S1_n_level ... 1910! ************************************************************************************************** 1911 SUBROUTINE allocate_mat_3c_gw(mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 1912 gw_corr_lev_tot, num_3c_repl, num_cells_dm, num_cells_R1_plus_S2, num_integ_points, & 1913 dimen_RI, nmo, RI_blk_sizes, prim_blk_sizes, matrix_s, cfm_mat_Q, & 1914 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 1915 mat_contr_gf_virt_re, mat_contr_gf_virt_im, mat_contr_W_re, mat_contr_W_im, & 1916 mat_I_muP_occ_re, mat_I_muP_virt_re, mat_I_muP_occ_im, mat_I_muP_virt_im, & 1917 has_3c_blocks_re, has_3c_blocks_im, cycle_R1_S2_n_level, & 1918 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, cycle_R1_plus_S2_n_level, & 1919 are_flops_I_T_R1_plus_S2_S1_n_level) 1920 1921 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_gw_kp_re, & 1922 mat_3c_overl_int_gw_kp_im 1923 INTEGER :: gw_corr_lev_tot, num_3c_repl, & 1924 num_cells_dm, num_cells_R1_plus_S2, & 1925 num_integ_points, dimen_RI, nmo 1926 INTEGER, DIMENSION(:), POINTER :: RI_blk_sizes, prim_blk_sizes 1927 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 1928 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 1929 TYPE(dbcsr_type), POINTER :: mat_contr_gf_occ_re, mat_contr_gf_occ_im, mat_contr_gf_virt_re, & 1930 mat_contr_gf_virt_im, mat_contr_W_re, mat_contr_W_im 1931 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_I_muP_occ_re, mat_I_muP_virt_re, & 1932 mat_I_muP_occ_im, mat_I_muP_virt_im 1933 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: has_3c_blocks_re, has_3c_blocks_im 1934 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :) :: cycle_R1_S2_n_level 1935 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: has_blocks_mat_dm_occ_S, & 1936 has_blocks_mat_dm_virt_S 1937 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: cycle_R1_plus_S2_n_level 1938 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: are_flops_I_T_R1_plus_S2_S1_n_level 1939 1940 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_mat_3c_gw', & 1941 routineP = moduleN//':'//routineN 1942 1943 INTEGER :: handle, i_cell, j_cell, n_level_gw 1944 TYPE(cp_fm_struct_type), POINTER :: fm_struct 1945 1946 CALL timeset(routineN, handle) 1947 1948 NULLIFY (mat_3c_overl_int_gw_kp_re) 1949 CALL dbcsr_allocate_matrix_set(mat_3c_overl_int_gw_kp_re, gw_corr_lev_tot, num_3c_repl, num_3c_repl) 1950 1951 NULLIFY (mat_3c_overl_int_gw_kp_im) 1952 CALL dbcsr_allocate_matrix_set(mat_3c_overl_int_gw_kp_im, gw_corr_lev_tot, num_3c_repl, num_3c_repl) 1953 1954 DO n_level_gw = 1, gw_corr_lev_tot 1955 1956 DO i_cell = 1, num_3c_repl 1957 1958 DO j_cell = 1, num_3c_repl 1959 1960 ALLOCATE (mat_3c_overl_int_gw_kp_re(n_level_gw, i_cell, j_cell)%matrix) 1961 CALL dbcsr_create(matrix=mat_3c_overl_int_gw_kp_re(n_level_gw, i_cell, j_cell)%matrix, & 1962 template=matrix_s(1)%matrix, & 1963 matrix_type=dbcsr_type_no_symmetry, & 1964 row_blk_size=RI_blk_sizes, & 1965 col_blk_size=prim_blk_sizes) 1966 1967 ALLOCATE (mat_3c_overl_int_gw_kp_im(n_level_gw, i_cell, j_cell)%matrix) 1968 CALL dbcsr_create(matrix=mat_3c_overl_int_gw_kp_im(n_level_gw, i_cell, j_cell)%matrix, & 1969 template=matrix_s(1)%matrix, & 1970 matrix_type=dbcsr_type_no_symmetry, & 1971 row_blk_size=RI_blk_sizes, & 1972 col_blk_size=prim_blk_sizes) 1973 1974 END DO 1975 1976 END DO 1977 1978 END DO 1979 1980 NULLIFY (mat_contr_gf_occ_re) 1981 CALL dbcsr_init_p(mat_contr_gf_occ_re) 1982 CALL dbcsr_create(matrix=mat_contr_gf_occ_re, & 1983 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 1984 CALL dbcsr_reserve_all_blocks(mat_contr_gf_occ_re) 1985 1986 NULLIFY (mat_contr_gf_occ_im) 1987 CALL dbcsr_init_p(mat_contr_gf_occ_im) 1988 CALL dbcsr_create(matrix=mat_contr_gf_occ_im, & 1989 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 1990 CALL dbcsr_reserve_all_blocks(mat_contr_gf_occ_im) 1991 1992 NULLIFY (mat_contr_gf_virt_re) 1993 CALL dbcsr_init_p(mat_contr_gf_virt_re) 1994 CALL dbcsr_create(matrix=mat_contr_gf_virt_re, & 1995 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 1996 CALL dbcsr_reserve_all_blocks(mat_contr_gf_virt_re) 1997 1998 NULLIFY (mat_contr_gf_virt_im) 1999 CALL dbcsr_init_p(mat_contr_gf_virt_im) 2000 CALL dbcsr_create(matrix=mat_contr_gf_virt_im, & 2001 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 2002 CALL dbcsr_reserve_all_blocks(mat_contr_gf_virt_im) 2003 2004 NULLIFY (mat_contr_W_re) 2005 CALL dbcsr_init_p(mat_contr_W_re) 2006 CALL dbcsr_create(matrix=mat_contr_W_re, & 2007 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 2008 2009 NULLIFY (mat_contr_W_im) 2010 CALL dbcsr_init_p(mat_contr_W_im) 2011 CALL dbcsr_create(matrix=mat_contr_W_im, & 2012 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 2013 2014 NULLIFY (mat_I_muP_occ_re) 2015 CALL dbcsr_allocate_matrix_set(mat_I_muP_occ_re, num_3c_repl) 2016 DO i_cell = 1, num_3c_repl 2017 ALLOCATE (mat_I_muP_occ_re(i_cell)%matrix) 2018 CALL dbcsr_create(matrix=mat_I_muP_occ_re(i_cell)%matrix, & 2019 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 2020 END DO 2021 2022 NULLIFY (mat_I_muP_virt_re) 2023 CALL dbcsr_allocate_matrix_set(mat_I_muP_virt_re, num_3c_repl) 2024 DO i_cell = 1, num_3c_repl 2025 ALLOCATE (mat_I_muP_virt_re(i_cell)%matrix) 2026 CALL dbcsr_create(matrix=mat_I_muP_virt_re(i_cell)%matrix, & 2027 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 2028 END DO 2029 2030 NULLIFY (mat_I_muP_occ_im) 2031 CALL dbcsr_allocate_matrix_set(mat_I_muP_occ_im, num_3c_repl) 2032 DO i_cell = 1, num_3c_repl 2033 ALLOCATE (mat_I_muP_occ_im(i_cell)%matrix) 2034 CALL dbcsr_create(matrix=mat_I_muP_occ_im(i_cell)%matrix, & 2035 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 2036 END DO 2037 2038 NULLIFY (mat_I_muP_virt_im) 2039 CALL dbcsr_allocate_matrix_set(mat_I_muP_virt_im, num_3c_repl) 2040 DO i_cell = 1, num_3c_repl 2041 ALLOCATE (mat_I_muP_virt_im(i_cell)%matrix) 2042 CALL dbcsr_create(matrix=mat_I_muP_virt_im(i_cell)%matrix, & 2043 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix) 2044 END DO 2045 2046 NULLIFY (fm_struct) 2047 CALL cp_fm_struct_create(fm_struct, context=cfm_mat_Q%matrix_struct%context, nrow_global=dimen_RI, & 2048 ncol_global=nmo, para_env=cfm_mat_Q%matrix_struct%para_env) 2049 2050 CALL cp_fm_struct_release(fm_struct) 2051 2052 ALLOCATE (has_3c_blocks_re(gw_corr_lev_tot, num_3c_repl, num_3c_repl)) 2053 has_3c_blocks_re = .FALSE. 2054 ALLOCATE (has_3c_blocks_im(gw_corr_lev_tot, num_3c_repl, num_3c_repl)) 2055 has_3c_blocks_im = .FALSE. 2056 ALLOCATE (has_blocks_mat_dm_occ_S(0:num_integ_points, num_cells_dm)) 2057 has_blocks_mat_dm_occ_S = .FALSE. 2058 ALLOCATE (has_blocks_mat_dm_virt_S(0:num_integ_points, num_cells_dm)) 2059 has_blocks_mat_dm_virt_S = .FALSE. 2060 ALLOCATE (cycle_R1_S2_n_level(num_cells_R1_plus_S2, num_3c_repl, gw_corr_lev_tot, 0:num_integ_points)) 2061 cycle_R1_S2_n_level = .FALSE. 2062 ALLOCATE (cycle_R1_plus_S2_n_level(num_cells_R1_plus_S2, gw_corr_lev_tot, 0:num_integ_points)) 2063 cycle_R1_plus_S2_n_level = .FALSE. 2064 ! 4 is for occ_re, occ_im, virt_re, virt_im 2065 ALLOCATE (are_flops_I_T_R1_plus_S2_S1_n_level(num_3c_repl, num_cells_R1_plus_S2, num_cells_dm, gw_corr_lev_tot, 4)) 2066 are_flops_I_T_R1_plus_S2_S1_n_level = .TRUE. 2067 CALL timestop(handle) 2068 2069 END SUBROUTINE 2070 2071! ************************************************************************************************** 2072!> \brief ... 2073!> \param mat_3c_overl_int_gw_kp_re ... 2074!> \param mat_3c_overl_int_gw_kp_im ... 2075!> \param index_to_cell_R2 ... 2076!> \param index_to_cell_R1_plus_S2 ... 2077!> \param mat_W_R ... 2078!> \param mat_dm_occ_S ... 2079!> \param mat_dm_virt_S ... 2080!> \param mat_contr_gf_occ_re ... 2081!> \param mat_contr_gf_virt_re ... 2082!> \param mat_contr_gf_occ_im ... 2083!> \param mat_contr_gf_virt_im ... 2084!> \param mat_contr_W_re ... 2085!> \param mat_contr_W_im ... 2086!> \param mat_I_muP_occ_re ... 2087!> \param mat_I_muP_virt_re ... 2088!> \param mat_I_muP_occ_im ... 2089!> \param mat_I_muP_virt_im ... 2090!> \param has_3c_blocks_re ... 2091!> \param has_3c_blocks_im ... 2092!> \param cycle_R1_S2_n_level ... 2093!> \param has_blocks_mat_dm_occ_S ... 2094!> \param has_blocks_mat_dm_virt_S ... 2095!> \param cycle_R1_plus_S2_n_level ... 2096!> \param are_flops_I_T_R1_plus_S2_S1_n_level ... 2097! ************************************************************************************************** 2098 SUBROUTINE clean_up_self_energy_kp(mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, & 2099 index_to_cell_R2, index_to_cell_R1_plus_S2, & 2100 mat_W_R, mat_dm_occ_S, mat_dm_virt_S, & 2101 mat_contr_gf_occ_re, mat_contr_gf_virt_re, & 2102 mat_contr_gf_occ_im, mat_contr_gf_virt_im, mat_contr_W_re, mat_contr_W_im, & 2103 mat_I_muP_occ_re, mat_I_muP_virt_re, mat_I_muP_occ_im, mat_I_muP_virt_im, & 2104 has_3c_blocks_re, has_3c_blocks_im, cycle_R1_S2_n_level, & 2105 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, cycle_R1_plus_S2_n_level, & 2106 are_flops_I_T_R1_plus_S2_S1_n_level) 2107 2108 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_gw_kp_re, & 2109 mat_3c_overl_int_gw_kp_im 2110 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_R2, & 2111 index_to_cell_R1_plus_S2 2112 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_W_R, mat_dm_occ_S, mat_dm_virt_S 2113 TYPE(dbcsr_type), POINTER :: mat_contr_gf_occ_re, mat_contr_gf_virt_re, mat_contr_gf_occ_im, & 2114 mat_contr_gf_virt_im, mat_contr_W_re, mat_contr_W_im 2115 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_I_muP_occ_re, mat_I_muP_virt_re, & 2116 mat_I_muP_occ_im, mat_I_muP_virt_im 2117 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: has_3c_blocks_re, has_3c_blocks_im 2118 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :) :: cycle_R1_S2_n_level 2119 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: has_blocks_mat_dm_occ_S, & 2120 has_blocks_mat_dm_virt_S 2121 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: cycle_R1_plus_S2_n_level 2122 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: are_flops_I_T_R1_plus_S2_S1_n_level 2123 2124 CHARACTER(LEN=*), PARAMETER :: routineN = 'clean_up_self_energy_kp', & 2125 routineP = moduleN//':'//routineN 2126 2127 INTEGER :: handle, imatrix, jmatrix 2128 2129 CALL timeset(routineN, handle) 2130 2131 CALL dbcsr_deallocate_matrix_set(mat_3c_overl_int_gw_kp_re) 2132 CALL dbcsr_deallocate_matrix_set(mat_3c_overl_int_gw_kp_im) 2133 DEALLOCATE (index_to_cell_R2) 2134 DEALLOCATE (index_to_cell_R1_plus_S2) 2135 2136! CALL dbcsr_deallocate_matrix_set(mat_dm_occ_S) 2137 DO imatrix = LBOUND(mat_dm_occ_S, 1), UBOUND(mat_dm_occ_S, 1) 2138 DO jmatrix = 1, SIZE(mat_dm_occ_S, 2) 2139 CALL dbcsr_deallocate_matrix(mat_dm_occ_S(imatrix, jmatrix)%matrix) 2140 END DO 2141 END DO 2142 DEALLOCATE (mat_dm_occ_S) 2143 2144 CALL dbcsr_deallocate_matrix_set(mat_dm_virt_S) 2145 2146! CALL dbcsr_deallocate_matrix_set(mat_W_R) 2147 DO imatrix = 1, SIZE(mat_W_R, 1) 2148 DO jmatrix = LBOUND(mat_W_R, 2), UBOUND(mat_W_R, 2) 2149 CALL dbcsr_deallocate_matrix(mat_W_R(imatrix, jmatrix)%matrix) 2150 END DO 2151 END DO 2152 DEALLOCATE (mat_W_R) 2153 2154 CALL dbcsr_release_P(mat_contr_gf_occ_re) 2155 CALL dbcsr_release_P(mat_contr_gf_virt_re) 2156 CALL dbcsr_release_P(mat_contr_gf_occ_im) 2157 CALL dbcsr_release_P(mat_contr_gf_virt_im) 2158 CALL dbcsr_release_P(mat_contr_W_re) 2159 CALL dbcsr_release_P(mat_contr_W_im) 2160 CALL dbcsr_deallocate_matrix_set(mat_I_muP_occ_re) 2161 CALL dbcsr_deallocate_matrix_set(mat_I_muP_virt_re) 2162 CALL dbcsr_deallocate_matrix_set(mat_I_muP_occ_im) 2163 CALL dbcsr_deallocate_matrix_set(mat_I_muP_virt_im) 2164 DEALLOCATE (has_3c_blocks_re) 2165 DEALLOCATE (has_3c_blocks_im) 2166 DEALLOCATE (cycle_R1_S2_n_level) 2167 DEALLOCATE (cycle_R1_plus_S2_n_level) 2168 DEALLOCATE (has_blocks_mat_dm_occ_S) 2169 DEALLOCATE (has_blocks_mat_dm_virt_S) 2170 DEALLOCATE (are_flops_I_T_R1_plus_S2_S1_n_level) 2171 2172 CALL timestop(handle) 2173 2174 END SUBROUTINE 2175 2176! ************************************************************************************************** 2177!> \brief ... 2178!> \param index_to_cell_R1_plus_S2 ... 2179!> \param cell_to_index_3c ... 2180!> \param num_cells_R1_plus_S2 ... 2181!> \param kpoints ... 2182!> \param unit_nr ... 2183!> \param maxcell ... 2184!> \param num_cells_dm ... 2185! ************************************************************************************************** 2186 SUBROUTINE compute_cell_vec_for_R1_plus_S2_or_R1(index_to_cell_R1_plus_S2, cell_to_index_3c, & 2187 num_cells_R1_plus_S2, kpoints, & 2188 unit_nr, maxcell, num_cells_dm) 2189 2190 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_R1_plus_S2 2191 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c 2192 INTEGER :: num_cells_R1_plus_S2 2193 TYPE(kpoint_type), POINTER :: kpoints 2194 INTEGER :: unit_nr, maxcell, num_cells_dm 2195 2196 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cell_vec_for_R1_plus_S2_or_R1', & 2197 routineP = moduleN//':'//routineN 2198 2199 INTEGER :: bound_1, bound_2, bound_3, bound_4, bound_5, bound_6, handle, i_cell_S_1, & 2200 size_set, x_cell_R1_plus_S2, x_cell_R1_plus_S2_minus_S1, y_cell_R1_plus_S2, & 2201 y_cell_R1_plus_S2_minus_S1, z_cell_R1_plus_S2, z_cell_R1_plus_S2_minus_S1 2202 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_tmp 2203 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 2204 LOGICAL :: already_there 2205 2206 CALL timeset(routineN, handle) 2207 2208 index_to_cell_dm => kpoints%index_to_cell 2209 2210 ALLOCATE (index_to_cell_R1_plus_S2(3, 1)) 2211 index_to_cell_R1_plus_S2 = 0 2212 2213 bound_1 = LBOUND(cell_to_index_3c, 1) 2214 bound_2 = UBOUND(cell_to_index_3c, 1) 2215 bound_3 = LBOUND(cell_to_index_3c, 2) 2216 bound_4 = UBOUND(cell_to_index_3c, 2) 2217 bound_5 = LBOUND(cell_to_index_3c, 3) 2218 bound_6 = UBOUND(cell_to_index_3c, 3) 2219 2220 maxcell = 8 2221 2222 DO x_cell_R1_plus_S2 = -maxcell, maxcell 2223 DO y_cell_R1_plus_S2 = -maxcell, maxcell 2224 DO z_cell_R1_plus_S2 = -maxcell, maxcell 2225 2226 already_there = .FALSE. 2227 2228 DO i_cell_S_1 = 1, num_cells_dm 2229 2230 IF (already_there) CYCLE 2231 2232 x_cell_R1_plus_S2_minus_S1 = x_cell_R1_plus_S2 - index_to_cell_dm(1, i_cell_S_1) 2233 y_cell_R1_plus_S2_minus_S1 = y_cell_R1_plus_S2 - index_to_cell_dm(2, i_cell_S_1) 2234 z_cell_R1_plus_S2_minus_S1 = z_cell_R1_plus_S2 - index_to_cell_dm(3, i_cell_S_1) 2235 2236 IF (x_cell_R1_plus_S2_minus_S1 .GE. bound_1 .AND. & 2237 x_cell_R1_plus_S2_minus_S1 .LE. bound_2 .AND. & 2238 y_cell_R1_plus_S2_minus_S1 .GE. bound_3 .AND. & 2239 y_cell_R1_plus_S2_minus_S1 .LE. bound_4 .AND. & 2240 z_cell_R1_plus_S2_minus_S1 .GE. bound_5 .AND. & 2241 z_cell_R1_plus_S2_minus_S1 .LE. bound_6) THEN 2242 2243 size_set = SIZE(index_to_cell_R1_plus_S2, 2) 2244 2245 IF (size_set == 1 .AND. & 2246 index_to_cell_R1_plus_S2(1, size_set) == 0 .AND. & 2247 index_to_cell_R1_plus_S2(2, size_set) == 0 .AND. & 2248 index_to_cell_R1_plus_S2(3, size_set) == 0) THEN 2249 2250 index_to_cell_R1_plus_S2(1, size_set) = x_cell_R1_plus_S2 2251 index_to_cell_R1_plus_S2(2, size_set) = y_cell_R1_plus_S2 2252 index_to_cell_R1_plus_S2(3, size_set) = z_cell_R1_plus_S2 2253 2254 ELSE 2255 2256 ALLOCATE (index_to_cell_tmp(3, size_set)) 2257 index_to_cell_tmp(1:3, 1:size_set) = index_to_cell_R1_plus_S2(1:3, 1:size_set) 2258 DEALLOCATE (index_to_cell_R1_plus_S2) 2259 ALLOCATE (index_to_cell_R1_plus_S2(3, size_set + 1)) 2260 index_to_cell_R1_plus_S2(1:3, 1:size_set) = index_to_cell_tmp(1:3, 1:size_set) 2261 index_to_cell_R1_plus_S2(1, size_set + 1) = x_cell_R1_plus_S2 2262 index_to_cell_R1_plus_S2(2, size_set + 1) = y_cell_R1_plus_S2 2263 index_to_cell_R1_plus_S2(3, size_set + 1) = z_cell_R1_plus_S2 2264 DEALLOCATE (index_to_cell_tmp) 2265 already_there = .TRUE. 2266 2267 END IF 2268 2269 END IF 2270 2271 END DO 2272 2273 END DO 2274 END DO 2275 END DO 2276 2277 num_cells_R1_plus_S2 = SIZE(index_to_cell_R1_plus_S2, 2) 2278 2279 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 2280 "GW_INFO| Number of periodic images for R_1+S_2 and R_2 sum in self-energy:", num_cells_R1_plus_S2 2281 2282 CALL timestop(handle) 2283 2284 END SUBROUTINE 2285 2286! ************************************************************************************************** 2287!> \brief ... 2288!> \param index_to_cell_R2 ... 2289!> \param cell_to_index_R2 ... 2290!> \param num_cells_R2 ... 2291!> \param unit_nr ... 2292!> \param index_to_cell_dm ... 2293!> \param qs_env ... 2294! ************************************************************************************************** 2295 SUBROUTINE compute_cell_vec_for_R2(index_to_cell_R2, cell_to_index_R2, num_cells_R2, unit_nr, & 2296 index_to_cell_dm, qs_env) 2297 2298 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_R2 2299 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_R2 2300 INTEGER :: num_cells_R2, unit_nr 2301 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 2302 TYPE(qs_environment_type), POINTER :: qs_env 2303 2304 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cell_vec_for_R2', & 2305 routineP = moduleN//':'//routineN 2306 2307 INTEGER :: bound_1, bound_2, bound_3, bound_4, & 2308 bound_5, bound_6, handle, icell, & 2309 num_cells_dm 2310 TYPE(kpoint_type), POINTER :: kpoints 2311 2312 CALL timeset(routineN, handle) 2313 2314 CALL get_qs_env(qs_env, kpoints=kpoints) 2315 2316 index_to_cell_dm => kpoints%index_to_cell 2317 2318 num_cells_dm = SIZE(index_to_cell_dm, 2) 2319 2320 ALLOCATE (index_to_cell_R2(3, num_cells_dm)) 2321 index_to_cell_R2(1:3, 1:num_cells_dm) = index_to_cell_dm(1:3, 1:num_cells_dm) 2322 2323 num_cells_R2 = SIZE(index_to_cell_R2, 2) 2324 2325 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 2326 "GW_INFO| Number of periodic images for R_2 W-sum in self-energy:", num_cells_R2 2327 2328 bound_1 = MINVAL(index_to_cell_R2(1, :)) 2329 bound_2 = MAXVAL(index_to_cell_R2(1, :)) 2330 bound_3 = MINVAL(index_to_cell_R2(2, :)) 2331 bound_4 = MAXVAL(index_to_cell_R2(2, :)) 2332 bound_5 = MINVAL(index_to_cell_R2(3, :)) 2333 bound_6 = MAXVAL(index_to_cell_R2(3, :)) 2334 2335 ALLOCATE (cell_to_index_R2(bound_1:bound_2, bound_3:bound_4, bound_5:bound_6)) 2336 cell_to_index_R2(:, :, :) = 0 2337 2338 DO icell = 1, SIZE(index_to_cell_R2, 2) 2339 2340 cell_to_index_R2(index_to_cell_R2(1, icell), & 2341 index_to_cell_R2(2, icell), & 2342 index_to_cell_R2(3, icell)) = icell 2343 2344 END DO 2345 2346 CALL timestop(handle) 2347 2348 END SUBROUTINE 2349 2350! ************************************************************************************************** 2351!> \brief ... 2352!> \param mat_3c_overl_int ... 2353!> \param ikp ... 2354!> \param mat_3c_overl_int_gw_kp_re ... 2355!> \param mat_3c_overl_int_gw_kp_im ... 2356!> \param kpoints ... 2357!> \param para_env ... 2358!> \param para_env_sub ... 2359!> \param matrix_s ... 2360!> \param mat_dm_virt_local ... 2361!> \param gw_corr_lev_occ ... 2362!> \param gw_corr_lev_virt ... 2363!> \param homo ... 2364!> \param nmo ... 2365!> \param cut_RI ... 2366!> \param num_3c_repl ... 2367!> \param row_from_LLL ... 2368!> \param my_group_L_starts_im_time ... 2369!> \param my_group_L_sizes_im_time ... 2370!> \param has_3c_blocks_re ... 2371!> \param has_3c_blocks_im ... 2372! ************************************************************************************************** 2373 SUBROUTINE get_mat_3c_gw_kp(mat_3c_overl_int, ikp, & 2374 mat_3c_overl_int_gw_kp_re, mat_3c_overl_int_gw_kp_im, kpoints, & 2375 para_env, para_env_sub, matrix_s, mat_dm_virt_local, & 2376 gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, cut_RI, num_3c_repl, & 2377 row_from_LLL, my_group_L_starts_im_time, & 2378 my_group_L_sizes_im_time, has_3c_blocks_re, has_3c_blocks_im) 2379 2380 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int 2381 INTEGER :: ikp 2382 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_gw_kp_re, & 2383 mat_3c_overl_int_gw_kp_im 2384 TYPE(kpoint_type), POINTER :: kpoints 2385 TYPE(cp_para_env_type), POINTER :: para_env, para_env_sub 2386 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 2387 TYPE(dbcsr_p_type) :: mat_dm_virt_local 2388 INTEGER :: gw_corr_lev_occ, gw_corr_lev_virt, homo, & 2389 nmo, cut_RI, num_3c_repl 2390 INTEGER, DIMENSION(:) :: row_from_LLL, my_group_L_starts_im_time, & 2391 my_group_L_sizes_im_time 2392 LOGICAL, DIMENSION(:, :, :) :: has_3c_blocks_re, has_3c_blocks_im 2393 2394 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_3c_gw_kp', & 2395 routineP = moduleN//':'//routineN 2396 2397 INTEGER :: handle, i_cell, i_cut_RI, icol_global, & 2398 irow_global, j_cell, n_level_gw, nblks 2399 REAL(KIND=dp) :: minus_one 2400 TYPE(cp_fm_type), POINTER :: fm_mat_mo_coeff_gw_im, & 2401 fm_mat_mo_coeff_gw_re, & 2402 fm_mat_mo_coeff_im, fm_mat_mo_coeff_re 2403 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_3c_overl_int_gw_for_mult_tmp 2404 TYPE(dbcsr_type), POINTER :: mat_mo_coeff_gw_local_tmp, & 2405 mat_mo_coeff_gw_tmp 2406 2407 CALL timeset(routineN, handle) 2408 2409 fm_mat_mo_coeff_re => kpoints%kp_env(ikp)%kpoint_env%mos(1, 1)%mo_set%mo_coeff 2410 fm_mat_mo_coeff_im => kpoints%kp_env(ikp)%kpoint_env%mos(2, 1)%mo_set%mo_coeff 2411 2412 NULLIFY (fm_mat_mo_coeff_gw_re) 2413 CALL cp_fm_create(fm_mat_mo_coeff_gw_re, fm_mat_mo_coeff_re%matrix_struct) 2414 CALL cp_fm_set_all(fm_mat_mo_coeff_gw_re, 0.0_dp) 2415 2416 NULLIFY (fm_mat_mo_coeff_gw_im) 2417 CALL cp_fm_create(fm_mat_mo_coeff_gw_im, fm_mat_mo_coeff_im%matrix_struct) 2418 CALL cp_fm_set_all(fm_mat_mo_coeff_gw_im, 0.0_dp) 2419 2420 CALL cp_fm_to_fm(fm_mat_mo_coeff_re, fm_mat_mo_coeff_gw_re) 2421 CALL cp_fm_to_fm(fm_mat_mo_coeff_im, fm_mat_mo_coeff_gw_im) 2422 2423 ! set MO coeffs to zero for non-GW corrected levels 2424 DO irow_global = 1, nmo 2425 DO icol_global = 1, homo - gw_corr_lev_occ 2426 CALL cp_fm_set_element(fm_mat_mo_coeff_gw_re, irow_global, icol_global, 0.0_dp) 2427 CALL cp_fm_set_element(fm_mat_mo_coeff_gw_im, irow_global, icol_global, 0.0_dp) 2428 END DO 2429 DO icol_global = homo + gw_corr_lev_virt + 1, nmo 2430 CALL cp_fm_set_element(fm_mat_mo_coeff_gw_re, irow_global, icol_global, 0.0_dp) 2431 CALL cp_fm_set_element(fm_mat_mo_coeff_gw_im, irow_global, icol_global, 0.0_dp) 2432 END DO 2433 END DO 2434 2435 NULLIFY (mat_mo_coeff_gw_tmp) 2436 CALL dbcsr_init_p(mat_mo_coeff_gw_tmp) 2437 CALL dbcsr_create(matrix=mat_mo_coeff_gw_tmp, & 2438 template=matrix_s(1)%matrix, & 2439 matrix_type=dbcsr_type_no_symmetry) 2440 2441 NULLIFY (mat_mo_coeff_gw_local_tmp) 2442 CALL dbcsr_init_p(mat_mo_coeff_gw_local_tmp) 2443 CALL dbcsr_create(matrix=mat_mo_coeff_gw_local_tmp, & 2444 template=mat_dm_virt_local%matrix, & 2445 matrix_type=dbcsr_type_no_symmetry) 2446 2447 NULLIFY (mat_3c_overl_int_gw_for_mult_tmp) 2448 CALL dbcsr_allocate_matrix_set(mat_3c_overl_int_gw_for_mult_tmp, cut_RI) 2449 DO i_cut_RI = 1, cut_RI 2450 ALLOCATE (mat_3c_overl_int_gw_for_mult_tmp(i_cut_RI)%matrix) 2451 CALL dbcsr_create(matrix=mat_3c_overl_int_gw_for_mult_tmp(i_cut_RI)%matrix, & 2452 template=mat_3c_overl_int(i_cut_RI, 1, 1)%matrix) 2453 END DO 2454 2455 ! ****************************** 1. REAL PART OF (nk_R muS P0) ************************************* 2456 CALL copy_fm_to_dbcsr(fm_mat_mo_coeff_gw_re, & 2457 mat_mo_coeff_gw_tmp, & 2458 keep_sparsity=.FALSE.) 2459 2460 ! just remove the blocks which have been set to zero 2461 CALL dbcsr_filter(mat_mo_coeff_gw_tmp, 1.0E-20_dp) 2462 2463 CALL dbcsr_set(mat_mo_coeff_gw_local_tmp, 0.0_dp) 2464 2465 CALL replicate_mat_to_subgroup_simple(para_env, para_env_sub, mat_mo_coeff_gw_tmp, nmo, & 2466 mat_mo_coeff_gw_local_tmp) 2467 2468 DO i_cell = 1, num_3c_repl 2469 DO j_cell = 1, num_3c_repl 2470 2471 DO i_cut_RI = 1, cut_RI 2472 CALL dbcsr_multiply("T", "N", 1.0_dp, mat_mo_coeff_gw_local_tmp, & 2473 mat_3c_overl_int(i_cut_RI, i_cell, j_cell)%matrix, & 2474 0.0_dp, mat_3c_overl_int_gw_for_mult_tmp(i_cut_RI)%matrix) 2475 END DO 2476 2477 CALL fill_mat_3c_overl_int_gw(mat_3c_overl_int_gw_kp_re(:, i_cell, j_cell), & 2478 mat_3c_overl_int_gw_for_mult_tmp, row_from_LLL, & 2479 my_group_L_starts_im_time, my_group_L_sizes_im_time, cut_RI, & 2480 para_env, gw_corr_lev_occ, gw_corr_lev_virt, homo) 2481 2482 END DO 2483 END DO 2484 2485 ! ****************************** 2. IMAG PART OF (nk_R muS P0) ************************************* 2486 CALL copy_fm_to_dbcsr(fm_mat_mo_coeff_gw_im, & 2487 mat_mo_coeff_gw_tmp, & 2488 keep_sparsity=.FALSE.) 2489 2490 ! just remove the blocks which have been set to zero 2491 CALL dbcsr_filter(mat_mo_coeff_gw_tmp, 1.0E-20_dp) 2492 2493 CALL dbcsr_set(mat_mo_coeff_gw_local_tmp, 0.0_dp) 2494 2495 CALL replicate_mat_to_subgroup_simple(para_env, para_env_sub, mat_mo_coeff_gw_tmp, nmo, & 2496 mat_mo_coeff_gw_local_tmp) 2497 2498 DO i_cell = 1, num_3c_repl 2499 DO j_cell = 1, num_3c_repl 2500 2501 DO i_cut_RI = 1, cut_RI 2502 ! minus one because MO coeffs are Hermitian conjugate 2503 minus_one = -1.0_dp 2504 CALL dbcsr_multiply("T", "N", minus_one, mat_mo_coeff_gw_local_tmp, & 2505 mat_3c_overl_int(i_cut_RI, i_cell, j_cell)%matrix, & 2506 0.0_dp, mat_3c_overl_int_gw_for_mult_tmp(i_cut_RI)%matrix) 2507 END DO 2508 2509 CALL fill_mat_3c_overl_int_gw(mat_3c_overl_int_gw_kp_im(:, i_cell, j_cell), & 2510 mat_3c_overl_int_gw_for_mult_tmp, row_from_LLL, & 2511 my_group_L_starts_im_time, my_group_L_sizes_im_time, cut_RI, & 2512 para_env, gw_corr_lev_occ, gw_corr_lev_virt, homo) 2513 2514 DO n_level_gw = 1, gw_corr_lev_occ + gw_corr_lev_virt 2515 2516 nblks = dbcsr_get_num_blocks(mat_3c_overl_int_gw_kp_re(n_level_gw, i_cell, j_cell)%matrix) 2517 CALL mp_sum(nblks, para_env%group) 2518 IF (nblks == 0) has_3c_blocks_re(n_level_gw, i_cell, j_cell) = .FALSE. 2519 IF (nblks > 0) has_3c_blocks_re(n_level_gw, i_cell, j_cell) = .TRUE. 2520 2521 nblks = dbcsr_get_num_blocks(mat_3c_overl_int_gw_kp_im(n_level_gw, i_cell, j_cell)%matrix) 2522 CALL mp_sum(nblks, para_env%group) 2523 IF (nblks == 0) has_3c_blocks_im(n_level_gw, i_cell, j_cell) = .FALSE. 2524 IF (nblks > 0) has_3c_blocks_im(n_level_gw, i_cell, j_cell) = .TRUE. 2525 2526 END DO 2527 2528 END DO 2529 END DO 2530 2531 CALL dbcsr_deallocate_matrix_set(mat_3c_overl_int_gw_for_mult_tmp) 2532 CALL dbcsr_release_p(mat_mo_coeff_gw_tmp) 2533 CALL dbcsr_release_p(mat_mo_coeff_gw_local_tmp) 2534 CALL cp_fm_release(fm_mat_mo_coeff_gw_re) 2535 CALL cp_fm_release(fm_mat_mo_coeff_gw_im) 2536 2537 CALL timestop(handle) 2538 2539 END SUBROUTINE 2540 2541! ************************************************************************************************** 2542!> \brief ... 2543!> \param index_to_cell_R2 ... 2544!> \param mat_W_R ... 2545!> \param mat_3c_overl_int_gw_kp_re ... 2546!> \param cfm_mat_W_kp_tau ... 2547!> \param kpoints ... 2548!> \param RI_blk_sizes ... 2549!> \param fm_mat_L ... 2550!> \param dimen_RI ... 2551!> \param qs_env ... 2552!> \param start_jquad ... 2553!> \param wkp_W ... 2554! ************************************************************************************************** 2555 SUBROUTINE trafo_W_from_k_to_R(index_to_cell_R2, mat_W_R, mat_3c_overl_int_gw_kp_re, & 2556 cfm_mat_W_kp_tau, kpoints, & 2557 RI_blk_sizes, fm_mat_L, dimen_RI, qs_env, start_jquad, wkp_W) 2558 2559 INTEGER, DIMENSION(:, :) :: index_to_cell_R2 2560 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_W_R 2561 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: mat_3c_overl_int_gw_kp_re 2562 TYPE(cp_cfm_p_type), DIMENSION(:, :), POINTER :: cfm_mat_W_kp_tau 2563 TYPE(kpoint_type), POINTER :: kpoints 2564 INTEGER, DIMENSION(:), POINTER :: RI_blk_sizes 2565 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L 2566 INTEGER :: dimen_RI 2567 TYPE(qs_environment_type), POINTER :: qs_env 2568 INTEGER :: start_jquad 2569 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 2570 2571 CHARACTER(LEN=*), PARAMETER :: routineN = 'trafo_W_from_k_to_R', & 2572 routineP = moduleN//':'//routineN 2573 2574 INTEGER :: handle, icell, ik, jquad, nkp, & 2575 num_cells_R2, num_integ_points, xcell, & 2576 ycell, zcell 2577 REAL(KIND=dp) :: arg, check, check_2, coskl, sinkl 2578 REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp 2579 TYPE(cp_fm_type), POINTER :: fm_tmp_im, fm_tmp_re 2580 TYPE(dbcsr_type), POINTER :: mat_work_im, mat_work_re 2581 2582 CALL timeset(routineN, handle) 2583 2584 NULLIFY (xkp) 2585 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp) 2586 2587 num_integ_points = SIZE(cfm_mat_W_kp_tau, 2) 2588 2589 IF (qs_env%mp2_env%ri_g0w0%print_exx == gw_read_exx) THEN 2590 ! we do not do HFX here, so we do not need jquad = 0 what corresponds to HFX 2591 start_jquad = 1 2592 ELSE 2593 start_jquad = 0 2594 END IF 2595 2596 num_cells_R2 = SIZE(index_to_cell_R2, 2) 2597 2598 NULLIFY (mat_W_R) 2599 ALLOCATE (mat_W_R(num_cells_R2, start_jquad:num_integ_points)) 2600 DO jquad = start_jquad, num_integ_points 2601 DO icell = 1, num_cells_R2 2602 ALLOCATE (mat_W_R(icell, jquad)%matrix) 2603 CALL dbcsr_create(matrix=mat_W_R(icell, jquad)%matrix, & 2604 template=mat_3c_overl_int_gw_kp_re(1, 1, 1)%matrix, & 2605 matrix_type=dbcsr_type_no_symmetry, & 2606 row_blk_size=RI_blk_sizes, & 2607 col_blk_size=RI_blk_sizes) 2608 CALL dbcsr_set(mat_W_R(icell, jquad)%matrix, 0.0_dp) 2609 END DO 2610 END DO 2611 2612 NULLIFY (fm_tmp_re) 2613 CALL cp_fm_create(fm_tmp_re, cfm_mat_W_kp_tau(1, 1)%matrix%matrix_struct) 2614 NULLIFY (fm_tmp_im) 2615 CALL cp_fm_create(fm_tmp_im, cfm_mat_W_kp_tau(1, 1)%matrix%matrix_struct) 2616 2617 NULLIFY (mat_work_re) 2618 CALL dbcsr_init_p(mat_work_re) 2619 CALL dbcsr_create(matrix=mat_work_re, & 2620 template=mat_W_R(1, 1)%matrix, & 2621 matrix_type=dbcsr_type_no_symmetry) 2622 2623 NULLIFY (mat_work_im) 2624 CALL dbcsr_init_p(mat_work_im) 2625 CALL dbcsr_create(matrix=mat_work_im, & 2626 template=mat_W_R(1, 1)%matrix, & 2627 matrix_type=dbcsr_type_no_symmetry) 2628 2629 check = 0.0_dp 2630 check_2 = 0.0_dp 2631 2632 DO jquad = start_jquad, num_integ_points 2633 2634 DO ik = 1, nkp 2635 2636 IF (jquad == 0) THEN 2637 ! jquad = 0 corresponds to exact exchange self-energy 2638 ! V(ik) = L(ik)*L^H(ik) 2639 2640 CALL cp_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, & 2641 fm_mat_L(ik, 1)%matrix, fm_mat_L(ik, 1)%matrix, & 2642 0.0_dp, fm_tmp_re) 2643 CALL cp_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, & 2644 fm_mat_L(ik, 2)%matrix, fm_mat_L(ik, 2)%matrix, & 2645 1.0_dp, fm_tmp_re) 2646 CALL cp_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, -1.0_dp, & 2647 fm_mat_L(ik, 1)%matrix, fm_mat_L(ik, 2)%matrix, & 2648 0.0_dp, fm_tmp_im) 2649 CALL cp_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, & 2650 fm_mat_L(ik, 2)%matrix, fm_mat_L(ik, 1)%matrix, & 2651 1.0_dp, fm_tmp_im) 2652 2653 ELSE 2654 2655 CALL cp_cfm_to_fm(cfm_mat_W_kp_tau(ik, jquad)%matrix, fm_tmp_re, fm_tmp_im) 2656 2657 END IF 2658 2659 CALL copy_fm_to_dbcsr(fm_tmp_re, mat_work_re, keep_sparsity=.FALSE.) 2660 CALL copy_fm_to_dbcsr(fm_tmp_im, mat_work_im, keep_sparsity=.FALSE.) 2661 2662 DO icell = 1, num_cells_R2 2663 2664 xcell = index_to_cell_R2(1, icell) 2665 ycell = index_to_cell_R2(2, icell) 2666 zcell = index_to_cell_R2(3, icell) 2667 2668 arg = REAL(xcell, dp)*xkp(1, ik) + REAL(ycell, dp)*xkp(2, ik) + REAL(zcell, dp)*xkp(3, ik) 2669 coskl = wkp_W(ik)*COS(twopi*arg) 2670 sinkl = wkp_W(ik)*SIN(twopi*arg) 2671 2672 CALL dbcsr_add(mat_W_R(icell, jquad)%matrix, mat_work_re, 1.0_dp, coskl) 2673 CALL dbcsr_add(mat_W_R(icell, jquad)%matrix, mat_work_im, 1.0_dp, sinkl) 2674 2675 END DO ! icell 2676 2677 END DO ! ik 2678 2679 END DO ! jquad 2680 2681 CALL cp_fm_release(fm_tmp_re) 2682 CALL cp_fm_release(fm_tmp_im) 2683 CALL dbcsr_release_p(mat_work_re) 2684 CALL dbcsr_release_p(mat_work_im) 2685 2686 CALL timestop(handle) 2687 2688 END SUBROUTINE 2689 2690END MODULE rpa_gw_kpoints 2691