40 dbcsr_type_no_symmetry
67 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
68 dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, &
69 dbt_get_block, dbt_get_info, dbt_iterator_blocks_left, dbt_iterator_next_block, &
70 dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, dbt_nblks_total, &
71 dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
148#include "./base/base_uses.f90"
154 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'rpa_gw'
193 num_integ_points, unit_nr, &
194 RI_blk_sizes, do_ic_model, &
195 para_env, fm_mat_W, fm_mat_Q, &
197 t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
198 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
199 starts_array_mc, ends_array_mc, &
200 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
201 matrix_s, mat_W, t_3c_overl_int, &
202 t_3c_O_compressed, t_3c_O_ind, &
205 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
206 INTEGER,
INTENT(IN) :: nmo, num_integ_points, unit_nr
207 INTEGER,
DIMENSION(:),
POINTER :: ri_blk_sizes
208 LOGICAL,
INTENT(IN) :: do_ic_model
210 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
211 INTENT(OUT) :: fm_mat_w
213 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: mo_coeff
214 TYPE(dbt_type) :: t_3c_overl_int_ao_mo
216 DIMENSION(:) :: t_3c_o_mo_compressed
218 DIMENSION(:),
INTENT(OUT) :: t_3c_o_mo_ind
219 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:), &
220 INTENT(INOUT) :: t_3c_overl_int_gw_ri, &
222 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc
223 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:), &
224 INTENT(INOUT) :: t_3c_overl_nnp_ic, &
225 t_3c_overl_nnp_ic_reflected
228 TYPE(dbt_type),
DIMENSION(:, :) :: t_3c_overl_int
233 CHARACTER(LEN=*),
PARAMETER :: routinen =
'allocate_matrices_gw_im_time'
235 INTEGER :: handle, jquad, nspins
236 LOGICAL :: my_open_shell
237 TYPE(dbt_type) :: t_3c_overl_int_ao_mo_beta
239 CALL timeset(routinen, handle)
242 my_open_shell = (nspins == 2)
244 ALLOCATE (t_3c_o_mo_ind(nspins), t_3c_overl_int_gw_ao(nspins), t_3c_overl_int_gw_ri(nspins), &
245 t_3c_overl_nnp_ic(nspins), t_3c_overl_nnp_ic_reflected(nspins), t_3c_o_mo_compressed(nspins))
247 t_3c_o_compressed, t_3c_o_ind, &
248 t_3c_overl_int_ao_mo, t_3c_o_mo_compressed(1), t_3c_o_mo_ind(1)%array, &
249 t_3c_overl_int_gw_ri(1), t_3c_overl_int_gw_ao(1), &
250 starts_array_mc, ends_array_mc, &
251 mo_coeff(1), matrix_s, &
252 gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), nmo, &
255 t_3c_overl_nnp_ic(1), t_3c_overl_nnp_ic_reflected(1), &
256 qs_env, unit_nr, do_alpha=.true.)
258 IF (my_open_shell)
THEN
261 t_3c_o_compressed, t_3c_o_ind, &
262 t_3c_overl_int_ao_mo_beta, t_3c_o_mo_compressed(2), t_3c_o_mo_ind(2)%array, &
263 t_3c_overl_int_gw_ri(2), t_3c_overl_int_gw_ao(2), &
264 starts_array_mc, ends_array_mc, &
265 mo_coeff(2), matrix_s, &
266 gw_corr_lev_occ(2), gw_corr_lev_virt(2), homo(2), nmo, &
269 t_3c_overl_nnp_ic(2), t_3c_overl_nnp_ic_reflected(2), &
270 qs_env, unit_nr, do_alpha=.false.)
272 IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma)
THEN
273 CALL dbt_destroy(t_3c_overl_int_ao_mo_beta)
278 ALLOCATE (fm_mat_w(num_integ_points))
280 DO jquad = 1, num_integ_points
282 CALL cp_fm_create(fm_mat_w(jquad), fm_mat_q%matrix_struct, set_zero=.true.)
291 template=matrix_s(1)%matrix, &
292 matrix_type=dbcsr_type_no_symmetry, &
293 row_blk_size=ri_blk_sizes, &
294 col_blk_size=ri_blk_sizes)
296 CALL timestop(handle)
340 gw_corr_lev_occ, gw_corr_lev_virt, homo, &
341 nmo, num_integ_group, num_integ_points, unit_nr, &
342 gw_corr_lev_tot, num_fit_points, omega_max_fit, &
343 do_minimax_quad, do_periodic, do_ri_Sigma_x, my_do_gw, &
344 first_cycle_periodic_correction, &
345 a_scaling, Eigenval, tj, vec_omega_fit_gw, vec_Sigma_x_gw, &
346 delta_corr, Eigenval_last, Eigenval_scf, vec_W_gw, &
347 fm_mat_S_gw, fm_mat_S_gw_work, &
348 para_env, mp2_env, kpoints, nkp, nkp_self_energy, &
349 do_kpoints_cubic_RPA, do_kpoints_from_Gamma)
351 COMPLEX(KIND=dp),
ALLOCATABLE, &
352 DIMENSION(:, :, :, :),
INTENT(OUT) :: vec_sigma_c_gw
353 INTEGER,
INTENT(IN) :: color_rpa_group, dimen_nm_gw
354 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
355 INTEGER,
INTENT(IN) :: nmo, num_integ_group, num_integ_points, &
357 INTEGER,
INTENT(INOUT) :: gw_corr_lev_tot, num_fit_points
358 REAL(kind=
dp) :: omega_max_fit
359 LOGICAL,
INTENT(IN) :: do_minimax_quad, do_periodic, &
360 do_ri_sigma_x, my_do_gw
361 LOGICAL,
INTENT(OUT) :: first_cycle_periodic_correction
362 REAL(kind=
dp),
INTENT(IN) :: a_scaling
363 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
364 INTENT(INOUT) :: eigenval
365 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
367 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
368 INTENT(OUT) :: vec_omega_fit_gw
369 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
370 INTENT(OUT) :: vec_sigma_x_gw
371 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
372 INTENT(INOUT) :: delta_corr
373 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
374 INTENT(OUT) :: eigenval_last, eigenval_scf
375 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
376 INTENT(OUT) :: vec_w_gw
377 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_s_gw
378 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
379 INTENT(INOUT) :: fm_mat_s_gw_work
383 INTEGER,
INTENT(OUT) :: nkp, nkp_self_energy
384 LOGICAL,
INTENT(IN) :: do_kpoints_cubic_rpa, &
385 do_kpoints_from_gamma
387 CHARACTER(LEN=*),
PARAMETER :: routinen =
'allocate_matrices_gw'
389 INTEGER :: handle, iquad, ispin, jquad, nspins
390 LOGICAL :: my_open_shell
391 REAL(kind=
dp) :: omega
392 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: vec_omega_gw
394 CALL timeset(routinen, handle)
396 nspins =
SIZE(eigenval, 3)
397 my_open_shell = (nspins == 2)
399 gw_corr_lev_tot = gw_corr_lev_occ(1) + gw_corr_lev_virt(1)
402 ALLOCATE (vec_omega_gw(num_integ_points))
403 vec_omega_gw = 0.0_dp
405 DO jquad = 1, num_integ_points
406 IF (do_minimax_quad)
THEN
409 omega = a_scaling/tan(tj(jquad))
411 vec_omega_gw(jquad) = omega
417 DO jquad = 1, num_integ_points
418 IF (vec_omega_gw(jquad) < omega_max_fit)
THEN
419 num_fit_points = num_fit_points + 1
424 IF (mp2_env%ri_g0w0%nparam_pade > num_fit_points)
THEN
425 IF (unit_nr > 0)
WRITE (unit=unit_nr, fmt=
"(T3,A)") &
426 "Pade approximation: more parameters than data points. Reset # of parameters."
427 mp2_env%ri_g0w0%nparam_pade = num_fit_points
428 IF (unit_nr > 0)
WRITE (unit=unit_nr, fmt=
"(T3,A,T74,I7)") &
429 "Number of pade parameters:", mp2_env%ri_g0w0%nparam_pade
434 ALLOCATE (vec_omega_fit_gw(num_fit_points))
438 DO jquad = 1, num_integ_points
439 IF (vec_omega_gw(jquad) < omega_max_fit)
THEN
441 vec_omega_fit_gw(iquad) = vec_omega_gw(jquad)
445 DEALLOCATE (vec_omega_gw)
447 IF (do_kpoints_cubic_rpa)
THEN
449 IF (mp2_env%ri_g0w0%do_gamma_only_sigma)
THEN
452 nkp_self_energy = nkp
454 ELSE IF (do_kpoints_from_gamma)
THEN
456 IF (mp2_env%ri_g0w0%do_kpoints_Sigma)
THEN
457 nkp_self_energy = mp2_env%ri_g0w0%nkp_self_energy
465 ALLOCATE (vec_sigma_c_gw(gw_corr_lev_tot, num_fit_points, nkp_self_energy, nspins))
468 ALLOCATE (eigenval_scf(nmo, nkp_self_energy, nspins))
469 eigenval_scf(:, :, :) = eigenval(:, :, :)
471 ALLOCATE (eigenval_last(nmo, nkp_self_energy, nspins))
472 eigenval_last(:, :, :) = eigenval(:, :, :)
474 IF (do_periodic)
THEN
476 ALLOCATE (delta_corr(1 + homo(1) - gw_corr_lev_occ(1):homo(1) + gw_corr_lev_virt(1)))
477 delta_corr(:) = 0.0_dp
479 first_cycle_periodic_correction = .true.
483 ALLOCATE (vec_sigma_x_gw(nmo, nkp_self_energy, nspins))
484 vec_sigma_x_gw = 0.0_dp
489 cpassert(.NOT. do_minimax_quad)
492 ALLOCATE (fm_mat_s_gw_work(nspins))
494 CALL cp_fm_create(fm_mat_s_gw_work(ispin), fm_mat_s_gw(ispin)%matrix_struct)
495 CALL cp_fm_set_all(matrix=fm_mat_s_gw_work(ispin), alpha=0.0_dp)
498 ALLOCATE (vec_w_gw(dimen_nm_gw, nspins))
502 IF (do_ri_sigma_x)
THEN
504 CALL get_vec_sigma_x(vec_sigma_x_gw(:, :, 1), nmo, fm_mat_s_gw(1), para_env, num_integ_group, color_rpa_group, &
505 homo(1), gw_corr_lev_occ(1), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, 1))
507 IF (my_open_shell)
THEN
508 CALL get_vec_sigma_x(vec_sigma_x_gw(:, :, 2), nmo, fm_mat_s_gw(2), para_env, num_integ_group, &
509 color_rpa_group, homo(2), gw_corr_lev_occ(2), &
510 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, 1))
517 CALL timestop(handle)
533 SUBROUTINE get_vec_sigma_x(vec_Sigma_x_gw, nmo, fm_mat_S_gw, para_env, num_integ_group, color_rpa_group, homo, &
534 gw_corr_lev_occ, vec_Sigma_x_minus_vxc_gw11)
536 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: vec_sigma_x_gw
537 INTEGER,
INTENT(IN) :: nmo
540 INTEGER,
INTENT(IN) :: num_integ_group, color_rpa_group, homo, &
542 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: vec_sigma_x_minus_vxc_gw11
544 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_vec_sigma_x'
546 INTEGER :: handle, iib, m_global, n_global, &
547 ncol_local, nm_global, nrow_local
548 INTEGER,
DIMENSION(:),
POINTER :: col_indices
550 CALL timeset(routinen, handle)
553 nrow_local=nrow_local, &
554 ncol_local=ncol_local, &
555 col_indices=col_indices)
560 DO iib = 1, ncol_local
563 IF (
modulo(1, num_integ_group) /= color_rpa_group) cycle
565 nm_global = col_indices(iib)
568 n_global = max(1, nm_global - 1)/nmo + 1
569 m_global = nm_global - (n_global - 1)*nmo
570 n_global = n_global + homo - gw_corr_lev_occ
572 IF (m_global <= homo)
THEN
575 vec_sigma_x_gw(n_global, 1) = &
576 vec_sigma_x_gw(n_global, 1) - &
577 dot_product(fm_mat_s_gw%local_data(:, iib), fm_mat_s_gw%local_data(:, iib))
585 CALL para_env%sum(vec_sigma_x_gw)
587 vec_sigma_x_minus_vxc_gw11(:) = &
588 vec_sigma_x_minus_vxc_gw11(:) + &
591 CALL timestop(handle)
593 END SUBROUTINE get_vec_sigma_x
612 vec_Sigma_x_minus_vxc_gw, Eigenval_last, &
613 Eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, kpoints, &
614 vec_Sigma_x_gw, my_do_gw)
616 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
617 INTENT(INOUT) :: fm_mat_s_gw_work
618 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
619 INTENT(INOUT) :: vec_w_gw
620 COMPLEX(KIND=dp),
ALLOCATABLE, &
621 DIMENSION(:, :, :, :),
INTENT(INOUT) :: vec_sigma_c_gw
622 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
623 INTENT(INOUT) :: vec_omega_fit_gw
624 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
625 INTENT(INOUT) :: vec_sigma_x_minus_vxc_gw, eigenval_last, &
627 LOGICAL,
INTENT(IN) :: do_periodic
628 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_re_mo_mo, &
629 matrix_berry_im_mo_mo
631 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
632 INTENT(INOUT) :: vec_sigma_x_gw
633 LOGICAL,
INTENT(IN) :: my_do_gw
635 CHARACTER(LEN=*),
PARAMETER :: routinen =
'deallocate_matrices_gw'
637 INTEGER :: handle, nspins
638 LOGICAL :: my_open_shell
640 CALL timeset(routinen, handle)
642 nspins =
SIZE(eigenval_last, 3)
643 my_open_shell = (nspins == 2)
647 DEALLOCATE (vec_sigma_x_minus_vxc_gw)
648 DEALLOCATE (vec_w_gw)
651 DEALLOCATE (vec_sigma_c_gw)
652 DEALLOCATE (vec_sigma_x_gw)
653 DEALLOCATE (vec_omega_fit_gw)
654 DEALLOCATE (eigenval_last)
655 DEALLOCATE (eigenval_scf)
657 IF (do_periodic)
THEN
663 CALL timestop(handle)
686 t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
687 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
688 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, mat_W, &
691 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
692 INTENT(INOUT) :: weights_cos_tf_w_to_t, &
693 weights_sin_tf_t_to_w
694 LOGICAL,
INTENT(IN) :: do_ic_model, do_kpoints_cubic_rpa
695 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
696 INTENT(INOUT) :: fm_mat_w
697 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_overl_int_ao_mo
699 DIMENSION(:) :: t_3c_o_mo_compressed
701 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:), &
702 INTENT(INOUT) :: t_3c_overl_int_gw_ri, &
703 t_3c_overl_int_gw_ao, &
705 t_3c_overl_nnp_ic_reflected
709 CHARACTER(LEN=*),
PARAMETER :: routinen =
'deallocate_matrices_gw_im_time'
711 INTEGER :: handle, ispin, nspins, unused
712 LOGICAL :: my_open_shell
714 CALL timeset(routinen, handle)
716 nspins =
SIZE(t_3c_overl_int_gw_ri)
717 my_open_shell = (nspins == 2)
719 IF (
ALLOCATED(weights_cos_tf_w_to_t))
DEALLOCATE (weights_cos_tf_w_to_t)
720 IF (
ALLOCATED(weights_sin_tf_t_to_w))
DEALLOCATE (weights_sin_tf_t_to_w)
722 IF (.NOT. do_kpoints_cubic_rpa)
THEN
728 CALL dbt_destroy(t_3c_overl_int_gw_ri(ispin))
729 CALL dbt_destroy(t_3c_overl_int_gw_ao(ispin))
731 DEALLOCATE (t_3c_overl_int_gw_ao, t_3c_overl_int_gw_ri)
732 IF (do_ic_model)
THEN
734 CALL dbt_destroy(t_3c_overl_nnp_ic(ispin))
735 CALL dbt_destroy(t_3c_overl_nnp_ic_reflected(ispin))
737 DEALLOCATE (t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected)
740 IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma)
THEN
742 DEALLOCATE (t_3c_o_mo_ind(ispin)%array)
745 DEALLOCATE (t_3c_o_mo_ind, t_3c_o_mo_compressed)
747 CALL dbt_destroy(t_3c_overl_int_ao_mo)
750 IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma)
THEN
752 CALL dbcsr_release(qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
753 DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
755 CALL dbcsr_release(qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
756 DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
758 DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc)
759 DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks)
762 CALL timestop(handle)
801 gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, &
802 do_im_time, do_periodic, &
803 first_cycle_periodic_correction, fermi_level_offset, &
804 omega, Eigenval, delta_corr, vec_omega_fit_gw, vec_W_gw, wj, &
805 fm_mat_Q, fm_mat_R_gw, fm_mat_S_gw, &
806 fm_mat_S_gw_work, mo_coeff, para_env, &
807 para_env_RPA, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, &
808 kpoints, qs_env, mp2_env)
810 COMPLEX(KIND=dp),
ALLOCATABLE, &
811 DIMENSION(:, :, :, :),
INTENT(INOUT) :: vec_sigma_c_gw
812 INTEGER,
INTENT(IN) :: dimen_nm_gw, dimen_ri
813 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
814 INTEGER,
INTENT(IN) :: jquad, nmo, num_fit_points
815 LOGICAL,
INTENT(IN) :: do_im_time, do_periodic
816 LOGICAL,
INTENT(INOUT) :: first_cycle_periodic_correction
817 REAL(kind=
dp),
INTENT(INOUT) :: fermi_level_offset, omega
818 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: eigenval
819 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
820 INTENT(INOUT) :: delta_corr
821 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
822 INTENT(IN) :: vec_omega_fit_gw
823 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
824 INTENT(INOUT) :: vec_w_gw
825 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
827 TYPE(
cp_fm_type),
INTENT(IN) :: fm_mat_q, fm_mat_r_gw
828 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_s_gw, fm_mat_s_gw_work
831 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_im_mo_mo, &
832 matrix_berry_re_mo_mo
837 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_GW_self_energy'
839 INTEGER :: handle, i_global, iib, ispin, j_global, &
840 jjb, ncol_local, nrow_local, nspins
841 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
843 CALL timeset(routinen, handle)
845 nspins =
SIZE(fm_mat_s_gw)
848 nrow_local=nrow_local, &
849 ncol_local=ncol_local, &
850 row_indices=row_indices, &
851 col_indices=col_indices)
853 IF (.NOT. do_im_time)
THEN
860 IF (do_periodic)
THEN
861 CALL calc_periodic_correction(delta_corr, qs_env, para_env, para_env_rpa, &
862 mp2_env%ri_g0w0%kp_grid, homo(1), nmo, gw_corr_lev_occ(1), &
863 gw_corr_lev_virt(1), omega, mo_coeff, eigenval(:, 1), &
864 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
865 first_cycle_periodic_correction, kpoints, &
866 mp2_env%ri_g0w0%do_mo_coeff_gamma, &
867 mp2_env%ri_g0w0%num_kp_grids, mp2_env%ri_g0w0%eps_kpoint, &
868 mp2_env%ri_g0w0%do_extra_kpoints, &
869 mp2_env%ri_g0w0%do_aux_bas_gw, mp2_env%ri_g0w0%frac_aux_mos)
877 DO jjb = 1, ncol_local
878 j_global = col_indices(jjb)
879 DO iib = 1, nrow_local
880 i_global = row_indices(iib)
881 IF (j_global == i_global .AND. i_global <= dimen_ri)
THEN
882 fm_mat_q%local_data(iib, jjb) = fm_mat_q%local_data(iib, jjb) - 1.0_dp
890 CALL compute_gw_self_energy_deep(vec_sigma_c_gw(:, :, :, ispin), dimen_nm_gw, dimen_ri, &
891 gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), &
892 homo(ispin), jquad, nmo, &
893 num_fit_points, do_periodic, fermi_level_offset, omega, &
894 eigenval(:, ispin), delta_corr, &
895 vec_omega_fit_gw, vec_w_gw(:, ispin), wj, fm_mat_q, &
896 fm_mat_s_gw(ispin), fm_mat_s_gw_work(ispin))
901 CALL timestop(handle)
914 REAL(kind=
dp),
INTENT(INOUT) :: fermi_level_offset
915 REAL(kind=
dp),
INTENT(IN) :: fermi_level_offset_input
916 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: eigenval
917 INTEGER,
DIMENSION(:),
INTENT(IN) :: homo
919 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_fermi_level_offset'
921 INTEGER :: handle, ispin, nspins
923 CALL timeset(routinen, handle)
925 nspins =
SIZE(eigenval, 2)
930 fermi_level_offset = fermi_level_offset_input
932 fermi_level_offset = min(fermi_level_offset, (eigenval(homo(ispin) + 1, ispin) - eigenval(homo(ispin), ispin))*0.5_dp)
935 CALL timestop(handle)
953 SUBROUTINE compute_w_cubic_gw(fm_mat_W, fm_mat_Q, fm_mat_work, dimen_RI, fm_mat_L, num_integ_points, &
954 tj, tau_tj, weights_cos_tf_w_to_t, jquad, omega)
955 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_w
956 TYPE(
cp_fm_type),
INTENT(IN) :: fm_mat_q, fm_mat_work
957 INTEGER,
INTENT(IN) :: dimen_ri
958 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: fm_mat_l
959 INTEGER,
INTENT(IN) :: num_integ_points
960 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
961 INTENT(IN) :: tj, tau_tj
962 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
963 INTENT(IN) :: weights_cos_tf_w_to_t
964 INTEGER,
INTENT(IN) :: jquad
965 REAL(kind=
dp),
INTENT(INOUT) :: omega
967 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_W_cubic_GW'
969 INTEGER :: handle, i_global, iib, iquad, j_global, &
970 jjb, ncol_local, nrow_local
971 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
972 REAL(kind=
dp) :: tau, weight
974 CALL timeset(routinen, handle)
977 nrow_local=nrow_local, &
978 ncol_local=ncol_local, &
979 row_indices=row_indices, &
980 col_indices=col_indices)
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
1001 CALL parallel_gemm(
'T',
'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_mat_l(1, 1), fm_mat_q, &
1002 0.0_dp, fm_mat_work)
1004 CALL parallel_gemm(
'N',
'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_mat_work, fm_mat_l(1, 1), &
1008 DO iquad = 1, num_integ_points
1012 weight = weights_cos_tf_w_to_t(iquad, jquad)*cos(tau*omega)
1014 IF (jquad == 1)
THEN
1020 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_w(iquad), beta=weight, matrix_b=fm_mat_q)
1024 CALL timestop(handle)
1050 SUBROUTINE compute_gw_self_energy_deep(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, &
1051 gw_corr_lev_occ, gw_corr_lev_virt, &
1052 homo, jquad, nmo, num_fit_points, &
1053 do_periodic, fermi_level_offset, omega, Eigenval, &
1054 delta_corr, vec_omega_fit_gw, vec_W_gw, &
1055 wj, fm_mat_Q, fm_mat_S_gw, fm_mat_S_gw_work)
1057 COMPLEX(KIND=dp),
DIMENSION(:, :, :), &
1058 INTENT(INOUT) :: vec_sigma_c_gw
1059 INTEGER,
INTENT(IN) :: dimen_nm_gw, dimen_ri, gw_corr_lev_occ, &
1060 gw_corr_lev_virt, homo, jquad, nmo, &
1062 LOGICAL,
INTENT(IN) :: do_periodic
1063 REAL(kind=
dp),
INTENT(IN) :: fermi_level_offset
1064 REAL(kind=
dp),
INTENT(INOUT) :: omega
1065 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: eigenval
1066 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: delta_corr, vec_omega_fit_gw
1067 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: vec_w_gw
1068 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: wj
1069 TYPE(
cp_fm_type),
INTENT(IN) :: fm_mat_q, fm_mat_s_gw, fm_mat_s_gw_work
1071 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_GW_self_energy_deep'
1073 INTEGER :: handle, iib, iquad, m_global, n_global, &
1074 ncol_local, nm_global
1075 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
1076 REAL(kind=
dp) :: delta_corr_nn, e_fermi, omega_i, &
1079 CALL timeset(routinen, handle)
1082 CALL parallel_gemm(transa=
"N", transb=
"N", m=dimen_ri, n=dimen_nm_gw, k=dimen_ri, alpha=1.0_dp, &
1083 matrix_a=fm_mat_q, matrix_b=fm_mat_s_gw, beta=0.0_dp, &
1084 matrix_c=fm_mat_s_gw_work)
1087 ncol_local=ncol_local, &
1088 row_indices=row_indices, &
1089 col_indices=col_indices)
1095 DO iib = 1, ncol_local
1096 nm_global = col_indices(iib)
1097 vec_w_gw(nm_global) = vec_w_gw(nm_global) + &
1098 dot_product(fm_mat_s_gw_work%local_data(:, iib), fm_mat_s_gw%local_data(:, iib))
1101 n_global = max(1, nm_global - 1)/nmo + 1
1102 m_global = nm_global - (n_global - 1)*nmo
1103 n_global = n_global + homo - gw_corr_lev_occ
1106 DO iquad = 1, num_fit_points
1109 IF (n_global <= homo)
THEN
1110 sign_occ_virt = -1.0_dp
1112 sign_occ_virt = 1.0_dp
1115 omega_i = vec_omega_fit_gw(iquad)*sign_occ_virt
1119 IF (n_global <= homo)
THEN
1120 e_fermi = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset
1122 e_fermi = minval(eigenval(homo + 1:homo + gw_corr_lev_virt)) - fermi_level_offset
1126 IF (do_periodic .AND. row_indices(1) == 1 .AND. n_global == m_global)
THEN
1127 delta_corr_nn = delta_corr(n_global)
1129 delta_corr_nn = 0.0_dp
1135 vec_sigma_c_gw(n_global - homo + gw_corr_lev_occ, iquad, 1) = &
1136 vec_sigma_c_gw(n_global - homo + gw_corr_lev_occ, iquad, 1) - &
1137 0.5_dp/
pi*wj(jquad)/2.0_dp*(vec_w_gw(nm_global) + delta_corr_nn)* &
1138 (1.0_dp/(
gaussi*(omega + omega_i) + e_fermi - eigenval(m_global)) + &
1139 1.0_dp/(
gaussi*(-omega + omega_i) + e_fermi - eigenval(m_global)))
1144 CALL timestop(handle)
1146 END SUBROUTINE compute_gw_self_energy_deep
1215 gw_corr_lev_tot, gw_corr_lev_virt, homo, &
1216 nmo, num_fit_points, num_integ_points, &
1217 unit_nr, do_apply_ic_corr_to_gw, do_im_time, &
1218 do_periodic, do_ri_Sigma_x, &
1219 first_cycle_periodic_correction, e_fermi, eps_filter, &
1220 fermi_level_offset, delta_corr, Eigenval, &
1221 Eigenval_last, Eigenval_scf, iter_sc_GW0, exit_ev_gw, tau_tj, tj, &
1222 vec_omega_fit_gw, vec_Sigma_x_gw, ic_corr_list, &
1223 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, &
1224 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_mo_coeff_occ, &
1225 fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
1226 mo_coeff, fm_mat_W, para_env, para_env_RPA, mat_dm, mat_MinvVMinv, &
1227 t_3c_O, t_3c_M, t_3c_overl_int_ao_mo, &
1228 t_3c_O_compressed, t_3c_O_mo_compressed, &
1229 t_3c_O_ind, t_3c_O_mo_ind, &
1230 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, matrix_berry_im_mo_mo, &
1231 matrix_berry_re_mo_mo, mat_W, matrix_s, &
1232 kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_RPA, &
1233 starts_array_mc, ends_array_mc)
1235 COMPLEX(KIND=dp),
DIMENSION(:, :, :, :), &
1236 INTENT(OUT) :: vec_sigma_c_gw
1237 INTEGER,
INTENT(IN) :: count_ev_sc_gw
1238 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ
1239 INTEGER,
INTENT(IN) :: gw_corr_lev_tot
1240 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_virt, homo
1241 INTEGER,
INTENT(IN) :: nmo, num_fit_points, num_integ_points, &
1243 LOGICAL,
INTENT(IN) :: do_apply_ic_corr_to_gw, do_im_time, &
1244 do_periodic, do_ri_sigma_x
1245 LOGICAL,
INTENT(INOUT) :: first_cycle_periodic_correction
1246 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: e_fermi
1247 REAL(kind=
dp),
INTENT(IN) :: eps_filter, fermi_level_offset
1248 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
1249 INTENT(INOUT) :: delta_corr
1250 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: eigenval
1251 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
1252 INTENT(INOUT) :: eigenval_last, eigenval_scf
1253 INTEGER,
INTENT(IN) :: iter_sc_gw0
1254 LOGICAL,
INTENT(INOUT) :: exit_ev_gw
1255 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
1256 INTENT(INOUT) :: tau_tj, tj, vec_omega_fit_gw
1257 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
1258 INTENT(INOUT) :: vec_sigma_x_gw
1260 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
1261 INTENT(IN) :: weights_cos_tf_t_to_w, &
1262 weights_sin_tf_t_to_w
1263 TYPE(
cp_fm_type),
INTENT(IN) :: fm_mo_coeff_occ_scaled, &
1264 fm_mo_coeff_virt_scaled
1265 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mo_coeff_occ, fm_mo_coeff_virt
1266 TYPE(
cp_fm_type),
INTENT(IN) :: fm_scaled_dm_occ_tau, &
1267 fm_scaled_dm_virt_tau, mo_coeff
1268 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
1269 INTENT(IN) :: fm_mat_w
1271 TYPE(
dbcsr_p_type),
INTENT(IN) :: mat_dm, mat_minvvminv
1272 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_o
1273 TYPE(dbt_type) :: t_3c_m, t_3c_overl_int_ao_mo
1275 DIMENSION(:, :, :),
INTENT(INOUT) :: t_3c_o_compressed
1278 DIMENSION(:, :, :),
INTENT(INOUT) :: t_3c_o_ind
1280 TYPE(dbt_type),
DIMENSION(:) :: t_3c_overl_int_gw_ri, &
1281 t_3c_overl_int_gw_ao
1282 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_im_mo_mo, &
1283 matrix_berry_re_mo_mo
1289 INTEGER,
INTENT(IN) :: nkp_self_energy
1290 LOGICAL,
INTENT(IN) :: do_kpoints_cubic_rpa
1291 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc
1293 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_QP_energies'
1295 INTEGER :: count_ev_sc_gw_print, count_sc_gw0, count_sc_gw0_print, crossing_search, handle, &
1296 idos, ikp, ispin, iunit, n_level_gw, ndos, nspins, num_points_corr, num_poles
1297 LOGICAL :: do_kpoints_sigma, my_open_shell
1298 REAL(kind=
dp) :: dos_lower_bound, dos_precision, dos_upper_bound, e_cbm_gw, e_cbm_gw_beta, &
1299 e_cbm_scf, e_cbm_scf_beta, e_vbm_gw, e_vbm_gw_beta, e_vbm_scf, e_vbm_scf_beta, stop_crit
1300 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: vec_gw_dos
1301 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: m_value, vec_gw_energ, z_value
1305 CALL timeset(routinen, handle)
1308 my_open_shell = (nspins == 2)
1310 do_kpoints_sigma = mp2_env%ri_g0w0%do_kpoints_Sigma
1312 DO count_sc_gw0 = 1, iter_sc_gw0
1315 IF (do_im_time .AND. .NOT. do_kpoints_cubic_rpa .AND. .NOT. do_kpoints_sigma)
THEN
1316 num_points_corr = mp2_env%ri_g0w0%num_omega_points
1318 DO ispin = 1, nspins
1319 CALL compute_self_energy_cubic_gw(num_integ_points, nmo, tau_tj, tj, &
1320 matrix_s, fm_mo_coeff_occ(ispin), &
1321 fm_mo_coeff_virt(ispin), fm_mo_coeff_occ_scaled, &
1322 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, &
1323 fm_scaled_dm_virt_tau, eigenval(:, 1, ispin), eps_filter, &
1324 e_fermi(ispin), fm_mat_w, &
1325 gw_corr_lev_tot, gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), homo(ispin), &
1326 count_ev_sc_gw, count_sc_gw0, &
1327 t_3c_overl_int_ao_mo, t_3c_o_mo_compressed(ispin), &
1328 t_3c_o_mo_ind(ispin)%array, &
1329 t_3c_overl_int_gw_ri(ispin), t_3c_overl_int_gw_ao(ispin), &
1330 mat_w, mat_minvvminv, mat_dm, &
1331 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_sigma_c_gw(:, :, :, ispin), &
1332 do_periodic, num_points_corr, delta_corr, qs_env, para_env, para_env_rpa, &
1333 mp2_env, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
1334 first_cycle_periodic_correction, kpoints, num_fit_points, mo_coeff, &
1335 do_ri_sigma_x, vec_sigma_x_gw(:, :, ispin), unit_nr, ispin)
1340 IF (do_kpoints_sigma)
THEN
1341 CALL compute_self_energy_cubic_gw_kpoints(num_integ_points, tau_tj, tj, &
1342 matrix_s, eigenval(:, :, :), 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_o, t_3c_m, t_3c_o_compressed, t_3c_o_ind, &
1346 mat_w, mat_minvvminv, &
1347 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_sigma_c_gw(:, :, :, :), &
1349 mp2_env, num_fit_points, mo_coeff, &
1350 do_ri_sigma_x, vec_sigma_x_gw(:, :, :), unit_nr, nspins, &
1351 starts_array_mc, ends_array_mc, eps_filter)
1355 IF (do_periodic .AND. mp2_env%ri_g0w0%do_average_deg_levels)
THEN
1357 DO ispin = 1, nspins
1358 CALL average_degenerate_levels(vec_sigma_c_gw(:, :, :, ispin), &
1359 eigenval(1 + homo(ispin) - gw_corr_lev_occ(ispin): &
1360 homo(ispin) + gw_corr_lev_virt(ispin), 1, ispin), &
1361 mp2_env%ri_g0w0%eps_eigenval)
1365 IF (.NOT. do_im_time)
THEN
1366 CALL para_env%sum(vec_sigma_c_gw)
1369 CALL para_env%sync()
1372 num_poles = mp2_env%ri_g0w0%num_poles
1373 crossing_search = mp2_env%ri_g0w0%crossing_search
1376 ALLOCATE (vec_gw_energ(gw_corr_lev_tot, nkp_self_energy, nspins))
1377 vec_gw_energ = 0.0_dp
1378 ALLOCATE (z_value(gw_corr_lev_tot, nkp_self_energy, nspins))
1380 ALLOCATE (m_value(gw_corr_lev_tot, nkp_self_energy, nspins))
1386 e_vbm_gw_beta = -1.0e3
1387 e_cbm_gw_beta = 1.0e3
1388 e_vbm_scf_beta = -1.0e3
1389 e_cbm_scf_beta = 1.0e3
1392 dos_precision = mp2_env%ri_g0w0%dos_prec
1393 dos_upper_bound = mp2_env%ri_g0w0%dos_upper
1394 dos_lower_bound = mp2_env%ri_g0w0%dos_lower
1396 IF (dos_lower_bound >= dos_upper_bound)
THEN
1397 CALL cp_abort(__location__,
"Invalid settings for GW_DOS calculation!")
1400 IF (dos_precision /= 0)
THEN
1401 ndos = int((dos_upper_bound - dos_lower_bound)/dos_precision)
1402 ALLOCATE (vec_gw_dos(ndos))
1407 DO ikp = 1, nkp_self_energy
1409 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
1412 DO n_level_gw = 1, gw_corr_lev_tot
1414 IF (
modulo(n_level_gw, para_env%num_pe) /= para_env%mepos) cycle
1416 SELECT CASE (mp2_env%ri_g0w0%analytic_continuation)
1418 CALL fit_and_continuation_2pole(vec_gw_energ(:, ikp, 1), vec_omega_fit_gw, &
1419 z_value(:, ikp, 1), m_value(:, ikp, 1), vec_sigma_c_gw(:, :, ikp, 1), &
1420 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), &
1421 eigenval(:, ikp, 1), eigenval_scf(:, ikp, 1), n_level_gw, &
1422 gw_corr_lev_occ(1), gw_corr_lev_virt(1), num_poles, &
1423 num_fit_points, crossing_search, homo(1), stop_crit, &
1424 fermi_level_offset, do_im_time)
1428 z_value(:, ikp, 1), m_value(:, ikp, 1), vec_sigma_c_gw(:, :, ikp, 1), &
1429 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), &
1430 eigenval(:, ikp, 1), eigenval_scf(:, ikp, 1), &
1431 mp2_env%ri_g0w0%do_hedin_shift, n_level_gw, &
1432 gw_corr_lev_occ(1), gw_corr_lev_virt(1), mp2_env%ri_g0w0%nparam_pade, &
1433 num_fit_points, crossing_search, homo(1), fermi_level_offset, &
1434 do_im_time, mp2_env%ri_g0w0%print_self_energy, count_ev_sc_gw, &
1435 vec_gw_dos, dos_lower_bound, dos_precision, ndos, &
1436 mp2_env%ri_g0w0%min_level_self_energy, &
1437 mp2_env%ri_g0w0%max_level_self_energy, mp2_env%ri_g0w0%dos_eta, &
1438 mp2_env%ri_g0w0%dos_min, mp2_env%ri_g0w0%dos_max)
1440 cpabort(
"Only two-model and Pade approximation are implemented.")
1443 IF (my_open_shell)
THEN
1444 SELECT CASE (mp2_env%ri_g0w0%analytic_continuation)
1446 CALL fit_and_continuation_2pole( &
1447 vec_gw_energ(:, ikp, 2), vec_omega_fit_gw, &
1448 z_value(:, ikp, 2), m_value(:, ikp, 2), vec_sigma_c_gw(:, :, ikp, 2), &
1449 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), &
1450 eigenval(:, ikp, 2), eigenval_scf(:, ikp, 2), n_level_gw, &
1451 gw_corr_lev_occ(2), gw_corr_lev_virt(2), num_poles, &
1452 num_fit_points, crossing_search, homo(2), stop_crit, &
1453 fermi_level_offset, do_im_time)
1456 z_value(:, ikp, 2), m_value(:, ikp, 2), vec_sigma_c_gw(:, :, ikp, 2), &
1457 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), &
1458 eigenval(:, ikp, 2), eigenval_scf(:, ikp, 2), &
1459 mp2_env%ri_g0w0%do_hedin_shift, n_level_gw, &
1460 gw_corr_lev_occ(2), gw_corr_lev_virt(2), mp2_env%ri_g0w0%nparam_pade, &
1461 num_fit_points, crossing_search, homo(2), &
1462 fermi_level_offset, do_im_time, &
1463 mp2_env%ri_g0w0%print_self_energy, count_ev_sc_gw, &
1464 vec_gw_dos, dos_lower_bound, dos_precision, ndos, &
1465 mp2_env%ri_g0w0%min_level_self_energy, &
1466 mp2_env%ri_g0w0%max_level_self_energy, mp2_env%ri_g0w0%dos_eta, &
1467 mp2_env%ri_g0w0%dos_min, mp2_env%ri_g0w0%dos_max)
1469 cpabort(
"Only two-pole model and Pade approximation are implemented.")
1476 CALL para_env%sum(vec_gw_energ)
1477 CALL para_env%sum(z_value)
1478 CALL para_env%sum(m_value)
1480 IF (dos_precision /= 0.0_dp)
THEN
1481 CALL para_env%sum(vec_gw_dos)
1484 CALL check_nan(vec_gw_energ, 0.0_dp)
1485 CALL check_nan(z_value, 1.0_dp)
1486 CALL check_nan(m_value, 0.0_dp)
1488 IF (do_im_time .OR. mp2_env%ri_g0w0%iter_sc_GW0 == 1)
THEN
1489 count_ev_sc_gw_print = count_ev_sc_gw
1490 count_sc_gw0_print = count_sc_gw0
1492 count_ev_sc_gw_print = count_sc_gw0
1493 count_sc_gw0_print = count_ev_sc_gw
1497 IF (my_open_shell)
THEN
1499 CALL print_and_update_for_ev_sc( &
1500 vec_gw_energ(:, ikp, 1), &
1501 z_value(:, ikp, 1), m_value(:, ikp, 1), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), &
1502 eigenval(:, ikp, 1), eigenval_last(:, ikp, 1), eigenval_scf(:, ikp, 1), &
1503 gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, &
1504 crossing_search, homo(1), unit_nr, count_ev_sc_gw_print, count_sc_gw0_print, &
1505 ikp, nkp_self_energy, kpoints_sigma, 1, e_vbm_gw, e_cbm_gw, e_vbm_scf, e_cbm_scf)
1507 CALL print_and_update_for_ev_sc( &
1508 vec_gw_energ(:, ikp, 2), &
1509 z_value(:, ikp, 2), m_value(:, ikp, 2), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), &
1510 eigenval(:, ikp, 2), eigenval_last(:, ikp, 2), eigenval_scf(:, ikp, 2), &
1511 gw_corr_lev_occ(2), gw_corr_lev_virt(2), gw_corr_lev_tot, &
1512 crossing_search, homo(2), unit_nr, count_ev_sc_gw_print, count_sc_gw0_print, &
1513 ikp, nkp_self_energy, kpoints_sigma, 2, e_vbm_gw_beta, e_cbm_gw_beta, e_vbm_scf_beta, e_cbm_scf_beta)
1515 IF (do_apply_ic_corr_to_gw .AND. count_ev_sc_gw == 1)
THEN
1517 CALL apply_ic_corr(eigenval(:, ikp, 1), eigenval_scf(:, ikp, 1), ic_corr_list(1)%array, &
1518 gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, &
1519 homo(1), nmo, unit_nr, do_alpha=.true.)
1521 CALL apply_ic_corr(eigenval(:, ikp, 2), eigenval_scf(:, ikp, 2), ic_corr_list(2)%array, &
1522 gw_corr_lev_occ(2), gw_corr_lev_virt(2), gw_corr_lev_tot, &
1523 homo(2), nmo, unit_nr, do_beta=.true.)
1529 CALL print_and_update_for_ev_sc( &
1530 vec_gw_energ(:, ikp, 1), &
1531 z_value(:, ikp, 1), m_value(:, ikp, 1), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), &
1532 eigenval(:, ikp, 1), eigenval_last(:, ikp, 1), eigenval_scf(:, ikp, 1), &
1533 gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, &
1534 crossing_search, homo(1), unit_nr, count_ev_sc_gw_print, count_sc_gw0_print, &
1535 ikp, nkp_self_energy, kpoints_sigma, 0, e_vbm_gw, e_cbm_gw, e_vbm_scf, e_cbm_scf)
1537 IF (do_apply_ic_corr_to_gw .AND. count_ev_sc_gw == 1)
THEN
1539 CALL apply_ic_corr(eigenval(:, ikp, 1), eigenval_scf(:, ikp, 1), ic_corr_list(1)%array, &
1540 gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, &
1541 homo(1), nmo, unit_nr)
1549 IF (nkp_self_energy > 1 .AND. unit_nr > 0)
THEN
1551 CALL print_gaps(e_vbm_scf, e_cbm_scf, e_vbm_scf_beta, e_cbm_scf_beta, &
1552 e_vbm_gw, e_cbm_gw, e_vbm_gw_beta, e_cbm_gw_beta, my_open_shell, unit_nr)
1558 IF (mp2_env%ri_g0w0%soc_type /=
soc_none)
THEN
1559 CALL calculate_and_print_soc(qs_env, eigenval_scf, eigenval_scf, gw_corr_lev_occ, gw_corr_lev_virt, &
1560 homo, unit_nr, do_soc_gw=.false., do_soc_scf=.true.)
1561 CALL calculate_and_print_soc(qs_env, eigenval, eigenval_scf, gw_corr_lev_occ, gw_corr_lev_virt, &
1562 homo, unit_nr, do_soc_gw=.true., do_soc_scf=.false.)
1566 IF (logger%para_env%is_source())
THEN
1572 IF (dos_precision /= 0.0_dp)
THEN
1574 CALL open_file(
'spectral.dat', unit_number=iunit, file_status=
"UNKNOWN", file_action=
"WRITE")
1578 WRITE (iunit,
'(E17.10, E17.10)') (dos_lower_bound + real(idos - 1, kind=
dp)*dos_precision)*
evolt, &
1583 DEALLOCATE (vec_gw_dos)
1586 DEALLOCATE (z_value)
1587 DEALLOCATE (m_value)
1588 DEALLOCATE (vec_gw_energ)
1590 exit_ev_gw = .false.
1593 IF (abs(eigenval(homo(1), 1, 1) - eigenval_last(homo(1), 1, 1) - &
1594 eigenval(homo(1) + 1, 1, 1) + eigenval_last(homo(1) + 1, 1, 1)) &
1595 < mp2_env%ri_g0w0%eps_iter)
THEN
1596 IF (count_sc_gw0 == 1) exit_ev_gw = .true.
1600 DO ispin = 1, nspins
1601 CALL shift_unshifted_levels(eigenval(:, 1, ispin), eigenval_last(:, 1, ispin), gw_corr_lev_occ(ispin), &
1602 gw_corr_lev_virt(ispin), homo(ispin), nmo)
1605 IF (do_im_time .AND. do_kpoints_sigma .AND. mp2_env%ri_g0w0%print_local_bandgap)
THEN
1606 CALL print_local_bandgap(qs_env, eigenval, gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1),
"GW")
1607 CALL print_local_bandgap(qs_env, eigenval_scf, gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1),
"DFT")
1611 IF (.NOT. do_im_time)
EXIT
1615 CALL timestop(handle)
1631 SUBROUTINE calculate_and_print_soc(qs_env, Eigenval, Eigenval_scf, gw_corr_lev_occ, gw_corr_lev_virt, &
1632 homo, unit_nr, do_soc_gw, do_soc_scf)
1634 REAL(kind=
dp),
DIMENSION(:, :, :) :: eigenval, eigenval_scf
1635 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
1637 LOGICAL :: do_soc_gw, do_soc_scf
1639 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_and_print_soc'
1641 INTEGER :: handle, i_dim, i_glob, i_row, ikp, j_col, j_glob, n_level_gw, nao, ncol_local, &
1642 nder, nkind, nkp_self_energy, nrow_local, periodic(3), size_real_space
1643 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: index0
1644 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
1645 LOGICAL :: calculate_forces, use_virial
1646 REAL(kind=
dp) :: avg_occ_qp_shift, avg_virt_qp_shift, e_cbm_gw_soc, e_gap_gw_soc, e_homo, &
1647 e_homo_gw_soc, e_i, e_j, e_lumo, e_lumo_gw_soc, e_vbm_gw_soc, e_window, eps_ppnl
1648 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eigenvalues_without_soc_sorted
1649 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eigenvalues
1652 TYPE(
cp_cfm_type) :: cfm_mat_h_double, cfm_mat_h_ks, &
1653 cfm_mat_s_double, cfm_mat_work_double, &
1654 cfm_mo_coeff, cfm_mo_coeff_double
1656 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, matrix_s_desymm
1657 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_vsoc_l_nosymm, mat_vsoc_lx_kp, &
1658 mat_vsoc_ly_kp, mat_vsoc_lz_kp, &
1659 matrix_dummy, matrix_l, &
1665 POINTER :: sab_orb, sap_ppnl
1668 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1672 CALL timeset(routinen, handle)
1674 cpassert(do_soc_gw .NEQV. do_soc_scf)
1677 matrix_s=matrix_s, &
1678 para_env=para_env, &
1679 qs_kind_set=qs_kind_set, &
1681 atomic_kind_set=atomic_kind_set, &
1682 particle_set=particle_set, &
1683 sap_ppnl=sap_ppnl, &
1684 dft_control=dft_control, &
1687 scf_control=scf_control)
1689 calculate_forces = .false.
1690 use_virial = .false.
1692 eps_ppnl = dft_control%qs_control%eps_ppnl
1694 CALL get_cell(cell=cell, periodic=periodic)
1696 size_real_space = 3**(periodic(1) + periodic(2) + periodic(3))
1701 ALLOCATE (matrix_l(i_dim, 1)%matrix)
1702 CALL dbcsr_create(matrix_l(i_dim, 1)%matrix, template=matrix_s(1)%matrix, &
1703 matrix_type=dbcsr_type_antisymmetric)
1705 CALL dbcsr_set(matrix_l(i_dim, 1)%matrix, 0.0_dp)
1708 NULLIFY (matrix_pot_dummy)
1710 ALLOCATE (matrix_pot_dummy(1, 1)%matrix)
1711 CALL dbcsr_create(matrix_pot_dummy(1, 1)%matrix, template=matrix_s(1)%matrix)
1713 CALL dbcsr_set(matrix_pot_dummy(1, 1)%matrix, 0.0_dp)
1715 CALL build_core_ppnl(matrix_pot_dummy, matrix_dummy, force, virial, calculate_forces, use_virial, nder, &
1716 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
1717 nimages=1, basis_type=
"ORB", matrix_l=matrix_l)
1719 CALL alloc_mat_set_2d(mat_vsoc_l_nosymm, 3, size_real_space, matrix_s(1)%matrix, explicitly_no_symmetry=.true.)
1721 CALL dbcsr_desymmetrize(matrix_l(i_dim, 1)%matrix, mat_vsoc_l_nosymm(i_dim, 1)%matrix)
1724 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
1726 CALL mat_kp_from_mat_gamma(qs_env, mat_vsoc_lx_kp, mat_vsoc_l_nosymm(1, 1)%matrix, kpoints_sigma, 1, .false.)
1727 CALL mat_kp_from_mat_gamma(qs_env, mat_vsoc_ly_kp, mat_vsoc_l_nosymm(2, 1)%matrix, kpoints_sigma, 1, .false.)
1728 CALL mat_kp_from_mat_gamma(qs_env, mat_vsoc_lz_kp, mat_vsoc_l_nosymm(3, 1)%matrix, kpoints_sigma, 1, .false.)
1730 nkp_self_energy = kpoints_sigma%nkp
1732 CALL get_mo_set(kpoints_sigma%kp_env(1)%kpoint_env%mos(1, 1), mo_coeff=rmos)
1734 CALL create_cfm_double_row_col_size(rmos, cfm_mat_h_double)
1735 CALL create_cfm_double_row_col_size(rmos, cfm_mat_s_double)
1736 CALL create_cfm_double_row_col_size(rmos, cfm_mo_coeff_double)
1737 CALL create_cfm_double_row_col_size(rmos, cfm_mat_work_double)
1746 NULLIFY (matrix_s_desymm)
1748 ALLOCATE (matrix_s_desymm(1)%matrix)
1749 CALL dbcsr_create(matrix=matrix_s_desymm(1)%matrix, template=matrix_s(1)%matrix, &
1750 matrix_type=dbcsr_type_no_symmetry)
1753 ALLOCATE (eigenvalues(2*nao))
1754 eigenvalues = 0.0_dp
1755 ALLOCATE (eigenvalues_without_soc_sorted(2*nao))
1757 e_window = qs_env%mp2_env%ri_g0w0%soc_energy_window
1758 IF (unit_nr > 0)
THEN
1759 WRITE (unit_nr,
'(T3,A)')
' '
1760 WRITE (unit_nr,
'(T3,A)')
'------------------------------------------------------------------------------'
1761 WRITE (unit_nr,
'(T3,A)')
' '
1762 WRITE (unit_nr,
'(T3,A,F42.1)')
'GW_SOC_INFO | SOC energy window (eV)', e_window*
evolt
1765 e_vbm_gw_soc = -1000.0_dp
1766 e_cbm_gw_soc = 1000.0_dp
1768 DO ikp = 1, nkp_self_energy
1770 CALL get_mo_set(kpoints_sigma%kp_env(ikp)%kpoint_env%mos(1, 1), mo_coeff=rmos)
1771 CALL get_mo_set(kpoints_sigma%kp_env(ikp)%kpoint_env%mos(2, 1), mo_coeff=imos)
1775 avg_occ_qp_shift = sum(eigenval(homo(1) - gw_corr_lev_occ(1) + 1:homo(1), ikp, 1) - &
1776 eigenval_scf(homo(1) - gw_corr_lev_occ(1) + 1:homo(1), ikp, 1))/gw_corr_lev_occ(1)
1777 avg_virt_qp_shift = sum(eigenval(homo(1):homo(1) + gw_corr_lev_virt(1), ikp, 1) - &
1778 eigenval_scf(homo(1):homo(1) + gw_corr_lev_virt(1), ikp, 1))/gw_corr_lev_virt(1)
1780 IF (gw_corr_lev_occ(1) < homo(1))
THEN
1781 eigenval(1:homo(1) - gw_corr_lev_occ(1), ikp, 1) = eigenval_scf(1:homo(1) - gw_corr_lev_occ(1), ikp, 1) &
1784 IF (gw_corr_lev_virt(1) < nao - homo(1) + 1)
THEN
1785 eigenval(homo(1) + gw_corr_lev_virt(1) + 1:nao, ikp, 1) = eigenval_scf(homo(1) + gw_corr_lev_virt(1) + 1:nao, ikp, 1) &
1790 CALL add_dbcsr_submatrix(cfm_mat_h_double, mat_vsoc_lx_kp(ikp, 1:2), cfm_mat_h_ks, nao + 1, 1,
z_one, .true.)
1791 CALL add_dbcsr_submatrix(cfm_mat_h_double, mat_vsoc_ly_kp(ikp, 1:2), cfm_mat_h_ks, nao + 1, 1,
gaussi, .true.)
1792 CALL add_dbcsr_submatrix(cfm_mat_h_double, mat_vsoc_lz_kp(ikp, 1:2), cfm_mat_h_ks, 1, 1,
z_one, .false.)
1793 CALL add_dbcsr_submatrix(cfm_mat_h_double, mat_vsoc_lz_kp(ikp, 1:2), cfm_mat_h_ks, nao + 1, nao + 1, -
z_one, .false.)
1796 cfm_mo_coeff_double%local_data =
z_zero
1797 CALL add_cfm_submatrix(cfm_mo_coeff_double, cfm_mo_coeff, 1, 1)
1798 CALL add_cfm_submatrix(cfm_mo_coeff_double, cfm_mo_coeff, nao + 1, nao + 1)
1801 nrow_local=nrow_local, &
1802 ncol_local=ncol_local, &
1803 row_indices=row_indices, &
1804 col_indices=col_indices)
1807 matrix_a=cfm_mat_h_double, matrix_b=cfm_mo_coeff_double, beta=
z_zero, &
1808 matrix_c=cfm_mat_work_double)
1811 matrix_a=cfm_mo_coeff_double, matrix_b=cfm_mat_work_double, beta=
z_zero, &
1812 matrix_c=cfm_mat_h_double)
1815 nrow_local=nrow_local, &
1816 ncol_local=ncol_local, &
1817 row_indices=row_indices, &
1818 col_indices=col_indices)
1822 e_homo = eigenval(homo(1), ikp, 1)
1823 e_lumo = eigenval(homo(1) + 1, ikp, 1)
1825 CALL para_env%sync()
1827 DO i_row = 1, nrow_local
1828 DO j_col = 1, ncol_local
1829 i_glob = row_indices(i_row)
1830 j_glob = col_indices(j_col)
1831 IF (i_glob <= nao)
THEN
1832 e_i = eigenval(i_glob, ikp, 1)
1834 e_i = eigenval(i_glob - nao, ikp, 1)
1836 IF (j_glob <= nao)
THEN
1837 e_j = eigenval(j_glob, ikp, 1)
1839 e_j = eigenval(j_glob - nao, ikp, 1)
1843 IF (i_glob == j_glob)
THEN
1844 cfm_mat_h_double%local_data(i_row, j_col) = cfm_mat_h_double%local_data(i_row, j_col) + e_i*
z_one
1845 cfm_mat_s_double%local_data(i_row, j_col) =
z_one
1847 IF (e_i < e_homo - 0.5_dp*e_window .OR. e_i > e_lumo + 0.5_dp*e_window .OR. &
1848 e_j < e_homo - 0.5_dp*e_window .OR. e_j > e_lumo + 0.5_dp*e_window)
THEN
1849 cfm_mat_h_double%local_data(i_row, j_col) =
z_zero
1856 CALL para_env%sync()
1858 eigenvalues = 0.0_dp
1859 CALL cp_cfm_geeig_canon(cfm_mat_h_double, cfm_mat_s_double, cfm_mo_coeff_double, eigenvalues, &
1860 cfm_mat_work_double, scf_control%eps_eigval)
1862 eigenvalues_without_soc_sorted(1:nao) = eigenval(:, ikp, 1)
1863 eigenvalues_without_soc_sorted(nao + 1:2*nao) = eigenval(:, ikp, 1)
1864 ALLOCATE (index0(2*nao))
1865 CALL sort(eigenvalues_without_soc_sorted, 2*nao, index0)
1868 e_homo_gw_soc = maxval(eigenvalues(2*homo(1) - 2*gw_corr_lev_occ(1) + 1:2*homo(1)))
1869 e_lumo_gw_soc = minval(eigenvalues(2*homo(1) + 1:2*homo(1) + 2*gw_corr_lev_virt(1)))
1870 e_gap_gw_soc = e_lumo_gw_soc - e_homo_gw_soc
1871 IF (e_homo_gw_soc > e_vbm_gw_soc) e_vbm_gw_soc = e_homo_gw_soc
1872 IF (e_lumo_gw_soc < e_cbm_gw_soc) e_cbm_gw_soc = e_lumo_gw_soc
1874 IF (unit_nr > 0)
THEN
1875 WRITE (unit_nr,
'(T3,A)')
' '
1876 WRITE (unit_nr,
'(T3,A7,I3,A3,I3,A8,3F7.3,A12,3F7.3)')
'Kpoint ', ikp,
' /', nkp_self_energy, &
1877 ' xkp =', kpoints_sigma%xkp(1, ikp), kpoints_sigma%xkp(2, ikp), kpoints_sigma%xkp(3, ikp), &
1878 ' and xkp =', -kpoints_sigma%xkp(1, ikp), -kpoints_sigma%xkp(2, ikp), -kpoints_sigma%xkp(3, ikp)
1879 WRITE (unit_nr,
'(T3,A)')
' '
1881 WRITE (unit_nr,
'(T3,A)')
' '
1882 WRITE (unit_nr,
'(T3,A,F13.4)')
'GW_SOC_INFO | Average GW shift of occupied levels compared to SCF', &
1883 avg_occ_qp_shift*
evolt
1884 WRITE (unit_nr,
'(T3,A,F11.4)')
'GW_SOC_INFO | Average GW shift of unoccupied levels compared to SCF', &
1885 avg_virt_qp_shift*
evolt
1886 WRITE (unit_nr,
'(T3,A)')
' '
1887 WRITE (unit_nr,
'(T3,2A)')
'Molecular orbital E_GW with SOC (eV) E_GW without SOC (eV) SOC shift (eV)'
1889 WRITE (unit_nr,
'(T3,2A)')
'Molecular orbital E_SCF with SOC (eV) E_SCF without SOC (eV) SOC shift (eV)'
1892 DO n_level_gw = 2*(homo(1) - gw_corr_lev_occ(1)) + 1, 2*homo(1)
1893 WRITE (unit_nr,
'(T3,I4,A,3F21.4)') n_level_gw,
' ( occ ) ', eigenvalues(n_level_gw)*
evolt, &
1894 eigenvalues_without_soc_sorted(n_level_gw)*
evolt, &
1895 (eigenvalues(n_level_gw) - eigenvalues_without_soc_sorted(n_level_gw))*
evolt
1897 DO n_level_gw = 2*homo(1) + 1, 2*(homo(1) + gw_corr_lev_virt(1))
1898 WRITE (unit_nr,
'(T3,I4,A,3F21.4)') n_level_gw,
' ( vir ) ', eigenvalues(n_level_gw)*
evolt, &
1899 eigenvalues_without_soc_sorted(n_level_gw)*
evolt, &
1900 (eigenvalues(n_level_gw) - eigenvalues_without_soc_sorted(n_level_gw))*
evolt
1902 WRITE (unit_nr,
'(T3,A)')
' '
1904 WRITE (unit_nr,
'(T3,A,F38.4)')
'GW+SOC direct gap at current kpoint (eV)', e_gap_gw_soc*
evolt
1906 WRITE (unit_nr,
'(T3,A,F37.4)')
'SCF+SOC direct gap at current kpoint (eV)', e_gap_gw_soc*
evolt
1908 WRITE (unit_nr,
'(T3,A)')
' '
1909 WRITE (unit_nr,
'(T3,A)')
'------------------------------------------------------------------------------'
1914 IF (unit_nr > 0)
THEN
1915 WRITE (unit_nr,
'(T3,A)')
' '
1917 WRITE (unit_nr,
'(T3,A,F46.4)')
'GW+SOC valence band maximum (eV)', e_vbm_gw_soc*
evolt
1918 WRITE (unit_nr,
'(T3,A,F43.4)')
'GW+SOC conduction band minimum (eV)', e_cbm_gw_soc*
evolt
1919 WRITE (unit_nr,
'(T3,A,F59.4)')
'GW+SOC bandgap (eV)', (e_cbm_gw_soc - e_vbm_gw_soc)*
evolt
1921 WRITE (unit_nr,
'(T3,A,F45.4)')
'SCF+SOC valence band maximum (eV)', e_vbm_gw_soc*
evolt
1922 WRITE (unit_nr,
'(T3,A,F42.4)')
'SCF+SOC conduction band minimum (eV)', e_cbm_gw_soc*
evolt
1923 WRITE (unit_nr,
'(T3,A,F58.4)')
'SCF+SOC bandgap (eV)', (e_cbm_gw_soc - e_vbm_gw_soc)*
evolt
1941 DEALLOCATE (eigenvalues)
1943 CALL timestop(handle)
1945 END SUBROUTINE calculate_and_print_soc
1957 SUBROUTINE add_dbcsr_submatrix(cfm_mat_target, mat_source, cfm_source_template, &
1958 nstart_row, nstart_col, factor, add_also_herm_conj)
1962 INTEGER :: nstart_row, nstart_col
1963 COMPLEX(KIND=dp) :: factor
1964 LOGICAL :: add_also_herm_conj
1966 CHARACTER(LEN=*),
PARAMETER :: routinen =
'add_dbcsr_submatrix'
1968 INTEGER :: handle, nao
1970 cfm_mat_work_double_2
1972 fm_mat_work_double_re, fm_mat_work_im, &
1975 CALL timeset(routinen, handle)
1977 CALL cp_fm_create(fm_mat_work_double_re, cfm_mat_target%matrix_struct)
1978 CALL cp_fm_create(fm_mat_work_double_im, cfm_mat_target%matrix_struct)
1982 CALL cp_cfm_create(cfm_mat_work_double, cfm_mat_target%matrix_struct)
1983 CALL cp_cfm_create(cfm_mat_work_double_2, cfm_mat_target%matrix_struct)
1987 CALL cp_fm_create(fm_mat_work_re, cfm_source_template%matrix_struct)
1988 CALL cp_fm_create(fm_mat_work_im, cfm_source_template%matrix_struct)
1996 nrow=nao, ncol=nao, &
1997 s_firstrow=1, s_firstcol=1, &
1998 t_firstrow=nstart_row, t_firstcol=nstart_col)
2001 nrow=nao, ncol=nao, &
2002 s_firstrow=1, s_firstcol=1, &
2003 t_firstrow=nstart_row, t_firstcol=nstart_col)
2012 IF (add_also_herm_conj)
THEN
2024 CALL timestop(handle)
2026 END SUBROUTINE add_dbcsr_submatrix
2035 SUBROUTINE add_cfm_submatrix(cfm_mat_target, cfm_mat_source, nstart_row, nstart_col)
2037 TYPE(
cp_cfm_type) :: cfm_mat_target, cfm_mat_source
2038 INTEGER :: nstart_row, nstart_col
2040 CHARACTER(LEN=*),
PARAMETER :: routinen =
'add_cfm_submatrix'
2042 INTEGER :: handle, nao
2044 fm_mat_work_double_re, fm_mat_work_im, &
2047 CALL timeset(routinen, handle)
2049 CALL cp_fm_create(fm_mat_work_double_re, cfm_mat_target%matrix_struct)
2050 CALL cp_fm_create(fm_mat_work_double_im, cfm_mat_target%matrix_struct)
2054 CALL cp_fm_create(fm_mat_work_re, cfm_mat_source%matrix_struct)
2055 CALL cp_fm_create(fm_mat_work_im, cfm_mat_source%matrix_struct)
2056 CALL cp_cfm_to_fm(cfm_mat_source, fm_mat_work_re, fm_mat_work_im)
2061 nrow=nao, ncol=nao, &
2062 s_firstrow=1, s_firstcol=1, &
2063 t_firstrow=nstart_row, t_firstcol=nstart_col)
2066 nrow=nao, ncol=nao, &
2067 s_firstrow=1, s_firstcol=1, &
2068 t_firstrow=nstart_row, t_firstcol=nstart_col)
2078 CALL timestop(handle)
2080 END SUBROUTINE add_cfm_submatrix
2087 SUBROUTINE create_cfm_double_row_col_size(fm_orig, cfm_double)
2091 CHARACTER(LEN=*),
PARAMETER :: routinen =
'create_cfm_double_row_col_size'
2093 INTEGER :: handle, ncol_global_orig, &
2097 CALL timeset(routinen, handle)
2099 CALL cp_fm_get_info(matrix=fm_orig, nrow_global=nrow_global_orig, ncol_global=ncol_global_orig)
2102 nrow_global=2*nrow_global_orig, &
2103 ncol_global=2*ncol_global_orig, &
2104 template_fmstruct=fm_orig%matrix_struct)
2110 CALL timestop(handle)
2112 END SUBROUTINE create_cfm_double_row_col_size
2127 SUBROUTINE print_gaps(E_VBM_SCF, E_CBM_SCF, E_VBM_SCF_beta, E_CBM_SCF_beta, &
2128 E_VBM_GW, E_CBM_GW, E_VBM_GW_beta, E_CBM_GW_beta, my_open_shell, unit_nr)
2130 REAL(kind=
dp) :: e_vbm_scf, e_cbm_scf, e_vbm_scf_beta, &
2131 e_cbm_scf_beta, e_vbm_gw, e_cbm_gw, &
2132 e_vbm_gw_beta, e_cbm_gw_beta
2133 LOGICAL :: my_open_shell
2136 IF (my_open_shell)
THEN
2137 WRITE (unit_nr,
'(T3,A)')
' '
2138 WRITE (unit_nr,
'(T3,A,F43.4)')
'Alpha SCF valence band maximum (eV)', e_vbm_scf*
evolt
2139 WRITE (unit_nr,
'(T3,A,F40.4)')
'Alpha SCF conduction band minimum (eV)', e_cbm_scf*
evolt
2140 WRITE (unit_nr,
'(T3,A,F56.4)')
'Alpha SCF bandgap (eV)', (e_cbm_scf - e_vbm_scf)*
evolt
2141 WRITE (unit_nr,
'(T3,A)')
' '
2142 WRITE (unit_nr,
'(T3,A,F44.4)')
'Beta SCF valence band maximum (eV)', e_vbm_scf_beta*
evolt
2143 WRITE (unit_nr,
'(T3,A,F41.4)')
'Beta SCF conduction band minimum (eV)', e_cbm_scf_beta*
evolt
2144 WRITE (unit_nr,
'(T3,A,F57.4)')
'Beta SCF bandgap (eV)', (e_cbm_scf_beta - e_vbm_scf_beta)*
evolt
2145 WRITE (unit_nr,
'(T3,A)')
' '
2146 WRITE (unit_nr,
'(T3,A,F44.4)')
'Alpha GW valence band maximum (eV)', e_vbm_gw*
evolt
2147 WRITE (unit_nr,
'(T3,A,F41.4)')
'Alpha GW conduction band minimum (eV)', e_cbm_gw*
evolt
2148 WRITE (unit_nr,
'(T3,A,F57.4)')
'Alpha GW bandgap (eV)', (e_cbm_gw - e_vbm_gw)*
evolt
2149 WRITE (unit_nr,
'(T3,A)')
' '
2150 WRITE (unit_nr,
'(T3,A,F45.4)')
'Beta GW valence band maximum (eV)', e_vbm_gw_beta*
evolt
2151 WRITE (unit_nr,
'(T3,A,F42.4)')
'Beta GW conduction band minimum (eV)', e_cbm_gw_beta*
evolt
2152 WRITE (unit_nr,
'(T3,A,F58.4)')
'Beta GW bandgap (eV)', (e_cbm_gw_beta - e_vbm_gw_beta)*
evolt
2154 WRITE (unit_nr,
'(T3,A)')
' '
2155 WRITE (unit_nr,
'(T3,A,F49.4)')
'SCF valence band maximum (eV)', e_vbm_scf*
evolt
2156 WRITE (unit_nr,
'(T3,A,F46.4)')
'SCF conduction band minimum (eV)', e_cbm_scf*
evolt
2157 WRITE (unit_nr,
'(T3,A,F62.4)')
'SCF bandgap (eV)', (e_cbm_scf - e_vbm_scf)*
evolt
2158 WRITE (unit_nr,
'(T3,A)')
' '
2159 WRITE (unit_nr,
'(T3,A,F50.4)')
'GW valence band maximum (eV)', e_vbm_gw*
evolt
2160 WRITE (unit_nr,
'(T3,A,F47.4)')
'GW conduction band minimum (eV)', e_cbm_gw*
evolt
2161 WRITE (unit_nr,
'(T3,A,F63.4)')
'GW bandgap (eV)', (e_cbm_gw - e_vbm_gw)*
evolt
2164 END SUBROUTINE print_gaps
2171 SUBROUTINE check_nan(array, real_value)
2172 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
2173 INTENT(INOUT) :: array
2174 REAL(kind=
dp),
INTENT(IN) :: real_value
2176 CHARACTER(LEN=*),
PARAMETER :: routinen =
'check_NaN'
2178 INTEGER :: handle, i, j, k
2180 CALL timeset(routinen, handle)
2182 DO i = 1,
SIZE(array, 1)
2183 DO j = 1,
SIZE(array, 2)
2184 DO k = 1,
SIZE(array, 3)
2187 IF (array(i, j, k) /= array(i, j, k)) array(i, j, k) = real_value
2193 CALL timestop(handle)
2195 END SUBROUTINE check_nan
2206 SUBROUTINE print_local_bandgap(qs_env, Eigenval, gw_corr_lev_occ, gw_corr_lev_virt, homo, dft_gw_char)
2208 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(IN) :: eigenval
2209 INTEGER :: gw_corr_lev_occ, gw_corr_lev_virt, homo
2210 CHARACTER(len=*) :: dft_gw_char
2212 CHARACTER(LEN=*),
PARAMETER :: routinen =
'print_local_bandgap'
2214 INTEGER :: handle, i_e
2220 CALL timeset(routinen, handle)
2222 CALL create_real_space_grids(e_gap_rspace, e_vbm_rspace, e_cbm_rspace, rho_g_dummy, ldos, auxbas_pw_pool, qs_env)
2224 CALL calculate_e_gap_rspace(e_gap_rspace, e_vbm_rspace, e_cbm_rspace, rho_g_dummy, &
2225 ldos, qs_env, eigenval, gw_corr_lev_occ, gw_corr_lev_virt, homo, dft_gw_char)
2227 CALL auxbas_pw_pool%give_back_pw(e_gap_rspace)
2228 CALL auxbas_pw_pool%give_back_pw(e_vbm_rspace)
2229 CALL auxbas_pw_pool%give_back_pw(e_cbm_rspace)
2230 CALL auxbas_pw_pool%give_back_pw(rho_g_dummy)
2231 DO i_e = 1,
SIZE(ldos)
2232 CALL auxbas_pw_pool%give_back_pw(ldos(i_e))
2236 CALL timestop(handle)
2238 END SUBROUTINE print_local_bandgap
2254 SUBROUTINE calculate_e_gap_rspace(E_gap_rspace, E_VBM_rspace, E_CBM_rspace, rho_g_dummy, &
2255 LDOS, qs_env, Eigenval, gw_corr_lev_occ, gw_corr_lev_virt, homo, dft_gw_char)
2260 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(IN) :: eigenval
2261 INTEGER :: gw_corr_lev_occ, gw_corr_lev_virt, homo
2262 CHARACTER(len=*) :: dft_gw_char
2264 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_E_gap_rspace'
2266 INTEGER :: handle, i_e, i_img, i_spin, i_x, i_y, i_z, ikp, imo, n_e, n_e_occ, n_x_end, &
2267 n_x_start, n_y_end, n_y_start, n_z_end, n_z_start, nimg, nkp, nkp_self_energy
2268 REAL(kind=
dp) :: avg_ldos_occ, avg_ldos_virt, d_e, e_cbm, &
2269 e_cbm_at_k, e_diff, e_vbm, e_vbm_at_k
2270 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: e_array
2271 REAL(kind=
dp),
DIMENSION(:),
POINTER :: occupation
2273 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: fm_work
2274 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, rho_ao
2275 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rho_ao_weighted
2288 CALL timeset(routinen, handle)
2290 CALL get_qs_env(qs_env=qs_env, para_env=para_env, mp2_env=mp2_env, ks_env=ks_env, matrix_s=matrix_s, &
2291 scf_env=scf_env, sab_orb=sab_orb, dft_control=dft_control, subsys=subsys)
2294 nkp =
SIZE(eigenval, 2)
2300 e_vbm_at_k = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo, ikp, 1))
2301 IF (e_vbm_at_k > e_vbm) e_vbm = e_vbm_at_k
2303 e_cbm_at_k = minval(eigenval(homo + 1:homo + gw_corr_lev_virt, ikp, 1))
2304 IF (e_cbm_at_k < e_cbm) e_cbm = e_cbm_at_k
2308 d_e = mp2_env%ri_g0w0%energy_spacing_print_loc_bandgap
2310 n_e = int(mp2_env%ri_g0w0%energy_window_print_loc_bandgap/d_e)
2313 ALLOCATE (e_array(n_e))
2315 e_array(i_e) = e_vbm - real(n_e_occ - i_e, kind=
dp)*d_e
2317 DO i_e = n_e_occ + 1, n_e
2318 e_array(i_e) = e_cbm + real(i_e - n_e_occ - 1, kind=
dp)*d_e
2321 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
2323 nkp_self_energy = kpoints_sigma%nkp
2324 cpassert(nkp == nkp_self_energy)
2326 kpoints_sigma%sab_nl => sab_orb
2328 DEALLOCATE (kpoints_sigma%cell_to_index)
2329 NULLIFY (kpoints_sigma%cell_to_index)
2332 nimg = maxval(kpoints_sigma%cell_to_index)
2334 NULLIFY (rho_ao_weighted)
2339 ALLOCATE (rho_ao_weighted(i_spin, i_img)%matrix)
2340 CALL dbcsr_create(matrix=rho_ao_weighted(i_spin, i_img)%matrix, template=matrix_s(1)%matrix)
2342 CALL dbcsr_set(rho_ao_weighted(i_spin, i_img)%matrix, 0.0_dp)
2346 ALLOCATE (fm_work(nimg))
2347 matrix_struct => kpoints_sigma%kp_env(1)%kpoint_env%mos(1, 1)%mo_coeff%matrix_struct
2356 CALL get_mo_set(kpoints_sigma%kp_env(ikp)%kpoint_env%mos(1, 1), &
2357 occupation_numbers=occupation)
2359 occupation(:) = 0.0_dp
2360 DO imo = homo - gw_corr_lev_occ + 1, homo + gw_corr_lev_virt
2361 e_diff = e_array(i_e) - eigenval(imo, ikp, 1)
2362 occupation(imo) = exp(-(e_diff/d_e)**2)
2367 CALL get_mo_set(kpoints_sigma%kp_env(1)%kpoint_env%mos(1, 1), &
2368 occupation_numbers=occupation)
2375 matrix_s(1)%matrix, sab_orb, fm_work)
2377 rho_ao => rho_ao_weighted(1, :)
2381 rho_gspace=rho_g_dummy, &
2386 CALL dbcsr_set(rho_ao_weighted(i_spin, i_img)%matrix, 0.0_dp)
2392 n_x_start = lbound(ldos(1)%array, 1)
2393 n_x_end = ubound(ldos(1)%array, 1)
2394 n_y_start = lbound(ldos(1)%array, 2)
2395 n_y_end = ubound(ldos(1)%array, 2)
2396 n_z_start = lbound(ldos(1)%array, 3)
2397 n_z_end = ubound(ldos(1)%array, 3)
2402 DO i_x = n_x_start, n_x_end
2403 DO i_y = n_y_start, n_y_end
2404 DO i_z = n_z_start, n_z_end
2406 avg_ldos_occ = 0.0_dp
2408 avg_ldos_occ = avg_ldos_occ + ldos(i_e)%array(i_x, i_y, i_z)
2410 avg_ldos_occ = avg_ldos_occ/real(n_e_occ, kind=
dp)
2412 avg_ldos_virt = 0.0_dp
2413 DO i_e = n_e_occ + 1, n_e
2414 avg_ldos_virt = avg_ldos_virt + ldos(i_e)%array(i_x, i_y, i_z)
2416 avg_ldos_virt = avg_ldos_virt/real(n_e - n_e_occ, kind=
dp)
2419 DO i_e = n_e_occ, 1, -1
2420 IF (ldos(i_e)%array(i_x, i_y, i_z) > mp2_env%ri_g0w0%ldos_thresh_print_loc_bandgap*avg_ldos_occ)
THEN
2421 e_vbm_rspace%array(i_x, i_y, i_z) = e_array(i_e)
2427 DO i_e = n_e_occ + 1, n_e
2428 IF (ldos(i_e)%array(i_x, i_y, i_z) > mp2_env%ri_g0w0%ldos_thresh_print_loc_bandgap*avg_ldos_virt)
THEN
2429 e_cbm_rspace%array(i_x, i_y, i_z) = e_array(i_e)
2441 CALL pw_copy(e_cbm_rspace, e_gap_rspace)
2442 CALL pw_axpy(e_vbm_rspace, e_gap_rspace, -1.0_dp)
2447 CALL print_file(e_gap_rspace, dft_gw_char//
"_Gap_in_eV", gw_section, particles, mp2_env)
2448 CALL print_file(e_vbm_rspace, dft_gw_char//
"_VBM_in_eV", gw_section, particles, mp2_env)
2449 CALL print_file(e_cbm_rspace, dft_gw_char//
"_CBM_in_eV", gw_section, particles, mp2_env)
2450 CALL print_file(ldos(n_e_occ), dft_gw_char//
"_LDOS_VBM_in_eV", gw_section, particles, mp2_env)
2451 CALL print_file(ldos(n_e_occ + 1), dft_gw_char//
"_LDOS_CBM_in_eV", gw_section, particles, mp2_env)
2457 DEALLOCATE (e_array)
2459 NULLIFY (kpoints_sigma%sab_nl)
2461 CALL timestop(handle)
2463 END SUBROUTINE calculate_e_gap_rspace
2473 SUBROUTINE print_file(pw_print, middle_name, gw_section, particles, mp2_env)
2475 CHARACTER(len=*) :: middle_name
2480 CHARACTER(LEN=*),
PARAMETER :: routinen =
'print_file'
2482 INTEGER :: handle, unit_nr_cube
2486 CALL timeset(routinen, handle)
2491 unit_nr_cube =
cp_print_key_unit_nr(logger, gw_section,
"PRINT%LOCAL_BANDGAP", extension=
".cube", &
2492 middle_name=middle_name, file_form=
"FORMATTED", mpi_io=mpi_io)
2493 CALL cp_pw_to_cube(pw_print, unit_nr_cube, middle_name, particles=particles, &
2494 stride=mp2_env%ri_g0w0%stride_loc_bandgap, mpi_io=mpi_io)
2496 "PRINT%LOCAL_BANDGAP", mpi_io=mpi_io)
2498 CALL timestop(handle)
2500 END SUBROUTINE print_file
2512 SUBROUTINE create_real_space_grids(E_gap_rspace, E_VBM_rspace, E_CBM_rspace, rho_g_dummy, LDOS, auxbas_pw_pool, qs_env)
2519 CHARACTER(LEN=*),
PARAMETER :: routinen =
'create_real_space_grids'
2521 INTEGER :: handle, i_e, n_e
2525 CALL timeset(routinen, handle)
2527 CALL get_qs_env(qs_env=qs_env, mp2_env=mp2_env, pw_env=pw_env)
2529 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2531 CALL auxbas_pw_pool%create_pw(e_gap_rspace)
2532 CALL auxbas_pw_pool%create_pw(e_vbm_rspace)
2533 CALL auxbas_pw_pool%create_pw(e_cbm_rspace)
2534 CALL auxbas_pw_pool%create_pw(rho_g_dummy)
2536 n_e = int(mp2_env%ri_g0w0%energy_window_print_loc_bandgap/ &
2537 mp2_env%ri_g0w0%energy_spacing_print_loc_bandgap)
2539 ALLOCATE (ldos(n_e))
2542 CALL auxbas_pw_pool%create_pw(ldos(i_e))
2545 CALL timestop(handle)
2547 END SUBROUTINE create_real_space_grids
2574 SUBROUTINE calc_periodic_correction(delta_corr, qs_env, para_env, para_env_RPA, kp_grid, homo, nmo, &
2575 gw_corr_lev_occ, gw_corr_lev_virt, omega, fm_mo_coeff, Eigenval, &
2576 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
2577 first_cycle_periodic_correction, kpoints, do_mo_coeff_Gamma_only, &
2578 num_kp_grids, eps_kpoint, do_extra_kpoints, do_aux_bas, frac_aux_mos)
2580 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
2581 INTENT(INOUT) :: delta_corr
2584 INTEGER,
DIMENSION(:),
POINTER :: kp_grid
2585 INTEGER,
INTENT(IN) :: homo, nmo, gw_corr_lev_occ, &
2587 REAL(kind=
dp),
INTENT(IN) :: omega
2589 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eigenval
2590 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_re_mo_mo, &
2591 matrix_berry_im_mo_mo
2592 LOGICAL,
INTENT(INOUT) :: first_cycle_periodic_correction
2594 LOGICAL,
INTENT(IN) :: do_mo_coeff_gamma_only
2595 INTEGER,
INTENT(IN) :: num_kp_grids
2596 REAL(kind=
dp),
INTENT(IN) :: eps_kpoint
2597 LOGICAL,
INTENT(IN) :: do_extra_kpoints, do_aux_bas
2598 REAL(kind=
dp),
INTENT(IN) :: frac_aux_mos
2600 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_periodic_correction'
2603 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eps_head, eps_inv_head
2604 REAL(kind=
dp),
DIMENSION(3, 3) :: h_inv
2606 CALL timeset(routinen, handle)
2608 IF (first_cycle_periodic_correction)
THEN
2610 CALL get_kpoints(qs_env, kpoints, kp_grid, num_kp_grids, para_env, h_inv, nmo, do_mo_coeff_gamma_only, &
2613 CALL get_berry_phase(qs_env, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, fm_mo_coeff, &
2614 para_env, do_mo_coeff_gamma_only, homo, nmo, gw_corr_lev_virt, eps_kpoint, do_aux_bas, &
2619 CALL compute_eps_head_berry(eps_head, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, para_env_rpa, &
2620 qs_env, homo, eigenval, omega)
2622 CALL compute_eps_inv_head(eps_inv_head, eps_head, kpoints)
2624 CALL kpoint_sum_for_eps_inv_head_berry(delta_corr, eps_inv_head, kpoints, qs_env, &
2625 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
2626 homo, gw_corr_lev_occ, gw_corr_lev_virt, para_env_rpa, &
2629 DEALLOCATE (eps_head, eps_inv_head)
2631 first_cycle_periodic_correction = .false.
2633 CALL timestop(handle)
2635 END SUBROUTINE calc_periodic_correction
2649 SUBROUTINE compute_eps_head_berry(eps_head, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, para_env_RPA, &
2650 qs_env, homo, Eigenval, omega)
2652 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
2653 INTENT(OUT) :: eps_head
2655 TYPE(
dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_berry_re_mo_mo, &
2656 matrix_berry_im_mo_mo
2659 INTEGER,
INTENT(IN) :: homo
2660 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eigenval
2661 REAL(kind=
dp),
INTENT(IN) :: omega
2663 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_eps_head_Berry'
2665 INTEGER :: col, col_end_in_block, col_offset, col_size, handle, i_col, i_row, ikp, nkp, nmo, &
2666 row, row_offset, row_size, row_start_in_block
2667 REAL(kind=
dp) :: abs_k_square, cell_volume, &
2668 correct_kpoint(3), cos_square, &
2669 eigen_diff, relative_kpoint(3), &
2671 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: p_head
2672 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_block
2676 CALL timeset(routinen, handle)
2679 CALL get_cell(cell=cell, deth=cell_volume)
2681 NULLIFY (data_block)
2685 nmo =
SIZE(eigenval)
2687 ALLOCATE (p_head(nkp))
2690 ALLOCATE (eps_head(nkp))
2691 eps_head(:) = 0.0_dp
2695 relative_kpoint(1:3) = matmul(cell%hmat, kpoints%xkp(1:3, ikp))
2697 correct_kpoint(1:3) =
twopi*kpoints%xkp(1:3, ikp)
2699 abs_k_square = (correct_kpoint(1))**2 + (correct_kpoint(2))**2 + (correct_kpoint(3))**2
2706 row_size=row_size, col_size=col_size, &
2707 row_offset=row_offset, col_offset=col_offset)
2709 IF (row_offset + row_size <= homo .OR. col_offset > homo) cycle
2711 IF (row_offset <= homo)
THEN
2712 row_start_in_block = homo - row_offset + 2
2714 row_start_in_block = 1
2717 IF (col_offset + col_size - 1 > homo)
THEN
2718 col_end_in_block = homo - col_offset + 1
2720 col_end_in_block = col_size
2723 DO i_row = row_start_in_block, min(row_size, nmo - row_offset + 1)
2725 DO i_col = 1, min(col_end_in_block, nmo - col_offset + 1)
2727 eigen_diff = eigenval(i_col + col_offset - 1) - eigenval(i_row + row_offset - 1)
2729 cos_square = (data_block(i_row, i_col))**2
2731 p_head(ikp) = p_head(ikp) + 2.0_dp*eigen_diff/(omega**2 + eigen_diff**2)*cos_square/abs_k_square
2746 row_size=row_size, col_size=col_size, &
2747 row_offset=row_offset, col_offset=col_offset)
2749 IF (row_offset + row_size <= homo .OR. col_offset > homo) cycle
2751 IF (row_offset <= homo)
THEN
2752 row_start_in_block = homo - row_offset + 2
2754 row_start_in_block = 1
2757 IF (col_offset + col_size - 1 > homo)
THEN
2758 col_end_in_block = homo - col_offset + 1
2760 col_end_in_block = col_size
2763 DO i_row = row_start_in_block, min(row_size, nmo - row_offset + 1)
2765 DO i_col = 1, min(col_end_in_block, nmo - col_offset + 1)
2767 eigen_diff = eigenval(i_col + col_offset - 1) - eigenval(i_row + row_offset - 1)
2769 sin_square = (data_block(i_row, i_col))**2
2771 p_head(ikp) = p_head(ikp) + 2.0_dp*eigen_diff/(omega**2 + eigen_diff**2)*sin_square/abs_k_square
2783 CALL para_env_rpa%sum(p_head)
2787 eps_head(:) = 1.0_dp - 2.0_dp*p_head(:)/cell_volume*
fourpi
2791 CALL timestop(handle)
2793 END SUBROUTINE compute_eps_head_berry
2811 SUBROUTINE get_berry_phase(qs_env, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, fm_mo_coeff, para_env, &
2812 do_mo_coeff_Gamma_only, homo, nmo, gw_corr_lev_virt, eps_kpoint, do_aux_bas, &
2816 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_re_mo_mo, &
2817 matrix_berry_im_mo_mo
2820 LOGICAL,
INTENT(IN) :: do_mo_coeff_gamma_only
2821 INTEGER,
INTENT(IN) :: homo, nmo, gw_corr_lev_virt
2822 REAL(kind=
dp),
INTENT(IN) :: eps_kpoint
2823 LOGICAL,
INTENT(IN) :: do_aux_bas
2824 REAL(kind=
dp),
INTENT(IN) :: frac_aux_mos
2826 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_berry_phase'
2828 INTEGER :: col_index, handle, i_col_local, ikind, &
2829 ikp, nao_aux, ncol_local, nkind, nkp, &
2831 INTEGER,
DIMENSION(:),
POINTER :: col_indices
2832 REAL(
dp) :: abs_kpoint, correct_kpoint(3), &
2834 REAL(kind=
dp),
DIMENSION(:),
POINTER :: evals_p, evals_p_sqrt_inv
2837 TYPE(
cp_fm_type) :: fm_mat_eigv_p, fm_mat_p, fm_mat_p_sqrt_inv, fm_mat_s_aux_aux_inv, &
2838 fm_mat_scaled_eigv_p, fm_mat_work_aux_aux
2839 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, matrix_s_aux_aux, &
2841 TYPE(
dbcsr_type),
POINTER :: cosmat, cosmat_desymm, mat_mo_coeff_aux, mat_mo_coeff_aux_2, &
2842 mat_mo_coeff_gamma_all, mat_mo_coeff_gamma_occ_and_gw, mat_mo_coeff_im, mat_mo_coeff_re, &
2843 mat_work_aux_orb, mat_work_aux_orb_2, matrix_p, matrix_p_sqrt, matrix_p_sqrt_inv, &
2844 matrix_s_inv_aux_aux, sinmat, sinmat_desymm, tmp
2848 POINTER :: sab_orb, sab_orb_mic, sgwgw_list, &
2850 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2854 CALL timeset(routinen, handle)
2858 NULLIFY (matrix_berry_re_mo_mo, matrix_s, cell, matrix_berry_im_mo_mo, sinmat, cosmat, tmp, &
2859 cosmat_desymm, sinmat_desymm, qs_kind_set, orb_basis_set_list, sab_orb_mic)
2863 matrix_s=matrix_s, &
2864 qs_kind_set=qs_kind_set, &
2869 ALLOCATE (orb_basis_set_list(nkind))
2875 NULLIFY (mat_mo_coeff_re)
2878 template=matrix_s(1)%matrix, &
2879 matrix_type=dbcsr_type_no_symmetry)
2881 NULLIFY (mat_mo_coeff_im)
2884 template=matrix_s(1)%matrix, &
2885 matrix_type=dbcsr_type_no_symmetry)
2887 NULLIFY (mat_mo_coeff_gamma_all)
2890 template=matrix_s(1)%matrix, &
2891 matrix_type=dbcsr_type_no_symmetry)
2893 CALL copy_fm_to_dbcsr(fm_mo_coeff, mat_mo_coeff_gamma_all, keep_sparsity=.false.)
2895 NULLIFY (mat_mo_coeff_gamma_occ_and_gw)
2897 CALL dbcsr_create(matrix=mat_mo_coeff_gamma_occ_and_gw, &
2898 template=matrix_s(1)%matrix, &
2899 matrix_type=dbcsr_type_no_symmetry)
2901 CALL copy_fm_to_dbcsr(fm_mo_coeff, mat_mo_coeff_gamma_occ_and_gw, keep_sparsity=.false.)
2903 IF (.NOT. do_aux_bas)
THEN
2911 CALL dbcsr_create(matrix=cosmat, template=matrix_s(1)%matrix)
2912 CALL dbcsr_create(matrix=sinmat, template=matrix_s(1)%matrix)
2914 template=matrix_s(1)%matrix, &
2915 matrix_type=dbcsr_type_no_symmetry)
2917 template=matrix_s(1)%matrix, &
2918 matrix_type=dbcsr_type_no_symmetry)
2920 template=matrix_s(1)%matrix, &
2921 matrix_type=dbcsr_type_no_symmetry)
2932 NULLIFY (gw_aux_basis_set_list)
2933 ALLOCATE (gw_aux_basis_set_list(nkind))
2937 NULLIFY (gw_aux_basis_set_list(ikind)%gto_basis_set)
2939 NULLIFY (basis_set_gw_aux)
2941 qs_kind => qs_kind_set(ikind)
2942 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_gw_aux, basis_type=
"AUX_GW")
2943 cpassert(
ASSOCIATED(basis_set_gw_aux))
2945 basis_set_gw_aux%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius
2947 gw_aux_basis_set_list(ikind)%gto_basis_set => basis_set_gw_aux
2952 NULLIFY (sgwgw_list, sgworb_list)
2954 CALL setup_neighbor_list(sgworb_list, gw_aux_basis_set_list, orb_basis_set_list, qs_env=qs_env)
2956 NULLIFY (matrix_s_aux_aux, matrix_s_aux_orb)
2960 gw_aux_basis_set_list, gw_aux_basis_set_list, sgwgw_list)
2963 gw_aux_basis_set_list, orb_basis_set_list, sgworb_list)
2965 CALL dbcsr_get_info(matrix_s_aux_aux(1)%matrix, nfullrows_total=nao_aux)
2967 nmo_for_aux_bas = floor(frac_aux_mos*real(nao_aux, kind=
dp))
2970 context=fm_mo_coeff%matrix_struct%context, &
2971 nrow_global=nao_aux, &
2972 ncol_global=nao_aux, &
2975 NULLIFY (mat_work_aux_orb)
2978 template=matrix_s_aux_orb(1)%matrix, &
2979 matrix_type=dbcsr_type_no_symmetry)
2981 NULLIFY (mat_work_aux_orb_2)
2984 template=matrix_s_aux_orb(1)%matrix, &
2985 matrix_type=dbcsr_type_no_symmetry)
2987 NULLIFY (mat_mo_coeff_aux)
2990 template=matrix_s_aux_orb(1)%matrix, &
2991 matrix_type=dbcsr_type_no_symmetry)
2993 NULLIFY (mat_mo_coeff_aux_2)
2996 template=matrix_s_aux_orb(1)%matrix, &
2997 matrix_type=dbcsr_type_no_symmetry)
2999 NULLIFY (matrix_s_inv_aux_aux)
3002 template=matrix_s_aux_aux(1)%matrix, &
3003 matrix_type=dbcsr_type_no_symmetry)
3008 template=matrix_s(1)%matrix, &
3009 matrix_type=dbcsr_type_no_symmetry)
3011 NULLIFY (matrix_p_sqrt)
3014 template=matrix_s(1)%matrix, &
3015 matrix_type=dbcsr_type_no_symmetry)
3017 NULLIFY (matrix_p_sqrt_inv)
3020 template=matrix_s(1)%matrix, &
3021 matrix_type=dbcsr_type_no_symmetry)
3023 CALL cp_fm_create(fm_mat_s_aux_aux_inv, fm_struct_aux_aux, name=
"inverse overlap mat")
3024 CALL cp_fm_create(fm_mat_work_aux_aux, fm_struct_aux_aux, name=
"work mat")
3026 CALL cp_fm_create(fm_mat_eigv_p, fm_mo_coeff%matrix_struct)
3027 CALL cp_fm_create(fm_mat_scaled_eigv_p, fm_mo_coeff%matrix_struct)
3028 CALL cp_fm_create(fm_mat_p_sqrt_inv, fm_mo_coeff%matrix_struct)
3031 ALLOCATE (evals_p(nmo))
3033 NULLIFY (evals_p_sqrt_inv)
3034 ALLOCATE (evals_p_sqrt_inv(nmo))
3043 CALL copy_fm_to_dbcsr(fm_mat_s_aux_aux_inv, matrix_s_inv_aux_aux, keep_sparsity=.false.)
3045 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, &
3046 filter_eps=1.0e-15_dp)
3048 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, &
3049 last_column=nmo_for_aux_bas, filter_eps=1.0e-15_dp)
3051 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, &
3052 filter_eps=1.0e-15_dp)
3054 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_aux_2, mat_work_aux_orb, 0.0_dp, matrix_p, &
3055 filter_eps=1.0e-15_dp)
3059 CALL cp_fm_syevd(fm_mat_p, fm_mat_eigv_p, evals_p)
3062 evals_p_sqrt_inv(1:nmo - nmo_for_aux_bas) = 0.0_dp
3063 evals_p_sqrt_inv(nmo - nmo_for_aux_bas + 1:nmo) = 1.0_dp/sqrt(evals_p(nmo - nmo_for_aux_bas + 1:nmo))
3065 CALL cp_fm_to_fm(fm_mat_eigv_p, fm_mat_scaled_eigv_p)
3068 ncol_local=ncol_local, &
3069 col_indices=col_indices)
3071 CALL para_env%sync()
3074 DO i_col_local = 1, ncol_local
3076 col_index = col_indices(i_col_local)
3078 fm_mat_scaled_eigv_p%local_data(:, i_col_local) = &
3079 fm_mat_scaled_eigv_p%local_data(:, i_col_local)*evals_p_sqrt_inv(col_index)
3083 CALL para_env%sync()
3085 CALL parallel_gemm(transa=
"N", transb=
"T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, &
3086 matrix_a=fm_mat_eigv_p, matrix_b=fm_mat_scaled_eigv_p, beta=0.0_dp, &
3087 matrix_c=fm_mat_p_sqrt_inv)
3089 CALL copy_fm_to_dbcsr(fm_mat_p_sqrt_inv, matrix_p_sqrt_inv, keep_sparsity=.false.)
3091 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, mat_mo_coeff_aux_2, matrix_p_sqrt_inv, 0.0_dp, mat_mo_coeff_aux, &
3092 filter_eps=1.0e-15_dp)
3100 CALL dbcsr_create(matrix=cosmat, template=matrix_s_aux_aux(1)%matrix)
3101 CALL dbcsr_create(matrix=sinmat, template=matrix_s_aux_aux(1)%matrix)
3103 template=matrix_s_aux_orb(1)%matrix, &
3104 matrix_type=dbcsr_type_no_symmetry)
3106 template=matrix_s_aux_aux(1)%matrix, &
3107 matrix_type=dbcsr_type_no_symmetry)
3109 template=matrix_s_aux_aux(1)%matrix, &
3110 matrix_type=dbcsr_type_no_symmetry)
3111 CALL dbcsr_copy(cosmat, matrix_s_aux_aux(1)%matrix)
3112 CALL dbcsr_copy(sinmat, matrix_s_aux_aux(1)%matrix)
3123 NULLIFY (mat_mo_coeff_gamma_all)
3126 template=matrix_s_aux_orb(1)%matrix, &
3127 matrix_type=dbcsr_type_no_symmetry)
3129 CALL dbcsr_copy(mat_mo_coeff_gamma_all, mat_mo_coeff_aux)
3131 NULLIFY (mat_mo_coeff_gamma_occ_and_gw)
3133 CALL dbcsr_create(matrix=mat_mo_coeff_gamma_occ_and_gw, &
3134 template=matrix_s_aux_orb(1)%matrix, &
3135 matrix_type=dbcsr_type_no_symmetry)
3137 CALL dbcsr_copy(mat_mo_coeff_gamma_occ_and_gw, mat_mo_coeff_aux)
3139 DEALLOCATE (evals_p, evals_p_sqrt_inv)
3143 CALL remove_unnecessary_blocks(mat_mo_coeff_gamma_occ_and_gw, homo, gw_corr_lev_virt)
3147 ALLOCATE (matrix_berry_re_mo_mo(ikp)%matrix)
3150 template=matrix_s(1)%matrix, &
3151 matrix_type=dbcsr_type_no_symmetry)
3153 CALL dbcsr_set(matrix_berry_re_mo_mo(ikp)%matrix, 0.0_dp)
3155 ALLOCATE (matrix_berry_im_mo_mo(ikp)%matrix)
3158 template=matrix_s(1)%matrix, &
3159 matrix_type=dbcsr_type_no_symmetry)
3161 CALL dbcsr_set(matrix_berry_im_mo_mo(ikp)%matrix, 0.0_dp)
3163 correct_kpoint(1:3) = -
twopi*kpoints%xkp(1:3, ikp)
3165 abs_kpoint = sqrt(correct_kpoint(1)**2 + correct_kpoint(2)**2 + correct_kpoint(3)**2)
3167 IF (abs_kpoint < eps_kpoint)
THEN
3169 scale_kpoint = eps_kpoint/abs_kpoint
3170 correct_kpoint(:) = correct_kpoint(:)*scale_kpoint
3175 IF (do_aux_bas)
THEN
3177 basis_type=
"AUX_GW")
3183 IF (do_mo_coeff_gamma_only)
THEN
3187 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_gamma_occ_and_gw, 0.0_dp, tmp, &
3188 filter_eps=1.0e-15_dp)
3190 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 0.0_dp, &
3191 matrix_berry_re_mo_mo(ikp)%matrix, filter_eps=1.0e-15_dp)
3195 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_gamma_occ_and_gw, 0.0_dp, tmp, &
3196 filter_eps=1.0e-15_dp)
3198 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 0.0_dp, &
3199 matrix_berry_im_mo_mo(ikp)%matrix, filter_eps=1.0e-15_dp)
3205 mat_mo_coeff_re, keep_sparsity=.false.)
3208 mat_mo_coeff_im, keep_sparsity=.false.)
3215 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_re, 0.0_dp, tmp)
3218 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 0.0_dp, &
3219 matrix_berry_re_mo_mo(ikp)%matrix)
3222 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_re, 0.0_dp, tmp)
3225 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 0.0_dp, &
3226 matrix_berry_im_mo_mo(ikp)%matrix)
3229 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_im, 0.0_dp, tmp)
3232 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 1.0_dp, &
3233 matrix_berry_im_mo_mo(ikp)%matrix)
3236 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_im, 0.0_dp, tmp)
3239 CALL dbcsr_multiply(
'T',
'N', -1.0_dp, mat_mo_coeff_gamma_all, tmp, 1.0_dp, &
3240 matrix_berry_re_mo_mo(ikp)%matrix)
3244 IF (abs_kpoint < eps_kpoint)
THEN
3246 CALL dbcsr_scale(matrix_berry_im_mo_mo(ikp)%matrix, 1.0_dp/scale_kpoint)
3247 CALL dbcsr_set(matrix_berry_re_mo_mo(ikp)%matrix, 0.0_dp)
3263 DEALLOCATE (orb_basis_set_list)
3267 IF (do_aux_bas)
THEN
3269 DEALLOCATE (gw_aux_basis_set_list)
3296 CALL timestop(handle)
3298 END SUBROUTINE get_berry_phase
3306 SUBROUTINE remove_unnecessary_blocks(mat_mo_coeff_Gamma_occ_and_GW, homo, gw_corr_lev_virt)
3308 TYPE(
dbcsr_type),
POINTER :: mat_mo_coeff_gamma_occ_and_gw
3309 INTEGER,
INTENT(IN) :: homo, gw_corr_lev_virt
3311 INTEGER :: col, col_offset, row
3312 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_block
3320 col_offset=col_offset)
3322 IF (col_offset > homo + gw_corr_lev_virt)
THEN
3332 CALL dbcsr_filter(mat_mo_coeff_gamma_occ_and_gw, 1.0e-15_dp)
3334 END SUBROUTINE remove_unnecessary_blocks
3350 SUBROUTINE kpoint_sum_for_eps_inv_head_berry(delta_corr, eps_inv_head, kpoints, qs_env, matrix_berry_re_mo_mo, &
3351 matrix_berry_im_mo_mo, homo, gw_corr_lev_occ, gw_corr_lev_virt, &
3352 para_env_RPA, do_extra_kpoints)
3354 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
3355 INTENT(INOUT) :: delta_corr
3356 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eps_inv_head
3359 TYPE(
dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_berry_re_mo_mo, &
3360 matrix_berry_im_mo_mo
3361 INTEGER,
INTENT(IN) :: homo, gw_corr_lev_occ, gw_corr_lev_virt
3363 LOGICAL,
INTENT(IN) :: do_extra_kpoints
3365 INTEGER :: col, col_offset, col_size, i_col, i_row, &
3366 ikp, m_level, n_level_gw, nkp, row, &
3367 row_offset, row_size
3368 REAL(kind=
dp) :: abs_k_square, cell_volume, &
3369 check_int_one_over_ksq, contribution, &
3371 REAL(kind=
dp),
DIMENSION(3) :: correct_kpoint
3372 REAL(kind=
dp),
DIMENSION(:),
POINTER :: delta_corr_extra
3373 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_block
3379 CALL get_cell(cell=cell, deth=cell_volume)
3385 IF (do_extra_kpoints)
THEN
3386 NULLIFY (delta_corr_extra)
3387 ALLOCATE (delta_corr_extra(1 + homo - gw_corr_lev_occ:homo + gw_corr_lev_virt))
3388 delta_corr_extra = 0.0_dp
3391 check_int_one_over_ksq = 0.0_dp
3395 weight = kpoints%wkp(ikp)
3397 correct_kpoint(1:3) =
twopi*kpoints%xkp(1:3, ikp)
3399 abs_k_square = (correct_kpoint(1))**2 + (correct_kpoint(2))**2 + (correct_kpoint(3))**2
3406 row_size=row_size, col_size=col_size, &
3407 row_offset=row_offset, col_offset=col_offset)
3409 DO i_col = 1, col_size
3411 DO n_level_gw = 1 + homo - gw_corr_lev_occ, homo + gw_corr_lev_virt
3413 IF (n_level_gw == i_col + col_offset - 1)
THEN
3415 DO i_row = 1, row_size
3417 contribution = weight*(eps_inv_head(ikp) - 1.0_dp)/abs_k_square*(data_block(i_row, i_col))**2
3419 m_level = i_row + row_offset - 1
3422 IF (m_level /= n_level_gw) cycle
3424 IF (.NOT. do_extra_kpoints)
THEN
3426 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution
3430 IF (ikp <= nkp*8/9)
THEN
3432 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution
3436 delta_corr_extra(n_level_gw) = delta_corr_extra(n_level_gw) + contribution
3459 row_size=row_size, col_size=col_size, &
3460 row_offset=row_offset, col_offset=col_offset)
3462 DO i_col = 1, col_size
3464 DO n_level_gw = 1 + homo - gw_corr_lev_occ, homo + gw_corr_lev_virt
3466 IF (n_level_gw == i_col + col_offset - 1)
THEN
3468 DO i_row = 1, row_size
3470 m_level = i_row + row_offset - 1
3472 contribution = weight*(eps_inv_head(ikp) - 1.0_dp)/abs_k_square*(data_block(i_row, i_col))**2
3475 IF (m_level /= n_level_gw) cycle
3477 IF (.NOT. do_extra_kpoints)
THEN
3479 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution
3483 IF (ikp <= nkp*8/9)
THEN
3485 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution
3489 delta_corr_extra(n_level_gw) = delta_corr_extra(n_level_gw) + contribution
3507 check_int_one_over_ksq = check_int_one_over_ksq + weight/abs_k_square
3512 delta_corr = delta_corr/cell_volume*
fourpi
3514 check_int_one_over_ksq = check_int_one_over_ksq/cell_volume
3516 CALL para_env_rpa%sum(delta_corr)
3518 IF (do_extra_kpoints)
THEN
3520 delta_corr_extra = delta_corr_extra/cell_volume*
fourpi
3522 CALL para_env_rpa%sum(delta_corr_extra)
3524 delta_corr(:) = delta_corr(:) + (delta_corr(:) - delta_corr_extra(:))
3526 DEALLOCATE (delta_corr_extra)
3530 END SUBROUTINE kpoint_sum_for_eps_inv_head_berry
3538 SUBROUTINE compute_eps_inv_head(eps_inv_head, eps_head, kpoints)
3539 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
3540 INTENT(OUT) :: eps_inv_head
3541 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eps_head
3544 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_eps_inv_head'
3546 INTEGER :: handle, ikp, nkp
3548 CALL timeset(routinen, handle)
3552 ALLOCATE (eps_inv_head(nkp))
3556 eps_inv_head(ikp) = 1.0_dp/eps_head(ikp)
3560 CALL timestop(handle)
3562 END SUBROUTINE compute_eps_inv_head
3576 SUBROUTINE get_kpoints(qs_env, kpoints, kp_grid, num_kp_grids, para_env, h_inv, nmo, &
3577 do_mo_coeff_Gamma_only, do_extra_kpoints)
3580 INTEGER,
DIMENSION(:),
POINTER :: kp_grid
3581 INTEGER,
INTENT(IN) :: num_kp_grids
3583 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT) :: h_inv
3584 INTEGER,
INTENT(IN) :: nmo
3585 LOGICAL,
INTENT(IN) :: do_mo_coeff_gamma_only, do_extra_kpoints
3587 INTEGER :: end_kp, i, i_grid_level, ix, iy, iz, &
3588 nkp_inner_grid, nkp_outer_grid, &
3590 INTEGER,
DIMENSION(3) :: outer_kp_grid
3591 REAL(kind=
dp) :: kpoint_weight_left, single_weight
3592 REAL(kind=
dp),
DIMENSION(3) :: kpt_latt, reducing_factor
3596 NULLIFY (kpoints, cell, particle_set)
3599 cpassert(mod(kp_grid(1)*kp_grid(2)*kp_grid(3), 2) == 0)
3600 IF (do_extra_kpoints)
THEN
3601 cpassert(do_mo_coeff_gamma_only)
3604 IF (do_mo_coeff_gamma_only)
THEN
3606 outer_kp_grid(1) = kp_grid(1) - 1
3607 outer_kp_grid(2) = kp_grid(2) - 1
3608 outer_kp_grid(3) = kp_grid(3) - 1
3610 CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set)
3616 kpoints%kp_scheme =
"GENERAL"
3617 kpoints%symmetry = .false.
3618 kpoints%verbose = .false.
3619 kpoints%full_grid = .false.
3620 kpoints%use_real_wfn = .false.
3621 kpoints%eps_geo = 1.e-6_dp
3622 npoints = kp_grid(1)*kp_grid(2)*kp_grid(3)/2 + &
3623 (num_kp_grids - 1)*((outer_kp_grid(1) + 1)/2*outer_kp_grid(2)*outer_kp_grid(3) - 1)
3625 IF (do_extra_kpoints)
THEN
3627 cpassert(num_kp_grids == 1)
3628 cpassert(mod(kp_grid(1), 4) == 0)
3629 cpassert(mod(kp_grid(2), 4) == 0)
3630 cpassert(mod(kp_grid(3), 4) == 0)
3634 IF (do_extra_kpoints)
THEN
3636 npoints = kp_grid(1)*kp_grid(2)*kp_grid(3)/2 + kp_grid(1)*kp_grid(2)*kp_grid(3)/2/8
3640 kpoints%full_grid = .true.
3641 kpoints%nkp = npoints
3642 ALLOCATE (kpoints%xkp(3, npoints), kpoints%wkp(npoints))
3643 kpoints%xkp = 0.0_dp
3644 kpoints%wkp = 0.0_dp
3646 nkp_outer_grid = outer_kp_grid(1)*outer_kp_grid(2)*outer_kp_grid(3)
3647 nkp_inner_grid = kp_grid(1)*kp_grid(2)*kp_grid(3)
3650 reducing_factor(:) = 1.0_dp
3651 kpoint_weight_left = 1.0_dp
3654 DO i_grid_level = 1, num_kp_grids - 1
3656 single_weight = kpoint_weight_left/real(nkp_outer_grid, kind=
dp)
3660 DO ix = 1, outer_kp_grid(1)
3661 DO iy = 1, outer_kp_grid(2)
3662 DO iz = 1, outer_kp_grid(3)
3665 IF (2*ix - outer_kp_grid(1) - 1 == 0 .AND. 2*iy - outer_kp_grid(2) - 1 == 0 .AND. &
3666 2*iz - outer_kp_grid(3) - 1 == 0) cycle
3669 IF (2*ix - outer_kp_grid(1) - 1 < 0) cycle
3672 kpt_latt(1) = real(2*ix - outer_kp_grid(1) - 1, kind=
dp)/(2._dp*real(outer_kp_grid(1), kind=
dp)) &
3674 kpt_latt(2) = real(2*iy - outer_kp_grid(2) - 1, kind=
dp)/(2._dp*real(outer_kp_grid(2), kind=
dp)) &
3676 kpt_latt(3) = real(2*iz - outer_kp_grid(3) - 1, kind=
dp)/(2._dp*real(outer_kp_grid(3), kind=
dp)) &
3678 kpoints%xkp(1:3, i) = matmul(transpose(h_inv), kpt_latt(:))
3680 IF (2*ix - outer_kp_grid(1) - 1 == 0)
THEN
3681 kpoints%wkp(i) = single_weight
3683 kpoints%wkp(i) = 2._dp*single_weight
3692 kpoint_weight_left = kpoint_weight_left - sum(kpoints%wkp(start_kp:end_kp))
3694 reducing_factor(1) = reducing_factor(1)/real(outer_kp_grid(1), kind=
dp)
3695 reducing_factor(2) = reducing_factor(2)/real(outer_kp_grid(2), kind=
dp)
3696 reducing_factor(3) = reducing_factor(3)/real(outer_kp_grid(3), kind=
dp)
3700 single_weight = kpoint_weight_left/real(nkp_inner_grid, kind=
dp)
3703 DO ix = 1, kp_grid(1)
3704 DO iy = 1, kp_grid(2)
3705 DO iz = 1, kp_grid(3)
3708 IF (2*ix - kp_grid(1) - 1 < 0) cycle
3711 kpt_latt(1) = real(2*ix - kp_grid(1) - 1, kind=
dp)/(2._dp*real(kp_grid(1), kind=
dp))*reducing_factor(1)
3712 kpt_latt(2) = real(2*iy - kp_grid(2) - 1, kind=
dp)/(2._dp*real(kp_grid(2), kind=
dp))*reducing_factor(2)
3713 kpt_latt(3) = real(2*iz - kp_grid(3) - 1, kind=
dp)/(2._dp*real(kp_grid(3), kind=
dp))*reducing_factor(3)
3715 kpoints%xkp(1:3, i) = matmul(transpose(h_inv), kpt_latt(:))
3717 kpoints%wkp(i) = 2._dp*single_weight
3723 IF (do_extra_kpoints)
THEN
3725 single_weight = kpoint_weight_left/real(kp_grid(1)*kp_grid(2)*kp_grid(3)/8, kind=
dp)
3727 DO ix = 1, kp_grid(1)/2
3728 DO iy = 1, kp_grid(2)/2
3729 DO iz = 1, kp_grid(3)/2
3732 IF (2*ix - kp_grid(1)/2 - 1 < 0) cycle
3735 kpt_latt(1) = real(2*ix - kp_grid(1)/2 - 1, kind=
dp)/(real(kp_grid(1), kind=
dp))
3736 kpt_latt(2) = real(2*iy - kp_grid(2)/2 - 1, kind=
dp)/(real(kp_grid(2), kind=
dp))
3737 kpt_latt(3) = real(2*iz - kp_grid(3)/2 - 1, kind=
dp)/(real(kp_grid(3), kind=
dp))
3739 kpoints%xkp(1:3, i) = matmul(transpose(h_inv), kpt_latt(:))
3741 kpoints%wkp(i) = 2._dp*single_weight
3750 ALLOCATE (kpoints%kp_sym(kpoints%nkp))
3751 DO i = 1, kpoints%nkp
3752 NULLIFY (kpoints%kp_sym(i)%kpoint_sym)
3762 CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set)
3764 CALL calculate_kp_orbitals(qs_env_kp_gamma_only, kpoints,
"MONKHORST-PACK", nadd=nmo, mp_grid=kp_grid(1:3), &
3765 group_size_ext=para_env%num_pe)
3768 DEALLOCATE (qs_env_kp_gamma_only)
3773 END SUBROUTINE get_kpoints
3781 PURE SUBROUTINE average_degenerate_levels(vec_Sigma_c_gw, Eigenval_DFT, eps_eigenval)
3782 COMPLEX(KIND=dp),
DIMENSION(:, :, :), &
3783 INTENT(INOUT) :: vec_sigma_c_gw
3784 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: eigenval_dft
3785 REAL(kind=dp),
INTENT(IN) :: eps_eigenval
3787 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: avg_self_energy
3788 INTEGER :: degeneracy, first_degenerate_level, i_deg_level, i_level_gw, j_deg_level, jquad, &
3789 num_deg_levels, num_integ_points, num_levels_gw
3790 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: list_degenerate_levels
3792 num_levels_gw =
SIZE(vec_sigma_c_gw, 1)
3794 ALLOCATE (list_degenerate_levels(num_levels_gw))
3795 list_degenerate_levels = 1
3797 num_integ_points =
SIZE(vec_sigma_c_gw, 2)
3799 ALLOCATE (avg_self_energy(num_integ_points))
3801 DO i_level_gw = 2, num_levels_gw
3803 IF (abs(eigenval_dft(i_level_gw) - eigenval_dft(i_level_gw - 1)) < eps_eigenval)
THEN
3805 list_degenerate_levels(i_level_gw) = list_degenerate_levels(i_level_gw - 1)
3809 list_degenerate_levels(i_level_gw) = list_degenerate_levels(i_level_gw - 1) + 1
3815 num_deg_levels = list_degenerate_levels(num_levels_gw)
3817 DO i_deg_level = 1, num_deg_levels
3821 DO i_level_gw = 1, num_levels_gw
3823 IF (degeneracy == 0 .AND. i_deg_level == list_degenerate_levels(i_level_gw))
THEN
3825 first_degenerate_level = i_level_gw
3829 IF (i_deg_level == list_degenerate_levels(i_level_gw))
THEN
3831 degeneracy = degeneracy + 1
3837 DO jquad = 1, num_integ_points
3839 avg_self_energy(jquad) = sum(vec_sigma_c_gw(first_degenerate_level:first_degenerate_level + degeneracy - 1, jquad, 1)) &
3840 /real(degeneracy, kind=dp)
3844 DO j_deg_level = 0, degeneracy - 1
3846 vec_sigma_c_gw(first_degenerate_level + j_deg_level, :, 1) = avg_self_energy(:)
3852 END SUBROUTINE average_degenerate_levels
3875 SUBROUTINE fit_and_continuation_2pole(vec_gw_energ, vec_omega_fit_gw, &
3876 z_value, m_value, vec_Sigma_c_gw, vec_Sigma_x_minus_vxc_gw, &
3877 Eigenval, Eigenval_scf, n_level_gw, &
3878 gw_corr_lev_occ, gw_corr_lev_vir, num_poles, &
3879 num_fit_points, crossing_search, homo, stop_crit, &
3880 fermi_level_offset, do_gw_im_time)
3882 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: vec_gw_energ, vec_omega_fit_gw, z_value, &
3884 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: vec_sigma_c_gw
3885 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_sigma_x_minus_vxc_gw, eigenval, &
3887 INTEGER,
INTENT(IN) :: n_level_gw, gw_corr_lev_occ, &
3888 gw_corr_lev_vir, num_poles, &
3889 num_fit_points, crossing_search, homo
3890 REAL(kind=dp),
INTENT(IN) :: stop_crit, fermi_level_offset
3891 LOGICAL,
INTENT(IN) :: do_gw_im_time
3893 CHARACTER(LEN=*),
PARAMETER :: routinen =
'fit_and_continuation_2pole'
3895 COMPLEX(KIND=dp) :: func_val, rho1
3896 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: dlambda, dlambda_2, lambda, &
3897 lambda_without_offset, vec_b_gw, &
3899 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: mat_a_gw, mat_b_gw
3900 INTEGER :: handle4, ierr, iii, iiter, info, &
3901 integ_range, jjj, jquad, kkk, &
3902 max_iter_fit, n_level_gw_ref, num_var, &
3904 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ipiv
3905 LOGICAL :: could_exit
3906 REAL(kind=dp) :: chi2, chi2_old, delta, deriv_val_real, e_fermi, gw_energ, ldown, &
3907 level_energ_gw, lup, range_step, scalparam, sign_occ_virt, stat_error
3908 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: lambda_im, lambda_re, stat_errors, &
3909 vec_n_gw, vec_omega_fit_gw_sign
3910 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: mat_n_gw
3912 max_iter_fit = 10000
3914 num_var = 2*num_poles + 1
3915 ALLOCATE (lambda(num_var))
3917 ALLOCATE (lambda_without_offset(num_var))
3918 lambda_without_offset = z_zero
3919 ALLOCATE (lambda_re(num_var))
3921 ALLOCATE (lambda_im(num_var))
3924 ALLOCATE (vec_omega_fit_gw_sign(num_fit_points))
3926 IF (n_level_gw <= gw_corr_lev_occ)
THEN
3927 sign_occ_virt = -1.0_dp
3929 sign_occ_virt = 1.0_dp
3932 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
3934 DO jquad = 1, num_fit_points
3935 vec_omega_fit_gw_sign(jquad) = abs(vec_omega_fit_gw(jquad))*sign_occ_virt
3939 range_step = (vec_omega_fit_gw_sign(num_fit_points) - vec_omega_fit_gw_sign(1))/(num_poles - 1)
3940 DO iii = 1, num_poles
3941 lambda_im(2*iii + 1) = vec_omega_fit_gw_sign(1) + (iii - 1)*range_step
3943 range_step = (vec_omega_fit_gw_sign(num_fit_points) - vec_omega_fit_gw_sign(1))/num_poles
3944 DO iii = 1, num_poles
3945 lambda_re(2*iii + 1) = abs(vec_omega_fit_gw_sign(1) + (iii - 0.5_dp)*range_step)
3949 lambda(iii) = lambda_re(iii) + gaussi*lambda_im(iii)
3952 CALL calc_chi2(chi2_old, lambda, vec_sigma_c_gw, vec_omega_fit_gw_sign, num_poles, &
3953 num_fit_points, n_level_gw)
3955 ALLOCATE (mat_a_gw(num_poles + 1, num_poles + 1))
3956 ALLOCATE (vec_b_gw(num_poles + 1))
3957 ALLOCATE (ipiv(num_poles + 1))
3961 mat_a_gw(1:num_poles + 1, 1) = z_one
3962 integ_range = num_fit_points/num_poles
3963 DO kkk = 1, num_poles + 1
3964 xpos = (kkk - 1)*integ_range + 1
3965 xpos = min(xpos, num_fit_points)
3967 DO iii = 1, num_poles
3969 func_val = z_one/(gaussi*vec_omega_fit_gw_sign(xpos) - &
3970 cmplx(lambda_re(jjj + 1), lambda_im(jjj + 1), kind=dp))
3971 mat_a_gw(kkk, iii + 1) = func_val
3973 vec_b_gw(kkk) = vec_sigma_c_gw(n_level_gw, xpos)
3977 CALL zgetrf(num_poles + 1, num_poles + 1, mat_a_gw, num_poles + 1, ipiv, info)
3979 CALL zgetrs(
'N', num_poles + 1, 1, mat_a_gw, num_poles + 1, ipiv, vec_b_gw, num_poles + 1, info)
3981 lambda_re(1) = real(vec_b_gw(1))
3982 lambda_im(1) = aimag(vec_b_gw(1))
3983 DO iii = 1, num_poles
3985 lambda_re(jjj) = real(vec_b_gw(iii + 1))
3986 lambda_im(jjj) = aimag(vec_b_gw(iii + 1))
3989 DEALLOCATE (mat_a_gw)
3990 DEALLOCATE (vec_b_gw)
3993 ALLOCATE (mat_a_gw(num_var*2, num_var*2))
3994 ALLOCATE (mat_b_gw(num_fit_points, num_var*2))
3995 ALLOCATE (dlambda(num_fit_points))
3996 ALLOCATE (dlambda_2(num_fit_points))
3997 ALLOCATE (vec_b_gw(num_var*2))
3998 ALLOCATE (vec_b_gw_copy(num_var*2))
3999 ALLOCATE (ipiv(num_var*2))
4004 could_exit = .false.
4007 DO iiter = 1, max_iter_fit
4009 CALL timeset(routinen//
"_fit_loop_1", handle4)
4013 lambda(iii) = lambda_re(iii) + gaussi*lambda_im(iii)
4017 DO kkk = 1, num_fit_points
4018 func_val = lambda(1)
4019 DO iii = 1, num_poles
4021 func_val = func_val + lambda(jjj)/(vec_omega_fit_gw_sign(kkk)*gaussi - lambda(jjj + 1))
4023 dlambda(kkk) = vec_sigma_c_gw(n_level_gw, kkk) - func_val
4025 rho1 = sum(dlambda*dlambda)
4029 DO iii = 1, num_fit_points
4030 mat_b_gw(iii, 1) = 1.0_dp
4031 mat_b_gw(iii, num_var + 1) = gaussi
4033 DO iii = 1, num_poles
4035 DO kkk = 1, num_fit_points
4036 mat_b_gw(kkk, jjj) = 1.0_dp/(gaussi*vec_omega_fit_gw_sign(kkk) - lambda(jjj + 1))
4037 mat_b_gw(kkk, jjj + num_var) = gaussi/(gaussi*vec_omega_fit_gw_sign(kkk) - lambda(jjj + 1))
4038 mat_b_gw(kkk, jjj + 1) = lambda(jjj)/(gaussi*vec_omega_fit_gw_sign(kkk) - lambda(jjj + 1))**2
4039 mat_b_gw(kkk, jjj + 1 + num_var) = (-lambda_im(jjj) + gaussi*lambda_re(jjj))/ &
4040 (gaussi*vec_omega_fit_gw_sign(kkk) - lambda(jjj + 1))**2
4044 CALL timestop(handle4)
4046 CALL timeset(routinen//
"_fit_matmul_1", handle4)
4048 CALL zgemm(
'C',
'N', num_var*2, num_var*2, num_fit_points, z_one, mat_b_gw, num_fit_points, mat_b_gw, num_fit_points, &
4049 z_zero, mat_a_gw, num_var*2)
4050 CALL timestop(handle4)
4052 CALL timeset(routinen//
"_fit_zgemv_1", handle4)
4053 CALL zgemv(
'C', num_fit_points, num_var*2, z_one, mat_b_gw, num_fit_points, dlambda, 1, &
4054 z_zero, vec_b_gw, 1)
4056 CALL timestop(handle4)
4059 DO iii = 1, num_var*2
4060 mat_a_gw(iii, iii) = mat_a_gw(iii, iii) + scalparam*mat_a_gw(iii, iii)
4067 CALL timeset(routinen//
"_fit_lin_eq_2", handle4)
4069 CALL zgetrf(2*num_var, 2*num_var, mat_a_gw, 2*num_var, ipiv, info)
4071 CALL zgetrs(
'N', 2*num_var, 1, mat_a_gw, 2*num_var, ipiv, vec_b_gw, 2*num_var, info)
4073 CALL timestop(handle4)
4076 lambda(iii) = lambda_re(iii) + gaussi*lambda_im(iii) + vec_b_gw(iii) + vec_b_gw(iii + num_var)
4080 CALL calc_chi2(chi2, lambda, vec_sigma_c_gw, vec_omega_fit_gw_sign, num_poles, &
4081 num_fit_points, n_level_gw)
4084 IF (chi2 < 1.0e-30_dp)
EXIT
4086 IF (chi2 < chi2_old)
THEN
4087 scalparam = max(scalparam/ldown, 1e-12_dp)
4089 lambda_re(iii) = lambda_re(iii) + real(vec_b_gw(iii) + vec_b_gw(iii + num_var))
4090 lambda_im(iii) = lambda_im(iii) + aimag(vec_b_gw(iii) + vec_b_gw(iii + num_var))
4092 IF (chi2_old/chi2 - 1.0_dp < stop_crit) could_exit = .true.
4095 scalparam = scalparam*lup
4097 IF (scalparam > 100.0_dp .AND. could_exit)
EXIT
4099 IF (scalparam > 1e+10_dp) scalparam = 1e-4_dp
4103 IF (.NOT. do_gw_im_time)
THEN
4107 func_val = lambda(1)
4108 DO iii = 1, num_poles
4111 func_val = func_val + lambda(jjj)/(-lambda(jjj + 1))
4114 lambda_re(1) = lambda_re(1) - real(func_val) + real(vec_sigma_c_gw(n_level_gw, num_fit_points))
4115 lambda_im(1) = lambda_im(1) - aimag(func_val) + aimag(vec_sigma_c_gw(n_level_gw, num_fit_points))
4119 lambda_without_offset(:) = lambda(:)
4122 lambda(iii) = cmplx(lambda_re(iii), lambda_im(iii), kind=dp)
4125 IF (do_gw_im_time)
THEN
4128 e_fermi = 0.5_dp*(eigenval(homo) + eigenval(homo + 1))
4132 IF (n_level_gw <= gw_corr_lev_occ)
THEN
4133 e_fermi = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset
4135 e_fermi = minval(eigenval(homo + 1:homo + gw_corr_lev_vir)) - fermi_level_offset
4140 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot .OR. &
4141 crossing_search == ri_rpa_g0w0_crossing_newton)
THEN
4144 func_val = lambda(1)
4145 z_value(n_level_gw) = 1.0_dp
4146 DO iii = 1, num_poles
4148 z_value(n_level_gw) = z_value(n_level_gw) + real(lambda(jjj)/ &
4149 (eigenval(n_level_gw_ref) - e_fermi - lambda(jjj + 1))**2)
4150 func_val = func_val + lambda(jjj)/(eigenval(n_level_gw_ref) - e_fermi - lambda(jjj + 1))
4153 m_value(n_level_gw) = 1.0_dp - z_value(n_level_gw)
4154 z_value(n_level_gw) = 1.0_dp/z_value(n_level_gw)
4155 gw_energ = real(func_val)
4156 vec_gw_energ(n_level_gw) = gw_energ
4159 IF (crossing_search == ri_rpa_g0w0_crossing_newton)
THEN
4161 level_energ_gw = (eigenval_scf(n_level_gw_ref) - &
4162 m_value(n_level_gw)*eigenval(n_level_gw_ref) + &
4163 vec_gw_energ(n_level_gw) + &
4164 vec_sigma_x_minus_vxc_gw(n_level_gw_ref))* &
4171 func_val = lambda(1)
4172 z_value(n_level_gw) = 1.0_dp
4173 DO iii = 1, num_poles
4175 func_val = func_val + lambda(jjj)/(level_energ_gw - e_fermi - lambda(jjj + 1))
4179 deriv_val_real = -1.0_dp
4180 DO iii = 1, num_poles
4182 deriv_val_real = deriv_val_real + real(lambda(jjj))/((abs(level_energ_gw - e_fermi - lambda(jjj + 1)))**2) &
4183 - (real(lambda(jjj))*(level_energ_gw - e_fermi) - real(lambda(jjj)*conjg(lambda(jjj + 1))))* &
4184 2.0_dp*(level_energ_gw - e_fermi - real(lambda(jjj + 1)))/ &
4185 ((abs(level_energ_gw - e_fermi - lambda(jjj + 1)))**2)
4189 delta = (eigenval_scf(n_level_gw_ref) + vec_sigma_x_minus_vxc_gw(n_level_gw_ref) + real(func_val) - level_energ_gw)/ &
4192 level_energ_gw = level_energ_gw - delta
4194 IF (abs(delta) < 1.0e-08)
EXIT
4200 vec_gw_energ(n_level_gw) = real(func_val)
4201 z_value(n_level_gw) = 1.0_dp
4202 m_value(n_level_gw) = 0.0_dp
4207 cpabort(
"Only NONE, ZSHOT and NEWTON implemented for 2-pole model")
4217 CALL calc_chi2(chi2, lambda_without_offset, vec_sigma_c_gw, vec_omega_fit_gw_sign, num_poles, &
4218 num_fit_points, n_level_gw)
4221 stat_error = sqrt(chi2/num_fit_points)
4224 ALLOCATE (vec_n_gw(num_var*2))
4227 ALLOCATE (mat_n_gw(num_var*2, num_var*2))
4230 DO iii = 1, num_var*2
4231 CALL calc_mat_n(vec_n_gw(iii), lambda_without_offset, vec_sigma_c_gw, vec_omega_fit_gw_sign, &
4232 iii, iii, num_poles, num_fit_points, n_level_gw, 0.001_dp)
4235 DO iii = 1, num_var*2
4236 DO jjj = 1, num_var*2
4237 CALL calc_mat_n(mat_n_gw(iii, jjj), lambda_without_offset, vec_sigma_c_gw, vec_omega_fit_gw_sign, &
4238 iii, jjj, num_poles, num_fit_points, n_level_gw, 0.001_dp)
4242 CALL dgetrf(2*num_var, 2*num_var, mat_n_gw, 2*num_var, ipiv, info)
4245 CALL dgetri(2*num_var, mat_n_gw, 2*num_var, ipiv, vec_b_gw, 2*num_var, info)
4247 ALLOCATE (stat_errors(2*num_var))
4248 stat_errors = 0.0_dp
4250 DO iii = 1, 2*num_var
4251 stat_errors(iii) = sqrt(abs(mat_n_gw(iii, iii)))*stat_error
4254 DEALLOCATE (mat_n_gw)
4255 DEALLOCATE (vec_n_gw)
4256 DEALLOCATE (mat_a_gw)
4257 DEALLOCATE (mat_b_gw)
4258 DEALLOCATE (stat_errors)
4259 DEALLOCATE (dlambda)
4260 DEALLOCATE (dlambda_2)
4261 DEALLOCATE (vec_b_gw)
4262 DEALLOCATE (vec_b_gw_copy)
4264 DEALLOCATE (vec_omega_fit_gw_sign)
4266 DEALLOCATE (lambda_without_offset)
4267 DEALLOCATE (lambda_re)
4268 DEALLOCATE (lambda_im)
4270 END SUBROUTINE fit_and_continuation_2pole
4306 z_value, m_value, vec_Sigma_c_gw, vec_Sigma_x_minus_vxc_gw, &
4307 Eigenval, Eigenval_scf, do_hedin_shift, n_level_gw, &
4308 gw_corr_lev_occ, gw_corr_lev_vir, &
4309 nparam_pade, num_fit_points, crossing_search, homo, &
4310 fermi_level_offset, do_gw_im_time, print_self_energy, count_ev_sc_GW, &
4311 vec_gw_dos, dos_lower_bound, dos_precision, ndos, &
4312 min_level_self_energy, max_level_self_energy, &
4313 dos_eta, dos_min, dos_max, e_fermi_ext)
4316 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: vec_gw_energ
4317 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_omega_fit_gw
4318 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: z_value, m_value
4319 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: vec_sigma_c_gw
4320 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_sigma_x_minus_vxc_gw, eigenval, &
4322 LOGICAL,
INTENT(IN) :: do_hedin_shift
4323 INTEGER,
INTENT(IN) :: n_level_gw, gw_corr_lev_occ, &
4324 gw_corr_lev_vir, nparam_pade, &
4325 num_fit_points, crossing_search, homo
4326 REAL(kind=dp),
INTENT(IN) :: fermi_level_offset
4327 LOGICAL,
INTENT(IN) :: do_gw_im_time, print_self_energy
4328 INTEGER,
INTENT(IN) :: count_ev_sc_gw
4329 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:),
OPTIONAL :: vec_gw_dos
4330 REAL(kind=dp),
OPTIONAL :: dos_lower_bound, dos_precision
4331 INTEGER,
INTENT(IN),
OPTIONAL :: ndos, min_level_self_energy, &
4332 max_level_self_energy
4333 REAL(kind=dp),
OPTIONAL :: dos_eta
4334 INTEGER,
INTENT(IN),
OPTIONAL :: dos_min, dos_max
4335 REAL(kind=dp),
OPTIONAL :: e_fermi_ext
4337 CHARACTER(LEN=*),
PARAMETER :: routinen =
'continuation_pade'
4339 CHARACTER(LEN=5) :: string_level
4340 CHARACTER(len=default_path_length) :: filename
4341 COMPLEX(KIND=dp) :: sigma_c_pade, sigma_c_pade_im_freq
4342 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: coeff_pade, omega_points_pade, &
4344 INTEGER :: handle, i_omega, idos, iunit, jquad, &
4345 n_level_gw_ref, num_omega
4346 REAL(kind=dp) :: e_fermi, energy_val, hedin_shift, &
4347 level_energ_gw_start, omega, &
4348 omega_dos, omega_dos_pade_eval, &
4350 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: vec_omega_fit_gw_sign, &
4351 vec_omega_fit_gw_sign_reorder, &
4352 vec_sigma_imag, vec_sigma_real
4353 TYPE(cp_logger_type),
POINTER :: logger
4355 CALL timeset(routinen, handle)
4357 ALLOCATE (vec_omega_fit_gw_sign(num_fit_points))
4359 IF (n_level_gw <= gw_corr_lev_occ)
THEN
4360 sign_occ_virt = -1.0_dp
4362 sign_occ_virt = 1.0_dp
4365 DO jquad = 1, num_fit_points
4366 vec_omega_fit_gw_sign(jquad) = abs(vec_omega_fit_gw(jquad))*sign_occ_virt
4369 IF (do_gw_im_time)
THEN
4372 e_fermi = 0.5_dp*(eigenval(homo) + eigenval(homo + 1))
4376 IF (n_level_gw <= gw_corr_lev_occ)
THEN
4377 e_fermi = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset
4379 e_fermi = minval(eigenval(homo + 1:homo + gw_corr_lev_vir)) - fermi_level_offset
4383 IF (
PRESENT(e_fermi_ext)) e_fermi = e_fermi_ext
4385 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
4388 ALLOCATE (sigma_c_gw_reorder(num_fit_points))
4389 ALLOCATE (vec_omega_fit_gw_sign_reorder(num_fit_points))
4391 IF (do_gw_im_time)
THEN
4392 DO jquad = 1, num_fit_points
4393 sigma_c_gw_reorder(jquad) = vec_sigma_c_gw(n_level_gw, jquad)
4394 vec_omega_fit_gw_sign_reorder(jquad) = vec_omega_fit_gw_sign(jquad)
4397 DO jquad = 1, num_fit_points
4398 sigma_c_gw_reorder(jquad) = vec_sigma_c_gw(n_level_gw, num_fit_points - jquad + 1)
4399 vec_omega_fit_gw_sign_reorder(jquad) = vec_omega_fit_gw_sign(num_fit_points - jquad + 1)
4404 ALLOCATE (coeff_pade(nparam_pade))
4405 ALLOCATE (omega_points_pade(nparam_pade))
4407 CALL get_pade_parameters(sigma_c_gw_reorder, vec_omega_fit_gw_sign_reorder, &
4408 num_fit_points, nparam_pade, omega_points_pade, coeff_pade)
4411 IF ((crossing_search == ri_rpa_g0w0_crossing_bisection) .OR. &
4412 (crossing_search == ri_rpa_g0w0_crossing_newton))
THEN
4413 energy_val = eigenval(n_level_gw_ref) - e_fermi
4414 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, &
4415 coeff_pade, sigma_c_pade)
4416 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, &
4417 coeff_pade, z_value(n_level_gw), m_value(n_level_gw))
4418 level_energ_gw_start = (eigenval_scf(n_level_gw_ref) - &
4419 m_value(n_level_gw)*eigenval(n_level_gw_ref) + &
4420 REAL(sigma_c_pade) + &
4421 vec_sigma_x_minus_vxc_gw(n_level_gw_ref))* &
4425 hedin_shift = 0.0_dp
4426 IF (do_hedin_shift) hedin_shift = real(sigma_c_pade) + &
4427 vec_sigma_x_minus_vxc_gw(n_level_gw_ref) &
4428 - eigenval(n_level_gw_ref) + eigenval_scf(n_level_gw_ref)
4431 IF (
PRESENT(min_level_self_energy) .AND.
PRESENT(max_level_self_energy))
THEN
4432 IF (n_level_gw_ref >= min_level_self_energy .AND. &
4433 n_level_gw_ref <= max_level_self_energy)
THEN
4434 ALLOCATE (vec_sigma_real(ndos))
4435 ALLOCATE (vec_sigma_imag(ndos))
4436 WRITE (string_level,
"(I4)") n_level_gw_ref
4437 string_level = adjustl(string_level)
4446 IF (
PRESENT(ndos))
THEN
4449 cpassert(.NOT. do_hedin_shift)
4450 logger => cp_get_default_logger()
4451 IF (logger%para_env%is_source())
THEN
4452 iunit = cp_logger_get_default_unit_nr()
4457 omega_dos = dos_lower_bound + real(idos - 1, kind=dp)*dos_precision
4458 omega_dos_pade_eval = omega_dos - e_fermi
4459 CALL evaluate_pade_function(omega_dos_pade_eval, nparam_pade, omega_points_pade, &
4460 coeff_pade, sigma_c_pade)
4462 IF (n_level_gw_ref >= min_level_self_energy .AND. &
4463 n_level_gw_ref <= max_level_self_energy .AND. iunit > 0)
THEN
4465 vec_sigma_real(idos) = (real(sigma_c_pade))
4466 vec_sigma_imag(idos) = (aimag(sigma_c_pade))
4470 IF (n_level_gw_ref >= dos_min .AND. &
4471 (n_level_gw_ref <= dos_max .OR. dos_max == 0))
THEN
4472 vec_gw_dos(idos) = vec_gw_dos(idos) + &
4473 (abs(aimag(sigma_c_pade)) + dos_eta) &
4475 (omega_dos - eigenval_scf(n_level_gw_ref) - &
4476 (real(sigma_c_pade) + vec_sigma_x_minus_vxc_gw(n_level_gw_ref)) &
4478 + (abs(aimag(sigma_c_pade)) + dos_eta)**2 &
4486 IF (
PRESENT(min_level_self_energy) .AND.
PRESENT(max_level_self_energy))
THEN
4487 logger => cp_get_default_logger()
4488 IF (logger%para_env%is_source())
THEN
4489 iunit = cp_logger_get_default_unit_nr()
4493 IF (n_level_gw_ref >= min_level_self_energy .AND. &
4494 n_level_gw_ref <= max_level_self_energy .AND. iunit > 0)
THEN
4496 CALL open_file(
'self_energy_re_'//trim(string_level)//
'.dat', unit_number=iunit, &
4497 file_status=
"UNKNOWN", file_action=
"WRITE")
4499 omega_dos = dos_lower_bound + real(idos - 1, kind=dp)*dos_precision
4500 WRITE (iunit,
'(F17.10, F17.10)') omega_dos*evolt, vec_sigma_real(idos)*evolt
4503 CALL close_file(iunit)
4505 CALL open_file(
'self_energy_im_'//trim(string_level)//
'.dat', unit_number=iunit, &
4506 file_status=
"UNKNOWN", file_action=
"WRITE")
4508 omega_dos = dos_lower_bound + real(idos - 1, kind=dp)*dos_precision
4509 WRITE (iunit,
'(F17.10, F17.10)') omega_dos*evolt, vec_sigma_imag(idos)*evolt
4512 CALL close_file(iunit)
4514 DEALLOCATE (vec_sigma_real)
4515 DEALLOCATE (vec_sigma_imag)
4520 SELECT CASE (crossing_search)
4521 CASE (ri_rpa_g0w0_crossing_z_shot)
4523 cpassert(.NOT. do_hedin_shift)
4524 energy_val = eigenval(n_level_gw_ref) - e_fermi
4525 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, &
4526 coeff_pade, sigma_c_pade)
4527 vec_gw_energ(n_level_gw) = real(sigma_c_pade)
4529 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, &
4530 coeff_pade, z_value(n_level_gw), m_value(n_level_gw))
4532 CASE (ri_rpa_g0w0_crossing_bisection)
4533 CALL get_sigma_c_bisection_pade(vec_gw_energ(n_level_gw), eigenval_scf(n_level_gw_ref), &
4534 vec_sigma_x_minus_vxc_gw(n_level_gw_ref), e_fermi, &
4535 nparam_pade, omega_points_pade, coeff_pade, &
4536 level_energ_gw_start, hedin_shift)
4537 z_value(n_level_gw) = 1.0_dp
4538 m_value(n_level_gw) = 0.0_dp
4540 CASE (ri_rpa_g0w0_crossing_newton)
4541 CALL get_sigma_c_newton_pade(vec_gw_energ(n_level_gw), eigenval_scf(n_level_gw_ref), &
4542 vec_sigma_x_minus_vxc_gw(n_level_gw_ref), e_fermi, &
4543 nparam_pade, omega_points_pade, coeff_pade, &
4544 level_energ_gw_start, hedin_shift)
4545 z_value(n_level_gw) = 1.0_dp
4546 m_value(n_level_gw) = 0.0_dp
4549 cpabort(
"Only Z_SHOT, NEWTON, and BISECTION crossing search implemented.")
4552 IF (print_self_energy)
THEN
4554 IF (count_ev_sc_gw == 1)
THEN
4556 IF (n_level_gw_ref < 10)
THEN
4557 WRITE (filename,
"(A26,I1)")
"G0W0_self_energy_level_000", n_level_gw_ref
4558 ELSE IF (n_level_gw_ref < 100)
THEN
4559 WRITE (filename,
"(A25,I2)")
"G0W0_self_energy_level_00", n_level_gw_ref
4560 ELSE IF (n_level_gw_ref < 1000)
THEN
4561 WRITE (filename,
"(A24,I3)")
"G0W0_self_energy_level_0", n_level_gw_ref
4563 WRITE (filename,
"(A23,I4)")
"G0W0_self_energy_level_", n_level_gw_ref
4568 IF (n_level_gw_ref < 10)
THEN
4569 WRITE (filename,
"(A11,I1,A22,I1)")
"evGW_cycle_", count_ev_sc_gw, &
4570 "_self_energy_level_000", n_level_gw_ref
4571 ELSE IF (n_level_gw_ref < 100)
THEN
4572 WRITE (filename,
"(A11,I1,A21,I2)")
"evGW_cycle_", count_ev_sc_gw, &
4573 "_self_energy_level_00", n_level_gw_ref
4574 ELSE IF (n_level_gw_ref < 1000)
THEN
4575 WRITE (filename,
"(A11,I1,A20,I3)")
"evGW_cycle_", count_ev_sc_gw, &
4576 "_self_energy_level_0", n_level_gw_ref
4578 WRITE (filename,
"(A11,I1,A19,I4)")
"evGW_cycle_", count_ev_sc_gw, &
4579 "_self_energy_level_", n_level_gw_ref
4584 logger => cp_get_default_logger()
4585 IF (logger%para_env%is_source())
THEN
4586 iunit = cp_logger_get_default_unit_nr()
4590 CALL open_file(trim(filename), unit_number=iunit, file_status=
"UNKNOWN", file_action=
"WRITE")
4594 WRITE (iunit,
"(2A42)")
" omega (eV) Sigma(omega) (eV) ", &
4595 " omega - e_n^DFT - Sigma_n^x - v_n^xc (eV)"
4597 DO i_omega = 0, num_omega
4599 omega = -50.0_dp/evolt + real(i_omega, kind=dp)/real(num_omega, kind=dp)*100.0_dp/evolt
4601 CALL evaluate_pade_function(omega - e_fermi, nparam_pade, omega_points_pade, &
4602 coeff_pade, sigma_c_pade)
4604 WRITE (iunit,
"(F12.2,2F17.5)") omega*evolt, real(sigma_c_pade)*evolt, &
4605 (omega - eigenval_scf(n_level_gw_ref) - vec_sigma_x_minus_vxc_gw(n_level_gw_ref))*evolt
4609 WRITE (iunit,
"(A51,A39)")
" w (eV) Re(Sigma(i*w)) (eV) Im(Sigma(i*w)) (eV) ", &
4610 " Re(Fit(i*w)) (eV) Im(Fit(iw)) (eV)"
4612 DO jquad = 1, num_fit_points
4614 CALL evaluate_pade_function(vec_omega_fit_gw_sign_reorder(jquad), &
4615 nparam_pade, omega_points_pade, &
4616 coeff_pade, sigma_c_pade_im_freq, do_imag_freq=.true.)
4618 WRITE (iunit,
"(F12.2,4F17.5)") vec_omega_fit_gw_sign_reorder(jquad)*evolt, &
4619 REAL(sigma_c_gw_reorder(jquad)*evolt), &
4620 aimag(sigma_c_gw_reorder(jquad)*evolt), &
4621 REAL(sigma_c_pade_im_freq*evolt), &
4622 aimag(sigma_c_pade_im_freq*evolt)
4626 CALL close_file(iunit)
4630 DEALLOCATE (vec_omega_fit_gw_sign)
4631 DEALLOCATE (sigma_c_gw_reorder)
4632 DEALLOCATE (vec_omega_fit_gw_sign_reorder)
4633 DEALLOCATE (coeff_pade, omega_points_pade)
4635 CALL timestop(handle)
4649 PURE SUBROUTINE get_pade_parameters(y, x, num_fit_points, nparam, xpoints, coeff)
4651 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: y
4652 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: x
4653 INTEGER,
INTENT(IN) :: num_fit_points, nparam
4654 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: xpoints, coeff
4656 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: ypoints
4657 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: g_mat
4658 INTEGER :: idat, iparam, nstep
4660 nstep = int(num_fit_points/(nparam - 1))
4662 ALLOCATE (ypoints(nparam))
4665 DO iparam = 1, nparam - 1
4666 xpoints(iparam) = gaussi*x(idat)
4667 ypoints(iparam) = y(idat)
4670 xpoints(nparam) = gaussi*x(num_fit_points)
4671 ypoints(nparam) = y(num_fit_points)
4675 ALLOCATE (g_mat(nparam, nparam))
4676 g_mat(:, 1) = ypoints(:)
4677 DO iparam = 2, nparam
4678 DO idat = iparam, nparam
4679 g_mat(idat, iparam) = (g_mat(iparam - 1, iparam - 1) - g_mat(idat, iparam - 1))/ &
4680 ((xpoints(idat) - xpoints(iparam - 1))*g_mat(idat, iparam - 1))
4684 DO iparam = 1, nparam
4685 coeff(iparam) = g_mat(iparam, iparam)
4688 DEALLOCATE (ypoints)
4691 END SUBROUTINE get_pade_parameters
4702 PURE SUBROUTINE evaluate_pade_function(x_val, nparam, xpoints, coeff, func_val, do_imag_freq)
4704 REAL(kind=dp),
INTENT(IN) :: x_val
4705 INTEGER,
INTENT(IN) :: nparam
4706 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: xpoints, coeff
4707 COMPLEX(KIND=dp),
INTENT(OUT) :: func_val
4708 LOGICAL,
INTENT(IN),
OPTIONAL :: do_imag_freq
4711 LOGICAL :: my_do_imag_freq
4713 my_do_imag_freq = .false.
4714 IF (
PRESENT(do_imag_freq)) my_do_imag_freq = do_imag_freq
4717 DO iparam = nparam, 2, -1
4718 IF (my_do_imag_freq)
THEN
4719 func_val = z_one + coeff(iparam)*(gaussi*x_val - xpoints(iparam - 1))/func_val
4721 func_val = z_one + coeff(iparam)*(x_val*z_one - xpoints(iparam - 1))/func_val
4725 func_val = coeff(1)/func_val
4727 END SUBROUTINE evaluate_pade_function
4738 PURE SUBROUTINE get_z_and_m_value_pade(x_val, nparam, xpoints, coeff, z_value, m_value)
4740 REAL(kind=dp),
INTENT(IN) :: x_val
4741 INTEGER,
INTENT(IN) :: nparam
4742 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: xpoints, coeff
4743 REAL(kind=dp),
INTENT(OUT),
OPTIONAL :: z_value, m_value
4745 COMPLEX(KIND=dp) :: denominator, dev_denominator, &
4746 dev_numerator, dev_val, func_val, &
4752 DO iparam = nparam, 2, -1
4753 numerator = coeff(iparam)*(x_val*z_one - xpoints(iparam - 1))
4754 dev_numerator = coeff(iparam)*z_one
4755 denominator = func_val
4756 dev_denominator = dev_val
4757 dev_val = dev_numerator/denominator - (numerator*dev_denominator)/(denominator**2)
4758 func_val = z_one + coeff(iparam)*(x_val*z_one - xpoints(iparam - 1))/func_val
4761 dev_val = -1.0_dp*coeff(1)/(func_val**2)*dev_val
4762 func_val = coeff(1)/func_val
4764 IF (
PRESENT(z_value))
THEN
4765 z_value = 1.0_dp - real(dev_val)
4766 z_value = 1.0_dp/z_value
4768 IF (
PRESENT(m_value)) m_value = real(dev_val)
4770 END SUBROUTINE get_z_and_m_value_pade
4784 SUBROUTINE get_sigma_c_bisection_pade(gw_energ, Eigenval_scf, Sigma_x_minus_vxc_gw, e_fermi, &
4785 nparam_pade, omega_points_pade, coeff_pade, start_val, &
4788 REAL(kind=dp),
INTENT(OUT) :: gw_energ
4789 REAL(kind=dp),
INTENT(IN) :: eigenval_scf, sigma_x_minus_vxc_gw, &
4791 INTEGER,
INTENT(IN) :: nparam_pade
4792 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: omega_points_pade, coeff_pade
4793 REAL(kind=dp),
INTENT(IN) :: start_val, hedin_shift
4795 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_sigma_c_bisection_pade'
4797 COMPLEX(KIND=dp) :: sigma_c
4798 INTEGER :: handle, icount
4799 REAL(kind=dp) :: delta, energy_val, qp_energy, &
4800 qp_energy_old, threshold
4802 CALL timeset(routinen, handle)
4804 threshold = 1.0e-7_dp
4806 qp_energy = start_val
4807 qp_energy_old = start_val
4811 DO WHILE (abs(delta) > threshold)
4813 qp_energy = qp_energy_old + 0.5_dp*delta
4814 qp_energy_old = qp_energy
4815 energy_val = qp_energy - e_fermi - hedin_shift
4816 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, &
4817 coeff_pade, sigma_c)
4818 qp_energy = eigenval_scf + real(sigma_c) + sigma_x_minus_vxc_gw
4819 delta = qp_energy - qp_energy_old
4821 IF (icount > 500)
EXIT
4824 gw_energ = real(sigma_c)
4826 CALL timestop(handle)
4828 END SUBROUTINE get_sigma_c_bisection_pade
4842 SUBROUTINE get_sigma_c_newton_pade(gw_energ, Eigenval_scf, Sigma_x_minus_vxc_gw, e_fermi, &
4843 nparam_pade, omega_points_pade, coeff_pade, start_val, &
4846 REAL(kind=dp),
INTENT(OUT) :: gw_energ
4847 REAL(kind=dp),
INTENT(IN) :: eigenval_scf, sigma_x_minus_vxc_gw, &
4849 INTEGER,
INTENT(IN) :: nparam_pade
4850 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: omega_points_pade, coeff_pade
4851 REAL(kind=dp),
INTENT(IN) :: start_val, hedin_shift
4853 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_sigma_c_newton_pade'
4855 COMPLEX(KIND=dp) :: sigma_c
4856 INTEGER :: handle, icount
4857 REAL(kind=dp) :: delta, energy_val, m_value, qp_energy, &
4858 qp_energy_old, threshold
4860 CALL timeset(routinen, handle)
4862 threshold = 1.0e-7_dp
4864 qp_energy = start_val
4865 qp_energy_old = start_val
4869 DO WHILE (abs(delta) > threshold)
4871 energy_val = qp_energy - e_fermi - hedin_shift
4872 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, &
4873 coeff_pade, sigma_c)
4875 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, &
4876 coeff_pade, m_value=m_value)
4877 qp_energy_old = qp_energy
4878 qp_energy = qp_energy - (eigenval_scf + sigma_x_minus_vxc_gw + real(sigma_c) - qp_energy)/ &
4880 delta = qp_energy - qp_energy_old
4882 IF (icount > 500)
EXIT
4885 gw_energ = real(sigma_c)
4887 CALL timestop(handle)
4889 END SUBROUTINE get_sigma_c_newton_pade
4918 SUBROUTINE print_and_update_for_ev_sc(vec_gw_energ, &
4919 z_value, m_value, vec_Sigma_x_minus_vxc_gw, Eigenval, &
4920 Eigenval_last, Eigenval_scf, &
4921 gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, &
4922 crossing_search, homo, unit_nr, count_ev_sc_GW, count_sc_GW0, &
4923 ikp, nkp_self_energy, kpoints, ispin, E_VBM_GW, E_CBM_GW, &
4924 E_VBM_SCF, E_CBM_SCF)
4926 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_gw_energ, z_value, m_value
4927 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: vec_sigma_x_minus_vxc_gw, eigenval, &
4928 eigenval_last, eigenval_scf
4929 INTEGER,
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, crossing_search, &
4930 homo, unit_nr, count_ev_sc_gw, count_sc_gw0, ikp, nkp_self_energy
4931 TYPE(kpoint_type),
INTENT(IN),
POINTER :: kpoints
4932 INTEGER,
INTENT(IN) :: ispin
4933 REAL(kind=dp),
INTENT(INOUT),
OPTIONAL :: e_vbm_gw, e_cbm_gw, e_vbm_scf, e_cbm_scf
4935 CHARACTER(LEN=*),
PARAMETER :: routinen =
'print_and_update_for_ev_sc'
4937 CHARACTER(4) :: occ_virt
4938 INTEGER :: handle, n_level_gw, n_level_gw_ref
4939 LOGICAL :: do_alpha, do_beta, do_closed_shell, &
4940 do_kpoints, is_energy_okay
4941 REAL(kind=dp) :: e_gap_gw, e_homo_gw, e_homo_scf, &
4942 e_lumo_gw, e_lumo_scf, new_energy
4944 CALL timeset(routinen, handle)
4946 do_alpha = (ispin == 1)
4947 do_beta = (ispin == 2)
4948 do_closed_shell = .NOT. (do_alpha .OR. do_beta)
4949 do_kpoints = (nkp_self_energy > 1)
4951 eigenval_last(:) = eigenval(:)
4953 IF (unit_nr > 0)
THEN
4955 IF (count_ev_sc_gw == 1 .AND. count_sc_gw0 == 1 .AND. ikp == 1)
THEN
4957 WRITE (unit_nr, *)
' '
4959 IF (do_alpha .OR. do_closed_shell)
THEN
4960 WRITE (unit_nr, *)
' '
4961 WRITE (unit_nr,
'(T3,A)')
'******************************************************************************'
4962 WRITE (unit_nr,
'(T3,A)')
'** **'
4963 WRITE (unit_nr,
'(T3,A)')
'** GW QUASIPARTICLE ENERGIES **'
4964 WRITE (unit_nr,
'(T3,A)')
'** **'
4965 WRITE (unit_nr,
'(T3,A)')
'******************************************************************************'
4966 WRITE (unit_nr,
'(T3,A)')
' '
4967 WRITE (unit_nr,
'(T3,A)')
' '
4968 WRITE (unit_nr,
'(T3,A)')
'The GW quasiparticle energies are calculated according to: '
4970 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot)
THEN
4971 WRITE (unit_nr,
'(T3,A)')
'E_GW = E_SCF + Z * ( Sigc(E_SCF) + Sigx - vxc )'
4973 WRITE (unit_nr,
'(T3,A)')
' '
4974 WRITE (unit_nr,
'(T3,A)')
' E_GW = E_SCF + Sigc(E_GW) + Sigx - vxc '
4975 WRITE (unit_nr,
'(T3,A)')
' '
4976 WRITE (unit_nr,
'(T3,A)')
'Upper equation is solved self-consistently for E_GW, see Eq. (12) in J. Phys.'
4977 WRITE (unit_nr,
'(T3,A)')
'Chem. Lett. 9, 306 (2018), doi: 10.1021/acs.jpclett.7b02740'
4979 WRITE (unit_nr, *)
' '
4980 WRITE (unit_nr, *)
' '
4981 WRITE (unit_nr,
'(T3,A)')
'------------'
4982 WRITE (unit_nr,
'(T3,A)')
'G0W0 results'
4983 WRITE (unit_nr,
'(T3,A)')
'------------'
4987 IF (.NOT. do_kpoints)
THEN
4989 WRITE (unit_nr, *)
' '
4990 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
4991 WRITE (unit_nr,
'(T3,A)')
'GW quasiparticle energies of alpha spins'
4992 WRITE (unit_nr,
'(T3,A)')
'----------------------------------------'
4993 ELSE IF (do_beta)
THEN
4994 WRITE (unit_nr, *)
' '
4995 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
4996 WRITE (unit_nr,
'(T3,A)')
'GW quasiparticle energies of beta spins'
4997 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
5003 IF (count_ev_sc_gw > 1)
THEN
5004 WRITE (unit_nr, *)
' '
5005 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
5006 WRITE (unit_nr,
'(T3,A,I4)')
'Eigenvalue-selfconsistency cycle: ', count_ev_sc_gw
5007 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
5010 IF (count_sc_gw0 > 1)
THEN
5011 WRITE (unit_nr,
'(T3,A)')
'----------------------------------'
5012 WRITE (unit_nr,
'(T3,A,I4)')
'scGW0 selfconsistency cycle: ', count_sc_gw0
5013 WRITE (unit_nr,
'(T3,A)')
'----------------------------------'
5016 IF (do_kpoints)
THEN
5017 WRITE (unit_nr, *)
' '
5018 WRITE (unit_nr,
'(T3,A7,I3,A3,I3,A8,3F7.3,A12,3F7.3)')
'Kpoint ', ikp,
' /', nkp_self_energy, &
5019 ' xkp =', kpoints%xkp(1, ikp), kpoints%xkp(2, ikp), kpoints%xkp(3, ikp), &
5020 ' and xkp =', -kpoints%xkp(1, ikp), -kpoints%xkp(2, ikp), -kpoints%xkp(3, ikp)
5021 WRITE (unit_nr,
'(T3,A72)')
'(Relative Brillouin zone size: [-0.5, 0.5] x [-0.5, 0.5] x [-0.5, 0.5])'
5022 WRITE (unit_nr, *)
' '
5024 WRITE (unit_nr,
'(T3,A)')
'GW quasiparticle energies of alpha spins:'
5025 ELSE IF (do_beta)
THEN
5026 WRITE (unit_nr,
'(T3,A)')
'GW quasiparticle energies of beta spins:'
5032 DO n_level_gw = 1, gw_corr_lev_tot
5034 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
5036 new_energy = (eigenval_scf(n_level_gw_ref) - &
5037 m_value(n_level_gw)*eigenval(n_level_gw_ref) + &
5038 vec_gw_energ(n_level_gw) + &
5039 vec_sigma_x_minus_vxc_gw(n_level_gw_ref))* &
5042 is_energy_okay = .true.
5044 IF (n_level_gw_ref > homo .AND. new_energy < eigenval(homo))
THEN
5045 is_energy_okay = .false.
5048 IF (is_energy_okay)
THEN
5049 eigenval(n_level_gw_ref) = new_energy
5054 IF (unit_nr > 0)
THEN
5055 WRITE (unit_nr,
'(T3,A)')
' '
5056 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot)
THEN
5057 WRITE (unit_nr,
'(T13,2A)')
'MO E_SCF (eV) Sigc (eV) Sigx-vxc (eV) Z E_GW (eV)'
5059 WRITE (unit_nr,
'(T3,2A)')
'Molecular orbital E_SCF (eV) Sigc (eV) Sigx-vxc (eV) E_GW (eV)'
5063 DO n_level_gw = 1, gw_corr_lev_tot
5064 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
5065 IF (n_level_gw <= gw_corr_lev_occ)
THEN
5071 IF (unit_nr > 0)
THEN
5072 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot)
THEN
5073 WRITE (unit_nr,
'(T3,I4,3A,5F13.4)') &
5074 n_level_gw_ref,
' ( ', occ_virt,
') ', &
5075 eigenval_last(n_level_gw_ref)*evolt, &
5076 vec_gw_energ(n_level_gw)*evolt, &
5077 vec_sigma_x_minus_vxc_gw(n_level_gw_ref)*evolt, &
5078 z_value(n_level_gw), &
5079 eigenval(n_level_gw_ref)*evolt
5081 WRITE (unit_nr,
'(T3,I4,3A,4F16.4)') &
5082 n_level_gw_ref,
' ( ', occ_virt,
') ', &
5083 eigenval_last(n_level_gw_ref)*evolt, &
5084 vec_gw_energ(n_level_gw)*evolt, &
5085 vec_sigma_x_minus_vxc_gw(n_level_gw_ref)*evolt, &
5086 eigenval(n_level_gw_ref)*evolt
5091 e_homo_scf = maxval(eigenval_last(homo - gw_corr_lev_occ + 1:homo))
5092 e_lumo_scf = minval(eigenval_last(homo + 1:homo + gw_corr_lev_virt))
5094 e_homo_gw = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo))
5095 e_lumo_gw = minval(eigenval(homo + 1:homo + gw_corr_lev_virt))
5096 e_gap_gw = e_lumo_gw - e_homo_gw
5098 IF (
PRESENT(e_vbm_scf) .AND.
PRESENT(e_cbm_scf) .AND. &
5099 PRESENT(e_vbm_gw) .AND.
PRESENT(e_cbm_gw))
THEN
5100 IF (e_homo_scf > e_vbm_scf) e_vbm_scf = e_homo_scf
5101 IF (e_lumo_scf < e_cbm_scf) e_cbm_scf = e_lumo_scf
5102 IF (e_homo_gw > e_vbm_gw) e_vbm_gw = e_homo_gw
5103 IF (e_lumo_gw < e_cbm_gw) e_cbm_gw = e_lumo_gw
5106 IF (unit_nr > 0)
THEN
5108 IF (do_kpoints)
THEN
5109 IF (do_closed_shell)
THEN
5110 WRITE (unit_nr,
'(T3,A)')
' '
5111 WRITE (unit_nr,
'(T3,A,F42.4)')
'GW direct gap at current kpoint (eV)', e_gap_gw*evolt
5112 ELSE IF (do_alpha)
THEN
5113 WRITE (unit_nr,
'(T3,A)')
' '
5114 WRITE (unit_nr,
'(T3,A,F36.4)')
'Alpha GW direct gap at current kpoint (eV)', &
5116 ELSE IF (do_beta)
THEN
5117 WRITE (unit_nr,
'(T3,A)')
' '
5118 WRITE (unit_nr,
'(T3,A,F37.4)')
'Beta GW direct gap at current kpoint (eV)', &
5122 IF (do_closed_shell)
THEN
5123 WRITE (unit_nr,
'(T3,A)')
' '
5124 IF (count_ev_sc_gw > 1)
THEN
5125 WRITE (unit_nr,
'(T3,A,I3,A,F39.4)')
'HOMO-LUMO gap in evGW iteration', &
5126 count_ev_sc_gw,
' (eV)', e_gap_gw*evolt
5127 ELSE IF (count_sc_gw0 > 1)
THEN
5128 WRITE (unit_nr,
'(T3,A,I3,A,F38.4)')
'HOMO-LUMO gap in evGW0 iteration', &
5129 count_sc_gw0,
' (eV)', e_gap_gw*evolt
5131 WRITE (unit_nr,
'(T3,A,F55.4)')
'G0W0 HOMO-LUMO gap (eV)', e_gap_gw*evolt
5133 ELSE IF (do_alpha)
THEN
5134 WRITE (unit_nr,
'(T3,A)')
' '
5135 WRITE (unit_nr,
'(T3,A,F51.4)')
'Alpha GW HOMO-LUMO gap (eV)', e_gap_gw*evolt
5136 ELSE IF (do_beta)
THEN
5137 WRITE (unit_nr,
'(T3,A)')
' '
5138 WRITE (unit_nr,
'(T3,A,F52.4)')
'Beta GW HOMO-LUMO gap (eV)', e_gap_gw*evolt
5143 IF (unit_nr > 0)
THEN
5144 WRITE (unit_nr, *)
' '
5145 WRITE (unit_nr,
'(T3,A)')
'------------------------------------------------------------------------------'
5148 CALL timestop(handle)
5150 END SUBROUTINE print_and_update_for_ev_sc
5161 PURE SUBROUTINE shift_unshifted_levels(Eigenval, Eigenval_last, gw_corr_lev_occ, gw_corr_lev_virt, &
5164 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: eigenval, eigenval_last
5165 INTEGER,
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo, &
5168 INTEGER :: n_level_gw, n_level_gw_ref
5169 REAL(kind=dp) :: eigen_diff
5173 IF (gw_corr_lev_occ < homo .AND. gw_corr_lev_occ > 0)
THEN
5178 DO n_level_gw = 1, gw_corr_lev_occ
5179 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
5180 eigen_diff = eigen_diff + eigenval(n_level_gw_ref) - eigenval_last(n_level_gw_ref)
5182 eigen_diff = eigen_diff/gw_corr_lev_occ
5185 DO n_level_gw = 1, homo - gw_corr_lev_occ
5186 eigenval(n_level_gw) = eigenval(n_level_gw) + eigen_diff
5192 IF (gw_corr_lev_virt < nmo - homo .AND. gw_corr_lev_virt > 0)
THEN
5196 DO n_level_gw = 1, gw_corr_lev_virt
5197 n_level_gw_ref = n_level_gw + homo
5198 eigen_diff = eigen_diff + eigenval(n_level_gw_ref) - eigenval_last(n_level_gw_ref)
5200 eigen_diff = eigen_diff/gw_corr_lev_virt
5203 DO n_level_gw = homo + gw_corr_lev_virt + 1, nmo
5204 eigenval(n_level_gw) = eigenval(n_level_gw) + eigen_diff
5209 END SUBROUTINE shift_unshifted_levels
5226 SUBROUTINE calc_mat_n(N_ij, Lambda, Sigma_c, vec_omega_fit_gw, i, j, &
5227 num_poles, num_fit_points, n_level_gw, h)
5228 REAL(kind=dp),
INTENT(OUT) :: n_ij
5229 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:), &
5230 INTENT(IN) :: lambda
5231 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: sigma_c
5232 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
5233 INTENT(IN) :: vec_omega_fit_gw
5234 INTEGER,
INTENT(IN) :: i, j, num_poles, num_fit_points, &
5236 REAL(kind=dp),
INTENT(IN) :: h
5238 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_mat_N'
5240 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: lambda_tmp
5241 INTEGER :: handle, num_var
5242 REAL(kind=dp) :: chi2, chi2_sum
5244 CALL timeset(routinen, handle)
5246 num_var = 2*num_poles + 1
5247 ALLOCATE (lambda_tmp(num_var))
5252 lambda_tmp(:) = lambda(:)
5253 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5254 num_fit_points, n_level_gw)
5257 lambda_tmp(:) = lambda(:)
5258 IF (
modulo(i, 2) == 0)
THEN
5259 lambda_tmp(i/2) = lambda_tmp(i/2) + h*z_one
5261 lambda_tmp((i + 1)/2) = lambda_tmp((i + 1)/2) + h*gaussi
5263 IF (
modulo(j, 2) == 0)
THEN
5264 lambda_tmp(j/2) = lambda_tmp(j/2) + h*z_one
5266 lambda_tmp((j + 1)/2) = lambda_tmp((j + 1)/2) + h*gaussi
5268 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5269 num_fit_points, n_level_gw)
5270 chi2_sum = chi2_sum + chi2
5272 IF (
modulo(i, 2) == 0)
THEN
5273 lambda_tmp(i/2) = lambda_tmp(i/2) - 2.0_dp*h*z_one
5275 lambda_tmp((i + 1)/2) = lambda_tmp((i + 1)/2) - 2.0_dp*h*gaussi
5277 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5278 num_fit_points, n_level_gw)
5279 chi2_sum = chi2_sum - chi2
5281 IF (
modulo(j, 2) == 0)
THEN
5282 lambda_tmp(j/2) = lambda_tmp(j/2) - 2.0_dp*h*z_one
5284 lambda_tmp((j + 1)/2) = lambda_tmp((j + 1)/2) - 2.0_dp*h*gaussi
5286 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5287 num_fit_points, n_level_gw)
5288 chi2_sum = chi2_sum + chi2
5290 IF (
modulo(i, 2) == 0)
THEN
5291 lambda_tmp(i/2) = lambda_tmp(i/2) + 2.0_dp*h*z_one
5293 lambda_tmp((i + 1)/2) = lambda_tmp((i + 1)/2) + 2.0_dp*h*gaussi
5295 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5296 num_fit_points, n_level_gw)
5297 chi2_sum = chi2_sum - chi2
5300 n_ij = 1.0_dp/2.0_dp*chi2_sum/(4.0_dp*h*h)
5302 DEALLOCATE (lambda_tmp)
5304 CALL timestop(handle)
5306 END SUBROUTINE calc_mat_n
5318 PURE SUBROUTINE calc_chi2(chi2, Lambda, Sigma_c, vec_omega_fit_gw, num_poles, &
5319 num_fit_points, n_level_gw)
5320 REAL(kind=dp),
INTENT(OUT) :: chi2
5321 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: lambda
5322 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: sigma_c
5323 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_omega_fit_gw
5324 INTEGER,
INTENT(IN) :: num_poles, num_fit_points, n_level_gw
5326 COMPLEX(KIND=dp) :: func_val
5327 INTEGER :: iii, jjj, kkk
5330 DO kkk = 1, num_fit_points
5331 func_val = lambda(1)
5332 DO iii = 1, num_poles
5335 func_val = func_val + lambda(jjj)/(gaussi*vec_omega_fit_gw(kkk) - lambda(jjj + 1))
5337 chi2 = chi2 + (abs(sigma_c(n_level_gw, kkk) - func_val))**2
5340 END SUBROUTINE calc_chi2
5394 SUBROUTINE compute_self_energy_cubic_gw(num_integ_points, nmo, tau_tj, tj, &
5395 matrix_s, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
5396 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, &
5397 fm_scaled_dm_virt_tau, Eigenval, eps_filter, &
5398 e_fermi, fm_mat_W, &
5399 gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, &
5400 count_ev_sc_GW, count_sc_GW0, &
5401 t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
5402 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
5403 mat_W, mat_MinvVMinv, mat_dm, &
5404 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_Sigma_c_gw, &
5405 do_periodic, num_points_corr, delta_corr, qs_env, para_env, para_env_RPA, &
5406 mp2_env, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
5407 first_cycle_periodic_correction, kpoints, num_fit_points, fm_mo_coeff, &
5408 do_ri_Sigma_x, vec_Sigma_x_gw, unit_nr, ispin)
5409 INTEGER,
INTENT(IN) :: num_integ_points, nmo
5410 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
5411 INTENT(IN) :: tau_tj, tj
5412 TYPE(dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_s
5413 TYPE(cp_fm_type),
INTENT(IN) :: fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
5414 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau
5415 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: eigenval
5416 REAL(kind=dp),
INTENT(IN) :: eps_filter
5417 REAL(kind=dp),
INTENT(INOUT) :: e_fermi
5418 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_w
5419 INTEGER,
INTENT(IN) :: gw_corr_lev_tot, gw_corr_lev_occ, &
5420 gw_corr_lev_virt, homo, &
5421 count_ev_sc_gw, count_sc_gw0
5422 TYPE(dbt_type) :: t_3c_overl_int_ao_mo
5423 TYPE(hfx_compression_type) :: t_3c_o_mo_compressed
5424 INTEGER,
DIMENSION(:, :) :: t_3c_o_mo_ind
5425 TYPE(dbt_type) :: t_3c_overl_int_gw_ri, &
5426 t_3c_overl_int_gw_ao
5427 TYPE(dbcsr_type),
INTENT(INOUT),
TARGET :: mat_w
5428 TYPE(dbcsr_p_type) :: mat_minvvminv, mat_dm
5429 REAL(kind=dp),
DIMENSION(:, :),
INTENT(IN) :: weights_cos_tf_t_to_w, &
5430 weights_sin_tf_t_to_w
5431 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(OUT) :: vec_sigma_c_gw
5432 LOGICAL,
INTENT(IN) :: do_periodic
5433 INTEGER,
INTENT(IN) :: num_points_corr
5434 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
5435 INTENT(INOUT) :: delta_corr
5436 TYPE(qs_environment_type),
POINTER :: qs_env
5437 TYPE(mp_para_env_type),
POINTER :: para_env, para_env_rpa
5438 TYPE(mp2_type),
INTENT(INOUT) :: mp2_env
5439 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_re_mo_mo, &
5440 matrix_berry_im_mo_mo
5441 LOGICAL,
INTENT(INOUT) :: first_cycle_periodic_correction
5442 TYPE(kpoint_type),
POINTER :: kpoints
5443 INTEGER,
INTENT(IN) :: num_fit_points
5444 TYPE(cp_fm_type),
INTENT(IN) :: fm_mo_coeff
5445 LOGICAL,
INTENT(IN) :: do_ri_sigma_x
5446 REAL(kind=dp),
DIMENSION(:, :),
INTENT(INOUT) :: vec_sigma_x_gw
5447 INTEGER,
INTENT(IN) :: unit_nr, ispin
5449 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_self_energy_cubic_gw'
5451 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: delta_corr_omega
5452 INTEGER :: gw_lev_end, gw_lev_start, handle, handle3, i, iblk_mo, iquad, jquad, mo_end, &
5453 mo_start, n_level_gw, n_level_gw_ref, nblk_mo, unit_nr_prv
5454 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_range_mo, dist1, dist2, mo_bsizes, &
5455 mo_offsets, sizes_ao, sizes_ri
5456 INTEGER,
DIMENSION(2) :: mo_bounds, pdims_2d
5457 LOGICAL :: memory_info
5458 REAL(kind=dp) :: ext_scaling, omega, omega_i, omega_sign, &
5459 sign_occ_virt, t_i_clenshaw, tau, &
5460 weight_cos, weight_i, weight_sin
5461 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: vec_sigma_c_gw_cos_omega, &
5462 vec_sigma_c_gw_cos_tau, vec_sigma_c_gw_neg_tau, vec_sigma_c_gw_pos_tau, &
5463 vec_sigma_c_gw_sin_omega, vec_sigma_c_gw_sin_tau
5464 TYPE(dbcsr_type),
TARGET :: mat_greens_fct_occ, mat_greens_fct_virt
5465 TYPE(dbt_pgrid_type) :: pgrid_2d
5466 TYPE(dbt_type) :: t_3c_ctr_ao, t_3c_ctr_ri, t_ao_tmp, &
5467 t_dm, t_greens_fct_occ, &
5468 t_greens_fct_virt, t_ri_tmp, &
5471 CALL timeset(routinen, handle)
5473 CALL decompress_tensor(t_3c_overl_int_ao_mo, t_3c_o_mo_ind, t_3c_o_mo_compressed, &
5474 mp2_env%ri_rpa_im_time%eps_compress)
5476 CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_gw_ri)
5477 CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_gw_ao, order=[2, 1, 3], move_data=.true.)
5479 memory_info = mp2_env%ri_rpa_im_time%memory_info
5480 IF (memory_info)
THEN
5481 unit_nr_prv = unit_nr
5486 mo_start = homo - gw_corr_lev_occ + 1
5487 mo_end = homo + gw_corr_lev_virt
5488 cpassert(mo_end - mo_start + 1 == gw_corr_lev_tot)
5490 vec_sigma_c_gw = z_zero
5491 ALLOCATE (vec_sigma_c_gw_pos_tau(gw_corr_lev_tot, num_integ_points))
5492 vec_sigma_c_gw_pos_tau = 0.0_dp
5493 ALLOCATE (vec_sigma_c_gw_neg_tau(gw_corr_lev_tot, num_integ_points))
5494 vec_sigma_c_gw_neg_tau = 0.0_dp
5495 ALLOCATE (vec_sigma_c_gw_cos_tau(gw_corr_lev_tot, num_integ_points))
5496 vec_sigma_c_gw_cos_tau = 0.0_dp
5497 ALLOCATE (vec_sigma_c_gw_sin_tau(gw_corr_lev_tot, num_integ_points))
5498 vec_sigma_c_gw_sin_tau = 0.0_dp
5500 ALLOCATE (vec_sigma_c_gw_cos_omega(gw_corr_lev_tot, num_integ_points))
5501 vec_sigma_c_gw_cos_omega = 0.0_dp
5502 ALLOCATE (vec_sigma_c_gw_sin_omega(gw_corr_lev_tot, num_integ_points))
5503 vec_sigma_c_gw_sin_omega = 0.0_dp
5505 ALLOCATE (delta_corr_omega(1 + homo - gw_corr_lev_occ:homo + gw_corr_lev_virt, num_integ_points))
5506 delta_corr_omega(:, :) = z_zero
5508 CALL dbcsr_create(matrix=mat_greens_fct_occ, &
5509 template=matrix_s(1)%matrix, &
5510 matrix_type=dbcsr_type_no_symmetry)
5512 CALL dbcsr_create(matrix=mat_greens_fct_virt, &
5513 template=matrix_s(1)%matrix, &
5514 matrix_type=dbcsr_type_no_symmetry)
5516 e_fermi = 0.5_dp*(eigenval(homo) + eigenval(homo + 1))
5518 nblk_mo = dbt_nblks_total(t_3c_overl_int_gw_ao, 3)
5519 ALLOCATE (mo_offsets(nblk_mo))
5520 ALLOCATE (mo_bsizes(nblk_mo))
5521 ALLOCATE (batch_range_mo(nblk_mo - 1))
5522 CALL dbt_get_info(t_3c_overl_int_gw_ao, blk_offset_3=mo_offsets, blk_size_3=mo_bsizes)
5525 CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d)
5526 ALLOCATE (sizes_ri(dbt_nblks_total(t_3c_overl_int_gw_ri, 1)))
5527 CALL dbt_get_info(t_3c_overl_int_gw_ri, blk_size_1=sizes_ri)
5529 CALL create_2c_tensor(t_w, dist1, dist2, pgrid_2d, sizes_ri, sizes_ri, name=
"(RI|RI)")
5531 DEALLOCATE (dist1, dist2)
5533 CALL dbt_create(mat_w, t_ri_tmp, name=
"(RI|RI)")
5535 CALL dbt_create(t_3c_overl_int_gw_ri, t_3c_ctr_ri)
5536 CALL dbt_create(t_3c_overl_int_gw_ao, t_3c_ctr_ao)
5538 ALLOCATE (sizes_ao(dbt_nblks_total(t_3c_overl_int_gw_ao, 1)))
5539 CALL dbt_get_info(t_3c_overl_int_gw_ao, blk_size_1=sizes_ao)
5540 CALL create_2c_tensor(t_greens_fct_occ, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5541 DEALLOCATE (dist1, dist2)
5542 CALL create_2c_tensor(t_greens_fct_virt, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5543 DEALLOCATE (dist1, dist2)
5545 DO jquad = 1, num_integ_points
5547 CALL compute_greens_function_time(mat_greens_fct_occ, mat_greens_fct_virt, &
5548 fm_mo_coeff_occ, fm_mo_coeff_virt, &
5549 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
5550 fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, eigenval, &
5551 nmo, eps_filter, e_fermi, tau_tj(jquad), para_env)
5553 CALL dbcsr_set(mat_w, 0.0_dp)
5554 CALL copy_fm_to_dbcsr(fm_mat_w(jquad), mat_w, keep_sparsity=.false.)
5556 IF (jquad == 1)
CALL dbt_create(mat_greens_fct_occ, t_ao_tmp, name=
"(AO|AO)")
5558 CALL dbt_copy_matrix_to_tensor(mat_w, t_ri_tmp)
5559 CALL dbt_copy(t_ri_tmp, t_w)
5560 CALL dbt_copy_matrix_to_tensor(mat_greens_fct_occ, t_ao_tmp)
5561 CALL dbt_copy(t_ao_tmp, t_greens_fct_occ)
5562 CALL dbt_copy_matrix_to_tensor(mat_greens_fct_virt, t_ao_tmp)
5563 CALL dbt_copy(t_ao_tmp, t_greens_fct_virt)
5565 batch_range_mo(:) = [(i, i=2, nblk_mo)]
5566 CALL dbt_batched_contract_init(t_3c_overl_int_gw_ao, batch_range_3=batch_range_mo)
5567 CALL dbt_batched_contract_init(t_3c_overl_int_gw_ri, batch_range_3=batch_range_mo)
5568 CALL dbt_batched_contract_init(t_3c_ctr_ao, batch_range_3=batch_range_mo)
5569 CALL dbt_batched_contract_init(t_3c_ctr_ri, batch_range_3=batch_range_mo)
5570 CALL dbt_batched_contract_init(t_w)
5571 CALL dbt_batched_contract_init(t_greens_fct_occ)
5572 CALL dbt_batched_contract_init(t_greens_fct_virt)
5576 DO iblk_mo = 2, nblk_mo - 1
5577 mo_bounds = [mo_offsets(iblk_mo), mo_offsets(iblk_mo) + mo_bsizes(iblk_mo) - 1]
5578 CALL contract_cubic_gw(t_3c_overl_int_gw_ao, t_3c_overl_int_gw_ri, &
5579 t_greens_fct_occ, t_w, [1.0_dp, -1.0_dp], &
5580 mo_bounds, unit_nr_prv, &
5581 t_3c_ctr_ri, t_3c_ctr_ao, calculate_ctr_ri=.true.)
5582 CALL trace_sigma_gw(t_3c_ctr_ao, t_3c_ctr_ri, vec_sigma_c_gw_neg_tau(:, jquad), mo_start, mo_bounds, para_env)
5584 CALL contract_cubic_gw(t_3c_overl_int_gw_ao, t_3c_overl_int_gw_ri, &
5585 t_greens_fct_virt, t_w, [1.0_dp, 1.0_dp], &
5586 mo_bounds, unit_nr_prv, &
5587 t_3c_ctr_ri, t_3c_ctr_ao, calculate_ctr_ri=.false.)
5589 CALL trace_sigma_gw(t_3c_ctr_ao, t_3c_ctr_ri, vec_sigma_c_gw_pos_tau(:, jquad), mo_start, mo_bounds, para_env)
5591 CALL dbt_batched_contract_finalize(t_3c_overl_int_gw_ao)
5592 CALL dbt_batched_contract_finalize(t_3c_overl_int_gw_ri)
5593 CALL dbt_batched_contract_finalize(t_3c_ctr_ao)
5594 CALL dbt_batched_contract_finalize(t_3c_ctr_ri)
5595 CALL dbt_batched_contract_finalize(t_w)
5596 CALL dbt_batched_contract_finalize(t_greens_fct_occ)
5597 CALL dbt_batched_contract_finalize(t_greens_fct_virt)
5599 CALL dbt_clear(t_3c_ctr_ao)
5600 CALL dbt_clear(t_3c_ctr_ri)
5602 vec_sigma_c_gw_cos_tau(:, jquad) = 0.5_dp*(vec_sigma_c_gw_pos_tau(:, jquad) + &
5603 vec_sigma_c_gw_neg_tau(:, jquad))
5605 vec_sigma_c_gw_sin_tau(:, jquad) = 0.5_dp*(vec_sigma_c_gw_pos_tau(:, jquad) - &
5606 vec_sigma_c_gw_neg_tau(:, jquad))
5609 CALL dbt_destroy(t_w)
5611 CALL dbt_destroy(t_greens_fct_occ)
5612 CALL dbt_destroy(t_greens_fct_virt)
5615 DO jquad = 1, num_fit_points
5617 DO iquad = 1, num_integ_points
5621 weight_cos = weights_cos_tf_t_to_w(jquad, iquad)*cos(omega*tau)
5622 weight_sin = weights_sin_tf_t_to_w(jquad, iquad)*sin(omega*tau)
5624 vec_sigma_c_gw_cos_omega(:, jquad) = vec_sigma_c_gw_cos_omega(:, jquad) + &
5625 weight_cos*vec_sigma_c_gw_cos_tau(:, iquad)
5627 vec_sigma_c_gw_sin_omega(:, jquad) = vec_sigma_c_gw_sin_omega(:, jquad) + &
5628 weight_sin*vec_sigma_c_gw_sin_tau(:, iquad)
5636 vec_sigma_c_gw_sin_omega(1:gw_corr_lev_occ, :) = -vec_sigma_c_gw_sin_omega(1:gw_corr_lev_occ, :)
5638 vec_sigma_c_gw(:, 1:num_fit_points, 1) = vec_sigma_c_gw_cos_omega(:, 1:num_fit_points) + &
5639 gaussi*vec_sigma_c_gw_sin_omega(:, 1:num_fit_points)
5641 CALL dbcsr_release(mat_greens_fct_occ)
5642 CALL dbcsr_release(mat_greens_fct_virt)
5644 IF (do_ri_sigma_x .AND. count_ev_sc_gw == 1 .AND. count_sc_gw0 == 1)
THEN
5646 CALL timeset(routinen//
"_RI_HFX_operation_1", handle3)
5649 CALL parallel_gemm(transa=
"N", transb=
"T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, &
5650 matrix_a=fm_mo_coeff_occ, matrix_b=fm_mo_coeff_occ, beta=0.0_dp, &
5651 matrix_c=fm_scaled_dm_occ_tau)
5653 CALL timestop(handle3)
5655 CALL timeset(routinen//
"_RI_HFX_operation_2", handle3)
5657 CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, &
5659 keep_sparsity=.false.)
5661 CALL timestop(handle3)
5663 CALL create_2c_tensor(t_dm, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5664 DEALLOCATE (dist1, dist2)
5666 CALL dbt_copy_matrix_to_tensor(mat_dm%matrix, t_ao_tmp)
5667 CALL dbt_copy(t_ao_tmp, t_dm)
5669 CALL create_2c_tensor(t_sinvvsinv, dist1, dist2, pgrid_2d, sizes_ri, sizes_ri, name=
"(RI|RI)")
5670 DEALLOCATE (dist1, dist2)
5672 CALL dbt_copy_matrix_to_tensor(mat_minvvminv%matrix, t_ri_tmp)
5673 CALL dbt_copy(t_ri_tmp, t_sinvvsinv)
5675 CALL dbt_batched_contract_init(t_3c_overl_int_gw_ao, batch_range_3=batch_range_mo)
5676 CALL dbt_batched_contract_init(t_3c_overl_int_gw_ri, batch_range_3=batch_range_mo)
5677 CALL dbt_batched_contract_init(t_3c_ctr_ri, batch_range_3=batch_range_mo)
5678 CALL dbt_batched_contract_init(t_3c_ctr_ao, batch_range_3=batch_range_mo)
5679 CALL dbt_batched_contract_init(t_dm)
5680 CALL dbt_batched_contract_init(t_sinvvsinv)
5682 DO iblk_mo = 2, nblk_mo - 1
5683 mo_bounds = [mo_offsets(iblk_mo), mo_offsets(iblk_mo) + mo_bsizes(iblk_mo) - 1]
5685 CALL contract_cubic_gw(t_3c_overl_int_gw_ao, t_3c_overl_int_gw_ri, &
5686 t_dm, t_sinvvsinv, [1.0_dp, -1.0_dp], &
5687 mo_bounds, unit_nr_prv, &
5688 t_3c_ctr_ri, t_3c_ctr_ao, calculate_ctr_ri=.true.)
5690 CALL trace_sigma_gw(t_3c_ctr_ao, t_3c_ctr_ri, vec_sigma_x_gw(mo_start:mo_end, 1), mo_start, mo_bounds, para_env)
5692 CALL dbt_batched_contract_finalize(t_3c_overl_int_gw_ao)
5693 CALL dbt_batched_contract_finalize(t_3c_overl_int_gw_ri)
5694 CALL dbt_batched_contract_finalize(t_dm)
5695 CALL dbt_batched_contract_finalize(t_sinvvsinv)
5696 CALL dbt_batched_contract_finalize(t_3c_ctr_ri)
5697 CALL dbt_batched_contract_finalize(t_3c_ctr_ao)
5699 CALL dbt_destroy(t_dm)
5700 CALL dbt_destroy(t_sinvvsinv)
5702 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, 1) = &
5703 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, 1) + &
5704 vec_sigma_x_gw(:, 1)
5708 CALL dbt_pgrid_destroy(pgrid_2d)
5710 CALL dbt_destroy(t_3c_ctr_ri)
5711 CALL dbt_destroy(t_3c_ctr_ao)
5712 CALL dbt_destroy(t_ao_tmp)
5713 CALL dbt_destroy(t_ri_tmp)
5716 IF (do_periodic)
THEN
5718 ext_scaling = 0.2_dp
5721 DO iquad = 1, num_points_corr
5724 t_i_clenshaw = iquad*pi/(2.0_dp*num_points_corr)
5725 omega_i = ext_scaling/tan(t_i_clenshaw)
5727 IF (iquad < num_points_corr)
THEN
5728 weight_i = ext_scaling*pi/(num_points_corr*sin(t_i_clenshaw)**2)
5730 weight_i = ext_scaling*pi/(2.0_dp*num_points_corr*sin(t_i_clenshaw)**2)
5733 CALL calc_periodic_correction(delta_corr, qs_env, para_env, para_env_rpa, &
5734 mp2_env%ri_g0w0%kp_grid, homo, nmo, gw_corr_lev_occ, &
5735 gw_corr_lev_virt, omega_i, fm_mo_coeff, eigenval, &
5736 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
5737 first_cycle_periodic_correction, kpoints, &
5738 mp2_env%ri_g0w0%do_mo_coeff_gamma, &
5739 mp2_env%ri_g0w0%num_kp_grids, mp2_env%ri_g0w0%eps_kpoint, &
5740 mp2_env%ri_g0w0%do_extra_kpoints, &
5741 mp2_env%ri_g0w0%do_aux_bas_gw, mp2_env%ri_g0w0%frac_aux_mos)
5743 DO n_level_gw = 1, gw_corr_lev_tot
5745 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
5747 IF (n_level_gw <= gw_corr_lev_occ)
THEN
5748 sign_occ_virt = -1.0_dp
5750 sign_occ_virt = 1.0_dp
5753 DO jquad = 1, num_integ_points
5755 omega_sign = tj(jquad)*sign_occ_virt
5757 delta_corr_omega(n_level_gw_ref, jquad) = &
5758 delta_corr_omega(n_level_gw_ref, jquad) - &
5759 0.5_dp/pi*weight_i/2.0_dp*delta_corr(n_level_gw_ref)* &
5760 (1.0_dp/(gaussi*(omega_i + omega_sign) + e_fermi - eigenval(n_level_gw_ref)) + &
5761 1.0_dp/(gaussi*(-omega_i + omega_sign) + e_fermi - eigenval(n_level_gw_ref)))
5769 gw_lev_start = 1 + homo - gw_corr_lev_occ
5770 gw_lev_end = homo + gw_corr_lev_virt
5773 vec_sigma_c_gw(1:gw_corr_lev_tot, :, 1) = vec_sigma_c_gw(1:gw_corr_lev_tot, :, 1) + &
5774 delta_corr_omega(gw_lev_start:gw_lev_end, 1:num_fit_points)
5778 DEALLOCATE (vec_sigma_c_gw_pos_tau)
5779 DEALLOCATE (vec_sigma_c_gw_neg_tau)
5780 DEALLOCATE (vec_sigma_c_gw_cos_tau)
5781 DEALLOCATE (vec_sigma_c_gw_sin_tau)
5782 DEALLOCATE (vec_sigma_c_gw_cos_omega)
5783 DEALLOCATE (vec_sigma_c_gw_sin_omega)
5784 DEALLOCATE (delta_corr_omega)
5786 CALL timestop(handle)
5788 END SUBROUTINE compute_self_energy_cubic_gw
5827 SUBROUTINE compute_self_energy_cubic_gw_kpoints(num_integ_points, tau_tj, tj, &
5828 matrix_s, Eigenval, e_fermi, fm_mat_W, &
5829 gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, &
5830 count_ev_sc_GW, count_sc_GW0, &
5831 t_3c_O, t_3c_M, t_3c_O_compressed, t_3c_O_ind, &
5832 mat_W, mat_MinvVMinv, &
5833 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_Sigma_c_gw, &
5835 mp2_env, num_fit_points, fm_mo_coeff, &
5836 do_ri_Sigma_x, vec_Sigma_x_gw, unit_nr, nspins, &
5837 starts_array_mc, ends_array_mc, eps_filter)
5839 INTEGER,
INTENT(IN) :: num_integ_points
5840 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
5841 INTENT(IN) :: tau_tj, tj
5842 TYPE(dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_s
5843 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(IN) :: eigenval
5844 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: e_fermi
5845 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_w
5846 INTEGER,
INTENT(IN) :: gw_corr_lev_tot
5847 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
5848 INTEGER,
INTENT(IN) :: count_ev_sc_gw, count_sc_gw0
5849 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_o
5850 TYPE(dbt_type) :: t_3c_m
5851 TYPE(hfx_compression_type),
ALLOCATABLE, &
5852 DIMENSION(:, :, :) :: t_3c_o_compressed
5853 TYPE(block_ind_type),
ALLOCATABLE, &
5854 DIMENSION(:, :, :),
INTENT(INOUT) :: t_3c_o_ind
5855 TYPE(dbcsr_type),
INTENT(INOUT),
TARGET :: mat_w
5856 TYPE(dbcsr_p_type) :: mat_minvvminv
5857 REAL(kind=dp),
DIMENSION(:, :),
INTENT(IN) :: weights_cos_tf_t_to_w, &
5858 weights_sin_tf_t_to_w
5859 COMPLEX(KIND=dp),
DIMENSION(:, :, :, :), &
5860 INTENT(OUT) :: vec_sigma_c_gw
5861 TYPE(qs_environment_type),
POINTER :: qs_env
5862 TYPE(mp_para_env_type),
POINTER :: para_env
5863 TYPE(mp2_type),
INTENT(INOUT) :: mp2_env
5864 INTEGER,
INTENT(IN) :: num_fit_points
5865 TYPE(cp_fm_type),
INTENT(IN) :: fm_mo_coeff
5866 LOGICAL,
INTENT(IN) :: do_ri_sigma_x
5867 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: vec_sigma_x_gw
5868 INTEGER,
INTENT(IN) :: unit_nr, nspins
5869 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc
5870 REAL(kind=dp),
INTENT(IN) :: eps_filter
5872 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_self_energy_cubic_gw_kpoints'
5874 INTEGER :: cut_memory, handle, handle2, i_mem, &
5875 iquad, ispin, j_mem, jquad, &
5876 nkp_self_energy, num_points, &
5878 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist1, dist2, sizes_ao, sizes_ri
5879 INTEGER,
DIMENSION(2) :: mo_end, mo_start, pdims_2d
5880 INTEGER,
DIMENSION(2, 1) :: bounds_ri_i
5881 INTEGER,
DIMENSION(2, 2) :: bounds_ao_ao_j
5882 INTEGER,
DIMENSION(3) :: dims_3c
5883 LOGICAL :: memory_info
5884 REAL(kind=dp) :: omega, t1, t2, tau, weight_cos, &
5886 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: vec_sigma_c_gw_cos_omega, &
5887 vec_sigma_c_gw_cos_tau, vec_sigma_c_gw_neg_tau, vec_sigma_c_gw_pos_tau, &
5888 vec_sigma_c_gw_sin_omega, vec_sigma_c_gw_sin_tau
5889 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_p_greens_fct_occ, &
5890 mat_p_greens_fct_virt
5891 TYPE(dbcsr_type),
TARGET :: mat_greens_fct_occ, mat_greens_fct_virt, mat_mo_coeff, &
5892 mat_self_energy_ao_ao_neg_tau, mat_self_energy_ao_ao_pos_tau
5893 TYPE(dbt_pgrid_type) :: pgrid_2d
5894 TYPE(dbt_type) :: t_3c_m_w_tmp, t_3c_o_all, t_3c_o_w, &
5895 t_ao_tmp, t_greens_fct_occ, &
5896 t_greens_fct_virt, t_ri_tmp, t_w
5898 CALL timeset(routinen, handle)
5900 memory_info = mp2_env%ri_rpa_im_time%memory_info
5901 IF (memory_info)
THEN
5902 unit_nr_prv = unit_nr
5907 cut_memory = mp2_env%ri_rpa_im_time%cut_memory
5909 DO ispin = 1, nspins
5910 mo_start(ispin) = homo(ispin) - gw_corr_lev_occ(ispin) + 1
5911 mo_end(ispin) = homo(ispin) + gw_corr_lev_virt(ispin)
5912 cpassert(mo_end(ispin) - mo_start(ispin) + 1 == gw_corr_lev_tot)
5915 nkp_self_energy = mp2_env%ri_g0w0%nkp_self_energy
5917 vec_sigma_c_gw = z_zero
5918 ALLOCATE (vec_sigma_c_gw_pos_tau(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5919 vec_sigma_c_gw_pos_tau = 0.0_dp
5920 ALLOCATE (vec_sigma_c_gw_neg_tau(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5921 vec_sigma_c_gw_neg_tau = 0.0_dp
5922 ALLOCATE (vec_sigma_c_gw_cos_tau(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5923 vec_sigma_c_gw_cos_tau = 0.0_dp
5924 ALLOCATE (vec_sigma_c_gw_sin_tau(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5925 vec_sigma_c_gw_sin_tau = 0.0_dp
5927 ALLOCATE (vec_sigma_c_gw_cos_omega(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5928 vec_sigma_c_gw_cos_omega = 0.0_dp
5929 ALLOCATE (vec_sigma_c_gw_sin_omega(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5930 vec_sigma_c_gw_sin_omega = 0.0_dp
5932 CALL dbcsr_create(matrix=mat_greens_fct_occ, &
5933 template=matrix_s(1)%matrix, &
5934 matrix_type=dbcsr_type_no_symmetry)
5936 CALL dbcsr_create(matrix=mat_greens_fct_virt, &
5937 template=matrix_s(1)%matrix, &
5938 matrix_type=dbcsr_type_no_symmetry)
5940 CALL dbcsr_create(matrix=mat_self_energy_ao_ao_neg_tau, &
5941 template=matrix_s(1)%matrix, &
5942 matrix_type=dbcsr_type_no_symmetry)
5944 CALL dbcsr_create(matrix=mat_self_energy_ao_ao_pos_tau, &
5945 template=matrix_s(1)%matrix, &
5946 matrix_type=dbcsr_type_no_symmetry)
5948 CALL dbcsr_create(matrix=mat_mo_coeff, &
5949 template=matrix_s(1)%matrix, &
5950 matrix_type=dbcsr_type_no_symmetry)
5952 CALL copy_fm_to_dbcsr(fm_mo_coeff, mat_mo_coeff, keep_sparsity=.false.)
5954 DO ispin = 1, nspins
5955 e_fermi(ispin) = 0.5_dp*(maxval(eigenval(homo, :, ispin)) + minval(eigenval(homo + 1, :, ispin)))
5959 CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d)
5960 ALLOCATE (sizes_ri(dbt_nblks_total(t_3c_o(1, 1), 1)))
5961 CALL dbt_get_info(t_3c_o(1, 1), blk_size_1=sizes_ri)
5963 CALL create_2c_tensor(t_w, dist1, dist2, pgrid_2d, sizes_ri, sizes_ri, name=
"(RI|RI)")
5964 DEALLOCATE (dist1, dist2)
5966 CALL dbt_create(mat_w, t_ri_tmp, name=
"(RI|RI)")
5968 ALLOCATE (sizes_ao(dbt_nblks_total(t_3c_o(1, 1), 2)))
5969 CALL dbt_get_info(t_3c_o(1, 1), blk_size_2=sizes_ao)
5970 CALL create_2c_tensor(t_greens_fct_occ, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5972 DEALLOCATE (dist1, dist2)
5973 CALL create_2c_tensor(t_greens_fct_virt, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5974 DEALLOCATE (dist1, dist2)
5976 CALL dbt_get_info(t_3c_m, nfull_total=dims_3c)
5978 CALL dbt_create(t_3c_o(1, 1), t_3c_o_all, name=
"O (RI AO | AO)")
5981 DO i_mem = 1, cut_memory
5982 CALL decompress_tensor(t_3c_o(1, 1), &
5983 t_3c_o_ind(1, 1, i_mem)%ind, &
5984 t_3c_o_compressed(1, 1, i_mem), &
5985 mp2_env%ri_rpa_im_time%eps_compress)
5986 CALL dbt_copy(t_3c_o(1, 1), t_3c_o_all, summation=.true., move_data=.true.)
5989 CALL dbt_create(t_3c_m, t_3c_m_w_tmp, name=
"M W (RI | AO AO)")
5990 CALL dbt_create(t_3c_o(1, 1), t_3c_o_w, name=
"M W (RI AO | AO)")
5992 CALL dbt_create(mat_greens_fct_occ, t_ao_tmp, name=
"(AO|AO)")
5994 IF (count_ev_sc_gw == 1 .AND. count_sc_gw0 == 1 .AND. do_ri_sigma_x)
THEN
5995 num_points = num_integ_points + 1
5997 num_points = num_integ_points
6000 DO jquad = 1, num_points
6004 IF (jquad <= num_integ_points)
THEN
6007 IF (unit_nr > 0)
WRITE (unit_nr,
'(/T3,A,1X,I3)') &
6008 'GW_INFO| Computing self-energy time point', jquad
6012 IF (unit_nr > 0)
WRITE (unit_nr,
'(/T3,A,1X,I3)') &
6013 'GW_INFO| Computing exchange self-energy'
6016 IF (jquad <= num_integ_points)
THEN
6017 CALL dbcsr_set(mat_w, 0.0_dp)
6018 CALL copy_fm_to_dbcsr(fm_mat_w(jquad), mat_w, keep_sparsity=.false.)
6019 CALL dbt_copy_matrix_to_tensor(mat_w, t_ri_tmp)
6021 CALL dbt_copy_matrix_to_tensor(mat_minvvminv%matrix, t_ri_tmp)
6024 CALL dbt_copy(t_ri_tmp, t_w)
6026 DO ispin = 1, nspins
6028 CALL compute_periodic_dm(mat_p_greens_fct_occ, qs_env, &
6029 ispin, num_points, jquad, e_fermi(ispin), tau, &
6030 remove_occ=.false., remove_virt=.true., &
6031 alloc_dm=(jquad == 1 .AND. ispin == 1))
6033 CALL compute_periodic_dm(mat_p_greens_fct_virt, qs_env, &
6034 ispin, num_points, jquad, e_fermi(ispin), tau, &
6035 remove_occ=.true., remove_virt=.false., &
6036 alloc_dm=(jquad == 1 .AND. ispin == 1))
6038 CALL dbcsr_set(mat_greens_fct_occ, 0.0_dp)
6039 CALL dbcsr_copy(mat_greens_fct_occ, mat_p_greens_fct_occ(jquad, 1)%matrix)
6041 CALL dbcsr_set(mat_greens_fct_virt, 0.0_dp)
6042 CALL dbcsr_copy(mat_greens_fct_virt, mat_p_greens_fct_virt(jquad, 1)%matrix)
6044 CALL dbt_copy_matrix_to_tensor(mat_greens_fct_occ, t_ao_tmp)
6045 CALL dbt_copy(t_ao_tmp, t_greens_fct_occ)
6047 CALL dbt_copy_matrix_to_tensor(mat_greens_fct_virt, t_ao_tmp)
6048 CALL dbt_copy(t_ao_tmp, t_greens_fct_virt)
6050 CALL dbcsr_set(mat_self_energy_ao_ao_neg_tau, 0.0_dp)
6051 CALL dbcsr_set(mat_self_energy_ao_ao_pos_tau, 0.0_dp)
6053 CALL dbt_copy(t_3c_o_all, t_3c_m)
6055 CALL dbt_batched_contract_init(t_3c_o_w)
6059 DO i_mem = 1, cut_memory
6065 bounds_ri_i(:, 1) = [qs_env%mp2_env%ri_rpa_im_time%starts_array_mc_RI(i_mem), &
6066 qs_env%mp2_env%ri_rpa_im_time%ends_array_mc_RI(i_mem)]
6068 DO j_mem = 1, cut_memory
6070 bounds_ao_ao_j(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
6071 bounds_ao_ao_j(:, 2) = [1, dims_3c(3)]
6073 CALL timeset(
"tensor_operation_3c_W", handle2)
6075 CALL dbt_contract(1.0_dp, t_w, t_3c_m, 0.0_dp, &
6077 contract_1=[2], notcontract_1=[1], &
6078 contract_2=[1], notcontract_2=[2, 3], &
6079 map_1=[1], map_2=[2, 3], &
6080 bounds_2=bounds_ri_i, &
6081 bounds_3=bounds_ao_ao_j, &
6082 filter_eps=eps_filter, &
6083 unit_nr=unit_nr_prv)
6085 CALL dbt_copy(t_3c_m_w_tmp, t_3c_o_w, order=[1, 2, 3], move_data=.true.)
6087 CALL timestop(handle2)
6089 CALL contract_to_self_energy(t_3c_o_all, t_greens_fct_occ, t_3c_o_w, &
6090 mat_self_energy_ao_ao_neg_tau, &
6091 bounds_ao_ao_j, bounds_ri_i, unit_nr_prv, &
6092 eps_filter, do_occ=.true., do_virt=.false.)
6094 CALL contract_to_self_energy(t_3c_o_all, t_greens_fct_virt, t_3c_o_w, &
6095 mat_self_energy_ao_ao_pos_tau, &
6096 bounds_ao_ao_j, bounds_ri_i, unit_nr_prv, &
6097 eps_filter, do_occ=.false., do_virt=.true.)
6107 CALL dbt_batched_contract_finalize(t_3c_o_w)
6111 IF (jquad <= num_integ_points)
THEN
6113 CALL trafo_to_mo_and_kpoints(qs_env, mat_self_energy_ao_ao_neg_tau, vec_sigma_c_gw_neg_tau(:, jquad, :, ispin), &
6114 homo(ispin), gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), ispin)
6116 CALL trafo_to_mo_and_kpoints(qs_env, mat_self_energy_ao_ao_pos_tau, vec_sigma_c_gw_pos_tau(:, jquad, :, ispin), &
6117 homo(ispin), gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), ispin)
6119 vec_sigma_c_gw_cos_tau(:, jquad, :, ispin) = 0.5_dp*(vec_sigma_c_gw_pos_tau(:, jquad, :, ispin) + &
6120 vec_sigma_c_gw_neg_tau(:, jquad, :, ispin))
6122 vec_sigma_c_gw_sin_tau(:, jquad, :, ispin) = 0.5_dp*(vec_sigma_c_gw_pos_tau(:, jquad, :, ispin) - &
6123 vec_sigma_c_gw_neg_tau(:, jquad, :, ispin))
6127 vec_sigma_x_gw(mo_start(ispin):mo_end(ispin), :, ispin), &
6128 homo(ispin), gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), ispin)
6136 IF (unit_nr > 0)
WRITE (unit_nr,
'(T6,A,T56,F25.1)')
'Execution time (s):', t2 - t1
6140 IF (count_ev_sc_gw == 1 .AND. count_sc_gw0 == 1)
THEN
6144 IF (do_ri_sigma_x)
THEN
6145 DO ispin = 1, nspins
6146 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, :) = mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, :) + &
6147 vec_sigma_x_gw(:, :, ispin)
6154 DO jquad = 1, num_fit_points
6156 DO iquad = 1, num_integ_points
6160 weight_cos = weights_cos_tf_t_to_w(jquad, iquad)*cos(omega*tau)
6161 weight_sin = weights_sin_tf_t_to_w(jquad, iquad)*sin(omega*tau)
6163 vec_sigma_c_gw_cos_omega(:, jquad, :, :) = vec_sigma_c_gw_cos_omega(:, jquad, :, :) + &
6164 weight_cos*vec_sigma_c_gw_cos_tau(:, iquad, :, :)
6166 vec_sigma_c_gw_sin_omega(:, jquad, :, :) = vec_sigma_c_gw_sin_omega(:, jquad, :, :) + &
6167 weight_sin*vec_sigma_c_gw_sin_tau(:, iquad, :, :)
6175 DO ispin = 1, nspins
6176 vec_sigma_c_gw_sin_omega(1:gw_corr_lev_occ(ispin), :, :, ispin) = &
6177 -vec_sigma_c_gw_sin_omega(1:gw_corr_lev_occ(ispin), :, :, ispin)
6180 vec_sigma_c_gw(:, 1:num_fit_points, :, :) = vec_sigma_c_gw_cos_omega(:, 1:num_fit_points, :, :) + &
6181 gaussi*vec_sigma_c_gw_sin_omega(:, 1:num_fit_points, :, :)
6183 CALL dbt_pgrid_destroy(pgrid_2d)
6185 CALL dbcsr_release(mat_greens_fct_occ)
6186 CALL dbcsr_release(mat_greens_fct_virt)
6187 CALL dbcsr_release(mat_self_energy_ao_ao_neg_tau)
6188 CALL dbcsr_release(mat_self_energy_ao_ao_pos_tau)
6189 CALL dbcsr_release(mat_mo_coeff)
6191 CALL dbcsr_deallocate_matrix_set(mat_p_greens_fct_occ)
6192 CALL dbcsr_deallocate_matrix_set(mat_p_greens_fct_virt)
6194 CALL dbt_destroy(t_w)
6195 CALL dbt_destroy(t_ri_tmp)
6196 CALL dbt_destroy(t_greens_fct_occ)
6197 CALL dbt_destroy(t_greens_fct_virt)
6198 CALL dbt_destroy(t_ao_tmp)
6199 CALL dbt_destroy(t_3c_o_all)
6200 CALL dbt_destroy(t_3c_m_w_tmp)
6201 CALL dbt_destroy(t_3c_o_w)
6203 DEALLOCATE (vec_sigma_c_gw_pos_tau)
6204 DEALLOCATE (vec_sigma_c_gw_neg_tau)
6205 DEALLOCATE (vec_sigma_c_gw_cos_tau)
6206 DEALLOCATE (vec_sigma_c_gw_sin_tau)
6207 DEALLOCATE (vec_sigma_c_gw_cos_omega)
6208 DEALLOCATE (vec_sigma_c_gw_sin_omega)
6210 CALL timestop(handle)
6212 END SUBROUTINE compute_self_energy_cubic_gw_kpoints
6219 TYPE(qs_environment_type),
POINTER :: qs_env
6221 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_minus_vxc_kpoints'
6223 INTEGER :: handle, ikp, ispin, nkp_self_energy, &
6225 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: diag_sigma_x_minus_vxc_mo_mo
6226 TYPE(cp_cfm_type) :: cfm_mo_coeff, ks_mat_ao_ao, &
6227 ks_mat_no_xc_ao_ao, vxc_ao_ao, &
6228 vxc_ao_mo, vxc_mo_mo
6229 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct
6230 TYPE(cp_fm_type) :: fm_dummy, fm_sigma_x_minus_vxc_mo_mo, &
6231 fm_tmp_im, fm_tmp_re
6232 TYPE(dft_control_type),
POINTER :: dft_control
6233 TYPE(kpoint_type),
POINTER :: kpoints_sigma, kpoints_sigma_no_xc
6234 TYPE(mp_para_env_type),
POINTER :: para_env
6236 CALL timeset(routinen, handle)
6238 CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control)
6240 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
6242 kpoints_sigma_no_xc => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc
6244 nkp_self_energy = kpoints_sigma%nkp
6246 nspins = dft_control%nspins
6248 matrix_struct => kpoints_sigma%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct
6250 CALL cp_cfm_create(ks_mat_ao_ao, matrix_struct)
6251 CALL cp_cfm_create(ks_mat_no_xc_ao_ao, matrix_struct)
6252 CALL cp_cfm_create(vxc_ao_ao, matrix_struct)
6253 CALL cp_cfm_create(vxc_ao_mo, matrix_struct)
6254 CALL cp_cfm_create(vxc_mo_mo, matrix_struct)
6255 CALL cp_cfm_create(cfm_mo_coeff, matrix_struct)
6256 CALL cp_fm_create(fm_sigma_x_minus_vxc_mo_mo, matrix_struct)
6257 CALL cp_fm_create(fm_tmp_re, matrix_struct)
6258 CALL cp_fm_create(fm_tmp_im, matrix_struct)
6260 CALL cp_cfm_get_info(cfm_mo_coeff, nrow_global=nmo)
6261 ALLOCATE (diag_sigma_x_minus_vxc_mo_mo(nmo))
6263 DEALLOCATE (qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw)
6265 ALLOCATE (qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(nmo, 2, nkp_self_energy))
6267 DO ikp = 1, nkp_self_energy
6269 DO ispin = 1, nspins
6271 associate(mos => kpoints_sigma%kp_env(ikp)%kpoint_env%mos)
6272 IF (
ASSOCIATED(mos(1, ispin)%mo_coeff))
THEN
6273 CALL cp_fm_copy_general(mos(1, ispin)%mo_coeff, fm_tmp_re, para_env)
6275 CALL cp_fm_copy_general(fm_dummy, fm_tmp_re, para_env)
6277 IF (
ASSOCIATED(mos(2, ispin)%mo_coeff))
THEN
6278 CALL cp_fm_copy_general(mos(2, ispin)%mo_coeff, fm_tmp_im, para_env)
6280 CALL cp_fm_copy_general(fm_dummy, fm_tmp_im, para_env)
6284 CALL cp_fm_to_cfm(fm_tmp_re, fm_tmp_im, cfm_mo_coeff)
6286 CALL cp_fm_to_cfm(kpoints_sigma%kp_env(ikp)%kpoint_env%wmat(1, ispin), &
6287 kpoints_sigma%kp_env(ikp)%kpoint_env%wmat(2, ispin), ks_mat_ao_ao)
6288 associate(wmat => kpoints_sigma_no_xc%kp_env(ikp)%kpoint_env%wmat)
6289 IF (
ASSOCIATED(wmat(1, ispin)%matrix_struct))
THEN
6290 CALL cp_fm_copy_general(wmat(1, ispin), fm_tmp_re, para_env)
6292 CALL cp_fm_copy_general(fm_dummy, fm_tmp_re, para_env)
6294 IF (
ASSOCIATED(wmat(2, ispin)%matrix_struct))
THEN
6295 CALL cp_fm_copy_general(wmat(2, ispin), fm_tmp_im, para_env)
6297 CALL cp_fm_copy_general(fm_dummy, fm_tmp_im, para_env)
6301 CALL cp_fm_to_cfm(fm_tmp_re, fm_tmp_im, vxc_ao_ao)
6303 CALL parallel_gemm(
'N',
'N', nmo, nmo, nmo, z_one, vxc_ao_ao, cfm_mo_coeff, z_zero, vxc_ao_mo)
6304 CALL parallel_gemm(
'C',
'N', nmo, nmo, nmo, z_one, cfm_mo_coeff, vxc_ao_mo, z_zero, vxc_mo_mo)
6306 CALL cp_cfm_to_fm(vxc_mo_mo, fm_sigma_x_minus_vxc_mo_mo)
6308 CALL cp_fm_get_diag(fm_sigma_x_minus_vxc_mo_mo, diag_sigma_x_minus_vxc_mo_mo)
6310 qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, ikp) = diag_sigma_x_minus_vxc_mo_mo(:)
6316 CALL cp_cfm_release(ks_mat_ao_ao)
6317 CALL cp_cfm_release(ks_mat_no_xc_ao_ao)
6318 CALL cp_cfm_release(vxc_ao_ao)
6319 CALL cp_cfm_release(vxc_ao_mo)
6320 CALL cp_cfm_release(vxc_mo_mo)
6321 CALL cp_cfm_release(cfm_mo_coeff)
6322 CALL cp_fm_release(fm_sigma_x_minus_vxc_mo_mo)
6323 CALL cp_fm_release(fm_tmp_re)
6324 CALL cp_fm_release(fm_tmp_im)
6326 DEALLOCATE (diag_sigma_x_minus_vxc_mo_mo)
6328 CALL timestop(handle)
6343 homo, gw_corr_lev_occ, gw_corr_lev_virt, ispin)
6344 TYPE(qs_environment_type),
POINTER :: qs_env
6345 TYPE(dbcsr_type),
TARGET :: mat_self_energy_ao_ao
6346 REAL(kind=dp),
DIMENSION(:, :) :: vec_sigma
6347 INTEGER :: homo, gw_corr_lev_occ, gw_corr_lev_virt, &
6350 CHARACTER(LEN=*),
PARAMETER :: routinen =
'trafo_to_mo_and_kpoints'
6352 INTEGER :: handle, ikp, nkp_self_energy, nmo, &
6353 periodic(3), size_real_space
6354 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: diag_self_energy
6355 TYPE(cell_type),
POINTER :: cell
6356 TYPE(cp_cfm_type) :: cfm_mo_coeff, cfm_self_energy_ao_ao, &
6357 cfm_self_energy_ao_mo, &
6358 cfm_self_energy_mo_mo
6359 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct
6360 TYPE(cp_fm_type) :: fm_self_energy_mo_mo
6361 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mat_self_energy_ao_ao_kp_im, &
6362 mat_self_energy_ao_ao_kp_re, mat_self_energy_ao_ao_real_space
6363 TYPE(kpoint_type),
POINTER :: kpoints_sigma
6364 TYPE(mp_para_env_type),
POINTER :: para_env
6366 CALL timeset(routinen, handle)
6368 CALL get_qs_env(qs_env, cell=cell, para_env=para_env)
6369 CALL get_cell(cell=cell, periodic=periodic)
6371 size_real_space = 3**(periodic(1) + periodic(2) + periodic(3))
6373 CALL alloc_mat_set(mat_self_energy_ao_ao_real_space, size_real_space, mat_self_energy_ao_ao)
6375 CALL dbcsr_copy(mat_self_energy_ao_ao_real_space(1)%matrix, mat_self_energy_ao_ao)
6377 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
6379 CALL get_mat_cell_t_from_mat_gamma(mat_self_energy_ao_ao_real_space, qs_env, kpoints_sigma, 0, 0)
6381 nkp_self_energy = kpoints_sigma%nkp
6383 CALL alloc_mat_set(mat_self_energy_ao_ao_kp_re, nkp_self_energy, mat_self_energy_ao_ao)
6384 CALL alloc_mat_set(mat_self_energy_ao_ao_kp_im, nkp_self_energy, mat_self_energy_ao_ao)
6386 CALL real_space_to_kpoint_transform_rpa(mat_self_energy_ao_ao_kp_re, mat_self_energy_ao_ao_kp_im, &
6387 mat_self_energy_ao_ao_real_space, kpoints_sigma, 1.0e-50_dp)
6389 CALL dbcsr_get_info(mat_self_energy_ao_ao, nfullrows_total=nmo)
6390 ALLOCATE (diag_self_energy(nmo))
6392 matrix_struct => kpoints_sigma%kp_env(1)%kpoint_env%mos(1, 1)%mo_coeff%matrix_struct
6394 CALL cp_cfm_create(cfm_self_energy_ao_ao, matrix_struct)
6395 CALL cp_cfm_create(cfm_self_energy_ao_mo, matrix_struct)
6396 CALL cp_cfm_create(cfm_self_energy_mo_mo, matrix_struct)
6397 CALL cp_cfm_set_all(cfm_self_energy_ao_ao, z_zero)
6398 CALL cp_cfm_set_all(cfm_self_energy_ao_mo, z_zero)
6399 CALL cp_cfm_set_all(cfm_self_energy_mo_mo, z_zero)
6401 CALL cp_fm_create(fm_self_energy_mo_mo, matrix_struct)
6402 CALL cp_cfm_create(cfm_mo_coeff, matrix_struct)
6404 DO ikp = 1, nkp_self_energy
6406 CALL dbcsr_to_cfm(mat_self_energy_ao_ao_kp_re(ikp)%matrix, &
6407 mat_self_energy_ao_ao_kp_im(ikp)%matrix, cfm_self_energy_ao_ao)
6409 CALL cp_fm_to_cfm(kpoints_sigma%kp_env(ikp)%kpoint_env%mos(1, ispin)%mo_coeff, &
6410 kpoints_sigma%kp_env(ikp)%kpoint_env%mos(2, ispin)%mo_coeff, cfm_mo_coeff)
6412 CALL parallel_gemm(
'N',
'N', nmo, nmo, nmo, z_one, cfm_self_energy_ao_ao, cfm_mo_coeff, &
6413 z_zero, cfm_self_energy_ao_mo)
6415 CALL parallel_gemm(
'C',
'N', nmo, nmo, nmo, z_one, cfm_mo_coeff, cfm_self_energy_ao_mo, &
6416 z_zero, cfm_self_energy_mo_mo)
6418 CALL cp_cfm_to_fm(cfm_self_energy_mo_mo, fm_self_energy_mo_mo)
6420 CALL cp_fm_get_diag(fm_self_energy_mo_mo, diag_self_energy)
6422 vec_sigma(:, ikp) = diag_self_energy(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt)
6426 CALL dbcsr_deallocate_matrix_set(mat_self_energy_ao_ao_real_space)
6427 CALL dbcsr_deallocate_matrix_set(mat_self_energy_ao_ao_kp_re)
6428 CALL dbcsr_deallocate_matrix_set(mat_self_energy_ao_ao_kp_im)
6430 CALL cp_cfm_release(cfm_self_energy_ao_ao)
6431 CALL cp_cfm_release(cfm_self_energy_ao_mo)
6432 CALL cp_cfm_release(cfm_self_energy_mo_mo)
6433 CALL cp_cfm_release(cfm_mo_coeff)
6434 CALL cp_fm_release(fm_self_energy_mo_mo)
6436 DEALLOCATE (diag_self_energy)
6438 CALL timestop(handle)
6448 SUBROUTINE dbcsr_to_cfm(dbcsr_re, dbcsr_im, cfm_mat)
6450 TYPE(dbcsr_type),
POINTER :: dbcsr_re, dbcsr_im
6451 TYPE(cp_cfm_type),
INTENT(IN) :: cfm_mat
6453 CHARACTER(LEN=*),
PARAMETER :: routinen =
'dbcsr_to_cfm'
6456 TYPE(cp_fm_type) :: fm_mat_im, fm_mat_re
6458 CALL timeset(routinen, handle)
6460 CALL cp_fm_create(fm_mat_re, cfm_mat%matrix_struct)
6461 CALL cp_fm_create(fm_mat_im, cfm_mat%matrix_struct)
6462 CALL cp_fm_set_all(fm_mat_re, 0.0_dp)
6463 CALL cp_fm_set_all(fm_mat_im, 0.0_dp)
6465 CALL copy_dbcsr_to_fm(dbcsr_re, fm_mat_re)
6466 CALL copy_dbcsr_to_fm(dbcsr_im, fm_mat_im)
6468 CALL cp_fm_to_cfm(fm_mat_re, fm_mat_im, cfm_mat)
6470 CALL cp_fm_release(fm_mat_re)
6471 CALL cp_fm_release(fm_mat_im)
6473 CALL timestop(handle)
6475 END SUBROUTINE dbcsr_to_cfm
6484 SUBROUTINE alloc_mat_set(mat_set, mat_size, template, explicitly_no_symmetry)
6485 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mat_set
6486 INTEGER,
INTENT(IN) :: mat_size
6487 TYPE(dbcsr_type),
TARGET :: template
6488 LOGICAL,
OPTIONAL :: explicitly_no_symmetry
6490 CHARACTER(LEN=*),
PARAMETER :: routinen =
'alloc_mat_set'
6492 INTEGER :: handle, i_size
6493 LOGICAL :: my_explicitly_no_symmetry
6495 CALL timeset(routinen, handle)
6497 my_explicitly_no_symmetry = .false.
6498 IF (
PRESENT(explicitly_no_symmetry)) my_explicitly_no_symmetry = explicitly_no_symmetry
6501 CALL dbcsr_allocate_matrix_set(mat_set, mat_size)
6502 DO i_size = 1, mat_size
6503 ALLOCATE (mat_set(i_size)%matrix)
6504 IF (my_explicitly_no_symmetry)
THEN
6505 CALL dbcsr_create(matrix=mat_set(i_size)%matrix, template=template, &
6506 matrix_type=dbcsr_type_no_symmetry)
6508 CALL dbcsr_create(matrix=mat_set(i_size)%matrix, template=template)
6510 CALL dbcsr_copy(mat_set(i_size)%matrix, template)
6511 CALL dbcsr_set(mat_set(i_size)%matrix, 0.0_dp)
6514 CALL timestop(handle)
6516 END SUBROUTINE alloc_mat_set
6526 SUBROUTINE alloc_mat_set_2d(mat_set, mat_size_1, mat_size_2, template, explicitly_no_symmetry)
6527 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_set
6528 INTEGER,
INTENT(IN) :: mat_size_1, mat_size_2
6529 TYPE(dbcsr_type),
TARGET :: template
6530 LOGICAL,
OPTIONAL :: explicitly_no_symmetry
6532 CHARACTER(LEN=*),
PARAMETER :: routinen =
'alloc_mat_set_2d'
6534 INTEGER :: handle, i_size, j_size
6535 LOGICAL :: my_explicitly_no_symmetry
6537 CALL timeset(routinen, handle)
6539 my_explicitly_no_symmetry = .false.
6540 IF (
PRESENT(explicitly_no_symmetry)) my_explicitly_no_symmetry = explicitly_no_symmetry
6543 CALL dbcsr_allocate_matrix_set(mat_set, mat_size_1, mat_size_2)
6544 DO i_size = 1, mat_size_1
6545 DO j_size = 1, mat_size_2
6546 ALLOCATE (mat_set(i_size, j_size)%matrix)
6547 IF (my_explicitly_no_symmetry)
THEN
6548 CALL dbcsr_create(matrix=mat_set(i_size, j_size)%matrix, template=template, &
6549 matrix_type=dbcsr_type_no_symmetry)
6551 CALL dbcsr_create(matrix=mat_set(i_size, j_size)%matrix, template=template)
6553 CALL dbcsr_copy(mat_set(i_size, j_size)%matrix, template)
6554 CALL dbcsr_set(mat_set(i_size, j_size)%matrix, 0.0_dp)
6558 CALL timestop(handle)
6560 END SUBROUTINE alloc_mat_set_2d
6575 SUBROUTINE contract_to_self_energy(t_3c_O_all, t_greens_fct, t_3c_O_W, &
6576 mat_self_energy_ao_ao, bounds_ao_ao_j, bounds_RI_i, &
6577 unit_nr, eps_filter, do_occ, do_virt)
6579 TYPE(dbt_type) :: t_3c_o_all, t_greens_fct, t_3c_o_w
6580 TYPE(dbcsr_type),
TARGET :: mat_self_energy_ao_ao
6581 INTEGER,
DIMENSION(2, 2) :: bounds_ao_ao_j
6582 INTEGER,
DIMENSION(2, 1) :: bounds_ri_i
6584 REAL(kind=dp) :: eps_filter
6585 LOGICAL :: do_occ, do_virt
6587 CHARACTER(LEN=*),
PARAMETER :: routinen =
'contract_to_self_energy'
6590 INTEGER,
DIMENSION(2, 1) :: bounds_ao_j
6591 INTEGER,
DIMENSION(2, 2) :: bounds_ao_all_ri_i, bounds_ri_i_ao_j
6592 REAL(kind=dp) :: sign_self_energy
6593 TYPE(dbt_type) :: t_3c_o_g, t_3c_o_g_tmp, t_self_energy, &
6596 CALL timeset(routinen, handle)
6598 cpassert(do_occ .EQV. (.NOT. do_virt))
6600 CALL dbt_create(t_3c_o_all, t_3c_o_g, name=
"M occ (RI AO | AO)")
6601 CALL dbt_create(t_3c_o_all, t_3c_o_g_tmp, name=
"M occ (RI AO | AO)")
6602 CALL dbt_create(t_greens_fct, t_self_energy, name=
"(AO|AO)")
6603 CALL dbt_create(mat_self_energy_ao_ao, t_self_energy_tmp)
6605 bounds_ao_j(:, 1) = bounds_ao_ao_j(:, 1)
6606 bounds_ao_all_ri_i(:, 1) = bounds_ri_i(:, 1)
6607 bounds_ao_all_ri_i(:, 2) = bounds_ao_ao_j(:, 2)
6609 CALL dbt_contract(1.0_dp, t_greens_fct, t_3c_o_all, 0.0_dp, &
6611 contract_1=[2], notcontract_1=[1], &
6612 contract_2=[3], notcontract_2=[1, 2], &
6613 map_1=[3], map_2=[1, 2], &
6614 bounds_2=bounds_ao_j, &
6615 bounds_3=bounds_ao_all_ri_i, &
6616 filter_eps=eps_filter, &
6619 CALL dbt_copy(t_3c_o_g_tmp, t_3c_o_g, order=[1, 3, 2], move_data=.true.)
6621 IF (do_occ) sign_self_energy = -1.0_dp
6622 IF (do_virt) sign_self_energy = 1.0_dp
6624 bounds_ri_i_ao_j(:, 1) = bounds_ri_i(:, 1)
6625 bounds_ri_i_ao_j(:, 2) = bounds_ao_ao_j(:, 1)
6627 CALL dbt_contract(sign_self_energy, t_3c_o_w, t_3c_o_g, 0.0_dp, &
6629 contract_1=[1, 2], notcontract_1=[3], &
6630 contract_2=[1, 2], notcontract_2=[3], &
6631 map_1=[1], map_2=[2], &
6632 bounds_1=bounds_ri_i_ao_j, &
6633 filter_eps=eps_filter, &
6636 CALL dbt_copy(t_self_energy, t_self_energy_tmp)
6637 CALL dbt_clear(t_self_energy)
6639 CALL dbt_copy_tensor_to_matrix(t_self_energy_tmp, mat_self_energy_ao_ao, summation=.true.)
6641 CALL dbt_destroy(t_3c_o_g)
6642 CALL dbt_destroy(t_3c_o_g_tmp)
6643 CALL dbt_destroy(t_self_energy)
6644 CALL dbt_destroy(t_self_energy_tmp)
6646 CALL timestop(handle)
6648 END SUBROUTINE contract_to_self_energy
6663 SUBROUTINE contract_cubic_gw(t_3c_overl_int_gw_AO, t_3c_overl_int_gw_RI, &
6664 t_AO, t_RI, prefac, &
6665 mo_bounds, unit_nr, &
6666 t_3c_ctr_RI, t_3c_ctr_AO, calculate_ctr_RI)
6667 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_overl_int_gw_ao, &
6668 t_3c_overl_int_gw_ri, t_ao, t_ri
6669 REAL(dp),
DIMENSION(2),
INTENT(IN) :: prefac
6670 INTEGER,
DIMENSION(2),
INTENT(IN) :: mo_bounds
6671 INTEGER,
INTENT(IN) :: unit_nr
6672 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_ctr_ri, t_3c_ctr_ao
6673 LOGICAL,
INTENT(IN) :: calculate_ctr_ri
6675 CHARACTER(LEN=*),
PARAMETER :: routinen =
'contract_cubic_gw'
6678 INTEGER,
DIMENSION(2, 2) :: ctr_bounds_mo
6679 INTEGER,
DIMENSION(3) :: bounds_3c
6681 CALL timeset(routinen, handle)
6683 IF (calculate_ctr_ri)
THEN
6684 CALL dbt_get_info(t_3c_overl_int_gw_ri, nfull_total=bounds_3c)
6685 ctr_bounds_mo(:, 1) = [1, bounds_3c(2)]
6686 ctr_bounds_mo(:, 2) = mo_bounds
6688 CALL dbt_contract(prefac(1), t_ri, t_3c_overl_int_gw_ri, 0.0_dp, &
6690 contract_1=[2], notcontract_1=[1], &
6691 contract_2=[1], notcontract_2=[2, 3], &
6692 map_1=[1], map_2=[2, 3], &
6693 bounds_3=ctr_bounds_mo, &
6698 CALL dbt_get_info(t_3c_overl_int_gw_ao, nfull_total=bounds_3c)
6699 ctr_bounds_mo(:, 1) = [1, bounds_3c(2)]
6700 ctr_bounds_mo(:, 2) = mo_bounds
6702 CALL dbt_contract(prefac(2), t_ao, t_3c_overl_int_gw_ao, 0.0_dp, &
6704 contract_1=[2], notcontract_1=[1], &
6705 contract_2=[1], notcontract_2=[2, 3], &
6706 map_1=[1], map_2=[2, 3], &
6707 bounds_3=ctr_bounds_mo, &
6710 CALL timestop(handle)
6712 END SUBROUTINE contract_cubic_gw
6723 SUBROUTINE trace_sigma_gw(t3c_1, t3c_2, vec_sigma, mo_offset, mo_bounds, para_env)
6724 TYPE(dbt_type),
INTENT(INOUT) :: t3c_1, t3c_2
6725 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: vec_sigma
6726 INTEGER,
INTENT(IN) :: mo_offset
6727 INTEGER,
DIMENSION(2),
INTENT(IN) :: mo_bounds
6728 TYPE(mp_para_env_type),
INTENT(IN) :: para_env
6730 CHARACTER(LEN=*),
PARAMETER :: routinen =
'trace_sigma_gw'
6732 INTEGER :: handle, n, n_end, n_end_block, n_start, &
6734 INTEGER,
DIMENSION(1) :: trace_shape
6735 INTEGER,
DIMENSION(2) :: mo_bounds_off
6736 INTEGER,
DIMENSION(3) :: boff, bsize, ind
6738 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: block_1, block_2
6740 DIMENSION(mo_bounds(2)-mo_bounds(1)+1) :: vec_sigma_prv
6741 TYPE(dbt_iterator_type) :: iter
6742 TYPE(dbt_type) :: t3c_1_redist
6744 CALL timeset(routinen, handle)
6746 CALL dbt_create(t3c_2, t3c_1_redist)
6747 CALL dbt_copy(t3c_1, t3c_1_redist, order=[2, 1, 3], move_data=.true.)
6749 vec_sigma_prv = 0.0_dp
6755 CALL dbt_iterator_start(iter, t3c_1_redist)
6756 DO WHILE (dbt_iterator_blocks_left(iter))
6757 CALL dbt_iterator_next_block(iter, ind, blk_size=bsize, blk_offset=boff)
6758 CALL dbt_get_block(t3c_1_redist, ind, block_1, found)
6760 CALL dbt_get_block(t3c_2, ind, block_2, found)
6761 IF (.NOT. found) cycle
6763 IF (boff(3) < mo_bounds(1))
THEN
6764 n_start_block = mo_bounds(1) - boff(3) + 1
6768 n_start = boff(3) - mo_bounds(1) + 1
6771 IF (boff(3) + bsize(3) - 1 > mo_bounds(2))
THEN
6772 n_end_block = mo_bounds(2) - boff(3) + 1
6773 n_end = mo_bounds(2) - mo_bounds(1) + 1
6775 n_end_block = bsize(3)
6776 n_end = boff(3) + bsize(3) - mo_bounds(1)
6779 trace_shape(1) =
SIZE(block_1, 1)*
SIZE(block_1, 2)
6780 vec_sigma_prv(n_start:n_end) = &
6781 vec_sigma_prv(n_start:n_end) + &
6782 [(dot_product(reshape(block_1(:, :, n), trace_shape), &
6783 reshape(block_2(:, :, n), trace_shape)), &
6784 n=n_start_block, n_end_block)]
6785 DEALLOCATE (block_1, block_2)
6787 CALL dbt_iterator_stop(iter)
6790 CALL dbt_destroy(t3c_1_redist)
6792 CALL para_env%sum(vec_sigma_prv)
6794 mo_bounds_off = mo_bounds - mo_offset + 1
6795 vec_sigma(mo_bounds_off(1):mo_bounds_off(2)) = &
6796 vec_sigma(mo_bounds_off(1):mo_bounds_off(2)) + vec_sigma_prv
6798 CALL timestop(handle)
6799 END SUBROUTINE trace_sigma_gw
6818 SUBROUTINE compute_greens_function_time(mat_greens_fct_occ, mat_greens_fct_virt, fm_mo_coeff_occ, fm_mo_coeff_virt, &
6819 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
6820 fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, Eigenval, nmo, &
6821 eps_filter, e_fermi, tau, para_env)
6823 TYPE(dbcsr_type),
INTENT(INOUT) :: mat_greens_fct_occ, mat_greens_fct_virt
6824 TYPE(cp_fm_type),
INTENT(IN) :: fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
6825 fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau
6826 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: eigenval
6827 INTEGER,
INTENT(IN) :: nmo
6828 REAL(kind=dp),
INTENT(IN) :: eps_filter, e_fermi, tau
6829 TYPE(mp_para_env_type),
INTENT(IN) :: para_env
6831 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_Greens_function_time'
6833 INTEGER :: handle, i_global, iib, jjb, ncol_local, &
6835 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
6836 REAL(kind=dp) :: stabilize_exp
6838 CALL timeset(routinen, handle)
6840 CALL para_env%sync()
6843 CALL cp_fm_get_info(matrix=fm_mo_coeff_occ, &
6844 nrow_local=nrow_local, &
6845 ncol_local=ncol_local, &
6846 row_indices=row_indices, &
6847 col_indices=col_indices)
6853 stabilize_exp = 70.0_dp
6856 DO jjb = 1, nrow_local
6857 DO iib = 1, ncol_local
6858 i_global = col_indices(iib)
6860 IF (abs(tau*0.5_dp*(eigenval(i_global) - e_fermi)) < stabilize_exp)
THEN
6861 fm_mo_coeff_occ_scaled%local_data(jjb, iib) = &
6862 fm_mo_coeff_occ%local_data(jjb, iib)*exp(tau*0.5_dp*(eigenval(i_global) - e_fermi))
6864 fm_mo_coeff_occ_scaled%local_data(jjb, iib) = 0.0_dp
6871 DO jjb = 1, nrow_local
6872 DO iib = 1, ncol_local
6873 i_global = col_indices(iib)
6875 IF (abs(tau*0.5_dp*(eigenval(i_global) - e_fermi)) < stabilize_exp)
THEN
6876 fm_mo_coeff_virt_scaled%local_data(jjb, iib) = &
6877 fm_mo_coeff_virt%local_data(jjb, iib)*exp(-tau*0.5_dp*(eigenval(i_global) - e_fermi))
6879 fm_mo_coeff_virt_scaled%local_data(jjb, iib) = 0.0_dp
6885 CALL para_env%sync()
6887 CALL parallel_gemm(transa=
"N", transb=
"T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, &
6888 matrix_a=fm_mo_coeff_occ_scaled, matrix_b=fm_mo_coeff_occ_scaled, beta=0.0_dp, &
6889 matrix_c=fm_scaled_dm_occ_tau)
6891 CALL parallel_gemm(transa=
"N", transb=
"T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, &
6892 matrix_a=fm_mo_coeff_virt_scaled, matrix_b=fm_mo_coeff_virt_scaled, beta=0.0_dp, &
6893 matrix_c=fm_scaled_dm_virt_tau)
6895 CALL dbcsr_set(mat_greens_fct_occ, 0.0_dp)
6897 CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, &
6898 mat_greens_fct_occ, &
6899 keep_sparsity=.false.)
6901 CALL dbcsr_filter(mat_greens_fct_occ, eps_filter)
6903 CALL dbcsr_set(mat_greens_fct_virt, 0.0_dp)
6905 CALL copy_fm_to_dbcsr(fm_scaled_dm_virt_tau, &
6906 mat_greens_fct_virt, &
6907 keep_sparsity=.false.)
6909 CALL dbcsr_filter(mat_greens_fct_virt, eps_filter)
6911 CALL timestop(handle)
6913 END SUBROUTINE compute_greens_function_time
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds, sdab, pab, force_a)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Define the atomic kind types and their sub types.
Handles all functions related to the CELL.
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltar, matrix_l, atcore)
...
Basic linear algebra operations for complex full matrices.
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
subroutine, public cp_cfm_transpose(matrix, trans, matrixt)
Transposes a BLACS distributed complex matrix.
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
used for collecting diagonalization schemes available for cp_cfm_type
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
subroutine, public cp_cfm_create(matrix, matrix_struct, name, set_zero)
Creates a new full matrix with the given structure.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_release_p(matrix)
...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Utility routines to open and close files. Tracking of preconnections.
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_invert(matrix, n, info_out)
used to replace the cholesky decomposition by the inverse
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, zero_tails, silent, mpi_io)
...
This is the start of a dbt_api, all publically needed functions are exported here....
Types and set/get functions for HFX.
subroutine, public dealloc_containers(data, memory_usage)
...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_path_length
Routines needed for kpoint calculation.
subroutine, public kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
generate real space density matrices in DBCSR format
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
subroutine, public kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
Calculate kpoint density matrices (rho(k), owned by kpoint groups)
Types and basic routines needed for a kpoint calculation.
subroutine, public kpoint_sym_create(kp_sym)
Create a single kpoint symmetry environment.
subroutine, public kpoint_release(kpoint)
Release a kpoint environment, deallocate all data.
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Machine interface based on Fortran 2003 and POSIX.
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
complex(kind=dp), parameter, public z_one
complex(kind=dp), parameter, public gaussi
real(kind=dp), parameter, public fourpi
real(kind=dp), parameter, public twopi
complex(kind=dp), parameter, public z_zero
Interface to the message passing library MPI.
Types needed for MP2 calculations.
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public evolt
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculation of band structures.
subroutine, public calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgeneral, group_size_ext)
diagonalize KS matrices at a set of kpoints
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public qs_env_release(qs_env)
releases the given qs_env (see doc/ReferenceCounting.html)
Initialize a qs_env for kpoint calculations starting from a gamma point qs_env.
subroutine, public create_kp_from_gamma(qs_env, qs_env_kp, with_xc_terms)
...
Some utility functions for the calculation of integrals.
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition and initialisation of the mo data type.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
subroutine, public build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
...
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Generate the atomic neighbor lists.
subroutine, public setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, mic, symmetric, molecular, operator_type)
Build a neighborlist.
Calculation of overlap matrix, its derivatives and forces.
subroutine, public build_overlap_matrix_simple(ks_env, matrix_s, basis_set_list_a, basis_set_list_b, sab_nl)
Calculation of the overlap matrix over Cartesian Gaussian functions.
module that contains the definitions of the scf types
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Utility methods to build 3-center integral tensors of various types.
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
Utility methods to build 3-center integral tensors of various types.
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
Routines to calculate image charge corrections.
subroutine, public apply_ic_corr(eigenval, eigenval_scf, ic_corr_list, gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, homo, nmo, unit_nr, do_alpha, do_beta)
...
Utility routines for GW with imaginary time.
subroutine, public get_tensor_3c_overl_int_gw(t_3c_overl_int, t_3c_o_compressed, t_3c_o_ind, t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, starts_array_mc, ends_array_mc, mo_coeff, matrix_s, gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, para_env, do_ic_model, t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, qs_env, unit_nr, do_alpha)
...
Routines treating GW and RPA calculations with kpoints.
subroutine, public get_mat_cell_t_from_mat_gamma(mat_p_omega, qs_env, kpoints, jquad, unit_nr)
...
subroutine, public real_space_to_kpoint_transform_rpa(real_mat_kp, imag_mat_kp, mat_real_space, kpoints, eps_filter_im_time, real_mat_real_space)
...
subroutine, public mat_kp_from_mat_gamma(qs_env, mat_kp, mat_gamma, kpoints, ispin, real_mat_real_space)
...
Routines for GW, continuous development [Jan Wilhelm].
subroutine, public deallocate_matrices_gw(fm_mat_s_gw_work, vec_w_gw, vec_sigma_c_gw, vec_omega_fit_gw, vec_sigma_x_minus_vxc_gw, eigenval_last, eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, kpoints, vec_sigma_x_gw, my_do_gw)
...
subroutine, public allocate_matrices_gw(vec_sigma_c_gw, color_rpa_group, dimen_nm_gw, gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, num_integ_group, num_integ_points, unit_nr, gw_corr_lev_tot, num_fit_points, omega_max_fit, do_minimax_quad, do_periodic, do_ri_sigma_x, my_do_gw, first_cycle_periodic_correction, a_scaling, eigenval, tj, vec_omega_fit_gw, vec_sigma_x_gw, delta_corr, eigenval_last, eigenval_scf, vec_w_gw, fm_mat_s_gw, fm_mat_s_gw_work, para_env, mp2_env, kpoints, nkp, nkp_self_energy, do_kpoints_cubic_rpa, do_kpoints_from_gamma)
...
subroutine, public compute_gw_self_energy(vec_sigma_c_gw, dimen_nm_gw, dimen_ri, gw_corr_lev_occ, gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, do_im_time, do_periodic, first_cycle_periodic_correction, fermi_level_offset, omega, eigenval, delta_corr, vec_omega_fit_gw, vec_w_gw, wj, fm_mat_q, fm_mat_r_gw, fm_mat_s_gw, fm_mat_s_gw_work, mo_coeff, para_env, para_env_rpa, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, kpoints, qs_env, mp2_env)
...
subroutine, public get_fermi_level_offset(fermi_level_offset, fermi_level_offset_input, eigenval, homo)
...
subroutine, public trafo_to_mo_and_kpoints(qs_env, mat_self_energy_ao_ao, vec_sigma, homo, gw_corr_lev_occ, gw_corr_lev_virt, ispin)
...
subroutine, public allocate_matrices_gw_im_time(gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, num_integ_points, unit_nr, ri_blk_sizes, do_ic_model, para_env, fm_mat_w, fm_mat_q, mo_coeff, t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, starts_array_mc, ends_array_mc, t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, matrix_s, mat_w, t_3c_overl_int, t_3c_o_compressed, t_3c_o_ind, qs_env)
...
subroutine, public deallocate_matrices_gw_im_time(weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, do_ic_model, do_kpoints_cubic_rpa, fm_mat_w, t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, mat_w, qs_env)
...
subroutine, public compute_minus_vxc_kpoints(qs_env)
...
subroutine, public compute_qp_energies(vec_sigma_c_gw, count_ev_sc_gw, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, homo, nmo, num_fit_points, num_integ_points, unit_nr, do_apply_ic_corr_to_gw, do_im_time, do_periodic, do_ri_sigma_x, first_cycle_periodic_correction, e_fermi, eps_filter, fermi_level_offset, delta_corr, eigenval, eigenval_last, eigenval_scf, iter_sc_gw0, exit_ev_gw, tau_tj, tj, vec_omega_fit_gw, vec_sigma_x_gw, ic_corr_list, weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, mo_coeff, fm_mat_w, para_env, para_env_rpa, mat_dm, mat_minvvminv, t_3c_o, t_3c_m, t_3c_overl_int_ao_mo, t_3c_o_compressed, t_3c_o_mo_compressed, t_3c_o_ind, t_3c_o_mo_ind, t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, mat_w, matrix_s, kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_rpa, starts_array_mc, ends_array_mc)
...
subroutine, public compute_w_cubic_gw(fm_mat_w, fm_mat_q, fm_mat_work, dimen_ri, fm_mat_l, num_integ_points, tj, tau_tj, weights_cos_tf_w_to_t, jquad, omega)
...
subroutine, public continuation_pade(vec_gw_energ, vec_omega_fit_gw, z_value, m_value, vec_sigma_c_gw, vec_sigma_x_minus_vxc_gw, eigenval, eigenval_scf, do_hedin_shift, n_level_gw, gw_corr_lev_occ, gw_corr_lev_vir, nparam_pade, num_fit_points, crossing_search, homo, fermi_level_offset, do_gw_im_time, print_self_energy, count_ev_sc_gw, vec_gw_dos, dos_lower_bound, dos_precision, ndos, min_level_self_energy, max_level_self_energy, dos_eta, dos_min, dos_max, e_fermi_ext)
perform analytic continuation with pade approximation
Routines for low-scaling RPA/GW with imaginary time.
subroutine, public compute_periodic_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, remove_occ, remove_virt, alloc_dm)
...
parameters that control an scf iteration
All kind of helpful little routines.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Represent a complex full matrix.
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information about kpoints.
stores all the informations relevant to an mpi environment
represent a list of objects
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...