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)
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 .LE. nao)
THEN
1832 e_i = eigenval(i_glob, ikp, 1)
1834 e_i = eigenval(i_glob - nao, ikp, 1)
1836 IF (j_glob .LE. 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)
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)
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) .NE. array(i, j, k)) array(i, j, k) = real_value
2193 CALL timestop(handle)
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 .NE. 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 .NE. 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)
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)
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_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
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.
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_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
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
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, 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 qs_env_release(qs_env)
releases the given qs_env (see doc/ReferenceCounting.html)
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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
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, 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, 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 ...