1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines 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_type 27 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr 28 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add 29 USE cp_fm_types, ONLY: cp_fm_copy_general,& 30 cp_fm_create,& 31 cp_fm_p_type,& 32 cp_fm_release,& 33 cp_fm_set_all,& 34 cp_fm_type 35 USE cp_para_types, ONLY: cp_para_env_type 36 USE dbcsr_api, ONLY: & 37 dbcsr_add, dbcsr_create, dbcsr_dot, dbcsr_get_block_p, dbcsr_get_num_blocks, & 38 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, & 39 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_release, dbcsr_set, & 40 dbcsr_type 41 USE kinds, ONLY: dp 42 USE kpoint_types, ONLY: get_kpoint_info,& 43 kpoint_type 44 USE mathconstants, ONLY: gaussi,& 45 twopi,& 46 z_one,& 47 z_zero 48 USE mathlib, ONLY: invmat 49 USE message_passing, ONLY: mp_sum 50 USE particle_methods, ONLY: get_particle_set 51 USE particle_types, ONLY: particle_type 52 USE qs_environment_types, ONLY: get_qs_env,& 53 qs_environment_type 54 USE qs_integral_utils, ONLY: basis_set_list_setup 55 USE qs_kind_types, ONLY: qs_kind_type 56 USE rpa_im_time, ONLY: compute_transl_dm 57#include "./base/base_uses.f90" 58 59 IMPLICIT NONE 60 61 PRIVATE 62 63 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_kpoints' 64 65 PUBLIC :: compute_Wc_real_space_tau_GW, compute_Wc_kp_tau_GW, & 66 compute_wkp_W 67 68CONTAINS 69 70! ************************************************************************************************** 71!> \brief ... 72!> \param fm_mat_W_tau ... 73!> \param cfm_mat_Q ... 74!> \param fm_mat_L_re ... 75!> \param fm_mat_L_im ... 76!> \param dimen_RI ... 77!> \param num_integ_points ... 78!> \param jquad ... 79!> \param ikp ... 80!> \param tj ... 81!> \param tau_tj ... 82!> \param weights_cos_tf_w_to_t ... 83!> \param ikp_local ... 84!> \param para_env ... 85!> \param kpoints ... 86!> \param qs_env ... 87!> \param wkp_W ... 88!> \param mat_SinvVSinv ... 89!> \param do_W_and_not_V ... 90! ************************************************************************************************** 91 SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm_mat_L_im, & 92 dimen_RI, num_integ_points, jquad, & 93 ikp, tj, tau_tj, weights_cos_tf_w_to_t, ikp_local, & 94 para_env, kpoints, qs_env, wkp_W, mat_SinvVSinv, do_W_and_not_V) 95 96 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fm_mat_W_tau 97 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 98 TYPE(cp_fm_type), POINTER :: fm_mat_L_re, fm_mat_L_im 99 INTEGER :: dimen_RI, num_integ_points, jquad, ikp 100 REAL(KIND=dp), DIMENSION(:) :: tj 101 REAL(KIND=dp), DIMENSION(0:num_integ_points) :: tau_tj 102 REAL(KIND=dp), DIMENSION(:, :) :: weights_cos_tf_w_to_t 103 INTEGER, DIMENSION(:) :: ikp_local 104 TYPE(cp_para_env_type), POINTER :: para_env 105 TYPE(kpoint_type), POINTER :: kpoints 106 TYPE(qs_environment_type), POINTER :: qs_env 107 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 108 TYPE(dbcsr_p_type) :: mat_SinvVSinv 109 LOGICAL :: do_W_and_not_V 110 111 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Wc_real_space_tau_GW', & 112 routineP = moduleN//':'//routineN 113 114 INTEGER :: handle, handle2, i_global, iatom, iatom_old, icell, iiB, iquad, irow, j_global, & 115 jatom, jatom_old, jcol, jjB, jkp, LLL, natom, ncol_local, nkind, nkp, nrow_local, & 116 num_cells, xcell, ycell, zcell 117 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_RI_index 118 INTEGER, DIMENSION(:), POINTER :: col_indices, row_blk_end, row_blk_start, & 119 row_indices 120 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell 121 LOGICAL :: do_V_and_not_W 122 REAL(KIND=dp) :: abs_rab_cell, arg, contribution, coskl, cutoff_exp, d_0, omega, sinkl, & 123 sum_exp, sum_exp_k_im, sum_exp_k_re, tau, weight, weight_im, weight_re 124 REAL(KIND=dp), DIMENSION(3) :: cell_vector, rab_cell_i 125 REAL(KIND=dp), DIMENSION(3, 3) :: hmat 126 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp 127 REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp 128 TYPE(cell_type), POINTER :: cell 129 TYPE(cp_cfm_type), POINTER :: cfm_mat_L, cfm_mat_work, cfm_mat_work_2 130 TYPE(cp_fm_type), POINTER :: fm_dummy, fm_mat_work_global, & 131 fm_mat_work_local 132 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_RI_tmp 133 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 134 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 135 136 CALL timeset(routineN, handle) 137 138 CALL timeset(routineN//"_1", handle2) 139 140 NULLIFY (cfm_mat_work) 141 CALL cp_cfm_create(cfm_mat_work, cfm_mat_Q%matrix_struct) 142 CALL cp_cfm_set_all(cfm_mat_work, z_zero) 143 144 NULLIFY (cfm_mat_work_2) 145 CALL cp_cfm_create(cfm_mat_work_2, cfm_mat_Q%matrix_struct) 146 CALL cp_cfm_set_all(cfm_mat_work_2, z_zero) 147 148 NULLIFY (cfm_mat_L) 149 CALL cp_cfm_create(cfm_mat_L, cfm_mat_Q%matrix_struct) 150 CALL cp_cfm_set_all(cfm_mat_L, z_zero) 151 152 ! Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L 153 CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_L, z_one, fm_mat_L_re) 154 CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_L, gaussi, fm_mat_L_im) 155 156 NULLIFY (fm_mat_work_global) 157 CALL cp_fm_create(fm_mat_work_global, fm_mat_W_tau(1)%matrix%matrix_struct) 158 CALL cp_fm_set_all(fm_mat_work_global, 0.0_dp) 159 160 NULLIFY (fm_mat_work_local) 161 CALL cp_fm_create(fm_mat_work_local, cfm_mat_Q%matrix_struct) 162 CALL cp_fm_set_all(fm_mat_work_local, 0.0_dp) 163 164 CALL timestop(handle2) 165 166 IF (do_W_and_not_V) THEN 167 168 CALL timeset(routineN//"_2", handle2) 169 170 ! calculate [1+Q(iw')]^-1 171 CALL cp_cfm_cholesky_invert(cfm_mat_Q) 172 173 ! symmetrize the result 174 CALL own_cfm_upper_to_full(cfm_mat_Q, cfm_mat_work) 175 176 ! subtract exchange part by subtracing identity matrix from epsilon 177 CALL cp_cfm_get_info(matrix=cfm_mat_Q, & 178 nrow_local=nrow_local, & 179 ncol_local=ncol_local, & 180 row_indices=row_indices, & 181 col_indices=col_indices) 182 183 DO jjB = 1, ncol_local 184 j_global = col_indices(jjB) 185 DO iiB = 1, nrow_local 186 i_global = row_indices(iiB) 187 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 188 cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB) - z_one 189 END IF 190 END DO 191 END DO 192 193 CALL timestop(handle2) 194 195 CALL timeset(routineN//"_3.1", handle2) 196 197 ! work = epsilon(iw,k)*L^H(k) 198 CALL cp_cfm_gemm('N', 'C', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_Q, cfm_mat_L, & 199 z_zero, cfm_mat_work) 200 201 ! W(iw,k) = L(k)*work 202 CALL cp_cfm_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_work, & 203 z_zero, cfm_mat_work_2) 204 205 CALL timestop(handle2) 206 207 ELSE 208 209 ! S^-1(k)V(k)S^-1(k) = L(k)*L(k)^H 210 CALL cp_cfm_gemm('N', 'C', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_L, & 211 z_zero, cfm_mat_work_2) 212 213 END IF 214 215 CALL timeset(routineN//"_4", handle2) 216 217 CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp) 218 index_to_cell => kpoints%index_to_cell 219 num_cells = SIZE(index_to_cell, 2) 220 d_0 = qs_env%mp2_env%ri_rpa_im_time%cutoff 221 cutoff_exp = 10000.0_dp 222 CALL cp_cfm_set_all(cfm_mat_work, z_zero) 223 224 NULLIFY (qs_kind_set, cell, particle_set) 225 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell, natom=natom, nkind=nkind, & 226 particle_set=particle_set) 227 228 ALLOCATE (row_blk_start(natom)) 229 ALLOCATE (row_blk_end(natom)) 230 ALLOCATE (basis_set_RI_tmp(nkind)) 231 CALL basis_set_list_setup(basis_set_RI_tmp, "RI_AUX", qs_kind_set) 232 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, & 233 basis=basis_set_RI_tmp) 234 DEALLOCATE (basis_set_RI_tmp) 235 ALLOCATE (atom_from_RI_index(dimen_RI)) 236 DO LLL = 1, dimen_RI 237 DO iatom = 1, natom 238 IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN 239 atom_from_RI_index(LLL) = iatom 240 END IF 241 END DO 242 END DO 243 CALL get_cell(cell=cell, h=hmat) 244 iatom_old = 0 245 jatom_old = 0 246 247 CALL cp_cfm_get_info(matrix=cfm_mat_Q, & 248 nrow_local=nrow_local, & 249 ncol_local=ncol_local, & 250 row_indices=row_indices, & 251 col_indices=col_indices) 252 253 DO irow = 1, nrow_local 254 DO jcol = 1, ncol_local 255 256 iatom = atom_from_RI_index(row_indices(irow)) 257 jatom = atom_from_RI_index(col_indices(jcol)) 258 259 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN 260 261 sum_exp = 0.0_dp 262 sum_exp_k_re = 0.0_dp 263 sum_exp_k_im = 0.0_dp 264 265 DO icell = 1, num_cells 266 267 xcell = index_to_cell(1, icell) 268 ycell = index_to_cell(2, icell) 269 zcell = index_to_cell(3, icell) 270 271 arg = REAL(xcell, dp)*xkp(1, ikp) + REAL(ycell, dp)*xkp(2, ikp) + REAL(zcell, dp)*xkp(3, ikp) 272 273 coskl = wkp_W(ikp)*COS(twopi*arg) 274 sinkl = wkp_W(ikp)*SIN(twopi*arg) 275 276 cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, icell), dp)) 277 278 rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - & 279 (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3)) 280 281 abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2) 282 283 IF (abs_rab_cell/d_0 < cutoff_exp) THEN 284 sum_exp = sum_exp + EXP(-abs_rab_cell/d_0) 285 sum_exp_k_re = sum_exp_k_re + EXP(-abs_rab_cell/d_0)*coskl 286 sum_exp_k_im = sum_exp_k_im + EXP(-abs_rab_cell/d_0)*sinkl 287 END IF 288 289 END DO 290 291 weight_re = sum_exp_k_re/sum_exp 292 weight_im = sum_exp_k_im/sum_exp 293 294 iatom_old = iatom 295 jatom_old = jatom 296 297 END IF 298 299 contribution = weight_re*REAL(cfm_mat_work_2%local_data(irow, jcol)) + & 300 weight_im*AIMAG(cfm_mat_work_2%local_data(irow, jcol)) 301 302 fm_mat_work_local%local_data(irow, jcol) = fm_mat_work_local%local_data(irow, jcol) + contribution 303 304 END DO 305 END DO 306 307 CALL timestop(handle2) 308 309 CALL timeset(routineN//"_5", handle2) 310 311 IF (do_W_and_not_V) THEN 312 313 IF (SUM(ikp_local) > nkp) THEN 314 315 CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env) 316 317 DO iquad = 1, num_integ_points 318 319 omega = tj(jquad) 320 tau = tau_tj(iquad) 321 weight = weights_cos_tf_w_to_t(iquad, jquad)*COS(tau*omega) 322 323 IF (jquad == 1 .AND. ikp == 1) THEN 324 CALL cp_fm_set_all(matrix=fm_mat_W_tau(iquad)%matrix, alpha=0.0_dp) 325 END IF 326 327 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) 328 329 END DO 330 331 ELSE 332 333 DO jkp = 1, nkp 334 335 IF (ANY(ikp_local(:) == jkp)) THEN 336 CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env) 337 ELSE 338 NULLIFY (fm_dummy) 339 CALL cp_fm_copy_general(fm_dummy, fm_mat_work_global, para_env) 340 END IF 341 342 DO iquad = 1, num_integ_points 343 344 omega = tj(jquad) 345 tau = tau_tj(iquad) 346 weight = weights_cos_tf_w_to_t(iquad, jquad)*COS(tau*omega) 347 348 IF (jquad == 1 .AND. jkp == 1) THEN 349 CALL cp_fm_set_all(matrix=fm_mat_W_tau(iquad)%matrix, alpha=0.0_dp) 350 END IF 351 352 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_W_tau(iquad)%matrix, beta=weight, & 353 matrix_b=fm_mat_work_global) 354 355 END DO 356 357 END DO 358 359 END IF 360 361 END IF 362 363 do_V_and_not_W = .NOT. do_W_and_not_V 364 IF (do_V_and_not_W) THEN 365 366 IF (SUM(ikp_local) > nkp) THEN 367 CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env) 368 CALL fm_mat_work_global_to_mat_SinvVSinv(mat_SinvVSinv, fm_mat_work_global) 369 ELSE 370 DO jkp = 1, nkp 371 IF (ANY(ikp_local(:) == jkp)) THEN 372 CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env) 373 ELSE 374 NULLIFY (fm_dummy) 375 CALL cp_fm_copy_general(fm_dummy, fm_mat_work_global, para_env) 376 END IF 377 CALL fm_mat_work_global_to_mat_SinvVSinv(mat_SinvVSinv, fm_mat_work_global) 378 END DO 379 END IF 380 END IF 381 382 CALL cp_cfm_release(cfm_mat_work) 383 CALL cp_cfm_release(cfm_mat_work_2) 384 CALL cp_cfm_release(cfm_mat_L) 385 CALL cp_fm_release(fm_mat_work_global) 386 CALL cp_fm_release(fm_mat_work_local) 387 DEALLOCATE (atom_from_RI_index) 388 DEALLOCATE (row_blk_start) 389 DEALLOCATE (row_blk_end) 390 391 CALL timestop(handle2) 392 393 CALL timestop(handle) 394 395 END SUBROUTINE 396 397! ************************************************************************************************** 398!> \brief ... 399!> \param mat_SinvVSinv ... 400!> \param fm_mat_work_global ... 401! ************************************************************************************************** 402 SUBROUTINE fm_mat_work_global_to_mat_SinvVSinv(mat_SinvVSinv, fm_mat_work_global) 403 404 TYPE(dbcsr_p_type) :: mat_SinvVSinv 405 TYPE(cp_fm_type), POINTER :: fm_mat_work_global 406 407 CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_mat_work_global_to_mat_SinvVSinv', & 408 routineP = moduleN//':'//routineN 409 410 INTEGER :: handle 411 TYPE(dbcsr_p_type) :: mat_work 412 413 CALL timeset(routineN, handle) 414 415 NULLIFY (mat_work%matrix) 416 ALLOCATE (mat_work%matrix) 417 CALL dbcsr_create(mat_work%matrix, template=mat_SinvVSinv%matrix) 418 419 CALL copy_fm_to_dbcsr(fm_mat_work_global, mat_work%matrix, keep_sparsity=.FALSE.) 420 421 CALL dbcsr_add(mat_SinvVSinv%matrix, mat_work%matrix, 1.0_dp, 1.0_dp) 422 423 CALL dbcsr_release(mat_work%matrix) 424 DEALLOCATE (mat_work%matrix) 425 426 CALL timestop(handle) 427 428 END SUBROUTINE 429 430! ************************************************************************************************** 431!> \brief ... 432!> \param cfm_mat_W_kp_tau ... 433!> \param cfm_mat_Q ... 434!> \param fm_mat_L_re ... 435!> \param fm_mat_L_im ... 436!> \param dimen_RI ... 437!> \param num_integ_points ... 438!> \param jquad ... 439!> \param ikp ... 440!> \param tj ... 441!> \param tau_tj ... 442!> \param weights_cos_tf_w_to_t ... 443! ************************************************************************************************** 444 SUBROUTINE compute_Wc_kp_tau_GW(cfm_mat_W_kp_tau, cfm_mat_Q, fm_mat_L_re, fm_mat_L_im, & 445 dimen_RI, num_integ_points, jquad, & 446 ikp, tj, tau_tj, weights_cos_tf_w_to_t) 447 448 TYPE(cp_cfm_p_type), DIMENSION(:, :), POINTER :: cfm_mat_W_kp_tau 449 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q 450 TYPE(cp_fm_type), POINTER :: fm_mat_L_re, fm_mat_L_im 451 INTEGER :: dimen_RI, num_integ_points, jquad, ikp 452 REAL(KIND=dp), DIMENSION(:) :: tj 453 REAL(KIND=dp), DIMENSION(0:num_integ_points) :: tau_tj 454 REAL(KIND=dp), DIMENSION(:, :) :: weights_cos_tf_w_to_t 455 456 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Wc_kp_tau_GW', & 457 routineP = moduleN//':'//routineN 458 459 INTEGER :: handle, handle2, i_global, iiB, iquad, & 460 j_global, jjB, ncol_local, nrow_local 461 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 462 REAL(KIND=dp) :: omega, tau, weight 463 TYPE(cp_cfm_type), POINTER :: cfm_mat_L, cfm_mat_work 464 465 CALL timeset(routineN, handle) 466 467 NULLIFY (cfm_mat_work) 468 CALL cp_cfm_create(cfm_mat_work, fm_mat_L_re%matrix_struct) 469 CALL cp_cfm_set_all(cfm_mat_work, z_zero) 470 471 NULLIFY (cfm_mat_L) 472 CALL cp_cfm_create(cfm_mat_L, fm_mat_L_re%matrix_struct) 473 CALL cp_cfm_set_all(cfm_mat_L, z_zero) 474 475 CALL timeset(routineN//"_cholesky_inv", handle2) 476 477 ! calculate [1+Q(iw')]^-1 478 CALL cp_cfm_cholesky_invert(cfm_mat_Q) 479 480 ! symmetrize the result 481 CALL own_cfm_upper_to_full(cfm_mat_Q, cfm_mat_work) 482 483 ! subtract exchange part by subtracing identity matrix from epsilon 484 CALL cp_cfm_get_info(matrix=cfm_mat_Q, & 485 nrow_local=nrow_local, & 486 ncol_local=ncol_local, & 487 row_indices=row_indices, & 488 col_indices=col_indices) 489 490 DO jjB = 1, ncol_local 491 j_global = col_indices(jjB) 492 DO iiB = 1, nrow_local 493 i_global = row_indices(iiB) 494 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 495 cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB) - z_one 496 END IF 497 END DO 498 END DO 499 500 CALL timestop(handle2) 501 502 ! Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L 503 CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_L, z_one, fm_mat_L_re) 504 CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_L, gaussi, fm_mat_L_im) 505 506 ! work = epsilon(iw,k)*L^H(k) 507 CALL cp_cfm_gemm('N', 'C', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_Q, cfm_mat_L, & 508 z_zero, cfm_mat_work) 509 510 ! W(iw,k) = L(k)*work 511 CALL cp_cfm_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_work, & 512 z_zero, cfm_mat_Q) 513 514 DO iquad = 1, num_integ_points 515 omega = tj(jquad) 516 tau = tau_tj(iquad) 517 weight = weights_cos_tf_w_to_t(iquad, jquad)*COS(tau*omega) 518 CALL cp_cfm_scale_and_add(alpha=z_one, matrix_a=cfm_mat_W_kp_tau(ikp, iquad)%matrix, & 519 beta=CMPLX(weight, KIND=dp), matrix_b=cfm_mat_Q) 520 END DO 521 522 CALL cp_cfm_release(cfm_mat_work) 523 CALL cp_cfm_release(cfm_mat_L) 524 525 CALL timestop(handle) 526 527 END SUBROUTINE 528 529! ************************************************************************************************** 530!> \brief ... 531!> \param wkp_W ... 532!> \param kpoints ... 533!> \param h_mat ... 534!> \param h_inv ... 535!> \param exp_kpoints ... 536!> \param periodic ... 537! ************************************************************************************************** 538 SUBROUTINE compute_wkp_W(wkp_W, kpoints, h_mat, h_inv, exp_kpoints, periodic) 539 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_W 540 TYPE(kpoint_type), POINTER :: kpoints 541 REAL(KIND=dp), DIMENSION(3, 3) :: h_mat, h_inv 542 REAL(KIND=dp) :: exp_kpoints 543 INTEGER, DIMENSION(3) :: periodic 544 545 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_wkp_W', routineP = moduleN//':'//routineN 546 547 INTEGER :: handle, i_dim, i_x, ikp, info, j_y, k_z, & 548 n_x, n_y, n_z, nkp, nsuperfine, & 549 num_lin_eqs 550 REAL(KIND=dp) :: a_vec_dot_k_vec, integral, k_sq, weight 551 REAL(KIND=dp), DIMENSION(3) :: a_vec, k_vec, x_vec 552 REAL(KIND=dp), DIMENSION(:), POINTER :: right_side, wkp 553 REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix_lin_eqs, xkp 554 555 CALL timeset(routineN, handle) 556 557 CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp) 558 559 ! we determine the kpoint weights of the Monkhors Pack mesh new 560 ! such that the functions 1/k^2, 1/k and const are integrated exactly 561 ! in the Brillouin zone 562 ! this is done by minimizing sum_i |w_i|^2 where w_i are the weights of 563 ! the i-th kpoint under the following constraints: 564 ! 1) 1/k^2, 1/k and const are integrated exactly 565 ! 2) the kpoint weights of kpoints with identical absolute value are 566 ! the same, of e.g. (1/8,3/8,3/8) same weight as for (3/8,1/8,3/8) 567 ! for 1d and 2d materials: we use normal Monkhorst-Pack mesh, checked 568 ! by SUM(periodic) == 3 569 570 IF (exp_kpoints < 2.0_dp .AND. SUM(periodic) == 3) THEN 571 572 ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid 573 nsuperfine = 500 574 integral = 0.0_dp 575 IF (exp_kpoints > 0.0_dp) exp_kpoints = -2.0_dp 576 577 ! actually, there is the factor *det_3x3(h_inv) missing to account for the 578 ! integration volume but for wkp det_3x3(h_inv) is needed 579 weight = 2.0_dp/(REAL(nsuperfine, dp))**3 580 DO i_x = 1, nsuperfine 581 DO j_y = 1, nsuperfine 582 DO k_z = 1, nsuperfine/2 583 584 x_vec = (/REAL(i_x - nsuperfine/2, dp) - 0.5_dp, & 585 REAL(j_y - nsuperfine/2, dp) - 0.5_dp, & 586 REAL(k_z - nsuperfine/2, dp) - 0.5_dp/)/ & 587 REAL(nsuperfine, dp) 588 k_vec = MATMUL(h_inv(1:3, 1:3), x_vec) 589 k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 590 integral = integral + weight*k_sq**(exp_kpoints*0.5_dp) 591 END DO 592 END DO 593 END DO 594 595 num_lin_eqs = nkp + 2 596 597 ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs)) 598 matrix_lin_eqs(:, :) = 0.0_dp 599 600 DO ikp = 1, nkp 601 602 k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp)) 603 k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 604 605 matrix_lin_eqs(ikp, ikp) = 2.0_dp 606 matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp 607 matrix_lin_eqs(ikp, nkp + 2) = 1.0_dp*k_sq**(exp_kpoints*0.5_dp) 608 609 matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp 610 matrix_lin_eqs(nkp + 2, ikp) = 1.0_dp*k_sq**(exp_kpoints*0.5_dp) 611 612 END DO 613 614 CALL invmat(matrix_lin_eqs, info) 615 ! check whether inversion was successful 616 CPASSERT(info == 0) 617 618 ALLOCATE (wkp_W(num_lin_eqs)) 619 620 ALLOCATE (right_side(num_lin_eqs)) 621 right_side = 0.0_dp 622 right_side(nkp + 1) = 1.0_dp 623 right_side(nkp + 2) = integral 624 625 wkp_W(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side) 626 627 DEALLOCATE (matrix_lin_eqs, right_side) 628 629 ELSE IF (exp_kpoints < 2.0_dp .AND. SUM(periodic) == 1) THEN 630 631 ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid 632 nsuperfine = 5000 633 integral = 0.0_dp 634 635 ! actually, there is the factor *det_3x3(h_inv) missing to account for the 636 ! integration volume but for wkp det_3x3(h_inv) is needed 637 weight = 1.0_dp/REAL(nsuperfine, dp) 638 IF (periodic(1) == 1) THEN 639 n_x = nsuperfine 640 ELSE 641 n_x = 1 642 END IF 643 IF (periodic(2) == 1) THEN 644 n_y = nsuperfine 645 ELSE 646 n_y = 1 647 END IF 648 IF (periodic(3) == 1) THEN 649 n_z = nsuperfine 650 ELSE 651 n_z = 1 652 END IF 653 654 a_vec = MATMUL(h_mat(1:3, 1:3), & 655 (/REAL(periodic(1), dp), REAL(periodic(2), dp), REAL(periodic(3), dp)/)) 656 657 DO i_x = 1, n_x 658 DO j_y = 1, n_y 659 DO k_z = 1, n_z 660 661 x_vec = (/REAL(i_x - nsuperfine/2, dp) - 0.5_dp, & 662 REAL(j_y - nsuperfine/2, dp) - 0.5_dp, & 663 REAL(k_z - nsuperfine/2, dp) - 0.5_dp/)/ & 664 REAL(nsuperfine, dp) 665 666 DO i_dim = 1, 3 667 IF (periodic(i_dim) == 0) THEN 668 x_vec(i_dim) = 0.0_dp 669 END IF 670 END DO 671 672 k_vec = MATMUL(h_inv(1:3, 1:3), x_vec) 673 a_vec_dot_k_vec = a_vec(1)*k_vec(1) + a_vec(2)*k_vec(2) + a_vec(3)*k_vec(3) 674 integral = integral + weight*LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) 675 END DO 676 END DO 677 END DO 678 679 num_lin_eqs = nkp + 2 680 681 ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs)) 682 matrix_lin_eqs(:, :) = 0.0_dp 683 684 DO ikp = 1, nkp 685 686 k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp)) 687 k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 688 689 matrix_lin_eqs(ikp, ikp) = 2.0_dp 690 matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp 691 692 a_vec_dot_k_vec = a_vec(1)*k_vec(1) + a_vec(2)*k_vec(2) + a_vec(3)*k_vec(3) 693 matrix_lin_eqs(ikp, nkp + 2) = LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) 694 695 matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp 696 matrix_lin_eqs(nkp + 2, ikp) = LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) 697 698 END DO 699 700 CALL invmat(matrix_lin_eqs, info) 701 ! check whether inversion was successful 702 CPASSERT(info == 0) 703 704 ALLOCATE (wkp_W(num_lin_eqs)) 705 706 ALLOCATE (right_side(num_lin_eqs)) 707 right_side = 0.0_dp 708 right_side(nkp + 1) = 1.0_dp 709 right_side(nkp + 2) = integral 710 711 wkp_W(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side) 712 713 DEALLOCATE (matrix_lin_eqs, right_side) 714 715 ELSE 716 717 ALLOCATE (wkp_W(nkp)) 718 wkp_W(:) = wkp(:) 719 720 END IF 721 722 CALL timestop(handle) 723 724 END SUBROUTINE 725 726! ************************************************************************************************** 727!> \brief ... 728!> \param cfm_mat_Q ... 729!> \param cfm_mat_work ... 730! ************************************************************************************************** 731 SUBROUTINE own_cfm_upper_to_full(cfm_mat_Q, cfm_mat_work) 732 733 TYPE(cp_cfm_type), POINTER :: cfm_mat_Q, cfm_mat_work 734 735 CHARACTER(LEN=*), PARAMETER :: routineN = 'own_cfm_upper_to_full', & 736 routineP = moduleN//':'//routineN 737 738 INTEGER :: handle, i_global, iiB, j_global, jjB, & 739 ncol_local, nrow_local 740 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 741 742 CALL timeset(routineN, handle) 743 744 ! get info of fm_mat_Q 745 CALL cp_cfm_get_info(matrix=cfm_mat_Q, & 746 nrow_local=nrow_local, & 747 ncol_local=ncol_local, & 748 row_indices=row_indices, & 749 col_indices=col_indices) 750 751 DO jjB = 1, ncol_local 752 j_global = col_indices(jjB) 753 DO iiB = 1, nrow_local 754 i_global = row_indices(iiB) 755 IF (j_global < i_global) THEN 756 cfm_mat_Q%local_data(iiB, jjB) = z_zero 757 END IF 758 IF (j_global == i_global) THEN 759 cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB)/(2.0_dp, 0.0_dp) 760 END IF 761 END DO 762 END DO 763 764 CALL cp_cfm_transpose(cfm_mat_Q, 'C', cfm_mat_work) 765 766 CALL cp_cfm_scale_and_add(z_one, cfm_mat_Q, z_one, cfm_mat_work) 767 768 CALL timestop(handle) 769 770 END SUBROUTINE 771 772! ************************************************************************************************** 773!> \brief ... 774!> \param mat_contr_W_re ... 775!> \param mat_contr_W_im ... 776!> \param n_level_gw ... 777!> \param jquad ... 778!> \param gw_corr_lev_occ ... 779!> \param homo ... 780!> \param mat_contr_gf_occ_re ... 781!> \param mat_contr_gf_occ_im ... 782!> \param mat_contr_gf_virt_re ... 783!> \param mat_contr_gf_virt_im ... 784!> \param vec_Sigma_c_gw_pos_tau ... 785!> \param vec_Sigma_c_gw_neg_tau ... 786!> \param vec_Sigma_x_gw ... 787!> \param ikp ... 788!> \param cycle_R1_S2_n_level ... 789!> \param cycle_R1_plus_S2_n_level ... 790!> \param eps_filter ... 791!> \param i_cell_R1_plus_S2 ... 792!> \param i_cell_S2 ... 793!> \param start_jquad ... 794! ************************************************************************************************** 795 SUBROUTINE trace_for_self_ener(mat_contr_W_re, mat_contr_W_im, n_level_gw, jquad, & 796 gw_corr_lev_occ, homo, & 797 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 798 mat_contr_gf_virt_re, mat_contr_gf_virt_im, & 799 vec_Sigma_c_gw_pos_tau, vec_Sigma_c_gw_neg_tau, & 800 vec_Sigma_x_gw, ikp, & 801 cycle_R1_S2_n_level, cycle_R1_plus_S2_n_level, & 802 eps_filter, i_cell_R1_plus_S2, i_cell_S2, & 803 start_jquad) 804 805 TYPE(dbcsr_type), POINTER :: mat_contr_W_re, mat_contr_W_im 806 INTEGER :: n_level_gw, jquad, gw_corr_lev_occ, homo 807 TYPE(dbcsr_type), POINTER :: mat_contr_gf_occ_re, & 808 mat_contr_gf_occ_im, & 809 mat_contr_gf_virt_re, & 810 mat_contr_gf_virt_im 811 COMPLEX(KIND=dp), DIMENSION(:, :) :: vec_Sigma_c_gw_pos_tau, & 812 vec_Sigma_c_gw_neg_tau 813 REAL(KIND=dp), DIMENSION(:) :: vec_Sigma_x_gw 814 INTEGER :: ikp 815 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :) :: cycle_R1_S2_n_level 816 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: cycle_R1_plus_S2_n_level 817 REAL(KIND=dp) :: eps_filter 818 INTEGER :: i_cell_R1_plus_S2, i_cell_S2, start_jquad 819 820 CHARACTER(LEN=*), PARAMETER :: routineN = 'trace_for_self_ener', & 821 routineP = moduleN//':'//routineN 822 823 INTEGER :: handle, n_level_gw_ref 824 REAL(KIND=dp) :: trace_neg_tau_im_1, trace_neg_tau_im_2, trace_neg_tau_re_1, & 825 trace_neg_tau_re_2, trace_pos_tau_im_1, trace_pos_tau_im_2, trace_pos_tau_re_1, & 826 trace_pos_tau_re_2 827 828 CALL timeset(routineN, handle) 829 830 CALL dbcsr_dot(mat_contr_gf_occ_re, & 831 mat_contr_W_re, & 832 trace_neg_tau_re_1) 833 834 CALL dbcsr_dot(mat_contr_gf_occ_im, & 835 mat_contr_W_im, & 836 trace_neg_tau_re_2) 837 838 CALL dbcsr_dot(mat_contr_gf_occ_re, & 839 mat_contr_W_im, & 840 trace_neg_tau_im_1) 841 842 CALL dbcsr_dot(mat_contr_gf_occ_im, & 843 mat_contr_W_re, & 844 trace_neg_tau_im_2) 845 846 IF (ABS(trace_neg_tau_re_1) + ABS(trace_neg_tau_re_2) + ABS(trace_neg_tau_im_1) + & 847 ABS(trace_neg_tau_im_2) < eps_filter) THEN 848 cycle_R1_S2_n_level(i_cell_R1_plus_S2, i_cell_S2, n_level_gw, jquad) = .TRUE. 849 END IF 850 851 IF (ikp == 1 .AND. jquad == start_jquad .AND. i_cell_S2 == 1) THEN 852 cycle_R1_plus_S2_n_level(i_cell_R1_plus_S2, n_level_gw, jquad) = .TRUE. 853 END IF 854 855 IF (ABS(trace_neg_tau_re_1) + ABS(trace_neg_tau_re_2) + ABS(trace_neg_tau_im_1) + & 856 ABS(trace_neg_tau_im_2) > eps_filter) THEN 857 cycle_R1_plus_S2_n_level(i_cell_R1_plus_S2, n_level_gw, jquad) = .FALSE. 858 END IF 859 860 IF (jquad == 0) THEN 861 862 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 863 864 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 865 866 ELSE 867 868 vec_Sigma_c_gw_neg_tau(n_level_gw, jquad) = vec_Sigma_c_gw_neg_tau(n_level_gw, jquad) + & 869 CMPLX(trace_neg_tau_re_1 - trace_neg_tau_re_2, & 870 trace_neg_tau_im_1 + trace_neg_tau_im_2, dp) 871 872 END IF 873 874 CALL dbcsr_dot(mat_contr_gf_virt_re, & 875 mat_contr_W_re, & 876 trace_pos_tau_re_1) 877 878 CALL dbcsr_dot(mat_contr_gf_virt_im, & 879 mat_contr_W_im, & 880 trace_pos_tau_re_2) 881 882 CALL dbcsr_dot(mat_contr_gf_virt_re, & 883 mat_contr_W_im, & 884 trace_pos_tau_im_1) 885 886 CALL dbcsr_dot(mat_contr_gf_virt_im, & 887 mat_contr_W_re, & 888 trace_pos_tau_im_2) 889 890 IF (jquad > 0) THEN 891 892 vec_Sigma_c_gw_pos_tau(n_level_gw, jquad) = vec_Sigma_c_gw_pos_tau(n_level_gw, jquad) + & 893 CMPLX(trace_pos_tau_re_1 - trace_pos_tau_re_2, & 894 trace_pos_tau_im_1 + trace_pos_tau_im_2, dp) 895 896 END IF 897 898 CALL timestop(handle) 899 900 END SUBROUTINE 901 902! ************************************************************************************************** 903!> \brief ... 904!> \param mat_I_muP_occ_re ... 905!> \param mat_I_muP_virt_re ... 906!> \param mat_I_muP_occ_im ... 907!> \param mat_I_muP_virt_im ... 908!> \param mat_contr_gf_occ_re ... 909!> \param mat_contr_gf_occ_im ... 910!> \param mat_contr_gf_virt_re ... 911!> \param mat_contr_gf_virt_im ... 912!> \param index_to_cell_3c ... 913!> \param kpoints ... 914!> \param ikp ... 915!> \param x_cell_R1 ... 916!> \param y_cell_R1 ... 917!> \param z_cell_R1 ... 918! ************************************************************************************************** 919 SUBROUTINE trafo_I_T_R1_plus_S2_to_M_R1_S2(mat_I_muP_occ_re, mat_I_muP_virt_re, & 920 mat_I_muP_occ_im, mat_I_muP_virt_im, & 921 mat_contr_gf_occ_re, mat_contr_gf_occ_im, & 922 mat_contr_gf_virt_re, mat_contr_gf_virt_im, & 923 index_to_cell_3c, kpoints, ikp, & 924 x_cell_R1, y_cell_R1, z_cell_R1) 925 926 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_I_muP_occ_re, mat_I_muP_virt_re, & 927 mat_I_muP_occ_im, mat_I_muP_virt_im 928 TYPE(dbcsr_type), POINTER :: mat_contr_gf_occ_re, & 929 mat_contr_gf_occ_im, & 930 mat_contr_gf_virt_re, & 931 mat_contr_gf_virt_im 932 INTEGER, DIMENSION(:, :) :: index_to_cell_3c 933 TYPE(kpoint_type), POINTER :: kpoints 934 INTEGER :: ikp, x_cell_R1, y_cell_R1, z_cell_R1 935 936 CHARACTER(LEN=*), PARAMETER :: routineN = 'trafo_I_T_R1_plus_S2_to_M_R1_S2', & 937 routineP = moduleN//':'//routineN 938 939 INTEGER :: handle, i_cell_T, num_cells_3c, xcell, & 940 ycell, zcell 941 REAL(KIND=dp) :: arg, coskl, sinkl 942 REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp 943 944 CALL timeset(routineN, handle) 945 946 num_cells_3c = SIZE(mat_I_muP_occ_re) 947 CALL get_kpoint_info(kpoints, xkp=xkp) 948 949 CALL dbcsr_set(mat_contr_gf_occ_re, 0.0_dp) 950 CALL dbcsr_set(mat_contr_gf_occ_im, 0.0_dp) 951 CALL dbcsr_set(mat_contr_gf_virt_re, 0.0_dp) 952 CALL dbcsr_set(mat_contr_gf_virt_im, 0.0_dp) 953 954 DO i_cell_T = 1, num_cells_3c 955 956 xcell = index_to_cell_3c(1, i_cell_T) - x_cell_R1 957 ycell = index_to_cell_3c(2, i_cell_T) - y_cell_R1 958 zcell = index_to_cell_3c(3, i_cell_T) - z_cell_R1 959 arg = REAL(xcell, dp)*xkp(1, ikp) + REAL(ycell, dp)*xkp(2, ikp) + REAL(zcell, dp)*xkp(3, ikp) 960 coskl = COS(twopi*arg) 961 sinkl = SIN(twopi*arg) 962 963 CALL dbcsr_scale_and_add_local(mat_contr_gf_occ_re, mat_I_muP_occ_re(i_cell_T)%matrix, coskl) 964 CALL dbcsr_scale_and_add_local(mat_contr_gf_occ_re, mat_I_muP_occ_im(i_cell_T)%matrix, -sinkl) 965 CALL dbcsr_scale_and_add_local(mat_contr_gf_occ_im, mat_I_muP_occ_re(i_cell_T)%matrix, sinkl) 966 CALL dbcsr_scale_and_add_local(mat_contr_gf_occ_im, mat_I_muP_occ_im(i_cell_T)%matrix, coskl) 967 968 CALL dbcsr_scale_and_add_local(mat_contr_gf_virt_re, mat_I_muP_virt_re(i_cell_T)%matrix, coskl) 969 CALL dbcsr_scale_and_add_local(mat_contr_gf_virt_re, mat_I_muP_virt_im(i_cell_T)%matrix, -sinkl) 970 CALL dbcsr_scale_and_add_local(mat_contr_gf_virt_im, mat_I_muP_virt_re(i_cell_T)%matrix, sinkl) 971 CALL dbcsr_scale_and_add_local(mat_contr_gf_virt_im, mat_I_muP_virt_im(i_cell_T)%matrix, coskl) 972 973 END DO 974 975 CALL timestop(handle) 976 977 END SUBROUTINE 978 979! ************************************************************************************************** 980!> \brief ... 981!> \param mat_A ... 982!> \param mat_B ... 983!> \param beta ... 984! ************************************************************************************************** 985 SUBROUTINE dbcsr_scale_and_add_local(mat_A, mat_B, beta) 986 TYPE(dbcsr_type), POINTER :: mat_A, mat_B 987 REAL(KIND=dp) :: beta 988 989 INTEGER :: col, row 990 LOGICAL :: found 991 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_to_compute, data_block 992 TYPE(dbcsr_iterator_type) :: iter 993 994 CALL dbcsr_iterator_start(iter, mat_B) 995 DO WHILE (dbcsr_iterator_blocks_left(iter)) 996 997 CALL dbcsr_iterator_next_block(iter, row, col, data_block) 998 999 NULLIFY (block_to_compute) 1000 CALL dbcsr_get_block_p(matrix=mat_A, & 1001 row=row, col=col, block=block_to_compute, found=found) 1002 1003 CPASSERT(found) 1004 1005 block_to_compute(:, :) = block_to_compute(:, :) + beta*data_block(:, :) 1006 1007 END DO 1008 1009 CALL dbcsr_iterator_stop(iter) 1010 1011 END SUBROUTINE 1012 1013! ************************************************************************************************** 1014!> \brief ... 1015!> \param mat_dm_occ_S ... 1016!> \param mat_dm_virt_S ... 1017!> \param qs_env ... 1018!> \param num_integ_points ... 1019!> \param e_fermi ... 1020!> \param eps_filter ... 1021!> \param tau_tj ... 1022!> \param num_cells_dm ... 1023!> \param index_to_cell_dm ... 1024!> \param has_blocks_mat_dm_occ_S ... 1025!> \param has_blocks_mat_dm_virt_S ... 1026!> \param para_env ... 1027! ************************************************************************************************** 1028 SUBROUTINE compute_G_real_space(mat_dm_occ_S, mat_dm_virt_S, qs_env, num_integ_points, & 1029 e_fermi, eps_filter, tau_tj, num_cells_dm, index_to_cell_dm, & 1030 has_blocks_mat_dm_occ_S, has_blocks_mat_dm_virt_S, para_env) 1031 1032 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_S, mat_dm_virt_S 1033 TYPE(qs_environment_type), POINTER :: qs_env 1034 INTEGER :: num_integ_points 1035 REAL(KIND=dp) :: e_fermi, eps_filter 1036 REAL(KIND=dp), DIMENSION(0:num_integ_points) :: tau_tj 1037 INTEGER :: num_cells_dm 1038 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 1039 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: has_blocks_mat_dm_occ_S, & 1040 has_blocks_mat_dm_virt_S 1041 TYPE(cp_para_env_type), POINTER :: para_env 1042 1043 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_G_real_space', & 1044 routineP = moduleN//':'//routineN 1045 1046 INTEGER :: handle, i_cell, ispin, jquad, nblks 1047 REAL(KIND=dp) :: tau 1048 1049 CALL timeset(routineN, handle) 1050 1051 ispin = 1 1052 1053 ! get denity matrix for exchange self-energy 1054 tau = 0.0_dp 1055 CALL compute_transl_dm(mat_dm_occ_S, qs_env, ispin, num_integ_points, 0, e_fermi, tau, & 1056 eps_filter, num_cells_dm, index_to_cell_dm, & 1057 remove_occ=.FALSE., remove_virt=.TRUE., first_jquad=0) 1058 1059 DO i_cell = 1, num_cells_dm 1060 1061 nblks = dbcsr_get_num_blocks(mat_dm_occ_S(0, i_cell)%matrix) 1062 CALL mp_sum(nblks, para_env%group) 1063 IF (nblks == 0) has_blocks_mat_dm_occ_S(0, i_cell) = .FALSE. 1064 IF (nblks > 0) has_blocks_mat_dm_occ_S(0, i_cell) = .TRUE. 1065 1066 END DO 1067 1068 DO jquad = 1, num_integ_points 1069 1070 tau = tau_tj(jquad) 1071 1072 CALL compute_transl_dm(mat_dm_occ_S, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & 1073 eps_filter, num_cells_dm, index_to_cell_dm, & 1074 remove_occ=.FALSE., remove_virt=.TRUE., first_jquad=0) 1075 1076 CALL compute_transl_dm(mat_dm_virt_S, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & 1077 eps_filter, num_cells_dm, index_to_cell_dm, & 1078 remove_occ=.TRUE., remove_virt=.FALSE., first_jquad=1) 1079 1080 DO i_cell = 1, num_cells_dm 1081 1082 nblks = dbcsr_get_num_blocks(mat_dm_occ_S(jquad, i_cell)%matrix) 1083 CALL mp_sum(nblks, para_env%group) 1084 IF (nblks == 0) has_blocks_mat_dm_occ_S(jquad, i_cell) = .FALSE. 1085 IF (nblks > 0) has_blocks_mat_dm_occ_S(jquad, i_cell) = .TRUE. 1086 1087 nblks = dbcsr_get_num_blocks(mat_dm_virt_S(jquad, i_cell)%matrix) 1088 CALL mp_sum(nblks, para_env%group) 1089 IF (nblks == 0) has_blocks_mat_dm_virt_S(jquad, i_cell) = .FALSE. 1090 IF (nblks > 0) has_blocks_mat_dm_virt_S(jquad, i_cell) = .TRUE. 1091 1092 END DO 1093 1094 END DO 1095 1096 CALL timestop(handle) 1097 1098 END SUBROUTINE 1099 1100! ************************************************************************************************** 1101!> \brief ... 1102!> \param index_to_cell_R1_plus_S2 ... 1103!> \param cell_to_index_3c ... 1104!> \param num_cells_R1_plus_S2 ... 1105!> \param kpoints ... 1106!> \param unit_nr ... 1107!> \param maxcell ... 1108!> \param num_cells_dm ... 1109! ************************************************************************************************** 1110 SUBROUTINE compute_cell_vec_for_R1_plus_S2_or_R1(index_to_cell_R1_plus_S2, cell_to_index_3c, & 1111 num_cells_R1_plus_S2, kpoints, & 1112 unit_nr, maxcell, num_cells_dm) 1113 1114 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_R1_plus_S2 1115 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c 1116 INTEGER :: num_cells_R1_plus_S2 1117 TYPE(kpoint_type), POINTER :: kpoints 1118 INTEGER :: unit_nr, maxcell, num_cells_dm 1119 1120 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cell_vec_for_R1_plus_S2_or_R1', & 1121 routineP = moduleN//':'//routineN 1122 1123 INTEGER :: bound_1, bound_2, bound_3, bound_4, bound_5, bound_6, handle, i_cell_S_1, & 1124 size_set, x_cell_R1_plus_S2, x_cell_R1_plus_S2_minus_S1, y_cell_R1_plus_S2, & 1125 y_cell_R1_plus_S2_minus_S1, z_cell_R1_plus_S2, z_cell_R1_plus_S2_minus_S1 1126 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_tmp 1127 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 1128 LOGICAL :: already_there 1129 1130 CALL timeset(routineN, handle) 1131 1132 index_to_cell_dm => kpoints%index_to_cell 1133 1134 ALLOCATE (index_to_cell_R1_plus_S2(3, 1)) 1135 index_to_cell_R1_plus_S2 = 0 1136 1137 bound_1 = LBOUND(cell_to_index_3c, 1) 1138 bound_2 = UBOUND(cell_to_index_3c, 1) 1139 bound_3 = LBOUND(cell_to_index_3c, 2) 1140 bound_4 = UBOUND(cell_to_index_3c, 2) 1141 bound_5 = LBOUND(cell_to_index_3c, 3) 1142 bound_6 = UBOUND(cell_to_index_3c, 3) 1143 1144 maxcell = 8 1145 1146 DO x_cell_R1_plus_S2 = -maxcell, maxcell 1147 DO y_cell_R1_plus_S2 = -maxcell, maxcell 1148 DO z_cell_R1_plus_S2 = -maxcell, maxcell 1149 1150 already_there = .FALSE. 1151 1152 DO i_cell_S_1 = 1, num_cells_dm 1153 1154 IF (already_there) CYCLE 1155 1156 x_cell_R1_plus_S2_minus_S1 = x_cell_R1_plus_S2 - index_to_cell_dm(1, i_cell_S_1) 1157 y_cell_R1_plus_S2_minus_S1 = y_cell_R1_plus_S2 - index_to_cell_dm(2, i_cell_S_1) 1158 z_cell_R1_plus_S2_minus_S1 = z_cell_R1_plus_S2 - index_to_cell_dm(3, i_cell_S_1) 1159 1160 IF (x_cell_R1_plus_S2_minus_S1 .GE. bound_1 .AND. & 1161 x_cell_R1_plus_S2_minus_S1 .LE. bound_2 .AND. & 1162 y_cell_R1_plus_S2_minus_S1 .GE. bound_3 .AND. & 1163 y_cell_R1_plus_S2_minus_S1 .LE. bound_4 .AND. & 1164 z_cell_R1_plus_S2_minus_S1 .GE. bound_5 .AND. & 1165 z_cell_R1_plus_S2_minus_S1 .LE. bound_6) THEN 1166 1167 size_set = SIZE(index_to_cell_R1_plus_S2, 2) 1168 1169 IF (size_set == 1 .AND. & 1170 index_to_cell_R1_plus_S2(1, size_set) == 0 .AND. & 1171 index_to_cell_R1_plus_S2(2, size_set) == 0 .AND. & 1172 index_to_cell_R1_plus_S2(3, size_set) == 0) THEN 1173 1174 index_to_cell_R1_plus_S2(1, size_set) = x_cell_R1_plus_S2 1175 index_to_cell_R1_plus_S2(2, size_set) = y_cell_R1_plus_S2 1176 index_to_cell_R1_plus_S2(3, size_set) = z_cell_R1_plus_S2 1177 1178 ELSE 1179 1180 ALLOCATE (index_to_cell_tmp(3, size_set)) 1181 index_to_cell_tmp(1:3, 1:size_set) = index_to_cell_R1_plus_S2(1:3, 1:size_set) 1182 DEALLOCATE (index_to_cell_R1_plus_S2) 1183 ALLOCATE (index_to_cell_R1_plus_S2(3, size_set + 1)) 1184 index_to_cell_R1_plus_S2(1:3, 1:size_set) = index_to_cell_tmp(1:3, 1:size_set) 1185 index_to_cell_R1_plus_S2(1, size_set + 1) = x_cell_R1_plus_S2 1186 index_to_cell_R1_plus_S2(2, size_set + 1) = y_cell_R1_plus_S2 1187 index_to_cell_R1_plus_S2(3, size_set + 1) = z_cell_R1_plus_S2 1188 DEALLOCATE (index_to_cell_tmp) 1189 already_there = .TRUE. 1190 1191 END IF 1192 1193 END IF 1194 1195 END DO 1196 1197 END DO 1198 END DO 1199 END DO 1200 1201 num_cells_R1_plus_S2 = SIZE(index_to_cell_R1_plus_S2, 2) 1202 1203 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 1204 "GW_INFO| Number of periodic images for R_1+S_2 and R_2 sum in self-energy:", num_cells_R1_plus_S2 1205 1206 CALL timestop(handle) 1207 1208 END SUBROUTINE 1209 1210! ************************************************************************************************** 1211!> \brief ... 1212!> \param index_to_cell_R2 ... 1213!> \param cell_to_index_R2 ... 1214!> \param num_cells_R2 ... 1215!> \param unit_nr ... 1216!> \param index_to_cell_dm ... 1217!> \param qs_env ... 1218! ************************************************************************************************** 1219 SUBROUTINE compute_cell_vec_for_R2(index_to_cell_R2, cell_to_index_R2, num_cells_R2, unit_nr, & 1220 index_to_cell_dm, qs_env) 1221 1222 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_R2 1223 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_R2 1224 INTEGER :: num_cells_R2, unit_nr 1225 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 1226 TYPE(qs_environment_type), POINTER :: qs_env 1227 1228 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cell_vec_for_R2', & 1229 routineP = moduleN//':'//routineN 1230 1231 INTEGER :: bound_1, bound_2, bound_3, bound_4, & 1232 bound_5, bound_6, handle, icell, & 1233 num_cells_dm 1234 TYPE(kpoint_type), POINTER :: kpoints 1235 1236 CALL timeset(routineN, handle) 1237 1238 CALL get_qs_env(qs_env, kpoints=kpoints) 1239 1240 index_to_cell_dm => kpoints%index_to_cell 1241 1242 num_cells_dm = SIZE(index_to_cell_dm, 2) 1243 1244 ALLOCATE (index_to_cell_R2(3, num_cells_dm)) 1245 index_to_cell_R2(1:3, 1:num_cells_dm) = index_to_cell_dm(1:3, 1:num_cells_dm) 1246 1247 num_cells_R2 = SIZE(index_to_cell_R2, 2) 1248 1249 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 1250 "GW_INFO| Number of periodic images for R_2 W-sum in self-energy:", num_cells_R2 1251 1252 bound_1 = MINVAL(index_to_cell_R2(1, :)) 1253 bound_2 = MAXVAL(index_to_cell_R2(1, :)) 1254 bound_3 = MINVAL(index_to_cell_R2(2, :)) 1255 bound_4 = MAXVAL(index_to_cell_R2(2, :)) 1256 bound_5 = MINVAL(index_to_cell_R2(3, :)) 1257 bound_6 = MAXVAL(index_to_cell_R2(3, :)) 1258 1259 ALLOCATE (cell_to_index_R2(bound_1:bound_2, bound_3:bound_4, bound_5:bound_6)) 1260 cell_to_index_R2(:, :, :) = 0 1261 1262 DO icell = 1, SIZE(index_to_cell_R2, 2) 1263 1264 cell_to_index_R2(index_to_cell_R2(1, icell), & 1265 index_to_cell_R2(2, icell), & 1266 index_to_cell_R2(3, icell)) = icell 1267 1268 END DO 1269 1270 CALL timestop(handle) 1271 1272 END SUBROUTINE 1273 1274END MODULE rpa_gw_kpoints 1275