39 dbcsr_type_no_symmetry
66 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
67 dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, &
68 dbt_get_block, dbt_get_info, dbt_iterator_blocks_left, dbt_iterator_next_block, &
69 dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, dbt_nblks_total, &
70 dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
151#include "./base/base_uses.f90"
157 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'rpa_gw'
196 num_integ_points, unit_nr, &
197 RI_blk_sizes, do_ic_model, &
198 para_env, fm_mat_W, fm_mat_Q, &
200 t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
201 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
202 starts_array_mc, ends_array_mc, &
203 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
204 matrix_s, mat_W, t_3c_overl_int, &
205 t_3c_O_compressed, t_3c_O_ind, &
208 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
209 INTEGER,
INTENT(IN) :: nmo, num_integ_points, unit_nr
210 INTEGER,
DIMENSION(:),
POINTER :: ri_blk_sizes
211 LOGICAL,
INTENT(IN) :: do_ic_model
213 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
214 INTENT(OUT) :: fm_mat_w
216 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: mo_coeff
217 TYPE(dbt_type) :: t_3c_overl_int_ao_mo
219 DIMENSION(:) :: t_3c_o_mo_compressed
221 DIMENSION(:),
INTENT(OUT) :: t_3c_o_mo_ind
222 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:), &
223 INTENT(INOUT) :: t_3c_overl_int_gw_ri, &
225 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc
226 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:), &
227 INTENT(INOUT) :: t_3c_overl_nnp_ic, &
228 t_3c_overl_nnp_ic_reflected
231 TYPE(dbt_type),
DIMENSION(:, :) :: t_3c_overl_int
236 CHARACTER(LEN=*),
PARAMETER :: routinen =
'allocate_matrices_gw_im_time'
238 INTEGER :: handle, jquad, nspins
239 LOGICAL :: my_open_shell
240 TYPE(dbt_type) :: t_3c_overl_int_ao_mo_beta
242 CALL timeset(routinen, handle)
245 my_open_shell = (nspins == 2)
247 ALLOCATE (t_3c_o_mo_ind(nspins), t_3c_overl_int_gw_ao(nspins), t_3c_overl_int_gw_ri(nspins), &
248 t_3c_overl_nnp_ic(nspins), t_3c_overl_nnp_ic_reflected(nspins), t_3c_o_mo_compressed(nspins))
250 t_3c_o_compressed, t_3c_o_ind, &
251 t_3c_overl_int_ao_mo, t_3c_o_mo_compressed(1), t_3c_o_mo_ind(1)%array, &
252 t_3c_overl_int_gw_ri(1), t_3c_overl_int_gw_ao(1), &
253 starts_array_mc, ends_array_mc, &
254 mo_coeff(1), matrix_s, &
255 gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), nmo, &
258 t_3c_overl_nnp_ic(1), t_3c_overl_nnp_ic_reflected(1), &
259 qs_env, unit_nr, do_alpha=.true.)
261 IF (my_open_shell)
THEN
264 t_3c_o_compressed, t_3c_o_ind, &
265 t_3c_overl_int_ao_mo_beta, t_3c_o_mo_compressed(2), t_3c_o_mo_ind(2)%array, &
266 t_3c_overl_int_gw_ri(2), t_3c_overl_int_gw_ao(2), &
267 starts_array_mc, ends_array_mc, &
268 mo_coeff(2), matrix_s, &
269 gw_corr_lev_occ(2), gw_corr_lev_virt(2), homo(2), nmo, &
272 t_3c_overl_nnp_ic(2), t_3c_overl_nnp_ic_reflected(2), &
273 qs_env, unit_nr, do_alpha=.false.)
275 IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma)
THEN
276 CALL dbt_destroy(t_3c_overl_int_ao_mo_beta)
281 ALLOCATE (fm_mat_w(num_integ_points))
283 DO jquad = 1, num_integ_points
285 CALL cp_fm_create(fm_mat_w(jquad), fm_mat_q%matrix_struct, set_zero=.true.)
292 template=matrix_s(1)%matrix, &
293 matrix_type=dbcsr_type_no_symmetry, &
294 row_blk_size=ri_blk_sizes, &
295 col_blk_size=ri_blk_sizes)
297 CALL timestop(handle)
341 gw_corr_lev_occ, gw_corr_lev_virt, homo, &
342 nmo, num_integ_group, num_integ_points, unit_nr, &
343 gw_corr_lev_tot, num_fit_points, omega_max_fit, &
344 do_minimax_quad, do_periodic, do_ri_Sigma_x, my_do_gw, &
345 first_cycle_periodic_correction, &
346 a_scaling, Eigenval, tj, vec_omega_fit_gw, vec_Sigma_x_gw, &
347 delta_corr, Eigenval_last, Eigenval_scf, vec_W_gw, &
348 fm_mat_S_gw, fm_mat_S_gw_work, &
349 para_env, mp2_env, kpoints, nkp, nkp_self_energy, &
350 do_kpoints_cubic_RPA, do_kpoints_from_Gamma)
352 COMPLEX(KIND=dp),
ALLOCATABLE, &
353 DIMENSION(:, :, :, :),
INTENT(OUT) :: vec_sigma_c_gw
354 INTEGER,
INTENT(IN) :: color_rpa_group, dimen_nm_gw
355 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
356 INTEGER,
INTENT(IN) :: nmo, num_integ_group, num_integ_points, &
358 INTEGER,
INTENT(INOUT) :: gw_corr_lev_tot, num_fit_points
359 REAL(kind=
dp) :: omega_max_fit
360 LOGICAL,
INTENT(IN) :: do_minimax_quad, do_periodic, &
361 do_ri_sigma_x, my_do_gw
362 LOGICAL,
INTENT(OUT) :: first_cycle_periodic_correction
363 REAL(kind=
dp),
INTENT(IN) :: a_scaling
364 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
365 INTENT(INOUT) :: eigenval
366 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
368 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
369 INTENT(OUT) :: vec_omega_fit_gw
370 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
371 INTENT(OUT) :: vec_sigma_x_gw
372 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
373 INTENT(INOUT) :: delta_corr
374 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
375 INTENT(OUT) :: eigenval_last, eigenval_scf
376 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
377 INTENT(OUT) :: vec_w_gw
378 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_s_gw
379 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
380 INTENT(INOUT) :: fm_mat_s_gw_work
384 INTEGER,
INTENT(OUT) :: nkp, nkp_self_energy
385 LOGICAL,
INTENT(IN) :: do_kpoints_cubic_rpa, &
386 do_kpoints_from_gamma
388 CHARACTER(LEN=*),
PARAMETER :: routinen =
'allocate_matrices_gw'
390 INTEGER :: handle, iquad, ispin, jquad, nspins
391 LOGICAL :: my_open_shell
392 REAL(kind=
dp) :: omega
393 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: vec_omega_gw
395 CALL timeset(routinen, handle)
397 nspins =
SIZE(eigenval, 3)
398 my_open_shell = (nspins == 2)
400 gw_corr_lev_tot = gw_corr_lev_occ(1) + gw_corr_lev_virt(1)
403 ALLOCATE (vec_omega_gw(num_integ_points))
404 vec_omega_gw = 0.0_dp
406 DO jquad = 1, num_integ_points
407 IF (do_minimax_quad)
THEN
410 omega = a_scaling/tan(tj(jquad))
412 vec_omega_gw(jquad) = omega
418 DO jquad = 1, num_integ_points
419 IF (vec_omega_gw(jquad) < omega_max_fit)
THEN
420 num_fit_points = num_fit_points + 1
425 IF (mp2_env%ri_g0w0%nparam_pade > num_fit_points)
THEN
426 IF (unit_nr > 0)
WRITE (unit=unit_nr, fmt=
"(T3,A)") &
427 "Pade approximation: more parameters than data points. Reset # of parameters."
428 mp2_env%ri_g0w0%nparam_pade = num_fit_points
429 IF (unit_nr > 0)
WRITE (unit=unit_nr, fmt=
"(T3,A,T74,I7)") &
430 "Number of pade parameters:", mp2_env%ri_g0w0%nparam_pade
435 ALLOCATE (vec_omega_fit_gw(num_fit_points))
439 DO jquad = 1, num_integ_points
440 IF (vec_omega_gw(jquad) < omega_max_fit)
THEN
442 vec_omega_fit_gw(iquad) = vec_omega_gw(jquad)
446 DEALLOCATE (vec_omega_gw)
448 IF (do_kpoints_cubic_rpa)
THEN
450 IF (mp2_env%ri_g0w0%do_gamma_only_sigma)
THEN
453 nkp_self_energy = nkp
455 ELSE IF (do_kpoints_from_gamma)
THEN
457 IF (mp2_env%ri_g0w0%do_kpoints_Sigma)
THEN
458 nkp_self_energy = mp2_env%ri_g0w0%nkp_self_energy
466 ALLOCATE (vec_sigma_c_gw(gw_corr_lev_tot, num_fit_points, nkp_self_energy, nspins))
469 ALLOCATE (eigenval_scf(nmo, nkp_self_energy, nspins))
470 eigenval_scf(:, :, :) = eigenval(:, :, :)
472 ALLOCATE (eigenval_last(nmo, nkp_self_energy, nspins))
473 eigenval_last(:, :, :) = eigenval(:, :, :)
475 IF (do_periodic)
THEN
477 ALLOCATE (delta_corr(1 + homo(1) - gw_corr_lev_occ(1):homo(1) + gw_corr_lev_virt(1)))
478 delta_corr(:) = 0.0_dp
480 first_cycle_periodic_correction = .true.
484 ALLOCATE (vec_sigma_x_gw(nmo, nkp_self_energy, nspins))
485 vec_sigma_x_gw = 0.0_dp
490 cpassert(.NOT. do_minimax_quad)
493 ALLOCATE (fm_mat_s_gw_work(nspins))
495 CALL cp_fm_create(fm_mat_s_gw_work(ispin), fm_mat_s_gw(ispin)%matrix_struct)
496 CALL cp_fm_set_all(matrix=fm_mat_s_gw_work(ispin), alpha=0.0_dp)
499 ALLOCATE (vec_w_gw(dimen_nm_gw, nspins))
503 IF (do_ri_sigma_x)
THEN
505 CALL get_vec_sigma_x(vec_sigma_x_gw(:, :, 1), nmo, fm_mat_s_gw(1), para_env, num_integ_group, color_rpa_group, &
506 homo(1), gw_corr_lev_occ(1), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, 1))
508 IF (my_open_shell)
THEN
509 CALL get_vec_sigma_x(vec_sigma_x_gw(:, :, 2), nmo, fm_mat_s_gw(2), para_env, num_integ_group, &
510 color_rpa_group, homo(2), gw_corr_lev_occ(2), &
511 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, 1))
518 CALL timestop(handle)
534 SUBROUTINE get_vec_sigma_x(vec_Sigma_x_gw, nmo, fm_mat_S_gw, para_env, num_integ_group, color_rpa_group, homo, &
535 gw_corr_lev_occ, vec_Sigma_x_minus_vxc_gw11)
537 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: vec_sigma_x_gw
538 INTEGER,
INTENT(IN) :: nmo
541 INTEGER,
INTENT(IN) :: num_integ_group, color_rpa_group, homo, &
543 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: vec_sigma_x_minus_vxc_gw11
545 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_vec_sigma_x'
547 INTEGER :: handle, iib, m_global, n_global, &
548 ncol_local, nm_global, nrow_local
549 INTEGER,
DIMENSION(:),
POINTER :: col_indices
551 CALL timeset(routinen, handle)
554 nrow_local=nrow_local, &
555 ncol_local=ncol_local, &
556 col_indices=col_indices)
561 DO iib = 1, ncol_local
564 IF (
modulo(1, num_integ_group) /= color_rpa_group) cycle
566 nm_global = col_indices(iib)
569 n_global = max(1, nm_global - 1)/nmo + 1
570 m_global = nm_global - (n_global - 1)*nmo
571 n_global = n_global + homo - gw_corr_lev_occ
573 IF (m_global <= homo)
THEN
576 vec_sigma_x_gw(n_global, 1) = &
577 vec_sigma_x_gw(n_global, 1) - &
578 dot_product(fm_mat_s_gw%local_data(:, iib), fm_mat_s_gw%local_data(:, iib))
586 CALL para_env%sum(vec_sigma_x_gw)
588 vec_sigma_x_minus_vxc_gw11(:) = &
589 vec_sigma_x_minus_vxc_gw11(:) + &
592 CALL timestop(handle)
594 END SUBROUTINE get_vec_sigma_x
613 vec_Sigma_x_minus_vxc_gw, Eigenval_last, &
614 Eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, kpoints, &
615 vec_Sigma_x_gw, my_do_gw)
617 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
618 INTENT(INOUT) :: fm_mat_s_gw_work
619 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
620 INTENT(INOUT) :: vec_w_gw
621 COMPLEX(KIND=dp),
ALLOCATABLE, &
622 DIMENSION(:, :, :, :),
INTENT(INOUT) :: vec_sigma_c_gw
623 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
624 INTENT(INOUT) :: vec_omega_fit_gw
625 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
626 INTENT(INOUT) :: vec_sigma_x_minus_vxc_gw, eigenval_last, &
628 LOGICAL,
INTENT(IN) :: do_periodic
629 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_re_mo_mo, &
630 matrix_berry_im_mo_mo
632 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
633 INTENT(INOUT) :: vec_sigma_x_gw
634 LOGICAL,
INTENT(IN) :: my_do_gw
636 CHARACTER(LEN=*),
PARAMETER :: routinen =
'deallocate_matrices_gw'
638 INTEGER :: handle, nspins
639 LOGICAL :: my_open_shell
641 CALL timeset(routinen, handle)
643 nspins =
SIZE(eigenval_last, 3)
644 my_open_shell = (nspins == 2)
648 DEALLOCATE (vec_sigma_x_minus_vxc_gw)
649 DEALLOCATE (vec_w_gw)
652 DEALLOCATE (vec_sigma_c_gw)
653 DEALLOCATE (vec_sigma_x_gw)
654 DEALLOCATE (vec_omega_fit_gw)
655 DEALLOCATE (eigenval_last)
656 DEALLOCATE (eigenval_scf)
658 IF (do_periodic)
THEN
664 CALL timestop(handle)
687 t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
688 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
689 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, mat_W, &
692 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
693 INTENT(INOUT) :: weights_cos_tf_w_to_t, &
694 weights_sin_tf_t_to_w
695 LOGICAL,
INTENT(IN) :: do_ic_model, do_kpoints_cubic_rpa
696 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
697 INTENT(INOUT) :: fm_mat_w
698 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_overl_int_ao_mo
700 DIMENSION(:) :: t_3c_o_mo_compressed
702 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:), &
703 INTENT(INOUT) :: t_3c_overl_int_gw_ri, &
704 t_3c_overl_int_gw_ao, &
706 t_3c_overl_nnp_ic_reflected
710 CHARACTER(LEN=*),
PARAMETER :: routinen =
'deallocate_matrices_gw_im_time'
712 INTEGER :: handle, ispin, nspins, unused
713 LOGICAL :: my_open_shell
715 CALL timeset(routinen, handle)
717 nspins =
SIZE(t_3c_overl_int_gw_ri)
718 my_open_shell = (nspins == 2)
720 IF (
ALLOCATED(weights_cos_tf_w_to_t))
DEALLOCATE (weights_cos_tf_w_to_t)
721 IF (
ALLOCATED(weights_sin_tf_t_to_w))
DEALLOCATE (weights_sin_tf_t_to_w)
723 IF (.NOT. do_kpoints_cubic_rpa)
THEN
729 CALL dbt_destroy(t_3c_overl_int_gw_ri(ispin))
730 CALL dbt_destroy(t_3c_overl_int_gw_ao(ispin))
732 DEALLOCATE (t_3c_overl_int_gw_ao, t_3c_overl_int_gw_ri)
733 IF (do_ic_model)
THEN
735 CALL dbt_destroy(t_3c_overl_nnp_ic(ispin))
736 CALL dbt_destroy(t_3c_overl_nnp_ic_reflected(ispin))
738 DEALLOCATE (t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected)
741 IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma)
THEN
743 DEALLOCATE (t_3c_o_mo_ind(ispin)%array)
746 DEALLOCATE (t_3c_o_mo_ind, t_3c_o_mo_compressed)
748 CALL dbt_destroy(t_3c_overl_int_ao_mo)
751 IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma)
THEN
753 CALL dbcsr_release(qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
754 DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
756 CALL dbcsr_release(qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
757 DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
759 DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc)
760 DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks)
763 CALL timestop(handle)
802 gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, &
803 do_im_time, do_periodic, &
804 first_cycle_periodic_correction, fermi_level_offset, &
805 omega, Eigenval, delta_corr, vec_omega_fit_gw, vec_W_gw, wj, &
806 fm_mat_Q, fm_mat_R_gw, fm_mat_S_gw, &
807 fm_mat_S_gw_work, mo_coeff, para_env, &
808 para_env_RPA, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, &
809 kpoints, qs_env, mp2_env)
811 COMPLEX(KIND=dp),
ALLOCATABLE, &
812 DIMENSION(:, :, :, :),
INTENT(INOUT) :: vec_sigma_c_gw
813 INTEGER,
INTENT(IN) :: dimen_nm_gw, dimen_ri
814 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
815 INTEGER,
INTENT(IN) :: jquad, nmo, num_fit_points
816 LOGICAL,
INTENT(IN) :: do_im_time, do_periodic
817 LOGICAL,
INTENT(INOUT) :: first_cycle_periodic_correction
818 REAL(kind=
dp),
INTENT(INOUT) :: fermi_level_offset, omega
819 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: eigenval
820 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
821 INTENT(INOUT) :: delta_corr
822 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
823 INTENT(IN) :: vec_omega_fit_gw
824 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
825 INTENT(INOUT) :: vec_w_gw
826 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
828 TYPE(
cp_fm_type),
INTENT(IN) :: fm_mat_q, fm_mat_r_gw
829 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_s_gw, fm_mat_s_gw_work
832 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_im_mo_mo, &
833 matrix_berry_re_mo_mo
838 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_GW_self_energy'
840 INTEGER :: handle, i_global, iib, ispin, j_global, &
841 jjb, ncol_local, nrow_local, nspins
842 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
844 CALL timeset(routinen, handle)
846 nspins =
SIZE(fm_mat_s_gw)
849 nrow_local=nrow_local, &
850 ncol_local=ncol_local, &
851 row_indices=row_indices, &
852 col_indices=col_indices)
854 IF (.NOT. do_im_time)
THEN
861 IF (do_periodic)
THEN
862 CALL calc_periodic_correction(delta_corr, qs_env, para_env, para_env_rpa, &
863 mp2_env%ri_g0w0%kp_grid, homo(1), nmo, gw_corr_lev_occ(1), &
864 gw_corr_lev_virt(1), omega, mo_coeff, eigenval(:, 1), &
865 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
866 first_cycle_periodic_correction, kpoints, &
867 mp2_env%ri_g0w0%do_mo_coeff_gamma, &
868 mp2_env%ri_g0w0%num_kp_grids, mp2_env%ri_g0w0%eps_kpoint, &
869 mp2_env%ri_g0w0%do_extra_kpoints, &
870 mp2_env%ri_g0w0%do_aux_bas_gw, mp2_env%ri_g0w0%frac_aux_mos)
873 CALL para_env_rpa%sync()
878 DO jjb = 1, ncol_local
879 j_global = col_indices(jjb)
880 DO iib = 1, nrow_local
881 i_global = row_indices(iib)
882 IF (j_global == i_global .AND. i_global <= dimen_ri)
THEN
883 fm_mat_q%local_data(iib, jjb) = fm_mat_q%local_data(iib, jjb) - 1.0_dp
888 CALL para_env_rpa%sync()
891 CALL compute_gw_self_energy_deep(vec_sigma_c_gw(:, :, :, ispin), dimen_nm_gw, dimen_ri, &
892 gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), &
893 homo(ispin), jquad, nmo, &
894 num_fit_points, do_periodic, fermi_level_offset, omega, &
895 eigenval(:, ispin), delta_corr, &
896 vec_omega_fit_gw, vec_w_gw(:, ispin), wj, fm_mat_q, &
897 fm_mat_s_gw(ispin), fm_mat_s_gw_work(ispin))
902 CALL timestop(handle)
915 REAL(kind=
dp),
INTENT(INOUT) :: fermi_level_offset
916 REAL(kind=
dp),
INTENT(IN) :: fermi_level_offset_input
917 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: eigenval
918 INTEGER,
DIMENSION(:),
INTENT(IN) :: homo
920 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_fermi_level_offset'
922 INTEGER :: handle, ispin, nspins
924 CALL timeset(routinen, handle)
926 nspins =
SIZE(eigenval, 2)
931 fermi_level_offset = fermi_level_offset_input
933 fermi_level_offset = min(fermi_level_offset, (eigenval(homo(ispin) + 1, ispin) - eigenval(homo(ispin), ispin))*0.5_dp)
936 CALL timestop(handle)
954 SUBROUTINE compute_w_cubic_gw(fm_mat_W, fm_mat_Q, fm_mat_work, dimen_RI, fm_mat_L, num_integ_points, &
955 tj, tau_tj, weights_cos_tf_w_to_t, jquad, omega)
956 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_w
957 TYPE(
cp_fm_type),
INTENT(IN) :: fm_mat_q, fm_mat_work
958 INTEGER,
INTENT(IN) :: dimen_ri
959 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: fm_mat_l
960 INTEGER,
INTENT(IN) :: num_integ_points
961 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
962 INTENT(IN) :: tj, tau_tj
963 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
964 INTENT(IN) :: weights_cos_tf_w_to_t
965 INTEGER,
INTENT(IN) :: jquad
966 REAL(kind=
dp),
INTENT(INOUT) :: omega
968 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_W_cubic_GW'
970 INTEGER :: handle, i_global, iib, iquad, j_global, &
971 jjb, ncol_local, nrow_local
972 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
973 REAL(kind=
dp) :: tau, weight
975 CALL timeset(routinen, handle)
978 nrow_local=nrow_local, &
979 ncol_local=ncol_local, &
980 row_indices=row_indices, &
981 col_indices=col_indices)
991 DO jjb = 1, ncol_local
992 j_global = col_indices(jjb)
993 DO iib = 1, nrow_local
994 i_global = row_indices(iib)
995 IF (j_global == i_global .AND. i_global <= dimen_ri)
THEN
996 fm_mat_q%local_data(iib, jjb) = fm_mat_q%local_data(iib, jjb) - 1.0_dp
1002 CALL parallel_gemm(
'T',
'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_mat_l(1, 1), fm_mat_q, &
1003 0.0_dp, fm_mat_work)
1005 CALL parallel_gemm(
'N',
'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_mat_work, fm_mat_l(1, 1), &
1009 DO iquad = 1, num_integ_points
1013 weight = weights_cos_tf_w_to_t(iquad, jquad)*cos(tau*omega)
1015 IF (jquad == 1)
THEN
1021 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_w(iquad), beta=weight, matrix_b=fm_mat_q)
1025 CALL timestop(handle)
1051 SUBROUTINE compute_gw_self_energy_deep(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, &
1052 gw_corr_lev_occ, gw_corr_lev_virt, &
1053 homo, jquad, nmo, num_fit_points, &
1054 do_periodic, fermi_level_offset, omega, Eigenval, &
1055 delta_corr, vec_omega_fit_gw, vec_W_gw, &
1056 wj, fm_mat_Q, fm_mat_S_gw, fm_mat_S_gw_work)
1058 COMPLEX(KIND=dp),
DIMENSION(:, :, :), &
1059 INTENT(INOUT) :: vec_sigma_c_gw
1060 INTEGER,
INTENT(IN) :: dimen_nm_gw, dimen_ri, gw_corr_lev_occ, &
1061 gw_corr_lev_virt, homo, jquad, nmo, &
1063 LOGICAL,
INTENT(IN) :: do_periodic
1064 REAL(kind=
dp),
INTENT(IN) :: fermi_level_offset
1065 REAL(kind=
dp),
INTENT(INOUT) :: omega
1066 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: eigenval
1067 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: delta_corr, vec_omega_fit_gw
1068 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: vec_w_gw
1069 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: wj
1070 TYPE(
cp_fm_type),
INTENT(IN) :: fm_mat_q, fm_mat_s_gw, fm_mat_s_gw_work
1072 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_GW_self_energy_deep'
1074 INTEGER :: handle, iib, iquad, m_global, n_global, &
1075 ncol_local, nm_global
1076 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
1077 REAL(kind=
dp) :: delta_corr_nn, e_fermi, omega_i, &
1080 CALL timeset(routinen, handle)
1083 CALL parallel_gemm(transa=
"N", transb=
"N", m=dimen_ri, n=dimen_nm_gw, k=dimen_ri, alpha=1.0_dp, &
1084 matrix_a=fm_mat_q, matrix_b=fm_mat_s_gw, beta=0.0_dp, &
1085 matrix_c=fm_mat_s_gw_work)
1088 ncol_local=ncol_local, &
1089 row_indices=row_indices, &
1090 col_indices=col_indices)
1096 DO iib = 1, ncol_local
1097 nm_global = col_indices(iib)
1098 vec_w_gw(nm_global) = vec_w_gw(nm_global) + &
1099 dot_product(fm_mat_s_gw_work%local_data(:, iib), fm_mat_s_gw%local_data(:, iib))
1102 n_global = max(1, nm_global - 1)/nmo + 1
1103 m_global = nm_global - (n_global - 1)*nmo
1104 n_global = n_global + homo - gw_corr_lev_occ
1107 DO iquad = 1, num_fit_points
1110 IF (n_global <= homo)
THEN
1111 sign_occ_virt = -1.0_dp
1113 sign_occ_virt = 1.0_dp
1116 omega_i = vec_omega_fit_gw(iquad)*sign_occ_virt
1120 IF (n_global <= homo)
THEN
1121 e_fermi = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset
1123 e_fermi = minval(eigenval(homo + 1:homo + gw_corr_lev_virt)) - fermi_level_offset
1127 IF (do_periodic .AND. row_indices(1) == 1 .AND. n_global == m_global)
THEN
1128 delta_corr_nn = delta_corr(n_global)
1130 delta_corr_nn = 0.0_dp
1136 vec_sigma_c_gw(n_global - homo + gw_corr_lev_occ, iquad, 1) = &
1137 vec_sigma_c_gw(n_global - homo + gw_corr_lev_occ, iquad, 1) - &
1138 0.5_dp/
pi*wj(jquad)/2.0_dp*(vec_w_gw(nm_global) + delta_corr_nn)* &
1139 (1.0_dp/(
gaussi*(omega + omega_i) + e_fermi - eigenval(m_global)) + &
1140 1.0_dp/(
gaussi*(-omega + omega_i) + e_fermi - eigenval(m_global)))
1145 CALL timestop(handle)
1147 END SUBROUTINE compute_gw_self_energy_deep
1211 gw_corr_lev_tot, gw_corr_lev_virt, homo, &
1212 nmo, num_fit_points, num_integ_points, &
1213 unit_nr, do_apply_ic_corr_to_gw, do_im_time, &
1214 do_periodic, do_ri_Sigma_x, &
1215 first_cycle_periodic_correction, e_fermi, eps_filter, &
1216 fermi_level_offset, delta_corr, Eigenval, &
1217 Eigenval_last, Eigenval_scf, iter_sc_GW0, exit_ev_gw, tau_tj, tj, &
1218 vec_omega_fit_gw, vec_Sigma_x_gw, ic_corr_list, &
1219 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, cfm_mo_coeff, mo_coeff, fm_mat_W, &
1220 para_env, para_env_RPA, mat_dm, mat_MinvVMinv, &
1221 t_3c_O, t_3c_M, t_3c_overl_int_ao_mo, &
1222 t_3c_O_compressed, t_3c_O_mo_compressed, &
1223 t_3c_O_ind, t_3c_O_mo_ind, &
1224 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, matrix_berry_im_mo_mo, &
1225 matrix_berry_re_mo_mo, mat_W, matrix_s, &
1226 kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_RPA, &
1227 starts_array_mc, ends_array_mc)
1229 COMPLEX(KIND=dp),
DIMENSION(:, :, :, :), &
1230 INTENT(OUT) :: vec_sigma_c_gw
1231 INTEGER,
INTENT(IN) :: count_ev_sc_gw
1232 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ
1233 INTEGER,
INTENT(IN) :: gw_corr_lev_tot
1234 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_virt, homo
1235 INTEGER,
INTENT(IN) :: nmo, num_fit_points, num_integ_points, &
1237 LOGICAL,
INTENT(IN) :: do_apply_ic_corr_to_gw, do_im_time, &
1238 do_periodic, do_ri_sigma_x
1239 LOGICAL,
INTENT(INOUT) :: first_cycle_periodic_correction
1240 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: e_fermi
1241 REAL(kind=
dp),
INTENT(IN) :: eps_filter, fermi_level_offset
1242 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
1243 INTENT(INOUT) :: delta_corr
1244 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: eigenval
1245 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
1246 INTENT(INOUT) :: eigenval_last, eigenval_scf
1247 INTEGER,
INTENT(IN) :: iter_sc_gw0
1248 LOGICAL,
INTENT(INOUT) :: exit_ev_gw
1249 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
1250 INTENT(INOUT) :: tau_tj, tj, vec_omega_fit_gw
1251 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
1252 INTENT(INOUT) :: vec_sigma_x_gw
1254 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
1255 INTENT(IN) :: weights_cos_tf_t_to_w, &
1256 weights_sin_tf_t_to_w
1257 TYPE(
cp_cfm_type),
DIMENSION(:),
INTENT(IN) :: cfm_mo_coeff
1259 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
1260 INTENT(IN) :: fm_mat_w
1262 TYPE(
dbcsr_p_type),
INTENT(IN) :: mat_dm, mat_minvvminv
1263 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_o
1264 TYPE(dbt_type) :: t_3c_m, t_3c_overl_int_ao_mo
1266 DIMENSION(:, :, :),
INTENT(INOUT) :: t_3c_o_compressed
1269 DIMENSION(:, :, :),
INTENT(INOUT) :: t_3c_o_ind
1271 TYPE(dbt_type),
DIMENSION(:) :: t_3c_overl_int_gw_ri, &
1272 t_3c_overl_int_gw_ao
1273 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_im_mo_mo, &
1274 matrix_berry_re_mo_mo
1280 INTEGER,
INTENT(IN) :: nkp_self_energy
1281 LOGICAL,
INTENT(IN) :: do_kpoints_cubic_rpa
1282 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc
1284 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_QP_energies'
1286 INTEGER :: count_ev_sc_gw_print, count_sc_gw0, count_sc_gw0_print, crossing_search, handle, &
1287 idos, ikp, ispin, iunit, n_level_gw, ndos, nspins, num_points_corr, num_poles
1288 LOGICAL :: do_kpoints_sigma, my_open_shell
1289 REAL(kind=
dp) :: dos_lower_bound, dos_precision, dos_upper_bound, e_cbm_gw, e_cbm_gw_beta, &
1290 e_cbm_scf, e_cbm_scf_beta, e_vbm_gw, e_vbm_gw_beta, e_vbm_scf, e_vbm_scf_beta, stop_crit
1291 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: vec_gw_dos
1292 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: m_value, vec_gw_energ, z_value
1296 CALL timeset(routinen, handle)
1299 my_open_shell = (nspins == 2)
1301 do_kpoints_sigma = mp2_env%ri_g0w0%do_kpoints_Sigma
1303 DO count_sc_gw0 = 1, iter_sc_gw0
1306 IF (do_im_time .AND. .NOT. do_kpoints_cubic_rpa .AND. .NOT. do_kpoints_sigma)
THEN
1307 num_points_corr = mp2_env%ri_g0w0%num_omega_points
1309 DO ispin = 1, nspins
1310 CALL compute_self_energy_cubic_gw(num_integ_points, nmo, tau_tj, tj, &
1311 matrix_s, cfm_mo_coeff(ispin), eigenval(:, 1, ispin), eps_filter, &
1312 e_fermi(ispin), fm_mat_w, &
1313 gw_corr_lev_tot, gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), homo(ispin), &
1314 count_ev_sc_gw, count_sc_gw0, &
1315 t_3c_overl_int_ao_mo, t_3c_o_mo_compressed(ispin), &
1316 t_3c_o_mo_ind(ispin)%array, &
1317 t_3c_overl_int_gw_ri(ispin), t_3c_overl_int_gw_ao(ispin), &
1318 mat_w, mat_minvvminv, mat_dm, &
1319 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_sigma_c_gw(:, :, :, ispin), &
1320 do_periodic, num_points_corr, delta_corr, qs_env, para_env, para_env_rpa, &
1321 mp2_env, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
1322 first_cycle_periodic_correction, kpoints, num_fit_points, mo_coeff, &
1323 do_ri_sigma_x, vec_sigma_x_gw(:, :, ispin), unit_nr, ispin)
1328 IF (do_kpoints_sigma)
THEN
1329 CALL compute_self_energy_cubic_gw_kpoints(num_integ_points, tau_tj, tj, &
1330 matrix_s, eigenval(:, :, :), e_fermi, fm_mat_w, &
1331 gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, &
1332 count_ev_sc_gw, count_sc_gw0, &
1333 t_3c_o, t_3c_m, t_3c_o_compressed, t_3c_o_ind, &
1334 mat_w, mat_minvvminv, &
1335 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_sigma_c_gw(:, :, :, :), &
1337 mp2_env, num_fit_points, mo_coeff, &
1338 do_ri_sigma_x, vec_sigma_x_gw(:, :, :), unit_nr, nspins, &
1339 starts_array_mc, ends_array_mc, eps_filter)
1343 IF (do_periodic .AND. mp2_env%ri_g0w0%do_average_deg_levels)
THEN
1345 DO ispin = 1, nspins
1346 CALL average_degenerate_levels(vec_sigma_c_gw(:, :, :, ispin), &
1347 eigenval(1 + homo(ispin) - gw_corr_lev_occ(ispin): &
1348 homo(ispin) + gw_corr_lev_virt(ispin), 1, ispin), &
1349 mp2_env%ri_g0w0%eps_eigenval)
1353 IF (.NOT. do_im_time)
THEN
1354 CALL para_env%sum(vec_sigma_c_gw)
1357 CALL para_env%sync()
1360 num_poles = mp2_env%ri_g0w0%num_poles
1361 crossing_search = mp2_env%ri_g0w0%crossing_search
1364 ALLOCATE (vec_gw_energ(gw_corr_lev_tot, nkp_self_energy, nspins))
1365 vec_gw_energ = 0.0_dp
1366 ALLOCATE (z_value(gw_corr_lev_tot, nkp_self_energy, nspins))
1368 ALLOCATE (m_value(gw_corr_lev_tot, nkp_self_energy, nspins))
1374 e_vbm_gw_beta = -1.0e3
1375 e_cbm_gw_beta = 1.0e3
1376 e_vbm_scf_beta = -1.0e3
1377 e_cbm_scf_beta = 1.0e3
1380 dos_precision = mp2_env%ri_g0w0%dos_prec
1381 dos_upper_bound = mp2_env%ri_g0w0%dos_upper
1382 dos_lower_bound = mp2_env%ri_g0w0%dos_lower
1384 IF (dos_lower_bound >= dos_upper_bound)
THEN
1385 CALL cp_abort(__location__,
"Invalid settings for GW_DOS calculation!")
1388 IF (dos_precision /= 0)
THEN
1389 ndos = int((dos_upper_bound - dos_lower_bound)/dos_precision)
1390 ALLOCATE (vec_gw_dos(ndos))
1395 DO ikp = 1, nkp_self_energy
1397 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
1400 DO n_level_gw = 1, gw_corr_lev_tot
1402 IF (
modulo(n_level_gw, para_env%num_pe) /= para_env%mepos) cycle
1404 SELECT CASE (mp2_env%ri_g0w0%analytic_continuation)
1406 CALL fit_and_continuation_2pole(vec_gw_energ(:, ikp, 1), vec_omega_fit_gw, &
1407 z_value(:, ikp, 1), m_value(:, ikp, 1), vec_sigma_c_gw(:, :, ikp, 1), &
1408 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), &
1409 eigenval(:, ikp, 1), eigenval_scf(:, ikp, 1), n_level_gw, &
1410 gw_corr_lev_occ(1), gw_corr_lev_virt(1), num_poles, &
1411 num_fit_points, crossing_search, homo(1), stop_crit, &
1412 fermi_level_offset, do_im_time)
1416 z_value(:, ikp, 1), m_value(:, ikp, 1), vec_sigma_c_gw(:, :, ikp, 1), &
1417 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), &
1418 eigenval(:, ikp, 1), eigenval_scf(:, ikp, 1), &
1419 mp2_env%ri_g0w0%do_hedin_shift, n_level_gw, &
1420 gw_corr_lev_occ(1), gw_corr_lev_virt(1), mp2_env%ri_g0w0%nparam_pade, &
1421 num_fit_points, crossing_search, homo(1), fermi_level_offset, &
1422 do_im_time, mp2_env%ri_g0w0%print_self_energy, count_ev_sc_gw, &
1423 vec_gw_dos, dos_lower_bound, dos_precision, ndos, &
1424 mp2_env%ri_g0w0%min_level_self_energy, &
1425 mp2_env%ri_g0w0%max_level_self_energy, mp2_env%ri_g0w0%dos_eta, &
1426 mp2_env%ri_g0w0%dos_min, mp2_env%ri_g0w0%dos_max)
1428 cpabort(
"Only two-model and Pade approximation are implemented.")
1431 IF (my_open_shell)
THEN
1432 SELECT CASE (mp2_env%ri_g0w0%analytic_continuation)
1434 CALL fit_and_continuation_2pole( &
1435 vec_gw_energ(:, ikp, 2), vec_omega_fit_gw, &
1436 z_value(:, ikp, 2), m_value(:, ikp, 2), vec_sigma_c_gw(:, :, ikp, 2), &
1437 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), &
1438 eigenval(:, ikp, 2), eigenval_scf(:, ikp, 2), n_level_gw, &
1439 gw_corr_lev_occ(2), gw_corr_lev_virt(2), num_poles, &
1440 num_fit_points, crossing_search, homo(2), stop_crit, &
1441 fermi_level_offset, do_im_time)
1444 z_value(:, ikp, 2), m_value(:, ikp, 2), vec_sigma_c_gw(:, :, ikp, 2), &
1445 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), &
1446 eigenval(:, ikp, 2), eigenval_scf(:, ikp, 2), &
1447 mp2_env%ri_g0w0%do_hedin_shift, n_level_gw, &
1448 gw_corr_lev_occ(2), gw_corr_lev_virt(2), mp2_env%ri_g0w0%nparam_pade, &
1449 num_fit_points, crossing_search, homo(2), &
1450 fermi_level_offset, do_im_time, &
1451 mp2_env%ri_g0w0%print_self_energy, count_ev_sc_gw, &
1452 vec_gw_dos, dos_lower_bound, dos_precision, ndos, &
1453 mp2_env%ri_g0w0%min_level_self_energy, &
1454 mp2_env%ri_g0w0%max_level_self_energy, mp2_env%ri_g0w0%dos_eta, &
1455 mp2_env%ri_g0w0%dos_min, mp2_env%ri_g0w0%dos_max)
1457 cpabort(
"Only two-pole model and Pade approximation are implemented.")
1464 CALL para_env%sum(vec_gw_energ)
1465 CALL para_env%sum(z_value)
1466 CALL para_env%sum(m_value)
1468 IF (dos_precision /= 0.0_dp)
THEN
1469 CALL para_env%sum(vec_gw_dos)
1472 CALL check_nan(vec_gw_energ, 0.0_dp)
1473 CALL check_nan(z_value, 1.0_dp)
1474 CALL check_nan(m_value, 0.0_dp)
1476 IF (do_im_time .OR. mp2_env%ri_g0w0%iter_sc_GW0 == 1)
THEN
1477 count_ev_sc_gw_print = count_ev_sc_gw
1478 count_sc_gw0_print = count_sc_gw0
1480 count_ev_sc_gw_print = count_sc_gw0
1481 count_sc_gw0_print = count_ev_sc_gw
1485 IF (my_open_shell)
THEN
1487 CALL print_and_update_for_ev_sc( &
1488 vec_gw_energ(:, ikp, 1), &
1489 z_value(:, ikp, 1), m_value(:, ikp, 1), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), &
1490 eigenval(:, ikp, 1), eigenval_last(:, ikp, 1), eigenval_scf(:, ikp, 1), &
1491 gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, &
1492 crossing_search, homo(1), unit_nr, count_ev_sc_gw_print, count_sc_gw0_print, &
1493 ikp, nkp_self_energy, kpoints_sigma, 1, e_vbm_gw, e_cbm_gw, e_vbm_scf, e_cbm_scf)
1495 CALL print_and_update_for_ev_sc( &
1496 vec_gw_energ(:, ikp, 2), &
1497 z_value(:, ikp, 2), m_value(:, ikp, 2), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), &
1498 eigenval(:, ikp, 2), eigenval_last(:, ikp, 2), eigenval_scf(:, ikp, 2), &
1499 gw_corr_lev_occ(2), gw_corr_lev_virt(2), gw_corr_lev_tot, &
1500 crossing_search, homo(2), unit_nr, count_ev_sc_gw_print, count_sc_gw0_print, &
1501 ikp, nkp_self_energy, kpoints_sigma, 2, e_vbm_gw_beta, e_cbm_gw_beta, e_vbm_scf_beta, e_cbm_scf_beta)
1503 IF (do_apply_ic_corr_to_gw .AND. count_ev_sc_gw == 1)
THEN
1505 CALL apply_ic_corr(eigenval(:, ikp, 1), eigenval_scf(:, ikp, 1), ic_corr_list(1)%array, &
1506 gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, &
1507 homo(1), nmo, unit_nr, do_alpha=.true.)
1509 CALL apply_ic_corr(eigenval(:, ikp, 2), eigenval_scf(:, ikp, 2), ic_corr_list(2)%array, &
1510 gw_corr_lev_occ(2), gw_corr_lev_virt(2), gw_corr_lev_tot, &
1511 homo(2), nmo, unit_nr, do_beta=.true.)
1517 CALL print_and_update_for_ev_sc( &
1518 vec_gw_energ(:, ikp, 1), &
1519 z_value(:, ikp, 1), m_value(:, ikp, 1), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), &
1520 eigenval(:, ikp, 1), eigenval_last(:, ikp, 1), eigenval_scf(:, ikp, 1), &
1521 gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, &
1522 crossing_search, homo(1), unit_nr, count_ev_sc_gw_print, count_sc_gw0_print, &
1523 ikp, nkp_self_energy, kpoints_sigma, 0, e_vbm_gw, e_cbm_gw, e_vbm_scf, e_cbm_scf)
1525 IF (do_apply_ic_corr_to_gw .AND. count_ev_sc_gw == 1)
THEN
1527 CALL apply_ic_corr(eigenval(:, ikp, 1), eigenval_scf(:, ikp, 1), ic_corr_list(1)%array, &
1528 gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, &
1529 homo(1), nmo, unit_nr)
1537 IF (nkp_self_energy > 1 .AND. unit_nr > 0)
THEN
1539 CALL print_gaps(e_vbm_scf, e_cbm_scf, e_vbm_scf_beta, e_cbm_scf_beta, &
1540 e_vbm_gw, e_cbm_gw, e_vbm_gw_beta, e_cbm_gw_beta, my_open_shell, unit_nr)
1546 IF (mp2_env%ri_g0w0%soc_type /=
soc_none)
THEN
1547 CALL calculate_and_print_soc(qs_env, eigenval_scf, eigenval_scf, gw_corr_lev_occ, gw_corr_lev_virt, &
1548 homo, unit_nr, do_soc_gw=.false., do_soc_scf=.true.)
1549 CALL calculate_and_print_soc(qs_env, eigenval, eigenval_scf, gw_corr_lev_occ, gw_corr_lev_virt, &
1550 homo, unit_nr, do_soc_gw=.true., do_soc_scf=.false.)
1554 IF (logger%para_env%is_source())
THEN
1560 IF (dos_precision /= 0.0_dp)
THEN
1562 CALL open_file(
'spectral.dat', unit_number=iunit, file_status=
"UNKNOWN", file_action=
"WRITE")
1566 WRITE (iunit,
'(E17.10, E17.10)') (dos_lower_bound + real(idos - 1, kind=
dp)*dos_precision)*
evolt, &
1571 DEALLOCATE (vec_gw_dos)
1574 DEALLOCATE (z_value)
1575 DEALLOCATE (m_value)
1576 DEALLOCATE (vec_gw_energ)
1578 exit_ev_gw = .false.
1581 IF (abs(eigenval(homo(1), 1, 1) - eigenval_last(homo(1), 1, 1) - &
1582 eigenval(homo(1) + 1, 1, 1) + eigenval_last(homo(1) + 1, 1, 1)) &
1583 < mp2_env%ri_g0w0%eps_iter)
THEN
1584 IF (count_sc_gw0 == 1) exit_ev_gw = .true.
1588 DO ispin = 1, nspins
1589 CALL shift_unshifted_levels(eigenval(:, 1, ispin), eigenval_last(:, 1, ispin), gw_corr_lev_occ(ispin), &
1590 gw_corr_lev_virt(ispin), homo(ispin), nmo)
1593 IF (do_im_time .AND. do_kpoints_sigma .AND. mp2_env%ri_g0w0%print_local_bandgap)
THEN
1594 CALL print_local_bandgap(qs_env, eigenval, gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1),
"GW")
1595 CALL print_local_bandgap(qs_env, eigenval_scf, gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1),
"DFT")
1599 IF (.NOT. do_im_time)
EXIT
1603 CALL timestop(handle)
1619 SUBROUTINE calculate_and_print_soc(qs_env, Eigenval, Eigenval_scf, gw_corr_lev_occ, gw_corr_lev_virt, &
1620 homo, unit_nr, do_soc_gw, do_soc_scf)
1622 REAL(kind=
dp),
DIMENSION(:, :, :) :: eigenval, eigenval_scf
1623 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
1625 LOGICAL :: do_soc_gw, do_soc_scf
1627 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_and_print_soc'
1629 INTEGER :: handle, i_dim, i_glob, i_row, ikp, j_col, j_glob, n_level_gw, nao, ncol_local, &
1630 nder, nkind, nkp_self_energy, nrow_local, periodic(3), size_real_space
1631 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: index0
1632 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
1633 LOGICAL :: calculate_forces, use_virial
1634 REAL(kind=
dp) :: avg_occ_qp_shift, avg_virt_qp_shift, e_cbm_gw_soc, e_gap_gw_soc, e_homo, &
1635 e_homo_gw_soc, e_i, e_j, e_lumo, e_lumo_gw_soc, e_vbm_gw_soc, e_window, eps_ppnl
1636 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eigenvalues_without_soc_sorted
1637 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eigenvalues
1640 TYPE(
cp_cfm_type) :: cfm_mat_h_double, cfm_mat_h_ks, &
1641 cfm_mat_s_double, cfm_mat_work_double, &
1642 cfm_mo_coeff, cfm_mo_coeff_double
1644 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, matrix_s_desymm
1645 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_vsoc_l_nosymm, mat_vsoc_lx_kp, &
1646 mat_vsoc_ly_kp, mat_vsoc_lz_kp, &
1647 matrix_dummy, matrix_l, &
1653 POINTER :: sab_orb, sap_ppnl
1656 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1660 CALL timeset(routinen, handle)
1662 cpassert(do_soc_gw .NEQV. do_soc_scf)
1665 matrix_s=matrix_s, &
1666 para_env=para_env, &
1667 qs_kind_set=qs_kind_set, &
1669 atomic_kind_set=atomic_kind_set, &
1670 particle_set=particle_set, &
1671 sap_ppnl=sap_ppnl, &
1672 dft_control=dft_control, &
1675 scf_control=scf_control)
1677 calculate_forces = .false.
1678 use_virial = .false.
1680 eps_ppnl = dft_control%qs_control%eps_ppnl
1682 CALL get_cell(cell=cell, periodic=periodic)
1684 size_real_space = 3**(periodic(1) + periodic(2) + periodic(3))
1689 ALLOCATE (matrix_l(i_dim, 1)%matrix)
1690 CALL dbcsr_create(matrix_l(i_dim, 1)%matrix, template=matrix_s(1)%matrix, &
1691 matrix_type=dbcsr_type_antisymmetric)
1693 CALL dbcsr_set(matrix_l(i_dim, 1)%matrix, 0.0_dp)
1696 NULLIFY (matrix_pot_dummy)
1698 ALLOCATE (matrix_pot_dummy(1, 1)%matrix)
1699 CALL dbcsr_create(matrix_pot_dummy(1, 1)%matrix, template=matrix_s(1)%matrix)
1701 CALL dbcsr_set(matrix_pot_dummy(1, 1)%matrix, 0.0_dp)
1703 CALL build_core_ppnl(matrix_pot_dummy, matrix_dummy, force, virial, calculate_forces, use_virial, nder, &
1704 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
1705 nimages=1, basis_type=
"ORB", matrix_l=matrix_l)
1707 CALL alloc_mat_set_2d(mat_vsoc_l_nosymm, 3, size_real_space, matrix_s(1)%matrix, explicitly_no_symmetry=.true.)
1709 CALL dbcsr_desymmetrize(matrix_l(i_dim, 1)%matrix, mat_vsoc_l_nosymm(i_dim, 1)%matrix)
1712 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
1714 CALL mat_kp_from_mat_gamma(qs_env, mat_vsoc_lx_kp, mat_vsoc_l_nosymm(1, 1)%matrix, kpoints_sigma, 1, .false.)
1715 CALL mat_kp_from_mat_gamma(qs_env, mat_vsoc_ly_kp, mat_vsoc_l_nosymm(2, 1)%matrix, kpoints_sigma, 1, .false.)
1716 CALL mat_kp_from_mat_gamma(qs_env, mat_vsoc_lz_kp, mat_vsoc_l_nosymm(3, 1)%matrix, kpoints_sigma, 1, .false.)
1718 nkp_self_energy = kpoints_sigma%nkp
1720 CALL get_mo_set(kpoints_sigma%kp_env(1)%kpoint_env%mos(1, 1), mo_coeff=rmos)
1722 CALL create_cfm_double_row_col_size(rmos, cfm_mat_h_double)
1723 CALL create_cfm_double_row_col_size(rmos, cfm_mat_s_double)
1724 CALL create_cfm_double_row_col_size(rmos, cfm_mo_coeff_double)
1725 CALL create_cfm_double_row_col_size(rmos, cfm_mat_work_double)
1734 NULLIFY (matrix_s_desymm)
1736 ALLOCATE (matrix_s_desymm(1)%matrix)
1737 CALL dbcsr_create(matrix=matrix_s_desymm(1)%matrix, template=matrix_s(1)%matrix, &
1738 matrix_type=dbcsr_type_no_symmetry)
1741 ALLOCATE (eigenvalues(2*nao))
1742 eigenvalues = 0.0_dp
1743 ALLOCATE (eigenvalues_without_soc_sorted(2*nao))
1745 e_window = qs_env%mp2_env%ri_g0w0%soc_energy_window
1746 IF (unit_nr > 0)
THEN
1747 WRITE (unit_nr,
'(T3,A)')
' '
1748 WRITE (unit_nr,
'(T3,A)')
'------------------------------------------------------------------------------'
1749 WRITE (unit_nr,
'(T3,A)')
' '
1750 WRITE (unit_nr,
'(T3,A,F42.1)')
'GW_SOC_INFO | SOC energy window (eV)', e_window*
evolt
1753 e_vbm_gw_soc = -1000.0_dp
1754 e_cbm_gw_soc = 1000.0_dp
1756 DO ikp = 1, nkp_self_energy
1758 CALL get_mo_set(kpoints_sigma%kp_env(ikp)%kpoint_env%mos(1, 1), mo_coeff=rmos)
1759 CALL get_mo_set(kpoints_sigma%kp_env(ikp)%kpoint_env%mos(2, 1), mo_coeff=imos)
1763 avg_occ_qp_shift = sum(eigenval(homo(1) - gw_corr_lev_occ(1) + 1:homo(1), ikp, 1) - &
1764 eigenval_scf(homo(1) - gw_corr_lev_occ(1) + 1:homo(1), ikp, 1))/gw_corr_lev_occ(1)
1765 avg_virt_qp_shift = sum(eigenval(homo(1):homo(1) + gw_corr_lev_virt(1), ikp, 1) - &
1766 eigenval_scf(homo(1):homo(1) + gw_corr_lev_virt(1), ikp, 1))/gw_corr_lev_virt(1)
1768 IF (gw_corr_lev_occ(1) < homo(1))
THEN
1769 eigenval(1:homo(1) - gw_corr_lev_occ(1), ikp, 1) = eigenval_scf(1:homo(1) - gw_corr_lev_occ(1), ikp, 1) &
1772 IF (gw_corr_lev_virt(1) < nao - homo(1) + 1)
THEN
1773 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) &
1778 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.)
1779 CALL add_dbcsr_submatrix(cfm_mat_h_double, mat_vsoc_ly_kp(ikp, 1:2), cfm_mat_h_ks, nao + 1, 1,
gaussi, .true.)
1780 CALL add_dbcsr_submatrix(cfm_mat_h_double, mat_vsoc_lz_kp(ikp, 1:2), cfm_mat_h_ks, 1, 1,
z_one, .false.)
1781 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.)
1784 cfm_mo_coeff_double%local_data =
z_zero
1785 CALL add_cfm_submatrix(cfm_mo_coeff_double, cfm_mo_coeff, 1, 1)
1786 CALL add_cfm_submatrix(cfm_mo_coeff_double, cfm_mo_coeff, nao + 1, nao + 1)
1789 nrow_local=nrow_local, &
1790 ncol_local=ncol_local, &
1791 row_indices=row_indices, &
1792 col_indices=col_indices)
1795 matrix_a=cfm_mat_h_double, matrix_b=cfm_mo_coeff_double, beta=
z_zero, &
1796 matrix_c=cfm_mat_work_double)
1799 matrix_a=cfm_mo_coeff_double, matrix_b=cfm_mat_work_double, beta=
z_zero, &
1800 matrix_c=cfm_mat_h_double)
1803 nrow_local=nrow_local, &
1804 ncol_local=ncol_local, &
1805 row_indices=row_indices, &
1806 col_indices=col_indices)
1810 e_homo = eigenval(homo(1), ikp, 1)
1811 e_lumo = eigenval(homo(1) + 1, ikp, 1)
1813 CALL para_env%sync()
1815 DO i_row = 1, nrow_local
1816 DO j_col = 1, ncol_local
1817 i_glob = row_indices(i_row)
1818 j_glob = col_indices(j_col)
1819 IF (i_glob <= nao)
THEN
1820 e_i = eigenval(i_glob, ikp, 1)
1822 e_i = eigenval(i_glob - nao, ikp, 1)
1824 IF (j_glob <= nao)
THEN
1825 e_j = eigenval(j_glob, ikp, 1)
1827 e_j = eigenval(j_glob - nao, ikp, 1)
1831 IF (i_glob == j_glob)
THEN
1832 cfm_mat_h_double%local_data(i_row, j_col) = cfm_mat_h_double%local_data(i_row, j_col) + e_i*
z_one
1833 cfm_mat_s_double%local_data(i_row, j_col) =
z_one
1835 IF (e_i < e_homo - 0.5_dp*e_window .OR. e_i > e_lumo + 0.5_dp*e_window .OR. &
1836 e_j < e_homo - 0.5_dp*e_window .OR. e_j > e_lumo + 0.5_dp*e_window)
THEN
1837 cfm_mat_h_double%local_data(i_row, j_col) =
z_zero
1844 CALL para_env%sync()
1846 eigenvalues = 0.0_dp
1847 CALL cp_cfm_geeig_canon(cfm_mat_h_double, cfm_mat_s_double, cfm_mo_coeff_double, eigenvalues, &
1848 cfm_mat_work_double, scf_control%eps_eigval)
1850 eigenvalues_without_soc_sorted(1:nao) = eigenval(:, ikp, 1)
1851 eigenvalues_without_soc_sorted(nao + 1:2*nao) = eigenval(:, ikp, 1)
1852 ALLOCATE (index0(2*nao))
1853 CALL sort(eigenvalues_without_soc_sorted, 2*nao, index0)
1856 e_homo_gw_soc = maxval(eigenvalues(2*homo(1) - 2*gw_corr_lev_occ(1) + 1:2*homo(1)))
1857 e_lumo_gw_soc = minval(eigenvalues(2*homo(1) + 1:2*homo(1) + 2*gw_corr_lev_virt(1)))
1858 e_gap_gw_soc = e_lumo_gw_soc - e_homo_gw_soc
1859 IF (e_homo_gw_soc > e_vbm_gw_soc) e_vbm_gw_soc = e_homo_gw_soc
1860 IF (e_lumo_gw_soc < e_cbm_gw_soc) e_cbm_gw_soc = e_lumo_gw_soc
1862 IF (unit_nr > 0)
THEN
1863 WRITE (unit_nr,
'(T3,A)')
' '
1864 WRITE (unit_nr,
'(T3,A7,I3,A3,I3,A8,3F7.3,A12,3F7.3)')
'Kpoint ', ikp,
' /', nkp_self_energy, &
1865 ' xkp =', kpoints_sigma%xkp(1, ikp), kpoints_sigma%xkp(2, ikp), kpoints_sigma%xkp(3, ikp), &
1866 ' and xkp =', -kpoints_sigma%xkp(1, ikp), -kpoints_sigma%xkp(2, ikp), -kpoints_sigma%xkp(3, ikp)
1867 WRITE (unit_nr,
'(T3,A)')
' '
1869 WRITE (unit_nr,
'(T3,A)')
' '
1870 WRITE (unit_nr,
'(T3,A,F13.4)')
'GW_SOC_INFO | Average GW shift of occupied levels compared to SCF', &
1871 avg_occ_qp_shift*
evolt
1872 WRITE (unit_nr,
'(T3,A,F11.4)')
'GW_SOC_INFO | Average GW shift of unoccupied levels compared to SCF', &
1873 avg_virt_qp_shift*
evolt
1874 WRITE (unit_nr,
'(T3,A)')
' '
1875 WRITE (unit_nr,
'(T3,2A)')
'Molecular orbital E_GW with SOC (eV) E_GW without SOC (eV) SOC shift (eV)'
1877 WRITE (unit_nr,
'(T3,2A)')
'Molecular orbital E_SCF with SOC (eV) E_SCF without SOC (eV) SOC shift (eV)'
1880 DO n_level_gw = 2*(homo(1) - gw_corr_lev_occ(1)) + 1, 2*homo(1)
1881 WRITE (unit_nr,
'(T3,I4,A,3F21.4)') n_level_gw,
' ( occ ) ', eigenvalues(n_level_gw)*
evolt, &
1882 eigenvalues_without_soc_sorted(n_level_gw)*
evolt, &
1883 (eigenvalues(n_level_gw) - eigenvalues_without_soc_sorted(n_level_gw))*
evolt
1885 DO n_level_gw = 2*homo(1) + 1, 2*(homo(1) + gw_corr_lev_virt(1))
1886 WRITE (unit_nr,
'(T3,I4,A,3F21.4)') n_level_gw,
' ( vir ) ', eigenvalues(n_level_gw)*
evolt, &
1887 eigenvalues_without_soc_sorted(n_level_gw)*
evolt, &
1888 (eigenvalues(n_level_gw) - eigenvalues_without_soc_sorted(n_level_gw))*
evolt
1890 WRITE (unit_nr,
'(T3,A)')
' '
1892 WRITE (unit_nr,
'(T3,A,F38.4)')
'GW+SOC direct gap at current kpoint (eV)', e_gap_gw_soc*
evolt
1894 WRITE (unit_nr,
'(T3,A,F37.4)')
'SCF+SOC direct gap at current kpoint (eV)', e_gap_gw_soc*
evolt
1896 WRITE (unit_nr,
'(T3,A)')
' '
1897 WRITE (unit_nr,
'(T3,A)')
'------------------------------------------------------------------------------'
1902 IF (unit_nr > 0)
THEN
1903 WRITE (unit_nr,
'(T3,A)')
' '
1905 WRITE (unit_nr,
'(T3,A,F46.4)')
'GW+SOC valence band maximum (eV)', e_vbm_gw_soc*
evolt
1906 WRITE (unit_nr,
'(T3,A,F43.4)')
'GW+SOC conduction band minimum (eV)', e_cbm_gw_soc*
evolt
1907 WRITE (unit_nr,
'(T3,A,F59.4)')
'GW+SOC bandgap (eV)', (e_cbm_gw_soc - e_vbm_gw_soc)*
evolt
1909 WRITE (unit_nr,
'(T3,A,F45.4)')
'SCF+SOC valence band maximum (eV)', e_vbm_gw_soc*
evolt
1910 WRITE (unit_nr,
'(T3,A,F42.4)')
'SCF+SOC conduction band minimum (eV)', e_cbm_gw_soc*
evolt
1911 WRITE (unit_nr,
'(T3,A,F58.4)')
'SCF+SOC bandgap (eV)', (e_cbm_gw_soc - e_vbm_gw_soc)*
evolt
1929 DEALLOCATE (eigenvalues)
1931 CALL timestop(handle)
1933 END SUBROUTINE calculate_and_print_soc
1945 SUBROUTINE add_dbcsr_submatrix(cfm_mat_target, mat_source, cfm_source_template, &
1946 nstart_row, nstart_col, factor, add_also_herm_conj)
1950 INTEGER :: nstart_row, nstart_col
1951 COMPLEX(KIND=dp) :: factor
1952 LOGICAL :: add_also_herm_conj
1954 CHARACTER(LEN=*),
PARAMETER :: routinen =
'add_dbcsr_submatrix'
1956 INTEGER :: handle, nao
1958 cfm_mat_work_double_2
1960 fm_mat_work_double_re, fm_mat_work_im, &
1963 CALL timeset(routinen, handle)
1965 CALL cp_fm_create(fm_mat_work_double_re, cfm_mat_target%matrix_struct)
1966 CALL cp_fm_create(fm_mat_work_double_im, cfm_mat_target%matrix_struct)
1970 CALL cp_cfm_create(cfm_mat_work_double, cfm_mat_target%matrix_struct)
1971 CALL cp_cfm_create(cfm_mat_work_double_2, cfm_mat_target%matrix_struct)
1975 CALL cp_fm_create(fm_mat_work_re, cfm_source_template%matrix_struct)
1976 CALL cp_fm_create(fm_mat_work_im, cfm_source_template%matrix_struct)
1984 nrow=nao, ncol=nao, &
1985 s_firstrow=1, s_firstcol=1, &
1986 t_firstrow=nstart_row, t_firstcol=nstart_col)
1989 nrow=nao, ncol=nao, &
1990 s_firstrow=1, s_firstcol=1, &
1991 t_firstrow=nstart_row, t_firstcol=nstart_col)
2000 IF (add_also_herm_conj)
THEN
2012 CALL timestop(handle)
2014 END SUBROUTINE add_dbcsr_submatrix
2023 SUBROUTINE add_cfm_submatrix(cfm_mat_target, cfm_mat_source, nstart_row, nstart_col)
2025 TYPE(
cp_cfm_type) :: cfm_mat_target, cfm_mat_source
2026 INTEGER :: nstart_row, nstart_col
2028 CHARACTER(LEN=*),
PARAMETER :: routinen =
'add_cfm_submatrix'
2030 INTEGER :: handle, nao
2032 fm_mat_work_double_re, fm_mat_work_im, &
2035 CALL timeset(routinen, handle)
2037 CALL cp_fm_create(fm_mat_work_double_re, cfm_mat_target%matrix_struct)
2038 CALL cp_fm_create(fm_mat_work_double_im, cfm_mat_target%matrix_struct)
2042 CALL cp_fm_create(fm_mat_work_re, cfm_mat_source%matrix_struct)
2043 CALL cp_fm_create(fm_mat_work_im, cfm_mat_source%matrix_struct)
2044 CALL cp_cfm_to_fm(cfm_mat_source, fm_mat_work_re, fm_mat_work_im)
2049 nrow=nao, ncol=nao, &
2050 s_firstrow=1, s_firstcol=1, &
2051 t_firstrow=nstart_row, t_firstcol=nstart_col)
2054 nrow=nao, ncol=nao, &
2055 s_firstrow=1, s_firstcol=1, &
2056 t_firstrow=nstart_row, t_firstcol=nstart_col)
2066 CALL timestop(handle)
2068 END SUBROUTINE add_cfm_submatrix
2075 SUBROUTINE create_cfm_double_row_col_size(fm_orig, cfm_double)
2079 CHARACTER(LEN=*),
PARAMETER :: routinen =
'create_cfm_double_row_col_size'
2081 INTEGER :: handle, ncol_global_orig, &
2085 CALL timeset(routinen, handle)
2087 CALL cp_fm_get_info(matrix=fm_orig, nrow_global=nrow_global_orig, ncol_global=ncol_global_orig)
2090 nrow_global=2*nrow_global_orig, &
2091 ncol_global=2*ncol_global_orig, &
2092 template_fmstruct=fm_orig%matrix_struct)
2098 CALL timestop(handle)
2100 END SUBROUTINE create_cfm_double_row_col_size
2115 SUBROUTINE print_gaps(E_VBM_SCF, E_CBM_SCF, E_VBM_SCF_beta, E_CBM_SCF_beta, &
2116 E_VBM_GW, E_CBM_GW, E_VBM_GW_beta, E_CBM_GW_beta, my_open_shell, unit_nr)
2118 REAL(kind=
dp) :: e_vbm_scf, e_cbm_scf, e_vbm_scf_beta, &
2119 e_cbm_scf_beta, e_vbm_gw, e_cbm_gw, &
2120 e_vbm_gw_beta, e_cbm_gw_beta
2121 LOGICAL :: my_open_shell
2124 IF (my_open_shell)
THEN
2125 WRITE (unit_nr,
'(T3,A)')
' '
2126 WRITE (unit_nr,
'(T3,A,F43.4)')
'Alpha SCF valence band maximum (eV)', e_vbm_scf*
evolt
2127 WRITE (unit_nr,
'(T3,A,F40.4)')
'Alpha SCF conduction band minimum (eV)', e_cbm_scf*
evolt
2128 WRITE (unit_nr,
'(T3,A,F56.4)')
'Alpha SCF bandgap (eV)', (e_cbm_scf - e_vbm_scf)*
evolt
2129 WRITE (unit_nr,
'(T3,A)')
' '
2130 WRITE (unit_nr,
'(T3,A,F44.4)')
'Beta SCF valence band maximum (eV)', e_vbm_scf_beta*
evolt
2131 WRITE (unit_nr,
'(T3,A,F41.4)')
'Beta SCF conduction band minimum (eV)', e_cbm_scf_beta*
evolt
2132 WRITE (unit_nr,
'(T3,A,F57.4)')
'Beta SCF bandgap (eV)', (e_cbm_scf_beta - e_vbm_scf_beta)*
evolt
2133 WRITE (unit_nr,
'(T3,A)')
' '
2134 WRITE (unit_nr,
'(T3,A,F44.4)')
'Alpha GW valence band maximum (eV)', e_vbm_gw*
evolt
2135 WRITE (unit_nr,
'(T3,A,F41.4)')
'Alpha GW conduction band minimum (eV)', e_cbm_gw*
evolt
2136 WRITE (unit_nr,
'(T3,A,F57.4)')
'Alpha GW bandgap (eV)', (e_cbm_gw - e_vbm_gw)*
evolt
2137 WRITE (unit_nr,
'(T3,A)')
' '
2138 WRITE (unit_nr,
'(T3,A,F45.4)')
'Beta GW valence band maximum (eV)', e_vbm_gw_beta*
evolt
2139 WRITE (unit_nr,
'(T3,A,F42.4)')
'Beta GW conduction band minimum (eV)', e_cbm_gw_beta*
evolt
2140 WRITE (unit_nr,
'(T3,A,F58.4)')
'Beta GW bandgap (eV)', (e_cbm_gw_beta - e_vbm_gw_beta)*
evolt
2142 WRITE (unit_nr,
'(T3,A)')
' '
2143 WRITE (unit_nr,
'(T3,A,F49.4)')
'SCF valence band maximum (eV)', e_vbm_scf*
evolt
2144 WRITE (unit_nr,
'(T3,A,F46.4)')
'SCF conduction band minimum (eV)', e_cbm_scf*
evolt
2145 WRITE (unit_nr,
'(T3,A,F62.4)')
'SCF bandgap (eV)', (e_cbm_scf - e_vbm_scf)*
evolt
2146 WRITE (unit_nr,
'(T3,A)')
' '
2147 WRITE (unit_nr,
'(T3,A,F50.4)')
'GW valence band maximum (eV)', e_vbm_gw*
evolt
2148 WRITE (unit_nr,
'(T3,A,F47.4)')
'GW conduction band minimum (eV)', e_cbm_gw*
evolt
2149 WRITE (unit_nr,
'(T3,A,F63.4)')
'GW bandgap (eV)', (e_cbm_gw - e_vbm_gw)*
evolt
2152 END SUBROUTINE print_gaps
2159 SUBROUTINE check_nan(array, real_value)
2160 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
2161 INTENT(INOUT) :: array
2162 REAL(kind=
dp),
INTENT(IN) :: real_value
2164 CHARACTER(LEN=*),
PARAMETER :: routinen =
'check_NaN'
2166 INTEGER :: handle, i, j, k
2168 CALL timeset(routinen, handle)
2170 DO i = 1,
SIZE(array, 1)
2171 DO j = 1,
SIZE(array, 2)
2172 DO k = 1,
SIZE(array, 3)
2175 IF (array(i, j, k) /= array(i, j, k)) array(i, j, k) = real_value
2181 CALL timestop(handle)
2183 END SUBROUTINE check_nan
2194 SUBROUTINE print_local_bandgap(qs_env, Eigenval, gw_corr_lev_occ, gw_corr_lev_virt, homo, dft_gw_char)
2196 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(IN) :: eigenval
2197 INTEGER :: gw_corr_lev_occ, gw_corr_lev_virt, homo
2198 CHARACTER(len=*) :: dft_gw_char
2200 CHARACTER(LEN=*),
PARAMETER :: routinen =
'print_local_bandgap'
2202 INTEGER :: handle, i_e
2208 CALL timeset(routinen, handle)
2210 CALL create_real_space_grids(e_gap_rspace, e_vbm_rspace, e_cbm_rspace, rho_g_dummy, ldos, auxbas_pw_pool, qs_env)
2212 CALL calculate_e_gap_rspace(e_gap_rspace, e_vbm_rspace, e_cbm_rspace, rho_g_dummy, &
2213 ldos, qs_env, eigenval, gw_corr_lev_occ, gw_corr_lev_virt, homo, dft_gw_char)
2215 CALL auxbas_pw_pool%give_back_pw(e_gap_rspace)
2216 CALL auxbas_pw_pool%give_back_pw(e_vbm_rspace)
2217 CALL auxbas_pw_pool%give_back_pw(e_cbm_rspace)
2218 CALL auxbas_pw_pool%give_back_pw(rho_g_dummy)
2219 DO i_e = 1,
SIZE(ldos)
2220 CALL auxbas_pw_pool%give_back_pw(ldos(i_e))
2224 CALL timestop(handle)
2226 END SUBROUTINE print_local_bandgap
2242 SUBROUTINE calculate_e_gap_rspace(E_gap_rspace, E_VBM_rspace, E_CBM_rspace, rho_g_dummy, &
2243 LDOS, qs_env, Eigenval, gw_corr_lev_occ, gw_corr_lev_virt, homo, dft_gw_char)
2248 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(IN) :: eigenval
2249 INTEGER :: gw_corr_lev_occ, gw_corr_lev_virt, homo
2250 CHARACTER(len=*) :: dft_gw_char
2252 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_E_gap_rspace'
2254 INTEGER :: handle, i_e, i_img, i_spin, i_x, i_y, i_z, ikp, imo, n_e, n_e_occ, n_x_end, &
2255 n_x_start, n_y_end, n_y_start, n_z_end, n_z_start, nimg, nkp, nkp_self_energy
2256 REAL(kind=
dp) :: avg_ldos_occ, avg_ldos_virt, d_e, e_cbm, &
2257 e_cbm_at_k, e_diff, e_vbm, e_vbm_at_k
2258 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: e_array
2259 REAL(kind=
dp),
DIMENSION(:),
POINTER :: occupation
2261 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: fm_work
2262 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, rho_ao
2263 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rho_ao_weighted
2276 CALL timeset(routinen, handle)
2278 CALL get_qs_env(qs_env=qs_env, para_env=para_env, mp2_env=mp2_env, ks_env=ks_env, matrix_s=matrix_s, &
2279 scf_env=scf_env, sab_orb=sab_orb, dft_control=dft_control, subsys=subsys)
2282 nkp =
SIZE(eigenval, 2)
2288 e_vbm_at_k = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo, ikp, 1))
2289 IF (e_vbm_at_k > e_vbm) e_vbm = e_vbm_at_k
2291 e_cbm_at_k = minval(eigenval(homo + 1:homo + gw_corr_lev_virt, ikp, 1))
2292 IF (e_cbm_at_k < e_cbm) e_cbm = e_cbm_at_k
2296 d_e = mp2_env%ri_g0w0%energy_spacing_print_loc_bandgap
2298 n_e = int(mp2_env%ri_g0w0%energy_window_print_loc_bandgap/d_e)
2301 ALLOCATE (e_array(n_e))
2303 e_array(i_e) = e_vbm - real(n_e_occ - i_e, kind=
dp)*d_e
2305 DO i_e = n_e_occ + 1, n_e
2306 e_array(i_e) = e_cbm + real(i_e - n_e_occ - 1, kind=
dp)*d_e
2309 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
2311 nkp_self_energy = kpoints_sigma%nkp
2312 cpassert(nkp == nkp_self_energy)
2314 kpoints_sigma%sab_nl => sab_orb
2316 DEALLOCATE (kpoints_sigma%cell_to_index)
2317 NULLIFY (kpoints_sigma%cell_to_index)
2320 nimg = maxval(kpoints_sigma%cell_to_index)
2322 NULLIFY (rho_ao_weighted)
2327 ALLOCATE (rho_ao_weighted(i_spin, i_img)%matrix)
2328 CALL dbcsr_create(matrix=rho_ao_weighted(i_spin, i_img)%matrix, template=matrix_s(1)%matrix)
2330 CALL dbcsr_set(rho_ao_weighted(i_spin, i_img)%matrix, 0.0_dp)
2334 ALLOCATE (fm_work(nimg))
2335 matrix_struct => kpoints_sigma%kp_env(1)%kpoint_env%mos(1, 1)%mo_coeff%matrix_struct
2344 CALL get_mo_set(kpoints_sigma%kp_env(ikp)%kpoint_env%mos(1, 1), &
2345 occupation_numbers=occupation)
2347 occupation(:) = 0.0_dp
2348 DO imo = homo - gw_corr_lev_occ + 1, homo + gw_corr_lev_virt
2349 e_diff = e_array(i_e) - eigenval(imo, ikp, 1)
2350 occupation(imo) = exp(-(e_diff/d_e)**2)
2355 CALL get_mo_set(kpoints_sigma%kp_env(1)%kpoint_env%mos(1, 1), &
2356 occupation_numbers=occupation)
2363 matrix_s(1)%matrix, sab_orb, fm_work)
2365 rho_ao => rho_ao_weighted(1, :)
2369 rho_gspace=rho_g_dummy, &
2374 CALL dbcsr_set(rho_ao_weighted(i_spin, i_img)%matrix, 0.0_dp)
2380 n_x_start = lbound(ldos(1)%array, 1)
2381 n_x_end = ubound(ldos(1)%array, 1)
2382 n_y_start = lbound(ldos(1)%array, 2)
2383 n_y_end = ubound(ldos(1)%array, 2)
2384 n_z_start = lbound(ldos(1)%array, 3)
2385 n_z_end = ubound(ldos(1)%array, 3)
2390 DO i_x = n_x_start, n_x_end
2391 DO i_y = n_y_start, n_y_end
2392 DO i_z = n_z_start, n_z_end
2394 avg_ldos_occ = 0.0_dp
2396 avg_ldos_occ = avg_ldos_occ + ldos(i_e)%array(i_x, i_y, i_z)
2398 avg_ldos_occ = avg_ldos_occ/real(n_e_occ, kind=
dp)
2400 avg_ldos_virt = 0.0_dp
2401 DO i_e = n_e_occ + 1, n_e
2402 avg_ldos_virt = avg_ldos_virt + ldos(i_e)%array(i_x, i_y, i_z)
2404 avg_ldos_virt = avg_ldos_virt/real(n_e - n_e_occ, kind=
dp)
2407 DO i_e = n_e_occ, 1, -1
2408 IF (ldos(i_e)%array(i_x, i_y, i_z) > mp2_env%ri_g0w0%ldos_thresh_print_loc_bandgap*avg_ldos_occ)
THEN
2409 e_vbm_rspace%array(i_x, i_y, i_z) = e_array(i_e)
2415 DO i_e = n_e_occ + 1, n_e
2416 IF (ldos(i_e)%array(i_x, i_y, i_z) > mp2_env%ri_g0w0%ldos_thresh_print_loc_bandgap*avg_ldos_virt)
THEN
2417 e_cbm_rspace%array(i_x, i_y, i_z) = e_array(i_e)
2429 CALL pw_copy(e_cbm_rspace, e_gap_rspace)
2430 CALL pw_axpy(e_vbm_rspace, e_gap_rspace, -1.0_dp)
2435 CALL print_file(e_gap_rspace, dft_gw_char//
"_Gap_in_eV", gw_section, particles, mp2_env)
2436 CALL print_file(e_vbm_rspace, dft_gw_char//
"_VBM_in_eV", gw_section, particles, mp2_env)
2437 CALL print_file(e_cbm_rspace, dft_gw_char//
"_CBM_in_eV", gw_section, particles, mp2_env)
2438 CALL print_file(ldos(n_e_occ), dft_gw_char//
"_LDOS_VBM_in_eV", gw_section, particles, mp2_env)
2439 CALL print_file(ldos(n_e_occ + 1), dft_gw_char//
"_LDOS_CBM_in_eV", gw_section, particles, mp2_env)
2445 DEALLOCATE (e_array)
2447 NULLIFY (kpoints_sigma%sab_nl)
2449 CALL timestop(handle)
2451 END SUBROUTINE calculate_e_gap_rspace
2461 SUBROUTINE print_file(pw_print, middle_name, gw_section, particles, mp2_env)
2463 CHARACTER(len=*) :: middle_name
2468 CHARACTER(LEN=*),
PARAMETER :: routinen =
'print_file'
2470 INTEGER :: handle, unit_nr_cube
2474 CALL timeset(routinen, handle)
2479 unit_nr_cube =
cp_print_key_unit_nr(logger, gw_section,
"PRINT%LOCAL_BANDGAP", extension=
".cube", &
2480 middle_name=middle_name, file_form=
"FORMATTED", mpi_io=mpi_io)
2481 CALL cp_pw_to_cube(pw_print, unit_nr_cube, middle_name, particles=particles, &
2482 stride=mp2_env%ri_g0w0%stride_loc_bandgap, mpi_io=mpi_io)
2484 "PRINT%LOCAL_BANDGAP", mpi_io=mpi_io)
2486 CALL timestop(handle)
2488 END SUBROUTINE print_file
2500 SUBROUTINE create_real_space_grids(E_gap_rspace, E_VBM_rspace, E_CBM_rspace, rho_g_dummy, LDOS, auxbas_pw_pool, qs_env)
2507 CHARACTER(LEN=*),
PARAMETER :: routinen =
'create_real_space_grids'
2509 INTEGER :: handle, i_e, n_e
2513 CALL timeset(routinen, handle)
2515 CALL get_qs_env(qs_env=qs_env, mp2_env=mp2_env, pw_env=pw_env)
2517 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2519 CALL auxbas_pw_pool%create_pw(e_gap_rspace)
2520 CALL auxbas_pw_pool%create_pw(e_vbm_rspace)
2521 CALL auxbas_pw_pool%create_pw(e_cbm_rspace)
2522 CALL auxbas_pw_pool%create_pw(rho_g_dummy)
2524 n_e = int(mp2_env%ri_g0w0%energy_window_print_loc_bandgap/ &
2525 mp2_env%ri_g0w0%energy_spacing_print_loc_bandgap)
2527 ALLOCATE (ldos(n_e))
2530 CALL auxbas_pw_pool%create_pw(ldos(i_e))
2533 CALL timestop(handle)
2535 END SUBROUTINE create_real_space_grids
2562 SUBROUTINE calc_periodic_correction(delta_corr, qs_env, para_env, para_env_RPA, kp_grid, homo, nmo, &
2563 gw_corr_lev_occ, gw_corr_lev_virt, omega, fm_mo_coeff, Eigenval, &
2564 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
2565 first_cycle_periodic_correction, kpoints, do_mo_coeff_Gamma_only, &
2566 num_kp_grids, eps_kpoint, do_extra_kpoints, do_aux_bas, frac_aux_mos)
2568 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
2569 INTENT(INOUT) :: delta_corr
2572 INTEGER,
DIMENSION(:),
POINTER :: kp_grid
2573 INTEGER,
INTENT(IN) :: homo, nmo, gw_corr_lev_occ, &
2575 REAL(kind=
dp),
INTENT(IN) :: omega
2577 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eigenval
2578 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_re_mo_mo, &
2579 matrix_berry_im_mo_mo
2580 LOGICAL,
INTENT(INOUT) :: first_cycle_periodic_correction
2582 LOGICAL,
INTENT(IN) :: do_mo_coeff_gamma_only
2583 INTEGER,
INTENT(IN) :: num_kp_grids
2584 REAL(kind=
dp),
INTENT(IN) :: eps_kpoint
2585 LOGICAL,
INTENT(IN) :: do_extra_kpoints, do_aux_bas
2586 REAL(kind=
dp),
INTENT(IN) :: frac_aux_mos
2588 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_periodic_correction'
2591 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eps_head, eps_inv_head
2592 REAL(kind=
dp),
DIMENSION(3, 3) :: h_inv
2594 CALL timeset(routinen, handle)
2596 IF (first_cycle_periodic_correction)
THEN
2598 CALL get_kpoints(qs_env, kpoints, kp_grid, num_kp_grids, para_env, h_inv, nmo, do_mo_coeff_gamma_only, &
2601 CALL get_berry_phase(qs_env, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, fm_mo_coeff, &
2602 para_env, do_mo_coeff_gamma_only, homo, nmo, gw_corr_lev_virt, eps_kpoint, do_aux_bas, &
2607 CALL compute_eps_head_berry(eps_head, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, para_env_rpa, &
2608 qs_env, homo, eigenval, omega)
2610 CALL compute_eps_inv_head(eps_inv_head, eps_head, kpoints)
2612 CALL kpoint_sum_for_eps_inv_head_berry(delta_corr, eps_inv_head, kpoints, qs_env, &
2613 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
2614 homo, gw_corr_lev_occ, gw_corr_lev_virt, para_env_rpa, &
2617 DEALLOCATE (eps_head, eps_inv_head)
2619 first_cycle_periodic_correction = .false.
2621 CALL timestop(handle)
2623 END SUBROUTINE calc_periodic_correction
2637 SUBROUTINE compute_eps_head_berry(eps_head, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, para_env_RPA, &
2638 qs_env, homo, Eigenval, omega)
2640 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
2641 INTENT(OUT) :: eps_head
2643 TYPE(
dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_berry_re_mo_mo, &
2644 matrix_berry_im_mo_mo
2647 INTEGER,
INTENT(IN) :: homo
2648 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eigenval
2649 REAL(kind=
dp),
INTENT(IN) :: omega
2651 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_eps_head_Berry'
2653 INTEGER :: col, col_end_in_block, col_offset, col_size, handle, i_col, i_row, ikp, nkp, nmo, &
2654 row, row_offset, row_size, row_start_in_block
2655 REAL(kind=
dp) :: abs_k_square, cell_volume, &
2656 correct_kpoint(3), cos_square, &
2657 eigen_diff, relative_kpoint(3), &
2659 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: p_head
2660 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_block
2664 CALL timeset(routinen, handle)
2667 CALL get_cell(cell=cell, deth=cell_volume)
2669 NULLIFY (data_block)
2673 nmo =
SIZE(eigenval)
2675 ALLOCATE (p_head(nkp))
2678 ALLOCATE (eps_head(nkp))
2679 eps_head(:) = 0.0_dp
2683 relative_kpoint(1:3) = matmul(cell%hmat, kpoints%xkp(1:3, ikp))
2685 correct_kpoint(1:3) =
twopi*kpoints%xkp(1:3, ikp)
2687 abs_k_square = (correct_kpoint(1))**2 + (correct_kpoint(2))**2 + (correct_kpoint(3))**2
2694 row_size=row_size, col_size=col_size, &
2695 row_offset=row_offset, col_offset=col_offset)
2697 IF (row_offset + row_size <= homo .OR. col_offset > homo) cycle
2699 IF (row_offset <= homo)
THEN
2700 row_start_in_block = homo - row_offset + 2
2702 row_start_in_block = 1
2705 IF (col_offset + col_size - 1 > homo)
THEN
2706 col_end_in_block = homo - col_offset + 1
2708 col_end_in_block = col_size
2711 DO i_row = row_start_in_block, min(row_size, nmo - row_offset + 1)
2713 DO i_col = 1, min(col_end_in_block, nmo - col_offset + 1)
2715 eigen_diff = eigenval(i_col + col_offset - 1) - eigenval(i_row + row_offset - 1)
2717 cos_square = (data_block(i_row, i_col))**2
2719 p_head(ikp) = p_head(ikp) + 2.0_dp*eigen_diff/(omega**2 + eigen_diff**2)*cos_square/abs_k_square
2734 row_size=row_size, col_size=col_size, &
2735 row_offset=row_offset, col_offset=col_offset)
2737 IF (row_offset + row_size <= homo .OR. col_offset > homo) cycle
2739 IF (row_offset <= homo)
THEN
2740 row_start_in_block = homo - row_offset + 2
2742 row_start_in_block = 1
2745 IF (col_offset + col_size - 1 > homo)
THEN
2746 col_end_in_block = homo - col_offset + 1
2748 col_end_in_block = col_size
2751 DO i_row = row_start_in_block, min(row_size, nmo - row_offset + 1)
2753 DO i_col = 1, min(col_end_in_block, nmo - col_offset + 1)
2755 eigen_diff = eigenval(i_col + col_offset - 1) - eigenval(i_row + row_offset - 1)
2757 sin_square = (data_block(i_row, i_col))**2
2759 p_head(ikp) = p_head(ikp) + 2.0_dp*eigen_diff/(omega**2 + eigen_diff**2)*sin_square/abs_k_square
2771 CALL para_env_rpa%sum(p_head)
2775 eps_head(:) = 1.0_dp - 2.0_dp*p_head(:)/cell_volume*
fourpi
2779 CALL timestop(handle)
2781 END SUBROUTINE compute_eps_head_berry
2799 SUBROUTINE get_berry_phase(qs_env, kpoints, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, fm_mo_coeff, para_env, &
2800 do_mo_coeff_Gamma_only, homo, nmo, gw_corr_lev_virt, eps_kpoint, do_aux_bas, &
2804 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_re_mo_mo, &
2805 matrix_berry_im_mo_mo
2808 LOGICAL,
INTENT(IN) :: do_mo_coeff_gamma_only
2809 INTEGER,
INTENT(IN) :: homo, nmo, gw_corr_lev_virt
2810 REAL(kind=
dp),
INTENT(IN) :: eps_kpoint
2811 LOGICAL,
INTENT(IN) :: do_aux_bas
2812 REAL(kind=
dp),
INTENT(IN) :: frac_aux_mos
2814 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_berry_phase'
2816 INTEGER :: col_index, handle, i_col_local, ikind, &
2817 ikp, nao_aux, ncol_local, nkind, nkp, &
2819 INTEGER,
DIMENSION(:),
POINTER :: col_indices
2820 REAL(
dp) :: abs_kpoint, correct_kpoint(3), &
2822 REAL(kind=
dp),
DIMENSION(:),
POINTER :: evals_p, evals_p_sqrt_inv
2825 TYPE(
cp_fm_type) :: fm_mat_eigv_p, fm_mat_p, fm_mat_p_sqrt_inv, fm_mat_s_aux_aux_inv, &
2826 fm_mat_scaled_eigv_p, fm_mat_work_aux_aux
2827 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, matrix_s_aux_aux, &
2829 TYPE(
dbcsr_type),
POINTER :: cosmat, cosmat_desymm, mat_mo_coeff_aux, mat_mo_coeff_aux_2, &
2830 mat_mo_coeff_gamma_all, mat_mo_coeff_gamma_occ_and_gw, mat_mo_coeff_im, mat_mo_coeff_re, &
2831 mat_work_aux_orb, mat_work_aux_orb_2, matrix_p, matrix_p_sqrt, matrix_p_sqrt_inv, &
2832 matrix_s_inv_aux_aux, sinmat, sinmat_desymm, tmp
2836 POINTER :: sab_orb, sab_orb_mic, sgwgw_list, &
2838 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2842 CALL timeset(routinen, handle)
2846 NULLIFY (matrix_berry_re_mo_mo, matrix_s, cell, matrix_berry_im_mo_mo, sinmat, cosmat, tmp, &
2847 cosmat_desymm, sinmat_desymm, qs_kind_set, orb_basis_set_list, sab_orb_mic)
2851 matrix_s=matrix_s, &
2852 qs_kind_set=qs_kind_set, &
2857 ALLOCATE (orb_basis_set_list(nkind))
2863 NULLIFY (mat_mo_coeff_re)
2866 template=matrix_s(1)%matrix, &
2867 matrix_type=dbcsr_type_no_symmetry)
2869 NULLIFY (mat_mo_coeff_im)
2872 template=matrix_s(1)%matrix, &
2873 matrix_type=dbcsr_type_no_symmetry)
2875 NULLIFY (mat_mo_coeff_gamma_all)
2878 template=matrix_s(1)%matrix, &
2879 matrix_type=dbcsr_type_no_symmetry)
2881 CALL copy_fm_to_dbcsr(fm_mo_coeff, mat_mo_coeff_gamma_all, keep_sparsity=.false.)
2883 NULLIFY (mat_mo_coeff_gamma_occ_and_gw)
2885 CALL dbcsr_create(matrix=mat_mo_coeff_gamma_occ_and_gw, &
2886 template=matrix_s(1)%matrix, &
2887 matrix_type=dbcsr_type_no_symmetry)
2889 CALL copy_fm_to_dbcsr(fm_mo_coeff, mat_mo_coeff_gamma_occ_and_gw, keep_sparsity=.false.)
2891 IF (.NOT. do_aux_bas)
THEN
2899 CALL dbcsr_create(matrix=cosmat, template=matrix_s(1)%matrix)
2900 CALL dbcsr_create(matrix=sinmat, template=matrix_s(1)%matrix)
2902 template=matrix_s(1)%matrix, &
2903 matrix_type=dbcsr_type_no_symmetry)
2905 template=matrix_s(1)%matrix, &
2906 matrix_type=dbcsr_type_no_symmetry)
2908 template=matrix_s(1)%matrix, &
2909 matrix_type=dbcsr_type_no_symmetry)
2918 NULLIFY (gw_aux_basis_set_list)
2919 ALLOCATE (gw_aux_basis_set_list(nkind))
2923 NULLIFY (gw_aux_basis_set_list(ikind)%gto_basis_set)
2925 NULLIFY (basis_set_gw_aux)
2927 qs_kind => qs_kind_set(ikind)
2928 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_gw_aux, basis_type=
"AUX_GW")
2929 cpassert(
ASSOCIATED(basis_set_gw_aux))
2931 basis_set_gw_aux%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius
2933 gw_aux_basis_set_list(ikind)%gto_basis_set => basis_set_gw_aux
2938 NULLIFY (sgwgw_list, sgworb_list)
2940 CALL setup_neighbor_list(sgworb_list, gw_aux_basis_set_list, orb_basis_set_list, qs_env=qs_env)
2942 NULLIFY (matrix_s_aux_aux, matrix_s_aux_orb)
2946 gw_aux_basis_set_list, gw_aux_basis_set_list, sgwgw_list)
2949 gw_aux_basis_set_list, orb_basis_set_list, sgworb_list)
2951 CALL dbcsr_get_info(matrix_s_aux_aux(1)%matrix, nfullrows_total=nao_aux)
2953 nmo_for_aux_bas = floor(frac_aux_mos*real(nao_aux, kind=
dp))
2956 context=fm_mo_coeff%matrix_struct%context, &
2957 nrow_global=nao_aux, &
2958 ncol_global=nao_aux, &
2961 NULLIFY (mat_work_aux_orb)
2964 template=matrix_s_aux_orb(1)%matrix, &
2965 matrix_type=dbcsr_type_no_symmetry)
2967 NULLIFY (mat_work_aux_orb_2)
2970 template=matrix_s_aux_orb(1)%matrix, &
2971 matrix_type=dbcsr_type_no_symmetry)
2973 NULLIFY (mat_mo_coeff_aux)
2976 template=matrix_s_aux_orb(1)%matrix, &
2977 matrix_type=dbcsr_type_no_symmetry)
2979 NULLIFY (mat_mo_coeff_aux_2)
2982 template=matrix_s_aux_orb(1)%matrix, &
2983 matrix_type=dbcsr_type_no_symmetry)
2985 NULLIFY (matrix_s_inv_aux_aux)
2988 template=matrix_s_aux_aux(1)%matrix, &
2989 matrix_type=dbcsr_type_no_symmetry)
2994 template=matrix_s(1)%matrix, &
2995 matrix_type=dbcsr_type_no_symmetry)
2997 NULLIFY (matrix_p_sqrt)
3000 template=matrix_s(1)%matrix, &
3001 matrix_type=dbcsr_type_no_symmetry)
3003 NULLIFY (matrix_p_sqrt_inv)
3006 template=matrix_s(1)%matrix, &
3007 matrix_type=dbcsr_type_no_symmetry)
3009 CALL cp_fm_create(fm_mat_s_aux_aux_inv, fm_struct_aux_aux, name=
"inverse overlap mat")
3010 CALL cp_fm_create(fm_mat_work_aux_aux, fm_struct_aux_aux, name=
"work mat")
3012 CALL cp_fm_create(fm_mat_eigv_p, fm_mo_coeff%matrix_struct)
3013 CALL cp_fm_create(fm_mat_scaled_eigv_p, fm_mo_coeff%matrix_struct)
3014 CALL cp_fm_create(fm_mat_p_sqrt_inv, fm_mo_coeff%matrix_struct)
3017 ALLOCATE (evals_p(nmo))
3019 NULLIFY (evals_p_sqrt_inv)
3020 ALLOCATE (evals_p_sqrt_inv(nmo))
3029 CALL copy_fm_to_dbcsr(fm_mat_s_aux_aux_inv, matrix_s_inv_aux_aux, keep_sparsity=.false.)
3031 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, &
3032 filter_eps=1.0e-15_dp)
3034 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, &
3035 last_column=nmo_for_aux_bas, filter_eps=1.0e-15_dp)
3037 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, &
3038 filter_eps=1.0e-15_dp)
3040 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_aux_2, mat_work_aux_orb, 0.0_dp, matrix_p, &
3041 filter_eps=1.0e-15_dp)
3045 CALL cp_fm_syevd(fm_mat_p, fm_mat_eigv_p, evals_p)
3048 evals_p_sqrt_inv(1:nmo - nmo_for_aux_bas) = 0.0_dp
3049 evals_p_sqrt_inv(nmo - nmo_for_aux_bas + 1:nmo) = 1.0_dp/sqrt(evals_p(nmo - nmo_for_aux_bas + 1:nmo))
3051 CALL cp_fm_to_fm(fm_mat_eigv_p, fm_mat_scaled_eigv_p)
3054 ncol_local=ncol_local, &
3055 col_indices=col_indices)
3057 CALL para_env%sync()
3060 DO i_col_local = 1, ncol_local
3062 col_index = col_indices(i_col_local)
3064 fm_mat_scaled_eigv_p%local_data(:, i_col_local) = &
3065 fm_mat_scaled_eigv_p%local_data(:, i_col_local)*evals_p_sqrt_inv(col_index)
3069 CALL para_env%sync()
3071 CALL parallel_gemm(transa=
"N", transb=
"T", m=nmo, n=nmo, k=nmo, alpha=1.0_dp, &
3072 matrix_a=fm_mat_eigv_p, matrix_b=fm_mat_scaled_eigv_p, beta=0.0_dp, &
3073 matrix_c=fm_mat_p_sqrt_inv)
3075 CALL copy_fm_to_dbcsr(fm_mat_p_sqrt_inv, matrix_p_sqrt_inv, keep_sparsity=.false.)
3077 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, mat_mo_coeff_aux_2, matrix_p_sqrt_inv, 0.0_dp, mat_mo_coeff_aux, &
3078 filter_eps=1.0e-15_dp)
3086 CALL dbcsr_create(matrix=cosmat, template=matrix_s_aux_aux(1)%matrix)
3087 CALL dbcsr_create(matrix=sinmat, template=matrix_s_aux_aux(1)%matrix)
3089 template=matrix_s_aux_orb(1)%matrix, &
3090 matrix_type=dbcsr_type_no_symmetry)
3092 template=matrix_s_aux_aux(1)%matrix, &
3093 matrix_type=dbcsr_type_no_symmetry)
3095 template=matrix_s_aux_aux(1)%matrix, &
3096 matrix_type=dbcsr_type_no_symmetry)
3097 CALL dbcsr_copy(cosmat, matrix_s_aux_aux(1)%matrix)
3098 CALL dbcsr_copy(sinmat, matrix_s_aux_aux(1)%matrix)
3107 NULLIFY (mat_mo_coeff_gamma_all)
3110 template=matrix_s_aux_orb(1)%matrix, &
3111 matrix_type=dbcsr_type_no_symmetry)
3113 CALL dbcsr_copy(mat_mo_coeff_gamma_all, mat_mo_coeff_aux)
3115 NULLIFY (mat_mo_coeff_gamma_occ_and_gw)
3117 CALL dbcsr_create(matrix=mat_mo_coeff_gamma_occ_and_gw, &
3118 template=matrix_s_aux_orb(1)%matrix, &
3119 matrix_type=dbcsr_type_no_symmetry)
3121 CALL dbcsr_copy(mat_mo_coeff_gamma_occ_and_gw, mat_mo_coeff_aux)
3123 DEALLOCATE (evals_p, evals_p_sqrt_inv)
3127 CALL remove_unnecessary_blocks(mat_mo_coeff_gamma_occ_and_gw, homo, gw_corr_lev_virt)
3131 ALLOCATE (matrix_berry_re_mo_mo(ikp)%matrix)
3134 template=matrix_s(1)%matrix, &
3135 matrix_type=dbcsr_type_no_symmetry)
3137 CALL dbcsr_set(matrix_berry_re_mo_mo(ikp)%matrix, 0.0_dp)
3139 ALLOCATE (matrix_berry_im_mo_mo(ikp)%matrix)
3142 template=matrix_s(1)%matrix, &
3143 matrix_type=dbcsr_type_no_symmetry)
3145 CALL dbcsr_set(matrix_berry_im_mo_mo(ikp)%matrix, 0.0_dp)
3147 correct_kpoint(1:3) = -
twopi*kpoints%xkp(1:3, ikp)
3149 abs_kpoint = sqrt(correct_kpoint(1)**2 + correct_kpoint(2)**2 + correct_kpoint(3)**2)
3151 IF (abs_kpoint < eps_kpoint)
THEN
3153 scale_kpoint = eps_kpoint/abs_kpoint
3154 correct_kpoint(:) = correct_kpoint(:)*scale_kpoint
3159 IF (do_aux_bas)
THEN
3161 basis_type=
"AUX_GW")
3167 IF (do_mo_coeff_gamma_only)
THEN
3171 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_gamma_occ_and_gw, 0.0_dp, tmp, &
3172 filter_eps=1.0e-15_dp)
3174 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 0.0_dp, &
3175 matrix_berry_re_mo_mo(ikp)%matrix, filter_eps=1.0e-15_dp)
3179 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_gamma_occ_and_gw, 0.0_dp, tmp, &
3180 filter_eps=1.0e-15_dp)
3182 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 0.0_dp, &
3183 matrix_berry_im_mo_mo(ikp)%matrix, filter_eps=1.0e-15_dp)
3189 mat_mo_coeff_re, keep_sparsity=.false.)
3192 mat_mo_coeff_im, keep_sparsity=.false.)
3199 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_re, 0.0_dp, tmp)
3202 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 0.0_dp, &
3203 matrix_berry_re_mo_mo(ikp)%matrix)
3206 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_re, 0.0_dp, tmp)
3209 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 0.0_dp, &
3210 matrix_berry_im_mo_mo(ikp)%matrix)
3213 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, cosmat_desymm, mat_mo_coeff_im, 0.0_dp, tmp)
3216 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, mat_mo_coeff_gamma_all, tmp, 1.0_dp, &
3217 matrix_berry_im_mo_mo(ikp)%matrix)
3220 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, sinmat_desymm, mat_mo_coeff_im, 0.0_dp, tmp)
3223 CALL dbcsr_multiply(
'T',
'N', -1.0_dp, mat_mo_coeff_gamma_all, tmp, 1.0_dp, &
3224 matrix_berry_re_mo_mo(ikp)%matrix)
3228 IF (abs_kpoint < eps_kpoint)
THEN
3230 CALL dbcsr_scale(matrix_berry_im_mo_mo(ikp)%matrix, 1.0_dp/scale_kpoint)
3231 CALL dbcsr_set(matrix_berry_re_mo_mo(ikp)%matrix, 0.0_dp)
3247 DEALLOCATE (orb_basis_set_list)
3251 IF (do_aux_bas)
THEN
3253 DEALLOCATE (gw_aux_basis_set_list)
3280 CALL timestop(handle)
3282 END SUBROUTINE get_berry_phase
3290 SUBROUTINE remove_unnecessary_blocks(mat_mo_coeff_Gamma_occ_and_GW, homo, gw_corr_lev_virt)
3292 TYPE(
dbcsr_type),
POINTER :: mat_mo_coeff_gamma_occ_and_gw
3293 INTEGER,
INTENT(IN) :: homo, gw_corr_lev_virt
3295 INTEGER :: col, col_offset, row
3296 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_block
3304 col_offset=col_offset)
3306 IF (col_offset > homo + gw_corr_lev_virt)
THEN
3316 CALL dbcsr_filter(mat_mo_coeff_gamma_occ_and_gw, 1.0e-15_dp)
3318 END SUBROUTINE remove_unnecessary_blocks
3334 SUBROUTINE kpoint_sum_for_eps_inv_head_berry(delta_corr, eps_inv_head, kpoints, qs_env, matrix_berry_re_mo_mo, &
3335 matrix_berry_im_mo_mo, homo, gw_corr_lev_occ, gw_corr_lev_virt, &
3336 para_env_RPA, do_extra_kpoints)
3338 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
3339 INTENT(INOUT) :: delta_corr
3340 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eps_inv_head
3343 TYPE(
dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_berry_re_mo_mo, &
3344 matrix_berry_im_mo_mo
3345 INTEGER,
INTENT(IN) :: homo, gw_corr_lev_occ, gw_corr_lev_virt
3347 LOGICAL,
INTENT(IN) :: do_extra_kpoints
3349 INTEGER :: col, col_offset, col_size, i_col, i_row, &
3350 ikp, m_level, n_level_gw, nkp, row, &
3351 row_offset, row_size
3352 REAL(kind=
dp) :: abs_k_square, cell_volume, &
3353 check_int_one_over_ksq, contribution, &
3355 REAL(kind=
dp),
DIMENSION(3) :: correct_kpoint
3356 REAL(kind=
dp),
DIMENSION(:),
POINTER :: delta_corr_extra
3357 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_block
3363 CALL get_cell(cell=cell, deth=cell_volume)
3369 IF (do_extra_kpoints)
THEN
3370 NULLIFY (delta_corr_extra)
3371 ALLOCATE (delta_corr_extra(1 + homo - gw_corr_lev_occ:homo + gw_corr_lev_virt))
3372 delta_corr_extra = 0.0_dp
3375 check_int_one_over_ksq = 0.0_dp
3379 weight = kpoints%wkp(ikp)
3381 correct_kpoint(1:3) =
twopi*kpoints%xkp(1:3, ikp)
3383 abs_k_square = (correct_kpoint(1))**2 + (correct_kpoint(2))**2 + (correct_kpoint(3))**2
3390 row_size=row_size, col_size=col_size, &
3391 row_offset=row_offset, col_offset=col_offset)
3393 DO i_col = 1, col_size
3395 DO n_level_gw = 1 + homo - gw_corr_lev_occ, homo + gw_corr_lev_virt
3397 IF (n_level_gw == i_col + col_offset - 1)
THEN
3399 DO i_row = 1, row_size
3401 contribution = weight*(eps_inv_head(ikp) - 1.0_dp)/abs_k_square*(data_block(i_row, i_col))**2
3403 m_level = i_row + row_offset - 1
3406 IF (m_level /= n_level_gw) cycle
3408 IF (.NOT. do_extra_kpoints)
THEN
3410 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution
3414 IF (ikp <= nkp*8/9)
THEN
3416 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution
3420 delta_corr_extra(n_level_gw) = delta_corr_extra(n_level_gw) + contribution
3443 row_size=row_size, col_size=col_size, &
3444 row_offset=row_offset, col_offset=col_offset)
3446 DO i_col = 1, col_size
3448 DO n_level_gw = 1 + homo - gw_corr_lev_occ, homo + gw_corr_lev_virt
3450 IF (n_level_gw == i_col + col_offset - 1)
THEN
3452 DO i_row = 1, row_size
3454 m_level = i_row + row_offset - 1
3456 contribution = weight*(eps_inv_head(ikp) - 1.0_dp)/abs_k_square*(data_block(i_row, i_col))**2
3459 IF (m_level /= n_level_gw) cycle
3461 IF (.NOT. do_extra_kpoints)
THEN
3463 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution
3467 IF (ikp <= nkp*8/9)
THEN
3469 delta_corr(n_level_gw) = delta_corr(n_level_gw) + contribution
3473 delta_corr_extra(n_level_gw) = delta_corr_extra(n_level_gw) + contribution
3491 check_int_one_over_ksq = check_int_one_over_ksq + weight/abs_k_square
3496 delta_corr = delta_corr/cell_volume*
fourpi
3498 check_int_one_over_ksq = check_int_one_over_ksq/cell_volume
3500 CALL para_env_rpa%sum(delta_corr)
3502 IF (do_extra_kpoints)
THEN
3504 delta_corr_extra = delta_corr_extra/cell_volume*
fourpi
3506 CALL para_env_rpa%sum(delta_corr_extra)
3508 delta_corr(:) = delta_corr(:) + (delta_corr(:) - delta_corr_extra(:))
3510 DEALLOCATE (delta_corr_extra)
3514 END SUBROUTINE kpoint_sum_for_eps_inv_head_berry
3522 SUBROUTINE compute_eps_inv_head(eps_inv_head, eps_head, kpoints)
3523 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
3524 INTENT(OUT) :: eps_inv_head
3525 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eps_head
3528 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_eps_inv_head'
3530 INTEGER :: handle, ikp, nkp
3532 CALL timeset(routinen, handle)
3536 ALLOCATE (eps_inv_head(nkp))
3540 eps_inv_head(ikp) = 1.0_dp/eps_head(ikp)
3544 CALL timestop(handle)
3546 END SUBROUTINE compute_eps_inv_head
3560 SUBROUTINE get_kpoints(qs_env, kpoints, kp_grid, num_kp_grids, para_env, h_inv, nmo, &
3561 do_mo_coeff_Gamma_only, do_extra_kpoints)
3564 INTEGER,
DIMENSION(:),
POINTER :: kp_grid
3565 INTEGER,
INTENT(IN) :: num_kp_grids
3567 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT) :: h_inv
3568 INTEGER,
INTENT(IN) :: nmo
3569 LOGICAL,
INTENT(IN) :: do_mo_coeff_gamma_only, do_extra_kpoints
3571 INTEGER :: end_kp, i, i_grid_level, ix, iy, iz, &
3572 nkp_inner_grid, nkp_outer_grid, &
3574 INTEGER,
DIMENSION(3) :: outer_kp_grid
3575 REAL(kind=
dp) :: kpoint_weight_left, single_weight
3576 REAL(kind=
dp),
DIMENSION(3) :: kpt_latt, reducing_factor
3580 NULLIFY (kpoints, cell, particle_set)
3583 cpassert(mod(kp_grid(1)*kp_grid(2)*kp_grid(3), 2) == 0)
3584 IF (do_extra_kpoints)
THEN
3585 cpassert(do_mo_coeff_gamma_only)
3588 IF (do_mo_coeff_gamma_only)
THEN
3590 outer_kp_grid(1) = kp_grid(1) - 1
3591 outer_kp_grid(2) = kp_grid(2) - 1
3592 outer_kp_grid(3) = kp_grid(3) - 1
3594 CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set)
3600 kpoints%kp_scheme =
"GENERAL"
3601 kpoints%symmetry = .false.
3602 kpoints%verbose = .false.
3603 kpoints%full_grid = .false.
3604 kpoints%use_real_wfn = .false.
3605 kpoints%eps_geo = 1.e-6_dp
3606 npoints = kp_grid(1)*kp_grid(2)*kp_grid(3)/2 + &
3607 (num_kp_grids - 1)*((outer_kp_grid(1) + 1)/2*outer_kp_grid(2)*outer_kp_grid(3) - 1)
3609 IF (do_extra_kpoints)
THEN
3611 cpassert(num_kp_grids == 1)
3612 cpassert(mod(kp_grid(1), 4) == 0)
3613 cpassert(mod(kp_grid(2), 4) == 0)
3614 cpassert(mod(kp_grid(3), 4) == 0)
3618 IF (do_extra_kpoints)
THEN
3620 npoints = kp_grid(1)*kp_grid(2)*kp_grid(3)/2 + kp_grid(1)*kp_grid(2)*kp_grid(3)/2/8
3624 kpoints%full_grid = .true.
3625 kpoints%nkp = npoints
3626 ALLOCATE (kpoints%xkp(3, npoints), kpoints%wkp(npoints))
3627 kpoints%xkp = 0.0_dp
3628 kpoints%wkp = 0.0_dp
3630 nkp_outer_grid = outer_kp_grid(1)*outer_kp_grid(2)*outer_kp_grid(3)
3631 nkp_inner_grid = kp_grid(1)*kp_grid(2)*kp_grid(3)
3634 reducing_factor(:) = 1.0_dp
3635 kpoint_weight_left = 1.0_dp
3638 DO i_grid_level = 1, num_kp_grids - 1
3640 single_weight = kpoint_weight_left/real(nkp_outer_grid, kind=
dp)
3644 DO ix = 1, outer_kp_grid(1)
3645 DO iy = 1, outer_kp_grid(2)
3646 DO iz = 1, outer_kp_grid(3)
3649 IF (2*ix - outer_kp_grid(1) - 1 == 0 .AND. 2*iy - outer_kp_grid(2) - 1 == 0 .AND. &
3650 2*iz - outer_kp_grid(3) - 1 == 0) cycle
3653 IF (2*ix - outer_kp_grid(1) - 1 < 0) cycle
3656 kpt_latt(1) = real(2*ix - outer_kp_grid(1) - 1, kind=
dp)/(2._dp*real(outer_kp_grid(1), kind=
dp)) &
3658 kpt_latt(2) = real(2*iy - outer_kp_grid(2) - 1, kind=
dp)/(2._dp*real(outer_kp_grid(2), kind=
dp)) &
3660 kpt_latt(3) = real(2*iz - outer_kp_grid(3) - 1, kind=
dp)/(2._dp*real(outer_kp_grid(3), kind=
dp)) &
3662 kpoints%xkp(1:3, i) = matmul(transpose(h_inv), kpt_latt(:))
3664 IF (2*ix - outer_kp_grid(1) - 1 == 0)
THEN
3665 kpoints%wkp(i) = single_weight
3667 kpoints%wkp(i) = 2._dp*single_weight
3676 kpoint_weight_left = kpoint_weight_left - sum(kpoints%wkp(start_kp:end_kp))
3678 reducing_factor(1) = reducing_factor(1)/real(outer_kp_grid(1), kind=
dp)
3679 reducing_factor(2) = reducing_factor(2)/real(outer_kp_grid(2), kind=
dp)
3680 reducing_factor(3) = reducing_factor(3)/real(outer_kp_grid(3), kind=
dp)
3684 single_weight = kpoint_weight_left/real(nkp_inner_grid, kind=
dp)
3687 DO ix = 1, kp_grid(1)
3688 DO iy = 1, kp_grid(2)
3689 DO iz = 1, kp_grid(3)
3692 IF (2*ix - kp_grid(1) - 1 < 0) cycle
3695 kpt_latt(1) = real(2*ix - kp_grid(1) - 1, kind=
dp)/(2._dp*real(kp_grid(1), kind=
dp))*reducing_factor(1)
3696 kpt_latt(2) = real(2*iy - kp_grid(2) - 1, kind=
dp)/(2._dp*real(kp_grid(2), kind=
dp))*reducing_factor(2)
3697 kpt_latt(3) = real(2*iz - kp_grid(3) - 1, kind=
dp)/(2._dp*real(kp_grid(3), kind=
dp))*reducing_factor(3)
3699 kpoints%xkp(1:3, i) = matmul(transpose(h_inv), kpt_latt(:))
3701 kpoints%wkp(i) = 2._dp*single_weight
3707 IF (do_extra_kpoints)
THEN
3709 single_weight = kpoint_weight_left/real(kp_grid(1)*kp_grid(2)*kp_grid(3)/8, kind=
dp)
3711 DO ix = 1, kp_grid(1)/2
3712 DO iy = 1, kp_grid(2)/2
3713 DO iz = 1, kp_grid(3)/2
3716 IF (2*ix - kp_grid(1)/2 - 1 < 0) cycle
3719 kpt_latt(1) = real(2*ix - kp_grid(1)/2 - 1, kind=
dp)/(real(kp_grid(1), kind=
dp))
3720 kpt_latt(2) = real(2*iy - kp_grid(2)/2 - 1, kind=
dp)/(real(kp_grid(2), kind=
dp))
3721 kpt_latt(3) = real(2*iz - kp_grid(3)/2 - 1, kind=
dp)/(real(kp_grid(3), kind=
dp))
3723 kpoints%xkp(1:3, i) = matmul(transpose(h_inv), kpt_latt(:))
3725 kpoints%wkp(i) = 2._dp*single_weight
3734 ALLOCATE (kpoints%kp_sym(kpoints%nkp))
3735 DO i = 1, kpoints%nkp
3736 NULLIFY (kpoints%kp_sym(i)%kpoint_sym)
3746 CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set)
3748 CALL calculate_kp_orbitals(qs_env_kp_gamma_only, kpoints,
"MONKHORST-PACK", nadd=nmo, mp_grid=kp_grid(1:3), &
3749 group_size_ext=para_env%num_pe)
3752 DEALLOCATE (qs_env_kp_gamma_only)
3757 END SUBROUTINE get_kpoints
3765 PURE SUBROUTINE average_degenerate_levels(vec_Sigma_c_gw, Eigenval_DFT, eps_eigenval)
3766 COMPLEX(KIND=dp),
DIMENSION(:, :, :), &
3767 INTENT(INOUT) :: vec_sigma_c_gw
3768 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: eigenval_dft
3769 REAL(kind=dp),
INTENT(IN) :: eps_eigenval
3771 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: avg_self_energy
3772 INTEGER :: degeneracy, first_degenerate_level, i_deg_level, i_level_gw, j_deg_level, jquad, &
3773 num_deg_levels, num_integ_points, num_levels_gw
3774 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: list_degenerate_levels
3776 num_levels_gw =
SIZE(vec_sigma_c_gw, 1)
3778 ALLOCATE (list_degenerate_levels(num_levels_gw))
3779 list_degenerate_levels = 1
3781 num_integ_points =
SIZE(vec_sigma_c_gw, 2)
3783 ALLOCATE (avg_self_energy(num_integ_points))
3785 DO i_level_gw = 2, num_levels_gw
3787 IF (abs(eigenval_dft(i_level_gw) - eigenval_dft(i_level_gw - 1)) < eps_eigenval)
THEN
3789 list_degenerate_levels(i_level_gw) = list_degenerate_levels(i_level_gw - 1)
3793 list_degenerate_levels(i_level_gw) = list_degenerate_levels(i_level_gw - 1) + 1
3799 num_deg_levels = list_degenerate_levels(num_levels_gw)
3801 DO i_deg_level = 1, num_deg_levels
3805 DO i_level_gw = 1, num_levels_gw
3807 IF (degeneracy == 0 .AND. i_deg_level == list_degenerate_levels(i_level_gw))
THEN
3809 first_degenerate_level = i_level_gw
3813 IF (i_deg_level == list_degenerate_levels(i_level_gw))
THEN
3815 degeneracy = degeneracy + 1
3821 DO jquad = 1, num_integ_points
3823 avg_self_energy(jquad) = sum(vec_sigma_c_gw(first_degenerate_level:first_degenerate_level + degeneracy - 1, jquad, 1)) &
3824 /real(degeneracy, kind=dp)
3828 DO j_deg_level = 0, degeneracy - 1
3830 vec_sigma_c_gw(first_degenerate_level + j_deg_level, :, 1) = avg_self_energy(:)
3836 END SUBROUTINE average_degenerate_levels
3859 SUBROUTINE fit_and_continuation_2pole(vec_gw_energ, vec_omega_fit_gw, &
3860 z_value, m_value, vec_Sigma_c_gw, vec_Sigma_x_minus_vxc_gw, &
3861 Eigenval, Eigenval_scf, n_level_gw, &
3862 gw_corr_lev_occ, gw_corr_lev_vir, num_poles, &
3863 num_fit_points, crossing_search, homo, stop_crit, &
3864 fermi_level_offset, do_gw_im_time)
3866 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: vec_gw_energ, vec_omega_fit_gw, z_value, &
3868 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: vec_sigma_c_gw
3869 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_sigma_x_minus_vxc_gw, eigenval, &
3871 INTEGER,
INTENT(IN) :: n_level_gw, gw_corr_lev_occ, &
3872 gw_corr_lev_vir, num_poles, &
3873 num_fit_points, crossing_search, homo
3874 REAL(kind=dp),
INTENT(IN) :: stop_crit, fermi_level_offset
3875 LOGICAL,
INTENT(IN) :: do_gw_im_time
3877 CHARACTER(LEN=*),
PARAMETER :: routinen =
'fit_and_continuation_2pole'
3879 COMPLEX(KIND=dp) :: func_val, rho1
3880 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: dlambda, dlambda_2, lambda, &
3881 lambda_without_offset, vec_b_gw, &
3883 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: mat_a_gw, mat_b_gw
3884 INTEGER :: handle4, ierr, iii, iiter, info, &
3885 integ_range, jjj, jquad, kkk, &
3886 max_iter_fit, n_level_gw_ref, num_var, &
3888 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ipiv
3889 LOGICAL :: could_exit
3890 REAL(kind=dp) :: chi2, chi2_old, delta, deriv_val_real, e_fermi, gw_energ, ldown, &
3891 level_energ_gw, lup, range_step, scalparam, sign_occ_virt, stat_error
3892 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: lambda_im, lambda_re, stat_errors, &
3893 vec_n_gw, vec_omega_fit_gw_sign
3894 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: mat_n_gw
3896 max_iter_fit = 10000
3898 num_var = 2*num_poles + 1
3899 ALLOCATE (lambda(num_var))
3901 ALLOCATE (lambda_without_offset(num_var))
3902 lambda_without_offset = z_zero
3903 ALLOCATE (lambda_re(num_var))
3905 ALLOCATE (lambda_im(num_var))
3908 ALLOCATE (vec_omega_fit_gw_sign(num_fit_points))
3910 IF (n_level_gw <= gw_corr_lev_occ)
THEN
3911 sign_occ_virt = -1.0_dp
3913 sign_occ_virt = 1.0_dp
3916 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
3918 DO jquad = 1, num_fit_points
3919 vec_omega_fit_gw_sign(jquad) = abs(vec_omega_fit_gw(jquad))*sign_occ_virt
3923 range_step = (vec_omega_fit_gw_sign(num_fit_points) - vec_omega_fit_gw_sign(1))/(num_poles - 1)
3924 DO iii = 1, num_poles
3925 lambda_im(2*iii + 1) = vec_omega_fit_gw_sign(1) + (iii - 1)*range_step
3927 range_step = (vec_omega_fit_gw_sign(num_fit_points) - vec_omega_fit_gw_sign(1))/num_poles
3928 DO iii = 1, num_poles
3929 lambda_re(2*iii + 1) = abs(vec_omega_fit_gw_sign(1) + (iii - 0.5_dp)*range_step)
3933 lambda(iii) = lambda_re(iii) + gaussi*lambda_im(iii)
3936 CALL calc_chi2(chi2_old, lambda, vec_sigma_c_gw, vec_omega_fit_gw_sign, num_poles, &
3937 num_fit_points, n_level_gw)
3939 ALLOCATE (mat_a_gw(num_poles + 1, num_poles + 1))
3940 ALLOCATE (vec_b_gw(num_poles + 1))
3941 ALLOCATE (ipiv(num_poles + 1))
3945 mat_a_gw(1:num_poles + 1, 1) = z_one
3946 integ_range = num_fit_points/num_poles
3947 DO kkk = 1, num_poles + 1
3948 xpos = (kkk - 1)*integ_range + 1
3949 xpos = min(xpos, num_fit_points)
3951 DO iii = 1, num_poles
3953 func_val = z_one/(gaussi*vec_omega_fit_gw_sign(xpos) - &
3954 cmplx(lambda_re(jjj + 1), lambda_im(jjj + 1), kind=dp))
3955 mat_a_gw(kkk, iii + 1) = func_val
3957 vec_b_gw(kkk) = vec_sigma_c_gw(n_level_gw, xpos)
3961 CALL zgetrf(num_poles + 1, num_poles + 1, mat_a_gw, num_poles + 1, ipiv, info)
3963 CALL zgetrs(
'N', num_poles + 1, 1, mat_a_gw, num_poles + 1, ipiv, vec_b_gw, num_poles + 1, info)
3965 lambda_re(1) = real(vec_b_gw(1))
3966 lambda_im(1) = aimag(vec_b_gw(1))
3967 DO iii = 1, num_poles
3969 lambda_re(jjj) = real(vec_b_gw(iii + 1))
3970 lambda_im(jjj) = aimag(vec_b_gw(iii + 1))
3973 DEALLOCATE (mat_a_gw)
3974 DEALLOCATE (vec_b_gw)
3977 ALLOCATE (mat_a_gw(num_var*2, num_var*2))
3978 ALLOCATE (mat_b_gw(num_fit_points, num_var*2))
3979 ALLOCATE (dlambda(num_fit_points))
3980 ALLOCATE (dlambda_2(num_fit_points))
3981 ALLOCATE (vec_b_gw(num_var*2))
3982 ALLOCATE (vec_b_gw_copy(num_var*2))
3983 ALLOCATE (ipiv(num_var*2))
3988 could_exit = .false.
3991 DO iiter = 1, max_iter_fit
3993 CALL timeset(routinen//
"_fit_loop_1", handle4)
3997 lambda(iii) = lambda_re(iii) + gaussi*lambda_im(iii)
4001 DO kkk = 1, num_fit_points
4002 func_val = lambda(1)
4003 DO iii = 1, num_poles
4005 func_val = func_val + lambda(jjj)/(vec_omega_fit_gw_sign(kkk)*gaussi - lambda(jjj + 1))
4007 dlambda(kkk) = vec_sigma_c_gw(n_level_gw, kkk) - func_val
4009 rho1 = sum(dlambda*dlambda)
4013 DO iii = 1, num_fit_points
4014 mat_b_gw(iii, 1) = 1.0_dp
4015 mat_b_gw(iii, num_var + 1) = gaussi
4017 DO iii = 1, num_poles
4019 DO kkk = 1, num_fit_points
4020 mat_b_gw(kkk, jjj) = 1.0_dp/(gaussi*vec_omega_fit_gw_sign(kkk) - lambda(jjj + 1))
4021 mat_b_gw(kkk, jjj + num_var) = gaussi/(gaussi*vec_omega_fit_gw_sign(kkk) - lambda(jjj + 1))
4022 mat_b_gw(kkk, jjj + 1) = lambda(jjj)/(gaussi*vec_omega_fit_gw_sign(kkk) - lambda(jjj + 1))**2
4023 mat_b_gw(kkk, jjj + 1 + num_var) = (-lambda_im(jjj) + gaussi*lambda_re(jjj))/ &
4024 (gaussi*vec_omega_fit_gw_sign(kkk) - lambda(jjj + 1))**2
4028 CALL timestop(handle4)
4030 CALL timeset(routinen//
"_fit_matmul_1", handle4)
4032 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, &
4033 z_zero, mat_a_gw, num_var*2)
4034 CALL timestop(handle4)
4036 CALL timeset(routinen//
"_fit_zgemv_1", handle4)
4037 CALL zgemv(
'C', num_fit_points, num_var*2, z_one, mat_b_gw, num_fit_points, dlambda, 1, &
4038 z_zero, vec_b_gw, 1)
4040 CALL timestop(handle4)
4043 DO iii = 1, num_var*2
4044 mat_a_gw(iii, iii) = mat_a_gw(iii, iii) + scalparam*mat_a_gw(iii, iii)
4051 CALL timeset(routinen//
"_fit_lin_eq_2", handle4)
4053 CALL zgetrf(2*num_var, 2*num_var, mat_a_gw, 2*num_var, ipiv, info)
4055 CALL zgetrs(
'N', 2*num_var, 1, mat_a_gw, 2*num_var, ipiv, vec_b_gw, 2*num_var, info)
4057 CALL timestop(handle4)
4060 lambda(iii) = lambda_re(iii) + gaussi*lambda_im(iii) + vec_b_gw(iii) + vec_b_gw(iii + num_var)
4064 CALL calc_chi2(chi2, lambda, vec_sigma_c_gw, vec_omega_fit_gw_sign, num_poles, &
4065 num_fit_points, n_level_gw)
4068 IF (chi2 < 1.0e-30_dp)
EXIT
4070 IF (chi2 < chi2_old)
THEN
4071 scalparam = max(scalparam/ldown, 1e-12_dp)
4073 lambda_re(iii) = lambda_re(iii) + real(vec_b_gw(iii) + vec_b_gw(iii + num_var))
4074 lambda_im(iii) = lambda_im(iii) + aimag(vec_b_gw(iii) + vec_b_gw(iii + num_var))
4076 IF (chi2_old/chi2 - 1.0_dp < stop_crit) could_exit = .true.
4079 scalparam = scalparam*lup
4081 IF (scalparam > 100.0_dp .AND. could_exit)
EXIT
4083 IF (scalparam > 1e+10_dp) scalparam = 1e-4_dp
4087 IF (.NOT. do_gw_im_time)
THEN
4091 func_val = lambda(1)
4092 DO iii = 1, num_poles
4095 func_val = func_val + lambda(jjj)/(-lambda(jjj + 1))
4098 lambda_re(1) = lambda_re(1) - real(func_val) + real(vec_sigma_c_gw(n_level_gw, num_fit_points))
4099 lambda_im(1) = lambda_im(1) - aimag(func_val) + aimag(vec_sigma_c_gw(n_level_gw, num_fit_points))
4103 lambda_without_offset(:) = lambda(:)
4106 lambda(iii) = cmplx(lambda_re(iii), lambda_im(iii), kind=dp)
4109 IF (do_gw_im_time)
THEN
4112 e_fermi = 0.5_dp*(eigenval(homo) + eigenval(homo + 1))
4116 IF (n_level_gw <= gw_corr_lev_occ)
THEN
4117 e_fermi = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset
4119 e_fermi = minval(eigenval(homo + 1:homo + gw_corr_lev_vir)) - fermi_level_offset
4124 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot .OR. &
4125 crossing_search == ri_rpa_g0w0_crossing_newton)
THEN
4128 func_val = lambda(1)
4129 z_value(n_level_gw) = 1.0_dp
4130 DO iii = 1, num_poles
4132 z_value(n_level_gw) = z_value(n_level_gw) + real(lambda(jjj)/ &
4133 (eigenval(n_level_gw_ref) - e_fermi - lambda(jjj + 1))**2)
4134 func_val = func_val + lambda(jjj)/(eigenval(n_level_gw_ref) - e_fermi - lambda(jjj + 1))
4137 m_value(n_level_gw) = 1.0_dp - z_value(n_level_gw)
4138 z_value(n_level_gw) = 1.0_dp/z_value(n_level_gw)
4139 gw_energ = real(func_val)
4140 vec_gw_energ(n_level_gw) = gw_energ
4143 IF (crossing_search == ri_rpa_g0w0_crossing_newton)
THEN
4145 level_energ_gw = (eigenval_scf(n_level_gw_ref) - &
4146 m_value(n_level_gw)*eigenval(n_level_gw_ref) + &
4147 vec_gw_energ(n_level_gw) + &
4148 vec_sigma_x_minus_vxc_gw(n_level_gw_ref))* &
4155 func_val = lambda(1)
4156 z_value(n_level_gw) = 1.0_dp
4157 DO iii = 1, num_poles
4159 func_val = func_val + lambda(jjj)/(level_energ_gw - e_fermi - lambda(jjj + 1))
4163 deriv_val_real = -1.0_dp
4164 DO iii = 1, num_poles
4166 deriv_val_real = deriv_val_real + real(lambda(jjj))/((abs(level_energ_gw - e_fermi - lambda(jjj + 1)))**2) &
4167 - (real(lambda(jjj))*(level_energ_gw - e_fermi) - real(lambda(jjj)*conjg(lambda(jjj + 1))))* &
4168 2.0_dp*(level_energ_gw - e_fermi - real(lambda(jjj + 1)))/ &
4169 ((abs(level_energ_gw - e_fermi - lambda(jjj + 1)))**2)
4173 delta = (eigenval_scf(n_level_gw_ref) + vec_sigma_x_minus_vxc_gw(n_level_gw_ref) + real(func_val) - level_energ_gw)/ &
4176 level_energ_gw = level_energ_gw - delta
4178 IF (abs(delta) < 1.0e-08)
EXIT
4184 vec_gw_energ(n_level_gw) = real(func_val)
4185 z_value(n_level_gw) = 1.0_dp
4186 m_value(n_level_gw) = 0.0_dp
4191 cpabort(
"Only NONE, ZSHOT and NEWTON implemented for 2-pole model")
4201 CALL calc_chi2(chi2, lambda_without_offset, vec_sigma_c_gw, vec_omega_fit_gw_sign, num_poles, &
4202 num_fit_points, n_level_gw)
4205 stat_error = sqrt(chi2/num_fit_points)
4208 ALLOCATE (vec_n_gw(num_var*2))
4211 ALLOCATE (mat_n_gw(num_var*2, num_var*2))
4214 DO iii = 1, num_var*2
4215 CALL calc_mat_n(vec_n_gw(iii), lambda_without_offset, vec_sigma_c_gw, vec_omega_fit_gw_sign, &
4216 iii, iii, num_poles, num_fit_points, n_level_gw, 0.001_dp)
4219 DO iii = 1, num_var*2
4220 DO jjj = 1, num_var*2
4221 CALL calc_mat_n(mat_n_gw(iii, jjj), lambda_without_offset, vec_sigma_c_gw, vec_omega_fit_gw_sign, &
4222 iii, jjj, num_poles, num_fit_points, n_level_gw, 0.001_dp)
4226 CALL dgetrf(2*num_var, 2*num_var, mat_n_gw, 2*num_var, ipiv, info)
4229 CALL dgetri(2*num_var, mat_n_gw, 2*num_var, ipiv, vec_b_gw, 2*num_var, info)
4231 ALLOCATE (stat_errors(2*num_var))
4232 stat_errors = 0.0_dp
4234 DO iii = 1, 2*num_var
4235 stat_errors(iii) = sqrt(abs(mat_n_gw(iii, iii)))*stat_error
4238 DEALLOCATE (mat_n_gw)
4239 DEALLOCATE (vec_n_gw)
4240 DEALLOCATE (mat_a_gw)
4241 DEALLOCATE (mat_b_gw)
4242 DEALLOCATE (stat_errors)
4243 DEALLOCATE (dlambda)
4244 DEALLOCATE (dlambda_2)
4245 DEALLOCATE (vec_b_gw)
4246 DEALLOCATE (vec_b_gw_copy)
4248 DEALLOCATE (vec_omega_fit_gw_sign)
4250 DEALLOCATE (lambda_without_offset)
4251 DEALLOCATE (lambda_re)
4252 DEALLOCATE (lambda_im)
4254 END SUBROUTINE fit_and_continuation_2pole
4290 z_value, m_value, vec_Sigma_c_gw, vec_Sigma_x_minus_vxc_gw, &
4291 Eigenval, Eigenval_scf, do_hedin_shift, n_level_gw, &
4292 gw_corr_lev_occ, gw_corr_lev_vir, &
4293 nparam_pade, num_fit_points, crossing_search, homo, &
4294 fermi_level_offset, do_gw_im_time, print_self_energy, count_ev_sc_GW, &
4295 vec_gw_dos, dos_lower_bound, dos_precision, ndos, &
4296 min_level_self_energy, max_level_self_energy, &
4297 dos_eta, dos_min, dos_max, e_fermi_ext)
4300 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: vec_gw_energ
4301 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_omega_fit_gw
4302 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: z_value, m_value
4303 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: vec_sigma_c_gw
4304 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_sigma_x_minus_vxc_gw, eigenval, &
4306 LOGICAL,
INTENT(IN) :: do_hedin_shift
4307 INTEGER,
INTENT(IN) :: n_level_gw, gw_corr_lev_occ, &
4308 gw_corr_lev_vir, nparam_pade, &
4309 num_fit_points, crossing_search, homo
4310 REAL(kind=dp),
INTENT(IN) :: fermi_level_offset
4311 LOGICAL,
INTENT(IN) :: do_gw_im_time, print_self_energy
4312 INTEGER,
INTENT(IN) :: count_ev_sc_gw
4313 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:),
OPTIONAL :: vec_gw_dos
4314 REAL(kind=dp),
OPTIONAL :: dos_lower_bound, dos_precision
4315 INTEGER,
INTENT(IN),
OPTIONAL :: ndos, min_level_self_energy, &
4316 max_level_self_energy
4317 REAL(kind=dp),
OPTIONAL :: dos_eta
4318 INTEGER,
INTENT(IN),
OPTIONAL :: dos_min, dos_max
4319 REAL(kind=dp),
OPTIONAL :: e_fermi_ext
4321 CHARACTER(LEN=*),
PARAMETER :: routinen =
'continuation_pade'
4323 CHARACTER(LEN=5) :: string_level
4324 CHARACTER(len=default_path_length) :: filename
4325 COMPLEX(KIND=dp) :: sigma_c_pade, sigma_c_pade_im_freq
4326 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: coeff_pade, omega_points_pade, &
4328 INTEGER :: handle, i_omega, idos, iunit, jquad, &
4329 n_level_gw_ref, num_omega
4330 REAL(kind=dp) :: e_fermi, energy_val, hedin_shift, &
4331 level_energ_gw_start, omega, &
4332 omega_dos, omega_dos_pade_eval, &
4334 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: vec_omega_fit_gw_sign, &
4335 vec_omega_fit_gw_sign_reorder, &
4336 vec_sigma_imag, vec_sigma_real
4337 TYPE(cp_logger_type),
POINTER :: logger
4339 CALL timeset(routinen, handle)
4341 ALLOCATE (vec_omega_fit_gw_sign(num_fit_points))
4343 IF (n_level_gw <= gw_corr_lev_occ)
THEN
4344 sign_occ_virt = -1.0_dp
4346 sign_occ_virt = 1.0_dp
4349 DO jquad = 1, num_fit_points
4350 vec_omega_fit_gw_sign(jquad) = abs(vec_omega_fit_gw(jquad))*sign_occ_virt
4353 IF (do_gw_im_time)
THEN
4356 e_fermi = 0.5_dp*(eigenval(homo) + eigenval(homo + 1))
4360 IF (n_level_gw <= gw_corr_lev_occ)
THEN
4361 e_fermi = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo)) + fermi_level_offset
4363 e_fermi = minval(eigenval(homo + 1:homo + gw_corr_lev_vir)) - fermi_level_offset
4367 IF (
PRESENT(e_fermi_ext)) e_fermi = e_fermi_ext
4369 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
4372 ALLOCATE (sigma_c_gw_reorder(num_fit_points))
4373 ALLOCATE (vec_omega_fit_gw_sign_reorder(num_fit_points))
4375 IF (do_gw_im_time)
THEN
4376 DO jquad = 1, num_fit_points
4377 sigma_c_gw_reorder(jquad) = vec_sigma_c_gw(n_level_gw, jquad)
4378 vec_omega_fit_gw_sign_reorder(jquad) = vec_omega_fit_gw_sign(jquad)
4381 DO jquad = 1, num_fit_points
4382 sigma_c_gw_reorder(jquad) = vec_sigma_c_gw(n_level_gw, num_fit_points - jquad + 1)
4383 vec_omega_fit_gw_sign_reorder(jquad) = vec_omega_fit_gw_sign(num_fit_points - jquad + 1)
4388 ALLOCATE (coeff_pade(nparam_pade))
4389 ALLOCATE (omega_points_pade(nparam_pade))
4391 CALL get_pade_parameters(sigma_c_gw_reorder, vec_omega_fit_gw_sign_reorder, &
4392 num_fit_points, nparam_pade, omega_points_pade, coeff_pade)
4395 IF ((crossing_search == ri_rpa_g0w0_crossing_bisection) .OR. &
4396 (crossing_search == ri_rpa_g0w0_crossing_newton))
THEN
4397 energy_val = eigenval(n_level_gw_ref) - e_fermi
4398 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, &
4399 coeff_pade, sigma_c_pade)
4400 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, &
4401 coeff_pade, z_value(n_level_gw), m_value(n_level_gw))
4402 level_energ_gw_start = (eigenval_scf(n_level_gw_ref) - &
4403 m_value(n_level_gw)*eigenval(n_level_gw_ref) + &
4404 REAL(sigma_c_pade) + &
4405 vec_sigma_x_minus_vxc_gw(n_level_gw_ref))* &
4409 hedin_shift = 0.0_dp
4410 IF (do_hedin_shift) hedin_shift = real(sigma_c_pade) + &
4411 vec_sigma_x_minus_vxc_gw(n_level_gw_ref) &
4412 - eigenval(n_level_gw_ref) + eigenval_scf(n_level_gw_ref)
4415 IF (
PRESENT(min_level_self_energy) .AND.
PRESENT(max_level_self_energy))
THEN
4416 IF (n_level_gw_ref >= min_level_self_energy .AND. &
4417 n_level_gw_ref <= max_level_self_energy)
THEN
4418 ALLOCATE (vec_sigma_real(ndos))
4419 ALLOCATE (vec_sigma_imag(ndos))
4420 WRITE (string_level,
"(I4)") n_level_gw_ref
4421 string_level = adjustl(string_level)
4430 IF (
PRESENT(ndos))
THEN
4433 cpassert(.NOT. do_hedin_shift)
4434 logger => cp_get_default_logger()
4435 IF (logger%para_env%is_source())
THEN
4436 iunit = cp_logger_get_default_unit_nr()
4441 omega_dos = dos_lower_bound + real(idos - 1, kind=dp)*dos_precision
4442 omega_dos_pade_eval = omega_dos - e_fermi
4443 CALL evaluate_pade_function(omega_dos_pade_eval, nparam_pade, omega_points_pade, &
4444 coeff_pade, sigma_c_pade)
4446 IF (n_level_gw_ref >= min_level_self_energy .AND. &
4447 n_level_gw_ref <= max_level_self_energy .AND. iunit > 0)
THEN
4449 vec_sigma_real(idos) = (real(sigma_c_pade))
4450 vec_sigma_imag(idos) = (aimag(sigma_c_pade))
4454 IF (n_level_gw_ref >= dos_min .AND. &
4455 (n_level_gw_ref <= dos_max .OR. dos_max == 0))
THEN
4456 vec_gw_dos(idos) = vec_gw_dos(idos) + &
4457 (abs(aimag(sigma_c_pade)) + dos_eta) &
4459 (omega_dos - eigenval_scf(n_level_gw_ref) - &
4460 (real(sigma_c_pade) + vec_sigma_x_minus_vxc_gw(n_level_gw_ref)) &
4462 + (abs(aimag(sigma_c_pade)) + dos_eta)**2 &
4470 IF (
PRESENT(min_level_self_energy) .AND.
PRESENT(max_level_self_energy))
THEN
4471 logger => cp_get_default_logger()
4472 IF (logger%para_env%is_source())
THEN
4473 iunit = cp_logger_get_default_unit_nr()
4477 IF (n_level_gw_ref >= min_level_self_energy .AND. &
4478 n_level_gw_ref <= max_level_self_energy .AND. iunit > 0)
THEN
4480 CALL open_file(
'self_energy_re_'//trim(string_level)//
'.dat', unit_number=iunit, &
4481 file_status=
"UNKNOWN", file_action=
"WRITE")
4483 omega_dos = dos_lower_bound + real(idos - 1, kind=dp)*dos_precision
4484 WRITE (iunit,
'(F17.10, F17.10)') omega_dos*evolt, vec_sigma_real(idos)*evolt
4487 CALL close_file(iunit)
4489 CALL open_file(
'self_energy_im_'//trim(string_level)//
'.dat', unit_number=iunit, &
4490 file_status=
"UNKNOWN", file_action=
"WRITE")
4492 omega_dos = dos_lower_bound + real(idos - 1, kind=dp)*dos_precision
4493 WRITE (iunit,
'(F17.10, F17.10)') omega_dos*evolt, vec_sigma_imag(idos)*evolt
4496 CALL close_file(iunit)
4498 DEALLOCATE (vec_sigma_real)
4499 DEALLOCATE (vec_sigma_imag)
4504 SELECT CASE (crossing_search)
4505 CASE (ri_rpa_g0w0_crossing_z_shot)
4507 cpassert(.NOT. do_hedin_shift)
4508 energy_val = eigenval(n_level_gw_ref) - e_fermi
4509 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, &
4510 coeff_pade, sigma_c_pade)
4511 vec_gw_energ(n_level_gw) = real(sigma_c_pade)
4513 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, &
4514 coeff_pade, z_value(n_level_gw), m_value(n_level_gw))
4516 CASE (ri_rpa_g0w0_crossing_bisection)
4517 CALL get_sigma_c_bisection_pade(vec_gw_energ(n_level_gw), eigenval_scf(n_level_gw_ref), &
4518 vec_sigma_x_minus_vxc_gw(n_level_gw_ref), e_fermi, &
4519 nparam_pade, omega_points_pade, coeff_pade, &
4520 level_energ_gw_start, hedin_shift)
4521 z_value(n_level_gw) = 1.0_dp
4522 m_value(n_level_gw) = 0.0_dp
4524 CASE (ri_rpa_g0w0_crossing_newton)
4525 CALL get_sigma_c_newton_pade(vec_gw_energ(n_level_gw), eigenval_scf(n_level_gw_ref), &
4526 vec_sigma_x_minus_vxc_gw(n_level_gw_ref), e_fermi, &
4527 nparam_pade, omega_points_pade, coeff_pade, &
4528 level_energ_gw_start, hedin_shift)
4529 z_value(n_level_gw) = 1.0_dp
4530 m_value(n_level_gw) = 0.0_dp
4533 cpabort(
"Only Z_SHOT, NEWTON, and BISECTION crossing search implemented.")
4536 IF (print_self_energy)
THEN
4538 IF (count_ev_sc_gw == 1)
THEN
4540 IF (n_level_gw_ref < 10)
THEN
4541 WRITE (filename,
"(A26,I1)")
"G0W0_self_energy_level_000", n_level_gw_ref
4542 ELSE IF (n_level_gw_ref < 100)
THEN
4543 WRITE (filename,
"(A25,I2)")
"G0W0_self_energy_level_00", n_level_gw_ref
4544 ELSE IF (n_level_gw_ref < 1000)
THEN
4545 WRITE (filename,
"(A24,I3)")
"G0W0_self_energy_level_0", n_level_gw_ref
4547 WRITE (filename,
"(A23,I4)")
"G0W0_self_energy_level_", n_level_gw_ref
4552 IF (n_level_gw_ref < 10)
THEN
4553 WRITE (filename,
"(A11,I1,A22,I1)")
"evGW_cycle_", count_ev_sc_gw, &
4554 "_self_energy_level_000", n_level_gw_ref
4555 ELSE IF (n_level_gw_ref < 100)
THEN
4556 WRITE (filename,
"(A11,I1,A21,I2)")
"evGW_cycle_", count_ev_sc_gw, &
4557 "_self_energy_level_00", n_level_gw_ref
4558 ELSE IF (n_level_gw_ref < 1000)
THEN
4559 WRITE (filename,
"(A11,I1,A20,I3)")
"evGW_cycle_", count_ev_sc_gw, &
4560 "_self_energy_level_0", n_level_gw_ref
4562 WRITE (filename,
"(A11,I1,A19,I4)")
"evGW_cycle_", count_ev_sc_gw, &
4563 "_self_energy_level_", n_level_gw_ref
4568 logger => cp_get_default_logger()
4569 IF (logger%para_env%is_source())
THEN
4570 iunit = cp_logger_get_default_unit_nr()
4574 CALL open_file(trim(filename), unit_number=iunit, file_status=
"UNKNOWN", file_action=
"WRITE")
4578 WRITE (iunit,
"(2A42)")
" omega (eV) Sigma(omega) (eV) ", &
4579 " omega - e_n^DFT - Sigma_n^x - v_n^xc (eV)"
4581 DO i_omega = 0, num_omega
4583 omega = -50.0_dp/evolt + real(i_omega, kind=dp)/real(num_omega, kind=dp)*100.0_dp/evolt
4585 CALL evaluate_pade_function(omega - e_fermi, nparam_pade, omega_points_pade, &
4586 coeff_pade, sigma_c_pade)
4588 WRITE (iunit,
"(F12.2,2F17.5)") omega*evolt, real(sigma_c_pade)*evolt, &
4589 (omega - eigenval_scf(n_level_gw_ref) - vec_sigma_x_minus_vxc_gw(n_level_gw_ref))*evolt
4593 WRITE (iunit,
"(A51,A39)")
" w (eV) Re(Sigma(i*w)) (eV) Im(Sigma(i*w)) (eV) ", &
4594 " Re(Fit(i*w)) (eV) Im(Fit(iw)) (eV)"
4596 DO jquad = 1, num_fit_points
4598 CALL evaluate_pade_function(vec_omega_fit_gw_sign_reorder(jquad), &
4599 nparam_pade, omega_points_pade, &
4600 coeff_pade, sigma_c_pade_im_freq, do_imag_freq=.true.)
4602 WRITE (iunit,
"(F12.2,4F17.5)") vec_omega_fit_gw_sign_reorder(jquad)*evolt, &
4603 REAL(sigma_c_gw_reorder(jquad)*evolt), &
4604 aimag(sigma_c_gw_reorder(jquad)*evolt), &
4605 REAL(sigma_c_pade_im_freq*evolt), &
4606 aimag(sigma_c_pade_im_freq*evolt)
4610 CALL close_file(iunit)
4614 DEALLOCATE (vec_omega_fit_gw_sign)
4615 DEALLOCATE (sigma_c_gw_reorder)
4616 DEALLOCATE (vec_omega_fit_gw_sign_reorder)
4617 DEALLOCATE (coeff_pade, omega_points_pade)
4619 CALL timestop(handle)
4633 PURE SUBROUTINE get_pade_parameters(y, x, num_fit_points, nparam, xpoints, coeff)
4635 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: y
4636 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: x
4637 INTEGER,
INTENT(IN) :: num_fit_points, nparam
4638 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: xpoints, coeff
4640 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: ypoints
4641 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: g_mat
4642 INTEGER :: idat, iparam, nstep
4644 nstep = int(num_fit_points/(nparam - 1))
4646 ALLOCATE (ypoints(nparam))
4649 DO iparam = 1, nparam - 1
4650 xpoints(iparam) = gaussi*x(idat)
4651 ypoints(iparam) = y(idat)
4654 xpoints(nparam) = gaussi*x(num_fit_points)
4655 ypoints(nparam) = y(num_fit_points)
4659 ALLOCATE (g_mat(nparam, nparam))
4660 g_mat(:, 1) = ypoints(:)
4661 DO iparam = 2, nparam
4662 DO idat = iparam, nparam
4663 g_mat(idat, iparam) = (g_mat(iparam - 1, iparam - 1) - g_mat(idat, iparam - 1))/ &
4664 ((xpoints(idat) - xpoints(iparam - 1))*g_mat(idat, iparam - 1))
4668 DO iparam = 1, nparam
4669 coeff(iparam) = g_mat(iparam, iparam)
4672 DEALLOCATE (ypoints)
4675 END SUBROUTINE get_pade_parameters
4686 PURE SUBROUTINE evaluate_pade_function(x_val, nparam, xpoints, coeff, func_val, do_imag_freq)
4688 REAL(kind=dp),
INTENT(IN) :: x_val
4689 INTEGER,
INTENT(IN) :: nparam
4690 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: xpoints, coeff
4691 COMPLEX(KIND=dp),
INTENT(OUT) :: func_val
4692 LOGICAL,
INTENT(IN),
OPTIONAL :: do_imag_freq
4695 LOGICAL :: my_do_imag_freq
4697 my_do_imag_freq = .false.
4698 IF (
PRESENT(do_imag_freq)) my_do_imag_freq = do_imag_freq
4701 DO iparam = nparam, 2, -1
4702 IF (my_do_imag_freq)
THEN
4703 func_val = z_one + coeff(iparam)*(gaussi*x_val - xpoints(iparam - 1))/func_val
4705 func_val = z_one + coeff(iparam)*(x_val*z_one - xpoints(iparam - 1))/func_val
4709 func_val = coeff(1)/func_val
4711 END SUBROUTINE evaluate_pade_function
4722 PURE SUBROUTINE get_z_and_m_value_pade(x_val, nparam, xpoints, coeff, z_value, m_value)
4724 REAL(kind=dp),
INTENT(IN) :: x_val
4725 INTEGER,
INTENT(IN) :: nparam
4726 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: xpoints, coeff
4727 REAL(kind=dp),
INTENT(OUT),
OPTIONAL :: z_value, m_value
4729 COMPLEX(KIND=dp) :: denominator, dev_denominator, &
4730 dev_numerator, dev_val, func_val, &
4736 DO iparam = nparam, 2, -1
4737 numerator = coeff(iparam)*(x_val*z_one - xpoints(iparam - 1))
4738 dev_numerator = coeff(iparam)*z_one
4739 denominator = func_val
4740 dev_denominator = dev_val
4741 dev_val = dev_numerator/denominator - (numerator*dev_denominator)/(denominator**2)
4742 func_val = z_one + coeff(iparam)*(x_val*z_one - xpoints(iparam - 1))/func_val
4745 dev_val = -1.0_dp*coeff(1)/(func_val**2)*dev_val
4746 func_val = coeff(1)/func_val
4748 IF (
PRESENT(z_value))
THEN
4749 z_value = 1.0_dp - real(dev_val)
4750 z_value = 1.0_dp/z_value
4752 IF (
PRESENT(m_value)) m_value = real(dev_val)
4754 END SUBROUTINE get_z_and_m_value_pade
4768 SUBROUTINE get_sigma_c_bisection_pade(gw_energ, Eigenval_scf, Sigma_x_minus_vxc_gw, e_fermi, &
4769 nparam_pade, omega_points_pade, coeff_pade, start_val, &
4772 REAL(kind=dp),
INTENT(OUT) :: gw_energ
4773 REAL(kind=dp),
INTENT(IN) :: eigenval_scf, sigma_x_minus_vxc_gw, &
4775 INTEGER,
INTENT(IN) :: nparam_pade
4776 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: omega_points_pade, coeff_pade
4777 REAL(kind=dp),
INTENT(IN) :: start_val, hedin_shift
4779 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_sigma_c_bisection_pade'
4781 COMPLEX(KIND=dp) :: sigma_c
4782 INTEGER :: handle, icount
4783 REAL(kind=dp) :: delta, energy_val, qp_energy, &
4784 qp_energy_old, threshold
4786 CALL timeset(routinen, handle)
4788 threshold = 1.0e-7_dp
4790 qp_energy = start_val
4791 qp_energy_old = start_val
4795 DO WHILE (abs(delta) > threshold)
4797 qp_energy = qp_energy_old + 0.5_dp*delta
4798 qp_energy_old = qp_energy
4799 energy_val = qp_energy - e_fermi - hedin_shift
4800 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, &
4801 coeff_pade, sigma_c)
4802 qp_energy = eigenval_scf + real(sigma_c) + sigma_x_minus_vxc_gw
4803 delta = qp_energy - qp_energy_old
4805 IF (icount > 500)
EXIT
4808 gw_energ = real(sigma_c)
4810 CALL timestop(handle)
4812 END SUBROUTINE get_sigma_c_bisection_pade
4826 SUBROUTINE get_sigma_c_newton_pade(gw_energ, Eigenval_scf, Sigma_x_minus_vxc_gw, e_fermi, &
4827 nparam_pade, omega_points_pade, coeff_pade, start_val, &
4830 REAL(kind=dp),
INTENT(OUT) :: gw_energ
4831 REAL(kind=dp),
INTENT(IN) :: eigenval_scf, sigma_x_minus_vxc_gw, &
4833 INTEGER,
INTENT(IN) :: nparam_pade
4834 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: omega_points_pade, coeff_pade
4835 REAL(kind=dp),
INTENT(IN) :: start_val, hedin_shift
4837 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_sigma_c_newton_pade'
4839 COMPLEX(KIND=dp) :: sigma_c
4840 INTEGER :: handle, icount
4841 REAL(kind=dp) :: delta, energy_val, m_value, qp_energy, &
4842 qp_energy_old, threshold
4844 CALL timeset(routinen, handle)
4846 threshold = 1.0e-7_dp
4848 qp_energy = start_val
4849 qp_energy_old = start_val
4853 DO WHILE (abs(delta) > threshold)
4855 energy_val = qp_energy - e_fermi - hedin_shift
4856 CALL evaluate_pade_function(energy_val, nparam_pade, omega_points_pade, &
4857 coeff_pade, sigma_c)
4859 CALL get_z_and_m_value_pade(energy_val, nparam_pade, omega_points_pade, &
4860 coeff_pade, m_value=m_value)
4861 qp_energy_old = qp_energy
4862 qp_energy = qp_energy - (eigenval_scf + sigma_x_minus_vxc_gw + real(sigma_c) - qp_energy)/ &
4864 delta = qp_energy - qp_energy_old
4866 IF (icount > 500)
EXIT
4869 gw_energ = real(sigma_c)
4871 CALL timestop(handle)
4873 END SUBROUTINE get_sigma_c_newton_pade
4902 SUBROUTINE print_and_update_for_ev_sc(vec_gw_energ, &
4903 z_value, m_value, vec_Sigma_x_minus_vxc_gw, Eigenval, &
4904 Eigenval_last, Eigenval_scf, &
4905 gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, &
4906 crossing_search, homo, unit_nr, count_ev_sc_GW, count_sc_GW0, &
4907 ikp, nkp_self_energy, kpoints, ispin, E_VBM_GW, E_CBM_GW, &
4908 E_VBM_SCF, E_CBM_SCF)
4910 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_gw_energ, z_value, m_value
4911 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: vec_sigma_x_minus_vxc_gw, eigenval, &
4912 eigenval_last, eigenval_scf
4913 INTEGER,
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, crossing_search, &
4914 homo, unit_nr, count_ev_sc_gw, count_sc_gw0, ikp, nkp_self_energy
4915 TYPE(kpoint_type),
INTENT(IN),
POINTER :: kpoints
4916 INTEGER,
INTENT(IN) :: ispin
4917 REAL(kind=dp),
INTENT(INOUT),
OPTIONAL :: e_vbm_gw, e_cbm_gw, e_vbm_scf, e_cbm_scf
4919 CHARACTER(LEN=*),
PARAMETER :: routinen =
'print_and_update_for_ev_sc'
4921 CHARACTER(4) :: occ_virt
4922 INTEGER :: handle, n_level_gw, n_level_gw_ref
4923 LOGICAL :: do_alpha, do_beta, do_closed_shell, &
4924 do_kpoints, is_energy_okay
4925 REAL(kind=dp) :: e_gap_gw, e_homo_gw, e_homo_scf, &
4926 e_lumo_gw, e_lumo_scf, new_energy
4928 CALL timeset(routinen, handle)
4930 do_alpha = (ispin == 1)
4931 do_beta = (ispin == 2)
4932 do_closed_shell = .NOT. (do_alpha .OR. do_beta)
4933 do_kpoints = (nkp_self_energy > 1)
4935 eigenval_last(:) = eigenval(:)
4937 IF (unit_nr > 0)
THEN
4939 IF (count_ev_sc_gw == 1 .AND. count_sc_gw0 == 1 .AND. ikp == 1)
THEN
4941 WRITE (unit_nr, *)
' '
4943 IF (do_alpha .OR. do_closed_shell)
THEN
4944 WRITE (unit_nr, *)
' '
4945 WRITE (unit_nr,
'(T3,A)')
'******************************************************************************'
4946 WRITE (unit_nr,
'(T3,A)')
'** **'
4947 WRITE (unit_nr,
'(T3,A)')
'** GW QUASIPARTICLE ENERGIES **'
4948 WRITE (unit_nr,
'(T3,A)')
'** **'
4949 WRITE (unit_nr,
'(T3,A)')
'******************************************************************************'
4950 WRITE (unit_nr,
'(T3,A)')
' '
4951 WRITE (unit_nr,
'(T3,A)')
' '
4952 WRITE (unit_nr,
'(T3,A)')
'The GW quasiparticle energies are calculated according to: '
4954 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot)
THEN
4955 WRITE (unit_nr,
'(T3,A)')
'E_GW = E_SCF + Z * ( Sigc(E_SCF) + Sigx - vxc )'
4957 WRITE (unit_nr,
'(T3,A)')
' '
4958 WRITE (unit_nr,
'(T3,A)')
' E_GW = E_SCF + Sigc(E_GW) + Sigx - vxc '
4959 WRITE (unit_nr,
'(T3,A)')
' '
4960 WRITE (unit_nr,
'(T3,A)')
'Upper equation is solved self-consistently for E_GW, see Eq. (12) in J. Phys.'
4961 WRITE (unit_nr,
'(T3,A)')
'Chem. Lett. 9, 306 (2018), doi: 10.1021/acs.jpclett.7b02740'
4963 WRITE (unit_nr, *)
' '
4964 WRITE (unit_nr, *)
' '
4965 WRITE (unit_nr,
'(T3,A)')
'------------'
4966 WRITE (unit_nr,
'(T3,A)')
'G0W0 results'
4967 WRITE (unit_nr,
'(T3,A)')
'------------'
4971 IF (.NOT. do_kpoints)
THEN
4973 WRITE (unit_nr, *)
' '
4974 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
4975 WRITE (unit_nr,
'(T3,A)')
'GW quasiparticle energies of alpha spins'
4976 WRITE (unit_nr,
'(T3,A)')
'----------------------------------------'
4977 ELSE IF (do_beta)
THEN
4978 WRITE (unit_nr, *)
' '
4979 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
4980 WRITE (unit_nr,
'(T3,A)')
'GW quasiparticle energies of beta spins'
4981 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
4987 IF (count_ev_sc_gw > 1)
THEN
4988 WRITE (unit_nr, *)
' '
4989 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
4990 WRITE (unit_nr,
'(T3,A,I4)')
'Eigenvalue-selfconsistency cycle: ', count_ev_sc_gw
4991 WRITE (unit_nr,
'(T3,A)')
'---------------------------------------'
4994 IF (count_sc_gw0 > 1)
THEN
4995 WRITE (unit_nr,
'(T3,A)')
'----------------------------------'
4996 WRITE (unit_nr,
'(T3,A,I4)')
'scGW0 selfconsistency cycle: ', count_sc_gw0
4997 WRITE (unit_nr,
'(T3,A)')
'----------------------------------'
5000 IF (do_kpoints)
THEN
5001 WRITE (unit_nr, *)
' '
5002 WRITE (unit_nr,
'(T3,A7,I3,A3,I3,A8,3F7.3,A12,3F7.3)')
'Kpoint ', ikp,
' /', nkp_self_energy, &
5003 ' xkp =', kpoints%xkp(1, ikp), kpoints%xkp(2, ikp), kpoints%xkp(3, ikp), &
5004 ' and xkp =', -kpoints%xkp(1, ikp), -kpoints%xkp(2, ikp), -kpoints%xkp(3, ikp)
5005 WRITE (unit_nr,
'(T3,A72)')
'(Relative Brillouin zone size: [-0.5, 0.5] x [-0.5, 0.5] x [-0.5, 0.5])'
5006 WRITE (unit_nr, *)
' '
5008 WRITE (unit_nr,
'(T3,A)')
'GW quasiparticle energies of alpha spins:'
5009 ELSE IF (do_beta)
THEN
5010 WRITE (unit_nr,
'(T3,A)')
'GW quasiparticle energies of beta spins:'
5016 DO n_level_gw = 1, gw_corr_lev_tot
5018 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
5020 new_energy = (eigenval_scf(n_level_gw_ref) - &
5021 m_value(n_level_gw)*eigenval(n_level_gw_ref) + &
5022 vec_gw_energ(n_level_gw) + &
5023 vec_sigma_x_minus_vxc_gw(n_level_gw_ref))* &
5026 is_energy_okay = .true.
5028 IF (n_level_gw_ref > homo .AND. new_energy < eigenval(homo))
THEN
5029 is_energy_okay = .false.
5032 IF (is_energy_okay)
THEN
5033 eigenval(n_level_gw_ref) = new_energy
5038 IF (unit_nr > 0)
THEN
5039 WRITE (unit_nr,
'(T3,A)')
' '
5040 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot)
THEN
5041 WRITE (unit_nr,
'(T13,2A)')
'MO E_SCF (eV) Sigc (eV) Sigx-vxc (eV) Z E_GW (eV)'
5043 WRITE (unit_nr,
'(T3,2A)')
'Molecular orbital E_SCF (eV) Sigc (eV) Sigx-vxc (eV) E_GW (eV)'
5047 DO n_level_gw = 1, gw_corr_lev_tot
5048 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
5049 IF (n_level_gw <= gw_corr_lev_occ)
THEN
5055 IF (unit_nr > 0)
THEN
5056 IF (crossing_search == ri_rpa_g0w0_crossing_z_shot)
THEN
5057 WRITE (unit_nr,
'(T3,I4,3A,5F13.4)') &
5058 n_level_gw_ref,
' ( ', occ_virt,
') ', &
5059 eigenval_last(n_level_gw_ref)*evolt, &
5060 vec_gw_energ(n_level_gw)*evolt, &
5061 vec_sigma_x_minus_vxc_gw(n_level_gw_ref)*evolt, &
5062 z_value(n_level_gw), &
5063 eigenval(n_level_gw_ref)*evolt
5065 WRITE (unit_nr,
'(T3,I4,3A,4F16.4)') &
5066 n_level_gw_ref,
' ( ', occ_virt,
') ', &
5067 eigenval_last(n_level_gw_ref)*evolt, &
5068 vec_gw_energ(n_level_gw)*evolt, &
5069 vec_sigma_x_minus_vxc_gw(n_level_gw_ref)*evolt, &
5070 eigenval(n_level_gw_ref)*evolt
5075 e_homo_scf = maxval(eigenval_last(homo - gw_corr_lev_occ + 1:homo))
5076 e_lumo_scf = minval(eigenval_last(homo + 1:homo + gw_corr_lev_virt))
5078 e_homo_gw = maxval(eigenval(homo - gw_corr_lev_occ + 1:homo))
5079 e_lumo_gw = minval(eigenval(homo + 1:homo + gw_corr_lev_virt))
5080 e_gap_gw = e_lumo_gw - e_homo_gw
5082 IF (
PRESENT(e_vbm_scf) .AND.
PRESENT(e_cbm_scf) .AND. &
5083 PRESENT(e_vbm_gw) .AND.
PRESENT(e_cbm_gw))
THEN
5084 IF (e_homo_scf > e_vbm_scf) e_vbm_scf = e_homo_scf
5085 IF (e_lumo_scf < e_cbm_scf) e_cbm_scf = e_lumo_scf
5086 IF (e_homo_gw > e_vbm_gw) e_vbm_gw = e_homo_gw
5087 IF (e_lumo_gw < e_cbm_gw) e_cbm_gw = e_lumo_gw
5090 IF (unit_nr > 0)
THEN
5092 IF (do_kpoints)
THEN
5093 IF (do_closed_shell)
THEN
5094 WRITE (unit_nr,
'(T3,A)')
' '
5095 WRITE (unit_nr,
'(T3,A,F42.4)')
'GW direct gap at current kpoint (eV)', e_gap_gw*evolt
5096 ELSE IF (do_alpha)
THEN
5097 WRITE (unit_nr,
'(T3,A)')
' '
5098 WRITE (unit_nr,
'(T3,A,F36.4)')
'Alpha GW direct gap at current kpoint (eV)', &
5100 ELSE IF (do_beta)
THEN
5101 WRITE (unit_nr,
'(T3,A)')
' '
5102 WRITE (unit_nr,
'(T3,A,F37.4)')
'Beta GW direct gap at current kpoint (eV)', &
5106 IF (do_closed_shell)
THEN
5107 WRITE (unit_nr,
'(T3,A)')
' '
5108 IF (count_ev_sc_gw > 1)
THEN
5109 WRITE (unit_nr,
'(T3,A,I3,A,F39.4)')
'HOMO-LUMO gap in evGW iteration', &
5110 count_ev_sc_gw,
' (eV)', e_gap_gw*evolt
5111 ELSE IF (count_sc_gw0 > 1)
THEN
5112 WRITE (unit_nr,
'(T3,A,I3,A,F38.4)')
'HOMO-LUMO gap in evGW0 iteration', &
5113 count_sc_gw0,
' (eV)', e_gap_gw*evolt
5115 WRITE (unit_nr,
'(T3,A,F55.4)')
'G0W0 HOMO-LUMO gap (eV)', e_gap_gw*evolt
5117 ELSE IF (do_alpha)
THEN
5118 WRITE (unit_nr,
'(T3,A)')
' '
5119 WRITE (unit_nr,
'(T3,A,F51.4)')
'Alpha GW HOMO-LUMO gap (eV)', e_gap_gw*evolt
5120 ELSE IF (do_beta)
THEN
5121 WRITE (unit_nr,
'(T3,A)')
' '
5122 WRITE (unit_nr,
'(T3,A,F52.4)')
'Beta GW HOMO-LUMO gap (eV)', e_gap_gw*evolt
5127 IF (unit_nr > 0)
THEN
5128 WRITE (unit_nr, *)
' '
5129 WRITE (unit_nr,
'(T3,A)')
'------------------------------------------------------------------------------'
5132 CALL timestop(handle)
5134 END SUBROUTINE print_and_update_for_ev_sc
5145 PURE SUBROUTINE shift_unshifted_levels(Eigenval, Eigenval_last, gw_corr_lev_occ, gw_corr_lev_virt, &
5148 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: eigenval, eigenval_last
5149 INTEGER,
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo, &
5152 INTEGER :: n_level_gw, n_level_gw_ref
5153 REAL(kind=dp) :: eigen_diff
5157 IF (gw_corr_lev_occ < homo .AND. gw_corr_lev_occ > 0)
THEN
5162 DO n_level_gw = 1, gw_corr_lev_occ
5163 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
5164 eigen_diff = eigen_diff + eigenval(n_level_gw_ref) - eigenval_last(n_level_gw_ref)
5166 eigen_diff = eigen_diff/gw_corr_lev_occ
5169 DO n_level_gw = 1, homo - gw_corr_lev_occ
5170 eigenval(n_level_gw) = eigenval(n_level_gw) + eigen_diff
5176 IF (gw_corr_lev_virt < nmo - homo .AND. gw_corr_lev_virt > 0)
THEN
5180 DO n_level_gw = 1, gw_corr_lev_virt
5181 n_level_gw_ref = n_level_gw + homo
5182 eigen_diff = eigen_diff + eigenval(n_level_gw_ref) - eigenval_last(n_level_gw_ref)
5184 eigen_diff = eigen_diff/gw_corr_lev_virt
5187 DO n_level_gw = homo + gw_corr_lev_virt + 1, nmo
5188 eigenval(n_level_gw) = eigenval(n_level_gw) + eigen_diff
5193 END SUBROUTINE shift_unshifted_levels
5210 SUBROUTINE calc_mat_n(N_ij, Lambda, Sigma_c, vec_omega_fit_gw, i, j, &
5211 num_poles, num_fit_points, n_level_gw, h)
5212 REAL(kind=dp),
INTENT(OUT) :: n_ij
5213 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:), &
5214 INTENT(IN) :: lambda
5215 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: sigma_c
5216 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
5217 INTENT(IN) :: vec_omega_fit_gw
5218 INTEGER,
INTENT(IN) :: i, j, num_poles, num_fit_points, &
5220 REAL(kind=dp),
INTENT(IN) :: h
5222 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_mat_N'
5224 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: lambda_tmp
5225 INTEGER :: handle, num_var
5226 REAL(kind=dp) :: chi2, chi2_sum
5228 CALL timeset(routinen, handle)
5230 num_var = 2*num_poles + 1
5231 ALLOCATE (lambda_tmp(num_var))
5236 lambda_tmp(:) = lambda(:)
5237 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5238 num_fit_points, n_level_gw)
5241 lambda_tmp(:) = lambda(:)
5242 IF (
modulo(i, 2) == 0)
THEN
5243 lambda_tmp(i/2) = lambda_tmp(i/2) + h*z_one
5245 lambda_tmp((i + 1)/2) = lambda_tmp((i + 1)/2) + h*gaussi
5247 IF (
modulo(j, 2) == 0)
THEN
5248 lambda_tmp(j/2) = lambda_tmp(j/2) + h*z_one
5250 lambda_tmp((j + 1)/2) = lambda_tmp((j + 1)/2) + h*gaussi
5252 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5253 num_fit_points, n_level_gw)
5254 chi2_sum = chi2_sum + chi2
5256 IF (
modulo(i, 2) == 0)
THEN
5257 lambda_tmp(i/2) = lambda_tmp(i/2) - 2.0_dp*h*z_one
5259 lambda_tmp((i + 1)/2) = lambda_tmp((i + 1)/2) - 2.0_dp*h*gaussi
5261 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5262 num_fit_points, n_level_gw)
5263 chi2_sum = chi2_sum - chi2
5265 IF (
modulo(j, 2) == 0)
THEN
5266 lambda_tmp(j/2) = lambda_tmp(j/2) - 2.0_dp*h*z_one
5268 lambda_tmp((j + 1)/2) = lambda_tmp((j + 1)/2) - 2.0_dp*h*gaussi
5270 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5271 num_fit_points, n_level_gw)
5272 chi2_sum = chi2_sum + chi2
5274 IF (
modulo(i, 2) == 0)
THEN
5275 lambda_tmp(i/2) = lambda_tmp(i/2) + 2.0_dp*h*z_one
5277 lambda_tmp((i + 1)/2) = lambda_tmp((i + 1)/2) + 2.0_dp*h*gaussi
5279 CALL calc_chi2(chi2, lambda_tmp, sigma_c, vec_omega_fit_gw, num_poles, &
5280 num_fit_points, n_level_gw)
5281 chi2_sum = chi2_sum - chi2
5284 n_ij = 1.0_dp/2.0_dp*chi2_sum/(4.0_dp*h*h)
5286 DEALLOCATE (lambda_tmp)
5288 CALL timestop(handle)
5290 END SUBROUTINE calc_mat_n
5302 PURE SUBROUTINE calc_chi2(chi2, Lambda, Sigma_c, vec_omega_fit_gw, num_poles, &
5303 num_fit_points, n_level_gw)
5304 REAL(kind=dp),
INTENT(OUT) :: chi2
5305 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: lambda
5306 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: sigma_c
5307 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: vec_omega_fit_gw
5308 INTEGER,
INTENT(IN) :: num_poles, num_fit_points, n_level_gw
5310 COMPLEX(KIND=dp) :: func_val
5311 INTEGER :: iii, jjj, kkk
5314 DO kkk = 1, num_fit_points
5315 func_val = lambda(1)
5316 DO iii = 1, num_poles
5319 func_val = func_val + lambda(jjj)/(gaussi*vec_omega_fit_gw(kkk) - lambda(jjj + 1))
5321 chi2 = chi2 + (abs(sigma_c(n_level_gw, kkk) - func_val))**2
5324 END SUBROUTINE calc_chi2
5373 SUBROUTINE compute_self_energy_cubic_gw(num_integ_points, nmo, tau_tj, tj, &
5374 matrix_s, cfm_mo_coeff, Eigenval, eps_filter, &
5375 e_fermi, fm_mat_W, &
5376 gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, &
5377 count_ev_sc_GW, count_sc_GW0, &
5378 t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
5379 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
5380 mat_W, mat_MinvVMinv, mat_dm, &
5381 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_Sigma_c_gw, &
5382 do_periodic, num_points_corr, delta_corr, qs_env, para_env, para_env_RPA, &
5383 mp2_env, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
5384 first_cycle_periodic_correction, kpoints, num_fit_points, fm_mo_coeff, &
5385 do_ri_Sigma_x, vec_Sigma_x_gw, unit_nr, ispin)
5386 INTEGER,
INTENT(IN) :: num_integ_points, nmo
5387 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
5388 INTENT(IN) :: tau_tj, tj
5389 TYPE(dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_s
5390 TYPE(cp_cfm_type),
INTENT(IN) :: cfm_mo_coeff
5391 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: eigenval
5392 REAL(kind=dp),
INTENT(IN) :: eps_filter
5393 REAL(kind=dp),
INTENT(INOUT) :: e_fermi
5394 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_w
5395 INTEGER,
INTENT(IN) :: gw_corr_lev_tot, gw_corr_lev_occ, &
5396 gw_corr_lev_virt, homo, &
5397 count_ev_sc_gw, count_sc_gw0
5398 TYPE(dbt_type) :: t_3c_overl_int_ao_mo
5399 TYPE(hfx_compression_type) :: t_3c_o_mo_compressed
5400 INTEGER,
DIMENSION(:, :) :: t_3c_o_mo_ind
5401 TYPE(dbt_type) :: t_3c_overl_int_gw_ri, &
5402 t_3c_overl_int_gw_ao
5403 TYPE(dbcsr_type),
INTENT(INOUT),
TARGET :: mat_w
5404 TYPE(dbcsr_p_type) :: mat_minvvminv, mat_dm
5405 REAL(kind=dp),
DIMENSION(:, :),
INTENT(IN) :: weights_cos_tf_t_to_w, &
5406 weights_sin_tf_t_to_w
5407 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(OUT) :: vec_sigma_c_gw
5408 LOGICAL,
INTENT(IN) :: do_periodic
5409 INTEGER,
INTENT(IN) :: num_points_corr
5410 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
5411 INTENT(INOUT) :: delta_corr
5412 TYPE(qs_environment_type),
POINTER :: qs_env
5413 TYPE(mp_para_env_type),
POINTER :: para_env, para_env_rpa
5414 TYPE(mp2_type),
INTENT(INOUT) :: mp2_env
5415 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_berry_re_mo_mo, &
5416 matrix_berry_im_mo_mo
5417 LOGICAL,
INTENT(INOUT) :: first_cycle_periodic_correction
5418 TYPE(kpoint_type),
POINTER :: kpoints
5419 INTEGER,
INTENT(IN) :: num_fit_points
5420 TYPE(cp_fm_type),
INTENT(IN) :: fm_mo_coeff
5421 LOGICAL,
INTENT(IN) :: do_ri_sigma_x
5422 REAL(kind=dp),
DIMENSION(:, :),
INTENT(INOUT) :: vec_sigma_x_gw
5423 INTEGER,
INTENT(IN) :: unit_nr, ispin
5425 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_self_energy_cubic_gw'
5427 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: delta_corr_omega
5428 INTEGER :: gw_lev_end, gw_lev_start, handle, handle3, i, iblk_mo, iquad, jquad, mo_end, &
5429 mo_start, n_level_gw, n_level_gw_ref, nao, nblk_mo, unit_nr_prv
5430 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_range_mo, dist1, dist2, mo_bsizes, &
5431 mo_offsets, sizes_ao, sizes_ri
5432 INTEGER,
DIMENSION(2) :: mo_bounds, pdims_2d
5433 INTEGER,
DIMENSION(3, 1) :: index_to_cell_zero
5434 LOGICAL :: memory_info
5435 REAL(kind=dp) :: ext_scaling, omega, omega_i, omega_sign, &
5436 sign_occ_virt, t_i_clenshaw, tau, &
5437 weight_cos, weight_i, weight_sin
5438 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: vec_sigma_c_gw_cos_omega, &
5439 vec_sigma_c_gw_cos_tau, vec_sigma_c_gw_neg_tau, vec_sigma_c_gw_pos_tau, &
5440 vec_sigma_c_gw_sin_omega, vec_sigma_c_gw_sin_tau
5441 TYPE(cp_cfm_type) :: cfm_scaled_dm_occ_tau
5442 TYPE(cp_fm_type) :: fm_scaled_dm_occ_tau
5443 TYPE(dbcsr_p_type),
DIMENSION(:, :, :),
POINTER :: greens_fct
5444 TYPE(dbcsr_type),
POINTER :: mat_greens_fct_occ, mat_greens_fct_virt
5445 TYPE(dbt_pgrid_type) :: pgrid_2d
5446 TYPE(dbt_type) :: t_3c_ctr_ao, t_3c_ctr_ri, t_ao_tmp, &
5447 t_dm, t_greens_fct_occ, &
5448 t_greens_fct_virt, t_ri_tmp, &
5451 CALL timeset(routinen, handle)
5453 CALL cp_cfm_get_info(cfm_mo_coeff, nrow_global=nao)
5455 CALL decompress_tensor(t_3c_overl_int_ao_mo, t_3c_o_mo_ind, t_3c_o_mo_compressed, &
5456 mp2_env%ri_rpa_im_time%eps_compress)
5458 CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_gw_ri)
5459 CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_gw_ao, order=[2, 1, 3], move_data=.true.)
5461 memory_info = mp2_env%ri_rpa_im_time%memory_info
5462 IF (memory_info)
THEN
5463 unit_nr_prv = unit_nr
5468 mo_start = homo - gw_corr_lev_occ + 1
5469 mo_end = homo + gw_corr_lev_virt
5470 cpassert(mo_end - mo_start + 1 == gw_corr_lev_tot)
5472 vec_sigma_c_gw = z_zero
5473 ALLOCATE (vec_sigma_c_gw_pos_tau(gw_corr_lev_tot, num_integ_points))
5474 vec_sigma_c_gw_pos_tau = 0.0_dp
5475 ALLOCATE (vec_sigma_c_gw_neg_tau(gw_corr_lev_tot, num_integ_points))
5476 vec_sigma_c_gw_neg_tau = 0.0_dp
5477 ALLOCATE (vec_sigma_c_gw_cos_tau(gw_corr_lev_tot, num_integ_points))
5478 vec_sigma_c_gw_cos_tau = 0.0_dp
5479 ALLOCATE (vec_sigma_c_gw_sin_tau(gw_corr_lev_tot, num_integ_points))
5480 vec_sigma_c_gw_sin_tau = 0.0_dp
5482 ALLOCATE (vec_sigma_c_gw_cos_omega(gw_corr_lev_tot, num_integ_points))
5483 vec_sigma_c_gw_cos_omega = 0.0_dp
5484 ALLOCATE (vec_sigma_c_gw_sin_omega(gw_corr_lev_tot, num_integ_points))
5485 vec_sigma_c_gw_sin_omega = 0.0_dp
5487 ALLOCATE (delta_corr_omega(1 + homo - gw_corr_lev_occ:homo + gw_corr_lev_virt, num_integ_points))
5488 delta_corr_omega(:, :) = z_zero
5490 e_fermi = 0.5_dp*(eigenval(homo) + eigenval(homo + 1))
5492 index_to_cell_zero = 0
5493 NULLIFY (greens_fct)
5494 CALL create_propagator_matrix_set(greens_fct, 1, matrix_s(1)%matrix, index_to_cell_zero)
5495 mat_greens_fct_occ => greens_fct(propagator_sector_occupied, 1, 1)%matrix
5496 mat_greens_fct_virt => greens_fct(propagator_sector_virtual, 1, 1)%matrix
5498 nblk_mo = dbt_nblks_total(t_3c_overl_int_gw_ao, 3)
5499 ALLOCATE (mo_offsets(nblk_mo))
5500 ALLOCATE (mo_bsizes(nblk_mo))
5501 ALLOCATE (batch_range_mo(nblk_mo - 1))
5502 CALL dbt_get_info(t_3c_overl_int_gw_ao, blk_offset_3=mo_offsets, blk_size_3=mo_bsizes)
5505 CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d)
5506 ALLOCATE (sizes_ri(dbt_nblks_total(t_3c_overl_int_gw_ri, 1)))
5507 CALL dbt_get_info(t_3c_overl_int_gw_ri, blk_size_1=sizes_ri)
5509 CALL create_2c_tensor(t_w, dist1, dist2, pgrid_2d, sizes_ri, sizes_ri, name=
"(RI|RI)")
5511 DEALLOCATE (dist1, dist2)
5513 CALL dbt_create(mat_w, t_ri_tmp, name=
"(RI|RI)")
5515 CALL dbt_create(t_3c_overl_int_gw_ri, t_3c_ctr_ri)
5516 CALL dbt_create(t_3c_overl_int_gw_ao, t_3c_ctr_ao)
5518 ALLOCATE (sizes_ao(dbt_nblks_total(t_3c_overl_int_gw_ao, 1)))
5519 CALL dbt_get_info(t_3c_overl_int_gw_ao, blk_size_1=sizes_ao)
5520 CALL create_2c_tensor(t_greens_fct_occ, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5521 DEALLOCATE (dist1, dist2)
5522 CALL create_2c_tensor(t_greens_fct_virt, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5523 DEALLOCATE (dist1, dist2)
5525 DO jquad = 1, num_integ_points
5527 CALL compute_gamma_propagator(greens_fct, 1, cfm_mo_coeff, homo, eigenval, nmo, &
5528 eps_filter, e_fermi, tau_tj(jquad), para_env)
5530 CALL dbcsr_set(mat_w, 0.0_dp)
5531 CALL copy_fm_to_dbcsr(fm_mat_w(jquad), mat_w, keep_sparsity=.false.)
5533 IF (jquad == 1)
CALL dbt_create(mat_greens_fct_occ, t_ao_tmp, name=
"(AO|AO)")
5535 CALL dbt_copy_matrix_to_tensor(mat_w, t_ri_tmp)
5536 CALL dbt_copy(t_ri_tmp, t_w)
5537 CALL dbt_copy_matrix_to_tensor(mat_greens_fct_occ, t_ao_tmp)
5538 CALL dbt_copy(t_ao_tmp, t_greens_fct_occ)
5539 CALL dbt_copy_matrix_to_tensor(mat_greens_fct_virt, t_ao_tmp)
5540 CALL dbt_copy(t_ao_tmp, t_greens_fct_virt)
5542 batch_range_mo(:) = [(i, i=2, nblk_mo)]
5543 CALL dbt_batched_contract_init(t_3c_overl_int_gw_ao, batch_range_3=batch_range_mo)
5544 CALL dbt_batched_contract_init(t_3c_overl_int_gw_ri, batch_range_3=batch_range_mo)
5545 CALL dbt_batched_contract_init(t_3c_ctr_ao, batch_range_3=batch_range_mo)
5546 CALL dbt_batched_contract_init(t_3c_ctr_ri, batch_range_3=batch_range_mo)
5547 CALL dbt_batched_contract_init(t_w)
5548 CALL dbt_batched_contract_init(t_greens_fct_occ)
5549 CALL dbt_batched_contract_init(t_greens_fct_virt)
5553 DO iblk_mo = 2, nblk_mo - 1
5554 mo_bounds = [mo_offsets(iblk_mo), mo_offsets(iblk_mo) + mo_bsizes(iblk_mo) - 1]
5555 CALL contract_cubic_gw(t_3c_overl_int_gw_ao, t_3c_overl_int_gw_ri, &
5556 t_greens_fct_occ, t_w, [1.0_dp, -1.0_dp], &
5557 mo_bounds, unit_nr_prv, &
5558 t_3c_ctr_ri, t_3c_ctr_ao, calculate_ctr_ri=.true.)
5559 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)
5561 CALL contract_cubic_gw(t_3c_overl_int_gw_ao, t_3c_overl_int_gw_ri, &
5562 t_greens_fct_virt, t_w, [1.0_dp, 1.0_dp], &
5563 mo_bounds, unit_nr_prv, &
5564 t_3c_ctr_ri, t_3c_ctr_ao, calculate_ctr_ri=.false.)
5566 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)
5568 CALL dbt_batched_contract_finalize(t_3c_overl_int_gw_ao)
5569 CALL dbt_batched_contract_finalize(t_3c_overl_int_gw_ri)
5570 CALL dbt_batched_contract_finalize(t_3c_ctr_ao)
5571 CALL dbt_batched_contract_finalize(t_3c_ctr_ri)
5572 CALL dbt_batched_contract_finalize(t_w)
5573 CALL dbt_batched_contract_finalize(t_greens_fct_occ)
5574 CALL dbt_batched_contract_finalize(t_greens_fct_virt)
5576 CALL dbt_clear(t_3c_ctr_ao)
5577 CALL dbt_clear(t_3c_ctr_ri)
5579 vec_sigma_c_gw_cos_tau(:, jquad) = 0.5_dp*(vec_sigma_c_gw_pos_tau(:, jquad) + &
5580 vec_sigma_c_gw_neg_tau(:, jquad))
5582 vec_sigma_c_gw_sin_tau(:, jquad) = 0.5_dp*(vec_sigma_c_gw_pos_tau(:, jquad) - &
5583 vec_sigma_c_gw_neg_tau(:, jquad))
5586 CALL dbt_destroy(t_w)
5588 CALL dbt_destroy(t_greens_fct_occ)
5589 CALL dbt_destroy(t_greens_fct_virt)
5592 DO jquad = 1, num_fit_points
5594 DO iquad = 1, num_integ_points
5598 weight_cos = weights_cos_tf_t_to_w(jquad, iquad)*cos(omega*tau)
5599 weight_sin = weights_sin_tf_t_to_w(jquad, iquad)*sin(omega*tau)
5601 vec_sigma_c_gw_cos_omega(:, jquad) = vec_sigma_c_gw_cos_omega(:, jquad) + &
5602 weight_cos*vec_sigma_c_gw_cos_tau(:, iquad)
5604 vec_sigma_c_gw_sin_omega(:, jquad) = vec_sigma_c_gw_sin_omega(:, jquad) + &
5605 weight_sin*vec_sigma_c_gw_sin_tau(:, iquad)
5613 vec_sigma_c_gw_sin_omega(1:gw_corr_lev_occ, :) = -vec_sigma_c_gw_sin_omega(1:gw_corr_lev_occ, :)
5615 vec_sigma_c_gw(:, 1:num_fit_points, 1) = vec_sigma_c_gw_cos_omega(:, 1:num_fit_points) + &
5616 gaussi*vec_sigma_c_gw_sin_omega(:, 1:num_fit_points)
5618 CALL dbcsr_deallocate_matrix_set(greens_fct)
5620 IF (do_ri_sigma_x .AND. count_ev_sc_gw == 1 .AND. count_sc_gw0 == 1)
THEN
5622 CALL timeset(routinen//
"_RI_HFX_operation_1", handle3)
5624 CALL cp_cfm_create(cfm_scaled_dm_occ_tau, cfm_mo_coeff%matrix_struct, nrow=nao, ncol=nao)
5627 CALL parallel_gemm(transa=
"N", transb=
"C", m=nao, n=nao, k=homo, alpha=(1.0_dp, 0.0_dp), &
5628 matrix_a=cfm_mo_coeff, matrix_b=cfm_mo_coeff, beta=(0.0_dp, 0.0_dp), &
5629 matrix_c=cfm_scaled_dm_occ_tau)
5631 CALL cp_fm_create(fm_scaled_dm_occ_tau, cfm_scaled_dm_occ_tau%matrix_struct)
5632 CALL cp_cfm_to_fm(cfm_scaled_dm_occ_tau, fm_scaled_dm_occ_tau)
5634 CALL timestop(handle3)
5636 CALL timeset(routinen//
"_RI_HFX_operation_2", handle3)
5638 CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, &
5640 keep_sparsity=.false.)
5642 CALL cp_fm_release(fm_scaled_dm_occ_tau)
5643 CALL cp_cfm_release(cfm_scaled_dm_occ_tau)
5645 CALL timestop(handle3)
5647 CALL create_2c_tensor(t_dm, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5648 DEALLOCATE (dist1, dist2)
5650 CALL dbt_copy_matrix_to_tensor(mat_dm%matrix, t_ao_tmp)
5651 CALL dbt_copy(t_ao_tmp, t_dm)
5653 CALL create_2c_tensor(t_sinvvsinv, dist1, dist2, pgrid_2d, sizes_ri, sizes_ri, name=
"(RI|RI)")
5654 DEALLOCATE (dist1, dist2)
5656 CALL dbt_copy_matrix_to_tensor(mat_minvvminv%matrix, t_ri_tmp)
5657 CALL dbt_copy(t_ri_tmp, t_sinvvsinv)
5659 CALL dbt_batched_contract_init(t_3c_overl_int_gw_ao, batch_range_3=batch_range_mo)
5660 CALL dbt_batched_contract_init(t_3c_overl_int_gw_ri, batch_range_3=batch_range_mo)
5661 CALL dbt_batched_contract_init(t_3c_ctr_ri, batch_range_3=batch_range_mo)
5662 CALL dbt_batched_contract_init(t_3c_ctr_ao, batch_range_3=batch_range_mo)
5663 CALL dbt_batched_contract_init(t_dm)
5664 CALL dbt_batched_contract_init(t_sinvvsinv)
5666 DO iblk_mo = 2, nblk_mo - 1
5667 mo_bounds = [mo_offsets(iblk_mo), mo_offsets(iblk_mo) + mo_bsizes(iblk_mo) - 1]
5669 CALL contract_cubic_gw(t_3c_overl_int_gw_ao, t_3c_overl_int_gw_ri, &
5670 t_dm, t_sinvvsinv, [1.0_dp, -1.0_dp], &
5671 mo_bounds, unit_nr_prv, &
5672 t_3c_ctr_ri, t_3c_ctr_ao, calculate_ctr_ri=.true.)
5674 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)
5676 CALL dbt_batched_contract_finalize(t_3c_overl_int_gw_ao)
5677 CALL dbt_batched_contract_finalize(t_3c_overl_int_gw_ri)
5678 CALL dbt_batched_contract_finalize(t_dm)
5679 CALL dbt_batched_contract_finalize(t_sinvvsinv)
5680 CALL dbt_batched_contract_finalize(t_3c_ctr_ri)
5681 CALL dbt_batched_contract_finalize(t_3c_ctr_ao)
5683 CALL dbt_destroy(t_dm)
5684 CALL dbt_destroy(t_sinvvsinv)
5686 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, 1) = &
5687 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, 1) + &
5688 vec_sigma_x_gw(:, 1)
5692 CALL dbt_pgrid_destroy(pgrid_2d)
5694 CALL dbt_destroy(t_3c_ctr_ri)
5695 CALL dbt_destroy(t_3c_ctr_ao)
5696 CALL dbt_destroy(t_ao_tmp)
5697 CALL dbt_destroy(t_ri_tmp)
5700 IF (do_periodic)
THEN
5702 ext_scaling = 0.2_dp
5705 DO iquad = 1, num_points_corr
5708 t_i_clenshaw = iquad*pi/(2.0_dp*num_points_corr)
5709 omega_i = ext_scaling/tan(t_i_clenshaw)
5711 IF (iquad < num_points_corr)
THEN
5712 weight_i = ext_scaling*pi/(num_points_corr*sin(t_i_clenshaw)**2)
5714 weight_i = ext_scaling*pi/(2.0_dp*num_points_corr*sin(t_i_clenshaw)**2)
5717 CALL calc_periodic_correction(delta_corr, qs_env, para_env, para_env_rpa, &
5718 mp2_env%ri_g0w0%kp_grid, homo, nmo, gw_corr_lev_occ, &
5719 gw_corr_lev_virt, omega_i, fm_mo_coeff, eigenval, &
5720 matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
5721 first_cycle_periodic_correction, kpoints, &
5722 mp2_env%ri_g0w0%do_mo_coeff_gamma, &
5723 mp2_env%ri_g0w0%num_kp_grids, mp2_env%ri_g0w0%eps_kpoint, &
5724 mp2_env%ri_g0w0%do_extra_kpoints, &
5725 mp2_env%ri_g0w0%do_aux_bas_gw, mp2_env%ri_g0w0%frac_aux_mos)
5727 DO n_level_gw = 1, gw_corr_lev_tot
5729 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
5731 IF (n_level_gw <= gw_corr_lev_occ)
THEN
5732 sign_occ_virt = -1.0_dp
5734 sign_occ_virt = 1.0_dp
5737 DO jquad = 1, num_integ_points
5739 omega_sign = tj(jquad)*sign_occ_virt
5741 delta_corr_omega(n_level_gw_ref, jquad) = &
5742 delta_corr_omega(n_level_gw_ref, jquad) - &
5743 0.5_dp/pi*weight_i/2.0_dp*delta_corr(n_level_gw_ref)* &
5744 (1.0_dp/(gaussi*(omega_i + omega_sign) + e_fermi - eigenval(n_level_gw_ref)) + &
5745 1.0_dp/(gaussi*(-omega_i + omega_sign) + e_fermi - eigenval(n_level_gw_ref)))
5753 gw_lev_start = 1 + homo - gw_corr_lev_occ
5754 gw_lev_end = homo + gw_corr_lev_virt
5757 vec_sigma_c_gw(1:gw_corr_lev_tot, :, 1) = vec_sigma_c_gw(1:gw_corr_lev_tot, :, 1) + &
5758 delta_corr_omega(gw_lev_start:gw_lev_end, 1:num_fit_points)
5762 DEALLOCATE (vec_sigma_c_gw_pos_tau)
5763 DEALLOCATE (vec_sigma_c_gw_neg_tau)
5764 DEALLOCATE (vec_sigma_c_gw_cos_tau)
5765 DEALLOCATE (vec_sigma_c_gw_sin_tau)
5766 DEALLOCATE (vec_sigma_c_gw_cos_omega)
5767 DEALLOCATE (vec_sigma_c_gw_sin_omega)
5768 DEALLOCATE (delta_corr_omega)
5770 CALL timestop(handle)
5772 END SUBROUTINE compute_self_energy_cubic_gw
5811 SUBROUTINE compute_self_energy_cubic_gw_kpoints(num_integ_points, tau_tj, tj, &
5812 matrix_s, Eigenval, e_fermi, fm_mat_W, &
5813 gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, &
5814 count_ev_sc_GW, count_sc_GW0, &
5815 t_3c_O, t_3c_M, t_3c_O_compressed, t_3c_O_ind, &
5816 mat_W, mat_MinvVMinv, &
5817 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, vec_Sigma_c_gw, &
5819 mp2_env, num_fit_points, fm_mo_coeff, &
5820 do_ri_Sigma_x, vec_Sigma_x_gw, unit_nr, nspins, &
5821 starts_array_mc, ends_array_mc, eps_filter)
5823 INTEGER,
INTENT(IN) :: num_integ_points
5824 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
5825 INTENT(IN) :: tau_tj, tj
5826 TYPE(dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_s
5827 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(IN) :: eigenval
5828 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: e_fermi
5829 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_w
5830 INTEGER,
INTENT(IN) :: gw_corr_lev_tot
5831 INTEGER,
DIMENSION(:),
INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
5832 INTEGER,
INTENT(IN) :: count_ev_sc_gw, count_sc_gw0
5833 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_o
5834 TYPE(dbt_type) :: t_3c_m
5835 TYPE(hfx_compression_type),
ALLOCATABLE, &
5836 DIMENSION(:, :, :) :: t_3c_o_compressed
5837 TYPE(block_ind_type),
ALLOCATABLE, &
5838 DIMENSION(:, :, :),
INTENT(INOUT) :: t_3c_o_ind
5839 TYPE(dbcsr_type),
INTENT(INOUT),
TARGET :: mat_w
5840 TYPE(dbcsr_p_type) :: mat_minvvminv
5841 REAL(kind=dp),
DIMENSION(:, :),
INTENT(IN) :: weights_cos_tf_t_to_w, &
5842 weights_sin_tf_t_to_w
5843 COMPLEX(KIND=dp),
DIMENSION(:, :, :, :), &
5844 INTENT(OUT) :: vec_sigma_c_gw
5845 TYPE(qs_environment_type),
POINTER :: qs_env
5846 TYPE(mp_para_env_type),
POINTER :: para_env
5847 TYPE(mp2_type),
INTENT(INOUT) :: mp2_env
5848 INTEGER,
INTENT(IN) :: num_fit_points
5849 TYPE(cp_fm_type),
INTENT(IN) :: fm_mo_coeff
5850 LOGICAL,
INTENT(IN) :: do_ri_sigma_x
5851 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: vec_sigma_x_gw
5852 INTEGER,
INTENT(IN) :: unit_nr, nspins
5853 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc
5854 REAL(kind=dp),
INTENT(IN) :: eps_filter
5856 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_self_energy_cubic_gw_kpoints'
5858 INTEGER :: cut_memory, handle, handle2, i_mem, &
5859 iquad, ispin, j_mem, jquad, &
5860 nkp_self_energy, num_points, &
5862 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist1, dist2, sizes_ao, sizes_ri
5863 INTEGER,
DIMENSION(2) :: mo_end, mo_start, pdims_2d
5864 INTEGER,
DIMENSION(2, 1) :: bounds_ri_i
5865 INTEGER,
DIMENSION(2, 2) :: bounds_ao_ao_j
5866 INTEGER,
DIMENSION(3) :: dims_3c
5867 LOGICAL :: memory_info
5868 REAL(kind=dp) :: omega, t1, t2, tau, weight_cos, &
5870 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: vec_sigma_c_gw_cos_omega, &
5871 vec_sigma_c_gw_cos_tau, vec_sigma_c_gw_neg_tau, vec_sigma_c_gw_pos_tau, &
5872 vec_sigma_c_gw_sin_omega, vec_sigma_c_gw_sin_tau
5873 TYPE(dbcsr_p_type),
DIMENSION(:, :, :),
POINTER :: propagator
5874 TYPE(dbcsr_type),
TARGET :: mat_greens_fct_occ, mat_greens_fct_virt, mat_mo_coeff, &
5875 mat_self_energy_ao_ao_neg_tau, mat_self_energy_ao_ao_pos_tau
5876 TYPE(dbt_pgrid_type) :: pgrid_2d
5877 TYPE(dbt_type) :: t_3c_m_w_tmp, t_3c_o_all, t_3c_o_w, &
5878 t_ao_tmp, t_greens_fct_occ, &
5879 t_greens_fct_virt, t_ri_tmp, t_w
5881 CALL timeset(routinen, handle)
5883 NULLIFY (propagator)
5885 memory_info = mp2_env%ri_rpa_im_time%memory_info
5886 IF (memory_info)
THEN
5887 unit_nr_prv = unit_nr
5892 cut_memory = mp2_env%ri_rpa_im_time%cut_memory
5894 DO ispin = 1, nspins
5895 mo_start(ispin) = homo(ispin) - gw_corr_lev_occ(ispin) + 1
5896 mo_end(ispin) = homo(ispin) + gw_corr_lev_virt(ispin)
5897 cpassert(mo_end(ispin) - mo_start(ispin) + 1 == gw_corr_lev_tot)
5900 nkp_self_energy = mp2_env%ri_g0w0%nkp_self_energy
5902 vec_sigma_c_gw = z_zero
5903 ALLOCATE (vec_sigma_c_gw_pos_tau(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5904 vec_sigma_c_gw_pos_tau = 0.0_dp
5905 ALLOCATE (vec_sigma_c_gw_neg_tau(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5906 vec_sigma_c_gw_neg_tau = 0.0_dp
5907 ALLOCATE (vec_sigma_c_gw_cos_tau(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5908 vec_sigma_c_gw_cos_tau = 0.0_dp
5909 ALLOCATE (vec_sigma_c_gw_sin_tau(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5910 vec_sigma_c_gw_sin_tau = 0.0_dp
5912 ALLOCATE (vec_sigma_c_gw_cos_omega(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5913 vec_sigma_c_gw_cos_omega = 0.0_dp
5914 ALLOCATE (vec_sigma_c_gw_sin_omega(gw_corr_lev_tot, num_integ_points, nkp_self_energy, nspins))
5915 vec_sigma_c_gw_sin_omega = 0.0_dp
5917 CALL dbcsr_create(matrix=mat_greens_fct_occ, &
5918 template=matrix_s(1)%matrix, &
5919 matrix_type=dbcsr_type_no_symmetry)
5921 CALL dbcsr_create(matrix=mat_greens_fct_virt, &
5922 template=matrix_s(1)%matrix, &
5923 matrix_type=dbcsr_type_no_symmetry)
5925 CALL dbcsr_create(matrix=mat_self_energy_ao_ao_neg_tau, &
5926 template=matrix_s(1)%matrix, &
5927 matrix_type=dbcsr_type_no_symmetry)
5929 CALL dbcsr_create(matrix=mat_self_energy_ao_ao_pos_tau, &
5930 template=matrix_s(1)%matrix, &
5931 matrix_type=dbcsr_type_no_symmetry)
5933 CALL dbcsr_create(matrix=mat_mo_coeff, &
5934 template=matrix_s(1)%matrix, &
5935 matrix_type=dbcsr_type_no_symmetry)
5937 CALL copy_fm_to_dbcsr(fm_mo_coeff, mat_mo_coeff, keep_sparsity=.false.)
5939 DO ispin = 1, nspins
5940 e_fermi(ispin) = 0.5_dp*(maxval(eigenval(homo, :, ispin)) + minval(eigenval(homo + 1, :, ispin)))
5944 CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d)
5945 ALLOCATE (sizes_ri(dbt_nblks_total(t_3c_o(1, 1), 1)))
5946 CALL dbt_get_info(t_3c_o(1, 1), blk_size_1=sizes_ri)
5948 CALL create_2c_tensor(t_w, dist1, dist2, pgrid_2d, sizes_ri, sizes_ri, name=
"(RI|RI)")
5949 DEALLOCATE (dist1, dist2)
5951 CALL dbt_create(mat_w, t_ri_tmp, name=
"(RI|RI)")
5953 ALLOCATE (sizes_ao(dbt_nblks_total(t_3c_o(1, 1), 2)))
5954 CALL dbt_get_info(t_3c_o(1, 1), blk_size_2=sizes_ao)
5955 CALL create_2c_tensor(t_greens_fct_occ, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5957 DEALLOCATE (dist1, dist2)
5958 CALL create_2c_tensor(t_greens_fct_virt, dist1, dist2, pgrid_2d, sizes_ao, sizes_ao, name=
"(AO|AO)")
5959 DEALLOCATE (dist1, dist2)
5961 CALL dbt_get_info(t_3c_m, nfull_total=dims_3c)
5963 CALL dbt_create(t_3c_o(1, 1), t_3c_o_all, name=
"O (RI AO | AO)")
5966 DO i_mem = 1, cut_memory
5967 CALL decompress_tensor(t_3c_o(1, 1), &
5968 t_3c_o_ind(1, 1, i_mem)%ind, &
5969 t_3c_o_compressed(1, 1, i_mem), &
5970 mp2_env%ri_rpa_im_time%eps_compress)
5971 CALL dbt_copy(t_3c_o(1, 1), t_3c_o_all, summation=.true., move_data=.true.)
5974 CALL dbt_create(t_3c_m, t_3c_m_w_tmp, name=
"M W (RI | AO AO)")
5975 CALL dbt_create(t_3c_o(1, 1), t_3c_o_w, name=
"M W (RI AO | AO)")
5977 CALL dbt_create(mat_greens_fct_occ, t_ao_tmp, name=
"(AO|AO)")
5979 IF (count_ev_sc_gw == 1 .AND. count_sc_gw0 == 1 .AND. do_ri_sigma_x)
THEN
5980 num_points = num_integ_points + 1
5982 num_points = num_integ_points
5985 DO jquad = 1, num_points
5989 IF (jquad <= num_integ_points)
THEN
5992 IF (unit_nr > 0)
WRITE (unit_nr,
'(/T3,A,1X,I3)') &
5993 'GW_INFO| Computing self-energy time point', jquad
5997 IF (unit_nr > 0)
WRITE (unit_nr,
'(/T3,A,1X,I3)') &
5998 'GW_INFO| Computing exchange self-energy'
6001 IF (jquad <= num_integ_points)
THEN
6002 CALL dbcsr_set(mat_w, 0.0_dp)
6003 CALL copy_fm_to_dbcsr(fm_mat_w(jquad), mat_w, keep_sparsity=.false.)
6004 CALL dbt_copy_matrix_to_tensor(mat_w, t_ri_tmp)
6006 CALL dbt_copy_matrix_to_tensor(mat_minvvminv%matrix, t_ri_tmp)
6009 CALL dbt_copy(t_ri_tmp, t_w)
6011 DO ispin = 1, nspins
6013 CALL compute_periodic_dm(propagator, qs_env, &
6014 ispin, num_points, jquad, e_fermi(ispin), tau, &
6015 sector=propagator_sector_occupied)
6017 CALL compute_periodic_dm(propagator, qs_env, &
6018 ispin, num_points, jquad, e_fermi(ispin), tau, &
6019 sector=propagator_sector_virtual)
6021 CALL dbcsr_set(mat_greens_fct_occ, 0.0_dp)
6022 CALL dbcsr_copy(mat_greens_fct_occ, &
6023 propagator(propagator_sector_occupied, jquad, 1)%matrix)
6025 CALL dbcsr_set(mat_greens_fct_virt, 0.0_dp)
6026 CALL dbcsr_copy(mat_greens_fct_virt, &
6027 propagator(propagator_sector_virtual, jquad, 1)%matrix)
6029 CALL dbt_copy_matrix_to_tensor(mat_greens_fct_occ, t_ao_tmp)
6030 CALL dbt_copy(t_ao_tmp, t_greens_fct_occ)
6032 CALL dbt_copy_matrix_to_tensor(mat_greens_fct_virt, t_ao_tmp)
6033 CALL dbt_copy(t_ao_tmp, t_greens_fct_virt)
6035 CALL dbcsr_set(mat_self_energy_ao_ao_neg_tau, 0.0_dp)
6036 CALL dbcsr_set(mat_self_energy_ao_ao_pos_tau, 0.0_dp)
6038 CALL dbt_copy(t_3c_o_all, t_3c_m)
6040 CALL dbt_batched_contract_init(t_3c_o_w)
6044 DO i_mem = 1, cut_memory
6050 bounds_ri_i(:, 1) = [qs_env%mp2_env%ri_rpa_im_time%starts_array_mc_RI(i_mem), &
6051 qs_env%mp2_env%ri_rpa_im_time%ends_array_mc_RI(i_mem)]
6053 DO j_mem = 1, cut_memory
6055 bounds_ao_ao_j(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
6056 bounds_ao_ao_j(:, 2) = [1, dims_3c(3)]
6058 CALL timeset(
"tensor_operation_3c_W", handle2)
6060 CALL dbt_contract(1.0_dp, t_w, t_3c_m, 0.0_dp, &
6062 contract_1=[2], notcontract_1=[1], &
6063 contract_2=[1], notcontract_2=[2, 3], &
6064 map_1=[1], map_2=[2, 3], &
6065 bounds_2=bounds_ri_i, &
6066 bounds_3=bounds_ao_ao_j, &
6067 filter_eps=eps_filter, &
6068 unit_nr=unit_nr_prv)
6070 CALL dbt_copy(t_3c_m_w_tmp, t_3c_o_w, order=[1, 2, 3], move_data=.true.)
6072 CALL timestop(handle2)
6074 CALL contract_to_self_energy(t_3c_o_all, t_greens_fct_occ, t_3c_o_w, &
6075 mat_self_energy_ao_ao_neg_tau, &
6076 bounds_ao_ao_j, bounds_ri_i, unit_nr_prv, &
6077 eps_filter, do_occ=.true., do_virt=.false.)
6079 CALL contract_to_self_energy(t_3c_o_all, t_greens_fct_virt, t_3c_o_w, &
6080 mat_self_energy_ao_ao_pos_tau, &
6081 bounds_ao_ao_j, bounds_ri_i, unit_nr_prv, &
6082 eps_filter, do_occ=.false., do_virt=.true.)
6092 CALL dbt_batched_contract_finalize(t_3c_o_w)
6096 IF (jquad <= num_integ_points)
THEN
6098 CALL trafo_to_mo_and_kpoints(qs_env, mat_self_energy_ao_ao_neg_tau, vec_sigma_c_gw_neg_tau(:, jquad, :, ispin), &
6099 homo(ispin), gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), ispin)
6101 CALL trafo_to_mo_and_kpoints(qs_env, mat_self_energy_ao_ao_pos_tau, vec_sigma_c_gw_pos_tau(:, jquad, :, ispin), &
6102 homo(ispin), gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), ispin)
6104 vec_sigma_c_gw_cos_tau(:, jquad, :, ispin) = 0.5_dp*(vec_sigma_c_gw_pos_tau(:, jquad, :, ispin) + &
6105 vec_sigma_c_gw_neg_tau(:, jquad, :, ispin))
6107 vec_sigma_c_gw_sin_tau(:, jquad, :, ispin) = 0.5_dp*(vec_sigma_c_gw_pos_tau(:, jquad, :, ispin) - &
6108 vec_sigma_c_gw_neg_tau(:, jquad, :, ispin))
6112 vec_sigma_x_gw(mo_start(ispin):mo_end(ispin), :, ispin), &
6113 homo(ispin), gw_corr_lev_occ(ispin), gw_corr_lev_virt(ispin), ispin)
6121 IF (unit_nr > 0)
WRITE (unit_nr,
'(T6,A,T56,F25.1)')
'Execution time (s):', t2 - t1
6125 IF (count_ev_sc_gw == 1 .AND. count_sc_gw0 == 1)
THEN
6129 IF (do_ri_sigma_x)
THEN
6130 DO ispin = 1, nspins
6131 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, :) = mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, :) + &
6132 vec_sigma_x_gw(:, :, ispin)
6139 DO jquad = 1, num_fit_points
6141 DO iquad = 1, num_integ_points
6145 weight_cos = weights_cos_tf_t_to_w(jquad, iquad)*cos(omega*tau)
6146 weight_sin = weights_sin_tf_t_to_w(jquad, iquad)*sin(omega*tau)
6148 vec_sigma_c_gw_cos_omega(:, jquad, :, :) = vec_sigma_c_gw_cos_omega(:, jquad, :, :) + &
6149 weight_cos*vec_sigma_c_gw_cos_tau(:, iquad, :, :)
6151 vec_sigma_c_gw_sin_omega(:, jquad, :, :) = vec_sigma_c_gw_sin_omega(:, jquad, :, :) + &
6152 weight_sin*vec_sigma_c_gw_sin_tau(:, iquad, :, :)
6160 DO ispin = 1, nspins
6161 vec_sigma_c_gw_sin_omega(1:gw_corr_lev_occ(ispin), :, :, ispin) = &
6162 -vec_sigma_c_gw_sin_omega(1:gw_corr_lev_occ(ispin), :, :, ispin)
6165 vec_sigma_c_gw(:, 1:num_fit_points, :, :) = vec_sigma_c_gw_cos_omega(:, 1:num_fit_points, :, :) + &
6166 gaussi*vec_sigma_c_gw_sin_omega(:, 1:num_fit_points, :, :)
6168 CALL dbt_pgrid_destroy(pgrid_2d)
6170 CALL dbcsr_release(mat_greens_fct_occ)
6171 CALL dbcsr_release(mat_greens_fct_virt)
6172 CALL dbcsr_release(mat_self_energy_ao_ao_neg_tau)
6173 CALL dbcsr_release(mat_self_energy_ao_ao_pos_tau)
6174 CALL dbcsr_release(mat_mo_coeff)
6176 CALL dbcsr_deallocate_matrix_set(propagator)
6178 CALL dbt_destroy(t_w)
6179 CALL dbt_destroy(t_ri_tmp)
6180 CALL dbt_destroy(t_greens_fct_occ)
6181 CALL dbt_destroy(t_greens_fct_virt)
6182 CALL dbt_destroy(t_ao_tmp)
6183 CALL dbt_destroy(t_3c_o_all)
6184 CALL dbt_destroy(t_3c_m_w_tmp)
6185 CALL dbt_destroy(t_3c_o_w)
6187 DEALLOCATE (vec_sigma_c_gw_pos_tau)
6188 DEALLOCATE (vec_sigma_c_gw_neg_tau)
6189 DEALLOCATE (vec_sigma_c_gw_cos_tau)
6190 DEALLOCATE (vec_sigma_c_gw_sin_tau)
6191 DEALLOCATE (vec_sigma_c_gw_cos_omega)
6192 DEALLOCATE (vec_sigma_c_gw_sin_omega)
6194 CALL timestop(handle)
6196 END SUBROUTINE compute_self_energy_cubic_gw_kpoints
6203 TYPE(qs_environment_type),
POINTER :: qs_env
6205 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_minus_vxc_kpoints'
6207 INTEGER :: handle, ikp, ispin, nkp_self_energy, &
6209 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: diag_sigma_x_minus_vxc_mo_mo
6210 TYPE(cp_cfm_type) :: cfm_mo_coeff, ks_mat_ao_ao, &
6211 ks_mat_no_xc_ao_ao, vxc_ao_ao, &
6212 vxc_ao_mo, vxc_mo_mo
6213 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct
6214 TYPE(cp_fm_type) :: fm_dummy, fm_sigma_x_minus_vxc_mo_mo, &
6215 fm_tmp_im, fm_tmp_re
6216 TYPE(dft_control_type),
POINTER :: dft_control
6217 TYPE(kpoint_type),
POINTER :: kpoints_sigma, kpoints_sigma_no_xc
6218 TYPE(mp_para_env_type),
POINTER :: para_env
6220 CALL timeset(routinen, handle)
6222 CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control)
6224 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
6226 kpoints_sigma_no_xc => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc
6228 nkp_self_energy = kpoints_sigma%nkp
6230 nspins = dft_control%nspins
6232 matrix_struct => kpoints_sigma%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct
6234 CALL cp_cfm_create(ks_mat_ao_ao, matrix_struct)
6235 CALL cp_cfm_create(ks_mat_no_xc_ao_ao, matrix_struct)
6236 CALL cp_cfm_create(vxc_ao_ao, matrix_struct)
6237 CALL cp_cfm_create(vxc_ao_mo, matrix_struct)
6238 CALL cp_cfm_create(vxc_mo_mo, matrix_struct)
6239 CALL cp_cfm_create(cfm_mo_coeff, matrix_struct)
6240 CALL cp_fm_create(fm_sigma_x_minus_vxc_mo_mo, matrix_struct)
6241 CALL cp_fm_create(fm_tmp_re, matrix_struct)
6242 CALL cp_fm_create(fm_tmp_im, matrix_struct)
6244 CALL cp_cfm_get_info(cfm_mo_coeff, nrow_global=nmo)
6245 ALLOCATE (diag_sigma_x_minus_vxc_mo_mo(nmo))
6247 DEALLOCATE (qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw)
6249 ALLOCATE (qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(nmo, 2, nkp_self_energy))
6251 DO ikp = 1, nkp_self_energy
6253 DO ispin = 1, nspins
6255 associate(mos => kpoints_sigma%kp_env(ikp)%kpoint_env%mos)
6256 IF (
ASSOCIATED(mos(1, ispin)%mo_coeff))
THEN
6257 CALL cp_fm_copy_general(mos(1, ispin)%mo_coeff, fm_tmp_re, para_env)
6259 CALL cp_fm_copy_general(fm_dummy, fm_tmp_re, para_env)
6261 IF (
ASSOCIATED(mos(2, ispin)%mo_coeff))
THEN
6262 CALL cp_fm_copy_general(mos(2, ispin)%mo_coeff, fm_tmp_im, para_env)
6264 CALL cp_fm_copy_general(fm_dummy, fm_tmp_im, para_env)
6268 CALL cp_fm_to_cfm(fm_tmp_re, fm_tmp_im, cfm_mo_coeff)
6270 CALL cp_fm_to_cfm(kpoints_sigma%kp_env(ikp)%kpoint_env%wmat(1, ispin), &
6271 kpoints_sigma%kp_env(ikp)%kpoint_env%wmat(2, ispin), ks_mat_ao_ao)
6272 associate(wmat => kpoints_sigma_no_xc%kp_env(ikp)%kpoint_env%wmat)
6273 IF (
ASSOCIATED(wmat(1, ispin)%matrix_struct))
THEN
6274 CALL cp_fm_copy_general(wmat(1, ispin), fm_tmp_re, para_env)
6276 CALL cp_fm_copy_general(fm_dummy, fm_tmp_re, para_env)
6278 IF (
ASSOCIATED(wmat(2, ispin)%matrix_struct))
THEN
6279 CALL cp_fm_copy_general(wmat(2, ispin), fm_tmp_im, para_env)
6281 CALL cp_fm_copy_general(fm_dummy, fm_tmp_im, para_env)
6285 CALL cp_fm_to_cfm(fm_tmp_re, fm_tmp_im, vxc_ao_ao)
6287 CALL parallel_gemm(
'N',
'N', nmo, nmo, nmo, z_one, vxc_ao_ao, cfm_mo_coeff, z_zero, vxc_ao_mo)
6288 CALL parallel_gemm(
'C',
'N', nmo, nmo, nmo, z_one, cfm_mo_coeff, vxc_ao_mo, z_zero, vxc_mo_mo)
6290 CALL cp_cfm_to_fm(vxc_mo_mo, fm_sigma_x_minus_vxc_mo_mo)
6292 CALL cp_fm_get_diag(fm_sigma_x_minus_vxc_mo_mo, diag_sigma_x_minus_vxc_mo_mo)
6294 qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, ikp) = diag_sigma_x_minus_vxc_mo_mo(:)
6300 CALL cp_cfm_release(ks_mat_ao_ao)
6301 CALL cp_cfm_release(ks_mat_no_xc_ao_ao)
6302 CALL cp_cfm_release(vxc_ao_ao)
6303 CALL cp_cfm_release(vxc_ao_mo)
6304 CALL cp_cfm_release(vxc_mo_mo)
6305 CALL cp_cfm_release(cfm_mo_coeff)
6306 CALL cp_fm_release(fm_sigma_x_minus_vxc_mo_mo)
6307 CALL cp_fm_release(fm_tmp_re)
6308 CALL cp_fm_release(fm_tmp_im)
6310 DEALLOCATE (diag_sigma_x_minus_vxc_mo_mo)
6312 CALL timestop(handle)
6327 homo, gw_corr_lev_occ, gw_corr_lev_virt, ispin)
6328 TYPE(qs_environment_type),
POINTER :: qs_env
6329 TYPE(dbcsr_type),
TARGET :: mat_self_energy_ao_ao
6330 REAL(kind=dp),
DIMENSION(:, :) :: vec_sigma
6331 INTEGER :: homo, gw_corr_lev_occ, gw_corr_lev_virt, &
6334 CHARACTER(LEN=*),
PARAMETER :: routinen =
'trafo_to_mo_and_kpoints'
6336 INTEGER :: handle, ikp, nkp_self_energy, nmo, &
6337 periodic(3), size_real_space
6338 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: diag_self_energy
6339 TYPE(cell_type),
POINTER :: cell
6340 TYPE(cp_cfm_type) :: cfm_mo_coeff, cfm_self_energy_ao_ao, &
6341 cfm_self_energy_ao_mo, &
6342 cfm_self_energy_mo_mo
6343 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct
6344 TYPE(cp_fm_type) :: fm_self_energy_mo_mo
6345 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mat_self_energy_ao_ao_kp_im, &
6346 mat_self_energy_ao_ao_kp_re, mat_self_energy_ao_ao_real_space
6347 TYPE(kpoint_type),
POINTER :: kpoints_sigma
6348 TYPE(mp_para_env_type),
POINTER :: para_env
6350 CALL timeset(routinen, handle)
6352 CALL get_qs_env(qs_env, cell=cell, para_env=para_env)
6353 CALL get_cell(cell=cell, periodic=periodic)
6355 size_real_space = 3**(periodic(1) + periodic(2) + periodic(3))
6357 CALL alloc_mat_set(mat_self_energy_ao_ao_real_space, size_real_space, mat_self_energy_ao_ao)
6359 CALL dbcsr_copy(mat_self_energy_ao_ao_real_space(1)%matrix, mat_self_energy_ao_ao)
6361 kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
6363 CALL get_mat_cell_t_from_mat_gamma(mat_self_energy_ao_ao_real_space, qs_env, kpoints_sigma, 0, 0)
6365 nkp_self_energy = kpoints_sigma%nkp
6367 CALL alloc_mat_set(mat_self_energy_ao_ao_kp_re, nkp_self_energy, mat_self_energy_ao_ao)
6368 CALL alloc_mat_set(mat_self_energy_ao_ao_kp_im, nkp_self_energy, mat_self_energy_ao_ao)
6370 CALL real_space_to_kpoint_transform_rpa(mat_self_energy_ao_ao_kp_re, mat_self_energy_ao_ao_kp_im, &
6371 mat_self_energy_ao_ao_real_space, kpoints_sigma, 1.0e-50_dp)
6373 CALL dbcsr_get_info(mat_self_energy_ao_ao, nfullrows_total=nmo)
6374 ALLOCATE (diag_self_energy(nmo))
6376 matrix_struct => kpoints_sigma%kp_env(1)%kpoint_env%mos(1, 1)%mo_coeff%matrix_struct
6378 CALL cp_cfm_create(cfm_self_energy_ao_ao, matrix_struct)
6379 CALL cp_cfm_create(cfm_self_energy_ao_mo, matrix_struct)
6380 CALL cp_cfm_create(cfm_self_energy_mo_mo, matrix_struct)
6381 CALL cp_cfm_set_all(cfm_self_energy_ao_ao, z_zero)
6382 CALL cp_cfm_set_all(cfm_self_energy_ao_mo, z_zero)
6383 CALL cp_cfm_set_all(cfm_self_energy_mo_mo, z_zero)
6385 CALL cp_fm_create(fm_self_energy_mo_mo, matrix_struct)
6386 CALL cp_cfm_create(cfm_mo_coeff, matrix_struct)
6388 DO ikp = 1, nkp_self_energy
6390 CALL dbcsr_to_cfm(mat_self_energy_ao_ao_kp_re(ikp)%matrix, &
6391 mat_self_energy_ao_ao_kp_im(ikp)%matrix, cfm_self_energy_ao_ao)
6393 CALL cp_fm_to_cfm(kpoints_sigma%kp_env(ikp)%kpoint_env%mos(1, ispin)%mo_coeff, &
6394 kpoints_sigma%kp_env(ikp)%kpoint_env%mos(2, ispin)%mo_coeff, cfm_mo_coeff)
6396 CALL parallel_gemm(
'N',
'N', nmo, nmo, nmo, z_one, cfm_self_energy_ao_ao, cfm_mo_coeff, &
6397 z_zero, cfm_self_energy_ao_mo)
6399 CALL parallel_gemm(
'C',
'N', nmo, nmo, nmo, z_one, cfm_mo_coeff, cfm_self_energy_ao_mo, &
6400 z_zero, cfm_self_energy_mo_mo)
6402 CALL cp_cfm_to_fm(cfm_self_energy_mo_mo, fm_self_energy_mo_mo)
6404 CALL cp_fm_get_diag(fm_self_energy_mo_mo, diag_self_energy)
6406 vec_sigma(:, ikp) = diag_self_energy(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt)
6410 CALL dbcsr_deallocate_matrix_set(mat_self_energy_ao_ao_real_space)
6411 CALL dbcsr_deallocate_matrix_set(mat_self_energy_ao_ao_kp_re)
6412 CALL dbcsr_deallocate_matrix_set(mat_self_energy_ao_ao_kp_im)
6414 CALL cp_cfm_release(cfm_self_energy_ao_ao)
6415 CALL cp_cfm_release(cfm_self_energy_ao_mo)
6416 CALL cp_cfm_release(cfm_self_energy_mo_mo)
6417 CALL cp_cfm_release(cfm_mo_coeff)
6418 CALL cp_fm_release(fm_self_energy_mo_mo)
6420 DEALLOCATE (diag_self_energy)
6422 CALL timestop(handle)
6432 SUBROUTINE dbcsr_to_cfm(dbcsr_re, dbcsr_im, cfm_mat)
6434 TYPE(dbcsr_type),
POINTER :: dbcsr_re, dbcsr_im
6435 TYPE(cp_cfm_type),
INTENT(IN) :: cfm_mat
6437 CHARACTER(LEN=*),
PARAMETER :: routinen =
'dbcsr_to_cfm'
6440 TYPE(cp_fm_type) :: fm_mat_im, fm_mat_re
6442 CALL timeset(routinen, handle)
6444 CALL cp_fm_create(fm_mat_re, cfm_mat%matrix_struct)
6445 CALL cp_fm_create(fm_mat_im, cfm_mat%matrix_struct)
6447 CALL copy_dbcsr_to_fm(dbcsr_re, fm_mat_re)
6448 CALL copy_dbcsr_to_fm(dbcsr_im, fm_mat_im)
6450 CALL cp_fm_to_cfm(fm_mat_re, fm_mat_im, cfm_mat)
6452 CALL cp_fm_release(fm_mat_re)
6453 CALL cp_fm_release(fm_mat_im)
6455 CALL timestop(handle)
6457 END SUBROUTINE dbcsr_to_cfm
6466 SUBROUTINE alloc_mat_set(mat_set, mat_size, template, explicitly_no_symmetry)
6467 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mat_set
6468 INTEGER,
INTENT(IN) :: mat_size
6469 TYPE(dbcsr_type),
TARGET :: template
6470 LOGICAL,
OPTIONAL :: explicitly_no_symmetry
6472 CHARACTER(LEN=*),
PARAMETER :: routinen =
'alloc_mat_set'
6474 INTEGER :: handle, i_size
6475 LOGICAL :: my_explicitly_no_symmetry
6477 CALL timeset(routinen, handle)
6479 my_explicitly_no_symmetry = .false.
6480 IF (
PRESENT(explicitly_no_symmetry)) my_explicitly_no_symmetry = explicitly_no_symmetry
6483 CALL dbcsr_allocate_matrix_set(mat_set, mat_size)
6484 DO i_size = 1, mat_size
6485 ALLOCATE (mat_set(i_size)%matrix)
6486 IF (my_explicitly_no_symmetry)
THEN
6487 CALL dbcsr_create(matrix=mat_set(i_size)%matrix, template=template, &
6488 matrix_type=dbcsr_type_no_symmetry)
6490 CALL dbcsr_create(matrix=mat_set(i_size)%matrix, template=template)
6492 CALL dbcsr_copy(mat_set(i_size)%matrix, template)
6493 CALL dbcsr_set(mat_set(i_size)%matrix, 0.0_dp)
6496 CALL timestop(handle)
6498 END SUBROUTINE alloc_mat_set
6508 SUBROUTINE alloc_mat_set_2d(mat_set, mat_size_1, mat_size_2, template, explicitly_no_symmetry)
6509 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_set
6510 INTEGER,
INTENT(IN) :: mat_size_1, mat_size_2
6511 TYPE(dbcsr_type),
TARGET :: template
6512 LOGICAL,
OPTIONAL :: explicitly_no_symmetry
6514 CHARACTER(LEN=*),
PARAMETER :: routinen =
'alloc_mat_set_2d'
6516 INTEGER :: handle, i_size, j_size
6517 LOGICAL :: my_explicitly_no_symmetry
6519 CALL timeset(routinen, handle)
6521 my_explicitly_no_symmetry = .false.
6522 IF (
PRESENT(explicitly_no_symmetry)) my_explicitly_no_symmetry = explicitly_no_symmetry
6525 CALL dbcsr_allocate_matrix_set(mat_set, mat_size_1, mat_size_2)
6526 DO i_size = 1, mat_size_1
6527 DO j_size = 1, mat_size_2
6528 ALLOCATE (mat_set(i_size, j_size)%matrix)
6529 IF (my_explicitly_no_symmetry)
THEN
6530 CALL dbcsr_create(matrix=mat_set(i_size, j_size)%matrix, template=template, &
6531 matrix_type=dbcsr_type_no_symmetry)
6533 CALL dbcsr_create(matrix=mat_set(i_size, j_size)%matrix, template=template)
6535 CALL dbcsr_copy(mat_set(i_size, j_size)%matrix, template)
6536 CALL dbcsr_set(mat_set(i_size, j_size)%matrix, 0.0_dp)
6540 CALL timestop(handle)
6542 END SUBROUTINE alloc_mat_set_2d
6557 SUBROUTINE contract_to_self_energy(t_3c_O_all, t_greens_fct, t_3c_O_W, &
6558 mat_self_energy_ao_ao, bounds_ao_ao_j, bounds_RI_i, &
6559 unit_nr, eps_filter, do_occ, do_virt)
6561 TYPE(dbt_type) :: t_3c_o_all, t_greens_fct, t_3c_o_w
6562 TYPE(dbcsr_type),
TARGET :: mat_self_energy_ao_ao
6563 INTEGER,
DIMENSION(2, 2) :: bounds_ao_ao_j
6564 INTEGER,
DIMENSION(2, 1) :: bounds_ri_i
6566 REAL(kind=dp) :: eps_filter
6567 LOGICAL :: do_occ, do_virt
6569 CHARACTER(LEN=*),
PARAMETER :: routinen =
'contract_to_self_energy'
6572 INTEGER,
DIMENSION(2, 1) :: bounds_ao_j
6573 INTEGER,
DIMENSION(2, 2) :: bounds_ao_all_ri_i, bounds_ri_i_ao_j
6574 REAL(kind=dp) :: sign_self_energy
6575 TYPE(dbt_type) :: t_3c_o_g, t_3c_o_g_tmp, t_self_energy, &
6578 CALL timeset(routinen, handle)
6580 cpassert(do_occ .EQV. (.NOT. do_virt))
6582 CALL dbt_create(t_3c_o_all, t_3c_o_g, name=
"M occ (RI AO | AO)")
6583 CALL dbt_create(t_3c_o_all, t_3c_o_g_tmp, name=
"M occ (RI AO | AO)")
6584 CALL dbt_create(t_greens_fct, t_self_energy, name=
"(AO|AO)")
6585 CALL dbt_create(mat_self_energy_ao_ao, t_self_energy_tmp)
6587 bounds_ao_j(:, 1) = bounds_ao_ao_j(:, 1)
6588 bounds_ao_all_ri_i(:, 1) = bounds_ri_i(:, 1)
6589 bounds_ao_all_ri_i(:, 2) = bounds_ao_ao_j(:, 2)
6591 CALL dbt_contract(1.0_dp, t_greens_fct, t_3c_o_all, 0.0_dp, &
6593 contract_1=[2], notcontract_1=[1], &
6594 contract_2=[3], notcontract_2=[1, 2], &
6595 map_1=[3], map_2=[1, 2], &
6596 bounds_2=bounds_ao_j, &
6597 bounds_3=bounds_ao_all_ri_i, &
6598 filter_eps=eps_filter, &
6601 CALL dbt_copy(t_3c_o_g_tmp, t_3c_o_g, order=[1, 3, 2], move_data=.true.)
6603 IF (do_occ) sign_self_energy = -1.0_dp
6604 IF (do_virt) sign_self_energy = 1.0_dp
6606 bounds_ri_i_ao_j(:, 1) = bounds_ri_i(:, 1)
6607 bounds_ri_i_ao_j(:, 2) = bounds_ao_ao_j(:, 1)
6609 CALL dbt_contract(sign_self_energy, t_3c_o_w, t_3c_o_g, 0.0_dp, &
6611 contract_1=[1, 2], notcontract_1=[3], &
6612 contract_2=[1, 2], notcontract_2=[3], &
6613 map_1=[1], map_2=[2], &
6614 bounds_1=bounds_ri_i_ao_j, &
6615 filter_eps=eps_filter, &
6618 CALL dbt_copy(t_self_energy, t_self_energy_tmp)
6619 CALL dbt_clear(t_self_energy)
6621 CALL dbt_copy_tensor_to_matrix(t_self_energy_tmp, mat_self_energy_ao_ao, summation=.true.)
6623 CALL dbt_destroy(t_3c_o_g)
6624 CALL dbt_destroy(t_3c_o_g_tmp)
6625 CALL dbt_destroy(t_self_energy)
6626 CALL dbt_destroy(t_self_energy_tmp)
6628 CALL timestop(handle)
6630 END SUBROUTINE contract_to_self_energy
6645 SUBROUTINE contract_cubic_gw(t_3c_overl_int_gw_AO, t_3c_overl_int_gw_RI, &
6646 t_AO, t_RI, prefac, &
6647 mo_bounds, unit_nr, &
6648 t_3c_ctr_RI, t_3c_ctr_AO, calculate_ctr_RI)
6649 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_overl_int_gw_ao, &
6650 t_3c_overl_int_gw_ri, t_ao, t_ri
6651 REAL(dp),
DIMENSION(2),
INTENT(IN) :: prefac
6652 INTEGER,
DIMENSION(2),
INTENT(IN) :: mo_bounds
6653 INTEGER,
INTENT(IN) :: unit_nr
6654 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_ctr_ri, t_3c_ctr_ao
6655 LOGICAL,
INTENT(IN) :: calculate_ctr_ri
6657 CHARACTER(LEN=*),
PARAMETER :: routinen =
'contract_cubic_gw'
6660 INTEGER,
DIMENSION(2, 2) :: ctr_bounds_mo
6661 INTEGER,
DIMENSION(3) :: bounds_3c
6663 CALL timeset(routinen, handle)
6665 IF (calculate_ctr_ri)
THEN
6666 CALL dbt_get_info(t_3c_overl_int_gw_ri, nfull_total=bounds_3c)
6667 ctr_bounds_mo(:, 1) = [1, bounds_3c(2)]
6668 ctr_bounds_mo(:, 2) = mo_bounds
6670 CALL dbt_contract(prefac(1), t_ri, t_3c_overl_int_gw_ri, 0.0_dp, &
6672 contract_1=[2], notcontract_1=[1], &
6673 contract_2=[1], notcontract_2=[2, 3], &
6674 map_1=[1], map_2=[2, 3], &
6675 bounds_3=ctr_bounds_mo, &
6680 CALL dbt_get_info(t_3c_overl_int_gw_ao, nfull_total=bounds_3c)
6681 ctr_bounds_mo(:, 1) = [1, bounds_3c(2)]
6682 ctr_bounds_mo(:, 2) = mo_bounds
6684 CALL dbt_contract(prefac(2), t_ao, t_3c_overl_int_gw_ao, 0.0_dp, &
6686 contract_1=[2], notcontract_1=[1], &
6687 contract_2=[1], notcontract_2=[2, 3], &
6688 map_1=[1], map_2=[2, 3], &
6689 bounds_3=ctr_bounds_mo, &
6692 CALL timestop(handle)
6694 END SUBROUTINE contract_cubic_gw
6705 SUBROUTINE trace_sigma_gw(t3c_1, t3c_2, vec_sigma, mo_offset, mo_bounds, para_env)
6706 TYPE(dbt_type),
INTENT(INOUT) :: t3c_1, t3c_2
6707 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: vec_sigma
6708 INTEGER,
INTENT(IN) :: mo_offset
6709 INTEGER,
DIMENSION(2),
INTENT(IN) :: mo_bounds
6710 TYPE(mp_para_env_type),
INTENT(IN) :: para_env
6712 CHARACTER(LEN=*),
PARAMETER :: routinen =
'trace_sigma_gw'
6714 INTEGER :: handle, n, n_end, n_end_block, n_start, &
6716 INTEGER,
DIMENSION(1) :: trace_shape
6717 INTEGER,
DIMENSION(2) :: mo_bounds_off
6718 INTEGER,
DIMENSION(3) :: boff, bsize, ind
6720 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: block_1, block_2
6722 DIMENSION(mo_bounds(2)-mo_bounds(1)+1) :: vec_sigma_prv
6723 TYPE(dbt_iterator_type) :: iter
6724 TYPE(dbt_type) :: t3c_1_redist
6726 CALL timeset(routinen, handle)
6728 CALL dbt_create(t3c_2, t3c_1_redist)
6729 CALL dbt_copy(t3c_1, t3c_1_redist, order=[2, 1, 3], move_data=.true.)
6731 vec_sigma_prv = 0.0_dp
6737 CALL dbt_iterator_start(iter, t3c_1_redist)
6738 DO WHILE (dbt_iterator_blocks_left(iter))
6739 CALL dbt_iterator_next_block(iter, ind, blk_size=bsize, blk_offset=boff)
6740 CALL dbt_get_block(t3c_1_redist, ind, block_1, found)
6742 CALL dbt_get_block(t3c_2, ind, block_2, found)
6743 IF (.NOT. found) cycle
6745 IF (boff(3) < mo_bounds(1))
THEN
6746 n_start_block = mo_bounds(1) - boff(3) + 1
6750 n_start = boff(3) - mo_bounds(1) + 1
6753 IF (boff(3) + bsize(3) - 1 > mo_bounds(2))
THEN
6754 n_end_block = mo_bounds(2) - boff(3) + 1
6755 n_end = mo_bounds(2) - mo_bounds(1) + 1
6757 n_end_block = bsize(3)
6758 n_end = boff(3) + bsize(3) - mo_bounds(1)
6761 trace_shape(1) =
SIZE(block_1, 1)*
SIZE(block_1, 2)
6762 vec_sigma_prv(n_start:n_end) = &
6763 vec_sigma_prv(n_start:n_end) + &
6764 [(dot_product(reshape(block_1(:, :, n), trace_shape), &
6765 reshape(block_2(:, :, n), trace_shape)), &
6766 n=n_start_block, n_end_block)]
6767 DEALLOCATE (block_1, block_2)
6769 CALL dbt_iterator_stop(iter)
6772 CALL dbt_destroy(t3c_1_redist)
6774 CALL para_env%sum(vec_sigma_prv)
6776 mo_bounds_off = mo_bounds - mo_offset + 1
6777 vec_sigma(mo_bounds_off(1):mo_bounds_off(2)) = &
6778 vec_sigma(mo_bounds_off(1):mo_bounds_off(2)) + vec_sigma_prv
6780 CALL timestop(handle)
6781 END SUBROUTINE trace_sigma_gw
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
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, nmo_retained)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
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)
...
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_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
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, nrow, ncol, set_zero)
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, zeff, stride, max_file_size_mb, 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, overlap_rs)
generate real space density matrices in DBCSR format
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
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 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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
subroutine, public kpoint_create(kpoint)
Create 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, kp_shift, gamma_centered)
diagonalize KS matrices at a set of kpoints
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, 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, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public qs_env_release(qs_env)
releases the given qs_env (see doc/ReferenceCounting.html)
Initialize a qs_env for kpoint calculations starting from a gamma point qs_env.
subroutine, public create_kp_from_gamma(qs_env, qs_env_kp, with_xc_terms)
...
Some utility functions for the calculation of integrals.
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition and initialisation of the mo data type.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>.
subroutine, public build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
...
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Generate the atomic neighbor lists.
subroutine, public setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, mic, symmetric, molecular, operator_type)
Build a neighborlist.
Calculation of overlap matrix, its derivatives and forces.
subroutine, public build_overlap_matrix_simple(ks_env, matrix_s, basis_set_list_a, basis_set_list_b, sab_nl, lcart)
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 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, cfm_mo_coeff, 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 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_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 create_propagator_matrix_set(propagator, ntime, matrix_template, index_to_cell)
Creates the sector-, time-, and cell-resolved propagator matrix set.
integer, parameter, public propagator_sector_virtual
subroutine, public compute_periodic_dm(propagator, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, sector)
...
integer, parameter, public propagator_sector_occupied
subroutine, public compute_gamma_propagator(propagator, jquad, cfm_mo_coeff, homo, eigenval, nmo, eps_filter, e_fermi, tau, para_env)
Builds the Gamma-point occupied and virtual propagators.
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 ...