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 GW 8!> \par History 9!> 03.2019 created [Frederick Stein] 10! ************************************************************************************************** 11MODULE rpa_gw 12 USE basis_set_types, ONLY: gto_basis_set_p_type,& 13 gto_basis_set_type 14 USE cell_types, ONLY: cell_type,& 15 get_cell 16 USE cp_cfm_types, ONLY: cp_cfm_p_type,& 17 cp_cfm_release 18 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 19 copy_fm_to_dbcsr,& 20 dbcsr_allocate_matrix_set,& 21 dbcsr_deallocate_matrix_set 22 USE cp_files, ONLY: close_file,& 23 open_file 24 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,& 25 cp_fm_upper_to_full 26 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& 27 cp_fm_cholesky_invert 28 USE cp_fm_diag, ONLY: cp_fm_syevd 29 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 30 cp_fm_struct_release,& 31 cp_fm_struct_type 32 USE cp_fm_types, ONLY: cp_fm_create,& 33 cp_fm_get_info,& 34 cp_fm_p_type,& 35 cp_fm_release,& 36 cp_fm_set_all,& 37 cp_fm_to_fm,& 38 cp_fm_type 39 USE cp_gemm_interface, ONLY: cp_gemm 40 USE cp_para_types, ONLY: cp_para_env_type 41 USE dbcsr_api, ONLY: & 42 dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_filter, & 43 dbcsr_get_info, dbcsr_init_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, & 44 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, & 45 dbcsr_p_type, dbcsr_release_p, dbcsr_scalar, dbcsr_scale, dbcsr_set, dbcsr_type, & 46 dbcsr_type_no_symmetry 47 USE dbcsr_tensor_api, ONLY: & 48 dbcsr_t_contract, dbcsr_t_copy, dbcsr_t_copy_matrix_to_tensor, dbcsr_t_create, & 49 dbcsr_t_destroy, dbcsr_t_get_block, dbcsr_t_get_info, dbcsr_t_iterator_blocks_left, & 50 dbcsr_t_iterator_next_block, dbcsr_t_iterator_start, dbcsr_t_iterator_stop, & 51 dbcsr_t_iterator_type, dbcsr_t_nblks_total, dbcsr_t_pgrid_create, dbcsr_t_pgrid_destroy, & 52 dbcsr_t_pgrid_type, dbcsr_t_type 53 USE input_constants, ONLY: gw_pade_approx,& 54 gw_two_pole_model,& 55 ri_rpa_g0w0_crossing_bisection,& 56 ri_rpa_g0w0_crossing_newton,& 57 ri_rpa_g0w0_crossing_z_shot 58 USE kinds, ONLY: default_path_length,& 59 dp 60 USE kpoint_types, ONLY: get_kpoint_info,& 61 kpoint_create,& 62 kpoint_release,& 63 kpoint_sym_create,& 64 kpoint_type 65 USE mathconstants, ONLY: fourpi,& 66 pi,& 67 twopi 68 USE message_passing, ONLY: mp_sum,& 69 mp_sync 70 USE mp2_types, ONLY: mp2_type 71 USE particle_types, ONLY: particle_type 72 USE physcon, ONLY: evolt 73 USE qs_band_structure, ONLY: calculate_kp_orbitals 74 USE qs_environment_types, ONLY: get_qs_env,& 75 qs_env_release,& 76 qs_environment_type 77 USE qs_gamma2kp, ONLY: create_kp_from_gamma 78 USE qs_integral_utils, ONLY: basis_set_list_setup 79 USE qs_kind_types, ONLY: get_qs_kind,& 80 qs_kind_type 81 USE qs_ks_types, ONLY: qs_ks_env_type 82 USE qs_mo_types, ONLY: get_mo_set,& 83 mo_set_type 84 USE qs_moments, ONLY: build_berry_moment_matrix 85 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,& 86 release_neighbor_list_sets 87 USE qs_neighbor_lists, ONLY: setup_neighbor_list 88 USE qs_overlap, ONLY: build_overlap_matrix_simple 89 USE qs_tensors_types, ONLY: create_2c_tensor 90 USE rpa_gw_im_time_util, ONLY: get_tensor_3c_overl_int_gw 91#include "./base/base_uses.f90" 92 93 IMPLICIT NONE 94 95 PRIVATE 96 97 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw' 98 99 PUBLIC :: allocate_matrices_gw_im_time, allocate_matrices_gw, GW_matrix_operations, compute_QP_energies, & 100 deallocate_matrices_gw_im_time, deallocate_matrices_gw 101 102CONTAINS 103 104! ************************************************************************************************** 105!> \brief ... 106!> \param gw_corr_lev_occ ... 107!> \param gw_corr_lev_virt ... 108!> \param homo ... 109!> \param nmo ... 110!> \param num_integ_points ... 111!> \param unit_nr ... 112!> \param RI_blk_sizes ... 113!> \param do_ic_model ... 114!> \param para_env ... 115!> \param fm_mat_W ... 116!> \param fm_mat_Q ... 117!> \param mo_coeff ... 118!> \param t_3c_overl_int_gw_RI ... 119!> \param t_3c_overl_int_gw_AO ... 120!> \param t_3c_overl_nnP_ic ... 121!> \param t_3c_overl_nnP_ic_reflected ... 122!> \param matrix_s ... 123!> \param mat_W ... 124!> \param t_3c_overl_int ... 125!> \param qs_env ... 126!> \param gw_corr_lev_occ_beta ... 127!> \param gw_corr_lev_virt_beta ... 128!> \param homo_beta ... 129!> \param mo_coeff_beta ... 130!> \param t_3c_overl_int_gw_RI_beta ... 131!> \param t_3c_overl_int_gw_AO_beta ... 132!> \param t_3c_overl_nnP_ic_beta ... 133!> \param t_3c_overl_nnP_ic_reflected_beta ... 134! ************************************************************************************************** 135 SUBROUTINE allocate_matrices_gw_im_time(gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, & 136 num_integ_points, unit_nr, & 137 RI_blk_sizes, do_ic_model, & 138 para_env, fm_mat_W, fm_mat_Q, & 139 mo_coeff, & 140 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, & 141 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, & 142 matrix_s, mat_W, t_3c_overl_int, & 143 qs_env, & 144 gw_corr_lev_occ_beta, gw_corr_lev_virt_beta, homo_beta, & 145 mo_coeff_beta, & 146 t_3c_overl_int_gw_RI_beta, t_3c_overl_int_gw_AO_beta, & 147 t_3c_overl_nnP_ic_beta, t_3c_overl_nnP_ic_reflected_beta) 148 149 INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo, & 150 nmo, num_integ_points, unit_nr 151 INTEGER, DIMENSION(:), POINTER :: RI_blk_sizes 152 LOGICAL, INTENT(IN) :: do_ic_model 153 TYPE(cp_para_env_type), POINTER :: para_env 154 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fm_mat_W 155 TYPE(cp_fm_type), POINTER :: fm_mat_Q, mo_coeff 156 TYPE(dbcsr_t_type) :: t_3c_overl_int_gw_RI, & 157 t_3c_overl_int_gw_AO, & 158 t_3c_overl_nnP_ic, & 159 t_3c_overl_nnP_ic_reflected 160 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 161 TYPE(dbcsr_type), POINTER :: mat_W 162 TYPE(dbcsr_t_type), DIMENSION(:, :) :: t_3c_overl_int 163 TYPE(qs_environment_type), POINTER :: qs_env 164 INTEGER, INTENT(IN), OPTIONAL :: gw_corr_lev_occ_beta, & 165 gw_corr_lev_virt_beta, homo_beta 166 TYPE(cp_fm_type), OPTIONAL, POINTER :: mo_coeff_beta 167 TYPE(dbcsr_t_type), OPTIONAL :: t_3c_overl_int_gw_RI_beta, t_3c_overl_int_gw_AO_beta, & 168 t_3c_overl_nnP_ic_beta, t_3c_overl_nnP_ic_reflected_beta 169 170 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_matrices_gw_im_time', & 171 routineP = moduleN//':'//routineN 172 173 INTEGER :: handle, jquad 174 LOGICAL :: my_open_shell 175 176 CALL timeset(routineN, handle) 177 178 my_open_shell = .FALSE. 179 IF (PRESENT(gw_corr_lev_occ_beta) .AND. PRESENT(gw_corr_lev_virt_beta) .AND. PRESENT(homo_beta) .AND. & 180 PRESENT(mo_coeff_beta) .AND. PRESENT(t_3c_overl_int_gw_RI_beta) .AND. PRESENT(t_3c_overl_int_gw_AO_beta) .AND. & 181 PRESENT(t_3c_overl_nnP_ic_beta) & 182 .AND. PRESENT(t_3c_overl_nnP_ic_reflected_beta)) THEN 183 my_open_shell = .TRUE. 184 185 END IF 186 187 CALL get_tensor_3c_overl_int_gw(t_3c_overl_int, & 188 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, & 189 mo_coeff, matrix_s, & 190 gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, & 191 para_env, & 192 do_ic_model, & 193 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, & 194 qs_env, unit_nr) 195 196 IF (my_open_shell) THEN 197 198 CALL get_tensor_3c_overl_int_gw(t_3c_overl_int, & 199 t_3c_overl_int_gw_RI_beta, t_3c_overl_int_gw_AO_beta, & 200 mo_coeff_beta, matrix_s, & 201 gw_corr_lev_occ_beta, gw_corr_lev_virt_beta, homo_beta, nmo, & 202 para_env, & 203 do_ic_model, & 204 t_3c_overl_nnP_ic_beta, t_3c_overl_nnP_ic_reflected_beta, & 205 qs_env, unit_nr) 206 207 END IF 208 209 NULLIFY (fm_mat_W) 210 ALLOCATE (fm_mat_W(num_integ_points)) 211 212 DO jquad = 1, num_integ_points 213 214 NULLIFY (fm_mat_W(jquad)%matrix) 215 CALL cp_fm_create(fm_mat_W(jquad)%matrix, fm_mat_Q%matrix_struct) 216 CALL cp_fm_to_fm(fm_mat_Q, fm_mat_W(jquad)%matrix) 217 CALL cp_fm_set_all(fm_mat_W(jquad)%matrix, 0.0_dp) 218 219 END DO 220 221 NULLIFY (mat_W) 222 CALL dbcsr_init_p(mat_W) 223 CALL dbcsr_create(matrix=mat_W, & 224 template=matrix_s(1)%matrix, & 225 matrix_type=dbcsr_type_no_symmetry, & 226 row_blk_size=RI_blk_sizes, & 227 col_blk_size=RI_blk_sizes) 228 229 CALL timestop(handle) 230 231 END SUBROUTINE allocate_matrices_gw_im_time 232 233! ************************************************************************************************** 234!> \brief ... 235!> \param vec_Sigma_c_gw ... 236!> \param color_rpa_group ... 237!> \param dimen_nm_gw ... 238!> \param gw_corr_lev_occ ... 239!> \param gw_corr_lev_virt ... 240!> \param homo ... 241!> \param nmo ... 242!> \param num_integ_group ... 243!> \param num_integ_points ... 244!> \param unit_nr ... 245!> \param gw_corr_lev_tot ... 246!> \param num_fit_points ... 247!> \param omega_max_fit ... 248!> \param do_minimax_quad ... 249!> \param do_periodic ... 250!> \param do_ri_Sigma_x ... 251!> \param my_do_gw ... 252!> \param first_cycle_periodic_correction ... 253!> \param a_scaling ... 254!> \param Eigenval ... 255!> \param tj ... 256!> \param vec_omega_fit_gw ... 257!> \param vec_Sigma_x_gw ... 258!> \param delta_corr ... 259!> \param Eigenval_last ... 260!> \param Eigenval_scf ... 261!> \param vec_W_gw ... 262!> \param fm_mat_S_gw ... 263!> \param fm_mat_S_gw_work ... 264!> \param para_env ... 265!> \param mp2_env ... 266!> \param kpoints ... 267!> \param nkp ... 268!> \param nkp_self_energy ... 269!> \param do_kpoints_cubic_RPA ... 270!> \param do_kpoints_from_Gamma ... 271!> \param vec_Sigma_c_gw_beta ... 272!> \param gw_corr_lev_occ_beta ... 273!> \param homo_beta ... 274!> \param Eigenval_beta ... 275!> \param vec_Sigma_x_gw_beta ... 276!> \param Eigenval_last_beta ... 277!> \param Eigenval_scf_beta ... 278!> \param vec_W_gw_beta ... 279!> \param fm_mat_S_gw_work_beta ... 280!> \param fm_mat_S_gw_beta ... 281! ************************************************************************************************** 282 SUBROUTINE allocate_matrices_gw(vec_Sigma_c_gw, color_rpa_group, dimen_nm_gw, & 283 gw_corr_lev_occ, gw_corr_lev_virt, homo, & 284 nmo, num_integ_group, num_integ_points, unit_nr, & 285 gw_corr_lev_tot, num_fit_points, omega_max_fit, & 286 do_minimax_quad, do_periodic, do_ri_Sigma_x, my_do_gw, & 287 first_cycle_periodic_correction, & 288 a_scaling, Eigenval, tj, vec_omega_fit_gw, vec_Sigma_x_gw, & 289 delta_corr, Eigenval_last, Eigenval_scf, vec_W_gw, & 290 fm_mat_S_gw, fm_mat_S_gw_work, & 291 para_env, mp2_env, kpoints, nkp, nkp_self_energy, & 292 do_kpoints_cubic_RPA, do_kpoints_from_Gamma, & 293 vec_Sigma_c_gw_beta, gw_corr_lev_occ_beta, homo_beta, Eigenval_beta, & 294 vec_Sigma_x_gw_beta, Eigenval_last_beta, Eigenval_scf_beta, vec_W_gw_beta, & 295 fm_mat_S_gw_work_beta, fm_mat_S_gw_beta) 296 297 COMPLEX(KIND=dp), ALLOCATABLE, & 298 DIMENSION(:, :, :), INTENT(OUT) :: vec_Sigma_c_gw 299 INTEGER, INTENT(IN) :: color_rpa_group, dimen_nm_gw, gw_corr_lev_occ, gw_corr_lev_virt, & 300 homo, nmo, num_integ_group, num_integ_points, unit_nr 301 INTEGER, INTENT(INOUT) :: gw_corr_lev_tot, num_fit_points 302 REAL(KIND=dp) :: omega_max_fit 303 LOGICAL, INTENT(IN) :: do_minimax_quad, do_periodic, & 304 do_ri_Sigma_x, my_do_gw 305 LOGICAL, INTENT(OUT) :: first_cycle_periodic_correction 306 REAL(KIND=dp), INTENT(IN) :: a_scaling 307 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: Eigenval 308 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 309 INTENT(IN) :: tj 310 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 311 INTENT(OUT) :: vec_omega_fit_gw 312 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 313 INTENT(OUT) :: vec_Sigma_x_gw 314 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 315 INTENT(INOUT) :: delta_corr, Eigenval_last, Eigenval_scf, & 316 vec_W_gw 317 TYPE(cp_fm_type), POINTER :: fm_mat_S_gw, fm_mat_S_gw_work 318 TYPE(cp_para_env_type), POINTER :: para_env 319 TYPE(mp2_type), POINTER :: mp2_env 320 TYPE(kpoint_type), POINTER :: kpoints 321 INTEGER, INTENT(OUT) :: nkp, nkp_self_energy 322 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA, & 323 do_kpoints_from_Gamma 324 COMPLEX(KIND=dp), ALLOCATABLE, & 325 DIMENSION(:, :, :), INTENT(OUT), OPTIONAL :: vec_Sigma_c_gw_beta 326 INTEGER, INTENT(IN), OPTIONAL :: gw_corr_lev_occ_beta, homo_beta 327 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), & 328 OPTIONAL :: Eigenval_beta 329 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 330 INTENT(OUT), OPTIONAL :: vec_Sigma_x_gw_beta 331 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 332 INTENT(INOUT), OPTIONAL :: Eigenval_last_beta, Eigenval_scf_beta, & 333 vec_W_gw_beta 334 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_mat_S_gw_work_beta, fm_mat_S_gw_beta 335 336 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_matrices_gw', & 337 routineP = moduleN//':'//routineN 338 339 INTEGER :: handle, iquad, jquad 340 LOGICAL :: my_open_shell 341 REAL(KIND=dp) :: omega 342 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vec_omega_gw 343 344 CALL timeset(routineN, handle) 345 346 my_open_shell = .FALSE. 347 IF (PRESENT(vec_Sigma_c_gw_beta) .AND. PRESENT(gw_corr_lev_occ_beta) .AND. PRESENT(homo_beta) .AND. & 348 PRESENT(Eigenval_beta) .AND. PRESENT(vec_Sigma_x_gw_beta) .AND. PRESENT(Eigenval_last_beta) .AND. & 349 PRESENT(Eigenval_scf_beta) .AND. PRESENT(vec_W_gw_beta) .AND. PRESENT(fm_mat_S_gw_work_beta) .AND. & 350 PRESENT(fm_mat_S_gw_beta)) THEN 351 my_open_shell = .TRUE. 352 END IF 353 354 gw_corr_lev_tot = gw_corr_lev_occ + gw_corr_lev_virt 355 356 ! fill the omega_frequency vector 357 ALLOCATE (vec_omega_gw(num_integ_points)) 358 vec_omega_gw = 0.0_dp 359 360 DO jquad = 1, num_integ_points 361 IF (do_minimax_quad) THEN 362 omega = tj(jquad) 363 ELSE 364 omega = a_scaling/TAN(tj(jquad)) 365 END IF 366 vec_omega_gw(jquad) = omega 367 END DO 368 369 ! determine number of fit points in the interval [0,w_max] for virt, or [-w_max,0] for occ 370 num_fit_points = 0 371 372 DO jquad = 1, num_integ_points 373 IF (vec_omega_gw(jquad) < omega_max_fit) THEN 374 num_fit_points = num_fit_points + 1 375 END IF 376 END DO 377 378 IF (mp2_env%ri_g0w0%analytic_continuation == gw_pade_approx) THEN 379 IF (mp2_env%ri_g0w0%nparam_pade > num_fit_points) THEN 380 IF (unit_nr > 0) & 381 CPWARN("Pade approximation: more parameters than data points. Reset # of parameters.") 382 mp2_env%ri_g0w0%nparam_pade = num_fit_points 383 IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T74,I7)") & 384 "Number of pade parameters:", mp2_env%ri_g0w0%nparam_pade 385 ENDIF 386 ENDIF 387 388 ! create new arrays containing omega values at which we calculate vec_Sigma_c_gw 389 ALLOCATE (vec_omega_fit_gw(num_fit_points)) 390 391 ! fill the omega vector with frequencies, where we calculate the self-energy 392 iquad = 0 393 DO jquad = 1, num_integ_points 394 IF (vec_omega_gw(jquad) < omega_max_fit) THEN 395 iquad = iquad + 1 396 vec_omega_fit_gw(iquad) = vec_omega_gw(jquad) 397 END IF 398 END DO 399 400 DEALLOCATE (vec_omega_gw) 401 402 IF (do_kpoints_cubic_RPA) THEN 403 CALL get_kpoint_info(kpoints, nkp=nkp) 404 IF (mp2_env%ri_g0w0%do_gamma_only_sigma) THEN 405 nkp_self_energy = 1 406 ELSE 407 nkp_self_energy = nkp 408 END IF 409 ELSE IF (do_kpoints_from_Gamma) THEN 410 CALL get_kpoint_info(kpoints, nkp=nkp) 411 nkp_self_energy = 1 412 ELSE 413 nkp = 1 414 nkp_self_energy = 1 415 END IF 416 ALLOCATE (vec_Sigma_c_gw(gw_corr_lev_tot, num_fit_points, nkp_self_energy)) 417 vec_Sigma_c_gw = (0.0_dp, 0.0_dp) 418 419 IF (my_open_shell) THEN 420 ALLOCATE (vec_Sigma_c_gw_beta(gw_corr_lev_tot, num_fit_points, nkp_self_energy)) 421 vec_Sigma_c_gw_beta = (0.0_dp, 0.0_dp) 422 END IF 423 424 ALLOCATE (Eigenval_scf(nmo)) 425 Eigenval_scf(:) = Eigenval(:) 426 427 ALLOCATE (Eigenval_last(nmo)) 428 Eigenval_last(:) = Eigenval(:) 429 430 ! Eigenval for beta 431 IF (my_open_shell) THEN 432 ALLOCATE (Eigenval_scf_beta(nmo)) 433 Eigenval_scf_beta(:) = Eigenval_beta(:) 434 435 ALLOCATE (Eigenval_last_beta(nmo)) 436 Eigenval_last_beta(:) = Eigenval_beta(:) 437 END IF 438 439 IF (do_periodic) THEN 440 441 ALLOCATE (delta_corr(1 + homo - gw_corr_lev_occ:homo + gw_corr_lev_virt)) 442 delta_corr(:) = 0.0_dp 443 444 first_cycle_periodic_correction = .TRUE. 445 446 END IF 447 448 IF (do_ri_Sigma_x) THEN 449 ALLOCATE (vec_Sigma_x_gw(nmo, nkp_self_energy)) 450 vec_Sigma_x_gw = 0.0_dp 451 452 IF (my_open_shell) THEN 453 ALLOCATE (vec_Sigma_x_gw_beta(nmo, nkp_self_energy)) 454 vec_Sigma_x_gw_beta = 0.0_dp 455 END IF 456 END IF 457 458 IF (my_do_gw) THEN 459 460 ! minimax grids not implemented for O(N^4) GW 461 CPASSERT(.NOT. do_minimax_quad) 462 463 ! create temporary matrix to store B*([1+Q(iw')]^-1-1), has the same size as B 464 NULLIFY (fm_mat_S_gw_work) 465 CALL cp_fm_create(fm_mat_S_gw_work, fm_mat_S_gw%matrix_struct) 466 CALL cp_fm_set_all(matrix=fm_mat_S_gw_work, alpha=0.0_dp) 467 468 IF (my_open_shell) THEN 469 NULLIFY (fm_mat_S_gw_work_beta) 470 CALL cp_fm_create(fm_mat_S_gw_work_beta, fm_mat_S_gw%matrix_struct) 471 CALL cp_fm_set_all(matrix=fm_mat_S_gw_work_beta, alpha=0.0_dp) 472 END IF 473 474 ALLOCATE (vec_W_gw(dimen_nm_gw)) 475 vec_W_gw = 0.0_dp 476 477 IF (my_open_shell) THEN 478 ALLOCATE (vec_W_gw_beta(dimen_nm_gw)) 479 vec_W_gw_beta = 0.0_dp 480 END IF 481 482 ! in case we do RI for Sigma_x, we calculate Sigma_x right here 483 IF (do_ri_Sigma_x) THEN 484 485 CALL get_vec_sigma_x(vec_Sigma_x_gw, nmo, fm_mat_S_gw, para_env, num_integ_group, color_rpa_group, & 486 homo, gw_corr_lev_occ, mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, 1)) 487 488 IF (my_open_shell) THEN 489 CALL get_vec_sigma_x(vec_Sigma_x_gw_beta, nmo, fm_mat_S_gw_beta, para_env, num_integ_group, & 490 color_rpa_group, homo_beta, gw_corr_lev_occ_beta, & 491 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, 1)) 492 END IF 493 494 END IF 495 496 END IF 497 498 CALL timestop(handle) 499 500 END SUBROUTINE allocate_matrices_gw 501 502! ************************************************************************************************** 503!> \brief ... 504!> \param vec_Sigma_x_gw ... 505!> \param nmo ... 506!> \param fm_mat_S_gw ... 507!> \param para_env ... 508!> \param num_integ_group ... 509!> \param color_rpa_group ... 510!> \param homo ... 511!> \param gw_corr_lev_occ ... 512!> \param vec_Sigma_x_minus_vxc_gw11 ... 513! ************************************************************************************************** 514 SUBROUTINE get_vec_sigma_x(vec_Sigma_x_gw, nmo, fm_mat_S_gw, para_env, num_integ_group, color_rpa_group, homo, & 515 gw_corr_lev_occ, vec_Sigma_x_minus_vxc_gw11) 516 517 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vec_Sigma_x_gw 518 INTEGER, INTENT(IN) :: nmo 519 TYPE(cp_fm_type), POINTER :: fm_mat_S_gw 520 TYPE(cp_para_env_type), POINTER :: para_env 521 INTEGER, INTENT(IN) :: num_integ_group, color_rpa_group, homo, & 522 gw_corr_lev_occ 523 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec_Sigma_x_minus_vxc_gw11 524 525 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_vec_sigma_x', & 526 routineP = moduleN//':'//routineN 527 528 INTEGER :: handle, iiB, jjB, m_global, n_global, & 529 ncol_local, nm_global, nrow_local 530 INTEGER, DIMENSION(:), POINTER :: row_indices 531 532 CALL timeset(routineN, handle) 533 534 CALL cp_fm_get_info(matrix=fm_mat_S_gw, & 535 nrow_local=nrow_local, & 536 ncol_local=ncol_local, & 537 row_indices=row_indices) 538 539 CALL mp_sync(para_env%group) 540 541 ! loop over (nm) index 542 DO iiB = 1, nrow_local 543 544 ! this is needed for correct values within parallelization 545 IF (MODULO(1, num_integ_group) /= color_rpa_group) CYCLE 546 547 nm_global = row_indices(iiB) 548 549 ! transform the index nm to n and m, formulae copied from Mauro's code 550 n_global = MAX(1, nm_global - 1)/nmo + 1 551 m_global = nm_global - (n_global - 1)*nmo 552 n_global = n_global + homo - gw_corr_lev_occ 553 554 IF (m_global <= homo) THEN 555 556 ! loop over auxiliary basis functions 557 DO jjB = 1, ncol_local 558 559 ! Sigma_x_n = -sum_m^occ sum_P (B_(nm)^P)^2 560 vec_Sigma_x_gw(n_global, 1) = & 561 vec_Sigma_x_gw(n_global, 1) - & 562 fm_mat_S_gw%local_data(iiB, jjB)**2 563 564 END DO 565 566 END IF 567 568 END DO 569 570 CALL mp_sum(vec_Sigma_x_gw, para_env%group) 571 572 vec_Sigma_x_minus_vxc_gw11(:) = & 573 vec_Sigma_x_minus_vxc_gw11(:) + & 574 vec_Sigma_x_gw(:, 1) 575 576 CALL timestop(handle) 577 578 END SUBROUTINE get_vec_sigma_x 579 580! ************************************************************************************************** 581!> \brief ... 582!> \param fm_mat_S_gw_work ... 583!> \param vec_W_gw ... 584!> \param vec_Sigma_c_gw ... 585!> \param vec_omega_fit_gw ... 586!> \param vec_Sigma_x_minus_vxc_gw ... 587!> \param Eigenval_last ... 588!> \param Eigenval_scf ... 589!> \param do_periodic ... 590!> \param matrix_berry_re_mo_mo ... 591!> \param matrix_berry_im_mo_mo ... 592!> \param kpoints ... 593!> \param do_ri_Sigma_x ... 594!> \param vec_Sigma_x_gw ... 595!> \param my_do_gw ... 596!> \param fm_mat_S_gw_work_beta ... 597!> \param vec_W_gw_beta ... 598!> \param vec_Sigma_c_gw_beta ... 599!> \param Eigenval_last_beta ... 600!> \param Eigenval_scf_beta ... 601!> \param vec_Sigma_x_gw_beta ... 602! ************************************************************************************************** 603 SUBROUTINE deallocate_matrices_gw(fm_mat_S_gw_work, vec_W_gw, vec_Sigma_c_gw, vec_omega_fit_gw, & 604 vec_Sigma_x_minus_vxc_gw, Eigenval_last, & 605 Eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, kpoints, & 606 do_ri_Sigma_x, vec_Sigma_x_gw, my_do_gw, & 607 fm_mat_S_gw_work_beta, vec_W_gw_beta, & 608 vec_Sigma_c_gw_beta, Eigenval_last_beta, Eigenval_scf_beta, vec_Sigma_x_gw_beta) 609 610 TYPE(cp_fm_type), POINTER :: fm_mat_S_gw_work 611 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 612 INTENT(INOUT) :: vec_W_gw 613 COMPLEX(KIND=dp), ALLOCATABLE, & 614 DIMENSION(:, :, :), INTENT(INOUT) :: vec_Sigma_c_gw 615 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 616 INTENT(INOUT) :: vec_omega_fit_gw 617 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & 618 INTENT(INOUT) :: vec_Sigma_x_minus_vxc_gw 619 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 620 INTENT(INOUT) :: Eigenval_last, Eigenval_scf 621 LOGICAL, INTENT(IN) :: do_periodic 622 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_re_mo_mo, & 623 matrix_berry_im_mo_mo 624 TYPE(kpoint_type), POINTER :: kpoints 625 LOGICAL, INTENT(IN) :: do_ri_Sigma_x 626 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 627 INTENT(INOUT) :: vec_Sigma_x_gw 628 LOGICAL, INTENT(IN) :: my_do_gw 629 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_mat_S_gw_work_beta 630 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 631 INTENT(INOUT), OPTIONAL :: vec_W_gw_beta 632 COMPLEX(KIND=dp), ALLOCATABLE, & 633 DIMENSION(:, :, :), INTENT(INOUT), OPTIONAL :: vec_Sigma_c_gw_beta 634 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 635 INTENT(INOUT), OPTIONAL :: Eigenval_last_beta, Eigenval_scf_beta 636 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 637 INTENT(INOUT), OPTIONAL :: vec_Sigma_x_gw_beta 638 639 CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_matrices_gw', & 640 routineP = moduleN//':'//routineN 641 642 INTEGER :: handle 643 LOGICAL :: my_open_shell 644 645 CALL timeset(routineN, handle) 646 647 my_open_shell = .FALSE. 648 IF (PRESENT(fm_mat_S_gw_work_beta) .AND. PRESENT(vec_W_gw_beta) .AND. PRESENT(vec_Sigma_c_gw_beta) .AND. & 649 PRESENT(Eigenval_last_beta) .AND. PRESENT(Eigenval_scf_beta) .AND. PRESENT(vec_Sigma_x_gw_beta)) THEN 650 my_open_shell = .TRUE. 651 END IF 652 653 IF (my_do_gw) THEN 654 CALL cp_fm_release(fm_mat_S_gw_work) 655 DEALLOCATE (vec_Sigma_x_minus_vxc_gw) 656 DEALLOCATE (vec_W_gw) 657 IF (my_open_shell) THEN 658 CALL cp_fm_release(fm_mat_S_gw_work_beta) 659 DEALLOCATE (vec_W_gw_beta) 660 END IF 661 END IF 662 663 DEALLOCATE (vec_Sigma_c_gw) 664 DEALLOCATE (vec_omega_fit_gw) 665 DEALLOCATE (Eigenval_last) 666 DEALLOCATE (Eigenval_scf) 667 IF (my_open_shell) THEN 668 DEALLOCATE (vec_Sigma_c_gw_beta) 669 DEALLOCATE (Eigenval_last_beta) 670 DEALLOCATE (Eigenval_scf_beta) 671 END IF 672 673 IF (do_periodic) THEN 674 CALL dbcsr_deallocate_matrix_set(matrix_berry_re_mo_mo) 675 CALL dbcsr_deallocate_matrix_set(matrix_berry_im_mo_mo) 676 CALL kpoint_release(kpoints) 677 END IF 678 IF (do_ri_Sigma_x) THEN 679 DEALLOCATE (vec_Sigma_x_gw) 680 IF (my_open_shell) THEN 681 DEALLOCATE (vec_Sigma_x_gw_beta) 682 END IF 683 END IF 684 685 CALL timestop(handle) 686 687 END SUBROUTINE deallocate_matrices_gw 688 689! ************************************************************************************************** 690!> \brief ... 691!> \param weights_cos_tf_w_to_t ... 692!> \param weights_sin_tf_t_to_w ... 693!> \param do_ic_model ... 694!> \param do_kpoints_cubic_RPA ... 695!> \param fm_mat_W ... 696!> \param t_3c_overl_int_gw_RI ... 697!> \param t_3c_overl_int_gw_AO ... 698!> \param t_3c_overl_nnP_ic ... 699!> \param t_3c_overl_nnP_ic_reflected ... 700!> \param mat_W ... 701!> \param ikp_local ... 702!> \param cfm_mat_W_kp_tau ... 703!> \param t_3c_overl_int_gw_RI_beta ... 704!> \param t_3c_overl_int_gw_AO_beta ... 705!> \param t_3c_overl_nnP_ic_beta ... 706!> \param t_3c_overl_nnP_ic_reflected_beta ... 707! ************************************************************************************************** 708 SUBROUTINE deallocate_matrices_gw_im_time(weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, do_ic_model, do_kpoints_cubic_RPA, & 709 fm_mat_W, t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, & 710 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, mat_W, & 711 ikp_local, cfm_mat_W_kp_tau, & 712 t_3c_overl_int_gw_RI_beta, t_3c_overl_int_gw_AO_beta, & 713 t_3c_overl_nnP_ic_beta, t_3c_overl_nnP_ic_reflected_beta) 714 715 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 716 INTENT(INOUT) :: weights_cos_tf_w_to_t, & 717 weights_sin_tf_t_to_w 718 LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA 719 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fm_mat_W 720 TYPE(dbcsr_t_type) :: t_3c_overl_int_gw_RI, & 721 t_3c_overl_int_gw_AO, & 722 t_3c_overl_nnP_ic, & 723 t_3c_overl_nnP_ic_reflected 724 TYPE(dbcsr_type), POINTER :: mat_W 725 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: ikp_local 726 TYPE(cp_cfm_p_type), DIMENSION(:, :), POINTER :: cfm_mat_W_kp_tau 727 TYPE(dbcsr_t_type), OPTIONAL :: t_3c_overl_int_gw_RI_beta, t_3c_overl_int_gw_AO_beta, & 728 t_3c_overl_nnP_ic_beta, t_3c_overl_nnP_ic_reflected_beta 729 730 CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_matrices_gw_im_time', & 731 routineP = moduleN//':'//routineN 732 733 INTEGER :: handle, ikp, jquad 734 LOGICAL :: my_open_shell 735 736 CALL timeset(routineN, handle) 737 738 my_open_shell = .FALSE. 739 IF (PRESENT(t_3c_overl_int_gw_RI_beta) .AND. PRESENT(t_3c_overl_int_gw_AO_beta) .AND. PRESENT(t_3c_overl_nnP_ic_beta) & 740 .AND. PRESENT(t_3c_overl_nnP_ic_reflected_beta)) THEN 741 my_open_shell = .TRUE. 742 END IF 743 744 IF (ALLOCATED(weights_cos_tf_w_to_t)) DEALLOCATE (weights_cos_tf_w_to_t) 745 IF (ALLOCATED(weights_sin_tf_t_to_w)) DEALLOCATE (weights_sin_tf_t_to_w) 746 747 IF (.NOT. do_kpoints_cubic_RPA) THEN 748 749 DO jquad = 1, SIZE(fm_mat_W, 1) 750 CALL cp_fm_release(fm_mat_W(jquad)%matrix) 751 END DO 752 753 DEALLOCATE (fm_mat_W) 754 755 CALL dbcsr_release_P(mat_W) 756 757 ELSE 758 DO jquad = 1, SIZE(cfm_mat_W_kp_tau, 2) 759 DO ikp = 1, SIZE(cfm_mat_W_kp_tau, 1) 760 IF (.NOT. (ANY(ikp_local(:) == ikp))) CYCLE 761 CALL cp_cfm_release(cfm_mat_W_kp_tau(ikp, jquad)%matrix) 762 END DO 763 END DO 764 DEALLOCATE (cfm_mat_W_kp_tau) 765 766 END IF 767 768 CALL dbcsr_t_destroy(t_3c_overl_int_gw_RI) 769 CALL dbcsr_t_destroy(t_3c_overl_int_gw_AO) 770 IF (PRESENT(t_3c_overl_int_gw_RI_beta)) CALL dbcsr_t_destroy(t_3c_overl_int_gw_RI_beta) 771 IF (PRESENT(t_3c_overl_int_gw_AO_beta)) CALL dbcsr_t_destroy(t_3c_overl_int_gw_AO_beta) 772 IF (do_ic_model) THEN 773 CALL dbcsr_t_destroy(t_3c_overl_nnP_ic) 774 CALL dbcsr_t_destroy(t_3c_overl_nnP_ic_reflected) 775 IF (PRESENT(t_3c_overl_nnP_ic_beta)) CALL dbcsr_t_destroy(t_3c_overl_nnP_ic_beta) 776 IF (PRESENT(t_3c_overl_nnP_ic_reflected_beta)) CALL dbcsr_t_destroy(t_3c_overl_nnP_ic_reflected_beta) 777 ENDIF 778 779 CALL timestop(handle) 780 781 END SUBROUTINE deallocate_matrices_gw_im_time 782 783! ************************************************************************************************** 784!> \brief ... 785!> \param vec_Sigma_c_gw ... 786!> \param dimen_nm_gw ... 787!> \param dimen_RI ... 788!> \param gw_corr_lev_occ ... 789!> \param gw_corr_lev_virt ... 790!> \param homo ... 791!> \param jquad ... 792!> \param nmo ... 793!> \param num_fit_points ... 794!> \param num_integ_points ... 795!> \param do_bse ... 796!> \param do_im_time ... 797!> \param do_periodic ... 798!> \param first_cycle_periodic_correction ... 799!> \param fermi_level_offset ... 800!> \param fermi_level_offset_input ... 801!> \param omega ... 802!> \param Eigenval ... 803!> \param delta_corr ... 804!> \param tau_tj ... 805!> \param tj ... 806!> \param vec_omega_fit_gw ... 807!> \param vec_W_gw ... 808!> \param wj ... 809!> \param weights_cos_tf_w_to_t ... 810!> \param fm_mat_W ... 811!> \param fm_mat_L ... 812!> \param fm_mat_Q ... 813!> \param fm_mat_Q_static_bse ... 814!> \param fm_mat_R_gw ... 815!> \param fm_mat_S_gw ... 816!> \param fm_mat_S_gw_work ... 817!> \param fm_mat_work ... 818!> \param mo_coeff ... 819!> \param para_env ... 820!> \param para_env_RPA ... 821!> \param matrix_berry_im_mo_mo ... 822!> \param matrix_berry_re_mo_mo ... 823!> \param kpoints ... 824!> \param qs_env ... 825!> \param mp2_env ... 826!> \param do_kpoints_cubic_RPA ... 827!> \param do_kpoints_from_Gamma ... 828!> \param vec_Sigma_c_gw_beta ... 829!> \param gw_corr_lev_occ_beta ... 830!> \param homo_beta ... 831!> \param Eigenval_beta ... 832!> \param vec_W_gw_beta ... 833!> \param fm_mat_S_gw_beta ... 834!> \param fm_mat_S_gw_work_beta ... 835! ************************************************************************************************** 836 SUBROUTINE GW_matrix_operations(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw_corr_lev_occ, & 837 gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, num_integ_points, & 838 do_bse, do_im_time, do_periodic, first_cycle_periodic_correction, & 839 fermi_level_offset, fermi_level_offset_input, & 840 omega, Eigenval, delta_corr, tau_tj, tj, vec_omega_fit_gw, & 841 vec_W_gw, wj, weights_cos_tf_w_to_t, fm_mat_W, fm_mat_L, & 842 fm_mat_Q, fm_mat_Q_static_bse, fm_mat_R_gw, fm_mat_S_gw, & 843 fm_mat_S_gw_work, fm_mat_work, mo_coeff, para_env, & 844 para_env_RPA, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, & 845 kpoints, qs_env, mp2_env, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, & 846 vec_Sigma_c_gw_beta, gw_corr_lev_occ_beta, homo_beta, Eigenval_beta, & 847 vec_W_gw_beta, fm_mat_S_gw_beta, fm_mat_S_gw_work_beta) 848 849 COMPLEX(KIND=dp), ALLOCATABLE, & 850 DIMENSION(:, :, :), INTENT(INOUT) :: vec_Sigma_c_gw 851 INTEGER, INTENT(IN) :: dimen_nm_gw, dimen_RI, gw_corr_lev_occ, & 852 gw_corr_lev_virt, homo, jquad, nmo, & 853 num_fit_points, num_integ_points 854 LOGICAL, INTENT(IN) :: do_bse, do_im_time, do_periodic 855 LOGICAL, INTENT(INOUT) :: first_cycle_periodic_correction 856 REAL(KIND=dp), INTENT(INOUT) :: fermi_level_offset 857 REAL(KIND=dp), INTENT(IN) :: fermi_level_offset_input 858 REAL(KIND=dp), INTENT(INOUT) :: omega 859 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: Eigenval 860 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 861 INTENT(INOUT) :: delta_corr 862 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 863 INTENT(IN) :: tau_tj, tj, vec_omega_fit_gw 864 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 865 INTENT(INOUT) :: vec_W_gw 866 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 867 INTENT(IN) :: wj 868 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 869 INTENT(IN) :: weights_cos_tf_w_to_t 870 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fm_mat_W 871 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L 872 TYPE(cp_fm_type), POINTER :: fm_mat_Q, fm_mat_Q_static_bse, & 873 fm_mat_R_gw, fm_mat_S_gw, & 874 fm_mat_S_gw_work, fm_mat_work, mo_coeff 875 TYPE(cp_para_env_type), POINTER :: para_env, para_env_RPA 876 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_im_mo_mo, & 877 matrix_berry_re_mo_mo 878 TYPE(kpoint_type), POINTER :: kpoints 879 TYPE(qs_environment_type), POINTER :: qs_env 880 TYPE(mp2_type), POINTER :: mp2_env 881 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA, & 882 do_kpoints_from_Gamma 883 COMPLEX(KIND=dp), ALLOCATABLE, & 884 DIMENSION(:, :, :), INTENT(INOUT), OPTIONAL :: vec_Sigma_c_gw_beta 885 INTEGER, INTENT(IN), OPTIONAL :: gw_corr_lev_occ_beta, homo_beta 886 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), & 887 OPTIONAL :: Eigenval_beta 888 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 889 INTENT(INOUT), OPTIONAL :: vec_W_gw_beta 890 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_mat_S_gw_beta, fm_mat_S_gw_work_beta 891 892 CHARACTER(LEN=*), PARAMETER :: routineN = 'GW_matrix_operations', & 893 routineP = moduleN//':'//routineN 894 895 INTEGER :: handle, i_global, iiB, iquad, j_global, & 896 jjB, ncol_local, nrow_local 897 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 898 LOGICAL :: my_open_shell 899 REAL(KIND=dp) :: tau, weight 900 901 CALL timeset(routineN, handle) 902 903 my_open_shell = .FALSE. 904 IF (PRESENT(vec_Sigma_c_gw_beta) .AND. PRESENT(gw_corr_lev_occ_beta) .AND. PRESENT(homo_beta) & 905 .AND. PRESENT(Eigenval_beta) .AND. PRESENT(vec_W_gw_beta) .AND. PRESENT(fm_mat_S_gw_beta) & 906 .AND. PRESENT(fm_mat_S_gw_work_beta)) THEN 907 my_open_shell = .TRUE. 908 END IF 909 910 ! Fermi level offset should have a maximum such that the Fermi level of occupied orbitals 911 ! is always closer to occupied orbitals than to virtual orbitals and vice versa 912 ! that means, the Fermi level offset is at most as big as half the bandgap 913 IF (fermi_level_offset_input > 0.5_dp*(Eigenval(homo + 1) - Eigenval(homo))) THEN 914 fermi_level_offset = (Eigenval(homo + 1) - Eigenval(homo))*0.5_dp 915 ELSE 916 fermi_level_offset = fermi_level_offset_input 917 END IF 918 IF (my_open_shell) THEN 919 IF (fermi_level_offset > 0.5_dp*(Eigenval_beta(homo_beta + 1) - Eigenval_beta(homo_beta))) THEN 920 fermi_level_offset = (Eigenval_beta(homo_beta + 1) - Eigenval_beta(homo_beta))*0.5_dp 921 END IF 922 END IF 923 924 CALL cp_fm_get_info(matrix=fm_mat_Q, & 925 nrow_local=nrow_local, & 926 ncol_local=ncol_local, & 927 row_indices=row_indices, & 928 col_indices=col_indices) 929 930 IF (.NOT. do_im_time) THEN 931 ! calculate [1+Q(iw')]^-1 932 CALL cp_fm_cholesky_invert(fm_mat_Q) 933 ! symmetrize the result, fm_mat_R_gw is only temporary work matrix 934 CALL cp_fm_upper_to_full(fm_mat_Q, fm_mat_R_gw) 935 936 IF (do_bse .AND. jquad == 1) THEN 937 CALL cp_fm_to_fm(fm_mat_Q, fm_mat_Q_static_bse) 938 END IF 939 940 ! periodic correction for GW 941 IF (do_periodic) THEN 942 CALL calc_periodic_correction(delta_corr, qs_env, para_env, para_env_RPA, & 943 mp2_env%ri_g0w0%kp_grid, homo, nmo, gw_corr_lev_occ, & 944 gw_corr_lev_virt, omega, mo_coeff, Eigenval, & 945 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, & 946 first_cycle_periodic_correction, kpoints, & 947 mp2_env%ri_g0w0%do_mo_coeff_gamma, & 948 mp2_env%ri_g0w0%num_kp_grids, mp2_env%ri_g0w0%eps_kpoint, & 949 mp2_env%ri_g0w0%do_extra_kpoints, & 950 mp2_env%ri_g0w0%do_aux_bas_gw, mp2_env%ri_g0w0%frac_aux_mos) 951 END IF 952 953 ! subtract 1 from the diagonal to get rid of exchange self-energy 954!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 955!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_mat_Q,dimen_RI) 956 DO jjB = 1, ncol_local 957 j_global = col_indices(jjB) 958 DO iiB = 1, nrow_local 959 i_global = row_indices(iiB) 960 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 961 fm_mat_Q%local_data(iiB, jjB) = fm_mat_Q%local_data(iiB, jjB) - 1.0_dp 962 END IF 963 END DO 964 END DO 965 966 CALL compute_self_energy_gw(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw_corr_lev_occ, homo, jquad, nmo, num_fit_points, & 967 do_periodic, fermi_level_offset, omega, Eigenval, delta_corr, vec_omega_fit_gw, vec_W_gw, & 968 wj, fm_mat_Q, fm_mat_S_gw, fm_mat_S_gw_work) 969 970 IF (my_open_shell) THEN 971 CALL compute_self_energy_gw(vec_Sigma_c_gw_beta, dimen_nm_gw, dimen_RI, gw_corr_lev_occ_beta, homo_beta, jquad, nmo, & 972 num_fit_points, do_periodic, fermi_level_offset, omega, Eigenval_beta, delta_corr, & 973 vec_omega_fit_gw, vec_W_gw_beta, wj, fm_mat_Q, fm_mat_S_gw_beta, fm_mat_S_gw_work_beta) 974 END IF 975 976 END IF ! GW 977 978 ! cubic scaling GW calculation for molecules 979 IF (do_im_time .AND. .NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN 980 981 ! calculate [1+Q(iw')]^-1 982 CALL cp_fm_cholesky_invert(fm_mat_Q) 983 984 ! symmetrize the result 985 CALL cp_fm_upper_to_full(fm_mat_Q, fm_mat_work) 986 987 ! subtract 1 from the diagonal to get rid of exchange self-energy 988!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 989!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_mat_Q,dimen_RI) 990 DO jjB = 1, ncol_local 991 j_global = col_indices(jjB) 992 DO iiB = 1, nrow_local 993 i_global = row_indices(iiB) 994 IF (j_global == i_global .AND. i_global <= dimen_RI) THEN 995 fm_mat_Q%local_data(iiB, jjB) = fm_mat_Q%local_data(iiB, jjB) - 1.0_dp 996 END IF 997 END DO 998 END DO 999 1000 ! multiply with L from the left and the right to get the screened Coulomb interaction 1001 CALL cp_gemm('T', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_mat_L(1, 1)%matrix, fm_mat_Q, & 1002 0.0_dp, fm_mat_work) 1003 1004 CALL cp_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_mat_work, fm_mat_L(1, 1)%matrix, & 1005 0.0_dp, fm_mat_Q) 1006 1007 ! Fourier transform from w to t 1008 DO iquad = 1, num_integ_points 1009 1010 omega = tj(jquad) 1011 tau = tau_tj(iquad) 1012 weight = weights_cos_tf_w_to_t(iquad, jquad)*COS(tau*omega) 1013 1014 IF (jquad == 1) THEN 1015 1016 CALL cp_fm_set_all(matrix=fm_mat_W(iquad)%matrix, alpha=0.0_dp) 1017 1018 END IF 1019 1020 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_W(iquad)%matrix, beta=weight, matrix_b=fm_mat_Q) 1021 1022 END DO 1023 1024 END IF 1025 1026 CALL timestop(handle) 1027 1028 END SUBROUTINE GW_matrix_operations 1029 1030! ************************************************************************************************** 1031!> \brief ... 1032!> \param vec_Sigma_c_gw ... 1033!> \param dimen_nm_gw ... 1034!> \param dimen_RI ... 1035!> \param gw_corr_lev_occ ... 1036!> \param homo ... 1037!> \param jquad ... 1038!> \param nmo ... 1039!> \param num_fit_points ... 1040!> \param do_periodic ... 1041!> \param fermi_level_offset ... 1042!> \param omega ... 1043!> \param Eigenval ... 1044!> \param delta_corr ... 1045!> \param vec_omega_fit_gw ... 1046!> \param vec_W_gw ... 1047!> \param wj ... 1048!> \param fm_mat_Q ... 1049!> \param fm_mat_S_gw ... 1050!> \param fm_mat_S_gw_work ... 1051! ************************************************************************************************** 1052 SUBROUTINE compute_self_energy_gw(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw_corr_lev_occ, homo, jquad, nmo, num_fit_points, & 1053 do_periodic, fermi_level_offset, omega, Eigenval, delta_corr, vec_omega_fit_gw, vec_W_gw, & 1054 wj, fm_mat_Q, fm_mat_S_gw, fm_mat_S_gw_work) 1055 1056 COMPLEX(KIND=dp), ALLOCATABLE, & 1057 DIMENSION(:, :, :), INTENT(INOUT) :: vec_Sigma_c_gw 1058 INTEGER, INTENT(IN) :: dimen_nm_gw, dimen_RI, gw_corr_lev_occ, & 1059 homo, jquad, nmo, num_fit_points 1060 LOGICAL, INTENT(IN) :: do_periodic 1061 REAL(KIND=dp), INTENT(IN) :: fermi_level_offset 1062 REAL(KIND=dp), INTENT(INOUT) :: omega 1063 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: Eigenval 1064 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: delta_corr, vec_omega_fit_gw 1065 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: vec_W_gw 1066 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wj 1067 TYPE(cp_fm_type), POINTER :: fm_mat_Q, fm_mat_S_gw, fm_mat_S_gw_work 1068 1069 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_self_energy_gw', & 1070 routineP = moduleN//':'//routineN 1071 COMPLEX(KIND=dp), PARAMETER :: im_unit = (0.0_dp, 1.0_dp) 1072 1073 INTEGER :: handle, iiB, iquad, jjB, m_global, & 1074 n_global, ncol_local, nm_global, & 1075 nrow_local 1076 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 1077 REAL(KIND=dp) :: delta_corr_nn, e_fermi, omega_i, & 1078 sign_occ_virt 1079 1080 CALL timeset(routineN, handle) 1081 1082 ! S_work_(nm)Q = B_(nm)P * ([1+Q]^-1-1)_PQ 1083 CALL cp_gemm(transa="N", transb="N", m=dimen_nm_gw, n=dimen_RI, k=dimen_RI, alpha=1.0_dp, & 1084 matrix_a=fm_mat_S_gw, matrix_b=fm_mat_Q, beta=0.0_dp, & 1085 matrix_c=fm_mat_S_gw_work) 1086 1087 CALL cp_fm_get_info(matrix=fm_mat_S_gw, & 1088 nrow_local=nrow_local, & 1089 ncol_local=ncol_local, & 1090 row_indices=row_indices, & 1091 col_indices=col_indices) 1092 1093 ! vector W_(nm) = S_work_(nm)Q * [B_(nm)Q]^T 1094 1095 vec_W_gw = 0.0_dp 1096 1097 DO iiB = 1, nrow_local 1098 nm_global = row_indices(iiB) 1099 DO jjB = 1, ncol_local 1100 vec_W_gw(nm_global) = vec_W_gw(nm_global) + & 1101 fm_mat_S_gw_work%local_data(iiB, jjB)*fm_mat_S_gw%local_data(iiB, jjB) 1102 END DO 1103 1104 ! transform the index nm of vec_W_gw back to n and m, formulae copied from Mauro's code 1105 n_global = MAX(1, nm_global - 1)/nmo + 1 1106 m_global = nm_global - (n_global - 1)*nmo 1107 n_global = n_global + homo - gw_corr_lev_occ 1108 1109 ! compute self-energy for imaginary frequencies 1110 DO iquad = 1, num_fit_points 1111 1112 ! for occ orbitals, we compute the self-energy for negative frequencies 1113 IF (n_global <= homo) THEN 1114 sign_occ_virt = -1.0_dp 1115 ELSE 1116 sign_occ_virt = 1.0_dp 1117 END IF 1118 1119 omega_i = vec_omega_fit_gw(iquad)*sign_occ_virt 1120 1121 ! set the Fermi energy for occ orbitals slightly above the HOMO and 1122 ! for virt orbitals slightly below the LUMO 1123 IF (n_global <= homo) THEN 1124 e_fermi = Eigenval(homo) + fermi_level_offset 1125 ELSE 1126 e_fermi = Eigenval(homo + 1) - fermi_level_offset 1127 END IF 1128 1129 ! add here the periodic correction 1130 IF (do_periodic .AND. col_indices(1) == 1 .AND. n_global == m_global) THEN 1131 delta_corr_nn = delta_corr(n_global) 1132 ELSE 1133 delta_corr_nn = 0.0_dp 1134 END IF 1135 1136 ! update the self-energy (use that vec_W_gw(iw) is symmetric), divide the integration 1137 ! weight by 2, because the integration is from -infty to +infty and not just 0 to +infty 1138 ! as for RPA, also we need for virtual orbitals a complex conjugate 1139 vec_Sigma_c_gw(n_global - homo + gw_corr_lev_occ, iquad, 1) = & 1140 vec_Sigma_c_gw(n_global - homo + gw_corr_lev_occ, iquad, 1) - & 1141 0.5_dp/pi*wj(jquad)/2.0_dp*(vec_W_gw(nm_global) + delta_corr_nn)* & 1142 (1.0_dp/(im_unit*(omega + omega_i) + e_fermi - Eigenval(m_global)) + & 1143 1.0_dp/(im_unit*(-omega + omega_i) + e_fermi - Eigenval(m_global))) 1144 END DO 1145 1146 END DO 1147 1148 CALL timestop(handle) 1149 1150 END SUBROUTINE compute_self_energy_gw 1151 1152! ************************************************************************************************** 1153!> \brief ... 1154!> \param vec_Sigma_c_gw ... 1155!> \param count_ev_sc_GW ... 1156!> \param gw_corr_lev_occ ... 1157!> \param gw_corr_lev_tot ... 1158!> \param gw_corr_lev_virt ... 1159!> \param homo ... 1160!> \param nmo ... 1161!> \param num_fit_points ... 1162!> \param num_integ_points ... 1163!> \param num_points_corr ... 1164!> \param unit_nr ... 1165!> \param do_apply_ic_corr_to_gw ... 1166!> \param do_im_time ... 1167!> \param do_periodic ... 1168!> \param do_ri_Sigma_x ... 1169!> \param first_cycle_periodic_correction ... 1170!> \param e_fermi ... 1171!> \param eps_filter ... 1172!> \param fermi_level_offset ... 1173!> \param delta_corr ... 1174!> \param Eigenval ... 1175!> \param Eigenval_last ... 1176!> \param Eigenval_scf ... 1177!> \param iter_sc_GW0 ... 1178!> \param exit_ev_gw ... 1179!> \param tau_tj ... 1180!> \param tj ... 1181!> \param vec_omega_fit_gw ... 1182!> \param vec_Sigma_x_gw ... 1183!> \param ic_corr_list ... 1184!> \param weights_cos_tf_t_to_w ... 1185!> \param weights_sin_tf_t_to_w ... 1186!> \param fm_mo_coeff_occ_scaled ... 1187!> \param fm_mo_coeff_virt_scaled ... 1188!> \param fm_mo_coeff_occ ... 1189!> \param fm_mo_coeff_virt ... 1190!> \param fm_scaled_dm_occ_tau ... 1191!> \param fm_scaled_dm_virt_tau ... 1192!> \param mo_coeff ... 1193!> \param fm_mat_W ... 1194!> \param para_env ... 1195!> \param para_env_RPA ... 1196!> \param mat_dm ... 1197!> \param mat_SinvVSinv ... 1198!> \param t_3c_overl_int_gw_RI ... 1199!> \param t_3c_overl_int_gw_AO ... 1200!> \param matrix_berry_im_mo_mo ... 1201!> \param matrix_berry_re_mo_mo ... 1202!> \param mat_W ... 1203!> \param matrix_s ... 1204!> \param kpoints ... 1205!> \param mp2_env ... 1206!> \param qs_env ... 1207!> \param nkp_self_energy ... 1208!> \param do_kpoints_cubic_RPA ... 1209!> \param Eigenval_kp ... 1210!> \param Eigenval_scf_kp ... 1211!> \param vec_Sigma_c_gw_beta ... 1212!> \param gw_corr_lev_occ_beta ... 1213!> \param gw_corr_lev_virt_beta ... 1214!> \param homo_beta ... 1215!> \param e_fermi_beta ... 1216!> \param Eigenval_beta ... 1217!> \param Eigenval_last_beta ... 1218!> \param Eigenval_scf_beta ... 1219!> \param vec_Sigma_x_gw_beta ... 1220!> \param ic_corr_list_beta ... 1221!> \param fm_mo_coeff_occ_beta ... 1222!> \param fm_mo_coeff_virt_beta ... 1223!> \param t_3c_overl_int_gw_RI_beta ... 1224!> \param t_3c_overl_int_gw_AO_beta ... 1225! ************************************************************************************************** 1226 SUBROUTINE compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, & 1227 gw_corr_lev_tot, gw_corr_lev_virt, homo, & 1228 nmo, num_fit_points, num_integ_points, & 1229 num_points_corr, unit_nr, do_apply_ic_corr_to_gw, do_im_time, & 1230 do_periodic, do_ri_Sigma_x, & 1231 first_cycle_periodic_correction, e_fermi, eps_filter, & 1232 fermi_level_offset, delta_corr, Eigenval, & 1233 Eigenval_last, Eigenval_scf, iter_sc_GW0, exit_ev_gw, tau_tj, tj, & 1234 vec_omega_fit_gw, vec_Sigma_x_gw, ic_corr_list, & 1235 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, & 1236 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_mo_coeff_occ, & 1237 fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, & 1238 mo_coeff, fm_mat_W, para_env, para_env_RPA, mat_dm, mat_SinvVSinv, & 1239 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, matrix_berry_im_mo_mo, & 1240 matrix_berry_re_mo_mo, mat_W, matrix_s, & 1241 kpoints, mp2_env, qs_env, & 1242 nkp_self_energy, do_kpoints_cubic_RPA, Eigenval_kp, Eigenval_scf_kp, & 1243 vec_Sigma_c_gw_beta, gw_corr_lev_occ_beta, gw_corr_lev_virt_beta, homo_beta, & 1244 e_fermi_beta, Eigenval_beta, Eigenval_last_beta, Eigenval_scf_beta, & 1245 vec_Sigma_x_gw_beta, ic_corr_list_beta, fm_mo_coeff_occ_beta, fm_mo_coeff_virt_beta, & 1246 t_3c_overl_int_gw_RI_beta, t_3c_overl_int_gw_AO_beta) 1247 1248 COMPLEX(KIND=dp), ALLOCATABLE, & 1249 DIMENSION(:, :, :), INTENT(INOUT) :: vec_Sigma_c_gw 1250 INTEGER, INTENT(IN) :: count_ev_sc_GW, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, & 1251 homo, nmo, num_fit_points, num_integ_points 1252 INTEGER :: num_points_corr 1253 INTEGER, INTENT(IN) :: unit_nr 1254 LOGICAL, INTENT(IN) :: do_apply_ic_corr_to_gw, do_im_time, & 1255 do_periodic, do_ri_Sigma_x 1256 LOGICAL, INTENT(INOUT) :: first_cycle_periodic_correction 1257 REAL(KIND=dp), INTENT(INOUT) :: e_fermi 1258 REAL(KIND=dp), INTENT(IN) :: eps_filter, fermi_level_offset 1259 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1260 INTENT(INOUT) :: delta_corr 1261 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: Eigenval 1262 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1263 INTENT(INOUT) :: Eigenval_last, Eigenval_scf 1264 INTEGER, INTENT(IN) :: iter_sc_GW0 1265 LOGICAL, INTENT(INOUT) :: exit_ev_gw 1266 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1267 INTENT(INOUT) :: tau_tj, tj, vec_omega_fit_gw 1268 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1269 INTENT(INOUT) :: vec_Sigma_x_gw 1270 REAL(KIND=dp), DIMENSION(:), POINTER :: ic_corr_list 1271 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1272 INTENT(IN) :: weights_cos_tf_t_to_w, & 1273 weights_sin_tf_t_to_w 1274 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, & 1275 fm_mo_coeff_occ, fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, mo_coeff 1276 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fm_mat_W 1277 TYPE(cp_para_env_type), POINTER :: para_env, para_env_RPA 1278 TYPE(dbcsr_p_type), INTENT(IN) :: mat_dm, mat_SinvVSinv 1279 TYPE(dbcsr_t_type) :: t_3c_overl_int_gw_RI, & 1280 t_3c_overl_int_gw_AO 1281 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_im_mo_mo, & 1282 matrix_berry_re_mo_mo 1283 TYPE(dbcsr_type), POINTER :: mat_W 1284 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 1285 TYPE(kpoint_type), POINTER :: kpoints 1286 TYPE(mp2_type), POINTER :: mp2_env 1287 TYPE(qs_environment_type), POINTER :: qs_env 1288 INTEGER, INTENT(IN) :: nkp_self_energy 1289 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA 1290 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1291 INTENT(INOUT) :: Eigenval_kp, Eigenval_scf_kp 1292 COMPLEX(KIND=dp), ALLOCATABLE, & 1293 DIMENSION(:, :, :), INTENT(INOUT), OPTIONAL :: vec_Sigma_c_gw_beta 1294 INTEGER, INTENT(IN), OPTIONAL :: gw_corr_lev_occ_beta, & 1295 gw_corr_lev_virt_beta, homo_beta 1296 REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: e_fermi_beta 1297 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), & 1298 OPTIONAL :: Eigenval_beta 1299 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1300 INTENT(INOUT), OPTIONAL :: Eigenval_last_beta, Eigenval_scf_beta 1301 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1302 INTENT(INOUT), OPTIONAL :: vec_Sigma_x_gw_beta 1303 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: ic_corr_list_beta 1304 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_mo_coeff_occ_beta, & 1305 fm_mo_coeff_virt_beta 1306 TYPE(dbcsr_t_type), OPTIONAL :: t_3c_overl_int_gw_RI_beta, & 1307 t_3c_overl_int_gw_AO_beta 1308 1309 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies', & 1310 routineP = moduleN//':'//routineN 1311 1312 INTEGER :: count_ev_sc_GW_print, count_sc_GW0, & 1313 count_sc_GW0_print, crossing_search, & 1314 handle, ikp, n_level_gw, num_poles 1315 LOGICAL :: my_open_shell 1316 REAL(KIND=dp) :: stop_crit 1317 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: m_value, m_value_beta, vec_gw_energ, & 1318 vec_gw_energ_beta, z_value, & 1319 z_value_beta 1320 1321 CALL timeset(routineN, handle) 1322 1323 my_open_shell = .FALSE. 1324 IF (PRESENT(vec_Sigma_c_gw_beta) .AND. PRESENT(gw_corr_lev_occ_beta) .AND. PRESENT(gw_corr_lev_virt_beta) .AND. & 1325 PRESENT(homo_beta) .AND. PRESENT(e_fermi_beta) .AND. PRESENT(Eigenval_beta) .AND. PRESENT(Eigenval_last_beta) & 1326 .AND. PRESENT(Eigenval_scf_beta) .AND. PRESENT(vec_Sigma_x_gw_beta) .AND. & 1327 PRESENT(ic_corr_list_beta) .AND. PRESENT(fm_mo_coeff_occ_beta) .AND. PRESENT(fm_mo_coeff_virt_beta) .AND. & 1328 PRESENT(t_3c_overl_int_gw_RI_beta) .AND. PRESENT(t_3c_overl_int_gw_AO_beta)) THEN 1329 my_open_shell = .TRUE. ! todo: this should be refactored (and any similar occurrence of my_open_shell) 1330 END IF 1331 1332 DO count_sc_GW0 = 1, iter_sc_GW0 1333 1334 ! postprocessing for cubic scaling GW calculation 1335 IF (do_im_time .AND. .NOT. do_kpoints_cubic_RPA) THEN 1336 num_points_corr = mp2_env%ri_g0w0%num_omega_points 1337 1338 CALL compute_self_energy_cubic_gw(num_integ_points, nmo, tau_tj, tj, & 1339 matrix_s, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, & 1340 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, & 1341 fm_scaled_dm_virt_tau, Eigenval, eps_filter, & 1342 e_fermi, fm_mat_W, & 1343 gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, & 1344 count_ev_sc_GW, count_sc_GW0, & 1345 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, & 1346 mat_W, mat_SinvVSinv, mat_dm, & 1347 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_Sigma_c_gw, & 1348 do_periodic, num_points_corr, delta_corr, qs_env, para_env, para_env_RPA, & 1349 mp2_env, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, & 1350 first_cycle_periodic_correction, kpoints, num_fit_points, mo_coeff, & 1351 do_ri_Sigma_x, vec_Sigma_x_gw, unit_nr) 1352 1353 IF (my_open_shell) THEN 1354 1355 CALL compute_self_energy_cubic_gw(num_integ_points, nmo, tau_tj, tj, & 1356 matrix_s, fm_mo_coeff_occ_beta, fm_mo_coeff_virt_beta, fm_mo_coeff_occ_scaled, & 1357 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, & 1358 fm_scaled_dm_virt_tau, Eigenval_beta, eps_filter, & 1359 e_fermi_beta, fm_mat_W, & 1360 gw_corr_lev_tot, gw_corr_lev_occ_beta, gw_corr_lev_virt_beta, homo_beta, & 1361 count_ev_sc_GW, count_sc_GW0, & 1362 t_3c_overl_int_gw_RI_beta, t_3c_overl_int_gw_AO_beta, & 1363 mat_W, mat_SinvVSinv, mat_dm, & 1364 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_Sigma_c_gw_beta, & 1365 do_periodic, num_points_corr, delta_corr, qs_env, para_env, para_env_RPA, & 1366 mp2_env, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, & 1367 first_cycle_periodic_correction, kpoints, num_fit_points, mo_coeff, & 1368 do_ri_Sigma_x, vec_Sigma_x_gw_beta, unit_nr, do_beta=.TRUE.) 1369 1370 END IF 1371 1372 END IF 1373 1374 IF (.NOT. do_im_time) THEN 1375 1376 CALL mp_sum(vec_Sigma_c_gw, para_env%group) 1377 1378 END IF 1379 1380 IF (do_periodic .AND. mp2_env%ri_g0w0%do_average_deg_levels) THEN 1381 1382 CALL average_degenerate_levels(vec_Sigma_c_gw, Eigenval(1 + homo - gw_corr_lev_occ:homo + gw_corr_lev_virt), & 1383 mp2_env%ri_g0w0%eps_eigenval) 1384 IF (my_open_shell) THEN 1385 CALL average_degenerate_levels(vec_Sigma_c_gw_beta, & 1386 Eigenval_beta(1 + homo_beta - gw_corr_lev_occ_beta: & 1387 homo_beta + gw_corr_lev_virt_beta), & 1388 mp2_env%ri_g0w0%eps_eigenval) 1389 END IF 1390 END IF 1391 1392 IF (my_open_shell .AND. .NOT. do_im_time) THEN 1393 CALL mp_sum(vec_Sigma_c_gw_beta, para_env%group) 1394 END IF 1395 1396 CALL mp_sync(para_env%group) 1397 1398 stop_crit = 1.0e-7 1399 num_poles = mp2_env%ri_g0w0%num_poles 1400 crossing_search = mp2_env%ri_g0w0%crossing_search 1401 1402 ! arrays storing the correlation self-energy, stat. error and z-shot value 1403 ALLOCATE (vec_gw_energ(gw_corr_lev_tot)) 1404 vec_gw_energ = 0.0_dp 1405 ALLOCATE (z_value(gw_corr_lev_tot)) 1406 z_value = 0.0_dp 1407 ALLOCATE (m_value(gw_corr_lev_tot)) 1408 m_value = 0.0_dp 1409 1410 ! the same for beta 1411 IF (my_open_shell) THEN 1412 ALLOCATE (vec_gw_energ_beta(gw_corr_lev_tot)) 1413 vec_gw_energ_beta = 0.0_dp 1414 ALLOCATE (z_value_beta(gw_corr_lev_tot)) 1415 z_value_beta = 0.0_dp 1416 ALLOCATE (m_value_beta(gw_corr_lev_tot)) 1417 m_value_beta = 0.0_dp 1418 END IF 1419 1420 ! for the normal code for molecules or Gamma only: nkp = 1 1421 DO ikp = 1, nkp_self_energy 1422 1423 IF (do_kpoints_cubic_RPA) THEN 1424 1425 vec_gw_energ = 0.0_dp 1426 z_value = 0.0_dp 1427 m_value = 0.0_dp 1428 1429 CALL get_eigenval_for_conti(Eigenval, Eigenval_scf, Eigenval_kp, Eigenval_scf_kp, kpoints, & 1430 ikp, my_open_shell) 1431 END IF 1432 1433 ! fit the self-energy on imaginary frequency axis and evaluate the fit on the MO energy of the SCF 1434 DO n_level_gw = 1, gw_corr_lev_tot 1435 ! processes perform different fits 1436 IF (MODULO(n_level_gw, para_env%num_pe) /= para_env%mepos) CYCLE 1437 1438 SELECT CASE (mp2_env%ri_g0w0%analytic_continuation) 1439 CASE (gw_two_pole_model) 1440 CALL fit_and_continuation_2pole(vec_gw_energ, vec_omega_fit_gw, & 1441 z_value, m_value, vec_Sigma_c_gw(:, :, ikp), & 1442 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), & 1443 Eigenval, Eigenval_scf, n_level_gw, gw_corr_lev_occ, num_poles, & 1444 num_fit_points, crossing_search, homo, stop_crit, & 1445 fermi_level_offset, do_im_time) 1446 1447 CASE (gw_pade_approx) 1448 CALL continuation_pade(vec_gw_energ, vec_omega_fit_gw, & 1449 z_value, m_value, vec_Sigma_c_gw(:, :, ikp), & 1450 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), & 1451 Eigenval, Eigenval_scf, n_level_gw, gw_corr_lev_occ, mp2_env%ri_g0w0%nparam_pade, & 1452 num_fit_points, crossing_search, homo, fermi_level_offset, & 1453 do_im_time, mp2_env%ri_g0w0%print_self_energy, count_ev_sc_GW) 1454 CASE DEFAULT 1455 CPABORT("Only two-model and Pade approximation are implemented.") 1456 END SELECT 1457 1458 IF (my_open_shell) THEN 1459 SELECT CASE (mp2_env%ri_g0w0%analytic_continuation) 1460 CASE (gw_two_pole_model) 1461 CALL fit_and_continuation_2pole( & 1462 vec_gw_energ_beta, vec_omega_fit_gw, & 1463 z_value_beta, m_value_beta, vec_Sigma_c_gw_beta(:, :, ikp), & 1464 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), & 1465 Eigenval_beta, Eigenval_scf_beta, n_level_gw, & 1466 gw_corr_lev_occ_beta, num_poles, & 1467 num_fit_points, crossing_search, homo_beta, stop_crit, & 1468 fermi_level_offset, do_im_time) 1469 CASE (gw_pade_approx) 1470 CALL continuation_pade(vec_gw_energ_beta, vec_omega_fit_gw, & 1471 z_value_beta, m_value_beta, vec_Sigma_c_gw_beta(:, :, ikp), & 1472 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), & 1473 Eigenval_beta, Eigenval_scf_beta, n_level_gw, & 1474 gw_corr_lev_occ_beta, mp2_env%ri_g0w0%nparam_pade, & 1475 num_fit_points, crossing_search, homo_beta, & 1476 fermi_level_offset, do_im_time, & 1477 mp2_env%ri_g0w0%print_self_energy, count_ev_sc_GW) 1478 CASE DEFAULT 1479 CPABORT("Only two-model and Pade approximation are implemented.") 1480 END SELECT 1481 1482 END IF 1483 1484 END DO ! n_level_gw 1485 1486 CALL mp_sum(vec_gw_energ, para_env%group) 1487 CALL mp_sum(z_value, para_env%group) 1488 CALL mp_sum(m_value, para_env%group) 1489 1490 IF (my_open_shell) THEN 1491 CALL mp_sum(vec_gw_energ_beta, para_env%group) 1492 CALL mp_sum(z_value_beta, para_env%group) 1493 CALL mp_sum(m_value_beta, para_env%group) 1494 END IF 1495 1496 IF (do_im_time .OR. mp2_env%ri_g0w0%iter_sc_GW0 == 1) THEN 1497 count_ev_sc_GW_print = count_ev_sc_GW 1498 count_sc_GW0_print = count_sc_GW0 1499 ELSE 1500 count_ev_sc_GW_print = count_sc_GW0 1501 count_sc_GW0_print = count_ev_sc_GW 1502 END IF 1503 1504 ! print the quasiparticle energies and update Eigenval in case you do eigenvalue self-consistent GW 1505 IF (my_open_shell) THEN 1506 1507 CALL print_and_update_for_ev_sc( & 1508 vec_gw_energ, & 1509 z_value, m_value, mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), Eigenval, & 1510 Eigenval_last, Eigenval_scf, gw_corr_lev_occ, gw_corr_lev_tot, & 1511 crossing_search, homo, unit_nr, count_ev_sc_GW_print, count_sc_GW0_print, & 1512 ikp, nkp_self_energy, kpoints, do_alpha=.TRUE.) 1513 1514 CALL print_and_update_for_ev_sc( & 1515 vec_gw_energ_beta, & 1516 z_value_beta, m_value_beta, mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), Eigenval_beta, & 1517 Eigenval_last_beta, Eigenval_scf_beta, gw_corr_lev_occ_beta, gw_corr_lev_tot, & 1518 crossing_search, homo_beta, unit_nr, count_ev_sc_GW_print, count_sc_GW0_print, & 1519 ikp, nkp_self_energy, kpoints, do_beta=.TRUE.) 1520 1521 IF (do_apply_ic_corr_to_gw .AND. count_ev_sc_GW == 1) THEN 1522 1523 CALL apply_ic_corr(Eigenval, Eigenval_scf, ic_corr_list, & 1524 gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, & 1525 homo, nmo, unit_nr, do_alpha=.TRUE.) 1526 1527 CALL apply_ic_corr(Eigenval_beta, Eigenval_scf_beta, ic_corr_list_beta, & 1528 gw_corr_lev_occ_beta, gw_corr_lev_virt_beta, gw_corr_lev_tot, & 1529 homo_beta, nmo, unit_nr, do_beta=.TRUE.) 1530 1531 END IF 1532 1533 ELSE 1534 1535 CALL print_and_update_for_ev_sc( & 1536 vec_gw_energ, & 1537 z_value, m_value, mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), Eigenval, & 1538 Eigenval_last, Eigenval_scf, gw_corr_lev_occ, gw_corr_lev_tot, & 1539 crossing_search, homo, unit_nr, count_ev_sc_GW_print, count_sc_GW0_print, & 1540 ikp, nkp_self_energy, kpoints) 1541 1542 IF (do_apply_ic_corr_to_gw .AND. count_ev_sc_GW == 1) THEN 1543 1544 CALL apply_ic_corr(Eigenval, Eigenval_scf, ic_corr_list, & 1545 gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, & 1546 homo, nmo, unit_nr) 1547 1548 END IF 1549 1550 END IF 1551 1552 END DO ! ikp 1553 1554 DEALLOCATE (z_value) 1555 DEALLOCATE (m_value) 1556 DEALLOCATE (vec_gw_energ) 1557 IF (my_open_shell) THEN 1558 DEALLOCATE (z_value_beta) 1559 DEALLOCATE (m_value_beta) 1560 DEALLOCATE (vec_gw_energ_beta) 1561 END IF 1562 1563 exit_ev_gw = .FALSE. 1564 1565 ! if HOMO-LUMO gap differs by less than mp2_env%ri_g0w0%eps_sc_iter, exit ev sc GW loop 1566 IF (ABS(Eigenval(homo) - Eigenval_last(homo) - Eigenval(homo + 1) + Eigenval_last(homo + 1)) & 1567 < mp2_env%ri_g0w0%eps_iter) THEN 1568 IF (count_sc_GW0 == 1) exit_ev_gw = .TRUE. 1569 EXIT 1570 END IF 1571 1572 CALL shift_unshifted_levels(Eigenval, Eigenval_last, gw_corr_lev_occ, gw_corr_lev_virt, & 1573 homo, nmo) 1574 IF (my_open_shell) THEN 1575 CALL shift_unshifted_levels(Eigenval_beta, Eigenval_last_beta, gw_corr_lev_occ_beta, & 1576 gw_corr_lev_virt_beta, homo_beta, nmo) 1577 END IF 1578 1579 ! in case of N^4 scaling GW, the scGW0 cycle is the eigenvalue sc cycle 1580 IF (.NOT. do_im_time) EXIT 1581 1582 END DO ! scGW0 1583 1584 CALL timestop(handle) 1585 1586 END SUBROUTINE compute_QP_energies 1587 1588! ************************************************************************************************** 1589!> \brief ... 1590!> \param delta_corr ... 1591!> \param qs_env ... 1592!> \param para_env ... 1593!> \param para_env_RPA ... 1594!> \param kp_grid ... 1595!> \param homo ... 1596!> \param nmo ... 1597!> \param gw_corr_lev_occ ... 1598!> \param gw_corr_lev_virt ... 1599!> \param omega ... 1600!> \param fm_mo_coeff ... 1601!> \param Eigenval ... 1602!> \param matrix_berry_re_mo_mo ... 1603!> \param matrix_berry_im_mo_mo ... 1604!> \param first_cycle_periodic_correction ... 1605!> \param kpoints ... 1606!> \param do_mo_coeff_Gamma_only ... 1607!> \param num_kp_grids ... 1608!> \param eps_kpoint ... 1609!> \param do_extra_kpoints ... 1610!> \param do_aux_bas ... 1611!> \param frac_aux_mos ... 1612! ************************************************************************************************** 1613 SUBROUTINE calc_periodic_correction(delta_corr, qs_env, para_env, para_env_RPA, kp_grid, homo, nmo, & 1614 gw_corr_lev_occ, gw_corr_lev_virt, omega, fm_mo_coeff, Eigenval, & 1615 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, & 1616 first_cycle_periodic_correction, kpoints, do_mo_coeff_Gamma_only, & 1617 num_kp_grids, eps_kpoint, do_extra_kpoints, do_aux_bas, frac_aux_mos) 1618 1619 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1620 INTENT(INOUT) :: delta_corr 1621 TYPE(qs_environment_type), POINTER :: qs_env 1622 TYPE(cp_para_env_type), POINTER :: para_env, para_env_RPA 1623 INTEGER, DIMENSION(:), POINTER :: kp_grid 1624 INTEGER, INTENT(IN) :: homo, nmo, gw_corr_lev_occ, & 1625 gw_corr_lev_virt 1626 REAL(KIND=dp), INTENT(IN) :: omega 1627 TYPE(cp_fm_type), POINTER :: fm_mo_coeff 1628 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 1629 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_re_mo_mo, & 1630 matrix_berry_im_mo_mo 1631 LOGICAL, INTENT(INOUT) :: first_cycle_periodic_correction 1632 TYPE(kpoint_type), POINTER :: kpoints 1633 LOGICAL, INTENT(IN) :: do_mo_coeff_Gamma_only 1634 INTEGER, INTENT(IN) :: num_kp_grids 1635 REAL(KIND=dp), INTENT(IN) :: eps_kpoint 1636 LOGICAL, INTENT(IN) :: do_extra_kpoints, do_aux_bas 1637 REAL(KIND=dp), INTENT(IN) :: frac_aux_mos 1638 1639 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_periodic_correction', & 1640 routineP = moduleN//':'//routineN 1641 1642 INTEGER :: handle 1643 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eps_head, eps_inv_head 1644 REAL(KIND=dp), DIMENSION(3, 3) :: h_inv 1645 1646 CALL timeset(routineN, handle) 1647 1648 IF (first_cycle_periodic_correction) THEN 1649 1650 CALL get_kpoints(qs_env, kpoints, kp_grid, num_kp_grids, para_env, h_inv, nmo, do_mo_coeff_Gamma_only, & 1651 do_extra_kpoints) 1652 1653 CALL get_berry_phase(qs_env, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, fm_mo_coeff, & 1654 para_env, do_mo_coeff_Gamma_only, homo, nmo, gw_corr_lev_virt, eps_kpoint, do_aux_bas, & 1655 frac_aux_mos) 1656 1657 END IF 1658 1659 CALL compute_eps_head_Berry(eps_head, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, para_env_RPA, & 1660 qs_env, homo, Eigenval, omega) 1661 1662 CALL compute_eps_inv_head(eps_inv_head, eps_head, kpoints) 1663 1664 CALL kpoint_sum_for_eps_inv_head_Berry(delta_corr, eps_inv_head, kpoints, qs_env, & 1665 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, & 1666 homo, gw_corr_lev_occ, gw_corr_lev_virt, para_env_RPA, & 1667 do_extra_kpoints) 1668 1669 DEALLOCATE (eps_head, eps_inv_head) 1670 1671 first_cycle_periodic_correction = .FALSE. 1672 1673 CALL timestop(handle) 1674 1675 END SUBROUTINE calc_periodic_correction 1676 1677! ************************************************************************************************** 1678!> \brief ... 1679!> \param eps_head ... 1680!> \param kpoints ... 1681!> \param matrix_berry_re_mo_mo ... 1682!> \param matrix_berry_im_mo_mo ... 1683!> \param para_env_RPA ... 1684!> \param qs_env ... 1685!> \param homo ... 1686!> \param Eigenval ... 1687!> \param omega ... 1688! ************************************************************************************************** 1689 SUBROUTINE compute_eps_head_Berry(eps_head, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, para_env_RPA, & 1690 qs_env, homo, Eigenval, omega) 1691 1692 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1693 INTENT(OUT) :: eps_head 1694 TYPE(kpoint_type), POINTER :: kpoints 1695 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_re_mo_mo, & 1696 matrix_berry_im_mo_mo 1697 TYPE(cp_para_env_type), POINTER :: para_env_RPA 1698 TYPE(qs_environment_type), POINTER :: qs_env 1699 INTEGER, INTENT(IN) :: homo 1700 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 1701 REAL(KIND=dp), INTENT(IN) :: omega 1702 1703 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_eps_head_Berry', & 1704 routineP = moduleN//':'//routineN 1705 1706 INTEGER :: col, col_end_in_block, col_offset, col_size, handle, i_col, i_row, ikp, nkp, row, & 1707 row_offset, row_size, row_start_in_block 1708 REAL(KIND=dp) :: abs_k_square, cell_volume, & 1709 correct_kpoint(3), cos_square, & 1710 eigen_diff, relative_kpoint(3), & 1711 sin_square 1712 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: P_head 1713 REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block 1714 TYPE(cell_type), POINTER :: cell 1715 TYPE(dbcsr_iterator_type) :: iter 1716 1717 CALL timeset(routineN, handle) 1718 1719 CALL get_qs_env(qs_env=qs_env, cell=cell) 1720 CALL get_cell(cell=cell, deth=cell_volume) 1721 1722 NULLIFY (data_block) 1723 1724 nkp = kpoints%nkp 1725 1726 ALLOCATE (P_head(nkp)) 1727 P_head(:) = 0.0_dp 1728 1729 ALLOCATE (eps_head(nkp)) 1730 eps_head(:) = 0.0_dp 1731 1732 DO ikp = 1, nkp 1733 1734 relative_kpoint(1:3) = MATMUL(cell%hmat, kpoints%xkp(1:3, ikp)) 1735 1736 correct_kpoint(1:3) = twopi*kpoints%xkp(1:3, ikp) 1737 1738 abs_k_square = (correct_kpoint(1))**2 + (correct_kpoint(2))**2 + (correct_kpoint(3))**2 1739 1740 ! real part of the Berry phase 1741 CALL dbcsr_iterator_start(iter, matrix_berry_re_mo_mo(ikp)%matrix) 1742 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1743 1744 CALL dbcsr_iterator_next_block(iter, row, col, data_block, & 1745 row_size=row_size, col_size=col_size, & 1746 row_offset=row_offset, col_offset=col_offset) 1747 1748 IF (row_offset + row_size <= homo .OR. col_offset > homo) CYCLE 1749 1750 IF (row_offset <= homo) THEN 1751 row_start_in_block = homo - row_offset + 2 1752 ELSE 1753 row_start_in_block = 1 1754 END IF 1755 1756 IF (col_offset + col_size - 1 > homo) THEN 1757 col_end_in_block = homo - col_offset + 1 1758 ELSE 1759 col_end_in_block = col_size 1760 END IF 1761 1762 DO i_row = row_start_in_block, row_size 1763 1764 DO i_col = 1, col_end_in_block 1765 1766 eigen_diff = Eigenval(i_col + col_offset - 1) - Eigenval(i_row + row_offset - 1) 1767 1768 cos_square = (data_block(i_row, i_col))**2 1769 1770 P_head(ikp) = P_head(ikp) + 2.0_dp*eigen_diff/(omega**2 + eigen_diff**2)*cos_square/abs_k_square 1771 1772 END DO 1773 1774 END DO 1775 1776 END DO 1777 1778 CALL dbcsr_iterator_stop(iter) 1779 1780 ! imaginary part of the Berry phase 1781 CALL dbcsr_iterator_start(iter, matrix_berry_im_mo_mo(ikp)%matrix) 1782 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1783 1784 CALL dbcsr_iterator_next_block(iter, row, col, data_block, & 1785 row_size=row_size, col_size=col_size, & 1786 row_offset=row_offset, col_offset=col_offset) 1787 1788 IF (row_offset + row_size <= homo .OR. col_offset > homo) CYCLE 1789 1790 IF (row_offset <= homo) THEN 1791 row_start_in_block = homo - row_offset + 2 1792 ELSE 1793 row_start_in_block = 1 1794 END IF 1795 1796 IF (col_offset + col_size - 1 > homo) THEN 1797 col_end_in_block = homo - col_offset + 1 1798 ELSE 1799 col_end_in_block = col_size 1800 END IF 1801 1802 DO i_row = row_start_in_block, row_size 1803 1804 DO i_col = 1, col_end_in_block 1805 1806 eigen_diff = Eigenval(i_col + col_offset - 1) - Eigenval(i_row + row_offset - 1) 1807 1808 sin_square = (data_block(i_row, i_col))**2 1809 1810 P_head(ikp) = P_head(ikp) + 2.0_dp*eigen_diff/(omega**2 + eigen_diff**2)*sin_square/abs_k_square 1811 1812 END DO 1813 1814 END DO 1815 1816 END DO 1817 1818 CALL dbcsr_iterator_stop(iter) 1819 1820 END DO 1821 1822 CALL mp_sum(P_head, para_env_RPA%group) 1823 1824 ! normalize eps_head 1825 ! 2.0_dp due to closed shell 1826 eps_head(:) = 1.0_dp - 2.0_dp*P_head(:)/cell_volume*fourpi 1827 1828 DEALLOCATE (P_head) 1829 1830 CALL timestop(handle) 1831 1832 END SUBROUTINE compute_eps_head_Berry 1833 1834! ************************************************************************************************** 1835!> \brief ... 1836!> \param qs_env ... 1837!> \param kpoints ... 1838!> \param matrix_berry_re_mo_mo ... 1839!> \param matrix_berry_im_mo_mo ... 1840!> \param fm_mo_coeff ... 1841!> \param para_env ... 1842!> \param do_mo_coeff_Gamma_only ... 1843!> \param homo ... 1844!> \param nmo ... 1845!> \param gw_corr_lev_virt ... 1846!> \param eps_kpoint ... 1847!> \param do_aux_bas ... 1848!> \param frac_aux_mos ... 1849! ************************************************************************************************** 1850 SUBROUTINE get_berry_phase(qs_env, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, fm_mo_coeff, para_env, & 1851 do_mo_coeff_Gamma_only, homo, nmo, gw_corr_lev_virt, eps_kpoint, do_aux_bas, & 1852 frac_aux_mos) 1853 TYPE(qs_environment_type), POINTER :: qs_env 1854 TYPE(kpoint_type), POINTER :: kpoints 1855 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_re_mo_mo, & 1856 matrix_berry_im_mo_mo 1857 TYPE(cp_fm_type), POINTER :: fm_mo_coeff 1858 TYPE(cp_para_env_type), POINTER :: para_env 1859 LOGICAL, INTENT(IN) :: do_mo_coeff_Gamma_only 1860 INTEGER, INTENT(IN) :: homo, nmo, gw_corr_lev_virt 1861 REAL(KIND=dp), INTENT(IN) :: eps_kpoint 1862 LOGICAL, INTENT(IN) :: do_aux_bas 1863 REAL(KIND=dp), INTENT(IN) :: frac_aux_mos 1864 1865 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_berry_phase', & 1866 routineP = moduleN//':'//routineN 1867 1868 INTEGER :: col_index, handle, i_col_local, ikind, & 1869 ikp, nao_aux, ncol_local, nkind, nkp, & 1870 nmo_for_aux_bas 1871 INTEGER, DIMENSION(:), POINTER :: col_indices 1872 REAL(dp) :: abs_kpoint, correct_kpoint(3), & 1873 scale_kpoint 1874 REAL(KIND=dp), DIMENSION(:), POINTER :: evals_P, evals_P_sqrt_inv 1875 TYPE(cell_type), POINTER :: cell 1876 TYPE(cp_fm_struct_type), POINTER :: fm_struct_aux_aux 1877 TYPE(cp_fm_type), POINTER :: fm_mat_eigv_P, fm_mat_P, fm_mat_P_sqrt_inv, & 1878 fm_mat_s_aux_aux_inv, fm_mat_scaled_eigv_P, fm_mat_work_aux_aux 1879 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_aux_aux, & 1880 matrix_s_aux_orb 1881 TYPE(dbcsr_type), POINTER :: cosmat, cosmat_desymm, mat_mo_coeff_aux, mat_mo_coeff_aux_2, & 1882 mat_mo_coeff_Gamma_all, mat_mo_coeff_Gamma_occ_and_GW, mat_mo_coeff_im, mat_mo_coeff_re, & 1883 mat_work_aux_orb, mat_work_aux_orb_2, matrix_P, matrix_P_sqrt, matrix_P_sqrt_inv, & 1884 matrix_s_inv_aux_aux, sinmat, sinmat_desymm, tmp 1885 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: gw_aux_basis_set_list, orb_basis_set_list 1886 TYPE(gto_basis_set_type), POINTER :: basis_set_gw_aux 1887 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1888 POINTER :: sab_orb, sab_orb_mic, sgwgw_list, & 1889 sgworb_list 1890 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1891 TYPE(qs_kind_type), POINTER :: qs_kind 1892 TYPE(qs_ks_env_type), POINTER :: ks_env 1893 1894 CALL timeset(routineN, handle) 1895 1896 nkp = kpoints%nkp 1897 1898 NULLIFY (matrix_berry_re_mo_mo, matrix_s, cell, matrix_berry_im_mo_mo, sinmat, cosmat, tmp, & 1899 cosmat_desymm, sinmat_desymm, qs_kind_set, orb_basis_set_list, sab_orb_mic) 1900 1901 CALL get_qs_env(qs_env=qs_env, & 1902 cell=cell, & 1903 matrix_s=matrix_s, & 1904 qs_kind_set=qs_kind_set, & 1905 nkind=nkind, & 1906 ks_env=ks_env, & 1907 sab_orb=sab_orb) 1908 1909 ALLOCATE (orb_basis_set_list(nkind)) 1910 CALL basis_set_list_setup(orb_basis_set_list, "ORB", qs_kind_set) 1911 1912 CALL setup_neighbor_list(sab_orb_mic, orb_basis_set_list, qs_env=qs_env, mic=.FALSE.) 1913 1914 ! create dbcsr matrix of mo_coeff for multiplcation 1915 NULLIFY (mat_mo_coeff_re) 1916 CALL dbcsr_init_p(mat_mo_coeff_re) 1917 CALL dbcsr_create(matrix=mat_mo_coeff_re, & 1918 template=matrix_s(1)%matrix, & 1919 matrix_type=dbcsr_type_no_symmetry) 1920 1921 NULLIFY (mat_mo_coeff_im) 1922 CALL dbcsr_init_p(mat_mo_coeff_im) 1923 CALL dbcsr_create(matrix=mat_mo_coeff_im, & 1924 template=matrix_s(1)%matrix, & 1925 matrix_type=dbcsr_type_no_symmetry) 1926 1927 NULLIFY (mat_mo_coeff_Gamma_all) 1928 CALL dbcsr_init_p(mat_mo_coeff_Gamma_all) 1929 CALL dbcsr_create(matrix=mat_mo_coeff_Gamma_all, & 1930 template=matrix_s(1)%matrix, & 1931 matrix_type=dbcsr_type_no_symmetry) 1932 1933 CALL copy_fm_to_dbcsr(fm_mo_coeff, mat_mo_coeff_Gamma_all, keep_sparsity=.FALSE.) 1934 1935 NULLIFY (mat_mo_coeff_Gamma_occ_and_GW) 1936 CALL dbcsr_init_p(mat_mo_coeff_Gamma_occ_and_GW) 1937 CALL dbcsr_create(matrix=mat_mo_coeff_Gamma_occ_and_GW, & 1938 template=matrix_s(1)%matrix, & 1939 matrix_type=dbcsr_type_no_symmetry) 1940 1941 CALL copy_fm_to_dbcsr(fm_mo_coeff, mat_mo_coeff_Gamma_occ_and_GW, keep_sparsity=.FALSE.) 1942 1943 IF (.NOT. do_aux_bas) THEN 1944 1945 ! allocate intermediate matrices 1946 CALL dbcsr_init_p(cosmat) 1947 CALL dbcsr_init_p(sinmat) 1948 CALL dbcsr_init_p(tmp) 1949 CALL dbcsr_init_p(cosmat_desymm) 1950 CALL dbcsr_init_p(sinmat_desymm) 1951 CALL dbcsr_create(matrix=cosmat, template=matrix_s(1)%matrix) 1952 CALL dbcsr_create(matrix=sinmat, template=matrix_s(1)%matrix) 1953 CALL dbcsr_create(matrix=tmp, & 1954 template=matrix_s(1)%matrix, & 1955 matrix_type=dbcsr_type_no_symmetry) 1956 CALL dbcsr_create(matrix=cosmat_desymm, & 1957 template=matrix_s(1)%matrix, & 1958 matrix_type=dbcsr_type_no_symmetry) 1959 CALL dbcsr_create(matrix=sinmat_desymm, & 1960 template=matrix_s(1)%matrix, & 1961 matrix_type=dbcsr_type_no_symmetry) 1962 CALL dbcsr_copy(cosmat, matrix_s(1)%matrix) 1963 CALL dbcsr_copy(sinmat, matrix_s(1)%matrix) 1964 CALL dbcsr_set(cosmat, 0.0_dp) 1965 CALL dbcsr_set(sinmat, 0.0_dp) 1966 1967 CALL dbcsr_allocate_matrix_set(matrix_berry_re_mo_mo, nkp) 1968 CALL dbcsr_allocate_matrix_set(matrix_berry_im_mo_mo, nkp) 1969 1970 END IF 1971 1972 IF (do_aux_bas) THEN 1973 1974 NULLIFY (gw_aux_basis_set_list) 1975 ALLOCATE (gw_aux_basis_set_list(nkind)) 1976 1977 DO ikind = 1, nkind 1978 1979 NULLIFY (gw_aux_basis_set_list(ikind)%gto_basis_set) 1980 1981 NULLIFY (basis_set_gw_aux) 1982 1983 qs_kind => qs_kind_set(ikind) 1984 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_gw_aux, basis_type="AUX_GW") 1985 CPASSERT(ASSOCIATED(basis_set_gw_aux)) 1986 1987 basis_set_gw_aux%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius 1988 1989 gw_aux_basis_set_list(ikind)%gto_basis_set => basis_set_gw_aux 1990 1991 END DO 1992 1993 ! neighbor lists 1994 NULLIFY (sgwgw_list, sgworb_list) 1995 CALL setup_neighbor_list(sgwgw_list, gw_aux_basis_set_list, qs_env=qs_env) 1996 CALL setup_neighbor_list(sgworb_list, gw_aux_basis_set_list, orb_basis_set_list, qs_env=qs_env) 1997 1998 NULLIFY (matrix_s_aux_aux, matrix_s_aux_orb) 1999 2000 ! build overlap matrix in gw aux basis and the mixed gw aux basis-orb basis 2001 CALL build_overlap_matrix_simple(ks_env, matrix_s_aux_aux, & 2002 gw_aux_basis_set_list, gw_aux_basis_set_list, sgwgw_list) 2003 2004 CALL build_overlap_matrix_simple(ks_env, matrix_s_aux_orb, & 2005 gw_aux_basis_set_list, orb_basis_set_list, sgworb_list) 2006 2007 CALL dbcsr_get_info(matrix_s_aux_aux(1)%matrix, nfullrows_total=nao_aux) 2008 2009 nmo_for_aux_bas = FLOOR(frac_aux_mos*REAL(nao_aux, KIND=dp)) 2010 2011 CALL cp_fm_struct_create(fm_struct_aux_aux, & 2012 context=fm_mo_coeff%matrix_struct%context, & 2013 nrow_global=nao_aux, & 2014 ncol_global=nao_aux, & 2015 para_env=para_env) 2016 2017 NULLIFY (mat_work_aux_orb) 2018 CALL dbcsr_init_p(mat_work_aux_orb) 2019 CALL dbcsr_create(matrix=mat_work_aux_orb, & 2020 template=matrix_s_aux_orb(1)%matrix, & 2021 matrix_type=dbcsr_type_no_symmetry) 2022 2023 NULLIFY (mat_work_aux_orb_2) 2024 CALL dbcsr_init_p(mat_work_aux_orb_2) 2025 CALL dbcsr_create(matrix=mat_work_aux_orb_2, & 2026 template=matrix_s_aux_orb(1)%matrix, & 2027 matrix_type=dbcsr_type_no_symmetry) 2028 2029 NULLIFY (mat_mo_coeff_aux) 2030 CALL dbcsr_init_p(mat_mo_coeff_aux) 2031 CALL dbcsr_create(matrix=mat_mo_coeff_aux, & 2032 template=matrix_s_aux_orb(1)%matrix, & 2033 matrix_type=dbcsr_type_no_symmetry) 2034 2035 NULLIFY (mat_mo_coeff_aux_2) 2036 CALL dbcsr_init_p(mat_mo_coeff_aux_2) 2037 CALL dbcsr_create(matrix=mat_mo_coeff_aux_2, & 2038 template=matrix_s_aux_orb(1)%matrix, & 2039 matrix_type=dbcsr_type_no_symmetry) 2040 2041 NULLIFY (matrix_s_inv_aux_aux) 2042 CALL dbcsr_init_p(matrix_s_inv_aux_aux) 2043 CALL dbcsr_create(matrix=matrix_s_inv_aux_aux, & 2044 template=matrix_s_aux_aux(1)%matrix, & 2045 matrix_type=dbcsr_type_no_symmetry) 2046 2047 NULLIFY (matrix_P) 2048 CALL dbcsr_init_p(matrix_P) 2049 CALL dbcsr_create(matrix=matrix_P, & 2050 template=matrix_s(1)%matrix, & 2051 matrix_type=dbcsr_type_no_symmetry) 2052 2053 NULLIFY (matrix_P_sqrt) 2054 CALL dbcsr_init_p(matrix_P_sqrt) 2055 CALL dbcsr_create(matrix=matrix_P_sqrt, & 2056 template=matrix_s(1)%matrix, & 2057 matrix_type=dbcsr_type_no_symmetry) 2058 2059 NULLIFY (matrix_P_sqrt_inv) 2060 CALL dbcsr_init_p(matrix_P_sqrt_inv) 2061 CALL dbcsr_create(matrix=matrix_P_sqrt_inv, & 2062 template=matrix_s(1)%matrix, & 2063 matrix_type=dbcsr_type_no_symmetry) 2064 2065 NULLIFY (fm_mat_s_aux_aux_inv) 2066 CALL cp_fm_create(fm_mat_s_aux_aux_inv, fm_struct_aux_aux, name="inverse overlap mat") 2067 2068 NULLIFY (fm_mat_work_aux_aux) 2069 CALL cp_fm_create(fm_mat_work_aux_aux, fm_struct_aux_aux, name="work mat") 2070 2071 NULLIFY (fm_mat_P) 2072 CALL cp_fm_create(fm_mat_P, fm_mo_coeff%matrix_struct) 2073 2074 NULLIFY (fm_mat_eigv_P) 2075 CALL cp_fm_create(fm_mat_eigv_P, fm_mo_coeff%matrix_struct) 2076 2077 NULLIFY (fm_mat_scaled_eigv_P) 2078 CALL cp_fm_create(fm_mat_scaled_eigv_P, fm_mo_coeff%matrix_struct) 2079 2080 NULLIFY (fm_mat_P_sqrt_inv) 2081 CALL cp_fm_create(fm_mat_P_sqrt_inv, fm_mo_coeff%matrix_struct) 2082 2083 NULLIFY (evals_P) 2084 ALLOCATE (evals_P(nmo)) 2085 2086 NULLIFY (evals_P_sqrt_inv) 2087 ALLOCATE (evals_P_sqrt_inv(nmo)) 2088 2089 CALL copy_dbcsr_to_fm(matrix_s_aux_aux(1)%matrix, fm_mat_s_aux_aux_inv) 2090 ! Calculate S_inverse 2091 CALL cp_fm_cholesky_decompose(fm_mat_s_aux_aux_inv) 2092 CALL cp_fm_cholesky_invert(fm_mat_s_aux_aux_inv) 2093 ! Symmetrize the guy 2094 CALL cp_fm_upper_to_full(fm_mat_s_aux_aux_inv, fm_mat_work_aux_aux) 2095 2096 CALL copy_fm_to_dbcsr(fm_mat_s_aux_aux_inv, matrix_s_inv_aux_aux, keep_sparsity=.FALSE.) 2097 2098 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_inv_aux_aux, matrix_s_aux_orb(1)%matrix, 0.0_dp, mat_work_aux_orb, & 2099 filter_eps=1.0E-15_dp) 2100 2101 CALL dbcsr_multiply('N', 'N', 1.0_dp, mat_work_aux_orb, mat_mo_coeff_Gamma_all, 0.0_dp, mat_mo_coeff_aux_2, & 2102 last_column=nmo_for_aux_bas, filter_eps=1.0E-15_dp) 2103 2104 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s_aux_aux(1)%matrix, mat_mo_coeff_aux_2, 0.0_dp, mat_work_aux_orb, & 2105 filter_eps=1.0E-15_dp) 2106 2107 CALL dbcsr_multiply('T', 'N', 1.0_dp, mat_mo_coeff_aux_2, mat_work_aux_orb, 0.0_dp, matrix_P, & 2108 filter_eps=1.0E-15_dp) 2109 2110 CALL copy_dbcsr_to_fm(matrix_P, fm_mat_P) 2111 2112 CALL cp_fm_syevd(fm_mat_P, fm_mat_eigv_P, evals_P) 2113 2114 ! only invert the eigenvalues which correspond to the MOs used in the aux. basis 2115 evals_P_sqrt_inv(1:nmo - nmo_for_aux_bas) = 0.0_dp 2116 evals_P_sqrt_inv(nmo - nmo_for_aux_bas + 1:nmo) = 1.0_dp/SQRT(evals_P(nmo - nmo_for_aux_bas + 1:nmo)) 2117 2118 CALL cp_fm_to_fm(fm_mat_eigv_P, fm_mat_scaled_eigv_P) 2119 2120 CALL cp_fm_get_info(matrix=fm_mat_scaled_eigv_P, & 2121 ncol_local=ncol_local, & 2122 col_indices=col_indices) 2123 2124 ! multiply eigenvectors with inverse sqrt of eigenvalues 2125 DO i_col_local = 1, ncol_local 2126 2127 col_index = col_indices(i_col_local) 2128 2129 fm_mat_scaled_eigv_P%local_data(:, i_col_local) = & 2130 fm_mat_scaled_eigv_P%local_data(:, i_col_local)*evals_P_sqrt_inv(col_index) 2131 2132 END DO 2133 2134 CALL cp_gemm(transa="N", transb="T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, & 2135 matrix_a=fm_mat_eigv_P, matrix_b=fm_mat_scaled_eigv_P, beta=0.0_dp, & 2136 matrix_c=fm_mat_P_sqrt_inv) 2137 2138 CALL copy_fm_to_dbcsr(fm_mat_P_sqrt_inv, matrix_P_sqrt_inv, keep_sparsity=.FALSE.) 2139 2140 CALL dbcsr_multiply('N', 'N', 1.0_dp, mat_mo_coeff_aux_2, matrix_P_sqrt_inv, 0.0_dp, mat_mo_coeff_aux, & 2141 filter_eps=1.0E-15_dp) 2142 2143 ! allocate intermediate matrices 2144 CALL dbcsr_init_p(cosmat) 2145 CALL dbcsr_init_p(sinmat) 2146 CALL dbcsr_init_p(tmp) 2147 CALL dbcsr_init_p(cosmat_desymm) 2148 CALL dbcsr_init_p(sinmat_desymm) 2149 CALL dbcsr_create(matrix=cosmat, template=matrix_s_aux_aux(1)%matrix) 2150 CALL dbcsr_create(matrix=sinmat, template=matrix_s_aux_aux(1)%matrix) 2151 CALL dbcsr_create(matrix=tmp, & 2152 template=matrix_s_aux_orb(1)%matrix, & 2153 matrix_type=dbcsr_type_no_symmetry) 2154 CALL dbcsr_create(matrix=cosmat_desymm, & 2155 template=matrix_s_aux_aux(1)%matrix, & 2156 matrix_type=dbcsr_type_no_symmetry) 2157 CALL dbcsr_create(matrix=sinmat_desymm, & 2158 template=matrix_s_aux_aux(1)%matrix, & 2159 matrix_type=dbcsr_type_no_symmetry) 2160 CALL dbcsr_copy(cosmat, matrix_s_aux_aux(1)%matrix) 2161 CALL dbcsr_copy(sinmat, matrix_s_aux_aux(1)%matrix) 2162 CALL dbcsr_set(cosmat, 0.0_dp) 2163 CALL dbcsr_set(sinmat, 0.0_dp) 2164 2165 CALL dbcsr_allocate_matrix_set(matrix_berry_re_mo_mo, nkp) 2166 CALL dbcsr_allocate_matrix_set(matrix_berry_im_mo_mo, nkp) 2167 2168 ! allocate the new MO coefficients in the aux basis 2169 CALL dbcsr_release_p(mat_mo_coeff_Gamma_all) 2170 CALL dbcsr_release_p(mat_mo_coeff_Gamma_occ_and_GW) 2171 2172 NULLIFY (mat_mo_coeff_Gamma_all) 2173 CALL dbcsr_init_p(mat_mo_coeff_Gamma_all) 2174 CALL dbcsr_create(matrix=mat_mo_coeff_Gamma_all, & 2175 template=matrix_s_aux_orb(1)%matrix, & 2176 matrix_type=dbcsr_type_no_symmetry) 2177 2178 CALL dbcsr_copy(mat_mo_coeff_Gamma_all, mat_mo_coeff_aux) 2179 2180 NULLIFY (mat_mo_coeff_Gamma_occ_and_GW) 2181 CALL dbcsr_init_p(mat_mo_coeff_Gamma_occ_and_GW) 2182 CALL dbcsr_create(matrix=mat_mo_coeff_Gamma_occ_and_GW, & 2183 template=matrix_s_aux_orb(1)%matrix, & 2184 matrix_type=dbcsr_type_no_symmetry) 2185 2186 CALL dbcsr_copy(mat_mo_coeff_Gamma_occ_and_GW, mat_mo_coeff_aux) 2187 2188 DEALLOCATE (evals_P, evals_P_sqrt_inv) 2189 2190 END IF 2191 2192 CALL remove_unnecessary_blocks(mat_mo_coeff_Gamma_occ_and_GW, homo, gw_corr_lev_virt) 2193 2194 DO ikp = 1, nkp 2195 2196 ALLOCATE (matrix_berry_re_mo_mo(ikp)%matrix) 2197 CALL dbcsr_init_p(matrix_berry_re_mo_mo(ikp)%matrix) 2198 CALL dbcsr_create(matrix_berry_re_mo_mo(ikp)%matrix, & 2199 template=matrix_s(1)%matrix, & 2200 matrix_type=dbcsr_type_no_symmetry) 2201 CALL dbcsr_desymmetrize(matrix_s(1)%matrix, matrix_berry_re_mo_mo(ikp)%matrix) 2202 CALL dbcsr_set(matrix_berry_re_mo_mo(ikp)%matrix, 0.0_dp) 2203 2204 ALLOCATE (matrix_berry_im_mo_mo(ikp)%matrix) 2205 CALL dbcsr_init_p(matrix_berry_im_mo_mo(ikp)%matrix) 2206 CALL dbcsr_create(matrix_berry_im_mo_mo(ikp)%matrix, & 2207 template=matrix_s(1)%matrix, & 2208 matrix_type=dbcsr_type_no_symmetry) 2209 CALL dbcsr_desymmetrize(matrix_s(1)%matrix, matrix_berry_im_mo_mo(ikp)%matrix) 2210 CALL dbcsr_set(matrix_berry_im_mo_mo(ikp)%matrix, 0.0_dp) 2211 2212 correct_kpoint(1:3) = -twopi*kpoints%xkp(1:3, ikp) 2213 2214 abs_kpoint = SQRT(correct_kpoint(1)**2 + correct_kpoint(2)**2 + correct_kpoint(3)**2) 2215 2216 IF (abs_kpoint < eps_kpoint) THEN 2217 2218 scale_kpoint = eps_kpoint/abs_kpoint 2219 correct_kpoint(:) = correct_kpoint(:)*scale_kpoint 2220 2221 END IF 2222 2223 ! get the Berry phase 2224 IF (do_aux_bas) THEN 2225 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, correct_kpoint, sab_orb_external=sab_orb_mic, & 2226 basis_type="AUX_GW") 2227 ELSE 2228 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, correct_kpoint, sab_orb_external=sab_orb_mic, & 2229 basis_type="ORB") 2230 END IF 2231 2232 IF (do_mo_coeff_Gamma_only) THEN 2233 2234 CALL dbcsr_desymmetrize(cosmat, cosmat_desymm) 2235 2236 CALL dbcsr_multiply('N', 'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_Gamma_occ_and_GW, 0.0_dp, tmp, & 2237 filter_eps=1.0E-15_dp) 2238 2239 CALL dbcsr_multiply('T', 'N', 1.0_dp, mat_mo_coeff_Gamma_all, tmp, 0.0_dp, & 2240 matrix_berry_re_mo_mo(ikp)%matrix, filter_eps=1.0E-15_dp) 2241 2242 CALL dbcsr_desymmetrize(sinmat, sinmat_desymm) 2243 2244 CALL dbcsr_multiply('N', 'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_Gamma_occ_and_GW, 0.0_dp, tmp, & 2245 filter_eps=1.0E-15_dp) 2246 2247 CALL dbcsr_multiply('T', 'N', 1.0_dp, mat_mo_coeff_Gamma_all, tmp, 0.0_dp, & 2248 matrix_berry_im_mo_mo(ikp)%matrix, filter_eps=1.0E-15_dp) 2249 2250 ELSE 2251 2252 ! get mo coeff at the ikp 2253 CALL copy_fm_to_dbcsr(kpoints%kp_env(ikp)%kpoint_env%mos(1, 1)%mo_set%mo_coeff, & 2254 mat_mo_coeff_re, keep_sparsity=.FALSE.) 2255 2256 CALL copy_fm_to_dbcsr(kpoints%kp_env(ikp)%kpoint_env%mos(2, 1)%mo_set%mo_coeff, & 2257 mat_mo_coeff_im, keep_sparsity=.FALSE.) 2258 2259 CALL dbcsr_desymmetrize(cosmat, cosmat_desymm) 2260 2261 CALL dbcsr_desymmetrize(sinmat, sinmat_desymm) 2262 2263 ! I. 2264 CALL dbcsr_multiply('N', 'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_re, 0.0_dp, tmp) 2265 2266 ! I.1 2267 CALL dbcsr_multiply('T', 'N', 1.0_dp, mat_mo_coeff_Gamma_all, tmp, 0.0_dp, & 2268 matrix_berry_re_mo_mo(ikp)%matrix) 2269 2270 ! II. 2271 CALL dbcsr_multiply('N', 'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_re, 0.0_dp, tmp) 2272 2273 ! II.5 2274 CALL dbcsr_multiply('T', 'N', 1.0_dp, mat_mo_coeff_Gamma_all, tmp, 0.0_dp, & 2275 matrix_berry_im_mo_mo(ikp)%matrix) 2276 2277 ! III. 2278 CALL dbcsr_multiply('N', 'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_im, 0.0_dp, tmp) 2279 2280 ! III.7 2281 CALL dbcsr_multiply('T', 'N', 1.0_dp, mat_mo_coeff_Gamma_all, tmp, 1.0_dp, & 2282 matrix_berry_im_mo_mo(ikp)%matrix) 2283 2284 ! IV. 2285 CALL dbcsr_multiply('N', 'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_im, 0.0_dp, tmp) 2286 2287 ! IV.3 2288 CALL dbcsr_multiply('T', 'N', -1.0_dp, mat_mo_coeff_Gamma_all, tmp, 1.0_dp, & 2289 matrix_berry_re_mo_mo(ikp)%matrix) 2290 2291 END IF 2292 2293 IF (abs_kpoint < eps_kpoint) THEN 2294 2295 CALL dbcsr_scale(matrix_berry_im_mo_mo(ikp)%matrix, 1.0_dp/scale_kpoint) 2296 CALL dbcsr_set(matrix_berry_re_mo_mo(ikp)%matrix, 0.0_dp) 2297 CALL dbcsr_add_on_diag(matrix_berry_re_mo_mo(ikp)%matrix, 1.0_dp) 2298 2299 END IF 2300 2301 END DO 2302 2303 CALL dbcsr_release_p(cosmat) 2304 CALL dbcsr_release_p(sinmat) 2305 CALL dbcsr_release_p(mat_mo_coeff_re) 2306 CALL dbcsr_release_p(mat_mo_coeff_im) 2307 CALL dbcsr_release_p(mat_mo_coeff_Gamma_all) 2308 CALL dbcsr_release_p(mat_mo_coeff_Gamma_occ_and_GW) 2309 CALL dbcsr_release_p(tmp) 2310 CALL dbcsr_release_p(cosmat_desymm) 2311 CALL dbcsr_release_p(sinmat_desymm) 2312 DEALLOCATE (orb_basis_set_list) 2313 2314 CALL release_neighbor_list_sets(sab_orb_mic) 2315 2316 IF (do_aux_bas) THEN 2317 2318 DEALLOCATE (gw_aux_basis_set_list) 2319 CALL dbcsr_deallocate_matrix_set(matrix_s_aux_aux) 2320 CALL dbcsr_deallocate_matrix_set(matrix_s_aux_orb) 2321 CALL dbcsr_release_p(mat_work_aux_orb) 2322 CALL dbcsr_release_p(mat_work_aux_orb_2) 2323 CALL dbcsr_release_p(mat_mo_coeff_aux) 2324 CALL dbcsr_release_p(mat_mo_coeff_aux_2) 2325 CALL dbcsr_release_p(matrix_s_inv_aux_aux) 2326 CALL dbcsr_release_p(matrix_P) 2327 CALL dbcsr_release_p(matrix_P_sqrt) 2328 CALL dbcsr_release_p(matrix_P_sqrt_inv) 2329 2330 CALL cp_fm_struct_release(fm_struct_aux_aux) 2331 2332 CALL cp_fm_release(fm_mat_s_aux_aux_inv) 2333 CALL cp_fm_release(fm_mat_work_aux_aux) 2334 CALL cp_fm_release(fm_mat_P) 2335 CALL cp_fm_release(fm_mat_eigv_P) 2336 CALL cp_fm_release(fm_mat_scaled_eigv_P) 2337 CALL cp_fm_release(fm_mat_P_sqrt_inv) 2338 2339 ! Deallocate the neighbor list structure 2340 CALL release_neighbor_list_sets(sgwgw_list) 2341 CALL release_neighbor_list_sets(sgworb_list) 2342 2343 END IF 2344 2345 CALL timestop(handle) 2346 2347 END SUBROUTINE get_berry_phase 2348 2349! ************************************************************************************************** 2350!> \brief ... 2351!> \param mat_mo_coeff_Gamma_occ_and_GW ... 2352!> \param homo ... 2353!> \param gw_corr_lev_virt ... 2354! ************************************************************************************************** 2355 SUBROUTINE remove_unnecessary_blocks(mat_mo_coeff_Gamma_occ_and_GW, homo, gw_corr_lev_virt) 2356 2357 TYPE(dbcsr_type), POINTER :: mat_mo_coeff_Gamma_occ_and_GW 2358 INTEGER, INTENT(IN) :: homo, gw_corr_lev_virt 2359 2360 INTEGER :: col, col_offset, row 2361 REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block 2362 TYPE(dbcsr_iterator_type) :: iter 2363 2364 CALL dbcsr_iterator_start(iter, mat_mo_coeff_Gamma_occ_and_GW) 2365 2366 DO WHILE (dbcsr_iterator_blocks_left(iter)) 2367 2368 CALL dbcsr_iterator_next_block(iter, row, col, data_block, & 2369 col_offset=col_offset) 2370 2371 IF (col_offset > homo + gw_corr_lev_virt) THEN 2372 2373 data_block = 0.0_dp 2374 2375 END IF 2376 2377 END DO 2378 2379 CALL dbcsr_iterator_stop(iter) 2380 2381 CALL dbcsr_filter(mat_mo_coeff_Gamma_occ_and_GW, 1.0E-15_dp) 2382 2383 END SUBROUTINE remove_unnecessary_blocks 2384 2385! ************************************************************************************************** 2386!> \brief ... 2387!> \param delta_corr ... 2388!> \param eps_inv_head ... 2389!> \param kpoints ... 2390!> \param qs_env ... 2391!> \param matrix_berry_re_mo_mo ... 2392!> \param matrix_berry_im_mo_mo ... 2393!> \param homo ... 2394!> \param gw_corr_lev_occ ... 2395!> \param gw_corr_lev_virt ... 2396!> \param para_env_RPA ... 2397!> \param do_extra_kpoints ... 2398! ************************************************************************************************** 2399 SUBROUTINE kpoint_sum_for_eps_inv_head_Berry(delta_corr, eps_inv_head, kpoints, qs_env, matrix_berry_re_mo_mo, & 2400 matrix_berry_im_mo_mo, homo, gw_corr_lev_occ, gw_corr_lev_virt, & 2401 para_env_RPA, do_extra_kpoints) 2402 2403 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 2404 INTENT(INOUT) :: delta_corr 2405 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 2406 INTENT(IN) :: eps_inv_head 2407 TYPE(kpoint_type), POINTER :: kpoints 2408 TYPE(qs_environment_type), POINTER :: qs_env 2409 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_re_mo_mo, & 2410 matrix_berry_im_mo_mo 2411 INTEGER, INTENT(IN) :: homo, gw_corr_lev_occ, gw_corr_lev_virt 2412 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env_RPA 2413 LOGICAL, INTENT(IN) :: do_extra_kpoints 2414 2415 INTEGER :: col, col_offset, col_size, i_col, i_row, & 2416 ikp, m_level, n_level_gw, nkp, row, & 2417 row_offset, row_size 2418 REAL(KIND=dp) :: abs_k_square, cell_volume, & 2419 check_int_one_over_ksq, contribution, & 2420 weight 2421 REAL(KIND=dp), DIMENSION(3) :: correct_kpoint 2422 REAL(KIND=dp), DIMENSION(:), POINTER :: delta_corr_extra 2423 REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block 2424 TYPE(cell_type), POINTER :: cell 2425 TYPE(dbcsr_iterator_type) :: iter, iter_new 2426 2427 CALL get_qs_env(qs_env=qs_env, cell=cell) 2428 2429 CALL get_cell(cell=cell, deth=cell_volume) 2430 2431 nkp = kpoints%nkp 2432 2433 delta_corr = 0.0_dp 2434 2435 IF (do_extra_kpoints) THEN 2436 NULLIFY (delta_corr_extra) 2437 ALLOCATE (delta_corr_extra(1 + homo - gw_corr_lev_occ:homo + gw_corr_lev_virt)) 2438 delta_corr_extra = 0.0_dp 2439 END IF 2440 2441 check_int_one_over_ksq = 0.0_dp 2442 2443 DO ikp = 1, nkp 2444 2445 weight = kpoints%wkp(ikp) 2446 2447 correct_kpoint(1:3) = twopi*kpoints%xkp(1:3, ikp) 2448 2449 abs_k_square = (correct_kpoint(1))**2 + (correct_kpoint(2))**2 + (correct_kpoint(3))**2 2450 2451 ! cos part of the Berry phase 2452 CALL dbcsr_iterator_start(iter, matrix_berry_re_mo_mo(ikp)%matrix) 2453 DO WHILE (dbcsr_iterator_blocks_left(iter)) 2454 2455 CALL dbcsr_iterator_next_block(iter, row, col, data_block, & 2456 row_size=row_size, col_size=col_size, & 2457 row_offset=row_offset, col_offset=col_offset) 2458 2459 DO i_col = 1, col_size 2460 2461 DO n_level_gw = 1 + homo - gw_corr_lev_occ, homo + gw_corr_lev_virt 2462 2463 IF (n_level_gw == i_col + col_offset - 1) THEN 2464 2465 DO i_row = 1, row_size 2466 2467 contribution = weight*(eps_inv_head(ikp) - 1.0_dp)/abs_k_square*(data_block(i_row, i_col))**2 2468 2469 m_level = i_row + row_offset - 1 2470 2471 ! we only compute the correction for n=m 2472 IF (m_level .NE. n_level_gw) CYCLE 2473 2474 IF (.NOT. do_extra_kpoints) THEN 2475 2476 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution 2477 2478 ELSE 2479 2480 IF (ikp <= nkp*8/9) THEN 2481 2482 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution 2483 2484 ELSE 2485 2486 delta_corr_extra(n_level_gw) = delta_corr_extra(n_level_gw) + contribution 2487 2488 END IF 2489 2490 END IF 2491 2492 END DO 2493 2494 END IF 2495 2496 END DO 2497 2498 END DO 2499 2500 END DO 2501 2502 CALL dbcsr_iterator_stop(iter) 2503 2504 ! the same for the im. part of the Berry phase 2505 CALL dbcsr_iterator_start(iter_new, matrix_berry_im_mo_mo(ikp)%matrix) 2506 DO WHILE (dbcsr_iterator_blocks_left(iter_new)) 2507 2508 CALL dbcsr_iterator_next_block(iter_new, row, col, data_block, & 2509 row_size=row_size, col_size=col_size, & 2510 row_offset=row_offset, col_offset=col_offset) 2511 2512 DO i_col = 1, col_size 2513 2514 DO n_level_gw = 1 + homo - gw_corr_lev_occ, homo + gw_corr_lev_virt 2515 2516 IF (n_level_gw == i_col + col_offset - 1) THEN 2517 2518 DO i_row = 1, row_size 2519 2520 m_level = i_row + row_offset - 1 2521 2522 contribution = weight*(eps_inv_head(ikp) - 1.0_dp)/abs_k_square*(data_block(i_row, i_col))**2 2523 2524 ! we only compute the correction for n=m 2525 IF (m_level .NE. n_level_gw) CYCLE 2526 2527 IF (.NOT. do_extra_kpoints) THEN 2528 2529 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution 2530 2531 ELSE 2532 2533 IF (ikp <= nkp*8/9) THEN 2534 2535 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution 2536 2537 ELSE 2538 2539 delta_corr_extra(n_level_gw) = delta_corr_extra(n_level_gw) + contribution 2540 2541 END IF 2542 2543 END IF 2544 2545 END DO 2546 2547 END IF 2548 2549 END DO 2550 2551 END DO 2552 2553 END DO 2554 2555 CALL dbcsr_iterator_stop(iter_new) 2556 2557 check_int_one_over_ksq = check_int_one_over_ksq + weight/abs_k_square 2558 2559 END DO 2560 2561 ! normalize by the cell volume 2562 delta_corr = delta_corr/cell_volume*fourpi 2563 2564 check_int_one_over_ksq = check_int_one_over_ksq/cell_volume 2565 2566 CALL mp_sum(delta_corr, para_env_RPA%group) 2567 2568 IF (do_extra_kpoints) THEN 2569 2570 delta_corr_extra = delta_corr_extra/cell_volume*fourpi 2571 2572 CALL mp_sum(delta_corr_extra, para_env_RPA%group) 2573 2574 delta_corr(:) = delta_corr(:) + (delta_corr(:) - delta_corr_extra(:)) 2575 2576 DEALLOCATE (delta_corr_extra) 2577 2578 END IF 2579 2580 END SUBROUTINE kpoint_sum_for_eps_inv_head_Berry 2581 2582! ************************************************************************************************** 2583!> \brief ... 2584!> \param eps_inv_head ... 2585!> \param eps_head ... 2586!> \param kpoints ... 2587! ************************************************************************************************** 2588 SUBROUTINE compute_eps_inv_head(eps_inv_head, eps_head, kpoints) 2589 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 2590 INTENT(OUT) :: eps_inv_head 2591 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 2592 INTENT(IN) :: eps_head 2593 TYPE(kpoint_type), POINTER :: kpoints 2594 2595 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_eps_inv_head', & 2596 routineP = moduleN//':'//routineN 2597 2598 INTEGER :: handle, ikp, nkp 2599 2600 CALL timeset(routineN, handle) 2601 2602 nkp = kpoints%nkp 2603 2604 ALLOCATE (eps_inv_head(nkp)) 2605 2606 DO ikp = 1, nkp 2607 2608 eps_inv_head(ikp) = 1.0_dp/eps_head(ikp) 2609 2610 END DO 2611 2612 CALL timestop(handle) 2613 2614 END SUBROUTINE compute_eps_inv_head 2615 2616! ************************************************************************************************** 2617!> \brief ... 2618!> \param qs_env ... 2619!> \param kpoints ... 2620!> \param kp_grid ... 2621!> \param num_kp_grids ... 2622!> \param para_env ... 2623!> \param h_inv ... 2624!> \param nmo ... 2625!> \param do_mo_coeff_Gamma_only ... 2626!> \param do_extra_kpoints ... 2627! ************************************************************************************************** 2628 SUBROUTINE get_kpoints(qs_env, kpoints, kp_grid, num_kp_grids, para_env, h_inv, nmo, & 2629 do_mo_coeff_Gamma_only, do_extra_kpoints) 2630 TYPE(qs_environment_type), POINTER :: qs_env 2631 TYPE(kpoint_type), POINTER :: kpoints 2632 INTEGER, DIMENSION(:), POINTER :: kp_grid 2633 INTEGER, INTENT(IN) :: num_kp_grids 2634 TYPE(cp_para_env_type), POINTER :: para_env 2635 REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_inv 2636 INTEGER, INTENT(IN) :: nmo 2637 LOGICAL, INTENT(IN) :: do_mo_coeff_Gamma_only, do_extra_kpoints 2638 2639 INTEGER :: end_kp, i, i_grid_level, ix, iy, iz, & 2640 nkp_inner_grid, nkp_outer_grid, & 2641 npoints, start_kp 2642 INTEGER, DIMENSION(3) :: outer_kp_grid 2643 REAL(KIND=dp) :: kpoint_weight_left, single_weight 2644 REAL(KIND=dp), DIMENSION(3) :: kpt_latt, reducing_factor 2645 TYPE(cell_type), POINTER :: cell 2646 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 2647 TYPE(qs_environment_type), POINTER :: qs_env_kp_Gamma_only 2648 2649 NULLIFY (kpoints, cell, particle_set, qs_env_kp_Gamma_only) 2650 2651 ! check whether kp_grid includes the Gamma point. If so, abort. 2652 CPASSERT(MOD(kp_grid(1)*kp_grid(2)*kp_grid(3), 2) == 0) 2653 IF (do_extra_kpoints) THEN 2654 CPASSERT(do_mo_coeff_Gamma_only) 2655 END IF 2656 2657 IF (do_mo_coeff_Gamma_only) THEN 2658 2659 outer_kp_grid(1) = kp_grid(1) - 1 2660 outer_kp_grid(2) = kp_grid(2) - 1 2661 outer_kp_grid(3) = kp_grid(3) - 1 2662 2663 CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set) 2664 2665 CALL get_cell(cell, h_inv=h_inv) 2666 2667 CALL kpoint_create(kpoints) 2668 2669 kpoints%kp_scheme = "GENERAL" 2670 kpoints%symmetry = .FALSE. 2671 kpoints%verbose = .FALSE. 2672 kpoints%full_grid = .FALSE. 2673 kpoints%use_real_wfn = .FALSE. 2674 kpoints%eps_geo = 1.e-6_dp 2675 npoints = kp_grid(1)*kp_grid(2)*kp_grid(3)/2 + & 2676 (num_kp_grids - 1)*((outer_kp_grid(1) + 1)/2*outer_kp_grid(2)*outer_kp_grid(3) - 1) 2677 2678 IF (do_extra_kpoints) THEN 2679 2680 CPASSERT(num_kp_grids == 1) 2681 CPASSERT(MOD(kp_grid(1), 4) == 0) 2682 CPASSERT(MOD(kp_grid(2), 4) == 0) 2683 CPASSERT(MOD(kp_grid(3), 4) == 0) 2684 2685 END IF 2686 2687 IF (do_extra_kpoints) THEN 2688 2689 npoints = kp_grid(1)*kp_grid(2)*kp_grid(3)/2 + kp_grid(1)*kp_grid(2)*kp_grid(3)/2/8 2690 2691 END IF 2692 2693 kpoints%full_grid = .TRUE. 2694 kpoints%nkp = npoints 2695 ALLOCATE (kpoints%xkp(3, npoints), kpoints%wkp(npoints)) 2696 kpoints%xkp = 0.0_dp 2697 kpoints%wkp = 0.0_dp 2698 2699 nkp_outer_grid = outer_kp_grid(1)*outer_kp_grid(2)*outer_kp_grid(3) 2700 nkp_inner_grid = kp_grid(1)*kp_grid(2)*kp_grid(3) 2701 2702 i = 0 2703 reducing_factor(:) = 1.0_dp 2704 kpoint_weight_left = 1.0_dp 2705 2706 ! the outer grids 2707 DO i_grid_level = 1, num_kp_grids - 1 2708 2709 single_weight = kpoint_weight_left/REAL(nkp_outer_grid, KIND=dp) 2710 2711 start_kp = i + 1 2712 2713 DO ix = 1, outer_kp_grid(1) 2714 DO iy = 1, outer_kp_grid(2) 2715 DO iz = 1, outer_kp_grid(3) 2716 2717 ! exclude Gamma 2718 IF (2*ix - outer_kp_grid(1) - 1 == 0 .AND. 2*iy - outer_kp_grid(2) - 1 == 0 .AND. & 2719 2*iz - outer_kp_grid(3) - 1 == 0) CYCLE 2720 2721 ! use time reversal symmetry k<->-k 2722 IF (2*ix - outer_kp_grid(1) - 1 < 0) CYCLE 2723 2724 i = i + 1 2725 kpt_latt(1) = REAL(2*ix - outer_kp_grid(1) - 1, KIND=dp)/(2._dp*REAL(outer_kp_grid(1), KIND=dp)) & 2726 *reducing_factor(1) 2727 kpt_latt(2) = REAL(2*iy - outer_kp_grid(2) - 1, KIND=dp)/(2._dp*REAL(outer_kp_grid(2), KIND=dp)) & 2728 *reducing_factor(2) 2729 kpt_latt(3) = REAL(2*iz - outer_kp_grid(3) - 1, KIND=dp)/(2._dp*REAL(outer_kp_grid(3), KIND=dp)) & 2730 *reducing_factor(3) 2731 kpoints%xkp(1:3, i) = MATMUL(TRANSPOSE(h_inv), kpt_latt(:)) 2732 2733 IF (2*ix - outer_kp_grid(1) - 1 == 0) THEN 2734 kpoints%wkp(i) = single_weight 2735 ELSE 2736 kpoints%wkp(i) = 2._dp*single_weight 2737 END IF 2738 2739 END DO 2740 END DO 2741 END DO 2742 2743 end_kp = i 2744 2745 kpoint_weight_left = kpoint_weight_left - SUM(kpoints%wkp(start_kp:end_kp)) 2746 2747 reducing_factor(1) = reducing_factor(1)/REAL(outer_kp_grid(1), KIND=dp) 2748 reducing_factor(2) = reducing_factor(2)/REAL(outer_kp_grid(2), KIND=dp) 2749 reducing_factor(3) = reducing_factor(3)/REAL(outer_kp_grid(3), KIND=dp) 2750 2751 END DO 2752 2753 single_weight = kpoint_weight_left/REAL(nkp_inner_grid, KIND=dp) 2754 2755 ! the inner grid 2756 DO ix = 1, kp_grid(1) 2757 DO iy = 1, kp_grid(2) 2758 DO iz = 1, kp_grid(3) 2759 2760 ! use time reversal symmetry k<->-k 2761 IF (2*ix - kp_grid(1) - 1 < 0) CYCLE 2762 2763 i = i + 1 2764 kpt_latt(1) = REAL(2*ix - kp_grid(1) - 1, KIND=dp)/(2._dp*REAL(kp_grid(1), KIND=dp))*reducing_factor(1) 2765 kpt_latt(2) = REAL(2*iy - kp_grid(2) - 1, KIND=dp)/(2._dp*REAL(kp_grid(2), KIND=dp))*reducing_factor(2) 2766 kpt_latt(3) = REAL(2*iz - kp_grid(3) - 1, KIND=dp)/(2._dp*REAL(kp_grid(3), KIND=dp))*reducing_factor(3) 2767 2768 kpoints%xkp(1:3, i) = MATMUL(TRANSPOSE(h_inv), kpt_latt(:)) 2769 2770 kpoints%wkp(i) = 2._dp*single_weight 2771 2772 END DO 2773 END DO 2774 END DO 2775 2776 IF (do_extra_kpoints) THEN 2777 2778 single_weight = kpoint_weight_left/REAL(kp_grid(1)*kp_grid(2)*kp_grid(3)/8, KIND=dp) 2779 2780 DO ix = 1, kp_grid(1)/2 2781 DO iy = 1, kp_grid(2)/2 2782 DO iz = 1, kp_grid(3)/2 2783 2784 ! use time reversal symmetry k<->-k 2785 IF (2*ix - kp_grid(1)/2 - 1 < 0) CYCLE 2786 2787 i = i + 1 2788 kpt_latt(1) = REAL(2*ix - kp_grid(1)/2 - 1, KIND=dp)/(REAL(kp_grid(1), KIND=dp)) 2789 kpt_latt(2) = REAL(2*iy - kp_grid(2)/2 - 1, KIND=dp)/(REAL(kp_grid(2), KIND=dp)) 2790 kpt_latt(3) = REAL(2*iz - kp_grid(3)/2 - 1, KIND=dp)/(REAL(kp_grid(3), KIND=dp)) 2791 2792 kpoints%xkp(1:3, i) = MATMUL(TRANSPOSE(h_inv), kpt_latt(:)) 2793 2794 kpoints%wkp(i) = 2._dp*single_weight 2795 2796 END DO 2797 END DO 2798 END DO 2799 2800 END IF 2801 2802 ! default: no symmetry settings 2803 ALLOCATE (kpoints%kp_sym(kpoints%nkp)) 2804 DO i = 1, kpoints%nkp 2805 NULLIFY (kpoints%kp_sym(i)%kpoint_sym) 2806 CALL kpoint_sym_create(kpoints%kp_sym(i)%kpoint_sym) 2807 END DO 2808 2809 ELSE 2810 2811 CALL create_kp_from_gamma(qs_env, qs_env_kp_Gamma_only) 2812 2813 CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set) 2814 2815 CALL calculate_kp_orbitals(qs_env_kp_Gamma_only, kpoints, "MONKHORST-PACK", nadd=nmo, mp_grid=kp_grid(1:3), & 2816 group_size_ext=para_env%num_pe) 2817 2818 CALL qs_env_release(qs_env_kp_Gamma_only) 2819 2820 END IF 2821 2822 END SUBROUTINE get_kpoints 2823 2824! ************************************************************************************************** 2825!> \brief ... 2826!> \param vec_Sigma_c_gw ... 2827!> \param Eigenval_DFT ... 2828!> \param eps_eigenval ... 2829! ************************************************************************************************** 2830 SUBROUTINE average_degenerate_levels(vec_Sigma_c_gw, Eigenval_DFT, eps_eigenval) 2831 COMPLEX(KIND=dp), ALLOCATABLE, & 2832 DIMENSION(:, :, :), INTENT(INOUT) :: vec_Sigma_c_gw 2833 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval_DFT 2834 REAL(KIND=dp), INTENT(IN) :: eps_eigenval 2835 2836 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: avg_self_energy 2837 INTEGER :: degeneracy, first_degenerate_level, i_deg_level, i_level_gw, j_deg_level, jquad, & 2838 num_deg_levels, num_integ_points, num_levels_gw 2839 INTEGER, ALLOCATABLE, DIMENSION(:) :: list_degenerate_levels 2840 2841 num_levels_gw = SIZE(vec_Sigma_c_gw, 1) 2842 2843 ALLOCATE (list_degenerate_levels(num_levels_gw)) 2844 list_degenerate_levels = 1 2845 2846 num_integ_points = SIZE(vec_Sigma_c_gw, 2) 2847 2848 ALLOCATE (avg_self_energy(num_integ_points)) 2849 2850 DO i_level_gw = 2, num_levels_gw 2851 2852 IF (ABS(Eigenval_DFT(i_level_gw) - Eigenval_DFT(i_level_gw - 1)) < eps_eigenval) THEN 2853 2854 list_degenerate_levels(i_level_gw) = list_degenerate_levels(i_level_gw - 1) 2855 2856 ELSE 2857 2858 list_degenerate_levels(i_level_gw) = list_degenerate_levels(i_level_gw - 1) + 1 2859 2860 END IF 2861 2862 END DO 2863 2864 num_deg_levels = list_degenerate_levels(num_levels_gw) 2865 2866 DO i_deg_level = 1, num_deg_levels 2867 2868 degeneracy = 0 2869 2870 DO i_level_gw = 1, num_levels_gw 2871 2872 IF (degeneracy == 0 .AND. i_deg_level == list_degenerate_levels(i_level_gw)) THEN 2873 2874 first_degenerate_level = i_level_gw 2875 2876 END IF 2877 2878 IF (i_deg_level == list_degenerate_levels(i_level_gw)) THEN 2879 2880 degeneracy = degeneracy + 1 2881 2882 END IF 2883 2884 END DO 2885 2886 DO jquad = 1, num_integ_points 2887 2888 avg_self_energy(jquad) = SUM(vec_Sigma_c_gw(first_degenerate_level:first_degenerate_level + degeneracy - 1, jquad, 1)) & 2889 /REAL(degeneracy, KIND=dp) 2890 2891 END DO 2892 2893 DO j_deg_level = 0, degeneracy - 1 2894 2895 vec_Sigma_c_gw(first_degenerate_level + j_deg_level, :, 1) = avg_self_energy(:) 2896 2897 END DO 2898 2899 END DO 2900 2901 END SUBROUTINE average_degenerate_levels 2902 2903! ************************************************************************************************** 2904!> \brief ... 2905!> \param Eigenval ... 2906!> \param Eigenval_scf ... 2907!> \param Eigenval_kp ... 2908!> \param Eigenval_scf_kp ... 2909!> \param kpoints ... 2910!> \param ikp ... 2911!> \param my_open_shell ... 2912! ************************************************************************************************** 2913 SUBROUTINE get_eigenval_for_conti(Eigenval, Eigenval_scf, Eigenval_kp, Eigenval_scf_kp, kpoints, ikp, my_open_shell) 2914 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: Eigenval, Eigenval_scf 2915 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 2916 INTENT(INOUT) :: Eigenval_kp, Eigenval_scf_kp 2917 TYPE(kpoint_type), POINTER :: kpoints 2918 INTEGER, INTENT(IN) :: ikp 2919 LOGICAL, INTENT(IN) :: my_open_shell 2920 2921 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_eigenval_for_conti', & 2922 routineP = moduleN//':'//routineN 2923 2924 INTEGER :: handle, ispin, jkp, nkp, nmo, nspin 2925 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues 2926 TYPE(mo_set_type), POINTER :: mo_set 2927 2928 CALL timeset(routineN, handle) 2929 2930 CALL get_kpoint_info(kpoints, nkp=nkp) 2931 2932 nmo = SIZE(Eigenval) 2933 2934 IF (my_open_shell) THEN 2935 nspin = 2 2936 ELSE 2937 nspin = 1 2938 END IF 2939 2940 IF (ikp == 1) THEN 2941 ALLOCATE (Eigenval_kp(SIZE(Eigenval), nkp)) 2942 ALLOCATE (Eigenval_scf_kp(SIZE(Eigenval), nkp)) 2943 2944 DO jkp = 1, nkp 2945 2946 DO ispin = 1, nspin 2947 2948 mo_set => kpoints%kp_env(jkp)%kpoint_env%mos(1, ispin)%mo_set 2949 2950 CALL get_mo_set(mo_set=mo_set, eigenvalues=mo_eigenvalues) 2951 2952 Eigenval_kp(1:nmo, jkp) = mo_eigenvalues(1:nmo) 2953 Eigenval_scf_kp(1:nmo, jkp) = mo_eigenvalues(1:nmo) 2954 2955 END DO 2956 2957 END DO 2958 2959 END IF 2960 2961 Eigenval(1:nmo) = Eigenval_kp(1:nmo, ikp) 2962 Eigenval_scf(1:nmo) = Eigenval_scf_kp(1:nmo, ikp) 2963 2964 CALL timestop(handle) 2965 2966 END SUBROUTINE 2967 2968! ************************************************************************************************** 2969!> \brief ... 2970!> \param vec_gw_energ ... 2971!> \param vec_omega_fit_gw ... 2972!> \param z_value ... 2973!> \param m_value ... 2974!> \param vec_Sigma_c_gw ... 2975!> \param vec_Sigma_x_minus_vxc_gw ... 2976!> \param Eigenval ... 2977!> \param Eigenval_scf ... 2978!> \param n_level_gw ... 2979!> \param gw_corr_lev_occ ... 2980!> \param num_poles ... 2981!> \param num_fit_points ... 2982!> \param crossing_search ... 2983!> \param homo ... 2984!> \param stop_crit ... 2985!> \param fermi_level_offset ... 2986!> \param do_gw_im_time ... 2987! ************************************************************************************************** 2988 SUBROUTINE fit_and_continuation_2pole(vec_gw_energ, vec_omega_fit_gw, & 2989 z_value, m_value, vec_Sigma_c_gw, vec_Sigma_x_minus_vxc_gw, & 2990 Eigenval, Eigenval_scf, n_level_gw, gw_corr_lev_occ, num_poles, & 2991 num_fit_points, crossing_search, homo, stop_crit, & 2992 fermi_level_offset, do_gw_im_time) 2993 2994 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 2995 INTENT(INOUT) :: vec_gw_energ, vec_omega_fit_gw, z_value, & 2996 m_value 2997 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vec_Sigma_c_gw 2998 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec_Sigma_x_minus_vxc_gw, Eigenval, & 2999 Eigenval_scf 3000 INTEGER, INTENT(IN) :: n_level_gw, gw_corr_lev_occ, num_poles, & 3001 num_fit_points, crossing_search, homo 3002 REAL(KIND=dp), INTENT(IN) :: stop_crit, fermi_level_offset 3003 LOGICAL, INTENT(IN) :: do_gw_im_time 3004 3005 CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_and_continuation_2pole', & 3006 routineP = moduleN//':'//routineN 3007 3008 COMPLEX(KIND=dp) :: func_val, im_unit, one, re_unit, rho1, & 3009 zero 3010 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dLambda, dLambda_2, Lambda, & 3011 Lambda_without_offset, vec_b_gw, & 3012 vec_b_gw_copy 3013 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A_gw, mat_B_gw 3014 INTEGER :: handle4, ierr, iii, iiter, info, & 3015 integ_range, jjj, jquad, kkk, & 3016 max_iter_fit, n_level_gw_ref, num_var, & 3017 xpos 3018 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv 3019 LOGICAL :: could_exit 3020 REAL(KIND=dp) :: chi2, chi2_old, delta, deriv_val_real, e_fermi, gw_energ, Ldown, & 3021 level_energ_GW, Lup, range_step, ScalParam, sign_occ_virt, stat_error 3022 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Lambda_Im, Lambda_Re, stat_errors, & 3023 vec_N_gw, vec_omega_fit_gw_sign 3024 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_N_gw 3025 3026 im_unit = (0.0_dp, 1.0_dp) 3027 re_unit = (1.0_dp, 0.0_dp) 3028 3029 max_iter_fit = 10000 3030 3031 num_var = 2*num_poles + 1 3032 ALLOCATE (Lambda(num_var)) 3033 Lambda = (0.0_dp, 0.0_dp) 3034 ALLOCATE (Lambda_without_offset(num_var)) 3035 Lambda_without_offset = (0.0_dp, 0.0_dp) 3036 ALLOCATE (Lambda_Re(num_var)) 3037 Lambda_Re = 0.0_dp 3038 ALLOCATE (Lambda_Im(num_var)) 3039 Lambda_Im = 0.0_dp 3040 3041 ALLOCATE (vec_omega_fit_gw_sign(num_fit_points)) 3042 3043 IF (n_level_gw <= gw_corr_lev_occ) THEN 3044 sign_occ_virt = -1.0_dp 3045 ELSE 3046 sign_occ_virt = 1.0_dp 3047 END IF 3048 3049 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 3050 3051 DO jquad = 1, num_fit_points 3052 vec_omega_fit_gw_sign(jquad) = ABS(vec_omega_fit_gw(jquad))*sign_occ_virt 3053 END DO 3054 3055 ! initial guess 3056 range_step = (vec_omega_fit_gw_sign(num_fit_points) - vec_omega_fit_gw_sign(1))/(num_poles - 1) 3057 DO iii = 1, num_poles 3058 Lambda_Im(2*iii + 1) = vec_omega_fit_gw_sign(1) + (iii - 1)*range_step 3059 END DO 3060 range_step = (vec_omega_fit_gw_sign(num_fit_points) - vec_omega_fit_gw_sign(1))/num_poles 3061 DO iii = 1, num_poles 3062 Lambda_Re(2*iii + 1) = ABS(vec_omega_fit_gw_sign(1) + (iii - 0.5_dp)*range_step) 3063 END DO 3064 3065 DO iii = 1, num_var 3066 Lambda(iii) = Lambda_Re(iii) + im_unit*Lambda_Im(iii) 3067 END DO 3068 3069 CALL calc_chi2(chi2_old, Lambda, vec_Sigma_c_gw, vec_omega_fit_gw_sign, num_poles, & 3070 num_fit_points, n_level_gw) 3071 3072 ALLOCATE (mat_A_gw(num_poles + 1, num_poles + 1)) 3073 ALLOCATE (vec_b_gw(num_poles + 1)) 3074 ALLOCATE (ipiv(num_poles + 1)) 3075 mat_A_gw = (0.0_dp, 0.0_dp) 3076 vec_b_gw = 0.0_dp 3077 3078 DO iii = 1, num_poles + 1 3079 mat_A_gw(iii, 1) = (1.0_dp, 0.0_dp) 3080 END DO 3081 integ_range = num_fit_points/num_poles 3082 DO kkk = 1, num_poles + 1 3083 xpos = (kkk - 1)*integ_range + 1 3084 xpos = MIN(xpos, num_fit_points) 3085 ! calculate coefficient at this point 3086 DO iii = 1, num_poles 3087 jjj = iii*2 3088 func_val = (1.0_dp, 0.0_dp)/(im_unit*vec_omega_fit_gw_sign(xpos) - & 3089 CMPLX(Lambda_Re(jjj + 1), Lambda_Im(jjj + 1), KIND=dp)) 3090 mat_A_gw(kkk, iii + 1) = func_val 3091 END DO 3092 vec_b_gw(kkk) = vec_Sigma_c_gw(n_level_gw, xpos) 3093 END DO 3094 3095 ! Solve system of linear equations 3096 CALL ZGETRF(num_poles + 1, num_poles + 1, mat_A_gw, num_poles + 1, ipiv, info) 3097 3098 CALL ZGETRS('N', num_poles + 1, 1, mat_A_gw, num_poles + 1, ipiv, vec_b_gw, num_poles + 1, info) 3099 3100 Lambda_Re(1) = REAL(vec_b_gw(1)) 3101 Lambda_Im(1) = AIMAG(vec_b_gw(1)) 3102 DO iii = 1, num_poles 3103 jjj = iii*2 3104 Lambda_Re(jjj) = REAL(vec_b_gw(iii + 1)) 3105 Lambda_Im(jjj) = AIMAG(vec_b_gw(iii + 1)) 3106 END DO 3107 3108 DEALLOCATE (mat_A_gw) 3109 DEALLOCATE (vec_b_gw) 3110 DEALLOCATE (ipiv) 3111 3112 ALLOCATE (mat_A_gw(num_var*2, num_var*2)) 3113 ALLOCATE (mat_B_gw(num_fit_points, num_var*2)) 3114 ALLOCATE (dLambda(num_fit_points)) 3115 ALLOCATE (dLambda_2(num_fit_points)) 3116 ALLOCATE (vec_b_gw(num_var*2)) 3117 ALLOCATE (vec_b_gw_copy(num_var*2)) 3118 ALLOCATE (ipiv(num_var*2)) 3119 3120 ScalParam = 0.01_dp 3121 Ldown = 1.5_dp 3122 Lup = 10.0_dp 3123 could_exit = .FALSE. 3124 3125 ! iteration loop for fitting 3126 DO iiter = 1, max_iter_fit 3127 3128 CALL timeset(routineN//"_fit_loop_1", handle4) 3129 3130 ! calc delta lambda 3131 DO iii = 1, num_var 3132 Lambda(iii) = Lambda_Re(iii) + im_unit*Lambda_Im(iii) 3133 END DO 3134 dLambda = (0.0_dp, 0.0_dp) 3135 3136 DO kkk = 1, num_fit_points 3137 func_val = Lambda(1) 3138 DO iii = 1, num_poles 3139 jjj = iii*2 3140 func_val = func_val + Lambda(jjj)/(vec_omega_fit_gw_sign(kkk)*im_unit - Lambda(jjj + 1)) 3141 END DO 3142 dLambda(kkk) = vec_Sigma_c_gw(n_level_gw, kkk) - func_val 3143 END DO 3144 rho1 = SUM(dLambda*dLambda) 3145 3146 ! fill matrix 3147 mat_B_gw = (0.0_dp, 0.0_dp) 3148 DO iii = 1, num_fit_points 3149 mat_B_gw(iii, 1) = 1.0_dp 3150 mat_B_gw(iii, num_var + 1) = im_unit 3151 END DO 3152 DO iii = 1, num_poles 3153 jjj = iii*2 3154 DO kkk = 1, num_fit_points 3155 mat_B_gw(kkk, jjj) = 1.0_dp/(im_unit*vec_omega_fit_gw_sign(kkk) - Lambda(jjj + 1)) 3156 mat_B_gw(kkk, jjj + num_var) = im_unit/(im_unit*vec_omega_fit_gw_sign(kkk) - Lambda(jjj + 1)) 3157 mat_B_gw(kkk, jjj + 1) = Lambda(jjj)/(im_unit*vec_omega_fit_gw_sign(kkk) - Lambda(jjj + 1))**2 3158 mat_B_gw(kkk, jjj + 1 + num_var) = (-Lambda_Im(jjj) + im_unit*Lambda_Re(jjj))/ & 3159 (im_unit*vec_omega_fit_gw_sign(kkk) - Lambda(jjj + 1))**2 3160 END DO 3161 END DO 3162 3163 CALL timestop(handle4) 3164 3165 CALL timeset(routineN//"_fit_matmul_1", handle4) 3166 3167 one = (1.0_dp, 0.0_dp) 3168 zero = (0.0_dp, 0.0_dp) 3169 CALL zgemm('C', 'N', num_var*2, num_var*2, num_fit_points, one, mat_B_gw, num_fit_points, mat_B_gw, num_fit_points, & 3170 zero, mat_A_gw, num_var*2) 3171 CALL timestop(handle4) 3172 3173 CALL timeset(routineN//"_fit_zgemv_1", handle4) 3174 CALL zgemv('C', num_fit_points, num_var*2, one, mat_B_gw, num_fit_points, dLambda, 1, & 3175 zero, vec_b_gw, 1) 3176 3177 CALL timestop(handle4) 3178 3179 ! scale diagonal elements of a_mat 3180 DO iii = 1, num_var*2 3181 mat_A_gw(iii, iii) = mat_A_gw(iii, iii) + ScalParam*mat_A_gw(iii, iii) 3182 END DO 3183 3184 ! solve linear system 3185 ierr = 0 3186 ipiv = 0 3187 3188 CALL timeset(routineN//"_fit_lin_eq_2", handle4) 3189 3190 CALL ZGETRF(2*num_var, 2*num_var, mat_A_gw, 2*num_var, ipiv, info) 3191 3192 CALL ZGETRS('N', 2*num_var, 1, mat_A_gw, 2*num_var, ipiv, vec_b_gw, 2*num_var, info) 3193 3194 CALL timestop(handle4) 3195 3196 DO iii = 1, num_var 3197 Lambda(iii) = Lambda_Re(iii) + im_unit*Lambda_Im(iii) + vec_b_gw(iii) + vec_b_gw(iii + num_var) 3198 END DO 3199 3200 ! calculate chi2 3201 CALL calc_chi2(chi2, Lambda, vec_Sigma_c_gw, vec_omega_fit_gw_sign, num_poles, & 3202 num_fit_points, n_level_gw) 3203 3204 ! if the fit is already super accurate, exit. otherwise maybe issues when dividing by 0 3205 IF (chi2 < 1.0E-30_dp) EXIT 3206 3207 IF (chi2 < chi2_old) THEN 3208 ScalParam = MAX(ScalParam/Ldown, 1E-12_dp) 3209 DO iii = 1, num_var 3210 Lambda_Re(iii) = Lambda_Re(iii) + REAL(vec_b_gw(iii) + vec_b_gw(iii + num_var)) 3211 Lambda_Im(iii) = Lambda_Im(iii) + AIMAG(vec_b_gw(iii) + vec_b_gw(iii + num_var)) 3212 END DO 3213 IF (chi2_old/chi2 - 1.0_dp < stop_crit) could_exit = .TRUE. 3214 chi2_old = chi2 3215 ELSE 3216 ScalParam = ScalParam*Lup 3217 END IF 3218 IF (ScalParam > 100.0_dp .AND. could_exit) EXIT 3219 3220 IF (ScalParam > 1E+10_dp) ScalParam = 1E-4_dp 3221 3222 END DO 3223 3224 IF (.NOT. do_gw_im_time) THEN 3225 3226 ! change a_0 [Lambda(1)], so that Sigma(i0) = Fit(i0) 3227 ! do not do this for imaginary time since we do not have many fit points and the fit should be perfect 3228 func_val = Lambda(1) 3229 DO iii = 1, num_poles 3230 jjj = iii*2 3231 ! calculate value of the fit function 3232 func_val = func_val + Lambda(jjj)/(-Lambda(jjj + 1)) 3233 END DO 3234 3235 Lambda_Re(1) = Lambda_Re(1) - REAL(func_val) + REAL(vec_Sigma_c_gw(n_level_gw, num_fit_points)) 3236 Lambda_Im(1) = Lambda_Im(1) - AIMAG(func_val) + AIMAG(vec_Sigma_c_gw(n_level_gw, num_fit_points)) 3237 3238 END IF 3239 3240 Lambda_without_offset(:) = Lambda(:) 3241 3242 DO iii = 1, num_var 3243 Lambda(iii) = CMPLX(Lambda_Re(iii), Lambda_Im(iii), KIND=dp) 3244 END DO 3245 3246 IF (do_gw_im_time) THEN 3247 ! for cubic-scaling GW, we have one Green's function for occ and virt states with the Fermi level 3248 ! in the middle of homo and lumo 3249 e_fermi = 0.5_dp*(Eigenval(homo) + Eigenval(homo + 1)) 3250 ELSE 3251 ! in case of O(N^4) GW, we have the Fermi level differently for occ and virt states, see 3252 ! Fig. 1 in JCTC 12, 3623-3635 (2016) 3253 IF (n_level_gw <= gw_corr_lev_occ) THEN 3254 e_fermi = Eigenval(homo) + fermi_level_offset 3255 ELSE 3256 e_fermi = Eigenval(homo + 1) - fermi_level_offset 3257 END IF 3258 END IF 3259 3260 ! either Z-shot or Newton/bisection crossing search for evaluating Sigma_c 3261 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot .OR. & 3262 crossing_search == ri_rpa_g0w0_crossing_newton) THEN 3263 3264 ! calculate Sigma_c_fit(e_n) and Z 3265 func_val = Lambda(1) 3266 z_value(n_level_gw) = 1.0_dp 3267 DO iii = 1, num_poles 3268 jjj = iii*2 3269 z_value(n_level_gw) = z_value(n_level_gw) + REAL(Lambda(jjj)/ & 3270 (Eigenval(n_level_gw_ref) - e_fermi - Lambda(jjj + 1))**2) 3271 func_val = func_val + Lambda(jjj)/(Eigenval(n_level_gw_ref) - e_fermi - Lambda(jjj + 1)) 3272 END DO 3273 ! m is the slope of the correl self-energy 3274 m_value(n_level_gw) = 1.0_dp - z_value(n_level_gw) 3275 z_value(n_level_gw) = 1.0_dp/z_value(n_level_gw) 3276 gw_energ = REAL(func_val) 3277 vec_gw_energ(n_level_gw) = gw_energ 3278 3279 ! in case one wants to do Newton-Raphson on top of the Z-shot 3280 IF (crossing_search == ri_rpa_g0w0_crossing_newton) THEN 3281 3282 level_energ_GW = (Eigenval_scf(n_level_gw_ref) - & 3283 m_value(n_level_gw)*Eigenval(n_level_gw_ref) + & 3284 vec_gw_energ(n_level_gw) + & 3285 vec_Sigma_x_minus_vxc_gw(n_level_gw_ref))* & 3286 z_value(n_level_gw) 3287 3288 ! Newton-Raphson iteration 3289 DO kkk = 1, 1000 3290 3291 ! calculate the value of the fit function for level_energ_GW 3292 func_val = Lambda(1) 3293 z_value(n_level_gw) = 1.0_dp 3294 DO iii = 1, num_poles 3295 jjj = iii*2 3296 func_val = func_val + Lambda(jjj)/(level_energ_GW - e_fermi - Lambda(jjj + 1)) 3297 END DO 3298 3299 ! calculate the derivative of the fit function for level_energ_GW 3300 deriv_val_real = -1.0_dp 3301 DO iii = 1, num_poles 3302 jjj = iii*2 3303 deriv_val_real = deriv_val_real + REAL(Lambda(jjj))/((ABS(level_energ_GW - e_fermi - Lambda(jjj + 1)))**2) & 3304 - (REAL(Lambda(jjj))*(level_energ_GW - e_fermi) - REAL(Lambda(jjj)*CONJG(Lambda(jjj + 1))))* & 3305 2.0_dp*(level_energ_GW - e_fermi - REAL(Lambda(jjj + 1)))/ & 3306 ((ABS(level_energ_GW - e_fermi - Lambda(jjj + 1)))**2) 3307 3308 END DO 3309 3310 delta = (Eigenval_scf(n_level_gw_ref) + vec_Sigma_x_minus_vxc_gw(n_level_gw_ref) + REAL(func_val) - level_energ_GW)/ & 3311 deriv_val_real 3312 3313 level_energ_GW = level_energ_GW - delta 3314 3315 IF (ABS(delta) < 1.0E-08) EXIT 3316 3317 END DO 3318 3319 ! update the GW-energy by Newton-Raphson and set the Z-value to 1 3320 3321 vec_gw_energ(n_level_gw) = REAL(func_val) 3322 z_value(n_level_gw) = 1.0_dp 3323 m_value(n_level_gw) = 0.0_dp 3324 3325 END IF ! Newton-Raphson on top of Z-shot 3326 3327 ELSE 3328 CPABORT("Only NONE, ZSHOT and NEWTON implemented for 2-pole model") 3329 END IF ! decision crossing search none, Z-shot 3330 3331 ! -------------------------------------------- 3332 ! | calculate statistical error due to fitting | 3333 ! -------------------------------------------- 3334 3335 ! estimate the statistical error of the calculated Sigma_c(i*omega) 3336 ! by sqrt(chi2/n), where n is the number of fit points 3337 3338 CALL calc_chi2(chi2, Lambda_without_offset, vec_Sigma_c_gw, vec_omega_fit_gw_sign, num_poles, & 3339 num_fit_points, n_level_gw) 3340 3341 ! Estimate the statistical error of every fit point 3342 stat_error = SQRT(chi2/num_fit_points) 3343 3344 ! allocate N array containing the second derivatives of chi^2 3345 ALLOCATE (vec_N_gw(num_var*2)) 3346 vec_N_gw = 0.0_dp 3347 3348 ALLOCATE (mat_N_gw(num_var*2, num_var*2)) 3349 mat_N_gw = 0.0_dp 3350 3351 DO iii = 1, num_var*2 3352 CALL calc_mat_N(vec_N_gw(iii), Lambda_without_offset, vec_Sigma_c_gw, vec_omega_fit_gw_sign, & 3353 iii, iii, num_poles, num_fit_points, n_level_gw, 0.001_dp) 3354 END DO 3355 3356 DO iii = 1, num_var*2 3357 DO jjj = 1, num_var*2 3358 CALL calc_mat_N(mat_N_gw(iii, jjj), Lambda_without_offset, vec_Sigma_c_gw, vec_omega_fit_gw_sign, & 3359 iii, jjj, num_poles, num_fit_points, n_level_gw, 0.001_dp) 3360 END DO 3361 END DO 3362 3363 CALL DGETRF(2*num_var, 2*num_var, mat_N_gw, 2*num_var, ipiv, info) 3364 3365 ! vec_b_gw is only working array 3366 CALL DGETRI(2*num_var, mat_N_gw, 2*num_var, ipiv, vec_b_gw, 2*num_var, info) 3367 3368 ALLOCATE (stat_errors(2*num_var)) 3369 stat_errors = 0.0_dp 3370 3371 DO iii = 1, 2*num_var 3372 stat_errors(iii) = SQRT(ABS(mat_N_gw(iii, iii)))*stat_error 3373 END DO 3374 3375 DEALLOCATE (mat_N_gw) 3376 DEALLOCATE (vec_N_gw) 3377 DEALLOCATE (mat_A_gw) 3378 DEALLOCATE (mat_B_gw) 3379 DEALLOCATE (stat_errors) 3380 DEALLOCATE (dLambda) 3381 DEALLOCATE (dLambda_2) 3382 DEALLOCATE (vec_b_gw) 3383 DEALLOCATE (vec_b_gw_copy) 3384 DEALLOCATE (ipiv) 3385 DEALLOCATE (vec_omega_fit_gw_sign) 3386 DEALLOCATE (Lambda) 3387 DEALLOCATE (Lambda_without_offset) 3388 DEALLOCATE (Lambda_Re) 3389 DEALLOCATE (Lambda_Im) 3390 3391 END SUBROUTINE fit_and_continuation_2pole 3392 3393! ************************************************************************************************** 3394!> \brief perform analytic continuation with pade approximation 3395!> \param vec_gw_energ real Sigma_c 3396!> \param vec_omega_fit_gw frequency points for Sigma_c(iomega) 3397!> \param z_value 1/(1-dev) 3398!> \param m_value derivative of real Sigma_c 3399!> \param vec_Sigma_c_gw complex Sigma_c(iomega) 3400!> \param vec_Sigma_x_minus_vxc_gw ... 3401!> \param Eigenval quasiparticle energy during ev self-consistent GW 3402!> \param Eigenval_scf KS/HF eigenvalue 3403!> \param n_level_gw ... 3404!> \param gw_corr_lev_occ ... 3405!> \param nparam_pade number of pade parameters 3406!> \param num_fit_points number of fit points for Sigma_c(iomega) 3407!> \param crossing_search type ofr cross search to find quasiparticle energies 3408!> \param homo ... 3409!> \param fermi_level_offset ... 3410!> \param do_gw_im_time ... 3411!> \param print_self_energy ... 3412!> \param count_ev_sc_GW ... 3413! ************************************************************************************************** 3414 SUBROUTINE continuation_pade(vec_gw_energ, vec_omega_fit_gw, & 3415 z_value, m_value, vec_Sigma_c_gw, vec_Sigma_x_minus_vxc_gw, & 3416 Eigenval, Eigenval_scf, n_level_gw, gw_corr_lev_occ, nparam_pade, & 3417 num_fit_points, crossing_search, homo, & 3418 fermi_level_offset, do_gw_im_time, print_self_energy, count_ev_sc_GW) 3419 3420 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec_gw_energ 3421 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec_omega_fit_gw 3422 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: z_value, m_value 3423 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vec_Sigma_c_gw 3424 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec_Sigma_x_minus_vxc_gw, Eigenval, & 3425 Eigenval_scf 3426 INTEGER, INTENT(IN) :: n_level_gw, gw_corr_lev_occ, & 3427 nparam_pade, num_fit_points, & 3428 crossing_search, homo 3429 REAL(KIND=dp), INTENT(IN) :: fermi_level_offset 3430 LOGICAL, INTENT(IN) :: do_gw_im_time, print_self_energy 3431 INTEGER, INTENT(IN) :: count_ev_sc_GW 3432 3433 CHARACTER(LEN=*), PARAMETER :: routineN = 'continuation_pade', & 3434 routineP = moduleN//':'//routineN 3435 3436 CHARACTER(len=default_path_length) :: filename 3437 COMPLEX(KIND=dp) :: im_unit, re_unit, sigma_c_pade 3438 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeff_pade, omega_points_pade, & 3439 Sigma_c_gw_reorder 3440 INTEGER :: handle, i_omega, iunit, jquad, & 3441 n_level_gw_ref, num_omega 3442 REAL(KIND=dp) :: e_fermi, energy_val, level_energ_GW, & 3443 omega, sign_occ_virt 3444 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vec_omega_fit_gw_sign, & 3445 vec_omega_fit_gw_sign_reorder 3446 3447 CALL timeset(routineN, handle) 3448 3449 im_unit = (0.0_dp, 1.0_dp) 3450 re_unit = (1.0_dp, 0.0_dp) 3451 3452 ALLOCATE (vec_omega_fit_gw_sign(num_fit_points)) 3453 3454 IF (n_level_gw <= gw_corr_lev_occ) THEN 3455 sign_occ_virt = -1.0_dp 3456 ELSE 3457 sign_occ_virt = 1.0_dp 3458 END IF 3459 3460 DO jquad = 1, num_fit_points 3461 vec_omega_fit_gw_sign(jquad) = ABS(vec_omega_fit_gw(jquad))*sign_occ_virt 3462 END DO 3463 3464 IF (do_gw_im_time) THEN 3465 ! for cubic-scaling GW, we have one Green's function for occ and virt states with the Fermi level 3466 ! in the middle of homo and lumo 3467 e_fermi = 0.5_dp*(Eigenval(homo) + Eigenval(homo + 1)) 3468 ELSE 3469 ! in case of O(N^4) GW, we have the Fermi level differently for occ and virt states, see 3470 ! Fig. 1 in JCTC 12, 3623-3635 (2016) 3471 IF (n_level_gw <= gw_corr_lev_occ) THEN 3472 e_fermi = Eigenval(homo) + fermi_level_offset 3473 ELSE 3474 e_fermi = Eigenval(homo + 1) - fermi_level_offset 3475 END IF 3476 END IF 3477 3478 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 3479 3480 !*** reorder, such that omega=i*0 is first entry 3481 ALLOCATE (Sigma_c_gw_reorder(num_fit_points)) 3482 ALLOCATE (vec_omega_fit_gw_sign_reorder(num_fit_points)) 3483 ! for cubic scaling GW fit points are ordered differently than in N^4 GW 3484 IF (do_gw_im_time) THEN 3485 DO jquad = 1, num_fit_points 3486 Sigma_c_gw_reorder(jquad) = vec_Sigma_c_gw(n_level_gw, jquad) 3487 vec_omega_fit_gw_sign_reorder(jquad) = vec_omega_fit_gw_sign(jquad) 3488 ENDDO 3489 ELSE 3490 DO jquad = 1, num_fit_points 3491 Sigma_c_gw_reorder(jquad) = vec_Sigma_c_gw(n_level_gw, num_fit_points - jquad + 1) 3492 vec_omega_fit_gw_sign_reorder(jquad) = vec_omega_fit_gw_sign(num_fit_points - jquad + 1) 3493 ENDDO 3494 ENDIF 3495 3496 !*** evaluate parameters for pade approximation 3497 ALLOCATE (coeff_pade(nparam_pade)) 3498 ALLOCATE (omega_points_pade(nparam_pade)) 3499 coeff_pade = 0.0_dp 3500 CALL get_pade_parameters(Sigma_c_gw_reorder, vec_omega_fit_gw_sign_reorder, & 3501 num_fit_points, nparam_pade, omega_points_pade, coeff_pade) 3502 3503 !*** calculate start_value for iterative cross-searching methods 3504 IF ((crossing_search == ri_rpa_g0w0_crossing_bisection) .OR. & 3505 (crossing_search == ri_rpa_g0w0_crossing_newton)) THEN 3506 energy_val = Eigenval(n_level_gw_ref) - e_fermi 3507 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, & 3508 coeff_pade, sigma_c_pade) 3509 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, & 3510 coeff_pade, z_value(n_level_gw), m_value(n_level_gw)) 3511 level_energ_GW = (Eigenval_scf(n_level_gw_ref) - & 3512 m_value(n_level_gw)*Eigenval(n_level_gw_ref) + & 3513 REAL(sigma_c_pade) + & 3514 vec_Sigma_x_minus_vxc_gw(n_level_gw_ref))* & 3515 z_value(n_level_gw) 3516 ENDIF 3517 3518 !*** perform crossing search 3519 SELECT CASE (crossing_search) 3520 CASE (ri_rpa_g0w0_crossing_z_shot) 3521 energy_val = Eigenval(n_level_gw_ref) - e_fermi 3522 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, & 3523 coeff_pade, sigma_c_pade) 3524 vec_gw_energ(n_level_gw) = REAL(sigma_c_pade) 3525 3526 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, & 3527 coeff_pade, z_value(n_level_gw), m_value(n_level_gw)) 3528 3529 CASE (ri_rpa_g0w0_crossing_bisection) 3530 CALL get_sigma_c_bisection_pade(vec_gw_energ(n_level_gw), Eigenval_scf(n_level_gw_ref), & 3531 vec_Sigma_x_minus_vxc_gw(n_level_gw_ref), e_fermi, & 3532 nparam_pade, omega_points_pade, coeff_pade, & 3533 n_level_gw_ref, start_val=level_energ_GW) 3534 z_value(n_level_gw) = 1.0_dp 3535 m_value(n_level_gw) = 0.0_dp 3536 3537 CASE (ri_rpa_g0w0_crossing_newton) 3538 CALL get_sigma_c_newton_pade(vec_gw_energ(n_level_gw), Eigenval_scf(n_level_gw_ref), & 3539 vec_Sigma_x_minus_vxc_gw(n_level_gw_ref), e_fermi, & 3540 nparam_pade, omega_points_pade, coeff_pade, & 3541 n_level_gw_ref, start_val=level_energ_GW) 3542 z_value(n_level_gw) = 1.0_dp 3543 m_value(n_level_gw) = 0.0_dp 3544 3545 CASE DEFAULT 3546 CPABORT("Only Z_SHOT, NEWTON, and BISECTION crossing search implemented.") 3547 END SELECT 3548 3549 IF (print_self_energy) THEN 3550 3551 IF (count_ev_sc_GW == 1) THEN 3552 3553 IF (n_level_gw_ref < 10) THEN 3554 WRITE (filename, "(A26,I1)") "G0W0_self_energy_level_000", n_level_gw_ref 3555 ELSE IF (n_level_gw_ref < 100) THEN 3556 WRITE (filename, "(A25,I2)") "G0W0_self_energy_level_00", n_level_gw_ref 3557 ELSE IF (n_level_gw_ref < 1000) THEN 3558 WRITE (filename, "(A24,I3)") "G0W0_self_energy_level_0", n_level_gw_ref 3559 ELSE 3560 WRITE (filename, "(A23,I4)") "G0W0_self_energy_level_", n_level_gw_ref 3561 END IF 3562 3563 ELSE 3564 3565 IF (n_level_gw_ref < 10) THEN 3566 WRITE (filename, "(A11,I1,A22,I1)") "evGW_cycle_", count_ev_sc_GW, & 3567 "_self_energy_level_000", n_level_gw_ref 3568 ELSE IF (n_level_gw_ref < 100) THEN 3569 WRITE (filename, "(A11,I1,A21,I2)") "evGW_cycle_", count_ev_sc_GW, & 3570 "_self_energy_level_00", n_level_gw_ref 3571 ELSE IF (n_level_gw_ref < 1000) THEN 3572 WRITE (filename, "(A11,I1,A20,I3)") "evGW_cycle_", count_ev_sc_GW, & 3573 "_self_energy_level_0", n_level_gw_ref 3574 ELSE 3575 WRITE (filename, "(A11,I1,A19,I4)") "evGW_cycle_", count_ev_sc_GW, & 3576 "_self_energy_level_", n_level_gw_ref 3577 END IF 3578 3579 END IF 3580 3581 CALL open_file(TRIM(filename), unit_number=iunit, file_status="UNKNOWN", file_action="WRITE") 3582 3583 num_omega = 10000 3584 3585 WRITE (iunit, "(2A42)") " omega (eV) Sigma(omega) (eV) ", & 3586 " omega - e_n^DFT - Sigma_n^x - v_n^xc (eV)" 3587 3588 DO i_omega = 0, num_omega 3589 3590 omega = -50.0_dp/evolt + REAL(i_omega, KIND=dp)/REAL(num_omega, KIND=dp)*100.0_dp/evolt 3591 3592 CALL evaluate_pade_function(omega - e_fermi, nparam_pade, omega_points_pade, & 3593 coeff_pade, sigma_c_pade) 3594 3595 WRITE (iunit, "(F12.2,2F17.5)") omega*evolt, REAL(sigma_c_pade)*evolt, & 3596 (omega - Eigenval_scf(n_level_gw_ref) - vec_Sigma_x_minus_vxc_gw(n_level_gw_ref))*evolt 3597 3598 END DO 3599 3600 CALL close_file(iunit) 3601 3602 END IF 3603 3604 DEALLOCATE (vec_omega_fit_gw_sign) 3605 DEALLOCATE (Sigma_c_gw_reorder) 3606 DEALLOCATE (vec_omega_fit_gw_sign_reorder) 3607 DEALLOCATE (coeff_pade, omega_points_pade) 3608 3609 CALL timestop(handle) 3610 3611 END SUBROUTINE continuation_pade 3612 3613! ************************************************************************************************** 3614!> \brief calculate pade parameter recursively as in Eq. (A2) in J. Low Temp. Phys., Vol. 29, 3615!> 1977, pp. 179 3616!> \param y f(x), here: Sigma_c(iomega) 3617!> \param x the frequency points omega 3618!> \param num_fit_points ... 3619!> \param nparam number of pade parameters 3620!> \param xpoints set of points used in pade approximation, selection of x 3621!> \param coeff pade coefficients 3622! ************************************************************************************************** 3623 SUBROUTINE get_pade_parameters(y, x, num_fit_points, nparam, xpoints, coeff) 3624 3625 COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: y 3626 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: x 3627 INTEGER, INTENT(IN) :: num_fit_points, nparam 3628 COMPLEX(KIND=dp), DIMENSION(:), INTENT(INOUT) :: xpoints, coeff 3629 3630 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_pade_parameters', & 3631 routineP = moduleN//':'//routineN 3632 3633 COMPLEX(KIND=dp) :: im_unit 3634 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ypoints 3635 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: g_mat 3636 INTEGER :: handle, idat, iparam, nstep 3637 3638 CALL timeset(routineN, handle) 3639 3640 im_unit = (0.0_dp, 1.0_dp) 3641 3642 nstep = INT(num_fit_points/(nparam - 1)) 3643 CPASSERT(LBOUND(x, 1) == 1) 3644 CPASSERT(LBOUND(y, 1) == 1) 3645 3646 ALLOCATE (ypoints(nparam)) 3647 !omega=i0 is in element x(1) 3648 idat = 1 3649 DO iparam = 1, nparam - 1 3650 xpoints(iparam) = im_unit*x(idat) 3651 ypoints(iparam) = y(idat) 3652 idat = idat + nstep 3653 ENDDO 3654 xpoints(nparam) = im_unit*x(num_fit_points) 3655 ypoints(nparam) = y(num_fit_points) 3656 3657 !*** generate parameters recursively 3658 3659 ALLOCATE (g_mat(nparam, nparam)) 3660 g_mat(:, 1) = ypoints(:) 3661 DO iparam = 2, nparam 3662 DO idat = iparam, nparam 3663 g_mat(idat, iparam) = (g_mat(iparam - 1, iparam - 1) - g_mat(idat, iparam - 1))/ & 3664 ((xpoints(idat) - xpoints(iparam - 1))*g_mat(idat, iparam - 1)) 3665 ENDDO 3666 ENDDO 3667 3668 DO iparam = 1, nparam 3669 coeff(iparam) = g_mat(iparam, iparam) 3670 ENDDO 3671 3672 DEALLOCATE (ypoints) 3673 DEALLOCATE (g_mat) 3674 3675 CALL timestop(handle) 3676 3677 END SUBROUTINE get_pade_parameters 3678 3679! ************************************************************************************************** 3680!> \brief evaluate pade function for a real value x_val 3681!> \param x_val real value 3682!> \param nparam number of pade parameters 3683!> \param xpoints selection of points of the original complex function, i.e. here of Sigma_c(iomega) 3684!> \param coeff pade coefficients 3685!> \param func_val function value 3686! ************************************************************************************************** 3687 SUBROUTINE evaluate_pade_function(x_val, nparam, xpoints, coeff, func_val) 3688 3689 REAL(KIND=dp), INTENT(IN) :: x_val 3690 INTEGER, INTENT(IN) :: nparam 3691 COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: xpoints, coeff 3692 COMPLEX(KIND=dp), INTENT(OUT) :: func_val 3693 3694 CHARACTER(LEN=*), PARAMETER :: routineN = 'evaluate_pade_function', & 3695 routineP = moduleN//':'//routineN 3696 3697 COMPLEX(KIND=dp) :: im_unit, re_unit 3698 INTEGER :: handle, iparam 3699 3700 CALL timeset(routineN, handle) 3701 3702 im_unit = (0.0_dp, 1.0_dp) 3703 re_unit = (1.0_dp, 0.0_dp) 3704 3705 func_val = re_unit 3706 DO iparam = nparam, 2, -1 3707 func_val = re_unit + coeff(iparam)*(re_unit*x_val - xpoints(iparam - 1))/func_val 3708 ENDDO 3709 3710 func_val = coeff(1)/func_val 3711 3712 CALL timestop(handle) 3713 3714 END SUBROUTINE evaluate_pade_function 3715 3716! ************************************************************************************************** 3717!> \brief get the z-value and the m-value (derivative) of the pade function 3718!> \param x_val real value 3719!> \param nparam number of pade parameters 3720!> \param xpoints selection of points of the original complex function, i.e. here of Sigma_c(iomega) 3721!> \param coeff pade coefficients 3722!> \param z_value 1/(1-dev) 3723!> \param m_value derivative 3724! ************************************************************************************************** 3725 SUBROUTINE get_z_and_m_value_pade(x_val, nparam, xpoints, coeff, z_value, m_value) 3726 3727 REAL(KIND=dp), INTENT(IN) :: x_val 3728 INTEGER, INTENT(IN) :: nparam 3729 COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: xpoints, coeff 3730 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: z_value, m_value 3731 3732 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_z_and_m_value_pade', & 3733 routineP = moduleN//':'//routineN 3734 3735 COMPLEX(KIND=dp) :: denominator, dev_denominator, & 3736 dev_numerator, dev_val, func_val, & 3737 im_unit, numerator, re_unit 3738 INTEGER :: iparam 3739 3740 im_unit = (0.0_dp, 1.0_dp) 3741 re_unit = (1.0_dp, 0.0_dp) 3742 3743 func_val = re_unit 3744 dev_val = (0.0_dp, 0.0_dp) 3745 DO iparam = nparam, 2, -1 3746 numerator = coeff(iparam)*(re_unit*x_val - xpoints(iparam - 1)) 3747 dev_numerator = coeff(iparam)*re_unit 3748 denominator = func_val 3749 dev_denominator = dev_val 3750 dev_val = dev_numerator/denominator - (numerator*dev_denominator)/(denominator**2) 3751 func_val = re_unit + coeff(iparam)*(re_unit*x_val - xpoints(iparam - 1))/func_val 3752 ENDDO 3753 3754 dev_val = -1.0_dp*coeff(1)/(func_val**2)*dev_val 3755 func_val = coeff(1)/func_val 3756 3757 IF (PRESENT(z_value)) THEN 3758 z_value = 1.0_dp - REAL(dev_val) 3759 z_value = 1.0_dp/z_value 3760 ENDIF 3761 IF (PRESENT(m_value)) m_value = REAL(dev_val) 3762 3763 END SUBROUTINE get_z_and_m_value_pade 3764 3765! ************************************************************************************************** 3766!> \brief crossing search using the bisection method to find the quasiparticle energy 3767!> \param gw_energ real Sigma_c 3768!> \param Eigenval_scf Eigenvalue from the SCF 3769!> \param Sigma_x_minus_vxc_gw ... 3770!> \param e_fermi fermi level 3771!> \param nparam_pade number of pade parameters 3772!> \param omega_points_pade selection of frequency points of Sigma_c(iomega) 3773!> \param coeff_pade pade coefficients 3774!> \param n_level_gw_ref ... 3775!> \param start_val start value for the quasiparticle iteration 3776! ************************************************************************************************** 3777 SUBROUTINE get_sigma_c_bisection_pade(gw_energ, Eigenval_scf, Sigma_x_minus_vxc_gw, e_fermi, & 3778 nparam_pade, omega_points_pade, coeff_pade, n_level_gw_ref, start_val) 3779 3780 REAL(KIND=dp), INTENT(OUT) :: gw_energ 3781 REAL(KIND=dp), INTENT(IN) :: Eigenval_scf, Sigma_x_minus_vxc_gw, & 3782 e_fermi 3783 INTEGER, INTENT(IN) :: nparam_pade 3784 COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: omega_points_pade, coeff_pade 3785 INTEGER, INTENT(IN) :: n_level_gw_ref 3786 REAL(KIND=dp), INTENT(IN), OPTIONAL :: start_val 3787 3788 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_sigma_c_bisection_pade', & 3789 routineP = moduleN//':'//routineN 3790 3791 CHARACTER(LEN=512) :: error_msg 3792 CHARACTER(LEN=64) :: n_level_gw_ref_char 3793 COMPLEX(KIND=dp) :: sigma_c 3794 INTEGER :: handle, icount 3795 REAL(KIND=dp) :: delta, energy_val, my_start_val, & 3796 qp_energy, qp_energy_old, threshold 3797 3798 CALL timeset(routineN, handle) 3799 3800 threshold = 1.0E-7_dp 3801 3802 IF (PRESENT(start_val)) THEN 3803 my_start_val = start_val 3804 ELSE 3805 my_start_val = Eigenval_scf 3806 ENDIF 3807 3808 qp_energy = my_start_val 3809 qp_energy_old = my_start_val 3810 delta = 1.0E-3_dp 3811 3812 icount = 0 3813 DO WHILE (ABS(delta) > threshold) 3814 icount = icount + 1 3815 qp_energy = qp_energy_old + 0.5_dp*delta 3816 qp_energy_old = qp_energy 3817 energy_val = qp_energy - e_fermi 3818 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, & 3819 coeff_pade, sigma_c) 3820 qp_energy = Eigenval_scf + REAL(sigma_c) + Sigma_x_minus_vxc_gw 3821 delta = qp_energy - qp_energy_old 3822 IF (icount > 500) THEN 3823 WRITE (n_level_gw_ref_char, '(I10)') n_level_gw_ref 3824 WRITE (error_msg, '(A,A,A)') " Self-consistent quasi-particle solution of "// & 3825 "MO ", TRIM(n_level_gw_ref_char), " has not been found." 3826 CPWARN(error_msg) 3827 EXIT 3828 ENDIF 3829 ENDDO 3830 3831 gw_energ = REAL(sigma_c) 3832 3833 CALL timestop(handle) 3834 3835 END SUBROUTINE get_sigma_c_bisection_pade 3836 3837! ************************************************************************************************** 3838!> \brief crossing search using the Newton method to find the quasiparticle energy 3839!> \param gw_energ real Sigma_c 3840!> \param Eigenval_scf Eigenvalue from the SCF 3841!> \param Sigma_x_minus_vxc_gw ... 3842!> \param e_fermi fermi level 3843!> \param nparam_pade number of pade parameters 3844!> \param omega_points_pade selection of frequency points of Sigma_c(iomega) 3845!> \param coeff_pade pade coefficients 3846!> \param n_level_gw_ref ... 3847!> \param start_val start value for the quasiparticle iteration 3848! ************************************************************************************************** 3849 SUBROUTINE get_sigma_c_newton_pade(gw_energ, Eigenval_scf, Sigma_x_minus_vxc_gw, e_fermi, & 3850 nparam_pade, omega_points_pade, coeff_pade, n_level_gw_ref, start_val) 3851 3852 REAL(KIND=dp), INTENT(OUT) :: gw_energ 3853 REAL(KIND=dp), INTENT(IN) :: Eigenval_scf, Sigma_x_minus_vxc_gw, & 3854 e_fermi 3855 INTEGER, INTENT(IN) :: nparam_pade 3856 COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: omega_points_pade, coeff_pade 3857 INTEGER, INTENT(IN) :: n_level_gw_ref 3858 REAL(KIND=dp), INTENT(IN), OPTIONAL :: start_val 3859 3860 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_sigma_c_newton_pade', & 3861 routineP = moduleN//':'//routineN 3862 3863 CHARACTER(LEN=512) :: error_msg 3864 CHARACTER(LEN=64) :: n_level_gw_ref_char 3865 COMPLEX(KIND=dp) :: sigma_c 3866 INTEGER :: handle, icount 3867 REAL(KIND=dp) :: delta, energy_val, m_value, & 3868 my_start_val, qp_energy, & 3869 qp_energy_old, threshold 3870 3871 CALL timeset(routineN, handle) 3872 3873 threshold = 1.0E-7_dp 3874 3875 IF (PRESENT(start_val)) THEN 3876 my_start_val = start_val 3877 ELSE 3878 my_start_val = Eigenval_scf 3879 ENDIF 3880 3881 qp_energy = my_start_val 3882 qp_energy_old = my_start_val 3883 delta = 1.0E-3_dp 3884 3885 icount = 0 3886 DO WHILE (ABS(delta) > threshold) 3887 icount = icount + 1 3888 energy_val = qp_energy - e_fermi 3889 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, & 3890 coeff_pade, sigma_c) 3891 !get m_value --> derivative of function 3892 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, & 3893 coeff_pade, m_value=m_value) 3894 qp_energy_old = qp_energy 3895 qp_energy = qp_energy - (Eigenval_scf + Sigma_x_minus_vxc_gw + REAL(sigma_c) - qp_energy)/ & 3896 (m_value - 1.0_dp) 3897 delta = qp_energy - qp_energy_old 3898 IF (icount > 500) THEN 3899 WRITE (n_level_gw_ref_char, '(I10)') n_level_gw_ref 3900 WRITE (error_msg, '(A,A,A)') " Self-consistent quasi-particle solution of "// & 3901 "MO ", TRIM(n_level_gw_ref_char), " has not been found." 3902 CPWARN(error_msg) 3903 EXIT 3904 ENDIF 3905 ENDDO 3906 3907 gw_energ = REAL(sigma_c) 3908 3909 CALL timestop(handle) 3910 3911 END SUBROUTINE get_sigma_c_newton_pade 3912 3913! ************************************************************************************************** 3914!> \brief Prints the GW stuff to the output and optinally to an external file. 3915!> Also updates the eigenvalues for eigenvalue-self-consistent GW 3916!> \param vec_gw_energ ... 3917!> \param z_value ... 3918!> \param m_value ... 3919!> \param vec_Sigma_x_minus_vxc_gw ... 3920!> \param Eigenval ... 3921!> \param Eigenval_last ... 3922!> \param Eigenval_scf ... 3923!> \param gw_corr_lev_occ ... 3924!> \param gw_corr_lev_tot ... 3925!> \param crossing_search ... 3926!> \param homo ... 3927!> \param unit_nr ... 3928!> \param count_ev_sc_GW ... 3929!> \param count_sc_GW0 ... 3930!> \param ikp ... 3931!> \param nkp_self_energy ... 3932!> \param kpoints ... 3933!> \param do_alpha ... 3934!> \param do_beta ... 3935! ************************************************************************************************** 3936 SUBROUTINE print_and_update_for_ev_sc(vec_gw_energ, & 3937 z_value, m_value, vec_Sigma_x_minus_vxc_gw, Eigenval, & 3938 Eigenval_last, Eigenval_scf, gw_corr_lev_occ, gw_corr_lev_tot, & 3939 crossing_search, homo, unit_nr, count_ev_sc_GW, count_sc_GW0, & 3940 ikp, nkp_self_energy, kpoints, do_alpha, do_beta) 3941 3942 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 3943 INTENT(IN) :: vec_gw_energ, z_value, m_value 3944 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec_Sigma_x_minus_vxc_gw, Eigenval 3945 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 3946 INTENT(INOUT) :: Eigenval_last, Eigenval_scf 3947 INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_tot, crossing_search, homo, unit_nr, & 3948 count_ev_sc_GW, count_sc_GW0, ikp, nkp_self_energy 3949 TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints 3950 LOGICAL, INTENT(IN), OPTIONAL :: do_alpha, do_beta 3951 3952 CHARACTER(LEN=*), PARAMETER :: routineN = 'print_and_update_for_ev_sc', & 3953 routineP = moduleN//':'//routineN 3954 3955 CHARACTER(4) :: occ_virt 3956 INTEGER :: handle, n_level_gw, n_level_gw_ref 3957 LOGICAL :: do_closed_shell, do_kpoints, & 3958 is_energy_okay, my_do_alpha, my_do_beta 3959 REAL(KIND=dp) :: new_energy 3960 3961 CALL timeset(routineN, handle) 3962 3963 IF (PRESENT(do_alpha)) THEN 3964 my_do_alpha = do_alpha 3965 ELSE 3966 my_do_alpha = .FALSE. 3967 END IF 3968 3969 IF (PRESENT(do_beta)) THEN 3970 my_do_beta = do_beta 3971 ELSE 3972 my_do_beta = .FALSE. 3973 END IF 3974 3975 do_closed_shell = .NOT. (my_do_alpha .OR. my_do_beta) 3976 do_kpoints = (nkp_self_energy > 1) 3977 3978 Eigenval_last(:) = Eigenval(:) 3979 3980 IF (unit_nr > 0) THEN 3981 3982 WRITE (unit_nr, *) ' ' 3983 3984 IF (count_ev_sc_GW == 1 .AND. count_sc_GW0 == 1) THEN 3985 3986 IF (do_closed_shell) THEN 3987 WRITE (unit_nr, *) ' ' 3988 WRITE (unit_nr, '(T3,A)') '******************************************************************************' 3989 WRITE (unit_nr, '(T3,A)') '** **' 3990 WRITE (unit_nr, '(T3,A)') '** GW QUASIPARTICLE ENERGIES **' 3991 WRITE (unit_nr, '(T3,A)') '** **' 3992 WRITE (unit_nr, '(T3,A)') '******************************************************************************' 3993 3994 ELSE IF (my_do_alpha) THEN 3995 WRITE (unit_nr, *) ' ' 3996 WRITE (unit_nr, '(T3,A)') '---------------------------------------' 3997 WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies of alpha spins' 3998 WRITE (unit_nr, '(T3,A)') '----------------------------------------' 3999 ELSE IF (my_do_beta) THEN 4000 WRITE (unit_nr, *) ' ' 4001 WRITE (unit_nr, '(T3,A)') '---------------------------------------' 4002 WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies of beta spins' 4003 WRITE (unit_nr, '(T3,A)') '---------------------------------------' 4004 END IF 4005 4006 WRITE (unit_nr, '(T3,A)') ' ' 4007 WRITE (unit_nr, '(T3,A)') ' ' 4008 WRITE (unit_nr, '(T3,A)') 'The GW quasiparticle energies are calculated according to: ' 4009 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot) THEN 4010 WRITE (unit_nr, '(T3,A)') 'E_GW = E_SCF + Z * ( Sigc(E_SCF) + Sigx - vxc )' 4011 ELSE 4012 WRITE (unit_nr, '(T3,A)') ' ' 4013 WRITE (unit_nr, '(T3,A)') ' E_GW = E_SCF + Sigc(E_GW) + Sigx - vxc ' 4014 WRITE (unit_nr, '(T3,A)') ' ' 4015 WRITE (unit_nr, '(T3,A)') 'Upper equation is solved self-consistently for E_GW, see Eq. (12) in J. Phys.' 4016 WRITE (unit_nr, '(T3,A)') 'Chem. Lett. 9, 306 (2018), doi: 10.1021/acs.jpclett.7b02740' 4017 END IF 4018 4019 WRITE (unit_nr, *) ' ' 4020 WRITE (unit_nr, *) ' ' 4021 WRITE (unit_nr, '(T3,A)') '------------' 4022 WRITE (unit_nr, '(T3,A)') 'G0W0 results' 4023 WRITE (unit_nr, '(T3,A)') '------------' 4024 4025 END IF 4026 4027 IF (count_ev_sc_GW > 1) THEN 4028 WRITE (unit_nr, *) ' ' 4029 WRITE (unit_nr, '(T3,A)') '---------------------------------------' 4030 WRITE (unit_nr, '(T3,A,I4)') 'Eigenvalue-selfconsistency cycle: ', count_ev_sc_GW 4031 WRITE (unit_nr, '(T3,A)') '---------------------------------------' 4032 END IF 4033 4034 IF (count_sc_GW0 > 1) THEN 4035 WRITE (unit_nr, '(T3,A)') '----------------------------------' 4036 WRITE (unit_nr, '(T3,A,I4)') 'scGW0 selfconsistency cycle: ', count_sc_GW0 4037 WRITE (unit_nr, '(T3,A)') '----------------------------------' 4038 END IF 4039 4040 IF (do_kpoints) THEN 4041 WRITE (unit_nr, *) ' ' 4042 WRITE (unit_nr, '(T3,A7,I3,A3,I3,A8,3F7.3,A12,3F7.3)') 'Kpoint ', ikp, ' /', nkp_self_energy, & 4043 ' xkp =', kpoints%xkp(1, ikp), kpoints%xkp(2, ikp), kpoints%xkp(3, ikp), & 4044 ' and xkp =', -kpoints%xkp(1, ikp), -kpoints%xkp(2, ikp), -kpoints%xkp(3, ikp) 4045 WRITE (unit_nr, '(T3,A72)') '(Relative Brillouin zone size: [-0.5, 0.5] x [-0.5, 0.5] x [-0.5, 0.5])' 4046 END IF 4047 4048 END IF 4049 4050 DO n_level_gw = 1, gw_corr_lev_tot 4051 4052 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 4053 4054 new_energy = (Eigenval_scf(n_level_gw_ref) - & 4055 m_value(n_level_gw)*Eigenval(n_level_gw_ref) + & 4056 vec_gw_energ(n_level_gw) + & 4057 vec_Sigma_x_minus_vxc_gw(n_level_gw_ref))* & 4058 z_value(n_level_gw) 4059 4060 is_energy_okay = .TRUE. 4061 4062 IF (n_level_gw_ref > homo .AND. new_energy < Eigenval(homo)) THEN 4063 is_energy_okay = .FALSE. 4064 END IF 4065 4066 IF (is_energy_okay) THEN 4067 Eigenval(n_level_gw_ref) = new_energy 4068 END IF 4069 4070 END DO 4071 4072 IF (unit_nr > 0) THEN 4073 WRITE (unit_nr, '(T3,A)') ' ' 4074 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot) THEN 4075 WRITE (unit_nr, '(T13,2A)') 'MO E_SCF (eV) Sigc (eV) Sigx-vxc (eV) Z E_GW (eV)' 4076 ELSE 4077 WRITE (unit_nr, '(T3,2A)') 'Molecular orbital E_SCF (eV) Sigc (eV) Sigx-vxc (eV) E_GW (eV)' 4078 END IF 4079 END IF 4080 4081 DO n_level_gw = 1, gw_corr_lev_tot 4082 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 4083 IF (n_level_gw <= gw_corr_lev_occ) THEN 4084 occ_virt = 'occ' 4085 ELSE 4086 occ_virt = 'vir' 4087 END IF 4088 4089 IF (unit_nr > 0) THEN 4090 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot) THEN 4091 WRITE (unit_nr, '(T3,I4,3A,5F13.3)') & 4092 n_level_gw_ref, ' ( ', occ_virt, ') ', & 4093 Eigenval_last(n_level_gw_ref)*evolt, & 4094 vec_gw_energ(n_level_gw)*evolt, & 4095 vec_Sigma_x_minus_vxc_gw(n_level_gw_ref)*evolt, & 4096 z_value(n_level_gw), & 4097 Eigenval(n_level_gw_ref)*evolt 4098 ELSE 4099 WRITE (unit_nr, '(T3,I4,3A,4F16.3)') & 4100 n_level_gw_ref, ' ( ', occ_virt, ') ', & 4101 Eigenval_last(n_level_gw_ref)*evolt, & 4102 vec_gw_energ(n_level_gw)*evolt, & 4103 vec_Sigma_x_minus_vxc_gw(n_level_gw_ref)*evolt, & 4104 Eigenval(n_level_gw_ref)*evolt 4105 END IF 4106 END IF 4107 END DO 4108 4109 IF (unit_nr > 0) THEN 4110 IF (do_closed_shell) THEN 4111 WRITE (unit_nr, '(T3,A)') ' ' 4112 WRITE (unit_nr, '(T3,A,F57.2)') 'GW HOMO-LUMO gap (eV)', (Eigenval(homo + 1) - Eigenval(homo))*evolt 4113 ELSE IF (my_do_alpha) THEN 4114 WRITE (unit_nr, '(T3,A)') ' ' 4115 WRITE (unit_nr, '(T3,A,F51.2)') 'Alpha GW HOMO-LUMO gap (eV)', (Eigenval(homo + 1) - Eigenval(homo))*evolt 4116 ELSE IF (my_do_beta) THEN 4117 WRITE (unit_nr, '(T3,A)') ' ' 4118 WRITE (unit_nr, '(T3,A,F52.2)') 'Beta GW HOMO-LUMO gap (eV)', (Eigenval(homo + 1) - Eigenval(homo))*evolt 4119 END IF 4120 END IF 4121 4122 IF (unit_nr > 0) THEN 4123 WRITE (unit_nr, *) ' ' 4124 END IF 4125 4126 CALL timestop(handle) 4127 4128 END SUBROUTINE print_and_update_for_ev_sc 4129 4130! ************************************************************************************************** 4131!> \brief ... 4132!> \param Eigenval ... 4133!> \param Eigenval_last ... 4134!> \param gw_corr_lev_occ ... 4135!> \param gw_corr_lev_virt ... 4136!> \param homo ... 4137!> \param nmo ... 4138! ************************************************************************************************** 4139 SUBROUTINE shift_unshifted_levels(Eigenval, Eigenval_last, gw_corr_lev_occ, gw_corr_lev_virt, & 4140 homo, nmo) 4141 4142 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: Eigenval 4143 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 4144 INTENT(INOUT) :: Eigenval_last 4145 INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo, & 4146 nmo 4147 4148 CHARACTER(LEN=*), PARAMETER :: routineN = 'shift_unshifted_levels', & 4149 routineP = moduleN//':'//routineN 4150 4151 INTEGER :: handle, n_level_gw, n_level_gw_ref 4152 REAL(KIND=dp) :: eigen_diff 4153 4154 CALL timeset(routineN, handle) 4155 4156 ! for eigenvalue self-consistent GW, all eigenvalues have to be corrected 4157 ! 1) the occupied; check if there are occupied MOs not being corrected by GW 4158 IF (gw_corr_lev_occ < homo .AND. gw_corr_lev_occ > 0) THEN 4159 4160 ! calculate average GW correction for occupied orbitals 4161 eigen_diff = 0.0_dp 4162 4163 DO n_level_gw = 1, gw_corr_lev_occ 4164 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 4165 eigen_diff = eigen_diff + Eigenval(n_level_gw_ref) - Eigenval_last(n_level_gw_ref) 4166 END DO 4167 eigen_diff = eigen_diff/gw_corr_lev_occ 4168 4169 ! correct the eigenvalues of the occupied orbitals which have not been corrected by GW 4170 DO n_level_gw = 1, homo - gw_corr_lev_occ 4171 Eigenval(n_level_gw) = Eigenval(n_level_gw) + eigen_diff 4172 END DO 4173 4174 END IF 4175 4176 ! 2) the virtual: check if there are virtual orbitals not being corrected by GW 4177 IF (gw_corr_lev_virt < nmo - homo .AND. gw_corr_lev_virt > 0) THEN 4178 4179 ! calculate average GW correction for virtual orbitals 4180 eigen_diff = 0.0_dp 4181 DO n_level_gw = 1, gw_corr_lev_virt 4182 n_level_gw_ref = n_level_gw + homo 4183 eigen_diff = eigen_diff + Eigenval(n_level_gw_ref) - Eigenval_last(n_level_gw_ref) 4184 END DO 4185 eigen_diff = eigen_diff/gw_corr_lev_virt 4186 4187 ! correct the eigenvalues of the virtual orbitals which have not been corrected by GW 4188 DO n_level_gw = homo + gw_corr_lev_virt + 1, nmo 4189 Eigenval(n_level_gw) = Eigenval(n_level_gw) + eigen_diff 4190 END DO 4191 4192 END IF 4193 4194 CALL timestop(handle) 4195 4196 END SUBROUTINE shift_unshifted_levels 4197 4198! ************************************************************************************************** 4199!> \brief Calculate the matrix mat_N_gw containing the second derivatives 4200!> with respect to the fitting parameters. The second derivatives are 4201!> calculated numerically by finite differences. 4202!> \param N_ij matrix element 4203!> \param Lambda fitting parameters 4204!> \param Sigma_c ... 4205!> \param vec_omega_fit_gw ... 4206!> \param i ... 4207!> \param j ... 4208!> \param num_poles ... 4209!> \param num_fit_points ... 4210!> \param n_level_gw ... 4211!> \param h ... 4212! ************************************************************************************************** 4213 SUBROUTINE calc_mat_N(N_ij, Lambda, Sigma_c, vec_omega_fit_gw, i, j, & 4214 num_poles, num_fit_points, n_level_gw, h) 4215 REAL(KIND=dp), INTENT(OUT) :: N_ij 4216 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:), & 4217 INTENT(IN) :: Lambda 4218 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Sigma_c 4219 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 4220 INTENT(IN) :: vec_omega_fit_gw 4221 INTEGER, INTENT(IN) :: i, j, num_poles, num_fit_points, & 4222 n_level_gw 4223 REAL(KIND=dp), INTENT(IN) :: h 4224 4225 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mat_N', routineP = moduleN//':'//routineN 4226 4227 COMPLEX(KIND=dp) :: im_unit, re_unit 4228 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Lambda_tmp 4229 INTEGER :: handle, num_var 4230 REAL(KIND=dp) :: chi2, chi2_sum 4231 4232 CALL timeset(routineN, handle) 4233 4234 num_var = 2*num_poles + 1 4235 ALLOCATE (Lambda_tmp(num_var)) 4236 Lambda_tmp = (0.0_dp, 0.0_dp) 4237 chi2_sum = 0.0_dp 4238 re_unit = (1.0_dp, 0.0_dp) 4239 im_unit = (0.0_dp, 1.0_dp) 4240 4241 !test 4242 Lambda_tmp(:) = Lambda(:) 4243 CALL calc_chi2(chi2, Lambda_tmp, Sigma_c, vec_omega_fit_gw, num_poles, & 4244 num_fit_points, n_level_gw) 4245 4246 ! Fitting parameters with offset h 4247 Lambda_tmp(:) = Lambda(:) 4248 IF (MODULO(i, 2) == 0) THEN 4249 Lambda_tmp(i/2) = Lambda_tmp(i/2) + h*re_unit 4250 ELSE 4251 Lambda_tmp((i + 1)/2) = Lambda_tmp((i + 1)/2) + h*im_unit 4252 END IF 4253 IF (MODULO(j, 2) == 0) THEN 4254 Lambda_tmp(j/2) = Lambda_tmp(j/2) + h*re_unit 4255 ELSE 4256 Lambda_tmp((j + 1)/2) = Lambda_tmp((j + 1)/2) + h*im_unit 4257 END IF 4258 CALL calc_chi2(chi2, Lambda_tmp, Sigma_c, vec_omega_fit_gw, num_poles, & 4259 num_fit_points, n_level_gw) 4260 chi2_sum = chi2_sum + chi2 4261 4262 IF (MODULO(i, 2) == 0) THEN 4263 Lambda_tmp(i/2) = Lambda_tmp(i/2) - 2.0_dp*h*re_unit 4264 ELSE 4265 Lambda_tmp((i + 1)/2) = Lambda_tmp((i + 1)/2) - 2.0_dp*h*im_unit 4266 END IF 4267 CALL calc_chi2(chi2, Lambda_tmp, Sigma_c, vec_omega_fit_gw, num_poles, & 4268 num_fit_points, n_level_gw) 4269 chi2_sum = chi2_sum - chi2 4270 4271 IF (MODULO(j, 2) == 0) THEN 4272 Lambda_tmp(j/2) = Lambda_tmp(j/2) - 2.0_dp*h*re_unit 4273 ELSE 4274 Lambda_tmp((j + 1)/2) = Lambda_tmp((j + 1)/2) - 2.0_dp*h*im_unit 4275 END IF 4276 CALL calc_chi2(chi2, Lambda_tmp, Sigma_c, vec_omega_fit_gw, num_poles, & 4277 num_fit_points, n_level_gw) 4278 chi2_sum = chi2_sum + chi2 4279 4280 IF (MODULO(i, 2) == 0) THEN 4281 Lambda_tmp(i/2) = Lambda_tmp(i/2) + 2.0_dp*h*re_unit 4282 ELSE 4283 Lambda_tmp((i + 1)/2) = Lambda_tmp((i + 1)/2) + 2.0_dp*h*im_unit 4284 END IF 4285 CALL calc_chi2(chi2, Lambda_tmp, Sigma_c, vec_omega_fit_gw, num_poles, & 4286 num_fit_points, n_level_gw) 4287 chi2_sum = chi2_sum - chi2 4288 4289 ! Second derivative with symmetric difference quotient 4290 N_ij = 1.0_dp/2.0_dp*chi2_sum/(4.0_dp*h*h) 4291 4292 DEALLOCATE (Lambda_tmp) 4293 4294 CALL timestop(handle) 4295 4296 END SUBROUTINE calc_mat_N 4297 4298! ************************************************************************************************** 4299!> \brief Calculate chi2 4300!> \param chi2 ... 4301!> \param Lambda fitting parameters 4302!> \param Sigma_c ... 4303!> \param vec_omega_fit_gw ... 4304!> \param num_poles ... 4305!> \param num_fit_points ... 4306!> \param n_level_gw ... 4307! ************************************************************************************************** 4308 SUBROUTINE calc_chi2(chi2, Lambda, Sigma_c, vec_omega_fit_gw, num_poles, & 4309 num_fit_points, n_level_gw) 4310 REAL(KIND=dp), INTENT(INOUT) :: chi2 4311 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:), & 4312 INTENT(IN) :: Lambda 4313 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Sigma_c 4314 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 4315 INTENT(IN) :: vec_omega_fit_gw 4316 INTEGER, INTENT(IN) :: num_poles, num_fit_points, n_level_gw 4317 4318 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_chi2', routineP = moduleN//':'//routineN 4319 4320 COMPLEX(KIND=dp) :: func_val, im_unit 4321 INTEGER :: handle, iii, jjj, kkk 4322 4323 CALL timeset(routineN, handle) 4324 4325 im_unit = (0.0_dp, 1.0_dp) 4326 chi2 = 0.0_dp 4327 DO kkk = 1, num_fit_points 4328 func_val = Lambda(1) 4329 DO iii = 1, num_poles 4330 jjj = iii*2 4331 ! calculate value of the fit function 4332 func_val = func_val + Lambda(jjj)/(im_unit*vec_omega_fit_gw(kkk) - Lambda(jjj + 1)) 4333 END DO 4334 chi2 = chi2 + (ABS(Sigma_c(n_level_gw, kkk) - func_val))**2 4335 END DO 4336 4337 CALL timestop(handle) 4338 4339 END SUBROUTINE calc_chi2 4340 4341! ************************************************************************************************** 4342!> \brief ... 4343!> \param Eigenval ... 4344!> \param Eigenval_scf ... 4345!> \param ic_corr_list ... 4346!> \param gw_corr_lev_occ ... 4347!> \param gw_corr_lev_virt ... 4348!> \param gw_corr_lev_tot ... 4349!> \param homo ... 4350!> \param nmo ... 4351!> \param unit_nr ... 4352!> \param do_alpha ... 4353!> \param do_beta ... 4354! ************************************************************************************************** 4355 SUBROUTINE apply_ic_corr(Eigenval, Eigenval_scf, ic_corr_list, & 4356 gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, & 4357 homo, nmo, unit_nr, do_alpha, do_beta) 4358 4359 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: Eigenval, Eigenval_scf 4360 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: ic_corr_list 4361 INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, & 4362 gw_corr_lev_tot, homo, nmo, unit_nr 4363 LOGICAL, INTENT(IN), OPTIONAL :: do_alpha, do_beta 4364 4365 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_ic_corr', routineP = moduleN//':'//routineN 4366 4367 CHARACTER(4) :: occ_virt 4368 INTEGER :: handle, n_level_gw, n_level_gw_ref 4369 LOGICAL :: do_closed_shell, my_do_alpha, my_do_beta 4370 REAL(KIND=dp) :: eigen_diff 4371 4372 CALL timeset(routineN, handle) 4373 4374 IF (PRESENT(do_alpha)) THEN 4375 my_do_alpha = do_alpha 4376 ELSE 4377 my_do_alpha = .FALSE. 4378 END IF 4379 4380 IF (PRESENT(do_beta)) THEN 4381 my_do_beta = do_beta 4382 ELSE 4383 my_do_beta = .FALSE. 4384 END IF 4385 4386 do_closed_shell = .NOT. (my_do_alpha .OR. my_do_beta) 4387 4388 ! check the number of input image charge corrected levels 4389 CPASSERT(SIZE(ic_corr_list) == gw_corr_lev_tot) 4390 4391 IF (unit_nr > 0) THEN 4392 4393 WRITE (unit_nr, *) ' ' 4394 4395 IF (do_closed_shell) THEN 4396 WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies with image charge (ic) correction' 4397 WRITE (unit_nr, '(T3,A)') '-----------------------------------------------------------' 4398 ELSE IF (my_do_alpha) THEN 4399 WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies of alpha spins with image charge (ic) correction' 4400 WRITE (unit_nr, '(T3,A)') '--------------------------------------------------------------------------' 4401 ELSE IF (my_do_beta) THEN 4402 WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies of beta spins with image charge (ic) correction' 4403 WRITE (unit_nr, '(T3,A)') '-------------------------------------------------------------------------' 4404 END IF 4405 4406 WRITE (unit_nr, *) ' ' 4407 4408 DO n_level_gw = 1, gw_corr_lev_tot 4409 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 4410 IF (n_level_gw <= gw_corr_lev_occ) THEN 4411 occ_virt = 'occ' 4412 ELSE 4413 occ_virt = 'vir' 4414 END IF 4415 4416 WRITE (unit_nr, '(T4,I4,3A,3F21.3)') & 4417 n_level_gw_ref, ' ( ', occ_virt, ') ', & 4418 Eigenval(n_level_gw_ref)*evolt, & 4419 ic_corr_list(n_level_gw)*evolt, & 4420 (Eigenval(n_level_gw_ref) + ic_corr_list(n_level_gw))*evolt 4421 4422 END DO 4423 4424 WRITE (unit_nr, *) ' ' 4425 4426 END IF 4427 4428 Eigenval(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt) = Eigenval(homo - gw_corr_lev_occ + 1: & 4429 homo + gw_corr_lev_virt) & 4430 + ic_corr_list(1:gw_corr_lev_tot) 4431 4432 Eigenval_scf(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt) = Eigenval_scf(homo - gw_corr_lev_occ + 1: & 4433 homo + gw_corr_lev_virt) & 4434 + ic_corr_list(1:gw_corr_lev_tot) 4435 4436 IF (unit_nr > 0) THEN 4437 4438 IF (do_closed_shell) THEN 4439 WRITE (unit_nr, '(T3,A,F52.2)') 'G0W0 IC HOMO-LUMO gap (eV)', Eigenval(homo + 1) - Eigenval(homo) 4440 ELSE IF (my_do_alpha) THEN 4441 WRITE (unit_nr, '(T3,A,F46.2)') 'G0W0 Alpha IC HOMO-LUMO gap (eV)', Eigenval(homo + 1) - Eigenval(homo) 4442 ELSE IF (my_do_beta) THEN 4443 WRITE (unit_nr, '(T3,A,F47.2)') 'G0W0 Beta IC HOMO-LUMO gap (eV)', Eigenval(homo + 1) - Eigenval(homo) 4444 END IF 4445 4446 WRITE (unit_nr, *) ' ' 4447 4448 END IF 4449 4450 ! for eigenvalue self-consistent GW, all eigenvalues have to be corrected 4451 ! 1) the occupied; check if there are occupied MOs not being corrected by the IC model 4452 IF (gw_corr_lev_occ < homo .AND. gw_corr_lev_occ > 0) THEN 4453 4454 ! calculate average IC contribution for occupied orbitals 4455 eigen_diff = 0.0_dp 4456 4457 DO n_level_gw = 1, gw_corr_lev_occ 4458 eigen_diff = eigen_diff + ic_corr_list(n_level_gw) 4459 END DO 4460 eigen_diff = eigen_diff/gw_corr_lev_occ 4461 4462 ! correct the eigenvalues of the occupied orbitals which have not been corrected by the IC model 4463 DO n_level_gw = 1, homo - gw_corr_lev_occ 4464 Eigenval(n_level_gw) = Eigenval(n_level_gw) + eigen_diff 4465 Eigenval_scf(n_level_gw) = Eigenval_scf(n_level_gw) + eigen_diff 4466 END DO 4467 4468 END IF 4469 4470 ! 2) the virtual: check if there are virtual orbitals not being corrected by the IC model 4471 IF (gw_corr_lev_virt < nmo - homo .AND. gw_corr_lev_virt > 0) THEN 4472 4473 ! calculate average IC correction for virtual orbitals 4474 eigen_diff = 0.0_dp 4475 DO n_level_gw = gw_corr_lev_occ + 1, gw_corr_lev_tot 4476 eigen_diff = eigen_diff + ic_corr_list(n_level_gw) 4477 END DO 4478 eigen_diff = eigen_diff/gw_corr_lev_virt 4479 4480 ! correct the eigenvalues of the virtual orbitals which have not been corrected by the IC model 4481 DO n_level_gw = homo + gw_corr_lev_virt + 1, nmo 4482 Eigenval(n_level_gw) = Eigenval(n_level_gw) + eigen_diff 4483 Eigenval_scf(n_level_gw) = Eigenval_scf(n_level_gw) + eigen_diff 4484 END DO 4485 4486 END IF 4487 4488 CALL timestop(handle) 4489 4490 END SUBROUTINE apply_ic_corr 4491 4492! ************************************************************************************************** 4493!> \brief ... 4494!> \param num_integ_points ... 4495!> \param nmo ... 4496!> \param tau_tj ... 4497!> \param tj ... 4498!> \param matrix_s ... 4499!> \param fm_mo_coeff_occ ... 4500!> \param fm_mo_coeff_virt ... 4501!> \param fm_mo_coeff_occ_scaled ... 4502!> \param fm_mo_coeff_virt_scaled ... 4503!> \param fm_scaled_dm_occ_tau ... 4504!> \param fm_scaled_dm_virt_tau ... 4505!> \param Eigenval ... 4506!> \param eps_filter ... 4507!> \param e_fermi ... 4508!> \param fm_mat_W ... 4509!> \param gw_corr_lev_tot ... 4510!> \param gw_corr_lev_occ ... 4511!> \param gw_corr_lev_virt ... 4512!> \param homo ... 4513!> \param count_ev_sc_GW ... 4514!> \param count_sc_GW0 ... 4515!> \param t_3c_overl_int_gw_RI ... 4516!> \param t_3c_overl_int_gw_AO ... 4517!> \param mat_W ... 4518!> \param mat_SinvVSinv ... 4519!> \param mat_dm ... 4520!> \param weights_cos_tf_t_to_w ... 4521!> \param weights_sin_tf_t_to_w ... 4522!> \param vec_Sigma_c_gw ... 4523!> \param do_periodic ... 4524!> \param num_points_corr ... 4525!> \param delta_corr ... 4526!> \param qs_env ... 4527!> \param para_env ... 4528!> \param para_env_RPA ... 4529!> \param mp2_env ... 4530!> \param matrix_berry_re_mo_mo ... 4531!> \param matrix_berry_im_mo_mo ... 4532!> \param first_cycle_periodic_correction ... 4533!> \param kpoints ... 4534!> \param num_fit_points ... 4535!> \param fm_mo_coeff ... 4536!> \param do_ri_Sigma_x ... 4537!> \param vec_Sigma_x_gw ... 4538!> \param unit_nr ... 4539!> \param do_beta ... 4540! ************************************************************************************************** 4541 SUBROUTINE compute_self_energy_cubic_gw(num_integ_points, nmo, tau_tj, tj, & 4542 matrix_s, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, & 4543 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, & 4544 fm_scaled_dm_virt_tau, Eigenval, eps_filter, & 4545 e_fermi, fm_mat_W, & 4546 gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, & 4547 count_ev_sc_GW, count_sc_GW0, & 4548 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, & 4549 mat_W, mat_SinvVSinv, mat_dm, & 4550 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_Sigma_c_gw, & 4551 do_periodic, num_points_corr, delta_corr, qs_env, para_env, para_env_RPA, & 4552 mp2_env, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, & 4553 first_cycle_periodic_correction, kpoints, num_fit_points, fm_mo_coeff, & 4554 do_ri_Sigma_x, vec_Sigma_x_gw, unit_nr, do_beta) 4555 INTEGER, INTENT(IN) :: num_integ_points, nmo 4556 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 4557 INTENT(IN) :: tau_tj, tj 4558 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 4559 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, & 4560 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau 4561 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 4562 REAL(KIND=dp), INTENT(IN) :: eps_filter 4563 REAL(KIND=dp), INTENT(INOUT) :: e_fermi 4564 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fm_mat_W 4565 INTEGER, INTENT(IN) :: gw_corr_lev_tot, gw_corr_lev_occ, & 4566 gw_corr_lev_virt, homo, & 4567 count_ev_sc_GW, count_sc_GW0 4568 TYPE(dbcsr_t_type) :: t_3c_overl_int_gw_RI, & 4569 t_3c_overl_int_gw_AO 4570 TYPE(dbcsr_type), POINTER :: mat_W 4571 TYPE(dbcsr_p_type) :: mat_SinvVSinv, mat_dm 4572 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 4573 INTENT(IN) :: weights_cos_tf_t_to_w, & 4574 weights_sin_tf_t_to_w 4575 COMPLEX(KIND=dp), ALLOCATABLE, & 4576 DIMENSION(:, :, :), INTENT(INOUT) :: vec_Sigma_c_gw 4577 LOGICAL, INTENT(IN) :: do_periodic 4578 INTEGER, INTENT(IN) :: num_points_corr 4579 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 4580 INTENT(INOUT) :: delta_corr 4581 TYPE(qs_environment_type), POINTER :: qs_env 4582 TYPE(cp_para_env_type), POINTER :: para_env, para_env_RPA 4583 TYPE(mp2_type), POINTER :: mp2_env 4584 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_re_mo_mo, & 4585 matrix_berry_im_mo_mo 4586 LOGICAL, INTENT(INOUT) :: first_cycle_periodic_correction 4587 TYPE(kpoint_type), POINTER :: kpoints 4588 INTEGER, INTENT(IN) :: num_fit_points 4589 TYPE(cp_fm_type), POINTER :: fm_mo_coeff 4590 LOGICAL, INTENT(IN) :: do_ri_Sigma_x 4591 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 4592 INTENT(INOUT) :: vec_Sigma_x_gw 4593 INTEGER, INTENT(IN) :: unit_nr 4594 LOGICAL, INTENT(IN), OPTIONAL :: do_beta 4595 4596 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_self_energy_cubic_gw', & 4597 routineP = moduleN//':'//routineN 4598 4599 COMPLEX(KIND=dp) :: im_unit 4600 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta_corr_omega 4601 INTEGER :: gw_lev_end, gw_lev_start, handle, & 4602 handle3, iquad, jquad, mo_end, & 4603 mo_start, n_level_gw, n_level_gw_ref, & 4604 unit_nr_prv 4605 LOGICAL :: memory_info, my_do_beta 4606 REAL(KIND=dp) :: ext_scaling, omega, omega_i, omega_sign, & 4607 sign_occ_virt, t_i_Clenshaw, tau, & 4608 weight_cos, weight_i, weight_sin 4609 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_Sigma_c_gw_cos_omega, & 4610 vec_Sigma_c_gw_cos_tau, vec_Sigma_c_gw_neg_tau, vec_Sigma_c_gw_pos_tau, & 4611 vec_Sigma_c_gw_sin_omega, vec_Sigma_c_gw_sin_tau 4612 TYPE(dbcsr_t_type) :: t_3c_ctr_RI 4613 TYPE(dbcsr_type), POINTER :: mat_greens_fct_occ, mat_greens_fct_virt 4614 4615 CALL timeset(routineN, handle) 4616 4617 memory_info = mp2_env%ri_rpa_im_time%memory_info 4618 IF (memory_info) THEN 4619 unit_nr_prv = unit_nr 4620 ELSE 4621 unit_nr_prv = 0 4622 ENDIF 4623 4624 my_do_beta = .FALSE. 4625 IF (PRESENT(do_beta)) my_do_beta = do_beta 4626 4627 im_unit = (0.0_dp, 1.0_dp) 4628 4629 mo_start = homo - gw_corr_lev_occ + 1 4630 mo_end = homo + gw_corr_lev_virt 4631 CPASSERT(mo_end - mo_start + 1 == gw_corr_lev_tot) 4632 4633 vec_Sigma_c_gw = (0.0_dp, 0.0_dp) 4634 ALLOCATE (vec_Sigma_c_gw_pos_tau(gw_corr_lev_tot, num_integ_points)) 4635 vec_Sigma_c_gw_pos_tau = 0.0_dp 4636 ALLOCATE (vec_Sigma_c_gw_neg_tau(gw_corr_lev_tot, num_integ_points)) 4637 vec_Sigma_c_gw_neg_tau = 0.0_dp 4638 ALLOCATE (vec_Sigma_c_gw_cos_tau(gw_corr_lev_tot, num_integ_points)) 4639 vec_Sigma_c_gw_cos_tau = 0.0_dp 4640 ALLOCATE (vec_Sigma_c_gw_sin_tau(gw_corr_lev_tot, num_integ_points)) 4641 vec_Sigma_c_gw_sin_tau = 0.0_dp 4642 4643 ALLOCATE (vec_Sigma_c_gw_cos_omega(gw_corr_lev_tot, num_integ_points)) 4644 vec_Sigma_c_gw_cos_omega = 0.0_dp 4645 ALLOCATE (vec_Sigma_c_gw_sin_omega(gw_corr_lev_tot, num_integ_points)) 4646 vec_Sigma_c_gw_sin_omega = 0.0_dp 4647 4648 ALLOCATE (delta_corr_omega(1 + homo - gw_corr_lev_occ:homo + gw_corr_lev_virt, num_integ_points)) 4649 delta_corr_omega(:, :) = (0.0_dp, 0.0_dp) 4650 4651 NULLIFY (mat_greens_fct_occ) 4652 CALL dbcsr_init_p(mat_greens_fct_occ) 4653 CALL dbcsr_create(matrix=mat_greens_fct_occ, & 4654 template=matrix_s(1)%matrix, & 4655 matrix_type=dbcsr_type_no_symmetry) 4656 4657 NULLIFY (mat_greens_fct_virt) 4658 CALL dbcsr_init_p(mat_greens_fct_virt) 4659 CALL dbcsr_create(matrix=mat_greens_fct_virt, & 4660 template=matrix_s(1)%matrix, & 4661 matrix_type=dbcsr_type_no_symmetry) 4662 4663 e_fermi = 0.5_dp*(Eigenval(homo) + Eigenval(homo + 1)) 4664 4665 DO jquad = 1, num_integ_points 4666 4667 CALL compute_Greens_function_time(mat_greens_fct_occ, mat_greens_fct_virt, & 4668 fm_mo_coeff_occ, fm_mo_coeff_virt, & 4669 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, & 4670 fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, Eigenval, & 4671 nmo, eps_filter, e_fermi, tau_tj(jquad)) 4672 4673 CALL dbcsr_set(mat_W, 0.0_dp) 4674 CALL copy_fm_to_dbcsr(fm_mat_W(jquad)%matrix, mat_W, keep_sparsity=.FALSE.) 4675 4676 CALL compute_sigma_gw(t_3c_overl_int_gw_AO, t_3c_overl_int_gw_RI, & 4677 mat_greens_fct_occ, mat_W, [1.0_dp, -1.0_dp], & 4678 vec_Sigma_c_gw_neg_tau(:, jquad), [mo_start, mo_end], para_env, unit_nr_prv, & 4679 t_3c_ctr_RI=t_3c_ctr_RI, t_3c_ctr_in=.FALSE.) 4680 4681 CALL compute_sigma_gw(t_3c_overl_int_gw_AO, t_3c_overl_int_gw_RI, & 4682 mat_greens_fct_virt, mat_W, [1.0_dp, 1.0_dp], & 4683 vec_Sigma_c_gw_pos_tau(:, jquad), [mo_start, mo_end], para_env, unit_nr_prv, & 4684 t_3c_ctr_RI=t_3c_ctr_RI, t_3c_ctr_in=.TRUE.) 4685 4686 CALL dbcsr_t_destroy(t_3c_ctr_RI) 4687 4688 vec_Sigma_c_gw_cos_tau(:, jquad) = 0.5_dp*(vec_Sigma_c_gw_pos_tau(:, jquad) + & 4689 vec_Sigma_c_gw_neg_tau(:, jquad)) 4690 4691 vec_Sigma_c_gw_sin_tau(:, jquad) = 0.5_dp*(vec_Sigma_c_gw_pos_tau(:, jquad) - & 4692 vec_Sigma_c_gw_neg_tau(:, jquad)) 4693 4694 END DO ! jquad (tau) 4695 4696 ! Fourier transform from time to frequency 4697 DO jquad = 1, num_fit_points 4698 4699 DO iquad = 1, num_integ_points 4700 4701 omega = tj(jquad) 4702 tau = tau_tj(iquad) 4703 weight_cos = weights_cos_tf_t_to_w(jquad, iquad)*COS(omega*tau) 4704 weight_sin = weights_sin_tf_t_to_w(jquad, iquad)*SIN(omega*tau) 4705 4706 vec_Sigma_c_gw_cos_omega(:, jquad) = vec_Sigma_c_gw_cos_omega(:, jquad) + & 4707 weight_cos*vec_Sigma_c_gw_cos_tau(:, iquad) 4708 4709 vec_Sigma_c_gw_sin_omega(:, jquad) = vec_Sigma_c_gw_sin_omega(:, jquad) + & 4710 weight_sin*vec_Sigma_c_gw_sin_tau(:, iquad) 4711 4712 END DO 4713 4714 END DO 4715 4716 ! for occupied levels, we need the correlation self-energy for negative omega. Therefore, weight_sin 4717 ! should be computed with -omega, which results in an additional minus for vec_Sigma_c_gw_sin_omega: 4718 vec_Sigma_c_gw_sin_omega(1:gw_corr_lev_occ, :) = -vec_Sigma_c_gw_sin_omega(1:gw_corr_lev_occ, :) 4719 4720 vec_Sigma_c_gw(:, 1:num_fit_points, 1) = vec_Sigma_c_gw_cos_omega(:, 1:num_fit_points) + & 4721 im_unit*vec_Sigma_c_gw_sin_omega(:, 1:num_fit_points) 4722 4723 CALL dbcsr_release_P(mat_greens_fct_occ) 4724 CALL dbcsr_release_P(mat_greens_fct_virt) 4725 4726 IF (do_ri_Sigma_x .AND. count_ev_sc_GW == 1 .AND. count_sc_GW0 == 1) THEN 4727 4728 CALL timeset(routineN//"_RI_HFX_operation_1", handle3) 4729 4730 ! get density matrix 4731 CALL cp_gemm(transa="N", transb="T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, & 4732 matrix_a=fm_mo_coeff_occ, matrix_b=fm_mo_coeff_occ, beta=0.0_dp, & 4733 matrix_c=fm_scaled_dm_occ_tau) 4734 4735 CALL timestop(handle3) 4736 4737 CALL timeset(routineN//"_RI_HFX_operation_2", handle3) 4738 4739 CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, & 4740 mat_dm%matrix, & 4741 keep_sparsity=.FALSE.) 4742 4743 CALL timestop(handle3) 4744 4745 CALL compute_sigma_gw(t_3c_overl_int_gw_AO, t_3c_overl_int_gw_RI, & 4746 mat_dm%matrix, mat_SinvVSinv%matrix, [1.0_dp, -1.0_dp], & 4747 vec_Sigma_x_gw(mo_start:mo_end, 1), [mo_start, mo_end], para_env, unit_nr_prv) 4748 4749 IF (my_do_beta) THEN 4750 4751 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, 1) = & 4752 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, 1) + & 4753 vec_Sigma_x_gw(:, 1) 4754 4755 ELSE 4756 4757 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, 1) = & 4758 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, 1) + & 4759 vec_Sigma_x_gw(:, 1) 4760 4761 END IF 4762 4763 END IF 4764 4765 ! compute and add the periodic correction 4766 IF (do_periodic) THEN 4767 4768 ext_scaling = 0.2_dp 4769 4770 ! loop over omega' (integration) 4771 DO iquad = 1, num_points_corr 4772 4773 ! use the Clenshaw-grid 4774 t_i_Clenshaw = iquad*pi/(2.0_dp*num_points_corr) 4775 omega_i = ext_scaling/TAN(t_i_Clenshaw) 4776 4777 IF (iquad < num_points_corr) THEN 4778 weight_i = ext_scaling*pi/(num_points_corr*SIN(t_i_Clenshaw)**2) 4779 ELSE 4780 weight_i = ext_scaling*pi/(2.0_dp*num_points_corr*SIN(t_i_Clenshaw)**2) 4781 END IF 4782 4783 CALL calc_periodic_correction(delta_corr, qs_env, para_env, para_env_RPA, & 4784 mp2_env%ri_g0w0%kp_grid, homo, nmo, gw_corr_lev_occ, & 4785 gw_corr_lev_virt, omega_i, fm_mo_coeff, Eigenval, & 4786 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, & 4787 first_cycle_periodic_correction, kpoints, & 4788 mp2_env%ri_g0w0%do_mo_coeff_gamma, & 4789 mp2_env%ri_g0w0%num_kp_grids, mp2_env%ri_g0w0%eps_kpoint, & 4790 mp2_env%ri_g0w0%do_extra_kpoints, & 4791 mp2_env%ri_g0w0%do_aux_bas_gw, mp2_env%ri_g0w0%frac_aux_mos) 4792 4793 DO n_level_gw = 1, gw_corr_lev_tot 4794 4795 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ 4796 4797 IF (n_level_gw <= gw_corr_lev_occ) THEN 4798 sign_occ_virt = -1.0_dp 4799 ELSE 4800 sign_occ_virt = 1.0_dp 4801 END IF 4802 4803 DO jquad = 1, num_integ_points 4804 4805 omega_sign = tj(jquad)*sign_occ_virt 4806 4807 delta_corr_omega(n_level_gw_ref, jquad) = & 4808 delta_corr_omega(n_level_gw_ref, jquad) - & 4809 0.5_dp/pi*weight_i/2.0_dp*delta_corr(n_level_gw_ref)* & 4810 (1.0_dp/(im_unit*(omega_i + omega_sign) + e_fermi - Eigenval(n_level_gw_ref)) + & 4811 1.0_dp/(im_unit*(-omega_i + omega_sign) + e_fermi - Eigenval(n_level_gw_ref))) 4812 4813 END DO 4814 4815 END DO 4816 4817 END DO 4818 4819 gw_lev_start = 1 + homo - gw_corr_lev_occ 4820 gw_lev_end = homo + gw_corr_lev_virt 4821 4822 ! add the periodic correction 4823 vec_Sigma_c_gw(1:gw_corr_lev_tot, :, 1) = vec_Sigma_c_gw(1:gw_corr_lev_tot, :, 1) + & 4824 delta_corr_omega(gw_lev_start:gw_lev_end, 1:num_fit_points) 4825 4826 END IF 4827 4828 DEALLOCATE (vec_Sigma_c_gw_pos_tau) 4829 DEALLOCATE (vec_Sigma_c_gw_neg_tau) 4830 DEALLOCATE (vec_Sigma_c_gw_cos_tau) 4831 DEALLOCATE (vec_Sigma_c_gw_sin_tau) 4832 DEALLOCATE (vec_Sigma_c_gw_cos_omega) 4833 DEALLOCATE (vec_Sigma_c_gw_sin_omega) 4834 DEALLOCATE (delta_corr_omega) 4835 4836 CALL timestop(handle) 4837 4838 END SUBROUTINE compute_self_energy_cubic_gw 4839 4840! ************************************************************************************************** 4841!> \brief ... 4842!> \param t_3c_overl_int_gw_AO ... 4843!> \param t_3c_overl_int_gw_RI ... 4844!> \param mat_AO ... 4845!> \param mat_RI ... 4846!> \param prefac ... 4847!> \param vec_Sigma ... 4848!> \param mo_bounds ... 4849!> \param para_env ... 4850!> \param unit_nr ... 4851!> \param t_3c_ctr_RI ... 4852!> \param t_3c_ctr_in ... 4853! ************************************************************************************************** 4854 SUBROUTINE compute_sigma_gw(t_3c_overl_int_gw_AO, t_3c_overl_int_gw_RI, & 4855 mat_AO, mat_RI, prefac, & 4856 vec_Sigma, mo_bounds, para_env, unit_nr, & 4857 t_3c_ctr_RI, t_3c_ctr_in) 4858 TYPE(dbcsr_t_type), INTENT(INOUT) :: t_3c_overl_int_gw_AO, & 4859 t_3c_overl_int_gw_RI 4860 TYPE(dbcsr_type), INTENT(IN) :: mat_AO, mat_RI 4861 REAL(dp), DIMENSION(2), INTENT(IN) :: prefac 4862 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec_Sigma 4863 INTEGER, DIMENSION(2), INTENT(IN) :: mo_bounds 4864 TYPE(cp_para_env_type), POINTER :: para_env 4865 INTEGER, INTENT(IN) :: unit_nr 4866 TYPE(dbcsr_t_type), INTENT(INOUT), OPTIONAL :: t_3c_ctr_RI 4867 LOGICAL, INTENT(IN), OPTIONAL :: t_3c_ctr_in 4868 4869 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_sigma_gw', & 4870 routineP = moduleN//':'//routineN 4871 4872 INTEGER :: handle 4873 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist1, dist2, sizes_AO, sizes_RI 4874 INTEGER, DIMENSION(2) :: pdims_2d 4875 LOGICAL :: t_3c_ctr_in_prv 4876 TYPE(dbcsr_t_pgrid_type) :: pgrid_2d 4877 TYPE(dbcsr_t_type) :: t_3c_ctr_RI_prv, t_3c_overl_int_gw_copy, & 4878 t_AO, t_AO_tmp, t_RI, t_RI_tmp 4879 4880 CALL timeset(routineN, handle) 4881 4882 IF (PRESENT(t_3c_ctr_in)) THEN 4883 t_3c_ctr_in_prv = t_3c_ctr_in 4884 ELSE 4885 t_3c_ctr_in_prv = .FALSE. 4886 ENDIF 4887 4888 pdims_2d = 0 4889 CALL dbcsr_t_pgrid_create(para_env%group, pdims_2d, pgrid_2d) 4890 4891 IF (t_3c_ctr_in_prv) THEN 4892 CPASSERT(PRESENT(t_3c_ctr_RI)) 4893 t_3c_ctr_RI_prv = t_3c_ctr_RI 4894 ELSE 4895 4896 ALLOCATE (sizes_RI(dbcsr_t_nblks_total(t_3c_overl_int_gw_RI, 1))) 4897 CALL dbcsr_t_get_info(t_3c_overl_int_gw_RI, blk_size_1=sizes_RI) 4898 4899 CALL create_2c_tensor(t_RI, dist1, dist2, pgrid_2d, sizes_RI, sizes_RI, name="(RI|RI)") 4900 DEALLOCATE (dist1, dist2, sizes_RI) 4901 4902 CALL dbcsr_t_create(mat_RI, t_RI_tmp, name="(RI|RI)") 4903 CALL dbcsr_t_copy_matrix_to_tensor(mat_RI, t_RI_tmp) 4904 CALL dbcsr_t_copy(t_RI_tmp, t_RI) 4905 CALL dbcsr_t_destroy(t_RI_tmp) 4906 4907 CALL dbcsr_t_create(t_3c_overl_int_gw_RI, t_3c_overl_int_gw_copy) 4908 4909 CALL dbcsr_t_contract(dbcsr_scalar(prefac(1)), t_RI, t_3c_overl_int_gw_RI, dbcsr_scalar(0.0_dp), & 4910 t_3c_overl_int_gw_copy, & 4911 contract_1=[2], notcontract_1=[1], & 4912 contract_2=[1], notcontract_2=[2, 3], & 4913 map_1=[1], map_2=[2, 3], & 4914 unit_nr=unit_nr) 4915 4916 CALL dbcsr_t_create(t_3c_overl_int_gw_AO, t_3c_ctr_RI_prv) 4917 CALL dbcsr_t_copy(t_3c_overl_int_gw_copy, t_3c_ctr_RI_prv, order=[2, 1, 3]) 4918 CALL dbcsr_t_destroy(t_3c_overl_int_gw_copy) 4919 CALL dbcsr_t_destroy(t_RI) 4920 4921 ENDIF 4922 4923 ALLOCATE (sizes_AO(dbcsr_t_nblks_total(t_3c_overl_int_gw_AO, 1))) 4924 CALL dbcsr_t_get_info(t_3c_overl_int_gw_AO, blk_size_1=sizes_AO) 4925 CALL create_2c_tensor(t_AO, dist1, dist2, pgrid_2d, sizes_AO, sizes_AO, name="(AO|AO)") 4926 DEALLOCATE (dist1, dist2, sizes_AO) 4927 4928 CALL dbcsr_t_create(mat_AO, t_AO_tmp, name="(AO|AO)") 4929 CALL dbcsr_t_copy_matrix_to_tensor(mat_AO, t_AO_tmp) 4930 CALL dbcsr_t_copy(t_AO_tmp, t_AO) 4931 CALL dbcsr_t_destroy(t_AO_tmp) 4932 4933 CALL dbcsr_t_create(t_3c_overl_int_gw_AO, t_3c_overl_int_gw_copy) 4934 CALL dbcsr_t_contract(dbcsr_scalar(prefac(2)), t_AO, t_3c_overl_int_gw_AO, dbcsr_scalar(0.0_dp), & 4935 t_3c_overl_int_gw_copy, & 4936 contract_1=[2], notcontract_1=[1], & 4937 contract_2=[1], notcontract_2=[2, 3], & 4938 map_1=[1], map_2=[2, 3], & 4939 unit_nr=unit_nr) 4940 4941 CALL trace_sigma_gw(t_3c_ctr_RI_prv, t_3c_overl_int_gw_copy, vec_sigma, mo_bounds, para_env) 4942 CALL dbcsr_t_destroy(t_3c_overl_int_gw_copy) 4943 4944 IF (PRESENT(t_3c_ctr_in)) THEN 4945 IF (.NOT. t_3c_ctr_in) THEN 4946 CPASSERT(PRESENT(t_3c_ctr_RI)) 4947 t_3c_ctr_RI = t_3c_ctr_RI_prv 4948 ENDIF 4949 ELSE 4950 CALL dbcsr_t_destroy(t_3c_ctr_RI_prv) 4951 ENDIF 4952 4953 CALL dbcsr_t_pgrid_destroy(pgrid_2d) 4954 CALL dbcsr_t_destroy(t_AO) 4955 4956 CALL timestop(handle) 4957 4958 END SUBROUTINE 4959 4960! ************************************************************************************************** 4961!> \brief ... 4962!> \param t3c_1 ... 4963!> \param t3c_2 ... 4964!> \param vec_sigma ... 4965!> \param mo_bounds ... 4966!> \param para_env ... 4967! ************************************************************************************************** 4968 SUBROUTINE trace_sigma_gw(t3c_1, t3c_2, vec_sigma, mo_bounds, para_env) 4969 TYPE(dbcsr_t_type), INTENT(INOUT) :: t3c_1, t3c_2 4970 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec_Sigma 4971 INTEGER, DIMENSION(2), INTENT(IN) :: mo_bounds 4972 TYPE(cp_para_env_type), INTENT(IN) :: para_env 4973 4974 CHARACTER(LEN=*), PARAMETER :: routineN = 'trace_sigma_gw', routineP = moduleN//':'//routineN 4975 4976 INTEGER :: blk, handle, n, n_end, n_end_block, & 4977 n_start, n_start_block 4978 INTEGER, DIMENSION(1) :: trace_shape 4979 INTEGER, DIMENSION(3) :: boff, bsize, ind 4980 LOGICAL :: found 4981 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: block_1, block_2 4982 TYPE(dbcsr_t_iterator_type) :: iter 4983 4984 CALL timeset(routineN, handle) 4985 4986 CALL dbcsr_t_iterator_start(iter, t3c_1) 4987 DO WHILE (dbcsr_t_iterator_blocks_left(iter)) 4988 CALL dbcsr_t_iterator_next_block(iter, ind, blk, blk_size=bsize, blk_offset=boff) 4989 CALL dbcsr_t_get_block(t3c_1, ind, block_1, found) 4990 CPASSERT(found) 4991 CALL dbcsr_t_get_block(t3c_2, ind, block_2, found) 4992 IF (.NOT. found) CYCLE 4993 4994 IF (boff(3) < mo_bounds(1)) THEN 4995 n_start_block = mo_bounds(1) - boff(3) + 1 4996 n_start = 1 4997 ELSE 4998 n_start_block = 1 4999 n_start = boff(3) - mo_bounds(1) + 1 5000 ENDIF 5001 5002 IF (boff(3) + bsize(3) - 1 > mo_bounds(2)) THEN 5003 n_end_block = mo_bounds(2) - boff(3) + 1 5004 n_end = mo_bounds(2) - mo_bounds(1) + 1 5005 ELSE 5006 n_end_block = bsize(3) 5007 n_end = boff(3) + bsize(3) - mo_bounds(1) 5008 ENDIF 5009 5010 trace_shape(1) = SIZE(block_1, 1)*SIZE(block_1, 2) 5011 vec_Sigma(n_start:n_end) = & 5012 vec_Sigma(n_start:n_end) + & 5013 (/(DOT_PRODUCT(RESHAPE(block_1(:, :, n), trace_shape), & 5014 RESHAPE(block_2(:, :, n), trace_shape)), & 5015 n=n_start_block, n_end_block)/) 5016 DEALLOCATE (block_1, block_2) 5017 ENDDO 5018 CALL dbcsr_t_iterator_stop(iter) 5019 5020 CALL mp_sum(vec_Sigma, para_env%group) 5021 5022 CALL timestop(handle) 5023 END SUBROUTINE 5024 5025! ************************************************************************************************** 5026!> \brief ... 5027!> \param mat_greens_fct_occ ... 5028!> \param mat_greens_fct_virt ... 5029!> \param fm_mo_coeff_occ ... 5030!> \param fm_mo_coeff_virt ... 5031!> \param fm_mo_coeff_occ_scaled ... 5032!> \param fm_mo_coeff_virt_scaled ... 5033!> \param fm_scaled_dm_occ_tau ... 5034!> \param fm_scaled_dm_virt_tau ... 5035!> \param Eigenval ... 5036!> \param nmo ... 5037!> \param eps_filter ... 5038!> \param e_fermi ... 5039!> \param tau ... 5040! ************************************************************************************************** 5041 SUBROUTINE compute_Greens_function_time(mat_greens_fct_occ, mat_greens_fct_virt, fm_mo_coeff_occ, fm_mo_coeff_virt, & 5042 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, & 5043 fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, Eigenval, nmo, & 5044 eps_filter, e_fermi, tau) 5045 5046 TYPE(dbcsr_type), POINTER :: mat_greens_fct_occ, mat_greens_fct_virt 5047 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, & 5048 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau 5049 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 5050 INTEGER, INTENT(IN) :: nmo 5051 REAL(KIND=dp), INTENT(IN) :: eps_filter, e_fermi, tau 5052 5053 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Greens_function_time', & 5054 routineP = moduleN//':'//routineN 5055 5056 INTEGER :: handle, i_global, iiB, jjB, ncol_local, & 5057 nrow_local 5058 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 5059 REAL(KIND=dp) :: stabilize_exp 5060 5061 CALL timeset(routineN, handle) 5062 5063 ! get info of fm_mo_coeff_occ 5064 CALL cp_fm_get_info(matrix=fm_mo_coeff_occ, & 5065 nrow_local=nrow_local, & 5066 ncol_local=ncol_local, & 5067 row_indices=row_indices, & 5068 col_indices=col_indices) 5069 5070 ! Multiply the occupied and the virtual MO coefficients with the factor exp((-e_i-e_F)*tau/2). 5071 ! Then, we simply get the sum over all occ states and virt. states by a simple matrix-matrix 5072 ! multiplication. 5073 5074 stabilize_exp = 70.0_dp 5075 5076 ! first, the occ 5077 DO jjB = 1, nrow_local 5078 DO iiB = 1, ncol_local 5079 i_global = col_indices(iiB) 5080 5081 IF (ABS(tau*0.5_dp*(Eigenval(i_global) - e_fermi)) < stabilize_exp) THEN 5082 fm_mo_coeff_occ_scaled%local_data(jjB, iiB) = & 5083 fm_mo_coeff_occ%local_data(jjB, iiB)*EXP(tau*0.5_dp*(Eigenval(i_global) - e_fermi)) 5084 ELSE 5085 fm_mo_coeff_occ_scaled%local_data(jjB, iiB) = 0.0_dp 5086 END IF 5087 5088 END DO 5089 END DO 5090 5091 ! the same for virt 5092 DO jjB = 1, nrow_local 5093 DO iiB = 1, ncol_local 5094 i_global = col_indices(iiB) 5095 5096 IF (ABS(tau*0.5_dp*(Eigenval(i_global) - e_fermi)) < stabilize_exp) THEN 5097 fm_mo_coeff_virt_scaled%local_data(jjB, iiB) = & 5098 fm_mo_coeff_virt%local_data(jjB, iiB)*EXP(-tau*0.5_dp*(Eigenval(i_global) - e_fermi)) 5099 ELSE 5100 fm_mo_coeff_virt_scaled%local_data(jjB, iiB) = 0.0_dp 5101 END IF 5102 5103 END DO 5104 END DO 5105 5106 CALL cp_gemm(transa="N", transb="T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, & 5107 matrix_a=fm_mo_coeff_occ_scaled, matrix_b=fm_mo_coeff_occ_scaled, beta=0.0_dp, & 5108 matrix_c=fm_scaled_dm_occ_tau) 5109 5110 CALL cp_gemm(transa="N", transb="T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, & 5111 matrix_a=fm_mo_coeff_virt_scaled, matrix_b=fm_mo_coeff_virt_scaled, beta=0.0_dp, & 5112 matrix_c=fm_scaled_dm_virt_tau) 5113 5114 CALL dbcsr_set(mat_greens_fct_occ, 0.0_dp) 5115 5116 CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, & 5117 mat_greens_fct_occ, & 5118 keep_sparsity=.FALSE.) 5119 5120 CALL dbcsr_filter(mat_greens_fct_occ, eps_filter) 5121 5122 CALL dbcsr_set(mat_greens_fct_virt, 0.0_dp) 5123 5124 CALL copy_fm_to_dbcsr(fm_scaled_dm_virt_tau, & 5125 mat_greens_fct_virt, & 5126 keep_sparsity=.FALSE.) 5127 5128 CALL dbcsr_filter(mat_greens_fct_virt, eps_filter) 5129 5130 CALL timestop(handle) 5131 5132 END SUBROUTINE compute_Greens_function_time 5133 5134! ************************************************************************************************** 5135!> \brief ... 5136!> \param mat_greens_fct_occ ... 5137!> \param mat_greens_fct_virt ... 5138!> \param fm_mo_coeff ... 5139!> \param fm_mo_coeff_occ_scaled ... 5140!> \param fm_mo_coeff_virt_scaled ... 5141!> \param fm_scaled_dm_occ_tau ... 5142!> \param fm_scaled_dm_virt_tau ... 5143!> \param Eigenval ... 5144!> \param fermi_level_offset ... 5145!> \param nmo ... 5146!> \param homo ... 5147!> \param eps_filter ... 5148!> \param omega ... 5149!> \param real_imaginary ... 5150! ************************************************************************************************** 5151 SUBROUTINE compute_Greens_function_freq(mat_greens_fct_occ, mat_greens_fct_virt, fm_mo_coeff, & 5152 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, & 5153 fm_scaled_dm_virt_tau, Eigenval, fermi_level_offset, nmo, & 5154 homo, eps_filter, omega, real_imaginary) 5155 5156 TYPE(dbcsr_type), POINTER :: mat_greens_fct_occ, mat_greens_fct_virt 5157 TYPE(cp_fm_type), POINTER :: fm_mo_coeff, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, & 5158 fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau 5159 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 5160 REAL(KIND=dp), INTENT(IN) :: fermi_level_offset 5161 INTEGER, INTENT(IN) :: nmo, homo 5162 REAL(KIND=dp), INTENT(IN) :: eps_filter, omega 5163 INTEGER, INTENT(IN) :: real_imaginary 5164 5165 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Greens_function_freq', & 5166 routineP = moduleN//':'//routineN 5167 5168 INTEGER :: handle, ncol_local 5169 INTEGER, DIMENSION(:), POINTER :: col_indices 5170 REAL(KIND=dp) :: e_fermi 5171 5172 CALL timeset(routineN, handle) 5173 5174 ! get info of fm_mo_coeff 5175 CALL cp_fm_get_info(matrix=fm_mo_coeff, & 5176 ncol_local=ncol_local, & 5177 col_indices=col_indices) 5178 5179 ! we have two different fermi levels for correcting occupied and virtual energy levels 5180 e_fermi = Eigenval(homo) + fermi_level_offset 5181 CALL scale_mo_coeff(fm_mo_coeff_occ_scaled, fm_mo_coeff, e_fermi, ncol_local, & 5182 Eigenval, omega, real_imaginary, col_indices) 5183 5184 e_fermi = Eigenval(homo + 1) - fermi_level_offset 5185 CALL scale_mo_coeff(fm_mo_coeff_virt_scaled, fm_mo_coeff, e_fermi, ncol_local, & 5186 Eigenval, omega, real_imaginary, col_indices) 5187 5188 CALL cp_gemm(transa="N", transb="T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, & 5189 matrix_a=fm_mo_coeff, matrix_b=fm_mo_coeff_occ_scaled, beta=0.0_dp, & 5190 matrix_c=fm_scaled_dm_occ_tau) 5191 5192 CALL cp_gemm(transa="N", transb="T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, & 5193 matrix_a=fm_mo_coeff, matrix_b=fm_mo_coeff_virt_scaled, beta=0.0_dp, & 5194 matrix_c=fm_scaled_dm_virt_tau) 5195 5196 CALL dbcsr_set(mat_greens_fct_occ, 0.0_dp) 5197 5198 CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, & 5199 mat_greens_fct_occ, & 5200 keep_sparsity=.FALSE.) 5201 5202 CALL dbcsr_filter(mat_greens_fct_occ, eps_filter) 5203 5204 CALL dbcsr_set(mat_greens_fct_virt, 0.0_dp) 5205 5206 CALL copy_fm_to_dbcsr(fm_scaled_dm_virt_tau, & 5207 mat_greens_fct_virt, & 5208 keep_sparsity=.FALSE.) 5209 5210 CALL dbcsr_filter(mat_greens_fct_virt, eps_filter) 5211 5212 CALL timestop(handle) 5213 5214 END SUBROUTINE compute_Greens_function_freq 5215 5216! ************************************************************************************************** 5217!> \brief ... 5218!> \param fm_mo_coeff_scaled ... 5219!> \param fm_mo_coeff ... 5220!> \param e_fermi ... 5221!> \param ncol_local ... 5222!> \param Eigenval ... 5223!> \param omega ... 5224!> \param real_imaginary ... 5225!> \param col_indices ... 5226! ************************************************************************************************** 5227 SUBROUTINE scale_mo_coeff(fm_mo_coeff_scaled, fm_mo_coeff, e_fermi, ncol_local, & 5228 Eigenval, omega, real_imaginary, col_indices) 5229 5230 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_scaled, fm_mo_coeff 5231 REAL(KIND=dp) :: e_fermi 5232 INTEGER, INTENT(IN) :: ncol_local 5233 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 5234 REAL(KIND=dp), INTENT(IN) :: omega 5235 INTEGER, INTENT(IN) :: real_imaginary 5236 INTEGER, DIMENSION(:), POINTER :: col_indices 5237 5238 CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_mo_coeff', routineP = moduleN//':'//routineN 5239 5240 INTEGER :: handle, i_global, iiB 5241 REAL(KIND=dp) :: energy_diff 5242 5243 CALL timeset(routineN, handle) 5244 5245 DO iiB = 1, ncol_local 5246 5247 i_global = col_indices(iiB) 5248 5249 energy_diff = e_fermi - Eigenval(i_global) 5250 5251 IF (real_imaginary == 1) THEN 5252 fm_mo_coeff_scaled%local_data(:, iiB) = fm_mo_coeff%local_data(:, iiB)* & 5253 energy_diff/(omega**2 + energy_diff**2) 5254 ELSE IF (real_imaginary == 2) THEN 5255 fm_mo_coeff_scaled%local_data(:, iiB) = fm_mo_coeff%local_data(:, iiB)* & 5256 (-omega/(omega**2 + energy_diff**2)) 5257 END IF 5258 5259 END DO 5260 5261 CALL timestop(handle) 5262 5263 END SUBROUTINE scale_mo_coeff 5264 5265END MODULE rpa_gw 5266 5267