From 024676661e34077435b59312aca1eae59daaf8ed Mon Sep 17 00:00:00 2001 From: Jan Wilhelm Date: Mon, 6 Jun 2022 11:03:09 +0200 Subject: [PATCH] 4-center Hartree-Fock and ADMM for exchange self-energy in GW+kpoints --- src/input_cp2k_mp2.F | 22 +- src/mp2.F | 3 +- src/mp2_integrals.F | 2 +- src/mp2_ri_2c.F | 18 +- src/mp2_setup.F | 5 + src/mp2_types.F | 9 +- src/qs_band_structure.F | 64 ++-- src/rpa_gw.F | 80 ++--- src/rpa_gw_sigma.F | 93 +++-- src/rpa_kpoints.F | 9 +- src/rpa_main.F | 322 ++++++++++++++---- src/xas_tdp_atom.F | 2 +- tests/QS/regtest-gw-cubic/TEST_FILES | 3 - .../G0W0_kpoints_from_Gamma.inp | 0 ...ints_from_Gamma_RI_regularization_1E-3.inp | 0 .../G0W0_kpoints_in_self_energy.inp | 0 ...elf_energy_Sigmax_from_four_center_HFX.inp | 105 ++++++ ...nergy_Sigmax_from_four_center_HFX_ADMM.inp | 130 +++++++ .../G0W0_kpoints_in_self_energy_at_HSE06.inp} | 59 +++- tests/QS/regtest-gw-kpoints/TEST_FILES | 7 + tests/QS/regtest-rpa-cubic-scaling/TEST_FILES | 5 +- tests/TEST_DIRS | 1 + 22 files changed, 752 insertions(+), 187 deletions(-) rename tests/QS/{regtest-gw-cubic => regtest-gw-kpoints}/G0W0_kpoints_from_Gamma.inp (100%) rename tests/QS/{regtest-gw-cubic => regtest-gw-kpoints}/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp (100%) rename tests/QS/{regtest-gw-cubic => regtest-gw-kpoints}/G0W0_kpoints_in_self_energy.inp (100%) create mode 100644 tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX.inp create mode 100644 tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX_ADMM.inp rename tests/QS/{regtest-rpa-cubic-scaling/RPA_kpoints_from_Gamma_H2O.inp => regtest-gw-kpoints/G0W0_kpoints_in_self_energy_at_HSE06.inp} (56%) create mode 100644 tests/QS/regtest-gw-kpoints/TEST_FILES diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F index e777b03aa9..d6a7eb7990 100644 --- a/src/input_cp2k_mp2.F +++ b/src/input_cp2k_mp2.F @@ -1224,7 +1224,17 @@ CONTAINS "means smaller expansion coefficients that leads to a more stable calculation at the price "// & "of a slightly worse RI approximation. In case the parameter 0.0 is chosen, ordinary RI is used.", & usage="REGULARIZATION_RI 1.0E-4", & - default_r_val=1.0E-2_dp) + default_r_val=0.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create( & + keyword, __LOCATION__, & + name="EPS_EIGVAL_S", & + description="Parameter to reduce the expansion coefficients in RI for periodic GW. Removes all "// & + "eigenvectors and eigenvalues of S_PQ(k) that are smaller than EPS_EIGVAL_S. ", & + usage="EPS_EIGVAL_S 1.0E-3", & + default_r_val=1.0E-5_dp) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -1251,6 +1261,16 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create( & + keyword, __LOCATION__, & + name="REL_CUTOFF_TRUNC_COLOUMB_RI_X", & + description="Only active in case TRUNC_COLOUMB_RI_X = True. Normally, relative cutoff = 0.5 is "// & + "good choice; still needs to be evaluated for RI schemes. ", & + usage="REL_CUTOFF_TRUNC_COLOUMB_RI_X 0.3", & + default_r_val=0.5_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create( & keyword, __LOCATION__, & name="KEEP_QUADRATURE", & diff --git a/src/mp2.F b/src/mp2.F index ac837f7904..334fecb155 100644 --- a/src/mp2.F +++ b/src/mp2.F @@ -697,7 +697,8 @@ CONTAINS DEALLOCATE (mos_mp2) ! if necessary reallocate hfx buffer - IF (free_hfx_buffer .AND. (.NOT. calc_forces)) THEN + IF (free_hfx_buffer .AND. (.NOT. calc_forces) .AND. & + (mp2_env%ri_g0w0%do_ri_Sigma_x .OR. .NOT. mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma)) THEN CALL timeset(routineN//"_alloc_hfx", handle2) DO irep = 1, n_rep_hf DO i_thread = 0, n_threads - 1 diff --git a/src/mp2_integrals.F b/src/mp2_integrals.F index 78c64186ac..bfa7c017b0 100644 --- a/src/mp2_integrals.F +++ b/src/mp2_integrals.F @@ -313,7 +313,7 @@ CONTAINS CPABORT("SVD not implemented for forces.!") END IF - do_kpoints_from_Gamma = SUM(qs_env%mp2_env%ri_rpa_im_time%kp_grid) > 0 + do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN CALL get_qs_env(qs_env=qs_env, & kpoints=kpoints) diff --git a/src/mp2_ri_2c.F b/src/mp2_ri_2c.F index e3617edf54..2d94ec9a09 100644 --- a/src/mp2_ri_2c.F +++ b/src/mp2_ri_2c.F @@ -260,7 +260,8 @@ CONTAINS do_kpoints, kpoints, put_mat_KS_env=.FALSE.) CALL inversion_of_S_and_mult_with_chol_dec_of_V(fm_matrix_L_RI_metric, fm_matrix_L_kpoints, & - fm_matrix_Sinv_Vtrunc_Sinv, dimen_RI, kpoints) + fm_matrix_Sinv_Vtrunc_Sinv, dimen_RI, kpoints, & + qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S) ELSE IF (calc_forces .AND. (.NOT. do_im_time)) THEN @@ -1404,15 +1405,18 @@ CONTAINS !> \param fm_matrix_Sinv_Vtrunc_Sinv ... !> \param dimen_RI ... !> \param kpoints ... +!> \param eps_eigval_S ... ! ************************************************************************************************** SUBROUTINE inversion_of_S_and_mult_with_chol_dec_of_V(fm_matrix_L_RI_metric, fm_matrix_L_kpoints, & - fm_matrix_Sinv_Vtrunc_Sinv, dimen_RI, kpoints) + fm_matrix_Sinv_Vtrunc_Sinv, dimen_RI, kpoints, & + eps_eigval_S) TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_matrix_L_RI_metric, & fm_matrix_L_kpoints, & fm_matrix_Sinv_Vtrunc_Sinv INTEGER, INTENT(IN) :: dimen_RI TYPE(kpoint_type), POINTER :: kpoints + REAL(KIND=dp), INTENT(IN) :: eps_eigval_S CHARACTER(LEN=*), PARAMETER :: routineN = 'inversion_of_S_and_mult_with_chol_dec_of_V' COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), & @@ -1449,7 +1453,7 @@ CONTAINS CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_Vtrunc_tmp, cone, fm_matrix_Sinv_Vtrunc_Sinv(ikp, 1)%matrix) CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_Vtrunc_tmp, ione, fm_matrix_Sinv_Vtrunc_Sinv(ikp, 2)%matrix) - CALL cp_cfm_robust_cholesky(cfm_matrix_S_tmp, cfm_matrix_S_inv_tmp, threshold=1.0E-12_dp, exponent=-0.5_dp) + CALL cp_cfm_robust_cholesky(cfm_matrix_S_tmp, cfm_matrix_S_inv_tmp, threshold=eps_eigval_S, exponent=-0.5_dp) ! get S^(-1) = U^(-H)U^(-1) CALL cp_cfm_gemm("C", "N", dimen_RI, dimen_RI, dimen_RI, cone, cfm_matrix_S_tmp, cfm_matrix_S_tmp, & @@ -1501,7 +1505,8 @@ CONTAINS INTEGER :: handle INTEGER, DIMENSION(3) :: periodic - REAL(KIND=dp) :: shortest_dist_cell_planes + REAL(KIND=dp) :: rel_cutoff_trunc_coulomb_ri_x, & + shortest_dist_cell_planes REAL(KIND=dp), DIMENSION(2) :: abs_cutoffs_chi_W TYPE(cell_type), POINTER :: cell @@ -1534,9 +1539,12 @@ CONTAINS qs_env%mp2_env%ri_rpa_im_time%abs_cutoffs_chi_W(1:2) = abs_cutoffs_chi_W(1:2) + rel_cutoff_trunc_coulomb_ri_x = qs_env%mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x + IF (PRESENT(trunc_coulomb)) THEN trunc_coulomb%potential_type = do_potential_truncated - trunc_coulomb%cutoff_radius = (abs_cutoffs_chi_W(1) + abs_cutoffs_chi_W(2))*0.5_dp + trunc_coulomb%cutoff_radius = (abs_cutoffs_chi_W(1) + abs_cutoffs_chi_W(2))* & + rel_cutoff_trunc_coulomb_ri_x trunc_coulomb%filename = "t_c_g.dat" ! dummy trunc_coulomb%omega = 0.0_dp diff --git a/src/mp2_setup.F b/src/mp2_setup.F index 7153cfa51d..b07cfcb0e3 100644 --- a/src/mp2_setup.F +++ b/src/mp2_setup.F @@ -238,6 +238,7 @@ CONTAINS l_val=mp2_env%ri_rpa_im_time%do_im_time_kpoints) CALL section_vals_val_get(low_scaling_section, "KPOINTS", & i_vals=mp2_env%ri_rpa_im_time%kp_grid) + mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma = SUM(mp2_env%ri_rpa_im_time%kp_grid) > 0 CALL section_vals_val_get(low_scaling_section, "KPOINT_WEIGHTS_W", & i_val=mp2_env%ri_rpa_im_time%kpoint_weights_W_method) CALL section_vals_val_get(low_scaling_section, "EXPONENT_TAILORED_WEIGHTS", & @@ -246,10 +247,14 @@ CONTAINS r_vals=mp2_env%ri_rpa_im_time%rel_cutoffs_chi_W) CALL section_vals_val_get(low_scaling_section, "REGULARIZATION_RI", & r_val=mp2_env%ri_rpa_im_time%regularization_RI) + CALL section_vals_val_get(low_scaling_section, "EPS_EIGVAL_S", & + r_val=mp2_env%ri_rpa_im_time%eps_eigval_S) CALL section_vals_val_get(low_scaling_section, "MAKE_CHI_POS_DEFINITE", & l_val=mp2_env%ri_rpa_im_time%make_chi_pos_definite) CALL section_vals_val_get(low_scaling_section, "TRUNC_COLOUMB_RI_X", & l_val=mp2_env%ri_rpa_im_time%trunc_coulomb_ri_x) + CALL section_vals_val_get(low_scaling_section, "REL_CUTOFF_TRUNC_COLOUMB_RI_X", & + r_val=mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x) CALL section_vals_val_get(low_scaling_section, "K_MESH_G_FACTOR", & i_val=mp2_env%ri_rpa_im_time%k_mesh_g_factor) CALL section_vals_val_get(low_scaling_section, "GREENS_FUNCTION", & diff --git a/src/mp2_types.F b/src/mp2_types.F index 18255b1546..9e0718aedd 100644 --- a/src/mp2_types.F +++ b/src/mp2_types.F @@ -153,9 +153,11 @@ MODULE mp2_types TYPE ri_rpa_im_time_type INTEGER :: cut_memory - LOGICAL :: memory_info, make_chi_pos_definite, trunc_coulomb_ri_x, keep_quad + LOGICAL :: memory_info, make_chi_pos_definite, trunc_coulomb_ri_x, keep_quad, & + do_kpoints_from_Gamma REAL(KIND=dp) :: eps_filter, & - eps_filter_factor, eps_compress, exp_tailored_weights, regularization_RI + eps_filter_factor, eps_compress, exp_tailored_weights, regularization_RI, & + eps_eigval_S, rel_cutoff_trunc_coulomb_ri_x REAL(KIND=dp), DIMENSION(:), POINTER :: rel_cutoffs_chi_W, tau_tj, tau_wj, tj, wj REAL(KIND=dp), DIMENSION(:, :), POINTER :: weights_cos_tf_t_to_w, weights_cos_tf_w_to_t REAL(KIND=dp), DIMENSION(2) :: abs_cutoffs_chi_W @@ -211,7 +213,8 @@ MODULE mp2_types INTEGER :: n_kp_in_kp_line, n_special_kp, nkp_self_energy, & nkp_self_energy_special_kp, nkp_self_energy_monkh_pack REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp_special_kp - TYPE(dbcsr_p_type), DIMENSION(:), ALLOCATABLE :: mat_exchange_for_kp_from_gamma + TYPE(dbcsr_p_type), DIMENSION(:), ALLOCATABLE :: & + matrix_sigma_x_minus_vxc, matrix_ks END TYPE TYPE ri_basis_opt diff --git a/src/qs_band_structure.F b/src/qs_band_structure.F index cba0b8e407..2c0fef8e4f 100644 --- a/src/qs_band_structure.F +++ b/src/qs_band_structure.F @@ -63,7 +63,7 @@ MODULE qs_band_structure CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_band_structure' - PUBLIC :: calculate_band_structure, calculate_kp_orbitals + PUBLIC :: calculate_band_structure, calculate_kp_orbitals, calculate_kpoints_for_bs ! ************************************************************************************************** @@ -310,7 +310,6 @@ CONTAINS OPTIONAL :: kpgeneral INTEGER, INTENT(IN), OPTIONAL :: group_size_ext - INTEGER :: i, ix, iy, iz, npoints TYPE(cp_blacs_env_type), POINTER :: blacs_env TYPE(cp_para_env_type), POINTER :: para_env TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s @@ -320,6 +319,46 @@ CONTAINS TYPE(qs_scf_env_type), POINTER :: scf_env TYPE(scf_control_type), POINTER :: scf_control + CALL calculate_kpoints_for_bs(kpoint, scheme, group_size_ext, mp_grid, kpgeneral) + + CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env) + kpoint%para_env => para_env + CALL cp_para_env_retain(para_env) + kpoint%blacs_env_all => blacs_env + CALL cp_blacs_env_retain(blacs_env) + CALL kpoint_env_initialize(kpoint) + + CALL kpoint_initialize_mos(kpoint, qs_env%mos, nadd) + CALL kpoint_initialize_mo_set(kpoint) + + CALL get_qs_env(qs_env, sab_kp=sab_nl, dft_control=dft_control) + CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control) + + CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, & + scf_env=scf_env, scf_control=scf_control) + CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE.) + + END SUBROUTINE calculate_kp_orbitals + +! ************************************************************************************************** +!> \brief ... +!> \param kpoint ... +!> \param scheme ... +!> \param group_size_ext ... +!> \param mp_grid ... +!> \param kpgeneral ... +! ************************************************************************************************** + SUBROUTINE calculate_kpoints_for_bs(kpoint, scheme, group_size_ext, mp_grid, kpgeneral) + + TYPE(kpoint_type), POINTER :: kpoint + CHARACTER(LEN=*), INTENT(IN) :: scheme + INTEGER, INTENT(IN), OPTIONAL :: group_size_ext + INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: mp_grid + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & + OPTIONAL :: kpgeneral + + INTEGER :: i, ix, iy, iz, npoints + CPASSERT(.NOT. ASSOCIATED(kpoint)) CALL kpoint_create(kpoint) @@ -389,25 +428,6 @@ CONTAINS CPABORT("Unknown kpoint scheme requested") END SELECT - CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env) - kpoint%para_env => para_env - CALL cp_para_env_retain(para_env) - kpoint%blacs_env_all => blacs_env - CALL cp_blacs_env_retain(blacs_env) - CALL kpoint_env_initialize(kpoint) - - CALL kpoint_initialize_mos(kpoint, qs_env%mos, nadd) - CALL kpoint_initialize_mo_set(kpoint) - - CALL get_qs_env(qs_env, sab_kp=sab_nl, dft_control=dft_control) - CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control) - - CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, & - scf_env=scf_env, scf_control=scf_control) - CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE.) - - END SUBROUTINE calculate_kp_orbitals - -! ************************************************************************************************** + END SUBROUTINE calculate_kpoints_for_bs END MODULE qs_band_structure diff --git a/src/rpa_gw.F b/src/rpa_gw.F index 23ee3355ba..bfb8f67e02 100644 --- a/src/rpa_gw.F +++ b/src/rpa_gw.F @@ -15,11 +15,15 @@ MODULE rpa_gw gto_basis_set_type USE cell_types, ONLY: cell_type,& get_cell - USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm,& - cp_cfm_scale_and_add - USE cp_cfm_types, ONLY: & - cp_cfm_create, cp_cfm_get_info, cp_cfm_p_type, cp_cfm_release, cp_cfm_set_all, & - cp_cfm_to_cfm, cp_cfm_to_fm, cp_cfm_type, cp_fm_to_cfm + USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm + USE cp_cfm_types, ONLY: cp_cfm_create,& + cp_cfm_get_info,& + cp_cfm_p_type,& + cp_cfm_release,& + cp_cfm_set_all,& + cp_cfm_to_fm,& + cp_cfm_type,& + cp_fm_to_cfm USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& copy_fm_to_dbcsr,& dbcsr_allocate_matrix_set,& @@ -1296,7 +1300,7 @@ CONTAINS qs_env, para_env, & mp2_env, num_fit_points, mo_coeff, & do_ri_Sigma_x, vec_Sigma_x_gw(:, :, 1), unit_nr, 1, & - starts_array_mc, ends_array_mc) + starts_array_mc, ends_array_mc, eps_filter) END IF @@ -4672,6 +4676,7 @@ CONTAINS !> \param ispin ... !> \param starts_array_mc ... !> \param ends_array_mc ... +!> \param eps_filter ... ! ************************************************************************************************** SUBROUTINE compute_self_energy_cubic_gw_kpoints(num_integ_points, tau_tj, tj, & matrix_s, Eigenval, e_fermi, fm_mat_W, & @@ -4684,7 +4689,7 @@ CONTAINS qs_env, para_env, & mp2_env, num_fit_points, fm_mo_coeff, & do_ri_Sigma_x, vec_Sigma_x_gw, unit_nr, ispin, & - starts_array_mc, ends_array_mc) + starts_array_mc, ends_array_mc, eps_filter) INTEGER, INTENT(IN) :: num_integ_points REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & @@ -4718,6 +4723,7 @@ CONTAINS REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vec_Sigma_x_gw INTEGER, INTENT(IN) :: unit_nr, ispin INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc + REAL(KIND=dp), INTENT(IN) :: eps_filter CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_self_energy_cubic_gw_kpoints' @@ -4921,6 +4927,7 @@ CONTAINS map_1=[1], map_2=[2, 3], & bounds_2=bounds_RI_i, & bounds_3=bounds_ao_ao_j, & + filter_eps=eps_filter, & unit_nr=unit_nr_prv) CALL dbt_copy(t_3c_M_W_tmp, t_3c_O_W, order=[1, 2, 3], move_data=.TRUE.) @@ -4930,12 +4937,12 @@ CONTAINS CALL contract_to_self_energy(t_3c_O_all, t_greens_fct_occ, t_3c_O_W, & mat_self_energy_ao_ao_neg_tau, & bounds_ao_ao_j, bounds_RI_i, unit_nr_prv, & - do_occ=.TRUE., do_virt=.FALSE.) + eps_filter, do_occ=.TRUE., do_virt=.FALSE.) CALL contract_to_self_energy(t_3c_O_all, t_greens_fct_virt, t_3c_O_W, & mat_self_energy_ao_ao_pos_tau, & bounds_ao_ao_j, bounds_RI_i, unit_nr_prv, & - do_occ=.FALSE., do_virt=.TRUE.) + eps_filter, do_occ=.FALSE., do_virt=.TRUE.) END DO ! j_mem @@ -4978,23 +4985,12 @@ CONTAINS IF (count_ev_sc_GW == 1 .AND. count_sc_GW0 == 1) THEN - IF (.NOT. do_ri_Sigma_x) THEN - - CALL mp_sync(para_env%group) - - CALL trafo_to_mo_and_kpoints(qs_env, qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(1)%matrix, & - vec_Sigma_x_gw(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt, :), & - homo, gw_corr_lev_occ, gw_corr_lev_virt) - - CALL dbcsr_release(qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(1)%matrix) - DEALLOCATE (qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(1)%matrix) - DEALLOCATE (qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma) - END IF - CALL compute_minus_vxc_kpoints(qs_env) - mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, :) = mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, :) + & - vec_Sigma_x_gw(:, :) + IF (do_ri_Sigma_x) THEN + mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, :) = mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, ispin, :) + & + vec_Sigma_x_gw(:, :) + END IF END IF @@ -5066,12 +5062,13 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_minus_vxc_kpoints' INTEGER :: handle, ikp, nkp_Sigma, nmo - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag_vxc + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag_Sigma_x_minus_vxc_mo_mo TYPE(cp_cfm_type), POINTER :: cfm_mo_coeff, ks_mat_ao_ao, & ks_mat_no_xc_ao_ao, vxc_ao_ao, & vxc_ao_mo, vxc_mo_mo TYPE(cp_fm_struct_type), POINTER :: matrix_struct - TYPE(cp_fm_type), POINTER :: fm_tmp_im, fm_tmp_re, fm_vxc_mo_mo + TYPE(cp_fm_type), POINTER :: fm_Sigma_x_minus_vxc_mo_mo, fm_tmp_im, & + fm_tmp_re TYPE(cp_para_env_type), POINTER :: para_env TYPE(kpoint_type), POINTER :: kpoints_Sigma, kpoints_Sigma_no_xc @@ -5087,7 +5084,7 @@ CONTAINS matrix_struct => kpoints_Sigma%kp_env(1)%kpoint_env%wmat(1, 1)%matrix%matrix_struct - NULLIFY (ks_mat_ao_ao, ks_mat_no_xc_ao_ao, vxc_ao_ao, vxc_ao_mo, vxc_mo_mo, cfm_mo_coeff, fm_vxc_mo_mo, & + NULLIFY (ks_mat_ao_ao, ks_mat_no_xc_ao_ao, vxc_ao_ao, vxc_ao_mo, vxc_mo_mo, cfm_mo_coeff, fm_Sigma_x_minus_vxc_mo_mo, & fm_tmp_re, fm_tmp_im) CALL cp_cfm_create(ks_mat_ao_ao, matrix_struct) CALL cp_cfm_create(ks_mat_no_xc_ao_ao, matrix_struct) @@ -5095,12 +5092,12 @@ CONTAINS CALL cp_cfm_create(vxc_ao_mo, matrix_struct) CALL cp_cfm_create(vxc_mo_mo, matrix_struct) CALL cp_cfm_create(cfm_mo_coeff, matrix_struct) - CALL cp_fm_create(fm_vxc_mo_mo, matrix_struct) + CALL cp_fm_create(fm_Sigma_x_minus_vxc_mo_mo, matrix_struct) CALL cp_fm_create(fm_tmp_re, matrix_struct) CALL cp_fm_create(fm_tmp_im, matrix_struct) CALL cp_cfm_get_info(cfm_mo_coeff, nrow_global=nmo) - ALLOCATE (diag_vxc(nmo)) + ALLOCATE (diag_Sigma_x_minus_vxc_mo_mo(nmo)) DEALLOCATE (qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw) @@ -5119,19 +5116,22 @@ CONTAINS CALL cp_fm_copy_general(kpoints_Sigma_no_xc%kp_env(ikp)%kpoint_env%wmat(1, 1)%matrix, fm_tmp_re, para_env) CALL cp_fm_copy_general(kpoints_Sigma_no_xc%kp_env(ikp)%kpoint_env%wmat(2, 1)%matrix, fm_tmp_im, para_env) - CALL cp_fm_to_cfm(fm_tmp_re, fm_tmp_im, ks_mat_no_xc_ao_ao) + CALL cp_fm_to_cfm(fm_tmp_re, fm_tmp_im, vxc_ao_ao) - CALL cp_cfm_scale_and_add(z_one, ks_mat_ao_ao, -z_one, ks_mat_no_xc_ao_ao) - CALL cp_cfm_to_cfm(ks_mat_ao_ao, vxc_ao_ao) +! CALL cp_fm_to_cfm(fm_tmp_re, fm_tmp_im, ks_mat_no_xc_ao_ao) + +! CALL cp_cfm_scale_and_add(z_one, ks_mat_ao_ao, -z_one, ks_mat_no_xc_ao_ao) +! CALL cp_cfm_to_cfm(ks_mat_ao_ao, vxc_ao_ao) CALL cp_cfm_gemm('N', 'N', nmo, nmo, nmo, z_one, vxc_ao_ao, cfm_mo_coeff, z_zero, vxc_ao_mo) CALL cp_cfm_gemm('C', 'N', nmo, nmo, nmo, z_one, cfm_mo_coeff, vxc_ao_mo, z_zero, vxc_mo_mo) - CALL cp_cfm_to_fm(vxc_mo_mo, fm_vxc_mo_mo) + CALL cp_cfm_to_fm(vxc_mo_mo, fm_Sigma_x_minus_vxc_mo_mo) - CALL cp_fm_get_diag(fm_vxc_mo_mo, diag_vxc) + CALL cp_fm_get_diag(fm_Sigma_x_minus_vxc_mo_mo, diag_Sigma_x_minus_vxc_mo_mo) - qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp) = -diag_vxc(:) +! qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp) = -diag_Sigma_x_minus_vxc_mo_mo(:) + qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp) = diag_Sigma_x_minus_vxc_mo_mo(:) END DO @@ -5141,11 +5141,11 @@ CONTAINS CALL cp_cfm_release(vxc_ao_mo) CALL cp_cfm_release(vxc_mo_mo) CALL cp_cfm_release(cfm_mo_coeff) - CALL cp_fm_release(fm_vxc_mo_mo) + CALL cp_fm_release(fm_Sigma_x_minus_vxc_mo_mo) CALL cp_fm_release(fm_tmp_re) CALL cp_fm_release(fm_tmp_im) - DEALLOCATE (diag_vxc) + DEALLOCATE (diag_Sigma_x_minus_vxc_mo_mo) CALL timestop(handle) @@ -5335,18 +5335,20 @@ CONTAINS !> \param bounds_ao_ao_j ... !> \param bounds_RI_i ... !> \param unit_nr ... +!> \param eps_filter ... !> \param do_occ ... !> \param do_virt ... ! ************************************************************************************************** SUBROUTINE contract_to_self_energy(t_3c_O_all, t_greens_fct, t_3c_O_W, & mat_self_energy_ao_ao, bounds_ao_ao_j, bounds_RI_i, & - unit_nr, do_occ, do_virt) + unit_nr, eps_filter, do_occ, do_virt) TYPE(dbt_type) :: t_3c_O_all, t_greens_fct, t_3c_O_W TYPE(dbcsr_type), TARGET :: mat_self_energy_ao_ao INTEGER, DIMENSION(2, 2) :: bounds_ao_ao_j INTEGER, DIMENSION(2, 1) :: bounds_RI_i INTEGER :: unit_nr + REAL(KIND=dp) :: eps_filter LOGICAL :: do_occ, do_virt CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_to_self_energy' @@ -5375,6 +5377,7 @@ CONTAINS contract_2=[3], notcontract_2=[1, 2], & map_1=[3], map_2=[1, 2], & bounds_2=bounds_ao_j, & + filter_eps=eps_filter, & unit_nr=unit_nr) CALL dbt_copy(t_3c_O_G_tmp, t_3c_O_G, order=[1, 3, 2], move_data=.TRUE.) @@ -5391,6 +5394,7 @@ CONTAINS contract_2=[1, 2], notcontract_2=[3], & map_1=[1], map_2=[2], & bounds_1=bounds_RI_i_ao_j, & + filter_eps=eps_filter, & unit_nr=unit_nr) CALL dbt_copy(t_self_energy, t_self_energy_tmp) diff --git a/src/rpa_gw_sigma.F b/src/rpa_gw_sigma.F index 7c766fd7dd..74ec7f4bc5 100644 --- a/src/rpa_gw_sigma.F +++ b/src/rpa_gw_sigma.F @@ -121,18 +121,19 @@ CONTAINS CHARACTER(LEN=40) :: line INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, handle, homo, i_img, & ikp, irep, ispin, iunit, myfun, myfun_aux, myfun_prim, n_level_gw, n_level_gw_ref, & - n_rep_hf, nkp, nkp_Sigma, nmo, ns, nspins, print_exx, virtual + n_rep_hf, nkp, nkp_Sigma, nmo, ns, nspins, print_exx LOGICAL :: charge_constrain_tmp, do_admm_rpa, do_hfx, do_kpoints_cubic_RPA, & do_kpoints_from_Gamma, do_ri_Sigma_x, really_read_line REAL(KIND=dp) :: E_GAP_GW, E_HOMO_GW, E_LUMO_GW, eh1, ehfx, eigval_dft, eigval_hf_at_dft, & energy_exc, energy_exc1, energy_exc1_aux_fit, energy_exc_aux_fit, energy_total, & - exx_minus_vxc, hfx_fraction, t1, t2, tmp + exx_minus_vxc, hfx_fraction, min_direct_HF_at_DFT_gap, t1, t2, tmp REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Eigenval_kp_HF_at_DFT, vec_Sigma_x REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_kp, vec_Sigma_x_minus_vxc_gw, & vec_Sigma_x_minus_vxc_gw_im TYPE(admm_type), POINTER :: admm_env TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: mat_exchange_for_kp_from_gamma TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, & matrix_ks_aux_fit_hfx, rho_ao TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_2d, matrix_ks_kp_im, & @@ -159,7 +160,8 @@ CONTAINS do_admm_rpa = mp2_env%ri_rpa%do_admm do_ri_Sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints - do_kpoints_from_Gamma = SUM(qs_env%mp2_env%ri_rpa_im_time%kp_grid) > 0 + do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma + print_exx = mp2_env%ri_g0w0%print_exx IF (do_kpoints_cubic_RPA) THEN CPASSERT(do_ri_Sigma_x) @@ -187,6 +189,8 @@ CONTAINS CPASSERT(admm_env%purification_method == do_admm_purify_none) CPASSERT(dft_control%admm_control%method == do_admm_basis_projection) + nkp = 1 + ELSE IF (do_kpoints_cubic_RPA) THEN @@ -235,6 +239,22 @@ CONTAINS END DO END IF + ! safe ks matrix for later: we will transform matrix_ks + ! to T-cell index and then to k-points for band structure calculation + IF (do_kpoints_from_Gamma) THEN + ! not yet there: open shell + ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(1)) + DO ispin = 1, 1 + NULLIFY (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix) + ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix) + CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix, & + template=matrix_ks(ispin)%matrix) + CALL dbcsr_desymmetrize(matrix_ks(1)%matrix, & + qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix) + + END DO + END IF + IF (do_kpoints_cubic_RPA) THEN CALL allocate_matrix_ks_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints) @@ -462,20 +482,21 @@ CONTAINS END IF END DO - IF (do_kpoints_from_Gamma) THEN - ! not yet there: open shell - ALLOCATE (qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(1)) + IF (do_kpoints_from_Gamma .AND. print_exx == gw_print_exx) THEN + ! JW not yet there: open shell + ALLOCATE (mat_exchange_for_kp_from_gamma(1)) + DO ispin = 1, 1 - NULLIFY (qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(ispin)%matrix) - ALLOCATE (qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(ispin)%matrix) - CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(ispin)%matrix, & - template=matrix_ks(ispin)%matrix) - CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, & - qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(ispin)%matrix) + NULLIFY (mat_exchange_for_kp_from_gamma(ispin)%matrix) + ALLOCATE (mat_exchange_for_kp_from_gamma(ispin)%matrix) + CALL dbcsr_create(mat_exchange_for_kp_from_gamma(ispin)%matrix, template=matrix_ks(ispin)%matrix) + CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, mat_exchange_for_kp_from_gamma(ispin)%matrix) END DO + END IF END IF + energy_ex = ehfx ! transform Fock-Matrix (calculated in integrate_four_center, written in matrix_ks_aux_fit in case @@ -488,6 +509,23 @@ CONTAINS CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, 1.0_dp) END DO + ! safe matrix_sigma_x_minus_vxc for later: for example, we will transform matrix_sigma_x_minus_vxc + ! to T-cell index and then to k-points for band structure calculation + IF (do_kpoints_from_Gamma) THEN + ! not yet there: open shell + ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(1)) + DO ispin = 1, 1 + NULLIFY (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix) + ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix) + CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix, & + template=matrix_ks(ispin)%matrix) + + CALL dbcsr_desymmetrize(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, & + qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix) + + END DO + END IF + CALL dbcsr_desymmetrize(matrix_ks(1)%matrix, mo_coeff_b) CALL dbcsr_set(mo_coeff_b, 0.0_dp) @@ -579,23 +617,22 @@ CONTAINS ALLOCATE (mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp)) - print_exx = mp2_env%ri_g0w0%print_exx - IF (print_exx == gw_print_exx) THEN IF (do_kpoints_from_Gamma) THEN gw_corr_lev_tot = gw_corr_lev_occ + gw_corr_lev_virt - virtual = nmo - homo CALL get_qs_env(qs_env=qs_env, & kpoints=kpoints) + CALL setup_abs_cutoffs_chi_and_trunc_coulomb_potential(qs_env) + CALL compute_kpoints(qs_env, kpoints, unit_nr) ALLOCATE (Eigenval_kp(nmo, 1, nspins)) - CALL get_bandstruc_and_k_dependent_MOs(qs_env, kpoints, virtual, Eigenval_kp) + CALL get_bandstruc_and_k_dependent_MOs(qs_env, kpoints, Eigenval_kp) CALL compute_minus_vxc_kpoints(qs_env) @@ -604,13 +641,15 @@ CONTAINS ALLOCATE (vec_Sigma_x(nmo, nkp_Sigma)) vec_Sigma_x(:, :) = 0.0_dp - CALL setup_abs_cutoffs_chi_and_trunc_coulomb_potential(qs_env) - CALL trafo_to_mo_and_kpoints(qs_env, & - qs_env%mp2_env%ri_g0w0%mat_exchange_for_kp_from_gamma(1)%matrix, & + mat_exchange_for_kp_from_gamma(1)%matrix, & vec_Sigma_x(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt, :), & homo, gw_corr_lev_occ, gw_corr_lev_virt) + CALL dbcsr_release(mat_exchange_for_kp_from_gamma(1)%matrix) + DEALLOCATE (mat_exchange_for_kp_from_gamma(1)%matrix) + DEALLOCATE (mat_exchange_for_kp_from_gamma) + DEALLOCATE (vec_Sigma_x_minus_vxc_gw) ALLOCATE (vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp_Sigma)) @@ -631,6 +670,8 @@ CONTAINS ALLOCATE (Eigenval_kp_HF_at_DFT(nmo, nkp_Sigma)) Eigenval_kp_HF_at_DFT(:, :) = Eigenval_kp(:, :, 1) + vec_Sigma_x_minus_vxc_gw(:, 1, :) + min_direct_HF_at_DFT_gap = 100.0_dp + WRITE (unit_nr, '(T3,A)') '' WRITE (unit_nr, '(T3,A)') 'Exchange energies' WRITE (unit_nr, '(T3,A)') '-----------------' @@ -665,21 +706,29 @@ CONTAINS E_HOMO_GW = MAXVAL(Eigenval_kp_HF_at_DFT(homo - gw_corr_lev_occ + 1:homo, ikp)) E_LUMO_GW = MINVAL(Eigenval_kp_HF_at_DFT(homo + 1:homo + gw_corr_lev_virt, ikp)) E_GAP_GW = E_LUMO_GW - E_HOMO_GW + IF (E_GAP_GW < min_direct_HF_at_DFT_gap) min_direct_HF_at_DFT_gap = E_GAP_GW WRITE (unit_nr, '(T3,A)') '' - WRITE (unit_nr, '(T3,A,F53.2)') 'HF@PBE HOMO-LUMO gap (eV)', E_GAP_GW*evolt + WRITE (unit_nr, '(T3,A,F53.2)') 'HF@DFT HOMO-LUMO gap (eV)', E_GAP_GW*evolt WRITE (unit_nr, '(T3,A)') '' END DO + WRITE (unit_nr, '(T3,A)') '' + WRITE (unit_nr, '(T3,A)') '' + WRITE (unit_nr, '(T3,A,F63.3)') 'HF@DFT direct bandgap (eV)', min_direct_HF_at_DFT_gap*evolt + WRITE (unit_nr, '(T3,A)') '' WRITE (unit_nr, '(T3,A)') 'End of exchange energies' WRITE (unit_nr, '(T3,A)') '------------------------' WRITE (unit_nr, '(T3,A)') '' + + CPABORT('Stop after printing exchange energies.') + + ELSE + CALL mp_sync(para_env%group) END IF END IF - CALL mp_sync(para_env%group) - IF (print_exx == gw_read_exx) THEN CALL open_file(unit_number=iunit, file_name="exx.out") diff --git a/src/rpa_kpoints.F b/src/rpa_kpoints.F index 61edd71ffa..d601121e5e 100644 --- a/src/rpa_kpoints.F +++ b/src/rpa_kpoints.F @@ -65,7 +65,8 @@ MODULE rpa_kpoints CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_kpoints' - PUBLIC :: RPA_postprocessing_kp, cp_cfm_robust_cholesky, real_space_to_kpoint_transform_rpa, get_mat_cell_T_from_mat_gamma + PUBLIC :: RPA_postprocessing_kp, cp_cfm_robust_cholesky, real_space_to_kpoint_transform_rpa, get_mat_cell_T_from_mat_gamma, & + transform_P_from_real_space_to_kpoints CONTAINS @@ -763,13 +764,15 @@ CONTAINS WRITE (unit_nr, '(T3,A,T72,F9.2)') 'GW_INFO| Relative cutoff 2: ', & qs_env%mp2_env%ri_rpa_im_time%rel_cutoffs_chi_W(2) WRITE (unit_nr, '(T3,A,T72,F9.2)') 'GW_INFO| Cutoff radius 1 (Angstroem):', & - cp_unit_from_cp2k(first_cutoff, "angstrom") + cp_unit_from_cp2k(first_cutoff, "Angstrom") WRITE (unit_nr, '(T3,A,T72,F9.2)') 'GW_INFO| Cutoff radius 2 (Angstroem):', & - cp_unit_from_cp2k(second_cutoff, "angstrom") + cp_unit_from_cp2k(second_cutoff, "Angstrom") WRITE (unit_nr, '(T3,A,T72,F9.2)') & 'GW_INFO| If Cutoff radius 1 = Cutoff radius 2 = 0.5: minimum image convention.' WRITE (unit_nr, '(T3,A,T66,ES15.2)') 'GW_INFO| RI regularization parameter: ', & qs_env%mp2_env%ri_rpa_im_time%regularization_RI + WRITE (unit_nr, '(T3,A,T66,ES15.2)') 'GW_INFO| eps_eigval_S: ', & + qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S IF (qs_env%mp2_env%ri_rpa_im_time%make_chi_pos_definite) THEN WRITE (unit_nr, '(T3,A,T81)') & 'GW_INFO| Make chi(iw,k) positive definite? TRUE' diff --git a/src/rpa_main.F b/src/rpa_main.F index c60b8ed4f9..e8b63fb0f5 100644 --- a/src/rpa_main.F +++ b/src/rpa_main.F @@ -15,18 +15,28 @@ !> 03.2019 Refactoring [Frederick Stein] ! ************************************************************************************************** MODULE rpa_main + USE admm_types, ONLY: admm_env_release USE bibliography, ONLY: & Bates2013, DelBen2013, DelBen2015, Ren2011, Ren2013, Wilhelm2016a, Wilhelm2016b, & Wilhelm2017, Wilhelm2018, cite_reference USE bse, ONLY: do_subspace_iterations,& mult_B_with_W_and_fill_local_3c_arrays + USE cell_types, ONLY: cell_type,& + get_cell USE cp_blacs_env, ONLY: cp_blacs_env_create,& cp_blacs_env_release,& + cp_blacs_env_retain,& cp_blacs_env_type - USE cp_cfm_types, ONLY: cp_cfm_p_type,& + USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm + USE cp_cfm_diag, ONLY: cp_cfm_geeig,& + cp_cfm_geeig_canon + USE cp_cfm_types, ONLY: cp_cfm_create,& + cp_cfm_p_type,& + cp_cfm_release,& + cp_cfm_to_fm,& cp_cfm_type - USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl - USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm + USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& + dbcsr_allocate_matrix_set USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add USE cp_fm_struct, ONLY: cp_fm_struct_create,& cp_fm_struct_release,& @@ -39,12 +49,13 @@ MODULE rpa_main cp_fm_to_fm,& cp_fm_type USE cp_para_env, ONLY: cp_para_env_release,& + cp_para_env_retain,& cp_para_env_split USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: & - dbcsr_add, dbcsr_clear, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, & - dbcsr_get_info, dbcsr_p_type, dbcsr_release, dbcsr_reserve_all_blocks, dbcsr_set, & - dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric + dbcsr_add, dbcsr_clear, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, & + dbcsr_desymmetrize, dbcsr_get_info, dbcsr_p_type, dbcsr_release, dbcsr_reserve_all_blocks, & + dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry USE dbt_api, ONLY: dbt_type USE group_dist_types, ONLY: create_group_dist,& get_group_dist,& @@ -52,14 +63,19 @@ MODULE rpa_main maxsize,& release_group_dist USE hfx_types, ONLY: block_ind_type,& - hfx_compression_type - USE input_constants, ONLY: gw_gf_gamma,& + hfx_compression_type,& + hfx_release + USE input_constants, ONLY: cholesky_off,& + gw_gf_gamma,& gw_gf_mic,& wfc_mm_style_gemm USE kinds, ONLY: dp,& int_8 - USE kpoint_methods, ONLY: rskp_transform + USE kpoint_methods, ONLY: kpoint_env_initialize,& + kpoint_initialize_mo_set,& + kpoint_initialize_mos USE kpoint_types, ONLY: get_kpoint_info,& + kpoint_env_type,& kpoint_type USE machine, ONLY: m_flush,& m_memory @@ -81,12 +97,11 @@ MODULE rpa_main three_dim_real_array,& two_dim_int_array,& two_dim_real_array - USE qs_band_structure, ONLY: calculate_kp_orbitals + USE qs_band_structure, ONLY: calculate_kpoints_for_bs USE qs_environment_types, ONLY: get_qs_env,& - qs_env_release,& qs_environment_type - USE qs_gamma2kp, ONLY: create_kp_from_gamma - USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type + USE qs_mo_types, ONLY: get_mo_set + USE qs_scf_types, ONLY: qs_scf_env_type USE rpa_axk, ONLY: compute_axk_ener USE rpa_gw, ONLY: GW_matrix_operations,& allocate_matrices_gw,& @@ -96,6 +111,7 @@ MODULE rpa_main deallocate_matrices_gw_im_time USE rpa_gw_ic, ONLY: calculate_ic_correction USE rpa_im_time, ONLY: compute_mat_P_omega,& + init_cell_index_rpa,& zero_mat_P_omega USE rpa_im_time_force_methods, ONLY: calc_laplace_loop_forces,& calc_post_loop_forces,& @@ -104,13 +120,16 @@ MODULE rpa_main keep_initial_quad USE rpa_im_time_force_types, ONLY: im_time_force_release,& im_time_force_type - USE rpa_kpoints, ONLY: RPA_postprocessing_kp + USE rpa_kpoints, ONLY: RPA_postprocessing_kp,& + get_mat_cell_T_from_mat_gamma,& + real_space_to_kpoint_transform_rpa USE rpa_util, ONLY: RPA_postprocessing_nokp,& RPA_postprocessing_start,& alloc_im_time,& calc_mat_Q,& contract_P_omega_with_mat_L,& dealloc_im_time + USE scf_control_types, ONLY: scf_control_type USE util, ONLY: get_limit #include "./base/base_uses.f90" @@ -655,9 +674,9 @@ CONTAINS END IF - do_kpoints_from_Gamma = SUM(mp2_env%ri_rpa_im_time%kp_grid) > 0 + do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma IF (do_kpoints_from_Gamma) THEN - CALL get_bandstruc_and_k_dependent_MOs(qs_env, kpoints, virtual(1), Eigenval_kp) + CALL get_bandstruc_and_k_dependent_MOs(qs_env, kpoints, Eigenval_kp) END IF ! Now start the RPA calculation @@ -725,13 +744,11 @@ CONTAINS !> \brief ... !> \param qs_env ... !> \param kpoints ... -!> \param virtual ... !> \param Eigenval_kp ... ! ************************************************************************************************** - SUBROUTINE get_bandstruc_and_k_dependent_MOs(qs_env, kpoints, virtual, Eigenval_kp) + SUBROUTINE get_bandstruc_and_k_dependent_MOs(qs_env, kpoints, Eigenval_kp) TYPE(qs_environment_type), POINTER :: qs_env TYPE(kpoint_type), POINTER :: kpoints - INTEGER :: virtual REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_kp CHARACTER(LEN=*), PARAMETER :: routineN = 'get_bandstruc_and_k_dependent_MOs' @@ -769,7 +786,7 @@ CONTAINS CALL get_qs_env(qs_env=qs_env, para_env=para_env) CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_G, & - "MONKHORST-PACK", virtual, para_env%num_pe, & + "MONKHORST-PACK", para_env%num_pe, & mp_grid=nkp_grid_G(1:3)) IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN @@ -778,11 +795,11 @@ CONTAINS CALL get_kpgeneral_for_Sigma_kpoints(qs_env, kpgeneral) CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, & - "GENERAL", virtual, para_env%num_pe, & + "GENERAL", para_env%num_pe, & kpgeneral=kpgeneral) CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc, & - "GENERAL", virtual, para_env%num_pe, & + "GENERAL", para_env%num_pe, & kpgeneral=kpgeneral, with_xc_terms=.FALSE.) kpoints_Sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma @@ -808,101 +825,270 @@ CONTAINS END IF + CALL dbcsr_release(qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(1)%matrix) + DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(1)%matrix) + DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc) + + CALL dbcsr_release(qs_env%mp2_env%ri_g0w0%matrix_ks(1)%matrix) + DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(1)%matrix) + DEALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks) + + CALL release_hfx_admm_stuff(qs_env) + CALL timestop(handle) END SUBROUTINE get_bandstruc_and_k_dependent_MOs +! ************************************************************************************************** +!> \brief releases part of the given qs_env in order to save memory +!> \param qs_env the object to release +! ************************************************************************************************** + SUBROUTINE release_hfx_admm_stuff(qs_env) + TYPE(qs_environment_type), POINTER :: qs_env + + IF (ASSOCIATED(qs_env%x_data) .AND. .NOT. qs_env%mp2_env%ri_g0w0%do_ri_Sigma_x) THEN + CALL hfx_release(qs_env%x_data) + END IF + IF (ASSOCIATED(qs_env%admm_env)) THEN + CALL admm_env_release(qs_env%admm_env) + END IF + + END SUBROUTINE release_hfx_admm_stuff + ! ************************************************************************************************** !> \brief ... !> \param qs_env ... !> \param kpoints ... !> \param scheme ... -!> \param nadd ... !> \param group_size_ext ... !> \param mp_grid ... !> \param kpgeneral ... !> \param with_xc_terms ... ! ************************************************************************************************** - SUBROUTINE create_kp_and_calc_kp_orbitals(qs_env, kpoints, scheme, nadd, & + SUBROUTINE create_kp_and_calc_kp_orbitals(qs_env, kpoints, scheme, & group_size_ext, mp_grid, kpgeneral, with_xc_terms) TYPE(qs_environment_type), POINTER :: qs_env TYPE(kpoint_type), POINTER :: kpoints CHARACTER(LEN=*), INTENT(IN) :: scheme - INTEGER :: nadd, group_size_ext + INTEGER :: group_size_ext INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: mp_grid REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & OPTIONAL :: kpgeneral LOGICAL, OPTIONAL :: with_xc_terms CHARACTER(LEN=*), PARAMETER :: routineN = 'create_kp_and_calc_kp_orbitals' + COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), & + czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp), ione = CMPLX(0.0_dp, 1.0_dp, KIND=dp) - INTEGER :: handle, ikp, nkp - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + INTEGER :: handle, i_dim, i_re_im, ikp, nkp + INTEGER, DIMENSION(3) :: cell_grid, periodic LOGICAL :: my_with_xc_terms - REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp - TYPE(dbcsr_type), POINTER :: cmat, rmat, tmpmat - TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_nl - TYPE(qs_environment_type), POINTER :: qs_env_kp_from_Gamma + REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues + TYPE(cell_type), POINTER :: cell + TYPE(cp_blacs_env_type), POINTER :: blacs_env + TYPE(cp_cfm_type), POINTER :: cksmat, cmos, csmat, cwork + TYPE(cp_fm_struct_type), POINTER :: matrix_struct + TYPE(cp_fm_type), POINTER :: fm_work, imos, rmos + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_desymm + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_ks_kp, mat_s_kp + TYPE(kpoint_env_type), POINTER :: kp + TYPE(qs_scf_env_type), POINTER :: scf_env + TYPE(scf_control_type), POINTER :: scf_control CALL timeset(routineN, handle) - NULLIFY (qs_env_kp_from_Gamma) - my_with_xc_terms = .TRUE. IF (PRESENT(with_xc_terms)) my_with_xc_terms = with_xc_terms - CALL create_kp_from_gamma(qs_env, qs_env_kp_from_Gamma, with_xc_terms=my_with_xc_terms) + CALL get_qs_env(qs_env, & + para_env=para_env, & + blacs_env=blacs_env, & + matrix_s=matrix_s, & + scf_env=scf_env, & + scf_control=scf_control, & + cell=cell) - ! k-dep. MO coeff. and band structure for Green's function - CALL calculate_kp_orbitals(qs_env_kp_from_Gamma, kpoints, scheme, nadd=nadd, kpgeneral=kpgeneral, & - mp_grid=mp_grid, group_size_ext=group_size_ext) + ! get kpoints + CALL calculate_kpoints_for_bs(kpoints, scheme, kpgeneral=kpgeneral, mp_grid=mp_grid, & + group_size_ext=group_size_ext) - ! Save k-dependent KS-matrix into wmatrix of kpoints - CALL get_qs_env(qs_env_kp_from_Gamma, matrix_ks_kp=matrix_ks_kp) + kpoints%para_env => para_env + CALL cp_para_env_retain(para_env) + kpoints%blacs_env_all => blacs_env + CALL cp_blacs_env_retain(blacs_env) + CALL kpoint_env_initialize(kpoints) - NULLIFY (xkp, sab_nl, cell_to_index) - CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, sab_nl=sab_nl, cell_to_index=cell_to_index) + ! calculate all MOs that are accessible in the given + ! Gaussian AO basis, therefore nadd=1E10 + CALL kpoint_initialize_mos(kpoints, qs_env%mos, 2000000000) + CALL kpoint_initialize_mo_set(kpoints) - ALLOCATE (rmat, cmat, tmpmat) - CALL dbcsr_create(rmat, template=matrix_ks_kp(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) - CALL dbcsr_create(cmat, template=matrix_ks_kp(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric) - CALL dbcsr_create(tmpmat, template=matrix_ks_kp(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry) - CALL cp_dbcsr_alloc_block_from_nbl(rmat, sab_nl) - CALL cp_dbcsr_alloc_block_from_nbl(cmat, sab_nl) - CALL dbcsr_reserve_all_blocks(tmpmat) + CALL get_cell(cell=cell, periodic=periodic) + + DO i_dim = 1, 3 + ! we have at most 3 neigboring cells per dimension and at least one because + ! the density response at Gamma is only divided to neighboring + IF (periodic(i_dim) == 1) THEN + cell_grid(i_dim) = MAX(MIN((kpoints%nkp_grid(i_dim)/2)*2 - 1, 1), 3) + ELSE + cell_grid(i_dim) = 1 + END IF + END DO + CALL init_cell_index_rpa(cell_grid, kpoints%cell_to_index, kpoints%index_to_cell, cell) + + ! get H(k) + IF (my_with_xc_terms) THEN + CALL mat_kp_from_mat_gamma(qs_env, mat_ks_kp, qs_env%mp2_env%ri_g0w0%matrix_ks(1)%matrix, kpoints) + ELSE + CALL mat_kp_from_mat_gamma(qs_env, mat_ks_kp, qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(1)%matrix, kpoints) + END IF + + ! get S(k) + CALL get_qs_env(qs_env, matrix_s=matrix_s, scf_env=scf_env, scf_control=scf_control) + + NULLIFY (matrix_s_desymm) + CALL dbcsr_allocate_matrix_set(matrix_s_desymm, 1) + ALLOCATE (matrix_s_desymm(1)%matrix) + CALL dbcsr_create(matrix=matrix_s_desymm(1)%matrix, template=matrix_s(1)%matrix, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_desymmetrize(matrix_s(1)%matrix, matrix_s_desymm(1)%matrix) + + CALL mat_kp_from_mat_gamma(qs_env, mat_s_kp, matrix_s_desymm(1)%matrix, kpoints) + + CALL get_kpoint_info(kpoints, nkp=nkp) + + matrix_struct => kpoints%kp_env(1)%kpoint_env%wmat(1, 1)%matrix%matrix_struct + + CALL cp_cfm_create(cksmat, matrix_struct) + CALL cp_cfm_create(csmat, matrix_struct) + CALL cp_cfm_create(cmos, matrix_struct) + CALL cp_cfm_create(cwork, matrix_struct) + CALL cp_fm_create(fm_work, matrix_struct) DO ikp = 1, nkp - CALL dbcsr_set(rmat, 0.0_dp) - CALL dbcsr_set(cmat, 0.0_dp) + CALL copy_dbcsr_to_fm(mat_ks_kp(ikp, 1)%matrix, kpoints%kp_env(ikp)%kpoint_env%wmat(1, 1)%matrix) + CALL cp_cfm_scale_and_add_fm(czero, cksmat, cone, kpoints%kp_env(ikp)%kpoint_env%wmat(1, 1)%matrix) - ! trafo matrix_ks from real space to k-space and write into wmat - CALL rskp_transform(rmatrix=rmat, cmatrix=cmat, & - rsmat=matrix_ks_kp, ispin=1, xkp=xkp(1:3, ikp), & - cell_to_index=cell_to_index, sab_nl=sab_nl) + CALL copy_dbcsr_to_fm(mat_ks_kp(ikp, 2)%matrix, kpoints%kp_env(ikp)%kpoint_env%wmat(2, 1)%matrix) + CALL cp_cfm_scale_and_add_fm(cone, cksmat, ione, kpoints%kp_env(ikp)%kpoint_env%wmat(2, 1)%matrix) - CALL dbcsr_set(tmpmat, 0.0_dp) - CALL dbcsr_desymmetrize(rmat, tmpmat) + CALL copy_dbcsr_to_fm(mat_s_kp(ikp, 1)%matrix, fm_work) + CALL cp_cfm_scale_and_add_fm(czero, csmat, cone, fm_work) - CALL cp_fm_set_all(kpoints%kp_env(ikp)%kpoint_env%wmat(1, 1)%matrix, 0.0_dp) - CALL copy_dbcsr_to_fm(tmpmat, kpoints%kp_env(ikp)%kpoint_env%wmat(1, 1)%matrix) + CALL copy_dbcsr_to_fm(mat_s_kp(ikp, 2)%matrix, fm_work) + CALL cp_cfm_scale_and_add_fm(cone, csmat, ione, fm_work) - CALL dbcsr_set(tmpmat, 0.0_dp) - CALL dbcsr_desymmetrize(cmat, tmpmat) + kp => kpoints%kp_env(ikp)%kpoint_env - CALL cp_fm_set_all(kpoints%kp_env(ikp)%kpoint_env%wmat(2, 1)%matrix, 0.0_dp) - CALL copy_dbcsr_to_fm(tmpmat, kpoints%kp_env(ikp)%kpoint_env%wmat(2, 1)%matrix) + CALL get_mo_set(kp%mos(1, 1)%mo_set, mo_coeff=rmos, eigenvalues=eigenvalues) + CALL get_mo_set(kp%mos(2, 1)%mo_set, mo_coeff=imos) + + IF (scf_env%cholesky_method == cholesky_off) THEN + CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, scf_control%eps_eigval) + ELSE + CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork) + END IF + + CALL cp_cfm_to_fm(cmos, rmos, imos) + + kp%mos(2, 1)%mo_set%eigenvalues = eigenvalues END DO - CALL qs_env_release(qs_env_kp_from_Gamma) + DO ikp = 1, nkp + DO i_re_im = 1, 2 + CALL dbcsr_deallocate_matrix(mat_ks_kp(ikp, i_re_im)%matrix) + END DO + END DO + DEALLOCATE (mat_ks_kp) - CALL dbcsr_deallocate_matrix(rmat) - CALL dbcsr_deallocate_matrix(cmat) - CALL dbcsr_deallocate_matrix(tmpmat) + DO ikp = 1, nkp + DO i_re_im = 1, 2 + CALL dbcsr_deallocate_matrix(mat_s_kp(ikp, i_re_im)%matrix) + END DO + END DO + DEALLOCATE (mat_s_kp) + + CALL dbcsr_deallocate_matrix(matrix_s_desymm(1)%matrix) + DEALLOCATE (matrix_s_desymm) + + CALL cp_cfm_release(cksmat) + CALL cp_cfm_release(csmat) + CALL cp_cfm_release(cwork) + CALL cp_cfm_release(cmos) + CALL cp_fm_release(fm_work) + + CALL timestop(handle) + + END SUBROUTINE create_kp_and_calc_kp_orbitals + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param mat_kp ... +!> \param mat_gamma ... +!> \param kpoints ... +! ************************************************************************************************** + SUBROUTINE mat_kp_from_mat_gamma(qs_env, mat_kp, mat_gamma, kpoints) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_kp + TYPE(dbcsr_type) :: mat_gamma + TYPE(kpoint_type), POINTER :: kpoints + + CHARACTER(LEN=*), PARAMETER :: routineN = 'mat_kp_from_mat_gamma' + + INTEGER :: handle, i_cell, i_re_im, ikp, nkp, & + num_cells + INTEGER, DIMENSION(3) :: periodic + INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp + TYPE(cell_type), POINTER :: cell + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_real_space + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env, cell=cell) + CALL get_cell(cell=cell, periodic=periodic) + num_cells = 3**(periodic(1) + periodic(2) + periodic(3)) + + NULLIFY (mat_real_space) + CALL dbcsr_allocate_matrix_set(mat_real_space, num_cells) + DO i_cell = 1, num_cells + ALLOCATE (mat_real_space(i_cell)%matrix) + CALL dbcsr_create(matrix=mat_real_space(i_cell)%matrix, & + template=mat_gamma) + CALL dbcsr_reserve_all_blocks(mat_real_space(i_cell)%matrix) + CALL dbcsr_set(mat_real_space(i_cell)%matrix, 0.0_dp) + END DO + + CALL dbcsr_copy(mat_real_space(1)%matrix, mat_gamma) + + CALL get_mat_cell_T_from_mat_gamma(mat_real_space, qs_env, kpoints, 2, 0) + + NULLIFY (xkp, cell_to_index) + CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, cell_to_index=cell_to_index) + + NULLIFY (mat_kp) + CALL dbcsr_allocate_matrix_set(mat_kp, nkp, 2) + DO ikp = 1, nkp + DO i_re_im = 1, 2 + ALLOCATE (mat_kp(ikp, i_re_im)%matrix) + CALL dbcsr_create(matrix=mat_kp(ikp, i_re_im)%matrix, template=mat_gamma) + CALL dbcsr_reserve_all_blocks(mat_kp(ikp, i_re_im)%matrix) + CALL dbcsr_set(mat_kp(ikp, i_re_im)%matrix, 0.0_dp) + END DO + END DO + + CALL real_space_to_kpoint_transform_rpa(mat_kp(:, 1), mat_kp(:, 2), mat_real_space, kpoints, 0.0_dp) + + DO i_cell = 1, num_cells + CALL dbcsr_deallocate_matrix(mat_real_space(i_cell)%matrix) + END DO + DEALLOCATE (mat_real_space) CALL timestop(handle) diff --git a/src/xas_tdp_atom.F b/src/xas_tdp_atom.F index 7144491dfd..e0576656df 100644 --- a/src/xas_tdp_atom.F +++ b/src/xas_tdp_atom.F @@ -3232,7 +3232,7 @@ CONTAINS maxso=maxso, npgf=npgf, nset=nset, zet=zet) ! Separate the functions into purely r and purely angular parts, compute them all -! and use matrix mutliplication for the integral. We use f for x derivative ang g for y +! and use matrix mutliplication for the integral. We use f for x derivative and g for y ! Separating the functions. Note that the radial part is the same for x and y derivatives ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3)) diff --git a/tests/QS/regtest-gw-cubic/TEST_FILES b/tests/QS/regtest-gw-cubic/TEST_FILES index 202ec0dc92..df904b6047 100644 --- a/tests/QS/regtest-gw-cubic/TEST_FILES +++ b/tests/QS/regtest-gw-cubic/TEST_FILES @@ -7,8 +7,5 @@ scGW0_and_evGW_H2O_PBE0_trunc_minimax.inp 11 1e-08 - scGW0_H2O_PBE0_trunc_minimax_RI_HFX.inp 11 1e-08 -13.191093644224157 G0W0_H2O_PBE_periodic.inp 78 1e-05 16.475 G0W0_OH_PBE.inp 79 1e-05 11.797 -G0W0_kpoints_from_Gamma.inp 78 1e-05 14.306 -G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp 78 1e-05 14.682 G0W0_OH_PBE_svd.inp 79 1e-05 11.797 -G0W0_kpoints_in_self_energy.inp 98 1e-05 12.091 #EOF diff --git a/tests/QS/regtest-gw-cubic/G0W0_kpoints_from_Gamma.inp b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_from_Gamma.inp similarity index 100% rename from tests/QS/regtest-gw-cubic/G0W0_kpoints_from_Gamma.inp rename to tests/QS/regtest-gw-kpoints/G0W0_kpoints_from_Gamma.inp diff --git a/tests/QS/regtest-gw-cubic/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp similarity index 100% rename from tests/QS/regtest-gw-cubic/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp rename to tests/QS/regtest-gw-kpoints/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp diff --git a/tests/QS/regtest-gw-cubic/G0W0_kpoints_in_self_energy.inp b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy.inp similarity index 100% rename from tests/QS/regtest-gw-cubic/G0W0_kpoints_in_self_energy.inp rename to tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy.inp diff --git a/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX.inp b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX.inp new file mode 100644 index 0000000000..f09c411dd2 --- /dev/null +++ b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX.inp @@ -0,0 +1,105 @@ +&GLOBAL + PROJECT G0W0_kpoints_from_Gamma + PRINT_LEVEL MEDIUM + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.01 + &END +&END GLOBAL +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + SORT_BASIS EXP + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &QS + METHOD GPW + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-15 + &END QS + &SCF + SCF_GUESS RESTART + EPS_SCF 1.0E-5 + MAX_SCF 100 + &PRINT + &RESTART ON + &END + &END + ADDED_MOS -1 + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &WF_CORRELATION + &INTEGRALS + SIZE_LATTICE_SUM 3 + &END INTEGRALS + &LOW_SCALING + KPOINTS 4 1 4 + &END LOW_SCALING + &RI_RPA + RPA_NUM_QUAD_POINTS 6 + &GW + CORR_OCC 1 + CORR_VIRT 1 + CROSSING_SEARCH NEWTON + ANALYTIC_CONTINUATION TWO_POLE + RI_SIGMA_X FALSE + KPOINTS_SELF_ENERGY 2 2 1 + &KPOINT_SET + SPECIAL_POINT 0.5 0.0 0.0 + SPECIAL_POINT 0.0 0.0 0.0 + NPOINTS 3 + &END + &END GW + &HF + FRACTION 1.0 + &SCREENING + EPS_SCHWARZ 1.0E-6 + SCREEN_ON_INITIAL_P TRUE + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 3.9 + T_C_G_DATA t_c_g.dat + &END + &MEMORY + ! In MB per MPI rank.. use as much as need to get in-core operation + MAX_MEMORY 0 + EPS_STORAGE_SCALING 0.1 + &END + &END + &END RI_RPA + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 8.000 8.000 8.000 + MULTIPLE_UNIT_CELL 1 1 1 + PERIODIC XZ + &END CELL + &KIND H + BASIS_SET ORB DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND O + BASIS_SET ORB DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q6 + &END KIND + &TOPOLOGY + MULTIPLE_UNIT_CELL 1 1 1 + &END TOPOLOGY + &COORD + H 0.0 -0.5 -4.5 + O 0.5 0.0 4.5 + H 0.0 0.5 -4.5 + &END COORD + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX_ADMM.inp b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX_ADMM.inp new file mode 100644 index 0000000000..9441913882 --- /dev/null +++ b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX_ADMM.inp @@ -0,0 +1,130 @@ +&GLOBAL + PROJECT G0W0_kpoints_from_Gamma + PRINT_LEVEL MEDIUM + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.01 + &END +&END GLOBAL +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + BASIS_SET_FILE_NAME BASIS_ADMM + SORT_BASIS EXP + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &QS + METHOD GPW + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-15 + &END QS + &SCF + SCF_GUESS RESTART + EPS_SCF 1.0E-5 + MAX_SCF 100 + &PRINT + &RESTART ON + &END + &END + ADDED_MOS -1 + &END SCF + &AUXILIARY_DENSITY_MATRIX_METHOD + METHOD BASIS_PROJECTION + ADMM_PURIFICATION_METHOD NONE + &END AUXILIARY_DENSITY_MATRIX_METHOD + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &HF + FRACTION 0.0 + &SCREENING + EPS_SCHWARZ 1.0E-3 + SCREEN_ON_INITIAL_P TRUE + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 3.0 + T_C_G_DATA t_c_g.dat + &END + &MEMORY + ! In MB per MPI rank.. use as much as need to get in-core operation + MAX_MEMORY 0 + EPS_STORAGE_SCALING 0.1 + &END + &END + &WF_CORRELATION + &INTEGRALS + SIZE_LATTICE_SUM 3 + &END INTEGRALS + &LOW_SCALING + KPOINTS 4 1 4 + &END LOW_SCALING + &RI_RPA + RPA_NUM_QUAD_POINTS 6 + ADMM + &GW + CORR_OCC 1 + CORR_VIRT 1 + CROSSING_SEARCH NEWTON + ANALYTIC_CONTINUATION TWO_POLE + RI_SIGMA_X FALSE + KPOINTS_SELF_ENERGY 2 2 1 + &KPOINT_SET + SPECIAL_POINT 0.5 0.0 0.0 + SPECIAL_POINT 0.0 0.0 0.0 + NPOINTS 3 + &END + &END GW + &HF + FRACTION 1.0 + &SCREENING + EPS_SCHWARZ 1.0E-6 + SCREEN_ON_INITIAL_P TRUE + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 3.9 + T_C_G_DATA t_c_g.dat + &END + &MEMORY + ! In MB per MPI rank.. use as much as need to get in-core operation + MAX_MEMORY 0 + EPS_STORAGE_SCALING 0.1 + &END + &END + &END RI_RPA + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 8.000 8.000 8.000 + MULTIPLE_UNIT_CELL 1 1 1 + PERIODIC XZ + &END CELL + &KIND H + BASIS_SET ORB DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + BASIS_SET AUX_FIT cFIT3 + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND O + BASIS_SET ORB DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + BASIS_SET AUX_FIT cFIT3 + POTENTIAL GTH-PBE-q6 + &END KIND + &TOPOLOGY + MULTIPLE_UNIT_CELL 1 1 1 + &END TOPOLOGY + &COORD + H 0.0 -0.5 -4.5 + O 0.5 0.0 4.5 + H 0.0 0.5 -4.5 + &END COORD + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rpa-cubic-scaling/RPA_kpoints_from_Gamma_H2O.inp b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_at_HSE06.inp similarity index 56% rename from tests/QS/regtest-rpa-cubic-scaling/RPA_kpoints_from_Gamma_H2O.inp rename to tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_at_HSE06.inp index 429f07fefb..289dfb731f 100644 --- a/tests/QS/regtest-rpa-cubic-scaling/RPA_kpoints_from_Gamma_H2O.inp +++ b/tests/QS/regtest-gw-kpoints/G0W0_kpoints_in_self_energy_at_HSE06.inp @@ -1,5 +1,5 @@ &GLOBAL - PROJECT RPA_H2O_kpoints + PROJECT kp_GW PRINT_LEVEL MEDIUM RUN_TYPE ENERGY &TIMINGS @@ -25,36 +25,63 @@ SCF_GUESS ATOMIC EPS_SCF 1.0E-5 MAX_SCF 100 - &PRINT - &RESTART OFF - &END - &END - ADDED_MOS 10000000 + ADDED_MOS -1 &END SCF &XC - &XC_FUNCTIONAL PBE + &XC_FUNCTIONAL + &XWPBE + SCALE_X -0.25 + SCALE_X0 1.0 + OMEGA 0.11 + &END + &PBE + SCALE_X 0.0 + SCALE_C 1.0 + &END PBE &END XC_FUNCTIONAL + &HF + &SCREENING + EPS_SCHWARZ 1.0E-10 + &END + &INTERACTION_POTENTIAL + POTENTIAL_TYPE SHORTRANGE + OMEGA 0.11 + &END + &MEMORY + MAX_MEMORY 10 + &END + FRACTION 0.25 + &END &WF_CORRELATION &INTEGRALS SIZE_LATTICE_SUM 3 &END INTEGRALS - MEMORY 200. - NUMBER_PROC 1 &LOW_SCALING - MEMORY_CUT 1 - KPOINTS 4 4 4 + KPOINTS 1 4 4 + REGULARIZATION_RI 1.0E-3 &END LOW_SCALING &RI_RPA RPA_NUM_QUAD_POINTS 6 + &GW + CORR_OCC 1 + CORR_VIRT 1 + RI_SIGMA_X + KPOINTS_SELF_ENERGY 2 2 1 + &KPOINT_SET + SPECIAL_POINT 0.5 0.0 0.0 + SPECIAL_POINT 0.0 0.0 0.0 + NPOINTS 3 + &END + &END GW &END RI_RPA &END &END XC &END DFT &SUBSYS &CELL - ABC [angstrom] 10.000 10.000 10.000 + ABC [angstrom] 8.000 8.000 8.000 MULTIPLE_UNIT_CELL 1 1 1 - PERIODIC XYZ + PERIODIC YZ &END CELL &KIND H BASIS_SET ORB DZVP-GTH @@ -70,9 +97,9 @@ MULTIPLE_UNIT_CELL 1 1 1 &END TOPOLOGY &COORD - H 0.0 -0.5 -5.5 - O 0.5 0.0 5.5 - H 0.0 0.5 -5.5 + H 0.0 -0.5 -4.5 + O 0.5 0.0 4.5 + H 0.0 0.5 -4.5 &END COORD &END SUBSYS &END FORCE_EVAL diff --git a/tests/QS/regtest-gw-kpoints/TEST_FILES b/tests/QS/regtest-gw-kpoints/TEST_FILES new file mode 100644 index 0000000000..b758036a4b --- /dev/null +++ b/tests/QS/regtest-gw-kpoints/TEST_FILES @@ -0,0 +1,7 @@ +G0W0_kpoints_from_Gamma.inp 78 1e-05 14.753 +G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp 78 1e-05 14.682 +G0W0_kpoints_in_self_energy.inp 98 1e-05 12.090 +G0W0_kpoints_in_self_energy_at_HSE06.inp 98 1e-05 12.569 +G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX.inp 98 1e-05 15.438 +G0W0_kpoints_in_self_energy_Sigmax_from_four_center_HFX_ADMM.inp 98 1e-05 15.484 +#EOF diff --git a/tests/QS/regtest-rpa-cubic-scaling/TEST_FILES b/tests/QS/regtest-rpa-cubic-scaling/TEST_FILES index 25a44a9c5c..241948e44e 100644 --- a/tests/QS/regtest-rpa-cubic-scaling/TEST_FILES +++ b/tests/QS/regtest-rpa-cubic-scaling/TEST_FILES @@ -3,9 +3,8 @@ Cubic_RPA_H2O_check_filtering.inp 11 1e-08 - Cubic_RPA_2x_H2_check_filtering.inp 11 1e-08 -2.119159936957426 Cubic_RPA_CH3.inp 11 3e-07 -7.414901056476346 Cubic_RPA_H2O_standard_svd.inp 11 1e-08 -17.160104846594532 -RPA_kpoints_H2O.inp 11 1e-08 -17.384092112689462 -RPA_kpoints_H2O_batched.inp 11 1e-08 -17.384092114445451 -RPA_kpoints_from_Gamma_H2O.inp 11 1e-08 -17.384093823208083 +RPA_kpoints_H2O.inp 11 1e-08 -17.384612021546573 +RPA_kpoints_H2O_batched.inp 11 1e-08 -17.384612040108440 Cubic_RPA_CH3_ri-hfx.inp 11 1e-08 -7.414775787676669 Cubic_RPA_H2O_ri-hfx.inp 11 1e-08 -17.159865017693086 #EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index d4965089cd..696fc606a4 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -118,6 +118,7 @@ QS/regtest-nmr-uks-1 QS/regtest-libxc libxc libint SE/regtest-3-2 QS/regtest-gw-cubic libint +QS/regtest-gw-kpoints libint QS/regtest-xc QS/regtest-almo-2 SE/regtest-2-2